← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/modules/vectorOperator/WMVectorOperator.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/vectorOperator/WMVectorOperator.cpp
 
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 <stdint.h>
 
26
 
 
27
#include <cmath>
 
28
#include <fstream>
 
29
#include <iostream>
 
30
#include <string>
 
31
#include <vector>
 
32
 
 
33
#include <boost/variant.hpp>
 
34
 
 
35
#include "core/common/math/linearAlgebra/WLinearAlgebra.h"
 
36
#include "core/common/WAssert.h"
 
37
#include "core/common/WProgress.h"
 
38
#include "core/common/WStringUtils.h"
 
39
#include "core/common/WTypeTraits.h"
 
40
#include "core/common/exceptions/WTypeMismatch.h"
 
41
#include "core/dataHandler/WGridRegular3D.h"
 
42
#include "core/dataHandler/WDataHandlerEnums.h"
 
43
#include "core/dataHandler/WDataHandler.h"
 
44
#include "core/dataHandler/exceptions/WDHValueSetMismatch.h"
 
45
#include "core/kernel/WKernel.h"
 
46
#include "WMVectorOperator.xpm"
 
47
#include "WMVectorOperator.h"
 
48
 
 
49
// This line is needed by the module loader to actually find your module.
 
50
W_LOADABLE_MODULE( WMVectorOperator )
 
51
 
 
52
WMVectorOperator::WMVectorOperator() :
 
53
    WModule()
 
54
{
 
55
    // initialize
 
56
}
 
57
 
 
58
WMVectorOperator::~WMVectorOperator()
 
59
{
 
60
    // cleanup
 
61
    removeConnectors();
 
62
}
 
63
 
 
64
boost::shared_ptr< WModule > WMVectorOperator::factory() const
 
65
{
 
66
    return boost::shared_ptr< WModule >( new WMVectorOperator() );
 
67
}
 
68
 
 
69
const char** WMVectorOperator::getXPMIcon() const
 
70
{
 
71
    return WMVectorOperator_xpm;
 
72
}
 
73
 
 
74
const std::string WMVectorOperator::getName() const
 
75
{
 
76
    return "Vector Operator";
 
77
}
 
78
 
 
79
const std::string WMVectorOperator::getDescription() const
 
80
{
 
81
    return "Applies an selected operator on a specified vector field.";
 
82
}
 
83
 
 
84
void WMVectorOperator::connectors()
 
85
{
 
86
    m_inputA = WModuleInputData< WDataSetVector >::createAndAdd( shared_from_this(), "operandA", "First operand of operation.." );
 
87
 
 
88
    m_output = WModuleOutputData< WDataSetScalar >::createAndAdd( shared_from_this(), "result", "Result of voxel-wise operation( A, B )." );
 
89
 
 
90
    // call WModules initialization
 
91
    WModule::connectors();
 
92
}
 
93
 
 
94
void WMVectorOperator::properties()
 
95
{
 
96
    m_propCondition = boost::shared_ptr< WCondition >( new WCondition() );
 
97
 
 
98
    // create a list of operations here
 
99
    m_operations = boost::shared_ptr< WItemSelection >( new WItemSelection() );
 
100
    m_operations->addItem( "Length", "Length of the vector." );          // NOTE: you can add XPM images here.
 
101
    m_operations->addItem( "Curvature", "Curvature at each voxel." );
 
102
 
 
103
    m_opSelection = m_properties->addProperty( "Operation", "The operation to apply on A and B.", m_operations->getSelectorFirst(),
 
104
                                               m_propCondition );
 
105
    WPropertyHelper::PC_SELECTONLYONE::addTo( m_opSelection );
 
106
    WPropertyHelper::PC_NOTEMPTY::addTo( m_opSelection );
 
107
 
 
108
    WModule::properties();
 
109
}
 
110
 
 
111
template< typename T >
 
112
T opLength( const WVector3d& vec,
 
113
            const WVector3d& /*mx*/, const WVector3d& /*px*/,
 
114
            const WVector3d& /*my*/, const WVector3d& /*py*/,
 
115
            const WVector3d& /*mz*/, const WVector3d& /*pz*/ )
 
116
{
 
117
    return length( vec );
 
118
}
 
119
 
 
120
template< typename T >
 
121
T opCurvature( const WVector3d& vec,
 
122
               const WVector3d& mx, const WVector3d& px,
 
123
               const WVector3d& my, const WVector3d& py,
 
124
               const WVector3d& mz, const WVector3d& pz )
 
