← Back to branch summary

~noskcaj/ubuntu/saucy/openwalnut/liberation

« back to all changes in this revision

Viewing changes to .pc/boost153/src/modules/anisotropicFiltering/WMAnisotropicFiltering.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:
.pc/boost153/src/modules/anisotropicFiltering/WMAnisotropicFiltering.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 <algorithm>
 
26
#include <string>
 
27
#include <sstream>
 
28
#include <vector>
 
29
 
 
30
#include "core/kernel/WKernel.h"
 
31
#include "core/dataHandler/WDataHandler.h"
 
32
#include "core/common/WPropertyHelper.h"
 
33
 
 
34
#include "WMAnisotropicFiltering.xpm"
 
35
#include "WMAnisotropicFiltering.h"
 
36
 
 
37
// This line is needed by the module loader to actually find your module.
 
38
W_LOADABLE_MODULE( WMAnisotropicFiltering )
 
39
 
 
40
WMAnisotropicFiltering::WMAnisotropicFiltering():
 
41
    WModule()
 
42
{
 
43
    m_k = 1.0;
 
44
}
 
45
 
 
46
WMAnisotropicFiltering::~WMAnisotropicFiltering()
 
47
{
 
48
}
 
49
 
 
50
boost::shared_ptr< WModule > WMAnisotropicFiltering::factory() const
 
51
{
 
52
    return boost::shared_ptr< WModule >( new WMAnisotropicFiltering() );
 
53
}
 
54
 
 
55
const std::string WMAnisotropicFiltering::getName() const
 
56
{
 
57
    return "Anisotropic Filter";
 
58
}
 
59
 
 
60
const std::string WMAnisotropicFiltering::getDescription() const
 
61
{
 
62
    return "Filters MRI data, preserving edges.";
 
63
}
 
64
 
 
65
const char** WMAnisotropicFiltering::getXPMIcon() const
 
66
{
 
67
    return WMAnisotropicFiltering_xpm;
 
68
}
 
69
 
 
70
void WMAnisotropicFiltering::connectors()
 
71
{
 
72
    m_input = boost::shared_ptr< WModuleInputData < WDataSetSingle  > >(
 
73
        new WModuleInputData< WDataSetSingle >( shared_from_this(), "in", "The input dataset." ) );
 
74
    addConnector( m_input );
 
75
 
 
76
    m_output = boost::shared_ptr< WModuleOutputData < WDataSetSingle  > >(
 
77
        new WModuleOutputData< WDataSetSingle >( shared_from_this(), "out", "The extracted image." ) );
 
78
    addConnector( m_output );
 
79
 
 
80
    WModule::connectors();
 
81
}
 
82
 
 
83
void WMAnisotropicFiltering::properties()
 
84
{
 
85
    m_propCondition = boost::shared_ptr< WCondition >( new WCondition() );
 
86
 
 
87
    m_iterations = m_properties->addProperty( "#Iterations", "Number of iterations.", 10, m_propCondition );
 
88
    m_iterations->setMax( 1000 );
 
89
 
 
90
    m_Kcoeff = m_properties->addProperty( "K", "The diffusion weighting coefficient. Increase this to better preserve edges.", 9.0, m_propCondition );
 
91
    m_Kcoeff->setMin( 0.1 );
 
92
    m_Kcoeff->setMax( 1000.0 );
 
93
 
 
94
    m_delta = m_properties->addProperty( "delta", "The time step for integration.", 0.1, m_propCondition );
 
95
    m_delta->setMax( 10.0 );
 
96
 
 
97
    WModule::properties();
 
98
}
 
99
 
 
100
void WMAnisotropicFiltering::activate()
 
101
{
 
102
    WModule::activate();
 
103
}
 
104
 
 
105
void WMAnisotropicFiltering::moduleMain()
 
