~noskcaj/ubuntu/saucy/openwalnut/liberation

« back to all changes in this revision

Viewing changes to .pc/boost153/src/modules/clusterSlicer/WMClusterSlicer.cpp

  • Committer: Package Import Robot
  • Author(s): Dmitrijs Ledkovs
  • Date: 2013-05-24 03:12:03 UTC
  • Revision ID: package-import@ubuntu.com-20130524031203-l5g1lzm1vd83fupi
Tags: 1.3.1+hg5849-1ubuntu1
Cherrypick boost1.53 pointer cast fixes.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
//---------------------------------------------------------------------------
 
2
//
 
3
// Project: OpenWalnut ( http://www.openwalnut.org )
 
4
//
 
5
// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
 
6
// For more information see http://www.openwalnut.org/copying
 
7
//
 
8
// This file is part of OpenWalnut.
 
9
//
 
10
// OpenWalnut is free software: you can redistribute it and/or modify
 
11
// it under the terms of the GNU Lesser General Public License as published by
 
12
// the Free Software Foundation, either version 3 of the License, or
 
13
// (at your option) any later version.
 
14
//
 
15
// OpenWalnut is distributed in the hope that it will be useful,
 
16
// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
17
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
18
// GNU Lesser General Public License for more details.
 
19
//
 
20
// You should have received a copy of the GNU Lesser General Public License
 
21
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
 
22
//
 
23
//---------------------------------------------------------------------------
 
24
 
 
25
#include <algorithm>
 
26
#include <functional>
 
27
#include <list>
 
28
#include <map>
 
29
#include <numeric>
 
30
#include <set>
 
31
#include <string>
 
32
#include <utility>
 
33
#include <vector>
 
34
 
 
35
#include <boost/shared_ptr.hpp>
 
36
 
 
37
#include <osg/Geode>
 
38
#include <osg/LightModel>
 
39
 
 
40
#include "core/common/math/WMath.h"
 
41
#include "core/common/math/WPlane.h"
 
42
#include "core/common/WAssert.h"
 
43
#include "core/common/WColor.h"
 
44
#include "core/common/WHistogramBasic.h"
 
45
#include "core/graphicsEngine/WGEGeodeUtils.h"
 
46
#include "core/graphicsEngine/WGEGeometryUtils.h"
 
47
#include "core/graphicsEngine/WGEUtils.h"
 
48
#include "core/graphicsEngine/WTriangleMesh.h"
 
49
#include "core/kernel/WKernel.h"
 
50
#include "WMClusterSlicer.h"
 
51
#include "WMClusterSlicer.xpm"
 
52
 
 
53
// This line is needed by the module loader to actually find your module.
 
54
W_LOADABLE_MODULE( WMClusterSlicer )
 
55
 
 
56
WMClusterSlicer::WMClusterSlicer()
 
57
    : WModule(),
 
58
      m_rootNode( osg::ref_ptr< WGEGroupNode >( new WGEGroupNode() ) ),
 
59
      m_fullUpdate( new WCondition ),
 
60
      m_maxMean( 0.0 ),
 
61
      m_minMean( 0.0 )
 
62
{
 
63
}
 
64
 
 
65
WMClusterSlicer::~WMClusterSlicer()
 
66
{
 
67
}
 
68
 
 
69
boost::shared_ptr< WModule > WMClusterSlicer::factory() const
 
70
{
 
71
    return boost::shared_ptr< WModule >( new WMClusterSlicer() );
 
72
}
 
73
 
 
74
const char** WMClusterSlicer::getXPMIcon() const
 
75
{
 
76
    return clusterSlicer_xpm;
 
77
}
 
78
 
 
79
const std::string WMClusterSlicer::getName() const
 
80
{
 
81
    return "Cluster Slicer";
 
82
}
 
83
 
 
84
const std::string WMClusterSlicer::getDescription() const
 
85
{
 
86
    return "Slices a cluster";
 
87
}
 
88
 
 
89
void WMClusterSlicer::connectors()
 
