← Back to branch summary

~noskcaj/ubuntu/saucy/openwalnut/liberation

« back to all changes in this revision

Viewing changes to .pc/boost153/src/modules/HARDIToSphericalHarmonics/WSphericalHarmonicsCoefficientsThread.h

  • 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:
.pc/boost153/src/modules/HARDIToSphericalHarmonics/WSphericalHarmonicsCoefficientsThread.h
 
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
#ifndef WSPHERICALHARMONICSCOEFFICIENTSTHREAD_H
 
26
#define WSPHERICALHARMONICSCOEFFICIENTSTHREAD_H
 
27
 
 
28
#include <cmath>
 
29
#include <fstream>
 
30
#include <iostream>
 
31
#include <map>
 
32
#include <string>
 
33
#include <utility> // for std::pair
 
34
#include <vector>
 
35
 
 
36
#include <boost/math/special_functions/spherical_harmonic.hpp>
 
37
#include <boost/thread/thread.hpp>
 
38
 
 
39
#include "core/common/WLimits.h"
 
40
#include "core/common/WAssert.h"
 
41
#include "core/common/WProgress.h"
 
42
#include "core/common/WThreadedRunner.h"
 
43
#include "core/common/math/WMatrix.h"
 
44
 
 
45
#include "core/dataHandler/WDataSetRawHARDI.h"
 
46
#include "core/dataHandler/WDataSetSphericalHarmonics.h"
 
47
#include "core/dataHandler/WGridRegular3D.h"
 
48
 
 
49
// forward declaration
 
50
template< typename T = void >
 
51
class WSphericalHarmonicsCoefficientsThread;
 
52
 
 
53
/**
 
54
 * A specialization for void.
 
55
 */
 
56
template<>
 
57
class WSphericalHarmonicsCoefficientsThread< void >
 
58
{
 
59
public:
 
60
    /**
 
61
    * This structure is a collection of parameter and pointer to input and output data needed by each thread.
 
62
    */
 
63
    struct ThreadParameter
 
64
    {
 
65
        /**
 
66
         * Constructor, we need to initialize a reference.
 
67
         *
 
68
         * \param flag A reference to a shutdown flag that we should listen to.
 
69
         */
 
70
        ThreadParameter( WBoolFlag const& flag ) // NOLINT no explicit
 
71
            : m_shutdownFlag( flag )
 
72
        {
 
73
        }
 
74
 
 
75
        /**
 
76
         * Pointer to the HARDI measurements
 
77
         */
 
78
        boost::shared_ptr< WValueSetBase > m_valueSet;
 
79
 
 
80
        /**
 
81
         * Indices of nonzero gradients
 
82
         */
 
83
        std::vector< size_t > m_validIndices;
 
84
 
 
85
        /**
 
86
         * Indices of zero gradients (this should be normal T2 measurement)
 
87
         */
 
88
        std::vector< size_t > m_S0Indexes;
 
89
 
 
90
        /**
 
91
         * Output data, the spherical harmonics coefficients
 
92
         */
 
93
        boost::shared_ptr< std::vector<double> > m_data;
 
94
 
 
95
        /**
 
96
         * The order of the calculated spherical harmonics
 
97
         */
 
98
        int m_order;
 
99
 
 
100
        /**
 
101
         * Transformation-Matrix for conversion from HARDI measurements to spherical harmonics coefficients
 
102
         * (see Descoteaux dissertation)
 
103
         */
 
104
        boost::shared_ptr< WMatrix< double > > m_TransformMatrix;
 
105
 
 
106
        /**
 
107
         * Gradients of all measurements (including )
 
108
         */
 
109
        std::vector< WVector3d > m_gradients;
 
110
 
 
111
        /**
 
112
         * Pointer to progess indicator
 
113
         */
 
114
        boost::shared_ptr< WProgress > m_progress;
 
115
 
 
116
        /**
 
117
         * Indicate if the is error calculation is done.
 
118
         */
 
119
        bool m_doErrorCalculation;
 
120
 
 
121
        /**
 
122
         * Indicate if the Funk-Radon-Transformation is done.
 
123
         */
 
124
        bool m_doFunkRadonTransformation;
 
125
 
 
126
        /**
 
127
         * Indicate if the residuals will be calculated.
 
128
         * The residuals are the reprojection error. This means the
 
129
         */
 
130
        bool m_doResidualCalculation;
 
131
 
 
132
        /**
 
133
         * The stored residuals.
 
134
         */
 
135
        boost::shared_ptr< std::vector<double> > m_dataResiduals;
 
136
 
 
137
        /**
 
138
         * The b-value used during the creation of the HARDI-data.
 
139
         */
 
140
        double m_bDiffusionWeightingFactor;
 
141
 
 
142
        /**
 
143
         * Indicate if the normalisation from 0 to 1 is done.
 
144
         */
 
145
        bool m_normalize;
 
146
 
 
147
        /**
 
148
         * Indicate if the constant solid angle reconstruction is done.
 
149
         */
 
150
        bool m_csa;
 
151
 
 
152
        /**
 
153
         * A shutdownFlag that may tell the thread to stop.
 
154
         */
 
155
        WBoolFlag const& m_shutdownFlag;
 
156
    };
 
157
};
 
