WMDetTractClustering.cpp 15.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
//---------------------------------------------------------------------------
//
// 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 <algorithm>
#include <iomanip>
#include <iostream>
28
#include <list>
29
#include <sstream>
30 31 32 33 34 35
#include <string>
#include <vector>

#include <osg/Geode>
#include <osg/Geometry>

36 37
#include "../../common/datastructures/WDXtLookUpTable.h"
#include "../../common/datastructures/WFiber.h"
38
#include "../../common/WAssert.h"
39
#include "../../common/WColor.h"
40
#include "../../common/WIOTools.h"
41
#include "../../common/WLogger.h"
42
#include "../../common/WProgress.h"
43
#include "../../common/WStringUtils.h"
44
#include "../../dataHandler/datastructures/WFiberCluster.h"
45
#include "../../dataHandler/exceptions/WDHIOFailure.h"
46 47
#include "../../dataHandler/io/WReaderLookUpTableVTK.h"
#include "../../dataHandler/io/WWriterLookUpTableVTK.h"
48 49
#include "../../dataHandler/WDataSetFiberVector.h"
#include "../../dataHandler/WSubject.h"
50
#include "../../graphicsEngine/WGEUtils.h"
51
#include "../../kernel/WKernel.h"
52 53
#include "detTractClustering.xpm"
#include "WMDetTractClustering.h"
54

55 56 57
// This line is needed by the module loader to actually find your module.
W_LOADABLE_MODULE( WMDetTractClustering )

58
WMDetTractClustering::WMDetTractClustering()
59
    : WModule(),
60
      m_lastTractsSize( 0 ),
61
      m_dLtTableExists( false ),
62
      m_update( new WCondition() )
63 64 65
{
}

66
WMDetTractClustering::~WMDetTractClustering()
67 68 69
{
}

70
boost::shared_ptr< WModule > WMDetTractClustering::factory() const
71
{
72
    return boost::shared_ptr< WModule >( new WMDetTractClustering() );
73 74
}

75
const char** WMDetTractClustering::getXPMIcon() const
76
{
77
    return detTractClustering_xpm;
78 79
}

80
void WMDetTractClustering::moduleMain()
81
{
82 83
    m_moduleState.setResetable( true, true ); // modules state remembers fired events while not waiting
    m_moduleState.add( m_tractInput->getDataChangedCondition() );
84
    m_moduleState.add( m_update );
Mathias Goldau's avatar
Mathias Goldau committed
85

86 87
    ready();

88
    while ( !m_shutdownFlag() )
89
    {
90 91
        m_moduleState.wait();

92
        if ( !m_tractInput->getData().get() ) // ok, the output has not yet sent data
Mathias Goldau's avatar
Mathias Goldau committed
93 94 95
        {
            continue;
        }
96
        if( m_rawTracts != m_tractInput->getData() ) // in case data has changed
97
        {
98
            m_rawTracts = m_tractInput->getData();
99 100
            boost::shared_ptr< WProgress > convertProgress( new WProgress( "Converting tracts", 1 ) );
            m_progress->addSubProgress( convertProgress );
101
            m_tracts = boost::shared_ptr< WDataSetFiberVector >( new WDataSetFiberVector( m_rawTracts ) );
102 103
            m_numTracts->set( static_cast< int32_t >( m_tracts->size() ) );
            convertProgress->finish();
104 105
        }

106
        if( m_run->get( true ) == WPVBaseTypes::PV_TRIGGER_TRIGGERED )
107
        {
108
            infoLog() << "Start processing tracts";
109
            update();
110 111
            infoLog() << "Processing finished";
            m_run->set( WPVBaseTypes::PV_TRIGGER_READY, false );
112
        }
113

114 115
        if( m_clusterOutputID->changed() )
        {
116
            m_clusterOutputID->get( true );
117
            updateOutput();
118 119 120 121
        }
    }
}