125
{
 
126
    // get partial differentiation in x direction:
 
127
    WVector3d dx = ( px - mx ) / 2.0;    // NOTE: step size h=2
 
128
    WVector3d dy = ( py - my ) / 2.0;
 
129
    WVector3d dz = ( pz - mz ) / 2.0;
 
130
 
 
131
    // get second derivative of tangent curve
 
132
    WVector3d L2 = ( vec.x() * dx ) + ( vec.y() * dy ) + ( vec.z() * dz );
 
133
 
 
134
    WVector3d::ValueType l = length( vec );
 
135
    return length( cross( vec, L2 ) ) / ( l * l * l );
 
136
}
 
137
 
 
138
/**
 
139
 * Get the ID of a voxel with the given coordinates.
 
140
 *
 
141
 * \param xDim number of voxels in x direction
 
142
 * \param yDim number of voxels in y direction
 
143
 * \param zDim number of voxels in z direction
 
144
 * \param x x coordinate of point to get index for
 
145
 * \param y y coordinate of point to get index for
 
146
 * \param z z coordinate of point to get index for
 
147
 * \param offset which coordinate? 0 = x, 3 = w
 
148
 * \param elements number of elements
 
149
 *
 
150
 * \return the index
 
151
 */
 
152
size_t getId( size_t xDim, size_t yDim, size_t /*zDim*/, size_t x, size_t y, size_t z, size_t offset = 0, size_t elements = 1 )
 
153
{
 
154
    return offset + ( elements * ( z * xDim * yDim + y * xDim + x ) );
 
155
}
 
156
 
 
157
/**
 
158
 * Visitor for discriminating the type of the first valueset.
 
159
 */
 
160
class VisitorVSetA: public boost::static_visitor< boost::shared_ptr< WValueSetBase > >
 
161
{
 
162
public:
 
163
    /**
 
164
     * Create visitor instance.
 
165
     *
 
166
     * \param opIdx The operator index.
 
167
     * \param grid the underlying grid
 
168
     */
 
169
    VisitorVSetA( boost::shared_ptr< WGridRegular3D > grid, size_t opIdx = 0 ):
 
170
        boost::static_visitor< result_type >(),
 
171
        m_grid( grid ),
 
172
        m_opIdx( opIdx )
 
173
    {
 
174
    }
 
175
 
 
176
    /**
 
177
     * Called by boost::varying during static visiting.
 
178
     *
 
179
     * \tparam T the real integral type of the first value set.
 
180
     * \param vsetA the first valueset currently visited.
 
181
     *
 
182
     * \return the result from the operation
 
183
     */
 
184
    template < typename T >
 
185
    result_type operator()( const WValueSet< T >* const& vsetA ) const             // NOLINT
 
186
    {
 
187
        // get some info
 
188
        std::vector< T > data;
 
189
        data.resize( vsetA->size() );
 
190
 
 
191
        // discriminate the right operation with the correct type. It would be nicer to use some kind of strategy pattern here, but the template
 
192
        // character of the operators forbids it as template methods can't be virtual. Besides this, at some point in the module main the
 
193
        // selector needs to be queried and its index mapped to a pointer. This is what we do here.
 
194
        boost::function< T( const WVector3d&,
 
195
                            const WVector3d&, const WVector3d&,
 
196
                            const WVector3d&, const WVector3d&,
 
197
                            const WVector3d&, const WVector3d& ) > op;
 
198
        switch ( m_opIdx )
 
199
        {
 
200
            case 1:
 
201
                op = &opCurvature< T >;
 
202
                break;
 
203
            case 0:
 
204
            default:
 
205
                op = &opLength< T >;
 
206
                break;
 
207
        }
 
208
 
 
209
        // some info needed for indexing the vector components
 
210
        size_t nX = m_grid->getNbCoordsX();
 
211
        size_t nY = m_grid->getNbCoordsY();
 
212
        size_t nZ = m_grid->getNbCoordsZ();
 
213
 
 
214
        // apply op to each value
 
215
        // iterate field
 
216
        for( size_t z = 1; z < nZ - 1; z++ )
 
217
        {
 
218
            for( size_t y = 1; y < nY - 1; y++ )
 
219
            {
 
220
                for( size_t x = 1; x < nX - 1; x++ )
 
221
                {
 
222
                    // this is ugly. We'll fix this crap with the upcoming new data handler
 
223
                    size_t idx = getId( nX, nY, nZ, x, y, z );
 
224
                    WVector3d vec = vsetA->getVector3D( idx );
 
225
 
 
226
                    // also get the neighbours
 
227
                    WVector3d mx = vsetA->getVector3D( getId( nX, nY, nZ, x - 1, y, z ) );
 
228
                    WVector3d px = vsetA->getVector3D( getId( nX, nY, nZ, x + 1, y, z ) );
 
229
                    WVector3d my = vsetA->getVector3D( getId( nX, nY, nZ, x, y - 1, z ) );
 
230
                    WVector3d py = vsetA->getVector3D( getId( nX, nY, nZ, x, y + 1, z ) );
 
231
                    WVector3d mz = vsetA->getVector3D( getId( nX, nY, nZ, x, y, z - 1 ) );
 
232
                    WVector3d pz = vsetA->getVector3D( getId( nX, nY, nZ, x, y, z + 1 ) );
 
233
 
 
234
                    data[ idx ] = op( vec, mx, px, my, py, mz, pz );
 
235
                }
 
236
            }
 
237
        }
 
238
 
 
239
        // init the border values
 
240
        for( size_t z = 0; z < 2; z++ )
 
241
        {
 
242
            for( size_t y = 1; y < 2; y++ )
 
243
            {
 
244
                for( size_t x = 1; x < 2; x++ )
 
245
                {
 
246
                    size_t idx = getId( nX, nY, nZ, x * ( nX - 1 ), y * ( nY - 1 ), z * ( nZ - 1 ) );
 
247
                    data[ idx ] = 0.0;
 
248
                }
 
249
            }
 
250
        }
 
251
 
 
252
        // create result value set
 
253
        return boost::shared_ptr< WValueSet< T > >( new WValueSet< T >( 0,
 
254
                                                                        1,
 
255
                                                                        boost::shared_ptr< std::vector< T > >(
 
256
                                                                            new std::vector< T >( data ) ),
 
257
                                                                        DataType< T >::type ) );
 
258
    }
 
259
 
 
260
    /**
 
261
     * The underlying grid.
 
262
     */
 
263
    boost::shared_ptr< WGridRegular3D > m_grid;
 
264
 
 
265
    /**
 
266
     * The operator index.
 
267
     */
 
268
    size_t m_opIdx;
 
269
};
 
