WMTransferCalc.cpp

//---------------------------------------------------------------------------
//
// 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
#define _USE_MATH_DEFINES

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <fstream>
#include <iostream>
#include <limits>
#include <sstream>
#include <string>
#include <vector>

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

#include <Eigen/Eigen>


#include "core/kernel/WKernel.h"
#include "core/common/WColor.h"
#include "core/common/WPathHelper.h"
#include "core/common/WPropertyHelper.h"
#include "core/common/math/WMatrix.h"
#include "core/common/math/linearAlgebra/WMatrixFixed.h"
#include "core/common/math/linearAlgebra/WVectorFixed.h"
#include "core/common/WItemSelection.h"
#include "core/common/WItemSelectionItem.h"
#include "core/common/WItemSelectionItemTyped.h"
#include "core/common/WItemSelector.h"
#include "core/graphicsEngine/WGEUtils.h"
#include "core/graphicsEngine/WGERequirement.h"

#include "../interaction/WViewEventHandler.h"
#include "../interaction/WGetMatrixCallback.h"

#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
{
    return "Transfer Function Calculation";
}

const std::string WMTransferCalc::getDescription() const
{
    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.";
}

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

    // 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." );

    // call WModules initialization
    WModule::connectors();
}

void WMTransferCalc::properties()
{
    m_propCondition   = boost::shared_ptr< WCondition >( new WCondition() );
    m_color           = m_properties->addProperty(  "Color", "Color of the ray.", WColor( 1.0, 0.0, 0.0, 1.0 ), m_propCondition );

    m_interval        = m_properties->addProperty( "Interval", "Interval for sampling rays.", 0.30, m_propCondition );
    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 );

    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." ) );

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


    m_RaySaveFilePath = m_properties->addProperty( "Save RayProfile to Folder:",
                                                   "Savefile for RayProfile.", WPathHelper::getAppPath(), m_propCondition );
    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 );
    WPropertyHelper::PC_NOTEMPTY::addTo( m_multiPlotDataSelection );

    WModule::properties();
}

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

bool WMTransferCalc::saveRayProfileTo( WPropFilename path, std::string filename, WRayProfile *profile )
{
    WRayProfile *prof;
    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;

    if( prof->size() == 0 )
    {
        return false; //  empty profiles will not be saved
    }

    // compose context of the soon-to-be-saved file which will be saved in given path with given filename
    std::ofstream savefile( ( path->get( true ).string() + "/" + filename ).c_str() /*, ios::app */ );
    // heading
    savefile << "# RayProfile data for OpenWalnut TransferCalcModule" << std::endl;
    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;
    // save data of each sample point
    for( size_t sample = 0; sample < prof->size(); sample++ )
    {
        savefile << sample + 1 << "\t";                         // print id
        savefile << (*prof)[sample].distance() << "\t";         // print t
        savefile << (*prof)[sample].value() << "\t";            // print value
        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
        savefile << (*prof)[sample].fracA() << "\t";            // print fractional anisotropy
        savefile << (*prof)[sample].angle() << "\t";            // print angle between gradient and ray
        savefile << (*prof)[sample].meanCurv() << "\t";         // print mean curvature
        savefile << (*prof)[sample].gaussCurv() << "\t";        // print gauss curvature
        savefile << std::endl;
    }
    savefile.close();
    return true;
}