106
{
 
107
    m_moduleState.setResetable( true, true );
 
108
    m_moduleState.add( m_input->getDataChangedCondition() );
 
109
    m_moduleState.add( m_propCondition );
 
110
 
 
111
    ready();
 
112
 
 
113
    while( !m_shutdownFlag() )
 
114
    {
 
115
        m_moduleState.wait();
 
116
 
 
117
        if( m_shutdownFlag() )
 
118
        {
 
119
            break;
 
120
        }
 
121
 
 
122
        boost::shared_ptr< WDataSetSingle > newDataSet = m_input->getData();
 
123
        bool dataChanged = ( m_dataSet != newDataSet );
 
124
        bool dataValid   = ( newDataSet );
 
125
 
 
126
        if( dataValid )
 
127
        {
 
128
            if( dataChanged || m_iterations->changed() || m_Kcoeff->changed() )
 
129
            {
 
130
                m_dataSet = newDataSet;
 
131
                WAssert( m_dataSet, "" );
 
132
                WAssert( m_dataSet->getValueSet(), "" );
 
133
 
 
134
                infoLog() << "Calculating.";
 
135
 
 
136
                m_k = m_Kcoeff->get( true );
 
137
                m_d = m_delta->get( true );
 
138
 
 
139
                calcSmoothedImages( m_iterations->get( true ) );
 
140
 
 
141
                infoLog() << "Calculation complete.";
 
142
            }
 
143
        }
 
144
    }
 
145
 
 
146
    debugLog() << "Shutting down...";
 
147
    debugLog() << "Finished! Good Bye!";
 
148
}
 
149
 
 
150
void WMAnisotropicFiltering::calcSmoothedImages( int iterations )
 
151
{
 
152
    if( iterations < 1 )
 
153
        return;
 
154
 
 
155
    std::size_t numImages = m_dataSet->getValueSet()->rawSize() / m_dataSet->getGrid()->size();
 
156
    infoLog() << "Images: " << numImages;
 
157
 
 
158
    boost::shared_ptr< std::vector< double > > smoothed( new std::vector< double >( m_dataSet->getValueSet()->rawSize() ) );
 
159
    boost::shared_ptr< WGridRegular3D > grid = boost::shared_dynamic_cast< WGridRegular3D >( m_dataSet->getGrid() );
 
160
 
 
161
    // fill in data from dataset
 
162
    copyData( smoothed, grid );
 
163
 
 
164
    // the 3 derivatives ( actually this is the gradient )
 
165
    std::vector< double > deriv( 3 * m_dataSet->getGrid()->size() );
 
166
 
 
167
    // the diffusion coeff
 
168
    std::vector< double > coeff( m_dataSet->getGrid()->size() );
 
169
 
 
170
    boost::shared_ptr< WProgress > prog( new WProgress( "Smoothing images", numImages ) );
 
171
    m_progress->addSubProgress( prog );
 
172
 
 
173
    for( std::size_t k = 0; k < numImages; ++k )
 
174
    {
 
175
        for( int i = 0; i < iterations; ++i )
 
176
        {
 
177
            calcDeriv( deriv, smoothed, grid, k, numImages );
 
178
            if( m_iterations->changed() || m_Kcoeff->changed() || m_delta->changed() || m_dataSet != m_input->getData() )
 
179
            {
 
180
                prog->finish();
 
181
                return;
 
182
            }
 
183
 
 
184
            calcCoeff( coeff, deriv, grid );
 
185
            if( m_iterations->changed() || m_Kcoeff->changed() || m_delta->changed() || m_dataSet != m_input->getData() )
 
186
            {
 
187
                prog->finish();
 
188
                return;
 
189
            }
 
190
 
 
191
            diffusion( deriv, coeff, smoothed, grid, k, numImages );
 
192
            if( m_iterations->changed() || m_Kcoeff->changed() || m_delta->changed() || m_dataSet != m_input->getData() )
 
193
            {
 
194
                prog->finish();
 
195
                return;
 
196
            }
 
197
        }
 
198
        ++*prog;
 
199
    }
 
200
 
 
201
    prog->finish();
 
202
 
 
203
    // create dataset and update connector
 
204
    boost::shared_ptr< WValueSet< double > > vs( new WValueSet< double >(
 
205
                                                    m_dataSet->getValueSet()->order(),
 
206
                                                    m_dataSet->getValueSet()->dimension(),
 
207
                                                    smoothed,
 
208
                                                    W_DT_DOUBLE ) );
 
209
    boost::shared_ptr< WDataSetSingle > ds = m_dataSet->clone( vs );
 
210
 
 
211
    m_output->updateData( ds );
 
212
}
 
213
 
 
214
std::size_t WMAnisotropicFiltering::coordsToIndex( boost::shared_ptr< WGridRegular3D > const& grid,
 
215
                                                   std::size_t x, std::size_t y, std::size_t z )
 
216
{
 
217
    return x + y * grid->getNbCoordsX() + z * grid->getNbCoordsX() * grid->getNbCoordsY();
 
218
}
 
219
 
 
220
void WMAnisotropicFiltering::copyData( boost::shared_ptr< std::vector< double > >& smoothed,  // NOLINT non-const ref
 
221
                                       boost::shared_ptr< WGridRegular3D > const& /* grid */ )
 