158
 
 
159
/**
 
160
 * Module for the creation of a spherical harmonic data set from raw HARDI data.
 
161
 * \ingroup modules
 
162
 */
 
163
template< typename T >
 
164
class WSphericalHarmonicsCoefficientsThread : public WThreadedRunner
 
165
{
 
166
public:
 
167
    /**
 
168
     * The constructor.
 
169
     * \param parameter collection of parameter and so on (description in more detail above)
 
170
     * \param range the range (start and end) of input data this thread should use
 
171
     */
 
172
    WSphericalHarmonicsCoefficientsThread( WSphericalHarmonicsCoefficientsThread<>::ThreadParameter parameter, std::pair< size_t, size_t > range );
 
173
 
 
174
    /**
 
175
     * The main function of the thread. Here the calculation of the spherical harmonics coefficients is done.
 
176
     */
 
177
    void threadMain();
 
178
 
 
179
    /**
 
180
     * Returns the average error of the thread.
 
181
     *
 
182
     * \return error measure
 
183
     */
 
184
    double getError() const;
 
185
 
 
186
private:
 
187
    double m_overallError;                                                  //!< accumulated error
 
188
    unsigned int m_errorCount;                                              //!< number of accumulated errors, necessary for average error calculation
 
189
    WSphericalHarmonicsCoefficientsThread<>::ThreadParameter m_parameter;   //!< collection of parameter and so on (description in more detail above)
 
190
    std::pair< size_t, size_t > m_range;                                    //!< the range (start and end) of input data this thread use
 
191
};
 
192
 
 
193
template< typename T >
 
194
WSphericalHarmonicsCoefficientsThread< T >::WSphericalHarmonicsCoefficientsThread(
 
195
    WSphericalHarmonicsCoefficientsThread<>::ThreadParameter parameter,
 
196
    std::pair< size_t, size_t > range )
 
197
    : m_parameter( parameter ),
 
198
      m_range( range )
 
199
{
 
200
}
 
201
 
 
202
template< typename T >
 
203
void WSphericalHarmonicsCoefficientsThread< T >::threadMain()
 