void WMTransferCalc::moduleMain()
{
    debugLog() << "moduleMain started...";

    // init observation of properties
    m_moduleState.setResetable( true, true );
    m_moduleState.add( m_inputData->getDataChangedCondition() );
    m_moduleState.add( m_inputDerivation->getDataChangedCondition() );
    m_moduleState.add( m_inputFA->getDataChangedCondition() );
    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 );

    // 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 );

    debugLog() << "Entering main loop";
    while( !m_shutdownFlag() )
    {
        debugLog() << "Waiting ...";
        m_moduleState.wait();
        debugLog() << "New Data Incoming ...";

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

        // activates if user pushes save button to save current RayProfile
        if( m_saveTrigger->get( true ) == WPVBaseTypes::PV_TRIGGER_TRIGGERED )
        {
            debugLog() << "User saves RayProfile(s).";

            // start progress for GUI
            boost::shared_ptr< WProgress > saving = boost::shared_ptr< WProgress >( new WProgress( "Saving calculated RayProfiles." ) );
            m_progress->addSubProgress( saving );

            // save data and plot file for every profile
            unsigned int nrProfiles = m_profiles.size();
            for( unsigned int id = 0; id < nrProfiles; id++ )
            {
                // generate header
                // first a general header for setting default options
                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;

                // generate filename of data file
                std::stringstream prepFilename;
                prepFilename << "RayProfile" << id + 1 << ".dat";
                std::string filename( prepFilename.str() );

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

                        // switch works with an enumeration
                        switch( s.at( i )->getAs< MyItemType >()->getValue() )
                        {
                            case VALUE:
                                prepGNUPlotScript << "'" << filename << "' using 2:3 title 'value' with linespoints";
                                if( count - 1 != 0 )
                                {
                                    prepGNUPlotScript << ", \\" << std::endl;
                                }
                                count--;
                                break;
                            case WEIGHT:
                                prepGNUPlotScript << "'" << filename << "' using 2:7 title 'gradient weight' with linespoints";
                                if( count - 1 != 0 )
                                {
                                    prepGNUPlotScript << ", \\" << std::endl;
                                }
                                count--;
                                break;
                            case FA:
                                prepGNUPlotScript << "'" << filename << "' using 2:($8*100) title 'FA value' with linespoints";
                                if( count - 1 != 0 )
                                {
                                    prepGNUPlotScript << ", \\" << std::endl;
                                }
                                count--;
                                break;
                            case ANGLE:
                                prepGNUPlotScript << "'" << filename << "' using 2:($9*(180/3.14159265)) title 'angle' with linespoints";
                                if( count - 1 != 0 )
                                {
                                    prepGNUPlotScript << ", \\" << std::endl;
                                }
                                count--;
                                break;
                            case MEAN:
                                prepGNUPlotScript << "'" << filename << "' using 2:($10*100) title 'mean curvature' with linespoints";
                                if( count - 1 != 0 )
                                {
                                    prepGNUPlotScript << ", \\" << std::endl;
                                }
                                count--;
                                break;
                            case GAUSS:
                                prepGNUPlotScript << "'" << filename << "' using 2:($11*100) title 'gaussian curvature' with linespoints";
                                if( count - 1 != 0 )
                                {
                                    prepGNUPlotScript << ", \\" << std::endl;
                                }
                                count--;
                                break;
                            default:
                                errorLog() << "No properties selected!";
                        }
                    }
                    // generate filename for plot script
                    std::stringstream prepPlotFilename;
                    prepPlotFilename << "/plot" << id + 1 << ".gsc";
                    std::string plotfilename( prepPlotFilename.str() );
                    // save plot script
                    std::ofstream plotFile( ( m_RaySaveFilePath->get( true ).string() + plotfilename ).c_str() );
                    plotFile << prepGNUPlotScript.str();
                    plotFile.close();
                }
            }
            // end progress
            saving->finish();
            m_progress->removeSubProgress( saving );
            // reactivate trigger
            m_saveTrigger->set( WPVBaseTypes::PV_TRIGGER_READY, false );
        }

        // check input connectors for new datasets
        bool dataChanged = ( m_inputData->getData() != m_dataSet );
        bool additionalData = ( m_inputFA->getData() != m_FAdataSet );
        bool newDerivative = ( m_inputDerivation->getData() != m_deriDataSet );
        bool dataValid  = static_cast< bool >( m_dataSet );
        m_DeriIsValid = static_cast< bool >( m_deriDataSet );
        m_FAisValid  = static_cast< bool >( m_FAdataSet );

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

            if( !dataValid )
            {
                debugLog() << "No valid data anymore. Cleaning up.";
                m_rootNode->clear();
                m_grid.reset();
            }
            else
            {
                // grab the grid
                m_grid = boost::dynamic_pointer_cast< WGridRegular3D >( m_dataSet->getGrid() );
                if( !m_grid )
                {
                    errorLog() << "The dataset does not provide a regular grid. Ignoring dataset.";
                    dataValid = false;
                }
            }
        }

        if( newDerivative ) // optional: dataset with derivated data
        {
            m_deriDataSet = m_inputDerivation->getData();
            m_DeriIsValid = static_cast< bool >( m_deriDataSet );
            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
                m_deriGrid = boost::dynamic_pointer_cast< WGridRegular3D >( m_deriDataSet->getGrid() );
                if( !m_deriGrid )
                {
                    errorLog() << "The derivated dataset does not provide a regular grid. Ignoring dataset.";
                    m_DeriIsValid = false;
                }
            }

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

        if( additionalData ) // optional: fractional anisotropy data
        {
            m_FAdataSet = m_inputFA->getData();
            m_FAisValid = static_cast< bool >( m_FAdataSet );
            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
                m_FAgrid = boost::dynamic_pointer_cast< WGridRegular3D >( m_FAdataSet->getGrid() );
                if( !m_FAgrid )
                {
                    errorLog() << "The FA dataset does not provide a regular grid. Ignoring dataset.";
                    m_FAisValid = false;
                }
            }
        }

        bool propsChanged = m_interval->changed() || m_rayNumber->changed() || m_radius->changed();

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

            // grid scale
            double x_scale = m_grid->getNbCoordsX();
            double y_scale = m_grid->getNbCoordsY();
            double z_scale = m_grid->getNbCoordsZ();

            // voxel scale
