WMTransferCalc.cpp 64.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenWalnut is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------

// Some rules to the inclusion of headers:
//  * Ordering:
//    * C Headers
//    * C++ Standard headers
//    * External Lib headers (like OSG or Boost headers)
//    * headers of other classes inside OpenWalnut
//    * your own header file
32
#define _USE_MATH_DEFINES
33

34
#include <math.h>
35
#include <stdio.h>
36
#include <stdlib.h>
37
#include <time.h>
38 39

#include <algorithm>
40
#include <cmath>
41 42
#include <cstdlib>
#include <ctime>
Franziska Jacob's avatar
Franziska Jacob committed
43
#include <fstream>
44 45 46 47 48
#include <iostream>
#include <limits>
#include <sstream>
#include <string>
#include <vector>
49 50 51 52 53 54

#include <osg/ShapeDrawable>
#include <osg/Group>
#include <osg/Geode>
#include <osg/Material>
#include <osg/StateAttribute>
55
#include <osg/Vec3>
56

57 58
#include <Eigen/Eigen>

59

60 61 62 63
#include "core/kernel/WKernel.h"
#include "core/common/WColor.h"
#include "core/common/WPathHelper.h"
#include "core/common/WPropertyHelper.h"
64
#include "core/common/math/WMatrix.h"
65 66
#include "core/common/math/linearAlgebra/WMatrixFixed.h"
#include "core/common/math/linearAlgebra/WVectorFixed.h"
67 68 69 70
#include "core/common/WItemSelection.h"
#include "core/common/WItemSelectionItem.h"
#include "core/common/WItemSelectionItemTyped.h"
#include "core/common/WItemSelector.h"
71 72 73
#include "core/graphicsEngine/WGEUtils.h"
#include "core/graphicsEngine/WGERequirement.h"

74 75 76
#include "../interaction/WViewEventHandler.h"
#include "../interaction/WGetMatrixCallback.h"

77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
#include "WMTransferCalc.xpm"
#include "WMTransferCalc.h"

WMTransferCalc::WMTransferCalc():
    WModule()
{
    // In the constructor, you can initialize your members and all this stuff. You must not initialize connectors or properties here! You also
    // should avoid doing computationally expensive stuff, since every module has its own thread which is intended to be used for such calculations.
    // Please keep in mind, that every member initialized here is also initialized in the prototype, which may be a problem if the member is large,
    // and therefore, wasting a lot of memory in your module's prototype instance.
}

WMTransferCalc::~WMTransferCalc()
{
    // Cleanup!
}

boost::shared_ptr< WModule > WMTransferCalc::factory() const
{
    return boost::shared_ptr< WModule >( new WMTransferCalc() );
}

const char** WMTransferCalc::getXPMIcon() const
{
    return WMTransferCalc_xpm;
}

const std::string WMTransferCalc::getName() const
{
106
    return "Transfer Function Calculation";
107 108 109 110
}

const std::string WMTransferCalc::getDescription() const
{
111 112 113
    return "Calculates optimal settings for the Transfer Function after klicking " \
            "by casting rays through the given data and analysing their profiles, " \
            "e. g. with regard to value, gradient or gradient weight.";
114 115 116 117
}

void WMTransferCalc::connectors()
{
118 119
    m_inputData = WModuleInputData< WDataSetScalar >::createAndAdd( shared_from_this(), "Compulsory Data",
                                                                    "Data which will be checked for displaying." );
120 121
    m_inputDerivation = WModuleInputData< WDataSetVector >::createAndAdd( shared_from_this(), "Derivated Data",
                                                                    "Derivated data by using the module 'Spatial Derivative'." );
122 123
    m_inputFA = WModuleInputData< WDataSetScalar >::createAndAdd( shared_from_this(), "Fractional Anisotropy Data (optional)",
                                                                    "Input for additional data concerning fractional anisotropy of a dataset." );
124

125 126 127 128 129 130
    // output data
    m_curveMeanOut = WModuleOutputData< WDataSetScalar >::createAndAdd( shared_from_this(),
                                                                        "(a+b)/2", "The calculated curvature dataset in mean form." );
    m_curveGaussOut = WModuleOutputData< WDataSetScalar >::createAndAdd( shared_from_this(),
                                                                         "a*b", "The calculated curvature dataset in Gauss form." );

131 132 133 134 135 136
    // call WModules initialization
    WModule::connectors();
}

