← Back to branch summary

~ubuntu-branches/ubuntu/precise/openwalnut/precise

« back to all changes in this revision

Viewing changes to src/modules/eigenSystem/WMEigenSystem.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Sebastian Eichelbaum
  • Date: 2011-06-21 10:26:54 UTC
  • Revision ID: james.westby@ubuntu.com-20110621102654-rq0zf436q949biih
Tags: upstream-1.2.5
Import upstream version 1.2.5

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/modules/eigenSystem/WMEigenSystem.cpp
 
1
//---------------------------------------------------------------------------
 
2
//
 
3
// Project: OpenWalnut ( http://www.openwalnut.org )
 
4
//
 
5
// Copyright 2009 OpenWalnut Community, BSV-Leipzig and CNCF-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 <string>
 
26
 
 
27
#include "core/dataHandler/WDataSetDTI.h"
 
28
#include "core/dataHandler/WDataSetVector.h"
 
29
#include "core/kernel/WKernel.h"
 
30
#include "core/kernel/WModuleInputData.h"
 
31
#include "core/kernel/WModuleOutputData.h"
 
32
#include "WMEigenSystem.xpm"
 
33
#include "WMEigenSystem.h"
 
34
 
 
35
// This line is needed by the module loader to actually find your module. Do not remove. Do NOT add a ";" here.
 
36
W_LOADABLE_MODULE( WMEigenSystem )
 
37
 
 
38
WMEigenSystem::WMEigenSystem():
 
39
    WModule()
 
40
{
 
41
}
 
42
 
 
43
WMEigenSystem::~WMEigenSystem()
 
44
{
 
45
    // Cleanup!
 
46
}
 
47
 
 
48
boost::shared_ptr< WModule > WMEigenSystem::factory() const
 
49
{
 
50
    // See "src/modules/template/" for an extensively documented example.
 
51
    return boost::shared_ptr< WModule >( new WMEigenSystem() );
 
52
}
 
53
 
 
54
const char** WMEigenSystem::getXPMIcon() const
 
55
{
 
56
    return WMEigenSystem_xpm;
 
57
}
 
58
const std::string WMEigenSystem::getName() const
 
59
{
 
60
    return "Eigen System";
 
61
}
 
62
 
 
63
const std::string WMEigenSystem::getDescription() const
 
64
{
 
65
    // Specify your module description here. Be detailed. This text is read by the user.
 
66
    // See "src/modules/template/" for an extensively documented example.
 
67
    return "Someone should add some documentation here. "
 
68
    "Probably the best person would be the modules's creator, i.e. \"math\"";
 
69
}
 
70
 
 
71
void WMEigenSystem::connectors()
 
72
{
 
73
    // Put the code for your connectors here. See "src/modules/template/" for an extensively documented example.
 
74
    m_tensorIC = WModuleInputData< WDataSetDTI >::createAndAdd( shared_from_this(), "tensorInput", "The tensor field" );
 
75
 
 
76
    m_evecOutputs.push_back( WModuleOutputData< WDataSetVector >::createAndAdd( shared_from_this(), "evec1Output", "The 1. eigenvector field" ) );
 
77
    m_evecOutputs.push_back( WModuleOutputData< WDataSetVector >::createAndAdd( shared_from_this(), "evec2Output", "The 2. eigenvector field" ) );
 
78
    m_evecOutputs.push_back( WModuleOutputData< WDataSetVector >::createAndAdd( shared_from_this(), "evec3Output", "The 3. eigenvector field" ) );
 
79
    m_evecOutputs.push_back( WModuleOutputData< WDataSetVector >::createAndAdd( shared_from_this(), "evalsOutput", "All three eigenvalues" ) );
 
80
 
 
81
    m_evalOutputs.push_back( WModuleOutputData< WDataSetScalar >::createAndAdd( shared_from_this(), "eval1Output", "The 1. eigenvalue field" ) );
 
82
    m_evalOutputs.push_back( WModuleOutputData< WDataSetScalar >::createAndAdd( shared_from_this(), "eval2Output", "The 2. eigenvalue field" ) );
 
83
    m_evalOutputs.push_back( WModuleOutputData< WDataSetScalar >::createAndAdd( shared_from_this(), "eval3Output", "The 3. eigenvalue field" ) );
 
84
 
 
85
    WModule::connectors();
 
86
}
 