//             double dist_x = m_grid->getOffsetX();
//             double dist_y = m_grid->getOffsetY();
//             double dist_z = m_grid->getOffsetZ();

            debugLog() << "x = " << x_scale << ", y = " << y_scale << ", z = " << z_scale;

            //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            // Calculating BoundingBox and OuterBounding
            //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            if( dataChanged )
            {
                std::vector< WVector4d > bounding_box( 8 );  // BoundingBox of the dataset

                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 );

                // 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];

                for( unsigned int k = 0; k < 8; k++ )
                {
                    // 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];
                }
            }
            //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            // Getting Data and setting Properties for GUI
            //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            m_interval->setMin( 0.0 );
            m_interval->setMax( 1.0 );
            //m_interval->setMax( length( WVector4d( dist_x, dist_y, dist_z, 0.0 ) ) );
            m_interval->setRecommendedValue( 0.30 );

            m_rayNumber->setMin( 1 );
            m_rayNumber->setRecommendedValue( 1 );

            m_radius->setMin( 1 );
            m_radius->setRecommendedValue( 2 );
            //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            // Drawing current position of the future ray from properties (debugging)
            //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            if( false ) //just for debugging
            {
                osg::ref_ptr< osg::Geode > rayGeode = new osg::Geode();
                rayGeode->addUpdateCallback( new SafeUpdateCallback( this ) );

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

                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 ) );
                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 );

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

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

                m_rootNode->clear();
                m_rootNode->insert( rayGeode );
            }
        }
    }
    // 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 )
{
/*    debugLog() << "Left Click at " << mousePos;
    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();
*/
    // get both matrices inverted
    WMatrix4d projectionMatrixInverted = invert( m_matrixCallback->getProjectionMatrix()->get() );
    WMatrix4d modelviewMatrixInverted = invert( m_matrixCallback->getModelViewMatrix()->get() );

    // 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,
                              0.0,
                              1.0 );
    WVector4d dirInClipSpace( 0.0,
                              0.0,
                              1.0,
                              0.0 );

    // unproject the clip-space vectors to the world space
    WVector4d pInWorldSpace   = projectionMatrixInverted * pInClipSpace;
    WVector4d dirInWorldSpace = projectionMatrixInverted * dirInClipSpace;
    // get back to model-space
    WVector4d pInObjectSpace   = modelviewMatrixInverted * pInWorldSpace;
    WVector4d dirInObjectSpace = modelviewMatrixInverted * dirInWorldSpace;

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

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

    // debug "vector"
    osg::ref_ptr< osg::Geode > debugGeode = new osg::Geode();
    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 ) ) );
    m_rootNode->clear();
    m_rootNode->insert( debugGeode );

    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    // Ray Casting
    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

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

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

    // seed for random
    std::srand( time( 0 ) );
    // cast additional rays if user defined more than one ray in properties
    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
        {
            // 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 ) );

            // create orthogonal plane to ray for calculating rays in vicinity
            // (the ray gets transformed and scaled before this step,
            // so we have to calculate this plane to get a uniform scale for the radius)
            double coord = std::max( ray.direction().x(), std::max( ray.direction().y(), ray.direction().z() ) );
            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 );
            // base vectors of the calculated plane
            WVector4d ax1_4( ax1.x(), ax1.y(), ax1.z(), 0.0 );
            WVector4d ax2_4( ax2.x(), ax2.y(), ax2.z(), 0.0 );

            // calculate random ray
            WVector4d vicPoint = ray.start() +
                                ( ax1_4 * sqrt( random_rad ) * std::cos( angle ) * radius ) +
                                ( ax2_4 * sqrt( random_rad ) * std::sin( angle ) * radius );
            WRay vicinity( vicPoint, ray.direction() );
            m_profiles.push_back( castRay( vicinity, interval, rayGeode ) );
        }
    }
    m_rootNode->clear();
    m_rootNode->insert( rayGeode );
}

