~ubuntu-branches/ubuntu/vivid/openwalnut/vivid

« back to all changes in this revision

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

  • Committer: Package Import Robot
  • Author(s): Sebastian Eichelbaum
  • Date: 2014-03-19 17:46:12 UTC
  • mfrom: (3.1.2 sid)
  • Revision ID: package-import@ubuntu.com-20140319174612-e4mgtr1avbq3f7ph
Tags: 1.4.0~rc1+hg3a3147463ee2-1
* Major functionality and stability improvements.
* Several bug fixes
* Changed ttf-liberation dependency to fonts-liberation (Closes: #722405)
* OpenWalnut now works properly with OpenSceneGraph 3.2 (Closes: #718381)
* See http://www.openwalnut.org/versions/2

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
 
}