90
{
 
91
    m_fiberClusterIC = WModuleInputData< WFiberCluster >::createAndAdd( shared_from_this(), "clusterInput", "A cluster of fibers" );
 
92
    m_voxelizedClusterIC = WModuleInputData< WDataSetScalar >::createAndAdd( shared_from_this(), "clusterDSInput", "DataSet from cluster" );
 
93
    m_paramIC = WModuleInputData< WDataSetScalar >::createAndAdd( shared_from_this(), "paramInput", "DataSet of the parameters" );
 
94
    m_triangleMeshIC = WModuleInputData< WTriangleMesh >::createAndAdd( shared_from_this(), "meshInput", "TrianglMesh" );
 
95
    m_colorMapOC = WModuleOutputData< WColoredVertices >::createAndAdd( shared_from_this(), "colorMapOutput", "VertexID and colors" );
 
96
    m_triangleMeshOC = WModuleOutputData< WTriangleMesh >::createAndAdd( shared_from_this(), "meshOutput", "The Mesh to forward it for rendering" );
 
97
 
 
98
    WModule::connectors();
 
99
}
 
100
 
 
101
void WMClusterSlicer::properties()
 
102
{
 
103
    m_drawIsoVoxels     = m_properties->addProperty( "Show|Hide Iso Voxels", "Show|Hide voxels withing a given isosurface.", true );
 
104
    m_drawSlices        = m_properties->addProperty( "Show|Hide Slices", "Show|Hide slices along center line", false );
 
105
    m_isoValue          = m_properties->addProperty( "Iso value", "", 0.01 );
 
106
    m_meanSelector      = m_properties->addProperty( "Mean Type", "Selects the mean type, must be on of:"
 
107
                                                                   " 0==arithmetic, 1==geometric, 2==median", 2, m_fullUpdate );
 
108
    m_planeNumX         = m_properties->addProperty( "Planes #X-SamplePoints", "#samplePoints in first direction", 40, m_fullUpdate );
 
109
    m_planeNumY         = m_properties->addProperty( "Planes #Y-SamplePoints", "#samplePoints in second direction", 40, m_fullUpdate );
 
110
    m_planeStepWidth    = m_properties->addProperty( "Planes Step Width", "Distance between sample points", 0.5, m_fullUpdate );
 
111
    m_centerLineScale   = m_properties->addProperty( "#Planes", "Scales the center line to have more or less samples", 1.0, m_fullUpdate );
 
112
    m_selectBiggestComponentOnly = m_properties->addProperty( "Biggest Component Only",
 
113
       "If true, first the mesh is decomposed into its components (expensive!) and the biggest will be drawn", false );
 
114
    m_alternateColoring = m_properties->addProperty( "Alternate Mesh Coloring", "En/Disables the alternative mesh colorer", true, m_fullUpdate );
 
115
 
 
116
    m_meanSelector->setMin( 0 );
 
117
    m_meanSelector->setMax( 4 );
 
118
    m_isoValue->setMin( 0.0 );
 
119
    m_isoValue->setMax( 100.0 );
 
120
    m_planeNumX->setMin( 1 );
 
121
    m_planeNumY->setMin( 1 );
 
122
    m_planeStepWidth->setMin( 0.0 );
 
123
    m_customScale       = m_properties->addProperty( "Custom Scale", "", true, m_fullUpdate );
 
124
    m_minScale          = m_properties->addProperty( "MinScale", "Mean threshold below which is mapped to 0", 0.1, m_fullUpdate );
 
125
    m_maxScale          = m_properties->addProperty( "MaxScale", "Mean threshold above which is mapped to 1", 0.9, m_fullUpdate );
 
126
    m_minScale->setMin( 0.0 );
 
127
    m_minScale->setMax( 1.0 );
 
128
    m_maxScale->setMin( 0.0 );
 
129
    m_maxScale->setMax( 1.0 );
 
130
    m_minScaleColor     = m_properties->addProperty( "MinScaleColor", "", WColor( 1.0, 0.0, 0.0, 1.0 ), m_fullUpdate );
 
131
    m_maxScaleColor     = m_properties->addProperty( "MaxScaleColor", "", WColor( 1.0, 0.0, 0.0, 1.0 ), m_fullUpdate );
 
132
 
 
133
    WModule::properties();
 
134
}
 
135
 
 
136
void WMClusterSlicer::moduleMain()
 