87
 
 
88
void WMEigenSystem::properties()
 
89
{
 
90
    // for selecting the strategy
 
91
    boost::shared_ptr< WItemSelection > strategies( new WItemSelection() );
 
92
    strategies->addItem( "LibEigen", "Eigensystem is computed via libEigen and its SelfAdjointEigenSolver" );
 
93
    strategies->addItem( "Jacobi", "Self implemented Jacobi iterative eigen decomposition" );
 
94
    m_strategySelector = m_properties->addProperty( "Strategy", "How the eigen system should be computed",
 
95
            strategies->getSelectorFirst() );
 
96
 
 
97
    WPropertyHelper::PC_SELECTONLYONE::addTo( m_strategySelector );
 
98
    WPropertyHelper::PC_NOTEMPTY::addTo( m_strategySelector );
 
99
 
 
100
    WModule::properties();
 
101
}
 
102
 
 
103
void WMEigenSystem::requirements()
 
104
{
 
105
    // Put the code for your requirements here. See "src/modules/template/" for an extensively documented example.
 
106
}
 
107
 
 
108
void WMEigenSystem::moduleMain()
 
109
{
 
110
    m_moduleState.setResetable( true, true );
 
111
    m_moduleState.add( m_tensorIC->getDataChangedCondition() );
 
112
    m_moduleState.add( m_strategySelector->getCondition() );
 
113
 
 
114
    ready();
 
115
 
 
116
    while( !m_shutdownFlag() )
 
117
    {
 
118
        infoLog() << "Waiting.";
 
119
        m_moduleState.wait();
 
120
 
 
121
        if( m_shutdownFlag() )
 
122
        {
 
123
            break;
 
124
        }
 
125
 
 
126
        boost::shared_ptr< WDataSetDTI > tensors = m_tensorIC->getData();
 
127
 
 
128
        if( !m_eigenPool && tensors )
 
129
        {
 
130
            // start the eigenvector computation if input is DTI double or float otherwise continue
 
131
            // when the computation finishes, we'll be notified by the threadspool's threadsDoneCondition
 
132
            resetEigenFunction( tensors );
 
133
            if( !m_eigenPool )
 
134
            {
 
135
                continue;
 
136
            }
 
137
            resetProgress( tensors->getValueSet()->size(), "Compute eigen system" );
 
138
            m_eigenPool->run();
 
139
            infoLog() << "Computing eigen systems...";
 
140
        }
 
141
        else if( m_eigenPool && m_eigenPool->status() == W_THREADS_FINISHED )
 
142
        {
 
143
            // the computation of the eigenvectors has finished we have a new field of eigen systems
 
144
            m_currentProgress->finish();
 
145
            WAssert( m_eigenOperationFloat || m_eigenOperationDouble, "Bug: No result is available. Checked double and float!" );
 
146
 
 
147
            if( tensors->getValueSet()->getDataType() == W_DT_DOUBLE )
 
148
            {
 
149
                updateOCs( m_eigenOperationDouble->getResult() );
 
150
            }
 
151
            else
 
152
            {
 
153
                updateOCs( m_eigenOperationFloat->getResult() );
 
154
            }
 
155
 
 
156
            m_eigenPool.reset();
 
157
            infoLog() << "Computing eigen systems done.";
 
158
        }
 
159
    }
 
160
}
 
161
 
 
162
void WMEigenSystem::updateOCs( boost::shared_ptr< const WDataSetSingle > es )
 