270
 
 
271
void WMVectorOperator::moduleMain()
 
272
{
 
273
    // let the main loop awake if the data changes or the properties changed.
 
274
    m_moduleState.setResetable( true, true );
 
275
    m_moduleState.add( m_inputA->getDataChangedCondition() );
 
276
    m_moduleState.add( m_propCondition );
 
277
 
 
278
    // signal ready state
 
279
    ready();
 
280
 
 
281
    // loop until the module container requests the module to quit
 
282
    while( !m_shutdownFlag() )
 
283
    {
 
284
        // Now, the moduleState variable comes into play. The module can wait for the condition, which gets fired whenever the input receives data
 
285
        // or an property changes. The main loop now waits until something happens.
 
286
        debugLog() << "Waiting ...";
 
287
        m_moduleState.wait();
 
288
 
 
289
        // woke up since the module is requested to finish
 
290
        if( m_shutdownFlag() )
 
291
        {
 
292
            break;
 
293
        }
 
294
 
 
295
        // has the data changed?
 
296
        if( m_opSelection->changed() || m_inputA->handledUpdate() )
 
297
        {
 
298
            boost::shared_ptr< WDataSetVector > dataSetA = m_inputA->getData();
 
299
 
 
300
            WItemSelector s = m_opSelection->get( true );
 
301
 
 
302
            // valid data?
 
303
            if( dataSetA )
 
304
            {
 
305
                boost::shared_ptr< WValueSetBase > valueSetA = dataSetA->getValueSet();
 
306
 
 
307
                // use a custom progress combiner
 
308
                boost::shared_ptr< WProgress > prog = boost::shared_ptr< WProgress >(
 
309
                    new WProgress( "Applying operator on data" ) );
 
310
                m_progress->addSubProgress( prog );
 
311
 
 
312
                // apply the operation to each voxel
 
313
                debugLog() << "Processing ...";
 
314
                boost::shared_ptr< WValueSetBase > newValueSet = valueSetA->applyFunction( VisitorVSetA(
 
315
                    boost::shared_dynamic_cast< WGridRegular3D >( dataSetA->getGrid() ), s )
 
316
                );
 
317
 
 
318
                // Create the new dataset and export it
 
319
                m_output->updateData( boost::shared_ptr<WDataSetScalar>( new WDataSetScalar( newValueSet, m_inputA->getData()->getGrid() ) ) );
 
320
 
 
321
                // done
 
322
                prog->finish();
 
323
                m_progress->removeSubProgress( prog );
 
324
            }
 
325
            else
 
326
            {
 
327
                debugLog() << "Resetting output.";
 
328
                m_output->reset();
 
329
            }
 
330
        }
 
331
    }
 
332
}
 
333