137
{
 
138
    WKernel::getRunningKernel()->getGraphicsEngine()->getScene()->insert( m_rootNode );
 
139
 
 
140
    m_moduleState.setResetable( true, true );
 
141
    m_moduleState.add( m_paramIC->getDataChangedCondition() );
 
142
    m_moduleState.add( m_triangleMeshIC->getDataChangedCondition() );
 
143
    m_moduleState.add( m_drawSlices->getCondition() );
 
144
    m_moduleState.add( m_drawIsoVoxels->getCondition() );
 
145
    m_moduleState.add( m_fullUpdate );
 
146
    m_moduleState.add( m_selectBiggestComponentOnly->getCondition() );
 
147
 
 
148
    ready();
 
149
 
 
150
    while( !m_shutdownFlag() )
 
151
    {
 
152
        debugLog() << "Waiting...";
 
153
        m_moduleState.wait();
 
154
 
 
155
        if( m_shutdownFlag() )
 
156
        {
 
157
            break;
 
158
        }
 
159
 
 
160
        boost::shared_ptr< WDataSetScalar > newClusterDS = m_voxelizedClusterIC->getData();
 
161
        boost::shared_ptr< WFiberCluster >  newCluster   = m_fiberClusterIC->getData();
 
162
        boost::shared_ptr< WDataSetScalar > newParamDS   = m_paramIC->getData();
 
163
        boost::shared_ptr< WTriangleMesh > newMesh      = m_triangleMeshIC->getData();
 
164
        bool meshChanged = ( m_mesh != newMesh );
 
165
        bool paramDSChanged = ( m_paramDS != newParamDS );
 
166
        bool clusterChanged = ( m_cluster != newCluster );
 
167
        bool clusterDSChanged = ( m_clusterDS != newClusterDS );
 
168
        bool dataChanged = clusterDSChanged || clusterChanged || paramDSChanged || meshChanged;
 
169
 
 
170
        m_clusterDS = newClusterDS;
 
171
        m_cluster = newCluster;
 
172
        m_paramDS = newParamDS;
 
173
        m_mesh = newMesh;
 
174
 
 
175
        if( !( m_clusterDS.get() && m_cluster.get() && m_paramDS.get() && m_mesh.get() ) )
 
176
        {
 
177
            debugLog() << "Invalid data. Waiting for data change again.";
 
178
            continue;
 
179
        }
 
180
 
 
181
        if( dataChanged )
 
182
        {
 
183
            infoLog() << "Coverage for isovalue: " << m_isoValue->get() << " is: " << countTractPointsInsideVolume( m_isoValue->get() );
 
184
            infoLog() << "Recommended isovalue for specified coverage: " << computeOptimalIsoValue();
 
185
        }
 
186
 
 
187
        if( dataChanged )
 
188
        {
 
189
            debugLog() << "Data changed...";
 
190
            if( meshChanged )
 
191
            {
 
192
                if( clusterDSChanged )
 
193
                {
 
194
                    debugLog() << "Building Join Tree";
 
195
                    m_joinTree = boost::shared_ptr< WJoinContourTree >( new WJoinContourTree( m_clusterDS ) );
 
196
                    m_joinTree->buildJoinTree();
 
197
                    debugLog() << "Finished building Join Tree";
 
198
                }
 
199
                // when the mesh has changed the isoValue must have had changed too
 
200
                m_isoVoxels = m_joinTree->getVolumeVoxelsEnclosedByIsoSurface( m_isoValue->get() );
 
201
                if( m_selectBiggestComponentOnly->get( true ) )
 
202
                {
 
203
                    debugLog() << "Start mesh decomposition";
 
204
                    m_components = tm_utils::componentDecomposition( *m_mesh );
 
205
                    debugLog() << "Decomposing mesh done";
 
206
                }
 
207
            }
 
208
        }
 
209
 
 
210
        if( meshChanged || paramDSChanged || m_meanSelector->changed() || m_planeNumX->changed() || m_alternateColoring->changed()
 
211
                        || m_planeNumY->changed() || m_planeStepWidth->changed() || m_centerLineScale->changed() || m_customScale->changed()
 
212
                        || m_minScale->changed() || m_maxScale->changed() || m_minScaleColor->changed() || m_maxScaleColor->changed() )
 
213
        {
 
214
            debugLog() << "Performing full update";
 
215
            generateSlices();
 
216
            sliceAndColorMesh( m_mesh );
 
217
            updateDisplay( true ); // force display update here (erasing invalid planes)
 
218
            debugLog() << "Full update done.";
 
219
        }
 
220
        else if( m_drawSlices->changed() || m_drawIsoVoxels->changed() )
 
221
        {
 
222
            updateDisplay();
 
223
        }
 
224
        else if( m_selectBiggestComponentOnly->changed() )
 
225
        {
 
226
            sliceAndColorMesh( m_mesh );
 
227
        }
 
228
    }
 
229
 
 
230
    WKernel::getRunningKernel()->getGraphicsEngine()->getScene()->remove( m_rootNode );
 
231
}
 
232
 
 
233
WValue< double > WMClusterSlicer::meanParameter( boost::shared_ptr< std::set< WPosition > > samplePoints ) const
 
