Commit 7b46ea43 authored by Mathias Goldau's avatar Mathias Goldau
Browse files

[MERGE]

parents 7b5a2898 a7b913b7
......@@ -7,10 +7,7 @@ the currect working directory to be project's main directory,
the following steps should do the trick:
> cd build
> ccmake ../src
PRESS 'c' for configuration
PRESS 'c' again for further configuration
PRESS 'g' to generate the build system including the Makefiles.
> cmake ../src
> make
Now the subdirectory "bin" of the "build" directory should contain an executable named walnut.
......
//---------------------------------------------------------------------------
//
// 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/>.
//
//---------------------------------------------------------------------------
#include <vector>
#include <map>
#include "WGeometryFunctions.h"
void wmath::tesselateIcosahedron( std::vector< WVector3D >* vertices, std::vector< unsigned int >* triangles, unsigned int level )
{
WAssert( vertices, "Missing input vector." );
WAssert( triangles, "Missing input vector." );
vertices->clear();
triangles->clear();
unsigned int nv = 12;
unsigned int nt = 20;
for( unsigned int t = 1; t <= level; ++t )
{
nv += 3 * nt / 2;
nt *= 4;
}
vertices->reserve( nv );
triangles->reserve( nt );
// add icosahedron vertices
double phi = 0.5 * ( 1.0 + sqrt( 5.0 ) );
std::vector< WVector3D > g;
vertices->push_back( WVector3D( phi, 1, 0 ) ); // 0
vertices->push_back( WVector3D( -phi, 1, 0 ) ); // 1
vertices->push_back( WVector3D( phi, -1, 0 ) ); // 2
vertices->push_back( WVector3D( -phi, -1, 0 ) ); // 3
vertices->push_back( WVector3D( 1, 0, phi ) ); // 4
vertices->push_back( WVector3D( -1, 0, phi ) ); // 5
vertices->push_back( WVector3D( 1, 0, -phi ) ); // 6
vertices->push_back( WVector3D( -1, 0, -phi ) ); // 7
vertices->push_back( WVector3D( 0, phi, 1 ) ); // 8
vertices->push_back( WVector3D( 0, -phi, 1 ) ); // 9
vertices->push_back( WVector3D( 0, phi, -1 ) ); // 10
vertices->push_back( WVector3D( 0, -phi, -1 ) ); // 11
for( std::vector< WVector3D >::iterator it = vertices->begin(); it != vertices->end(); ++it )
{
*it = it->normalized();
}
// add triangle indices
unsigned int inc[ 60 ] =
{
8, 5, 4,
8, 4, 0,
8, 0, 10,
8, 10, 1,
8, 1, 5,
5, 9, 4,
4, 2, 0,
0, 6, 10,
10, 7, 1,
1, 3, 5,
4, 9, 2,
0, 2, 6,
10, 6, 7,
1, 7, 3,
5, 3, 9,
9, 11, 2,
2, 11, 6,
6, 11, 7,
7, 11, 3,
3, 11, 9
};
triangles->assign( inc, inc + 60 );
std::map< utility::Edge, unsigned int > edgeVertices;
for( unsigned int t = 0; t < level; ++t )
{
// for every triangle
std::size_t numTriangles = triangles->size() / 3;
for( std::size_t k = 0; k < numTriangles; ++k )
{
unsigned int idx[ 3 ];
// generate a new vertex for every edge (if there is no vertex yet)
for( int i = 0; i < 3; ++i )
{
utility::Edge e( ( *triangles )[ 3 * k + i ], ( *triangles )[ 3 * k + ( i + 1 ) % 3 ] );
if( edgeVertices.find( e ) == edgeVertices.end() )
{
WVector3D v0 = vertices->at( e.first );
WVector3D v1 = vertices->at( e.second );
WVector3D v = v0 + v1;
v = v.normalized();
vertices->push_back( v );
edgeVertices[ e ] = vertices->size() - 1;
}
idx[ i ] = edgeVertices[ e ];
}
// make 4 triangles from the current one
// add 1st triangle
triangles->push_back( ( *triangles )[ 3 * k + 0 ] );
triangles->push_back( idx[ 0 ] );
triangles->push_back( idx[ 2 ] );
// add 2nd triangle
triangles->push_back( ( *triangles )[ 3 * k + 1 ] );
triangles->push_back( idx[ 1 ] );
triangles->push_back( idx[ 0 ] );
// add 3rd triangle
triangles->push_back( ( *triangles )[ 3 * k + 2 ] );
triangles->push_back( idx[ 2 ] );
triangles->push_back( idx[ 1 ] );
// changed indices of this triangle to the indices of the 4th triangle
( *triangles )[ 3 * k + 0 ] = idx[ 0 ];
( *triangles )[ 3 * k + 1 ] = idx[ 1 ];
( *triangles )[ 3 * k + 2 ] = idx[ 2 ];
}
edgeVertices.clear();
}
}
//---------------------------------------------------------------------------
//
// 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/>.
//
//---------------------------------------------------------------------------
#ifndef WGEOMETRYFUNCTIONS_H
#define WGEOMETRYFUNCTIONS_H
#include <utility>
#include <map>
#include <vector>
#include "WVector3D.h"
#include "../WAssert.h"
namespace wmath
{
namespace utility
{
/**
* \class Edge
*
* A helper class that is used to store edges as pairs of vertex indices. The indices are sorted.
*/
class Edge : public std::pair< unsigned int, unsigned int >
{
public:
/**
* Constructor that sorts the input indices.
*
* \param i the first index.
* \param j The second index.
*/
Edge( unsigned int i, unsigned int j )
: std::pair< unsigned int, unsigned int >( i < j ? i : j, i < j ? j : i )
{
}
/**
* Compare two edges. This operator defines a weak ordering on the edges.
*
* \param e The edge to compare to.
* \return True, iff this edge is 'smaller' than the given edge.
*/
bool operator < ( Edge const& e ) const
{
return first < e.first || ( first == e.first && second < e.second );
}
/**
* Compare two edges for equality.
*
* \param e The edge to compare to.
* \return True, iff this edge has the same vertex indices as the given edge.
*/
bool operator == ( Edge const& e ) const
{
return first == e.first && second == e.second;
}
};
} // namespace utility
/**
* Tesselates an icosahedron in order to generate a triangle-based approximation of a sphere.
* The content of the provided vectors will be cleared and replaced. Triangle vertices are stored as
* successive values in the triangle vector.
*
* \param[out] vertices The vertices of the mesh.
* \param[out] triangles The resulting triangles as a list of indices into the vertex vector.
* \param level The tesselation level.
*/
void OWCOMMON_EXPORT tesselateIcosahedron( std::vector< WVector3D >* vertices, std::vector< unsigned int >* triangles, unsigned int level );
} // namespace wmath
#endif // WGEOMETRYFUNCTIONS_H
......@@ -37,6 +37,7 @@
#include "../common/WSharedObject.h"
#include "../common/WSharedSequenceContainer.h"
#include "WDataSetSingle.h"
#include "WDataSetScalar.h"
#include "WValueSet.h"
#include "WDataHandlerEnums.h"
......@@ -223,10 +224,18 @@ void WThreadedPerVoxelOperation< Value_T, numValues, Output_T, numOutputs >::com
template< typename Value_T, std::size_t numValues, typename Output_T, std::size_t numOutputs >
boost::shared_ptr< WDataSetSingle > WThreadedPerVoxelOperation< Value_T, numValues, Output_T, numOutputs >::getResult()
{
boost::shared_ptr< OutValueSetType > values( new OutValueSetType( 1, numOutputs,
m_output.getReadTicket()->get(),
DataType< Output_T >::type ) );
return boost::shared_ptr< WDataSetSingle >( new WDataSetSingle( values, m_grid ) );
boost::shared_ptr< OutValueSetType > values;
switch( numOutputs )
{
case 1:
values = boost::shared_ptr< OutValueSetType >( new OutValueSetType( 0, 1, m_output.getReadTicket()->get(),
DataType< Output_T >::type ) );
return boost::shared_ptr< WDataSetScalar >( new WDataSetScalar( values, m_grid ) );
default:
values = boost::shared_ptr< OutValueSetType >( new OutValueSetType( 1, numOutputs, m_output.getReadTicket()->get(),
DataType< Output_T >::type ) );
return boost::shared_ptr< WDataSetSingle >( new WDataSetSingle( values, m_grid ) );
}
}
#endif // WTHREADEDPERVOXELOPERATION_H
......@@ -190,7 +190,7 @@ boost::shared_ptr< WDataSet > WReaderNIfTI::load()
else if( header->dim[ 5 ] > 1 )
{
WAssert( header->dim[ 4 ] == 1, "Only scalar datasets are supported for time series so far." );
wlog::debug( "WReaderNIfTI" ) << "Load as time series";
wlog::debug( "WReaderNIfTI" ) << "Load as WDataSetTimeSeries";
std::size_t numTimeSlices = header->dim[ 5 ];
float tw = header->pixdim[ 5 ];
WAssert( tw != 0.0f, "" );
......@@ -288,14 +288,18 @@ boost::shared_ptr< WDataSet > WReaderNIfTI::load()
{
if( vDim == 3 )
{
wlog::debug( "WReaderNIfTI" ) << "Load as WDataSetVector";
newDataSet = boost::shared_ptr< WDataSet >( new WDataSetVector( newValueSet, newGrid ) );
}
else if( vDim == 1 )
{
wlog::debug( "WReaderNIfTI" ) << "Load as WDataSetScalar";
newDataSet = boost::shared_ptr< WDataSet >( new WDataSetScalar( newValueSet, newGrid ) );
}
else if( vDim > 20 && header->dim[ 5 ] == 1 ) // hardi data, order 1
else if( vDim > 20 && header->dim[ 5 ] <= 1 ) // hardi data, order 1
{
wlog::debug( "WReaderNIfTI" ) << "Load as WDataSetRawHARDI";
std::string gradientFileName = m_fname;
using wiotools::getSuffix;
std::string suffix = getSuffix( m_fname );
......@@ -322,6 +326,7 @@ boost::shared_ptr< WDataSet > WReaderNIfTI::load()
}
// cannot find the appropriate gradient vectors, build a dataSetSingle instead of hardi
newDataSet = boost::shared_ptr< WDataSet >( new WDataSetSingle( newValueSet, newGrid ) );
wlog::debug( "WReaderNIfTI" ) << "Could not find corresponding gradients file, loading as WDataSetSingle instead.";
}
else
{
......@@ -348,6 +353,7 @@ boost::shared_ptr< WDataSet > WReaderNIfTI::load()
}
else
{
wlog::debug( "WReaderNIfTI" ) << "Load as WDataSetSingle";
newDataSet = boost::shared_ptr< WDataSet >( new WDataSetSingle( newValueSet, newGrid ) );
}
}
......
......@@ -248,15 +248,21 @@ void WPickHandler::pick( osgViewer::View* view, const osgGA::GUIEventAdapter& ea
}
}
}
}
else // if( intersetionsExist )
}// end of if( intersetionsExist )
else
{
pickInfo = WPickInfo( "nothing", m_viewerName, wmath::WPosition( 0.0, 0.0, 0.0 ), std::make_pair( x, y ),
m_startPick.getModifierKey(), m_mouseButton, wmath::WVector3D( 0.0, 0.0, 0.0 ) );
// if we found no intersection and we have noting pickecd befor
// we want to return "nothing" in order to provide the pixel coordinates
// even though we did not hit anything.
if( m_startPick.getName() == "" )
{
pickInfo = WPickInfo( "nothing", m_viewerName, wmath::WPosition( 0.0, 0.0, 0.0 ), std::make_pair( x, y ),
m_startPick.getModifierKey(), m_mouseButton, wmath::WVector3D( 0.0, 0.0, 0.0 ) );
m_hitResult = pickInfo;
m_pickSignal( getHitResult() );
return;
m_hitResult = pickInfo;
m_pickSignal( getHitResult() );
return;
}
}
// Set the new pickInfo if the previously picked is still in list or we have a pick in conjunction with previously no pick
......@@ -273,7 +279,7 @@ void WPickHandler::pick( osgViewer::View* view, const osgGA::GUIEventAdapter& ea
pickNormal[1] = hitr->getWorldIntersectNormal()[1];
pickNormal[2] = hitr->getWorldIntersectNormal()[2];
pickInfo = WPickInfo( extractSuitableName( hitr ), m_viewerName, pickPos, std::make_pair( x, y ),
WPickInfo::NONE, m_mouseButton, pickNormal );
pickInfo.getModifierKey(), m_mouseButton, pickNormal );
}
// Use the old PickInfo with updated pixel info if we have previously picked something but the old is not in list anymore
......
......@@ -33,7 +33,6 @@
#include "../../common/WLimits.h"
#include "../../dataHandler/WDataHandler.h"
#include "../../dataHandler/WDataTexture3D.h"
#include "../../dataHandler/WThreadedPerVoxelOperation.h"
#include "../../kernel/WKernel.h"
#include "bermanTracking.xpm"
#include "WMBermanTracking.h"
......@@ -86,7 +85,7 @@ void WMBermanTracking::connectors()
"inSH", "A spherical harmonics dataset." )
);
m_inputGFA = boost::shared_ptr< WModuleInputData< WDataSetSingle > >( new WModuleInputData< WDataSetSingle >( shared_from_this(),
m_inputGFA = boost::shared_ptr< WModuleInputData< WDataSetScalar > >( new WModuleInputData< WDataSetScalar >( shared_from_this(),
"inGFA", "GFA." )
);
......@@ -94,7 +93,7 @@ void WMBermanTracking::connectors()
"inResiduals", "The residual HARDI data." )
);
m_output = boost::shared_ptr< WModuleOutputData< WDataSetSingle > >( new WModuleOutputData< WDataSetSingle >( shared_from_this(),
m_output = boost::shared_ptr< WModuleOutputData< WDataSetScalar > >( new WModuleOutputData< WDataSetScalar >( shared_from_this(),
"outProb", "The probabilistic tracking result." )
);
......@@ -103,10 +102,6 @@ void WMBermanTracking::connectors()
"outFibers", "The probabilistic fibers." )
);
m_outputGFA = boost::shared_ptr< WModuleOutputData< WDataSetSingle > >( new WModuleOutputData< WDataSetSingle >( shared_from_this(),
"outGFA", "The generalized fractional anisotropy map." )
);
addConnector( m_input );
addConnector( m_inputResidual );
......@@ -116,7 +111,6 @@ void WMBermanTracking::connectors()
// temporary
addConnector( m_outputFibers );
addConnector( m_outputGFA );
// call WModules initialization
WModule::connectors();
......@@ -162,6 +156,7 @@ void WMBermanTracking::moduleMain()
m_moduleState.setResetable( true, true );
m_moduleState.add( m_input->getDataChangedCondition() );
m_moduleState.add( m_inputResidual->getDataChangedCondition() );
m_moduleState.add( m_inputGFA->getDataChangedCondition() );
m_moduleState.add( m_propCondition );
m_moduleState.add( m_exceptionCondition );
......@@ -182,32 +177,37 @@ void WMBermanTracking::moduleMain()
boost::shared_ptr< WDataSetSphericalHarmonics > inData = m_input->getData();
boost::shared_ptr< WDataSetRawHARDI > inData2 = m_inputResidual->getData();
bool dataChanged = ( m_dataSet != inData ) || ( m_dataSetResidual != inData2 );
boost::shared_ptr< WDataSetScalar > inData3 = m_inputGFA->getData();
bool dataChanged = ( m_dataSet != inData ) || ( m_dataSetResidual != inData2 ) || ( m_gfa != inData3 );
if( dataChanged || !m_dataSet || !m_dataSetResidual )
if( dataChanged )
{
if( dataChanged )
if( m_trackingPool )
{
if( m_trackingPool )
m_trackingPool->stop();
m_trackingPool->wait();
if( m_currentProgress )
{
m_trackingPool->stop();
m_trackingPool->wait();
if( m_currentProgress )
{
m_currentProgress->finish();
}
m_currentProgress->finish();
}
}
m_dataSet = inData;
m_dataSetResidual = inData2;
m_gfa = inData3;
}
if( !m_dataSet || !m_dataSetResidual )
{
continue;
}
else
if( m_startTrigger->get( true ) == WPVBaseTypes::PV_TRIGGER_TRIGGERED )
{
m_startTrigger->set( WPVBaseTypes::PV_TRIGGER_READY, false );
if( m_dataSet && m_dataSetResidual && m_gfa )
{
if( m_result )
{
WDataHandler::deregisterDataSet( m_result );
}
// calculate some matrices needed for the residuals
std::vector< wmath::WUnitSphereCoordinates > c;
for( std::size_t i = 0; i < m_dataSetResidual->getOrientations().size(); ++i )
......@@ -222,39 +222,20 @@ void WMBermanTracking::moduleMain()
m_HMat = m_BMat * wmath::WSymmetricSphericalHarmonic::getSHFittingMatrix( c, 4, 0.0, false );
WAssert( m_HMat.getNbCols() == m_HMat.getNbRows(), "" );
// calculate the generalized fractional anisotropy for the whole dataset
if( !m_inputGFA->getData() )
{
calcGFA();
}
else
{
m_gfa = m_inputGFA->getData();
}
}
}
if( m_dataSet && m_dataSetResidual && ( dataChanged || m_startTrigger->get( true ) == WPVBaseTypes::PV_TRIGGER_TRIGGERED ) )
{
m_startTrigger->set( WPVBaseTypes::PV_TRIGGER_READY, false );
if( m_result )
{
WDataHandler::deregisterDataSet( m_result );
// get current properties
m_currentMinFA = m_minFA->get( true );
m_currentMinPoints = static_cast< std::size_t >( m_minPoints->get( true ) );
m_currentMinCos = m_minCos->get( true );
m_currentProbabilistic = m_probabilistic->get( true );
m_currentRatio = m_ratio->get( true );
m_currentEpsImpr = m_epsImpr->get( true );
// start tracking
resetTracking();
m_hits = boost::shared_ptr< std::vector< float > >( new std::vector< float >( m_dataSet->getGrid()->size(), 0.0 ) );
m_trackingPool->run();
debugLog() << "Running tracking function.";
}
// get current properties
m_currentMinFA = m_minFA->get( true );
m_currentMinPoints = static_cast< std::size_t >( m_minPoints->get( true ) );
m_currentMinCos = m_minCos->get( true );
m_currentProbabilistic = m_probabilistic->get( true );
m_currentRatio = m_ratio->get( true );
m_currentEpsImpr = m_epsImpr->get( true );
// start tracking
resetTracking();
m_hits = boost::shared_ptr< std::vector< float > >( new std::vector< float >( m_dataSet->getGrid()->size(), 0.0 ) );
m_trackingPool->run();
debugLog() << "Running tracking function.";
}
else if( m_trackingPool && m_trackingPool->status() == W_THREADS_FINISHED )
{
......@@ -268,9 +249,9 @@ void WMBermanTracking::moduleMain()
m_outputFibers->updateData( m_fiberSet );
m_trackingPool = boost::shared_ptr< TrackingFuncType >();
boost::shared_ptr< WValueSet< float > > vs( new WValueSet< float >( 1, 1, *m_hits, DataType< float >::type ) );
boost::shared_ptr< WValueSet< float > > vs( new WValueSet< float >( 0, 1, *m_hits, DataType< float >::type ) );
m_result = boost::shared_ptr< WDataSetSingle >( new WDataSetSingle( vs, m_dataSet->getGrid() ) );
m_result = boost::shared_ptr< WDataSetScalar >( new WDataSetScalar( vs, m_dataSet->getGrid() ) );
m_result->setFileName( "Berman_prob_tracking_result" );
m_result->getTexture()->setThreshold( 0.05f );
m_result->getTexture()->setSelectedColormap( 2 );
......@@ -298,60 +279,6 @@ void WMBermanTracking::moduleMain()
}
}
void WMBermanTracking::calcGFA()
{
if( m_gfa )
{
WDataHandler::deregisterDataSet( m_gfa );
}
m_gfa = boost::shared_ptr< WDataSetSingle >();
typedef WThreadedPerVoxelOperation< double, 15, double, 1 > GFAFunc;
typedef WThreadedFunction< GFAFunc > ThreadPool;
boost::shared_ptr< GFAFunc > gfafunc( new GFAFunc( m_dataSet, boost::bind( &This::perVoxelGFAFunc, this, _1 ) ) );
debugLog() << "Starting GFA computation.";
resetProgress( m_dataSet->getGrid()->size() );
ThreadPool t( W_AUTOMATIC_NB_THREADS, gfafunc );
t.subscribeExceptionSignal( boost::bind( &This::handleException, this, _1 ) );
t.run();
t.wait();
m_currentProgress->finish();
if( m_lastException )
{
throw WException( *m_lastException );
}
WAssert( t.status() == W_THREADS_FINISHED, "" );
debugLog() << "GFA computation finished.";
m_gfa = gfafunc->getResult();
m_gfa->getTexture()->setInterpolation( false );
m_gfa->setFileName( "generalized fractional anisotropy" );
m_outputGFA->updateData( m_gfa );
WDataHandler::registerDataSet( m_gfa );
}
boost::array< double, 1 > WMBermanTracking::perVoxelGFAFunc( WValueSet< double >::SubArray const& s )
{
++*m_currentProgress;
boost::array< double, 1 > a;
wmath::WValue< double > w( 15 );
for( int i = 0; i < 15; ++i )
{
w[ i ] = s[ i ];
}
wmath::WSymmetricSphericalHarmonic h( w );
a[ 0 ] = h.calcGFA( m_BMat );