void WMTransferCalc::properties()
{
Franziska Jacob's avatar
Franziska Jacob committed
137
    m_propCondition   = boost::shared_ptr< WCondition >( new WCondition() );
138 139
    m_color           = m_properties->addProperty(  "Color", "Color of the ray.", WColor( 1.0, 0.0, 0.0, 1.0 ), m_propCondition );

140
    m_interval        = m_properties->addProperty( "Interval", "Interval for sampling rays.", 0.30, m_propCondition );
141 142 143 144
    m_rayNumber       = m_properties->addProperty( "# of Rays", "Number of rays being casted with one click.", 10, m_propCondition );
    m_radius          = m_properties->addProperty( "Radius", "Radius for circular vicinity in which the rays shall be casted.",
                                                   2, m_propCondition );

145 146 147 148 149 150 151
    m_plotDataSelection = WItemSelection::SPtr( new WItemSelection() );
    m_plotDataSelection->addItem( MyItemType::create( 3, "Plot value.", "Chooses 'value' data for plotting." ) );
    m_plotDataSelection->addItem( MyItemType::create( 7, "Plot gradient weight.", "Chooses 'gradient weight' data for plotting." ) );
    m_plotDataSelection->addItem( MyItemType::create( 8, "Plot FA value.", "Chooses 'fractional anisotropy' data for plotting." ) );
    m_plotDataSelection->addItem( MyItemType::create( 9, "Plot angle.", "Chooses 'angle between ray direction and gradient' for plotting." ) );
    m_plotDataSelection->addItem( MyItemType::create( 10, "Plot mean curvature.", "Chooses 'mean curvature' data for plotting." ) );
    m_plotDataSelection->addItem( MyItemType::create( 11, "Plot gaussian curvature.", "Chooses 'Gaussian curvature' data for plotting." ) );
152

153 154
    m_multiPlotDataSelection = m_properties->addProperty( "Choose plot data", "Choose one of these for data plotting.",
                                                              m_plotDataSelection->getSelectorFirst(), m_propCondition );
155

Franziska Jacob's avatar
Franziska Jacob committed
156

157 158
    m_RaySaveFilePath = m_properties->addProperty( "Save RayProfile to Folder:",
                                                   "Savefile for RayProfile.", WPathHelper::getAppPath(), m_propCondition );
Franziska Jacob's avatar
Franziska Jacob committed
159 160 161 162
    m_saveTrigger     = m_properties->addProperty( "Save", "Save RayProfile to given directory.", WPVBaseTypes::PV_TRIGGER_READY, m_propCondition );

    WPropertyHelper::PC_PATHEXISTS::addTo( m_RaySaveFilePath );
    WPropertyHelper::PC_ISDIRECTORY::addTo( m_RaySaveFilePath );
163
    WPropertyHelper::PC_NOTEMPTY::addTo( m_multiPlotDataSelection );
164 165 166 167 168 169 170 171 172

    WModule::properties();
}

void WMTransferCalc::requirements()
{
    m_requirements.push_back( new WGERequirement() );
}

Franziska Jacob's avatar
Franziska Jacob committed
173 174 175
bool WMTransferCalc::saveRayProfileTo( WPropFilename path, std::string filename, WRayProfile *profile )
{
    WRayProfile *prof;
176 177 178 179 180 181 182 183 184 185 186 187
    if( profile == NULL )
    {
        // if no profile is given, the main profile will be used as default
        prof = &m_mainProfile;
    }
    else
    {
        // use given profile
        prof = profile;
    }
    //( profile == NULL ) ? prof = &m_mainProfile : prof = profile;

188 189
    if( prof->size() == 0 )
    {
190
        return false; //  empty profiles will not be saved
191
    }
Franziska Jacob's avatar
Franziska Jacob committed
192

193
    // compose context of the soon-to-be-saved file which will be saved in given path with given filename
194
    std::ofstream savefile( ( path->get( true ).string() + "/" + filename ).c_str() /*, ios::app */ );
195
    // heading
Franziska Jacob's avatar
Franziska Jacob committed
196
    savefile << "# RayProfile data for OpenWalnut TransferCalcModule" << std::endl;
197
    savefile << "# ID \t dist \t value \t grad_x \t grad_y \t grad_z \t gradW \t FA \t angle \t mean \t gauss" << std::endl;
198
    // save data of each sample point
Franziska Jacob's avatar
Franziska Jacob committed
199 200
    for( size_t sample = 0; sample < prof->size(); sample++ )
    {
201
        savefile << sample + 1 << "\t";                         // print id
202 203
        savefile << (*prof)[sample].distance() << "\t";         // print t
        savefile << (*prof)[sample].value() << "\t";            // print value
204 205 206 207
        savefile << (*prof)[sample].gradient()[0] << "\t";      // print gradient x
        savefile << (*prof)[sample].gradient()[1] << "\t";      // print gradient y
        savefile << (*prof)[sample].gradient()[2] << "\t";      // print gradient z
        savefile << (*prof)[sample].gradWeight() << "\t";       // print weight of gradient
208
        savefile << (*prof)[sample].fracA() << "\t";            // print fractional anisotropy
209
        savefile << (*prof)[sample].angle() << "\t";            // print angle between gradient and ray
210 211
        savefile << (*prof)[sample].meanCurv() << "\t";         // print mean curvature
        savefile << (*prof)[sample].gaussCurv() << "\t";        // print gauss curvature
212
        savefile << std::endl;
Franziska Jacob's avatar
Franziska Jacob committed
213 214 215 216 217
    }
    savefile.close();
    return true;
}