163
{
 
164
    WAssert( es, "Bug: updating the output with no data makes no sense." );
 
165
    boost::shared_ptr< const WValueSet< double > > vs = boost::shared_dynamic_cast< WValueSet< double > >( es->getValueSet() );
 
166
    WAssert( vs, "Bug: invalid value-set from WThreadedPerVoxelOperation dataset" );
 
167
    boost::shared_ptr< WGrid > grid = es->getGrid();
 
168
 
 
169
    typedef boost::shared_ptr< std::vector< double > > DataPointer;
 
170
    std::vector< DataPointer > vecdata( 4, DataPointer( new std::vector< double >( vs->size() * 3 ) ) );
 
171
    std::vector< DataPointer > valdata( 3, DataPointer( new std::vector< double >( vs->size() ) ) );
 
172
 
 
173
    for( size_t i = 0; i < vs->size(); ++i )
 
174
    {
 
175
        WValue< double > sys = vs->getWValue( i );
 
176
        // eigenvalues
 
177
        ( *vecdata[0] )[ i * 3 ] = sys[0];
 
178
        ( *valdata[0] )[ i ] = sys[0];
 
179
        ( *vecdata[0] )[ i * 3 + 1] = sys[4];
 
180
        ( *valdata[1] )[ i ] = sys[4];
 
181
        ( *vecdata[0] )[ i * 3 + 2] = sys[8];
 
182
        ( *valdata[2] )[ i ] = sys[8];
 
183
        // first eigenvector
 
184
        ( *vecdata[1] )[ i * 3 ] = sys[1];
 
185
        ( *vecdata[1] )[ i * 3 + 1] = sys[2];
 
186
        ( *vecdata[1] )[ i * 3 + 2] = sys[3];
 
187
        // second eigenvector
 
188
        ( *vecdata[2] )[ i * 3 ] = sys[5];
 
189
        ( *vecdata[2] )[ i * 3 + 1] = sys[6];
 
190
        ( *vecdata[2] )[ i * 3 + 2] = sys[7];
 
191
        // third eigenvector
 
192
        ( *vecdata[3] )[ i * 3 ] = sys[9];
 
193
        ( *vecdata[3] )[ i * 3 + 1] = sys[10];
 
194
        ( *vecdata[3] )[ i * 3 + 2] = sys[11];
 
195
    }
 
196
    typedef boost::shared_ptr< WDataSetVector > PDSV;
 
197
    typedef boost::shared_ptr< WDataSetScalar > PDSS;
 
198
    typedef WValueSet< double > WVSDBL;
 
199
    typedef boost::shared_ptr< WVSDBL > PWVSDBL;
 
200
 
 
201
    m_evecOutputs[0]->updateData( PDSV( new WDataSetVector( PWVSDBL( new WVSDBL( 1, 3, vecdata[0], W_DT_DOUBLE ) ), grid ) ) );
 
202
    m_evecOutputs[1]->updateData( PDSV( new WDataSetVector( PWVSDBL( new WVSDBL( 1, 3, vecdata[1], W_DT_DOUBLE ) ), grid ) ) );
 
203
    m_evecOutputs[2]->updateData( PDSV( new WDataSetVector( PWVSDBL( new WVSDBL( 1, 3, vecdata[2], W_DT_DOUBLE ) ), grid ) ) );
 
204
    m_evecOutputs[3]->updateData( PDSV( new WDataSetVector( PWVSDBL( new WVSDBL( 1, 3, vecdata[3], W_DT_DOUBLE ) ), grid ) ) );
 
205
    m_evalOutputs[0]->updateData( PDSS( new WDataSetScalar( PWVSDBL( new WVSDBL( 0, 1, valdata[0], W_DT_DOUBLE ) ), grid ) ) );
 
206
    m_evalOutputs[1]->updateData( PDSS( new WDataSetScalar( PWVSDBL( new WVSDBL( 0, 1, valdata[1], W_DT_DOUBLE ) ), grid ) ) );
 
207
    m_evalOutputs[2]->updateData( PDSS( new WDataSetScalar( PWVSDBL( new WVSDBL( 0, 1, valdata[2], W_DT_DOUBLE ) ), grid ) ) );
 
208
}
 
209
 
 
210
boost::array< double, 12 > WMEigenSystem::computeEigenSystem( WTensorSym< 2, 3, double > const& m ) const
 
211
{
 
212
    RealEigenSystem sys;
 
213
    jacobiEigenvector3D( m, &sys );
 
214
    boost::array< double, 12 > result = { { sys[0].first, sys[0].second[0],
 
215
                                                          sys[0].second[1],
 
216
                                                          sys[0].second[2],
 
217
                                            sys[1].first, sys[1].second[0],
 
218
                                                          sys[1].second[1],
 
219
                                                          sys[1].second[2],
 
220
                                            sys[2].first, sys[2].second[0],
 
221
                                                          sys[2].second[1],
 
222
                                                          sys[2].second[2] } }; // NOLINT
 
223
    return result;
 
224
}
 
