1
//---------------------------------------------------------------------------
3
// Project: OpenWalnut ( http://www.openwalnut.org )
5
// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
6
// For more information see http://www.openwalnut.org/copying
8
// This file is part of OpenWalnut.
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.
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.
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/>.
23
//---------------------------------------------------------------------------
33
#include <boost/variant.hpp>
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"
49
// This line is needed by the module loader to actually find your module.
50
W_LOADABLE_MODULE( WMVectorOperator )
52
WMVectorOperator::WMVectorOperator() :
58
WMVectorOperator::~WMVectorOperator()
64
boost::shared_ptr< WModule > WMVectorOperator::factory() const
66
return boost::shared_ptr< WModule >( new WMVectorOperator() );
69
const char** WMVectorOperator::getXPMIcon() const
71
return WMVectorOperator_xpm;
74
const std::string WMVectorOperator::getName() const
76
return "Vector Operator";
79
const std::string WMVectorOperator::getDescription() const
81
return "Applies an selected operator on a specified vector field.";
84
void WMVectorOperator::connectors()
86
m_inputA = WModuleInputData< WDataSetVector >::createAndAdd( shared_from_this(), "operandA", "First operand of operation.." );
88
m_output = WModuleOutputData< WDataSetScalar >::createAndAdd( shared_from_this(), "result", "Result of voxel-wise operation( A, B )." );
90
// call WModules initialization
91
WModule::connectors();
94
void WMVectorOperator::properties()
96
m_propCondition = boost::shared_ptr< WCondition >( new WCondition() );
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." );
103
m_opSelection = m_properties->addProperty( "Operation", "The operation to apply on A and B.", m_operations->getSelectorFirst(),
105
WPropertyHelper::PC_SELECTONLYONE::addTo( m_opSelection );
106
WPropertyHelper::PC_NOTEMPTY::addTo( m_opSelection );
108
WModule::properties();
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*/ )
117
return length( vec );
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 )
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;
131
// get second derivative of tangent curve
132
WVector3d L2 = ( vec.x() * dx ) + ( vec.y() * dy ) + ( vec.z() * dz );
134
WVector3d::ValueType l = length( vec );
135
return length( cross( vec, L2 ) ) / ( l * l * l );
139
* Get the ID of a voxel with the given coordinates.
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
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 )
154
return offset + ( elements * ( z * xDim * yDim + y * xDim + x ) );
158
* Visitor for discriminating the type of the first valueset.
160
class VisitorVSetA: public boost::static_visitor< boost::shared_ptr< WValueSetBase > >
164
* Create visitor instance.
166
* \param opIdx The operator index.
167
* \param grid the underlying grid
169
VisitorVSetA( boost::shared_ptr< WGridRegular3D > grid, size_t opIdx = 0 ):
170
boost::static_visitor< result_type >(),
177
* Called by boost::varying during static visiting.
179
* \tparam T the real integral type of the first value set.
180
* \param vsetA the first valueset currently visited.
182
* \return the result from the operation
184
template < typename T >
185
result_type operator()( const WValueSet< T >* const& vsetA ) const // NOLINT
188
std::vector< T > data;
189
data.resize( vsetA->size() );
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;
201
op = &opCurvature< T >;
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();
214
// apply op to each value
216
for( size_t z = 1; z < nZ - 1; z++ )
218
for( size_t y = 1; y < nY - 1; y++ )
220
for( size_t x = 1; x < nX - 1; x++ )
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 );
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 ) );
234
data[ idx ] = op( vec, mx, px, my, py, mz, pz );
239
// init the border values
240
for( size_t z = 0; z < 2; z++ )
242
for( size_t y = 1; y < 2; y++ )
244
for( size_t x = 1; x < 2; x++ )
246
size_t idx = getId( nX, nY, nZ, x * ( nX - 1 ), y * ( nY - 1 ), z * ( nZ - 1 ) );
252
// create result value set
253
return boost::shared_ptr< WValueSet< T > >( new WValueSet< T >( 0,
255
boost::shared_ptr< std::vector< T > >(
256
new std::vector< T >( data ) ),
257
DataType< T >::type ) );
261
* The underlying grid.
263
boost::shared_ptr< WGridRegular3D > m_grid;
266
* The operator index.
271
void WMVectorOperator::moduleMain()
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 );
278
// signal ready state
281
// loop until the module container requests the module to quit
282
while( !m_shutdownFlag() )
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();
289
// woke up since the module is requested to finish
290
if( m_shutdownFlag() )
295
// has the data changed?
296
if( m_opSelection->changed() || m_inputA->handledUpdate() )
298
boost::shared_ptr< WDataSetVector > dataSetA = m_inputA->getData();
300
WItemSelector s = m_opSelection->get( true );
305
boost::shared_ptr< WValueSetBase > valueSetA = dataSetA->getValueSet();
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 );
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 )
318
// Create the new dataset and export it
319
m_output->updateData( boost::shared_ptr<WDataSetScalar>( new WDataSetScalar( newValueSet, m_inputA->getData()->getGrid() ) ) );
323
m_progress->removeSubProgress( prog );
327
debugLog() << "Resetting output.";