218 219 220 221
void WMTransferCalc::moduleMain()
{
    debugLog() << "moduleMain started...";

222
    // init observation of properties
223 224
    m_moduleState.setResetable( true, true );
    m_moduleState.add( m_inputData->getDataChangedCondition() );
225
    m_moduleState.add( m_inputDerivation->getDataChangedCondition() );
226
    m_moduleState.add( m_inputFA->getDataChangedCondition() );
227 228 229 230 231 232 233 234 235
    m_moduleState.add( m_propCondition );

    ready();
    debugLog() << "Module is now ready.";

    m_rootNode = new WGEManagedGroupNode( m_active );
    //m_rootNode->addUpdateCallback( new TranslateCallback( this ) );
    WKernel::getRunningKernel()->getGraphicsEngine()->getScene()->insert( m_rootNode );

236 237 238 239 240 241 242 243 244
    // Create and add an event handler for the main view as we are interested in click events
    WViewEventHandler::RefPtr handler( new WViewEventHandler() );
    WKernel::getRunningKernel()->getGraphicsEngine()->getViewer()->getView()->addEventHandler( handler );
    handler->onLeftClick( boost::bind( &WMTransferCalc::onClick, this, _1 ) );

    // as we need the projection and view matrices: add a callback which queries this
    m_matrixCallback = new WGetMatrixCallback();
    m_rootNode->addCullCallback( m_matrixCallback );

245 246 247 248 249
    debugLog() << "Entering main loop";
    while( !m_shutdownFlag() )
    {
        debugLog() << "Waiting ...";
        m_moduleState.wait();
250
        debugLog() << "New Data Incoming ...";
251 252 253 254 255 256 257

        // woke up since the module is requested to finish
        if( m_shutdownFlag() )
        {
            debugLog() << "Closing module...";
            break;
        }
258

Franziska Jacob's avatar
Franziska Jacob committed
259 260 261
        // activates if user pushes save button to save current RayProfile
        if( m_saveTrigger->get( true ) == WPVBaseTypes::PV_TRIGGER_TRIGGERED )
        {
262
            debugLog() << "User saves RayProfile(s).";
Franziska Jacob's avatar
Franziska Jacob committed
263

264
            // start progress for GUI
265
            boost::shared_ptr< WProgress > saving = boost::shared_ptr< WProgress >( new WProgress( "Saving calculated RayProfiles." ) );
Franziska Jacob's avatar
Franziska Jacob committed
266 267
            m_progress->addSubProgress( saving );

268
            // save data and plot file for every profile
269 270
            unsigned int nrProfiles = m_profiles.size();
            for( unsigned int id = 0; id < nrProfiles; id++ )
271
            {
272 273
                // generate header
                // first a general header for setting default options
274 275 276 277 278 279 280
                std::stringstream prepGNUPlotScript;
                prepGNUPlotScript << "# Measured RayProfile data for Profile "<< id + 1 << std::endl
                                << "#===========================================" << std::endl << std::endl
                                << "set title 'RayProfile " << id + 1 << "'" << std::endl
                                << "set pointsize 0.5" << std::endl
                                << "set grid" << std::endl
                                << "set xlabel 'uniform distance on ray'" << std::endl;
281

282
                // generate filename of data file
283
                std::stringstream prepFilename;
284
                prepFilename << "RayProfile" << id + 1 << ".dat";
285
                std::string filename( prepFilename.str() );
Sebastian Eichelbaum's avatar
[MERGE]  
Sebastian Eichelbaum committed
286

287
                // save data which will be plotted with the plot script with this method
288
                bool success = saveRayProfileTo( m_RaySaveFilePath, filename, &m_profiles[id] );
289 290
                if( !success )
                {
291
                    errorLog() << "RayProfile " << id + 1 << " was empty or could not be saved.";
292
                }
293
                else // saving successful
294
                {
295
                    // check selector for the chosen data to plot
296 297 298
                    WItemSelector s = m_multiPlotDataSelection->get( true );
                    unsigned int count = s.size();
                    prepGNUPlotScript << "plot ";
299
                    // each selected option gets an individual command line in the later script
300
                    for( size_t i = 0; i < s.size(); ++i )
301
                    {
302 303
//                         infoLog() << "The user chose " << s.at( i )->getName()
//                                 << " with the value " << s.at( i )->getAs< MyItemType >()->getValue();
304

305
                        // switch works with an enumeration
306 307 308
                        switch( s.at( i )->getAs< MyItemType >()->getValue() )
                        {
                            case VALUE:
309
                                prepGNUPlotScript << "'" << filename << "' using 2:3 title 'value' with linespoints";
310 311 312 313 314 315 316
                                if( count - 1 != 0 )
                                {
                                    prepGNUPlotScript << ", \\" << std::endl;
                                }
                                count--;
                                break;
                            case WEIGHT:
317
                                prepGNUPlotScript << "'" << filename << "' using 2:7 title 'gradient weight' with linespoints";
318 319 320 321 322 323 324
                                if( count - 1 != 0 )
                                {
                                    prepGNUPlotScript << ", \\" << std::endl;
                                }
                                count--;
                                break;
                            case FA:
325
                                prepGNUPlotScript << "'" << filename << "' using 2:($8*100) title 'FA value' with linespoints";
326 327 328 329 330 331 332
                                if( count - 1 != 0 )
                                {
                                    prepGNUPlotScript << ", \\" << std::endl;
                                }
                                count--;
                                break;
                            case ANGLE:
333
                                prepGNUPlotScript << "'" << filename << "' using 2:($9*(180/3.14159265)) title 'angle' with linespoints";
334 335 336 337 338 339 340
                                if( count - 1 != 0 )
                                {
                                    prepGNUPlotScript << ", \\" << std::endl;
                                }
                                count--;
                                break;
                            case MEAN:
341
                                prepGNUPlotScript << "'" << filename << "' using 2:($10*100) title 'mean curvature' with linespoints";
342 343 344 345 346 347 348
                                if( count - 1 != 0 )
                                {
                                    prepGNUPlotScript << ", \\" << std::endl;
                                }
                                count--;
                                break;
                            case GAUSS:
349
                                prepGNUPlotScript << "'" << filename << "' using 2:($11*100) title 'gaussian curvature' with linespoints";
350 351 352 353 354 355 356 357 358
                                if( count - 1 != 0 )
                                {
                                    prepGNUPlotScript << ", \\" << std::endl;
                                }
                                count--;
                                break;
                            default:
                                errorLog() << "No properties selected!";
                        }
359
                    }
360
                    // generate filename for plot script
361 362 363
                    std::stringstream prepPlotFilename;
                    prepPlotFilename << "/plot" << id + 1 << ".gsc";
                    std::string plotfilename( prepPlotFilename.str() );
364
                    // save plot script
365 366 367
                    std::ofstream plotFile( ( m_RaySaveFilePath->get( true ).string() + plotfilename ).c_str() );
                    plotFile << prepGNUPlotScript.str();
                    plotFile.close();
368
                }
369
            }
370
            // end progress
Franziska Jacob's avatar
Franziska Jacob committed
371 372
            saving->finish();
            m_progress->removeSubProgress( saving );
373
            // reactivate trigger
Franziska Jacob's avatar
Franziska Jacob committed
374 375 376
            m_saveTrigger->set( WPVBaseTypes::PV_TRIGGER_READY, false );
        }

377
        // check input connectors for new datasets
378
        bool dataChanged = ( m_inputData->getData() != m_dataSet );
379
        bool additionalData = ( m_inputFA->getData() != m_FAdataSet );
380
        bool newDerivative = ( m_inputDerivation->getData() != m_deriDataSet );
381 382 383
        bool dataValid  = static_cast< bool >( m_dataSet );
        m_DeriIsValid = static_cast< bool >( m_deriDataSet );
        m_FAisValid  = static_cast< bool >( m_FAdataSet );
384

385
        if( dataChanged ) // main data, compulsory
386
        {
387
            m_dataSet = m_inputData->getData();
388
            dataValid = static_cast< bool >( m_dataSet );
389
            debugLog() << "New Data!";
390

391 392 393 394 395 396 397 398
            if( !dataValid )
            {
                debugLog() << "No valid data anymore. Cleaning up.";
                m_rootNode->clear();
                m_grid.reset();
            }
            else
            {
Franziska Jacob's avatar
Franziska Jacob committed
399
                // grab the grid
400
                m_grid = boost::dynamic_pointer_cast< WGridRegular3D >( m_dataSet->getGrid() );
Franziska Jacob's avatar
Franziska Jacob committed
401 402 403 404 405 406 407
                if( !m_grid )
                {
                    errorLog() << "The dataset does not provide a regular grid. Ignoring dataset.";
                    dataValid = false;
                }
            }
        }
Franziska Jacob's avatar
Franziska Jacob committed
408

409
        if( newDerivative ) // optional: dataset with derivated data
410 411
        {
            m_deriDataSet = m_inputDerivation->getData();
412
            m_DeriIsValid = static_cast< bool >( m_deriDataSet );
413 414 415 416 417 418 419 420 421 422 423 424 425
            debugLog() << "New Derivated Data!";

            if( !m_DeriIsValid )
            {
                debugLog() << "No valid derivated data anymore. Cleaning up.";
                m_rootNode->clear();
                m_deriGrid.reset();
                m_curveMeanOut->reset();
                m_curveGaussOut->reset();
            }
            else
            {
                // grab the grid
426
                m_deriGrid = boost::dynamic_pointer_cast< WGridRegular3D >( m_deriDataSet->getGrid() );
427 428 429 430 431 432
                if( !m_deriGrid )
                {
                    errorLog() << "The derivated dataset does not provide a regular grid. Ignoring dataset.";
                    m_DeriIsValid = false;
                }
            }
433 434 435 436 437 438 439 440 441

            if( m_DeriIsValid )
            {
                // reset old calculations
                m_meanCurvDataSet.reset();
                m_gaussCurvDataSet.reset();
                // this may take some time, depending on the data
                calculateCurvature();
            }
442 443
        }

444
        if( additionalData ) // optional: fractional anisotropy data
445 446
        {
            m_FAdataSet = m_inputFA->getData();
447
            m_FAisValid = static_cast< bool >( m_FAdataSet );
448 449 450 451 452 453 454 455 456 457 458
            debugLog() << "New FA Data!";

            if( !m_FAisValid )
            {
                debugLog() << "No valid FA data anymore. Cleaning up.";
                m_rootNode->clear();
                m_FAgrid.reset();
            }
            else
            {
                // grab the grid
459
                m_FAgrid = boost::dynamic_pointer_cast< WGridRegular3D >( m_FAdataSet->getGrid() );
460 461 462 463 464 465 466 467
                if( !m_FAgrid )
                {
                    errorLog() << "The FA dataset does not provide a regular grid. Ignoring dataset.";
                    m_FAisValid = false;
                }
            }
        }

468
        bool propsChanged = m_interval->changed() || m_rayNumber->changed() || m_radius->changed();
Franziska Jacob's avatar
Franziska Jacob committed
469

470 471
        // observes properties and connectors
        // main calculation of the module starts in onClick function
472
        if( ( propsChanged || dataChanged ) && dataValid )
473
        {
474
            debugLog() << "Data or properties changed.";
475

476
            // grid scale
Franziska Jacob's avatar
Franziska Jacob committed
477 478 479 480
            double x_scale = m_grid->getNbCoordsX();
            double y_scale = m_grid->getNbCoordsY();
            double z_scale = m_grid->getNbCoordsZ();

481
            // voxel scale
482 483 484
//             double dist_x = m_grid->getOffsetX();
//             double dist_y = m_grid->getOffsetY();
//             double dist_z = m_grid->getOffsetZ();
485

Franziska Jacob's avatar
Franziska Jacob committed
486
            debugLog() << "x = " << x_scale << ", y = " << y_scale << ", z = " << z_scale;
Sebastian Eichelbaum's avatar
[MERGE]  
Sebastian Eichelbaum committed
487

488
            //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Franziska Jacob's avatar
Franziska Jacob committed
489
            // Calculating BoundingBox and OuterBounding
490
            //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Franziska Jacob's avatar
Franziska Jacob committed
491 492 493 494
            if( dataChanged )
            {
                std::vector< WVector4d > bounding_box( 8 );  // BoundingBox of the dataset

495 496 497 498 499 500 501 502
                bounding_box[0] = WVector4d( 0,       0,       0,       1 );
                bounding_box[1] = WVector4d( x_scale, 0,       0,       1 );
                bounding_box[2] = WVector4d( 0,       y_scale, 0,       1 );
                bounding_box[3] = WVector4d( x_scale, y_scale, 0,       1 );
                bounding_box[4] = WVector4d( 0,       0,       z_scale, 1 );
                bounding_box[5] = WVector4d( x_scale, 0,       z_scale, 1 );
                bounding_box[6] = WVector4d( 0,       y_scale, z_scale, 1 );
                bounding_box[7] = WVector4d( x_scale, y_scale, z_scale, 1 );
Sebastian Eichelbaum's avatar
[MERGE]  
Sebastian Eichelbaum committed
503

Franziska Jacob's avatar
Franziska Jacob committed
504 505 506 507
                // BoundingBox parallel to the screen (just min and max)
                m_outer_bounding.resize( 2 );
                m_outer_bounding[0] = bounding_box[0];
                m_outer_bounding[1] = bounding_box[0];
Sebastian Eichelbaum's avatar
[MERGE]  
Sebastian Eichelbaum committed
508

509
                for( unsigned int k = 0; k < 8; k++ )
Franziska Jacob's avatar
Franziska Jacob committed
510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525
                {
                    // minimal point in BoundingBox
                    m_outer_bounding[0][0] = std::min( bounding_box[k][0], m_outer_bounding[0][0] ); // x-min
                    m_outer_bounding[0][1] = std::min( bounding_box[k][1], m_outer_bounding[0][1] ); // y-min
                    m_outer_bounding[0][2] = std::min( bounding_box[k][2], m_outer_bounding[0][2] ); // z-min
                    m_outer_bounding[0][3] = 1;
    //                  debugLog() << "min-vec:" << m_outer_bounding[0];
                    // maximal point in BoundingBox
                    m_outer_bounding[1][0] = std::max( bounding_box[k][0], m_outer_bounding[1][0] ); // x-max
                    m_outer_bounding[1][1] = std::max( bounding_box[k][1], m_outer_bounding[1][1] ); // y-max
                    m_outer_bounding[1][2] = std::max( bounding_box[k][2], m_outer_bounding[1][2] ); // z-max
                    m_outer_bounding[1][3] = 1;
    //                  debugLog() << "max-vec:" << m_outer_bounding[1];
                }
            }
            //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
526
            // Getting Data and setting Properties for GUI
Franziska Jacob's avatar
Franziska Jacob committed
527
            //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
528 529 530
            m_interval->setMin( 0.0 );
            m_interval->setMax( 1.0 );
            //m_interval->setMax( length( WVector4d( dist_x, dist_y, dist_z, 0.0 ) ) );
531
            m_interval->setRecommendedValue( 0.30 );
532

533
            m_rayNumber->setMin( 1 );
534
            m_rayNumber->setRecommendedValue( 1 );
535

536
            m_radius->setMin( 1 );
537 538
            m_radius->setRecommendedValue( 2 );
            //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
539
            // Drawing current position of the future ray from properties (debugging)
540
            //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
541
            if( false ) //just for debugging
542 543 544 545
            {
                osg::ref_ptr< osg::Geode > rayGeode = new osg::Geode();
                rayGeode->addUpdateCallback( new SafeUpdateCallback( this ) );

546
                WVector4d cylStart = WVector4d( 0.5 * x_scale, 0.5 * y_scale, 0.0, 1.0 );
547

548 549
                osg::Cylinder* cylinder = new osg::Cylinder( getAs3D( cylStart + WVector4d( 0.0, 0.0, 0.5 * z_scale, 0.0 ) ),
                                                                0.5f, static_cast<float>( z_scale ) );
550 551 552 553 554
                osg::Quat rot;
                rot.makeRotate( osg::Vec3d( 0.0, 0.0, 1.0 ),
                                osg::Vec3d( 1.0, 0.0, 0.0 )   // <-- ray direction vector
                            );
                cylinder->setRotation( rot );
555

556 557
                osg::Cone* cone = new osg::Cone( getAs3D( cylStart + WVector4d( 0.0, 0.0, z_scale, 0.0 ) ), 1.0f, 5.0f );
                cone->setRotation( rot );
558

559 560
                rayGeode->addDrawable( new osg::ShapeDrawable( cylinder ) );
                rayGeode->addDrawable( new osg::ShapeDrawable( cone ) );
561

562 563 564
                m_rootNode->clear();
                m_rootNode->insert( rayGeode );
            }
565 566 567 568 569 570 571 572 573
        }
    }
    // we technically need to remove the event handler too but there is no method available to do this ...
    // WKernel::getRunningKernel()->getGraphicsEngine()->getViewer()->removeEventHandler( handler );
    WKernel::getRunningKernel()->getGraphicsEngine()->getScene()->remove( m_rootNode );
}