234
{
 
235
    std::vector< double > samples;
 
236
    samples.reserve( samplePoints->size() );
 
237
 
 
238
    boost::shared_ptr< WGridRegular3D > grid = boost::shared_dynamic_cast< WGridRegular3D >( m_clusterDS->getGrid() );
 
239
    assert( grid != 0 );
 
240
 
 
241
    for( std::set< WPosition >::iterator pos = samplePoints->begin(); pos != samplePoints->end(); )
 
242
    {
 
243
        // ATTENTION: We don't interpolate in m_clusterDS since it might be the wrong component!!
 
244
        int id = grid->getVoxelNum( *pos );
 
245
        if( id >= 0 ) // check if the position is inside of the main component (biggest volume)
 
246
        {
 
247
            if( m_isoVoxels->count( static_cast< size_t >( id ) ) == 1 ) // check if its in main component
 
248
            {
 
249
                bool inParamGrid = false;
 
250
                double isoValue = m_clusterDS->interpolate( *pos, &inParamGrid );
 
251
                if( inParamGrid )
 
252
                {
 
253
                    if( isoValue > m_isoValue->get() )
 
254
                    {
 
255
                        inParamGrid = false;
 
256
                        double value = m_paramDS->interpolate( *pos, &inParamGrid );
 
257
                        if( inParamGrid ) // check if its in paramDS
 
258
                        {
 
259
                            samples.push_back( value );
 
260
                            ++pos;
 
261
                            continue;
 
262
                        }
 
263
                    }
 
264
                }
 
265
            }
 
266
        }
 
267
        samplePoints->erase( pos++ ); // erase in case the pos is not in main component or in paramDS or in clusterDS
 
268
    }
 
269
    double max = *std::max_element( samples.begin(), samples.end() );
 
270
    double min = *std::min_element( samples.begin(), samples.end() );
 
271
 
 
272
    double arithmeticMean = std::accumulate( samples.begin(), samples.end(), 0.0 );
 
273
    arithmeticMean = arithmeticMean / ( samples.empty() ? 1.0 : samples.size() );
 
274
 
 
275
    std::nth_element( samples.begin(), samples.begin() + samples.size() / 2, samples.end() );
 
276
    double median = ( samples.empty() ? 0.0 : samples[ samples.size() / 2 ] );
 
277
 
 
278
    // discard first all elements <= 0.0 since then the geometric mean does not make any sense
 
279
 
 
280
    std::vector< double >::iterator newEnd = std::remove_if( samples.begin(), samples.end(), std::bind2nd( std::less_equal< double >(), 0.0 ) );
 
281
    double geometricMean = std::accumulate( samples.begin(), newEnd, 1.0, std::multiplies< double >() );
 
282
    geometricMean = std::pow( geometricMean, ( samples.empty() ? 0.0 : 1.0 / samples.size() ) );
 
283
 
 
284
    WValue< double > result( 5 );
 
285
    result[0] = arithmeticMean;
 
286
    result[1] = geometricMean;
 
287
    result[2] = median;
 
288
    result[3] = max;
 
289
    result[4] = min;
 
290
    return result;
 
291
}
 
292
 
 
293
void WMClusterSlicer::generateSlices()
 