WRayProfile WMTransferCalc::castRay( WRay ray, double interval, osg::ref_ptr< osg::Geode > rayGeode )
{
    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    // Calculating a single profile with all samples
    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    // start progress for GUI
    boost::shared_ptr< WProgress > prog = boost::shared_ptr< WProgress >( new WProgress( "Casting ray." ) );
    m_progress->addSubProgress( prog );

    //TODO(fjacob): the not set RaySamples may be deleted after collecting all profiles
    size_t max_nbSamples = ceil( length( m_outer_bounding[1] - m_outer_bounding[0] ) / interval );
    WRayProfile curProfile( max_nbSamples );
//   debugLog() << "Max samples: " << max_nbSamples;

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

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

    // helpers for visualization
    WVector3d geodeStartVec, geodeEndVec, cylinderVec;

    size_t sampleCount = 0;
    for( double sample_t = start_t; sample_t <= end_t; sample_t += interval )
    {
        WVector4d current = ray.getSpot( sample_t );
        WVector3d cur3D = getAs3D( current );

        // do not calculate anything for vectors outside of the data grid
        // (should not be neccessary anymore, though)
        if( !m_grid->encloses( cur3D ) )
        {
            //debugLog() << "continue...";
            continue;
        }
        // save position of the sample point in profile
        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 );

        // set start vector for drawing the ray
        if( sampleCount == 0 )
        {
            geodeStartVec = cur3D;
        }

        // get interpolated value at current sample point
        double val = interpolate( current, m_grid, m_dataSet );
        // save value
        curProfile[sampleCount].value() = val;
        // save distance t
        curProfile[sampleCount].distance() = sample_t;

        // if there is any FA data available, save it
        if( m_FAisValid )
        {
            if( m_FAgrid->encloses( cur3D ) )
            {
                // save fractional anisotropy
                curProfile[sampleCount].fracA() = interpolate( current, m_FAgrid, m_FAdataSet );
            }
        }
        // if there is any derivated data available, save calculated curvature
        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 );
                }
            }
        }

        // calculate gradient
        WVector4d grad = getGradient( current );
        double weight = length( grad );
        // save gradient weight
        curProfile[sampleCount].gradWeight() = weight;
        // normalize
        if( weight != 0 )
        {
            // save gradient
            curProfile[sampleCount].gradient() = grad * ( 1 / weight );
        }

        // 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] );
//         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 ) );

        // set helper for visualization
        geodeEndVec = cur3D;
        sampleCount++;
    }
    // draw ray (a cylinder and a cone)
    // (it is a little tricky because osg does not support
    // native rotation of its drawables on initialization)
    cylinderVec = geodeEndVec - geodeStartVec;

    // set cylinder center in the middle of the ray (will be rotation center)
    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 );

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

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

    // end progress
    prog->finish();
    m_progress->removeSubProgress( prog );

    return curProfile;
}

struct BoxIntersecParameter WMTransferCalc::rayIntersectsBox( WRay ray )
{
    // ray = start + t * dir
    double tnear = - std::numeric_limits<double>::max();
    double tfar  =   std::numeric_limits<double>::max();
    double t0, t1;

    struct BoxIntersecParameter result;
    result.isNull = true;

    for( unsigned int coord = 0; coord < 3; coord++ ) // for x,y,z
    {
        if( ray.direction()[coord] == 0 ) // ray parallel to the min/max coordinate planes
        {
            if( ray.start()[coord] < m_outer_bounding[0][coord] || ray.start()[coord] > m_outer_bounding[1][coord] )
            {
                return result; // ray not within min and max values of current coordinate
            }
        }
        else // ray not parallel, so calculate intersections
        {
            // time at which ray intersects min plane
            t0 = ( m_outer_bounding[0][coord] - ray.start()[coord] ) / ray.direction()[coord];
            // time at which ray intersects max plane