void WMTransferCalc::onClick( WVector2i mousePos )
{
574
/*    debugLog() << "Left Click at " << mousePos;
575 576 577 578
    debugLog() << "Projection Matrix: " << m_matrixCallback->getProjectionMatrix()->get();
    debugLog() << "ModelView Matrix: " << m_matrixCallback->getModelViewMatrix()->get();
    debugLog() << "Viewport: " << m_matrixCallback->getViewportX()->get() << ", " << m_matrixCallback->getViewportY()->get() << ", "
                               << m_matrixCallback->getViewportWidth()->get() << ", " << m_matrixCallback->getViewportHeight()->get();
579
*/
580
    // get both matrices inverted
581 582
    WMatrix4d projectionMatrixInverted = invert( m_matrixCallback->getProjectionMatrix()->get() );
    WMatrix4d modelviewMatrixInverted = invert( m_matrixCallback->getModelViewMatrix()->get() );
583

584 585 586
    // get the clip-space coordinate:
    WVector4d pInClipSpace( ( 2.0 * mousePos.x() / m_matrixCallback->getViewportWidth()->get() ) - 1.0,
                            ( 2.0 * mousePos.y() / m_matrixCallback->getViewportHeight()->get() ) - 1.0,
587 588 589 590
                              0.0,
                              1.0 );
    WVector4d dirInClipSpace( 0.0,
                              0.0,
591 592 593
                              1.0,
                              0.0 );

594 595 596
    // unproject the clip-space vectors to the world space
    WVector4d pInWorldSpace   = projectionMatrixInverted * pInClipSpace;
    WVector4d dirInWorldSpace = projectionMatrixInverted * dirInClipSpace;
597
    // get back to model-space
598 599
    WVector4d pInObjectSpace   = modelviewMatrixInverted * pInWorldSpace;
    WVector4d dirInObjectSpace = modelviewMatrixInverted * dirInWorldSpace;
600

601
    // for multilication with a vertex
602
    //WMatrix4d totalTransform = modelviewMatrixInverted * projectionMatrixInverted;
603

604 605
//     debugLog() << pInWorldSpace << " --- " << pInObjectSpace;
//     debugLog() << dirInWorldSpace << " --- " << dirInObjectSpace;
606

607
    // debug "vector"
608
    osg::ref_ptr< osg::Geode > debugGeode = new osg::Geode();
609 610 611
    debugGeode->addUpdateCallback( new SafeUpdateCallback( this ) );
    debugGeode->addDrawable( new osg::ShapeDrawable( new osg::Sphere( getAs3D( pInObjectSpace, true ), 2.5f ) ) );
    debugGeode->addDrawable( new osg::ShapeDrawable( new osg::Sphere( getAs3D( pInObjectSpace + dirInObjectSpace, true ), 5.0f ) ) );
612
    m_rootNode->clear();
613
    m_rootNode->insert( debugGeode );
614

615 616 617
    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    // Ray Casting
    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
618

619 620
    // ray object calculated from click position: ray = start + t * direction
    WRay ray( pInObjectSpace, dirInObjectSpace );
621
    // delete old profiles
622 623
    m_profiles.clear();

624
    // get user defined properties
625 626
    double interval = m_interval->get( true );
    unsigned int samplesInVicinity = static_cast<unsigned int>( m_rayNumber->get( true ) );
627
    unsigned int radius = static_cast<unsigned int>( m_radius->get( true ) );
628

629 630 631
    osg::ref_ptr< osg::Geode > rayGeode = new osg::Geode();
    rayGeode->addUpdateCallback( new SafeUpdateCallback( this ) );

632
    // seed for random
633
    std::srand( time( 0 ) );
634
    // cast additional rays if user defined more than one ray in properties
635 636 637 638 639 640 641 642 643
    for( unsigned int n = 0; n < samplesInVicinity; n++ )
    {
        if( n == 0 ) // the initial profile
        {
            m_mainProfile = castRay( ray, interval, rayGeode );
            m_profiles.push_back( m_mainProfile );
        }
        else // random profiles in vicinity
        {
644 645 646 647 648
            // calculate random angle in [0,2pi]
            double angle = 2 * M_PI * rand() / ( static_cast< double >( RAND_MAX ) );
            // calculate random radius in [0,1]
            double random_rad = rand() / ( static_cast< double >( RAND_MAX ) );

649
            // create orthogonal plane to ray for calculating rays in vicinity
650 651
            // (the ray gets transformed and scaled before this step,
            // so we have to calculate this plane to get a uniform scale for the radius)
652
            double coord = std::max( ray.direction().x(), std::max( ray.direction().y(), ray.direction().z() ) );
653 654 655 656 657
            WVector3d helper( coord, coord, coord );
            WVector3d ax1 = cross( getAs3D( ray.direction(), true ), helper );
            WVector3d ax2 = cross( getAs3D( ray.direction(), true ), ax1 );
            ax1 = normalize( ax1 );
            ax2 = normalize( ax2 );
658
            // base vectors of the calculated plane
659 660 661
            WVector4d ax1_4( ax1.x(), ax1.y(), ax1.z(), 0.0 );
            WVector4d ax2_4( ax2.x(), ax2.y(), ax2.z(), 0.0 );

662
            // calculate random ray
663 664
            WVector4d vicPoint = ray.start() +
                                ( ax1_4 * sqrt( random_rad ) * std::cos( angle ) * radius ) +
665
                                ( ax2_4 * sqrt( random_rad ) * std::sin( angle ) * radius );
666
            WRay vicinity( vicPoint, ray.direction() );
667 668 669 670 671
            m_profiles.push_back( castRay( vicinity, interval, rayGeode ) );
        }
    }
    m_rootNode->clear();
    m_rootNode->insert( rayGeode );
Franziska Jacob's avatar
Franziska Jacob committed
672
}
673