294
{
 
295
    debugLog() << "Generating Slices";
 
296
    WFiber centerLine( *m_cluster->getCenterLine() ); // copy the centerline
 
297
    if( centerLine.empty() )
 
298
    {
 
299
        errorLog() << "CenterLine of the bundle is empty => no slices are drawn";
 
300
        return;
 
301
    }
 
302
    centerLine.resampleByNumberOfPoints( static_cast< size_t >( m_centerLineScale->get( true ) * centerLine.size() ) );
 
303
 
 
304
    m_slices = boost::shared_ptr< std::vector< std::pair< double, WPlane > > >( new std::vector< std::pair< double, WPlane > > );
 
305
    m_maxMean = wlimits::MIN_DOUBLE;
 
306
    m_minMean = wlimits::MAX_DOUBLE;
 
307
    const int nbX = m_planeNumX->get( true );
 
308
    const int nbY = m_planeNumY->get( true );
 
309
    const double stepWidth = m_planeStepWidth->get( true );
 
310
    const int meanType = m_meanSelector->get( true );
 
311
    m_rootNode->remove( m_samplePointsGeode );
 
312
    m_samplePointsGeode = osg::ref_ptr< WGEGroupNode >( new WGEGroupNode ); // discard old geode
 
313
 
 
314
    WVector3d generator = ( centerLine.front() - midPoint( centerLine ) );
 
315
    generator = cross( generator, centerLine.back() - midPoint( centerLine ) );
 
316
    for( size_t i = 1; i < centerLine.size(); ++i )
 
317
    {
 
318
        WVector3d tangent = centerLine[i] - centerLine[i-1];
 
319
        WVector3d first = cross( tangent, generator );
 
320
        WVector3d second = cross( tangent, first );
 
321
        boost::shared_ptr< WPlane > p = boost::shared_ptr< WPlane >( new WPlane( tangent, centerLine[i-1], first, second ) );
 
322
 
 
323
        boost::shared_ptr< std::set< WPosition > > samplePoints = p->samplePoints( stepWidth, nbX, nbY  );
 
324
        double mean = meanParameter( samplePoints )[ meanType ];
 
325
        if( mean > m_maxMean )
 
326
        {
 
327
            m_maxMean = mean;
 
328
        }
 
329
        if( mean < m_minMean )
 
330
        {
 
331
            m_minMean = mean;
 
332
        }
 
333
        m_slices->push_back( std::make_pair( mean, WPlane( *p ) ) );
 
334
        if( m_drawSlices->get( true ) )
 
335
        {
 
336
            m_samplePointsGeode->insert( wge::genPointBlobs( samplePoints, 0.1 ) );
 
337
        }
 
338
    }
 
339
    infoLog() << "Max-Mean: " << m_maxMean;
 
340
    infoLog() << "Min-Mean: " << m_minMean;
 
341
    if( m_drawSlices->get( true ) )
 
342
    {
 
343
        m_rootNode->insert( m_samplePointsGeode );
 
344
    }
 
345
}
 
346
 
 
347
/**
 
348
 * Compares two WTrianglesMeshes on their size of vertices. This is private here since it makes sense only to this module ATM.
 
349
 */
 
350
struct WMeshSizeComp
 
351
{
 
352
    /**
 
353
     * Comparator on num vertex of two WTriangleMeshes
 
354
     *
 
355
     * \param m First Mesh
 
356
     * \param n Second Mesh
 
357
     *
 
358
     * \return True if and only if the first Mesh has less vertices as the second mesh.
 
359
     */
 
360
    bool operator()( const boost::shared_ptr< WTriangleMesh >& m, const boost::shared_ptr< WTriangleMesh >& n ) const
 
361
    {
 
362
        return m->vertSize() < n->vertSize();
 
363
    }
 
364
};
 
365
 
 
366
void WMClusterSlicer::sliceAndColorMesh( boost::shared_ptr< WTriangleMesh > mesh )
 