122
void WMDetTractClustering::properties()
123
{
124 125 126 127 128
    m_maxDistance_t   = m_properties->addProperty( "Max cluster distance", "Maximum distance of two tracts in one cluster.", 6.5 );
    m_proximity_t     = m_properties->addProperty( "Min point distance", "Min distance of points of two tracts which should be considered", 0.0 );
    m_minClusterSize  = m_properties->addProperty( "Min cluster size", "Minium of tracts per cluster", 10 );
    m_clusterOutputID = m_properties->addProperty( "Output cluster ID", "This cluster ID will be connected to the output.", 0, m_update );
    m_run             = m_properties->addProperty( "Start clustering", "Start", WPVBaseTypes::PV_TRIGGER_READY, m_update );
129 130 131 132 133 134 135 136 137 138 139 140 141 142

    // information properties
    m_numTracts = m_infoProperties->addProperty( "#Tracts:", "Number of tracts beeing processed", 0 );
    m_numTracts->setMin( 0 );
    m_numTracts->setMax( wlimits::MAX_INT32_T );
    m_numUsedTracts = m_infoProperties->addProperty( "#Tracts used:", "Number of tracts beeing finally used in the clustering", 0 );
    m_numUsedTracts->setMin( 0 );
    m_numUsedTracts->setMax( wlimits::MAX_INT32_T );
    m_numClusters = m_infoProperties->addProperty( "#Clusters:", "Number of clusters beeing computed", 0 );
    m_numClusters->setMin( 0 );
    m_numClusters->setMax( wlimits::MAX_INT32_T );
    m_numValidClusters = m_infoProperties->addProperty( "#Clusters used:", "Number of clusters beeing bigger than the given threshold", 0 );
    m_numValidClusters->setMin( 0 );
    m_numValidClusters->setMax( wlimits::MAX_INT32_T );
143
    m_clusterSizes = m_infoProperties->addProperty( "Cluster sizes:", "Size of each valid cluster", std::string() );
144
}
Mathias Goldau's avatar
Mathias Goldau committed
145

146
void WMDetTractClustering::updateOutput()
147 148 149
{
    if( m_clusters.empty() )
    {
150
        warnLog() << "There are no clusters for the output, leave output connecter untouched";
151 152 153
        return;
    }

154
    if( m_clusterOutputID->get() >= static_cast< int >( m_clusters.size() ) || m_clusterOutputID->get() < 0 )
155
    {
156 157
        warnLog() << "Invalid cluster ID for output selected: " << m_clusterOutputID->get() << " using default ID 0";
        m_clusterOutputID->set( 0, true );
Mathias Goldau's avatar
Mathias Goldau committed
158
    }
159
    m_output->updateData( boost::shared_ptr< WFiberCluster >( new WFiberCluster( m_clusters[ m_clusterOutputID->get() ] ) ) );
Mathias Goldau's avatar
Mathias Goldau committed
160 161
}

162
void WMDetTractClustering::update()
Mathias Goldau's avatar
Mathias Goldau committed
163
{
164
    if( !( m_dLtTableExists = dLtTableExists() ) )
165
    {
166
        debugLog() << "Consider old table as invalid.";
167
        m_dLtTable.reset( new WDXtLookUpTable( m_tracts->size() ) );
168 169
    }

170
    cluster();
171 172 173 174

    boost::shared_ptr< WProgress > saveProgress( new WProgress( "Saving tracts", 1 ) );
    m_progress->addSubProgress( saveProgress );
    if( !wiotools::fileExists( lookUpTableFileName() ) )
175
    {
176 177 178 179 180 181 182
        WWriterLookUpTableVTK w( lookUpTableFileName(), true );
        try
        {
            w.writeTable( m_dLtTable->getData(), m_lastTractsSize );
        }
        catch( const WDHIOFailure& e )
        {
183
            errorLog() << "Could not write dlt file: " << e.what();
184
        }
185
    }
186 187 188 189 190 191 192 193
    saveProgress->finish();

    // reset min max of the slider for the cluster ID
    m_clusterOutputID->setMin( 0 );
    m_clusterOutputID->setMax( m_clusters.size() - 1 );

    WKernel::getRunningKernel()->getGraphicsEngine()->getScene()->remove( m_osgNode );
    m_osgNode = paint();
194
    WKernel::getRunningKernel()->getGraphicsEngine()->getScene()->insert( m_osgNode );
195 196

    // TODO(math): For reasons of simplicity (no multiple input connectors possible) just forward one cluster to the voxelizer
197 198
    for( size_t i = 0; i < m_clusters.size(); ++i )
    {
199
        m_clusters[ i ].setDataSetReference( m_tracts );
200
        m_clusters[ i ].generateCenterLine();
201 202
    }

203
    updateOutput();
204 205
}