204
{
 
205
    m_errorCount = 0;
 
206
    m_overallError = 0.0;
 
207
 
 
208
    WMatrix<double> transformMatrix( *m_parameter.m_TransformMatrix );
 
209
    size_t l = ( m_parameter.m_order + 1 ) * ( m_parameter.m_order + 2 ) / 2;
 
210
 
 
211
    for( size_t i = m_range.first; i < m_range.second; i++ )
 
212
    {
 
213
        if( m_parameter.m_shutdownFlag() )
 
214
        {
 
215
            break;
 
216
        }
 
217
 
 
218
        // get measure vector
 
219
        boost::shared_ptr< WValueSet< T > > vs = boost::shared_dynamic_cast< WValueSet< T > >( m_parameter.m_valueSet );
 
220
        if( !vs )
 
221
        {
 
222
            throw WException( "Valueset pointer not valid." );
 
223
        }
 
224
 
 
225
        WValue< T > allMeasures( vs->getWValue( i ) );
 
226
 
 
227
        // extract measures for gradients != 0
 
228
        WValue< double > measures( m_parameter.m_validIndices.size() );
 
229
        unsigned int idx = 0;
 
230
 
 
231
        // find max S0 value
 
232
        double S0avg = 0.0;
 
233
        for( std::vector< size_t >::const_iterator it = m_parameter.m_S0Indexes.begin(); it != m_parameter.m_S0Indexes.end(); it++ )
 
234
        {
 
235
            S0avg += static_cast< double >( allMeasures[ *it ] );
 
236
        }
 
237
        S0avg /= m_parameter.m_S0Indexes.size();
 
238
 
 
239
        // to have a valid value for the average S0 signal
 
240
        if( S0avg <= 0.01 )
 
241
        {
 
242
            S0avg = 0.01;
 
243
        }
 
244
 
 
245
        // double minVal = 1e99;
 
246
        // double maxVal = -1e99;
 
247
        double thresholdDelta1 = 0.01;
 
248
        double thresholdDelta2 = 0.01;
 
249
        for( std::vector< size_t >::const_iterator it = m_parameter.m_validIndices.begin(); it != m_parameter.m_validIndices.end(); it++, idx++ )
 
250
        {
 
251
            if( m_parameter.m_csa )
 
252
            {
 
253
                double val = static_cast< double >( allMeasures[ *it ] ) / S0avg;
 
254
                if( val < 0.0)
 
255
                {
 
256
                    val = thresholdDelta1 / 2.0;
 
257
                }
 
258
                else if( val < thresholdDelta1 )
 
259
                {
 
260
                    val = thresholdDelta1 / 2.0 + val * val / ( 2.0 * thresholdDelta1 );
 
261
                }
 
262
                else if( val > 1.0 - thresholdDelta2 && val < 1.0 )
 
263
                {
 
264
                    val = 1.0 - thresholdDelta2 / 2.0 - std::pow( 1.0 - val, 2 ) / ( 2.0 * thresholdDelta2 );
 
265
                }
 
266
                else if( val >= 1.0 )
 
267
                {
 
268
                    val = 1.0 - thresholdDelta2 / 2.0;
 
269
                }
 
270
                measures[ idx ] = std::log( -std::log( val  ) );
 
271
            }
 
272
            else
 
273
            {
 
274
                measures[ idx ] = static_cast< double >( allMeasures[ *it ] ) / S0avg;
 
275
            }
 
276
        }
 
277
 
 
278
        WValue< double > coefficients( *m_parameter.m_TransformMatrix * measures );
 
279
 
 
280
        if( m_parameter.m_doResidualCalculation || m_parameter.m_doErrorCalculation )
 
281
        {
 
282
            WSymmetricSphericalHarmonic currentSphericalHarmonic( coefficients );
 
283
            for( idx = 0; idx < m_parameter.m_validIndices.size(); idx++ )
 
284
            {
 
285
                double error = static_cast< double >( measures[ idx ] )
 
286
                              - currentSphericalHarmonic.getValue( WUnitSphereCoordinates( m_parameter.m_gradients[ idx ] ) );
 
287
 
 
288
                if( m_parameter.m_doResidualCalculation )
 
289
                {
 
290
                    m_parameter.m_dataResiduals->operator[]( m_parameter.m_validIndices.size() * i + idx ) = error;
 
291
                }
 
292
                if( m_parameter.m_doErrorCalculation )
 
293
                {
 
294
                    m_overallError += fabs( error );
 
295
                    m_errorCount++;
 
296
                }
 
297
            }
 
298
        }
 
299
 
 
300
        // show progress
 
301
        ++( *( m_parameter.m_progress ) );
 
302
 
 
303
        // copy coefficients to output "data"
 
304
 
 
305
        // normalization
 
306
        double scale = 1.0;
 
307
        if( m_parameter.m_normalize )
 
308
        {
 
309
            scale *= std::sqrt( 4.0 * piDouble ) / coefficients[ 0 ];
 
310
        }
 
311
 
 
312
        if( m_parameter.m_csa )
 
313
        {
 
314
            coefficients[ 0 ] = 1.0 / ( 2.0 * std::sqrt( piDouble ) );
 
315
        }
 
316
 
 
317
        for( std::size_t j = 0; j < l; j++ )
 
318
        {
 
319
            m_parameter.m_data->operator[]( l * i + j ) = coefficients[ j ];
 
320
            // wlog::debug( "WSphericalHarmonicsCoefficientsThread" ) << "coefficients[" << j << "]: " << coefficients[ j ];
 
321
        }
 
322
    }
 
323
}
 
324
 
 
325
template< typename T >
 
326
double WSphericalHarmonicsCoefficientsThread< T >::getError() const
 
327
{
 
328
    if( m_errorCount == 0 ) return 0.0;
 
329
    return m_overallError / static_cast< double >( m_errorCount );
 
330
}
 
331
 
 
332
#endif  // WSPHERICALHARMONICSCOEFFICIENTSTHREAD_H