222
{
 
223
    for( std::size_t k = 0; k < m_dataSet->getValueSet()->rawSize(); ++k )
 
224
    {
 
225
        ( *smoothed )[ k ] = m_dataSet->getValueSet()->getScalarDouble( k );
 
226
    }
 
227
}
 
228
 
 
229
void WMAnisotropicFiltering::calcDeriv( std::vector< double >& deriv,  // NOLINT non-const ref
 
230
                                        boost::shared_ptr< std::vector< double > > const& smoothed,
 
231
                                        boost::shared_ptr< WGridRegular3D > const& grid, std::size_t image, std::size_t numImages )
 
232
{
 
233
    std::size_t s[] = { grid->getNbCoordsX(), grid->getNbCoordsY(), grid->getNbCoordsZ() };
 
234
    double d[] = { fabs( grid->getOffsetX() ), fabs( grid->getOffsetY() ), fabs( grid->getOffsetZ() ) };
 
235
 
 
236
    for( std::size_t x = 0; x < grid->getNbCoordsX(); ++x )
 
237
    {
 
238
        for( std::size_t y = 0; y < grid->getNbCoordsY(); ++y )
 
239
        {
 
240
            for( std::size_t z = 0; z < grid->getNbCoordsZ(); ++z )
 
241
            {
 
242
                double const& x1 = ( *smoothed )[ numImages * coordsToIndex( grid, ( x + 1 ) % s[ 0 ], y, z ) + image ];
 
243
                double const& x2 = ( *smoothed )[ numImages * coordsToIndex( grid, ( x + s[ 0 ] - 1 ) % s[ 0 ], y, z ) + image ];
 
244
                deriv[ 3 * coordsToIndex( grid, x, y, z ) + 0 ] = ( x1 - x2 ) / ( 2.0 * d[ 0 ] );
 
245
 
 
246
                double const& y1 = ( *smoothed )[ numImages * coordsToIndex( grid, x, ( y + 1 ) % s[ 1 ], z ) + image ];
 
247
                double const& y2 = ( *smoothed )[ numImages * coordsToIndex( grid, x, ( y + s[ 1 ] - 1 ) % s[ 1 ], z ) + image ];
 
248
                deriv[ 3 * coordsToIndex( grid, x, y, z ) + 1 ] = ( y1 - y2 ) / ( 2.0 * d[ 1 ] );
 
249
 
 
250
                double const& z1 = ( *smoothed )[ numImages * coordsToIndex( grid, x, y, ( z + 1 ) % s[ 2 ] ) + image ];
 
251
                double const& z2 = ( *smoothed )[ numImages * coordsToIndex( grid, x, y, ( z + s[ 2 ] - 1 ) % s[ 2 ] ) + image ];
 
252
                deriv[ 3 * coordsToIndex( grid, x, y, z ) + 2 ] = ( z1 - z2 ) / ( 2.0 * d[ 2 ] );
 
253
            }
 
254
        }
 
255
    }
 
256
}
 
257
 
 
258
void WMAnisotropicFiltering::calcCoeff( std::vector< double >& coeff,  // NOLINT non-const ref
 
259
                                        std::vector< double > const& deriv,
 
260
                                        boost::shared_ptr< WGridRegular3D > const& grid )
 
261
{
 
262
    for( std::size_t x = 0; x < grid->getNbCoordsX(); ++x )
 
263
    {
 
264
        for( std::size_t y = 0; y < grid->getNbCoordsY(); ++y )
 
265
        {
 
266
            for( std::size_t z = 0; z < grid->getNbCoordsZ(); ++z )
 
267
            {
 
268
                // coeff = exp( -sqr( |I|/K ) )
 
269
                double gradIAbsSquared = deriv[ 3 * coordsToIndex( grid, x, y, z ) + 0 ]
 
270
                                       * deriv[ 3 * coordsToIndex( grid, x, y, z ) + 0 ]
 
271
                                       + deriv[ 3 * coordsToIndex( grid, x, y, z ) + 1 ]
 
272
                                       * deriv[ 3 * coordsToIndex( grid, x, y, z ) + 1 ]
 
273
                                       + deriv[ 3 * coordsToIndex( grid, x, y, z ) + 2 ]
 
274
                                       * deriv[ 3 * coordsToIndex( grid, x, y, z ) + 2 ];
 
275
                coeff[ coordsToIndex( grid, x, y, z ) ] = 1.0 / exp( gradIAbsSquared / ( m_k * m_k ) );
 
276
            }
 
277
        }
 
278
    }
 
279
}
 