206
bool WMDetTractClustering::dLtTableExists()
207
{
208 209 210
    boost::shared_ptr< WProgress > readProgress( new WProgress( "Try to read dLt table", 1 ) );
    m_progress->addSubProgress( readProgress );

211 212 213
    std::string dLtFileName = lookUpTableFileName();

    if( wiotools::fileExists( dLtFileName ) )
214 215 216
    {
        try
        {
217 218
            debugLog() << "trying to read table from: " << dLtFileName;
            WReaderLookUpTableVTK r( dLtFileName );
219
            boost::shared_ptr< std::vector< double > > data( new std::vector< double >() );
220 221
            r.readTable( data );
            m_dLtTable.reset( new WDXtLookUpTable( static_cast< size_t >( data->back() ) ) );
222
            m_lastTractsSize = static_cast< size_t >( data->back() );
223 224 225 226 227

            // remove the dimension from data array since it's not representing any distance
            data->pop_back();

            m_dLtTable->setData( *data );
228
            readProgress->finish();
229
            return true;
230
        }
Mathias Goldau's avatar
Mathias Goldau committed
231
        catch( const WDHException& e )
232
        {
233
            debugLog() << e.what() << std::endl;
234 235
        }
    }
236
    readProgress->finish();
237
    return false;
238 239
}

240
void WMDetTractClustering::cluster()
241
{
242 243 244 245
    double proximity_t = m_proximity_t->get();
    double maxDistance_t = m_maxDistance_t->get();
    size_t minClusterSize = m_minClusterSize->get();

246
    size_t numTracts = m_tracts->size();
247

248 249 250 251
    infoLog() << "Start clustering with " << numTracts << " tracts.";

    m_clusters.clear();  // remove evtl. old clustering

252
    m_clusterIDs = std::vector< size_t >( numTracts, 0 );
253

254
    for( size_t i = 0; i < numTracts; ++i )
255
    {
256 257
        m_clusters.push_back( WFiberCluster( i ) );
        m_clusterIDs[i] = i;
258
    }
259

260
    boost::shared_ptr< WProgress > progress( new WProgress( "Tract clustering", numTracts ) );
261
    m_progress->addSubProgress( progress );
262

263 264 265 266
    boost::function< double ( const wmath::WFiber& q, const wmath::WFiber& r ) > dLt; // NOLINT
    const double proxSquare = proximity_t * proximity_t;
    dLt = boost::bind( wmath::WFiber::distDLT, proxSquare, _1, _2 );

267
    for( size_t q = 0; q < numTracts; ++q )  // loop over all "symmetric" tract pairs
268
    {
269
        for( size_t r = q + 1;  r < numTracts; ++r )
270
        {
271
            if( m_clusterIDs[q] != m_clusterIDs[r] )  // both tracts are in different clusters
272 273 274
            {
                if( !m_dLtTableExists )
                {
275
                    (*m_dLtTable)( q, r ) = dLt( (*m_tracts)[q], (*m_tracts)[r] );
276
                }
277
                if( (*m_dLtTable)( q, r ) < maxDistance_t )  // q and r provide an inter-cluster-link
278
                {
279
                    meld( m_clusterIDs[q], m_clusterIDs[r] );
280 281 282
                }
            }
        }
283
        ++*progress;
284
    }
285
    progress->finish();
286 287
    m_dLtTableExists = true;

288 289 290
    boost::shared_ptr< WProgress > eraseProgress( new WProgress( "Erasing clusters", 1 ) );
    m_progress->addSubProgress( eraseProgress );

291
    // remove empty clusters
292
    WFiberCluster emptyCluster;
293 294 295 296 297 298 299
    m_clusters.erase( std::remove( m_clusters.begin(), m_clusters.end(), emptyCluster ), m_clusters.end() );

    // determine #clusters and #small_clusters which are below a certain size
    size_t numSmallClusters = 0;
    for( size_t i = 0; i < m_clusters.size(); ++i )
    {
        m_clusters[i].sort();  // Refactor: why do we need sorting here?
300
        if( m_clusters[i].size() < minClusterSize )
301 302 303 304 305
        {
            m_clusters[i].clear();  // make small clusters empty to remove them easier
            ++numSmallClusters;
        }
    }
306
    m_numClusters->set( static_cast< int32_t >( m_clusters.size() ) );
307
    m_clusters.erase( std::remove( m_clusters.begin(), m_clusters.end(), emptyCluster ), m_clusters.end() );
308
    m_numValidClusters->set( static_cast< int32_t >( m_clusters.size() ) );
309

310
    m_lastTractsSize = m_tracts->size();
311
    eraseProgress->finish();
312 313
}