367
{
 
368
    debugLog() << "Selecting mesh component...";
 
369
    boost::shared_ptr< WTriangleMesh > renderMesh = mesh;
 
370
    assert( renderMesh.get() );
 
371
    if( m_selectBiggestComponentOnly->get( true ) )
 
372
    {
 
373
        if( !m_components.get() )
 
374
        {
 
375
            debugLog() << "Start mesh decomposition";
 
376
            m_components = tm_utils::componentDecomposition( *renderMesh );
 
377
            debugLog() << "Decomposing mesh done";
 
378
        }
 
379
        renderMesh = *std::max_element( m_components->begin(), m_components->end(), WMeshSizeComp() );
 
380
    }
 
381
 
 
382
    debugLog() << "Mesh selected";
 
383
 
 
384
    if( renderMesh != m_triangleMeshOC->getData() )
 
385
    {
 
386
        m_triangleMeshOC->updateData( renderMesh );
 
387
    }
 
388
 
 
389
    debugLog() << "Building mesh color map...";
 
390
    m_colorMap = boost::shared_ptr< WColoredVertices >( new WColoredVertices );
 
391
    std::map< size_t, WColor > cmData;
 
392
 
 
393
    if( !m_alternateColoring->get( true ) )
 
394
    {
 
395
//        std::map< size_t, std::pair< double, int > > cm;
 
396
//
 
397
//        for( std::vector< std::pair< double, WPlane > >::const_iterator slice = m_slices->begin(); slice != m_slices->end(); ++slice )
 
398
//        {
 
399
//            boost::shared_ptr< std::set< size_t > > coloredVertices = tm_utils::intersection( *renderMesh, slice->second );
 
400
//            double scaledMean = mapMeanOntoScale( slice->first );
 
401
//            for( std::set< size_t >::const_iterator coloredVertex = coloredVertices->begin(); coloredVertex != coloredVertices->end(); ++coloredVertex ) // NOLINT
 
402
//            {
 
403
//                if( cm.find( *coloredVertex ) != cm.end() )
 
404
//                {
 
405
//                    cm[ *coloredVertex ].first += scaledMean;
 
406
//                    cm[ *coloredVertex ].second++;
 
407
//                }
 
408
//                else
 
409
//                {
 
410
//                    cm[ *coloredVertex ].first  = scaledMean;
 
411
//                    cm[ *coloredVertex ].second = 1;
 
412
//                }
 
413
//            }
 
414
//        }
 
415
//
 
416
//        for( std::map< size_t, std::pair< double, int > >::const_iterator vertexColor = cm.begin(); vertexColor != cm.end(); ++vertexColor )
 
417
//        {
 
418
//            cmData[ vertexColor->first ] = WColor( 0.0, vertexColor->second.first / vertexColor->second.second, 1.0, 1.0 );
 
419
//        }
 
420
    }
 
421
    else
 
422
    {
 
423
        osg::ref_ptr< osg::Vec3Array > vertices = renderMesh->getVertexArray();
 
424
        for( size_t i = 0; i < vertices->size(); ++i )
 
425
        {
 
426
            WAssert( m_slices->size() > 2, "We only support alternative coloring with more than 2 slices" );
 
427
            WPosition vertex( ( *vertices)[i].x(), ( *vertices)[i].y(), ( *vertices)[i].z() );
 
428
            std::vector< PlanePair > planePairs = computeNeighbouringPlanePairs( vertex );
 
429
            if( !planePairs.empty() )
 
430
            {
 
431
                PlanePair closestPlanes = closestPlanePair( planePairs, vertex );
 
432
                if( closestPlanes.first != 0 || closestPlanes.second != 0 ) // if(0,0) then it may be a boundary vertex
 
433
                {
 
434
                    cmData[ i ] = colorFromPlanePair( vertex, closestPlanes );
 
435
                }
 
436
            }
 
437
        }
 
438
    }
 
439
 
 
440
    m_colorMap->setData( cmData );
 
441
    debugLog() << "Done with color map building";
 
442
    m_colorMapOC->updateData( m_colorMap );
 
443
}
 
444
 
 
445
double WMClusterSlicer::mapMeanOntoScale( double meanValue ) const
 
446
{
 
447
    if( m_customScale->get( true ) )
 
448
    {
 
449
        if( meanValue < m_minScale->get( true ) )
 
450
        {
 
451
            return 0.0;
 
452
        }
 
453
        else if( meanValue > m_maxScale->get( true ) )
 
454
        {
 
455
            return 1.0;
 
456
        }
 
457
        else
 
458
        {
 
459
            return ( m_maxScale->get() == m_minScale->get() ? 0.0 : ( meanValue - m_minScale->get() ) / ( m_maxScale->get() - m_minScale->get() ) );
 
460
        }
 
461
    }
 
462
    else
 
463
    {
 
464
        return ( m_maxMean == m_minMean ? 0.0 : ( meanValue - m_minMean ) / ( m_maxMean - m_minMean ) );
 
465
    }
 
466
}
 
467
 
 
468
bool WMClusterSlicer::isInBetween( const WPosition& vertex, const PlanePair& pp ) const
 
469
{
 
470
    const WPlane& p = ( *m_slices )[ pp.first ].second;
 
471
    const WPlane& q = ( *m_slices )[ pp.second ].second;
 
472
    double r = dot( p.getNormal(), vertex - p.getPosition() );
 
473
    double s = dot( q.getNormal(), vertex - q.getPosition() );
 
474
    if( signum( r ) != signum( s ) )
 
475
    {
 
476
        return true;
 
477
    }
 
478
    return false;
 
479
}
 
480
 
 
481
std::vector< WMClusterSlicer::PlanePair > WMClusterSlicer::computeNeighbouringPlanePairs( const WPosition& vertex ) const
 
482
{
 
483
    // assume base point/origin of every plane is on center line
 
484
    std::vector< PlanePair > result;
 
485
    for( size_t i = 1; i < m_slices->size(); ++i )
 
486
    {
 
487
        if( isInBetween( vertex, std::make_pair( i, i-1 ) ) )
 
488
        {
 
489
            result.push_back( std::make_pair( i, i-1 ) );
 
490
        }
 
491
    }
 
492
    return result;
 
493
}
 