280
 
 
281
void WMAnisotropicFiltering::diffusion( std::vector< double > const& deriv, std::vector< double > const& coeff,
 
282
                                        boost::shared_ptr< std::vector< double > >& smoothed,  // NOLINT non-const ref
 
283
                                        boost::shared_ptr< WGridRegular3D > const& grid, std::size_t image, std::size_t numImages )
 
284
{
 
285
    std::size_t s[] = { grid->getNbCoordsX(), grid->getNbCoordsY(), grid->getNbCoordsZ() };
 
286
    double d[] = { fabs( grid->getOffsetX() ), fabs( grid->getOffsetY() ), fabs( grid->getOffsetZ() ) };
 
287
 
 
288
    for( std::size_t x = 0; x < grid->getNbCoordsX(); ++x )
 
289
    {
 
290
        for( std::size_t y = 0; y < grid->getNbCoordsY(); ++y )
 
291
        {
 
292
            for( std::size_t z = 0; z < grid->getNbCoordsZ(); ++z )
 
293
            {
 
294
                // first deriv of the image intensity
 
295
                double const& dIx = deriv[ 3 * coordsToIndex( grid, x, y, z ) + 0 ];
 
296
                double const& dIy = deriv[ 3 * coordsToIndex( grid, x, y, z ) + 1 ];
 
297
                double const& dIz = deriv[ 3 * coordsToIndex( grid, x, y, z ) + 2 ];
 
298
 
 
299
                // first deriv of the diffusion coeff
 
300
                double dcx = ( coeff[ coordsToIndex( grid, ( x + 1 ) % s[ 0 ], y, z ) ]
 
301
                             - coeff[ coordsToIndex( grid, ( x + s[ 0 ] - 1 ) % s[ 0 ], y, z ) ] )
 
302
                              / ( 2.0 * d[ 0 ] );
 
303
                double dcy = ( coeff[ coordsToIndex( grid, x, ( y + 1 ) % s[ 1 ], z ) ]
 
304
                             - coeff[ coordsToIndex( grid, x, ( y + s[ 1 ] - 1 ) % s[ 1 ], z ) ] )
 
305
                              / ( 2.0 * d[ 1 ] );
 
306
                double dcz = ( coeff[ coordsToIndex( grid, x, y, ( z + 1 ) % s[ 2 ] ) ]
 
307
                             - coeff[ coordsToIndex( grid, x, y, ( z + s[ 2 ] - 1 ) % s[ 2 ] ) ] )
 
308
                              / ( 2.0 * d[ 2 ] );
 
309
 
 
310
                // diffusion coeff
 
311
                double const& c = coeff[ coordsToIndex( grid, x, y, z ) ];
 
312
 
 
313
                // 2nd derivatives in x, y, and z of the image intensity
 
314
                double ddIxx = ( deriv[ 3 * coordsToIndex( grid, ( x + 1 ) % s[ 0 ], y, z ) + 0 ]
 
315
                               - deriv[ 3 * coordsToIndex( grid, ( x + s[ 0 ] - 1 ) % s[ 0 ], y, z ) + 0 ] )
 
316
                                / ( 2 * d[ 0 ] );
 
317
                double ddIyy = ( deriv[ 3 * coordsToIndex( grid, x, ( y + 1 ) % s[ 1 ], z ) + 1 ]
 
318
                               - deriv[ 3 * coordsToIndex( grid, x, ( y + s[ 1 ] - 1 ) % s[ 1 ], z ) + 1 ] )
 
319
                                / ( 2 * d[ 1 ] );
 
320
                double ddIzz = ( deriv[ 3 * coordsToIndex( grid, x, y, ( z + 1 ) % s[ 2 ] ) + 2 ]
 
321
                               - deriv[ 3 * coordsToIndex( grid, x, y, ( z + s[ 2 ] - 1 ) % s[ 2 ] ) + 2 ] )
 
322
                                / ( 2 * d[ 2 ] );
 
323
 
 
324
                // d * ( grad I * grad c + c * grad grad I )
 
325
                ( *smoothed )[ numImages * coordsToIndex( grid, x, y, z ) + image ] += m_d * ( dIx * dcx + dIy * dcy + dIz * dcz
 
326
                                                                                             + c * ( ddIxx + ddIyy + ddIzz ) );
 
327
            }
 
328
        }
 
329
    }
 
330
}