314
osg::ref_ptr< osg::Geode > WMDetTractClustering::genTractGeode( const WFiberCluster &cluster, const WColor& color ) const
315 316 317 318 319 320 321 322
{
    using osg::ref_ptr;
    ref_ptr< osg::Vec3Array > vertices = ref_ptr< osg::Vec3Array >( new osg::Vec3Array );
    ref_ptr< osg::Geometry > geometry = ref_ptr< osg::Geometry >( new osg::Geometry );

    std::list< size_t >::const_iterator cit = cluster.getIndices().begin();
    for( ; cit !=  cluster.getIndices().end(); ++cit )
    {
323
        const wmath::WFiber &fib = (*m_tracts)[ *cit ];
324 325 326 327
        for( size_t i = 0; i < fib.size(); ++i )
        {
            vertices->push_back( osg::Vec3( fib[i][0], fib[i][1], fib[i][2] ) );
        }
328
        geometry->addPrimitiveSet( new osg::DrawArrays( osg::PrimitiveSet::LINE_STRIP, vertices->size() - fib.size(), fib.size() ) );
329 330 331 332 333
    }

    geometry->setVertexArray( vertices );

    ref_ptr< osg::Vec4Array > colors = ref_ptr< osg::Vec4Array >( new osg::Vec4Array );
334
    colors->push_back( wge::osgColor( color ) );
335 336 337 338
    geometry->setColorArray( colors );
    geometry->setColorBinding( osg::Geometry::BIND_OVERALL );
    osg::ref_ptr< osg::Geode > geode = osg::ref_ptr< osg::Geode >( new osg::Geode );
    geode->addDrawable( geometry.get() );
339

340 341 342 343
    return geode;
}