674
WRayProfile WMTransferCalc::castRay( WRay ray, double interval, osg::ref_ptr< osg::Geode > rayGeode )
Franziska Jacob's avatar
Franziska Jacob committed
675 676
{
    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
677
    // Calculating a single profile with all samples
Franziska Jacob's avatar
Franziska Jacob committed
678
    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
679
    // start progress for GUI
680 681
    boost::shared_ptr< WProgress > prog = boost::shared_ptr< WProgress >( new WProgress( "Casting ray." ) );
    m_progress->addSubProgress( prog );
682

683
    //TODO(fjacob): the not set RaySamples may be deleted after collecting all profiles
Franziska Jacob's avatar
Franziska Jacob committed
684 685 686
    size_t max_nbSamples = ceil( length( m_outer_bounding[1] - m_outer_bounding[0] ) / interval );
    WRayProfile curProfile( max_nbSamples );
//   debugLog() << "Max samples: " << max_nbSamples;
687

688
    // BoxIntersecParameter contains two values, the first and second intersection parameter with the grid
689
    struct BoxIntersecParameter bounds = rayIntersectsBox( ray );
690
    if( bounds.isNull ) //no intersection with data grid
691 692 693 694 695
    {
        prog->finish();
        m_progress->removeSubProgress( prog );
        return curProfile; //empty profile
    }
696 697
    double start_t = bounds.minimum_t;
    double end_t   = bounds.maximum_t;
698

699 700
//    debugLog() << "Start:" << start_t;
//    debugLog() << "End  :" << end_t;
701

702
    // helpers for visualization
703
    WVector3d geodeStartVec, geodeEndVec, cylinderVec;
Sebastian Eichelbaum's avatar
[MERGE]  
Sebastian Eichelbaum committed
704

705 706
    size_t sampleCount = 0;
    for( double sample_t = start_t; sample_t <= end_t; sample_t += interval )
Franziska Jacob's avatar
Franziska Jacob committed
707 708
    {
        WVector4d current = ray.getSpot( sample_t );
709
        WVector3d cur3D = getAs3D( current );
710

Franziska Jacob's avatar
Franziska Jacob committed
711
        // do not calculate anything for vectors outside of the data grid
712
        // (should not be neccessary anymore, though)
713
        if( !m_grid->encloses( cur3D ) )
Franziska Jacob's avatar
Franziska Jacob committed
714 715 716
        {
            //debugLog() << "continue...";
            continue;
717
        }
718
        // save position of the sample point in profile
719 720 721 722 723
        double vox_x = static_cast< double >( m_grid->getVoxelNum( cur3D ) );
        double vox_y = static_cast< double >( m_grid->getVoxelNum( cur3D ) );
        double vox_z = static_cast< double >( m_grid->getVoxelNum( cur3D ) );
        curProfile[sampleCount].position() = WVector4d( vox_x, vox_y, vox_z, 1 );

724
        // set start vector for drawing the ray
725
        if( sampleCount == 0 )
726
        {
727
            geodeStartVec = cur3D;
728 729
        }

730
        // get interpolated value at current sample point
731
        double val = interpolate( current, m_grid, m_dataSet );
732
        // save value
Franziska Jacob's avatar
Franziska Jacob committed
733
        curProfile[sampleCount].value() = val;
734
        // save distance t
735
        curProfile[sampleCount].distance() = sample_t;
736

737
        // if there is any FA data available, save it
738 739
        if( m_FAisValid )
        {
740
            if( m_FAgrid->encloses( cur3D ) )
741
            {
742
                // save fractional anisotropy
743 744 745
                curProfile[sampleCount].fracA() = interpolate( current, m_FAgrid, m_FAdataSet );
            }
        }
746
        // if there is any derivated data available, save calculated curvature
747 748 749 750 751 752 753 754 755 756 757 758
        if( m_DeriIsValid )
        {
            if( m_deriGrid->encloses( cur3D ) )
            {
                if( m_meanCurvDataSet )
                {
                    curProfile[sampleCount].meanCurv() = interpolate( current, m_deriGrid, m_meanCurvDataSet );
                }
                if( m_gaussCurvDataSet )
                {
                    curProfile[sampleCount].gaussCurv() = interpolate( current, m_deriGrid, m_gaussCurvDataSet );
                }
759 760 761
            }
        }

762
        // calculate gradient
763
        WVector4d grad = getGradient( current );
764
        double weight = length( grad );
765
        // save gradient weight
766
        curProfile[sampleCount].gradWeight() = weight;
767 768
        // normalize
        if( weight != 0 )
769
        {
770
            // save gradient
771 772 773
            curProfile[sampleCount].gradient() = grad * ( 1 / weight );
        }

774 775 776 777 778
        // calculate angle between gradient and ray
        // gradVec and rayVec are already normalized
        Eigen::Vector3d gradVec( curProfile[sampleCount].gradient()[0], curProfile[sampleCount].gradient()[1],
                                    curProfile[sampleCount].gradient()[2] );
        Eigen::Vector3d rayVec( ray.direction()[0], ray.direction()[1], ray.direction()[2] );
779 780 781
//         curProfile[sampleCount].angle() = std::min( std::acos( gradVec.dot( rayVec ) ),
//                                                     std::acos( gradVec.dot( rayVec * -1 ) ) );
        curProfile[sampleCount].angle() = std::acos( gradVec.dot( rayVec * -1 ) );
782

783
        // set helper for visualization
784
        geodeEndVec = cur3D;
Franziska Jacob's avatar
Franziska Jacob committed
785
        sampleCount++;
786
    }
787
    // draw ray (a cylinder and a cone)
788 789
    // (it is a little tricky because osg does not support
    // native rotation of its drawables on initialization)
790
    cylinderVec = geodeEndVec - geodeStartVec;
791

792
    // set cylinder center in the middle of the ray (will be rotation center)
793 794
    osg::Cylinder* cylinder = new osg::Cylinder( geodeStartVec + ( 0.5 * cylinderVec ), 0.5f, length( cylinderVec ) );
    osg::Cone* cone = new osg::Cone( geodeEndVec, 1.0f, 5.0f );
795

796 797 798
    // let osg calculate the neccessary rotation matrix and rotate the shapes
    osg::Quat rotation;
    rotation.makeRotate( osg::Vec3d( 0.0, 0.0, 1.0 ),
799
                    osg::Vec3d( ray.direction().x(), ray.direction().y(), ray.direction().z() )
800
                );
801 802
    cylinder->setRotation( rotation );
    cone->setRotation( rotation );
803

804
    // draw it
805 806 807
    rayGeode->addDrawable( new osg::ShapeDrawable( cylinder ) );
    rayGeode->addDrawable( new osg::ShapeDrawable( cone ) );

808
    // end progress
Franziska Jacob's avatar
Franziska Jacob committed
809 810
    prog->finish();
    m_progress->removeSubProgress( prog );
Franziska Jacob's avatar
Franziska Jacob committed
811

Franziska Jacob's avatar
Franziska Jacob committed
812 813
    return curProfile;
}
814