494
 
 
495
WMClusterSlicer::PlanePair WMClusterSlicer::closestPlanePair( const std::vector< PlanePair >& pairs, const WPosition& vertex ) const
 
496
{
 
497
    double minDistance = wlimits::MAX_DOUBLE;
 
498
    PlanePair result( 0, 0 );
 
499
    for( std::vector< PlanePair >::const_iterator pp = pairs.begin(); pp != pairs.end(); ++pp )
 
500
    {
 
501
        const WPlane& p = ( *m_slices )[ pp->first ].second;
 
502
        const WPlane& q = ( *m_slices )[ pp->second ].second;
 
503
        double sqDistance = std::min( length2( vertex - p.getPosition() ),  length2( vertex - q.getPosition() ) );
 
504
        if( minDistance > sqDistance )
 
505
        {
 
506
            minDistance = sqDistance;
 
507
            result = *pp;
 
508
        }
 
509
    }
 
510
    // check if coloring is necessary
 
511
    const WPlane& firstPlane = m_slices->front().second;
 
512
    const WPlane& lastPlane =  m_slices->back().second;
 
513
 
 
514
    double distanceToFirst = dot( firstPlane.getNormal(), vertex - firstPlane.getPosition() );
 
515
    double distanceToLast  = dot( lastPlane.getNormal(), vertex - lastPlane.getPosition() );
 
516
    if( ( distanceToFirst < 0 && length2( vertex - firstPlane.getPosition() ) < minDistance ) ||
 
517
        ( distanceToLast > 0 && length2( vertex - lastPlane.getPosition() ) < minDistance ) )
 
518
    {
 
519
        return std::make_pair( 0, 0 );
 
520
    }
 
521
 
 
522
    return result;
 
523
}
 
524
 
 
525
WColor WMClusterSlicer::colorFromPlanePair( const WPosition& vertex, const PlanePair& pp ) const
 
526
{
 
527
    const WPlane& p = ( *m_slices )[ pp.first ].second;
 
528
    const WPlane& q = ( *m_slices )[ pp.second ].second;
 
529
    double distanceToP = std::abs( dot( p.getNormal(), vertex - p.getPosition() ) );
 
530
    double distanceToQ = std::abs( dot( q.getNormal(), vertex - q.getPosition() ) );
 
531
    // std::cout << "distP, distQ: " << distanceToP << " " << distanceToQ << std::endl;
 
532
    double colorP = ( *m_slices )[ pp.first ].first;
 
533
    double colorQ = ( *m_slices )[ pp.second ].first;
 
534
    double vertexColor = colorQ * ( distanceToP / ( distanceToP + distanceToQ ) ) + colorP * ( distanceToQ / ( distanceToP + distanceToQ ) );
 
535
    // std::cout << "colorP, colorQ, vertexColor: " << colorP << " " << colorQ << " " << vertexColor << std::endl;
 
536
    return WColor( 0.0, mapMeanOntoScale( vertexColor ), 1.0, 1.0 );
 
537
}
 
538
 
 
539
void WMClusterSlicer::updateDisplay( bool force )
 
540
{
 
541
    debugLog() << "Forced updating display: " << force;
 
542
    if( m_drawIsoVoxels->changed() || force )
 
543
    {
 
544
        m_rootNode->remove( m_isoVoxelGeode );
 
545
        m_isoVoxelGeode = osg::ref_ptr< osg::Geode >( new osg::Geode() ); // discard old geode
 
546
        if( m_drawIsoVoxels->get( true ) )
 
547
        {
 
548
            assert( m_isoVoxels && "JoinTree cannot be valid since there is no valid m_clusterDS." );
 
549
            m_isoVoxelGeode = generateIsoVoxelGeode();
 
550
            m_rootNode->insert( m_isoVoxelGeode );
 
551
        }
 
552
    }
 
553
 
 
554
    if( m_drawSlices->changed() || force )
 
555
    {
 
556
        m_rootNode->remove( m_sliceGeode );
 
557
 
 
558
        generateSlices(); // for recomputation of sample point geode
 
559
 
 
560
        m_sliceGeode = osg::ref_ptr< WGEGroupNode >( new WGEGroupNode ); // discard old geode
 
561
        if( m_drawSlices->get( true ) ) // regenerate
 
562
        {
 
563
            const double width = m_planeNumX->get() * m_planeStepWidth->get();
 
564
            const double height = m_planeNumY->get() * m_planeStepWidth->get();
 
565
            for( std::vector< std::pair< double, WPlane > >::const_iterator cit = m_slices->begin(); cit != m_slices->end(); ++cit )
 
566
            {
 
567
                double scaledMean = mapMeanOntoScale( cit->first );
 
568
                WColor color( scaledMean, scaledMean, 1.0, 1.0 );
 
569
                m_sliceGeode->insert( wge::genFinitePlane( width, height, cit->second, color, true ) );
 
570
            }
 
571
            m_rootNode->insert( m_sliceGeode );
 
572
        }
 
573
    }
 
574
}
 