344
osg::ref_ptr< WGEManagedGroupNode > WMDetTractClustering::paint() const
345
{
346
    // get different colors via HSV color model for each cluster
347 348 349 350
    double hue = 0.0;
    double hue_increment = 1.0 / m_clusters.size();
    WColor color;

351 352
    osg::ref_ptr< WGEManagedGroupNode > result = osg::ref_ptr< WGEManagedGroupNode >( new WGEManagedGroupNode( m_active ) );
    size_t numUsedTracts = 0;
353
    std::stringstream clusterLog;
354 355 356
    for( size_t i = 0; i < m_clusters.size(); ++i, hue += hue_increment )
    {
        color.setHSV( hue, 1.0, 0.75 );
357
        result->insert( genTractGeode( m_clusters[i], color ).get() );
358
        clusterLog << m_clusters[i].size() << " ";
359
        numUsedTracts += m_clusters[i].size();
360
    }
361 362
    m_clusterSizes->set( std::string( clusterLog.str() ) );
    m_numUsedTracts->set( static_cast< int32_t >( numUsedTracts ) );
363 364
    result->getOrCreateStateSet()->setMode( GL_LIGHTING, osg::StateAttribute::OFF );
    return result;
365
}
366

367
void WMDetTractClustering::meld( size_t qClusterID, size_t rClusterID )
368 369 370 371 372 373 374 375 376 377
{
    // first choose the cluster with the smaller ID
    if( qClusterID > rClusterID )  // merge always to the cluster with the smaller id
    {
        std::swap( qClusterID, rClusterID );
    }

    WFiberCluster& qCluster = m_clusters[ qClusterID ];
    WFiberCluster& rCluster = m_clusters[ rClusterID ];

378
    WAssert( !qCluster.empty() && !rCluster.empty(), "At least one cluster was empty while trying to merge!" );
379 380 381 382 383 384 385 386 387 388 389 390

    // second update m_clusterIDs array
    std::list< size_t >::const_iterator cit = rCluster.getIndices().begin();
    std::list< size_t >::const_iterator cit_end = rCluster.getIndices().end();
    for( ; cit != cit_end; ++cit )
    {
        m_clusterIDs[*cit] = qClusterID;
    }

    // and at last merge them
    qCluster.merge( rCluster );

391
    WAssert( rCluster.empty(), "The right cluster was no empty after melting!" );
392
}
Mathias Goldau's avatar
Mathias Goldau committed
393

394
void WMDetTractClustering::connectors()
Mathias Goldau's avatar
Mathias Goldau committed
395
{
396
    typedef WModuleInputData< WDataSetFibers > InputData;  // just an alias
Mathias Goldau's avatar
Mathias Goldau committed
397
    typedef WModuleOutputData< WFiberCluster > OutputData; // -"-
Mathias Goldau's avatar
Mathias Goldau committed
398

399 400
    m_tractInput = boost::shared_ptr< InputData >( new InputData( shared_from_this(), "tractInput", "A deterministic tract dataset." ) );
    m_output = boost::shared_ptr< OutputData >( new OutputData( shared_from_this(), "clusterOutput", "A set of tract indices aka cluster" ) );
Mathias Goldau's avatar
Mathias Goldau committed
401

402
    addConnector( m_tractInput );
Mathias Goldau's avatar
Mathias Goldau committed
403
    addConnector( m_output );
Mathias Goldau's avatar
Mathias Goldau committed
404 405 406
    WModule::connectors();  // call WModules initialization
}

407
std::string WMDetTractClustering::lookUpTableFileName() const
408
{
409 410
    std::stringstream newExtension;
    newExtension << std::fixed << std::setprecision( 2 );
411
    newExtension << ".pt-" << m_proximity_t->get() << ".dlt";
412 413
    boost::filesystem::path tractFileName( m_tracts->getFileName() );
    return tractFileName.replace_extension( newExtension.str() ).string();
414
}
415 416


417
WMDetTractClustering::OutputIDBound::OutputIDBound( const std::vector< WFiberCluster >& clusters )
418 419 420 421
    : m_clusters( clusters )
{
}

422
bool WMDetTractClustering::OutputIDBound::accept( boost::shared_ptr< WPropertyVariable< WPVBaseTypes::PV_INT > > /* property */,
423 424 425 426
                                               WPVBaseTypes::PV_INT value )
{
    return ( value >= 0 ) && ( value < static_cast< int >( m_clusters.size() ) );
}