225
 
 
226
WMEigenSystem::TPVOFloat::OutTransmitType const WMEigenSystem::eigenFuncFloat( TPVOFloat::TransmitType const& input )
 
227
{
 
228
    WTensorSym< 2, 3, double > m;
 
229
    m( 0, 0 ) = static_cast< double >( input[ 0 ] );
 
230
    m( 0, 1 ) = static_cast< double >( input[ 1 ] );
 
231
    m( 0, 2 ) = static_cast< double >( input[ 2 ] );
 
232
    m( 1, 1 ) = static_cast< double >( input[ 3 ] );
 
233
    m( 1, 2 ) = static_cast< double >( input[ 4 ] );
 
234
    m( 2, 2 ) = static_cast< double >( input[ 5 ] );
 
235
 
 
236
    ++*m_currentProgress;
 
237
 
 
238
    return computeEigenSystem( m );
 
239
}
 
240
 
 
241
WMEigenSystem::TPVODouble::OutTransmitType const WMEigenSystem::eigenFuncDouble( TPVODouble::TransmitType const& input )
 
242
{
 
243
    WTensorSym< 2, 3, double > m;
 
244
    m( 0, 0 ) = static_cast< double >( input[ 0 ] );
 
245
    m( 0, 1 ) = static_cast< double >( input[ 1 ] );
 
246
    m( 0, 2 ) = static_cast< double >( input[ 2 ] );
 
247
    m( 1, 1 ) = static_cast< double >( input[ 3 ] );
 
248
    m( 1, 2 ) = static_cast< double >( input[ 4 ] );
 
249
    m( 2, 2 ) = static_cast< double >( input[ 5 ] );
 
250
 
 
251
    ++*m_currentProgress;
 
252
 
 
253
    return computeEigenSystem( m );
 
254
}
 
255
 
 
256
WMEigenSystem::TPVODouble::OutTransmitType const WMEigenSystem::eigenSolverDouble( TPVODouble::TransmitType const &input )
 
257
{
 
258
    Eigen::Matrix3d m;
 
259
    m << input[0], input[1], input[2],
 
260
         input[1], input[3], input[4],
 
261
         input[2], input[4], input[5];
 
262
 
 
263
    return applyEigenSolver( m );
 
264
}
 
265
 
 
266
WMEigenSystem::TPVOFloat::OutTransmitType const WMEigenSystem::eigenSolverFloat( TPVOFloat::TransmitType const &input )
 
267
{
 
268
    Eigen::Matrix3d m;
 
269
    m << input[0], input[1], input[2],
 
270
         input[1], input[3], input[4],
 
271
         input[2], input[4], input[5];
 
272
 
 
273
    return applyEigenSolver( m );
 
274
}
 
275
 
 
276
boost::array< double, 12 > WMEigenSystem::applyEigenSolver( const Eigen::Matrix3d &m ) const
 
277
{
 
278
    Eigen::SelfAdjointEigenSolver< Eigen::Matrix3d > es( m );
 
279
    const Eigen::Matrix3d &evecs = es.eigenvectors();
 
280
    const Eigen::Vector3d &evals = es.eigenvalues();
 
281
    boost::array< double, 12 > result = { { evals( 0 ), evecs( 0, 0 ),
 
282
                                                        evecs( 1, 0 ),
 
283
                                                        evecs( 2, 0 ),
 
284
                                            evals( 1 ), evecs( 0, 1 ),
 
285
                                                        evecs( 1, 1 ),
 
286
                                                        evecs( 2, 1 ),
 
287
                                            evals( 2 ), evecs( 0, 2 ),
 
288
                                                        evecs( 1, 2 ),
 
289
                                                        evecs( 2, 2 ) } }; // NOLINT curly braces
 
290
    return result;
 
291
}
 
292
 
 
293
void WMEigenSystem::resetProgress( std::size_t todo, std::string name )
 