575
 
 
576
osg::ref_ptr< osg::Geode > WMClusterSlicer::generateIsoVoxelGeode() const
 
577
{
 
578
    boost::shared_ptr< std::set< WPosition > > points( new std::set< WPosition > );
 
579
    boost::shared_ptr< WGridRegular3D > grid = boost::shared_dynamic_cast< WGridRegular3D >( m_clusterDS->getGrid() );
 
580
    assert( grid != 0 );
 
581
    for( std::set< size_t >::const_iterator cit = m_isoVoxels->begin(); cit != m_isoVoxels->end(); ++cit )
 
582
    {
 
583
        points->insert( grid->getPosition( *cit ) );
 
584
    }
 
585
    return wge::genPointBlobs< std::set< WPosition > >( points, grid->getOffsetX() );
 
586
}
 
587
 
 
588
 
 
589
double WMClusterSlicer::countTractPointsInsideVolume( double isoValue ) const
 
590
{
 
591
    const std::list< size_t >& indices = m_cluster->getIndices();
 
592
    boost::shared_ptr< const WDataSetFiberVector > tracts = m_cluster->getDataSetReference();
 
593
    std::list< size_t >::const_iterator iter;
 
594
    size_t pointCounter = 0;
 
595
    size_t pointsInside = 0;
 
596
    for( iter = indices.begin(); iter != indices.end(); ++iter )
 
597
    {
 
598
        const WFiber& tract = ( *tracts )[ *iter ];
 
599
        WFiber::const_iterator pos;
 
600
        for( pos = tract.begin(); pos != tract.end(); ++pos, ++pointCounter )
 
601
        {
 
602
            bool inGrid = false;
 
603
            double value = m_clusterDS->interpolate( *pos, &inGrid );
 
604
            if( inGrid && value >= isoValue )
 
605
            {
 
606
                ++pointsInside;
 
607
            }
 
608
        }
 
609
    }
 
610
    return static_cast< double >( pointsInside ) / static_cast< double >( pointCounter );
 
611
}
 
612
 
 
613
double WMClusterSlicer::computeOptimalIsoValue( double coverage ) const
 
614
{
 
615
    const std::list< size_t >& indices = m_cluster->getIndices();
 
616
    boost::shared_ptr< const WDataSetFiberVector > tracts = m_cluster->getDataSetReference();
 
617
    std::list< size_t >::const_iterator iter;
 
618
 
 
619
    WHistogramBasic h( m_clusterDS->getMin(), m_clusterDS->getMax() );
 
620
 
 
621
    for( iter = indices.begin(); iter != indices.end(); ++iter )
 
622
    {
 
623
        const WFiber& tract = ( *tracts )[ *iter ];
 
624
        WFiber::const_iterator pos;
 
625
        for( pos = tract.begin(); pos != tract.end(); ++pos )
 
626
        {
 
627
            bool inGrid = false;
 
628
            double value = m_clusterDS->interpolate( *pos, &inGrid );
 
629
            if( inGrid )
 
630
            {
 
631
                h.insert( value );
 
632
            }
 
633
        }
 
634
    }
 
635
    size_t valuesSoFar = 0;
 
636
    int64_t index = h.size() - 1;
 
637
    while( static_cast< double >( valuesSoFar ) / static_cast< double >( h.valuesSize() ) < coverage && index >= 0 )
 
638
    {
 
639
        valuesSoFar += h[index--];
 
640
    }
 
641
    return  m_clusterDS->getMin() + std::abs( m_clusterDS->getMin() - m_clusterDS->getMax() ) * index / static_cast< double >( h.size() );
 
642
}
 
643
 
 
644
void WMClusterSlicer::activate()
 
645
{
 
646
    if( m_rootNode )
 
647
    {
 
648
        if( m_active->get() )
 
649
        {
 
650
            m_rootNode->setNodeMask( 0xFFFFFFFF );
 
651
        }
 
652
        else
 
653
        {
 
654
            m_rootNode->setNodeMask( 0x0 );
 
655
        }
 
656
    }
 
657
 
 
658
    WModule::activate();
 
659
}