815
struct BoxIntersecParameter WMTransferCalc::rayIntersectsBox( WRay ray )
Franziska Jacob's avatar
Franziska Jacob committed
816
{
817
    // ray = start + t * dir
818 819
    double tnear = - std::numeric_limits<double>::max();
    double tfar  =   std::numeric_limits<double>::max();
Franziska Jacob's avatar
Franziska Jacob committed
820
    double t0, t1;
821

822 823 824
    struct BoxIntersecParameter result;
    result.isNull = true;

Franziska Jacob's avatar
Franziska Jacob committed
825 826
    for( unsigned int coord = 0; coord < 3; coord++ ) // for x,y,z
    {
827
        if( ray.direction()[coord] == 0 ) // ray parallel to the min/max coordinate planes
Franziska Jacob's avatar
Franziska Jacob committed
828
        {
829
            if( ray.start()[coord] < m_outer_bounding[0][coord] || ray.start()[coord] > m_outer_bounding[1][coord] )
830
            {
831
                return result; // ray not within min and max values of current coordinate
832
            }
Franziska Jacob's avatar
Franziska Jacob committed
833 834 835 836
        }
        else // ray not parallel, so calculate intersections
        {
            // time at which ray intersects min plane
837
            t0 = ( m_outer_bounding[0][coord] - ray.start()[coord] ) / ray.direction()[coord];
Franziska Jacob's avatar
Franziska Jacob committed
838
            // time at which ray intersects max plane
Franziska Jacob's avatar