294
{
 
295
    if( m_currentProgress )
 
296
    {
 
297
        m_currentProgress->finish();
 
298
    }
 
299
    m_currentProgress = boost::shared_ptr< WProgress >( new WProgress( name, todo ) );
 
300
    m_progress->addSubProgress( m_currentProgress );
 
301
}
 
302
 
 
303
void WMEigenSystem::resetEigenFunction( boost::shared_ptr< WDataSetDTI > tensors )
 
304
{
 
305
    // check if we need to stop the currently running eigencomputation
 
306
    if( m_eigenPool )
 
307
    {
 
308
        debugLog() << "Stopping eigencomputation.";
 
309
        WThreadedFunctionStatus s = m_eigenPool->status();
 
310
        if( s != W_THREADS_FINISHED && s != W_THREADS_ABORTED )
 
311
        {
 
312
            m_eigenPool->stop();
 
313
            m_eigenPool->wait();
 
314
            s = m_eigenPool->status();
 
315
            WAssert( s == W_THREADS_FINISHED || s == W_THREADS_ABORTED, "" );
 
316
        }
 
317
        m_moduleState.remove( m_eigenPool->getThreadsDoneCondition() );
 
318
    }
 
319
    // the threadpool should have finished computing by now
 
320
 
 
321
    m_eigenOperationFloat = boost::shared_ptr< TPVOFloat >();
 
322
    m_eigenOperationDouble = boost::shared_ptr< TPVODouble >();
 
323
 
 
324
    // create a new one
 
325
    if( tensors->getValueSet()->getDataType() == W_DT_DOUBLE )
 
326
    {
 
327
        if( m_strategySelector->get().at( 0 )->getName() == "LibEigen" )
 
328
        {
 
329
            m_eigenOperationDouble = boost::shared_ptr< TPVODouble >( new TPVODouble( tensors, boost::bind( &WMEigenSystem::eigenSolverDouble, this, _1 ) ) ); // NOLINT line length
 
330
        }
 
331
        else if( m_strategySelector->get().at( 0 )->getName() == "Jacobi" )
 
332
        {
 
333
            m_eigenOperationDouble = boost::shared_ptr< TPVODouble >( new TPVODouble( tensors, boost::bind( &WMEigenSystem::eigenFuncDouble, this, _1 ) ) ); // NOLINT line length
 
334
        }
 
335
        else
 
336
        {
 
337
            WAssert( 0, "Bug: Invalid strategy for eigen decomposition selected!" );
 
338
        }
 
339
        m_eigenPool = boost::shared_ptr< WThreadedFunctionBase >( new EigenFunctionTypeDouble( W_AUTOMATIC_NB_THREADS, m_eigenOperationDouble ) );
 
340
        m_moduleState.add( m_eigenPool->getThreadsDoneCondition() );
 
341
    }
 
342
    else if( tensors->getValueSet()->getDataType() == W_DT_FLOAT )
 
343
    {
 
344
        if( m_strategySelector->get().at( 0 )->getName() == "LibEigen" )
 
345
        {
 
346
            m_eigenOperationFloat = boost::shared_ptr< TPVOFloat >( new TPVOFloat( tensors, boost::bind( &WMEigenSystem::eigenSolverFloat, this, _1 ) ) ); // NOLINT line length
 
347
        }
 
348
        else if( m_strategySelector->get().at( 0 )->getName() == "Jacobi" )
 
349
        {
 
350
            m_eigenOperationFloat = boost::shared_ptr< TPVOFloat >( new TPVOFloat( tensors, boost::bind( &WMEigenSystem::eigenFuncFloat, this, _1 ) ) ); // NOLINT line length
 
351
        }
 
352
        else
 
353
        {
 
354
            WAssert( 0, "Bug: Invalid strategy for eigen decomposition selected!" );
 
355
        }
 
356
        m_eigenPool = boost::shared_ptr< WThreadedFunctionBase >( new EigenFunctionTypeFloat( W_AUTOMATIC_NB_THREADS, m_eigenOperationFloat ) );
 
357
        m_moduleState.add( m_eigenPool->getThreadsDoneCondition() );
 
358
    }
 
359
    else
 
360
    {
 
361
        errorLog() << "Input data does not contain floating point values, skipping.";
 
362
        m_eigenPool.reset();
 
363
    }
 
364
}