[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

multi_blockwise.hxx
1/************************************************************************/
2/* */
3/* Copyright 2015 by Thorsten Beier */
4/* */
5/* This file is part of the VIGRA computer vision library. */
6/* The VIGRA Website is */
7/* http://hci.iwr.uni-heidelberg.de/vigra/ */
8/* Please direct questions, bug reports, and contributions to */
9/* ullrich.koethe@iwr.uni-heidelberg.de or */
10/* vigra@informatik.uni-hamburg.de */
11/* */
12/* Permission is hereby granted, free of charge, to any person */
13/* obtaining a copy of this software and associated documentation */
14/* files (the "Software"), to deal in the Software without */
15/* restriction, including without limitation the rights to use, */
16/* copy, modify, merge, publish, distribute, sublicense, and/or */
17/* sell copies of the Software, and to permit persons to whom the */
18/* Software is furnished to do so, subject to the following */
19/* conditions: */
20/* */
21/* The above copyright notice and this permission notice shall be */
22/* included in all copies or substantial portions of the */
23/* Software. */
24/* */
25/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32/* OTHER DEALINGS IN THE SOFTWARE. */
33/* */
34/************************************************************************/
35
36
37#ifndef VIGRA_MULTI_BLOCKWISE_HXX
38#define VIGRA_MULTI_BLOCKWISE_HXX
39
40#include <cmath>
41#include <vector>
42#include "multi_blocking.hxx"
43#include "multi_convolution.hxx"
44#include "multi_tensorutilities.hxx"
45#include "threadpool.hxx"
46#include "array_vector.hxx"
47
48 namespace vigra{
49
50 /** Option base class for blockwise algorithms.
51
52 Attaches blockshape to ParallelOptions.
53 */
54 class BlockwiseOptions
55: public ParallelOptions
56{
57 public:
58 typedef ArrayVector<MultiArrayIndex> Shape;
59
60 BlockwiseOptions()
61 : ParallelOptions()
62 , blockShape_()
63 {}
64
65 /** Retrieve block shape as a std::vector.
66
67 If the returned vector is empty, a default block shape should be used.
68 If the returned vector has length 1, square blocks of size
69 <tt>getBlockShape()[0]</tt> should be used.
70 */
71 Shape const & getBlockShape() const
72 {
73 return blockShape_;
74 }
75
76 // for Python bindings
77 Shape readBlockShape() const
78 {
79 return blockShape_;
80 }
81
82 /** Retrieve block shape as a fixed-size vector.
83
84 Default shape specifications are appropriately expanded.
85 An exception is raised if the stored shape's length is
86 incompatible with dimension <tt>N</tt>.
87 */
88 template <int N>
89 TinyVector<MultiArrayIndex, N> getBlockShapeN() const
90 {
91 if(blockShape_.size() > 1)
92 {
93 vigra_precondition(blockShape_.size() == (size_t)N,
94 "BlockwiseOptions::getBlockShapeN(): dimension mismatch between N and stored block shape.");
95 return TinyVector<MultiArrayIndex, N>(blockShape_.data());
96 }
97 else if(blockShape_.size() == 1)
98 {
99 return TinyVector<MultiArrayIndex, N>(blockShape_[0]);
100 }
101 else
102 {
103 return detail::ChunkShape<N>::defaultShape();
104 }
105 }
106
107 /** Specify block shape as a std::vector of appropriate length.
108 If <tt>blockShape.size() == 0</tt>, the default shape is used.
109 If <tt>blockShape.size() == 1</tt>, square blocks of size
110 <tt>blockShape[0]</tt> are used.
111
112 Default: Use square blocks with side length <tt>VIGRA_DEFAULT_BLOCK_SHAPE</tt>.
113 */
114 BlockwiseOptions & blockShape(const Shape & blockShape){
115 blockShape_ = blockShape;
116 return *this;
117 }
118
119 // for Python bindings
120 void setBlockShape(const Shape & blockShape){
121 blockShape_ = blockShape;
122 }
123
124 /** Specify block shape by a fixed-size shape object.
125 */
126 template <class T, int N>
127 BlockwiseOptions & blockShape(const TinyVector<T, N> & blockShape){
128 Shape(blockShape.begin(), blockShape.end()).swap(blockShape_);
129 return *this;
130 }
131
132 /** Specify square block shape by its side length.
133 */
134 BlockwiseOptions & blockShape(MultiArrayIndex blockShape){
135 Shape(1, blockShape).swap(blockShape_);
136 return *this;
137 }
138
139 BlockwiseOptions & numThreads(const int n)
140 {
141 ParallelOptions::numThreads(n);
142 return *this;
143 }
144
145 void setNumThreads(const int n)
146 {
147 ParallelOptions::numThreads(n);
148 }
149
150 private:
151 Shape blockShape_;
152};
153
154 /** Option class for blockwise convolution algorithms.
155
156 Simply derives from \ref vigra::BlockwiseOptions and
157 \ref vigra::ConvolutionOptions to join their capabilities.
158 */
159 template<unsigned int N>
160 class BlockwiseConvolutionOptions
161: public BlockwiseOptions
162, public ConvolutionOptions<N>{
163 public:
164 BlockwiseConvolutionOptions()
165 : BlockwiseOptions(),
166 ConvolutionOptions<N>()
167 {}
168};
169
170
171 namespace blockwise{
172
173 /**
174 helper function to create blockwise parallel filters.
175 This implementation should be used if the filter functor
176 does not support the ROI/sub array options.
177 */
178 template<
179 unsigned int DIM,
180 class T_IN, class ST_IN,
181 class T_OUT, class ST_OUT,
182 class FILTER_FUNCTOR,
183 class C
184 >
185 void blockwiseCallerNoRoiApi(
186 const vigra::MultiArrayView<DIM, T_IN, ST_IN > & source,
187 const vigra::MultiArrayView<DIM, T_OUT, ST_OUT> & dest,
188 FILTER_FUNCTOR & functor,
189 const vigra::MultiBlocking<DIM, C> & blocking,
190 const typename vigra::MultiBlocking<DIM, C>::Shape & borderWidth,
191 const BlockwiseConvolutionOptions<DIM> & options
192 ){
193
194 typedef typename MultiBlocking<DIM, C>::BlockWithBorder BlockWithBorder;
195
196 auto beginIter = blocking.blockWithBorderBegin(borderWidth);
197 auto endIter = blocking.blockWithBorderEnd(borderWidth);
198
199 parallel_foreach(options.getNumThreads(),
200 beginIter, endIter,
201 [&](const int /*threadId*/, const BlockWithBorder bwb)
202 {
203 // get the input of the block as a view
204 vigra::MultiArrayView<DIM, T_IN, ST_IN> sourceSub = source.subarray(bwb.border().begin(),
205 bwb.border().end());
206 // get the output as NEW allocated array
207 vigra::MultiArray<DIM, T_OUT> destSub(sourceSub.shape());
208 // call the functor
209 functor(sourceSub, destSub);
210 // write the core global out
211 vigra::MultiArrayView<DIM, T_OUT, ST_OUT> destSubCore = destSub.subarray(bwb.localCore().begin(),
212 bwb.localCore().end());
213 // write the core global out
214 dest.subarray(bwb.core().begin()-blocking.roiBegin(),
215 bwb.core().end() -blocking.roiBegin() ) = destSubCore;
216 },
217 blocking.numBlocks()
218 );
219
220 }
221
222 /**
223 helper function to create blockwise parallel filters.
224 This implementation should be used if the filter functor
225 does support the ROI/sub array options.
226 */
227 template<
228 unsigned int DIM,
229 class T_IN, class ST_IN,
230 class T_OUT, class ST_OUT,
231 class FILTER_FUNCTOR,
232 class C
233 >
234 void blockwiseCaller(
235 const vigra::MultiArrayView<DIM, T_IN, ST_IN > & source,
236 const vigra::MultiArrayView<DIM, T_OUT, ST_OUT> & dest,
237 FILTER_FUNCTOR & functor,
238 const vigra::MultiBlocking<DIM, C> & blocking,
239 const typename vigra::MultiBlocking<DIM, C>::Shape & borderWidth,
240 const BlockwiseConvolutionOptions<DIM> & options
241 ){
242
243 typedef typename MultiBlocking<DIM, C>::BlockWithBorder BlockWithBorder;
244 //typedef typename MultiBlocking<DIM, C>::BlockWithBorderIter BlockWithBorderIter;
245 typedef typename MultiBlocking<DIM, C>::Block Block;
246
247
248 auto beginIter = blocking.blockWithBorderBegin(borderWidth);
249 auto endIter = blocking.blockWithBorderEnd(borderWidth);
250
251 parallel_foreach(options.getNumThreads(),
252 beginIter, endIter,
253 [&](const int /*threadId*/, const BlockWithBorder bwb)
254 {
255 // get the input of the block as a view
256 vigra::MultiArrayView<DIM, T_IN, ST_IN> sourceSub = source.subarray(bwb.border().begin(),
257 bwb.border().end());
258 // get the output of the blocks core as a view
259 vigra::MultiArrayView<DIM, T_OUT, ST_OUT> destCore = dest.subarray(bwb.core().begin(),
260 bwb.core().end());
261 const Block localCore = bwb.localCore();
262 // call the functor
263 functor(sourceSub, destCore, localCore.begin(), localCore.end());
264 },
265 blocking.numBlocks()
266 );
267
268
269 }
270
271 #define CONVOLUTION_FUNCTOR(FUNCTOR_NAME, FUNCTION_NAME) \
272 template<unsigned int DIM> \
273 class FUNCTOR_NAME{ \
274 public: \
275 typedef ConvolutionOptions<DIM> ConvOpt; \
276 FUNCTOR_NAME(const ConvOpt & convOpt) \
277 : sharedOpt_(convOpt){} \
278 template<class S, class D> \
279 void operator()(const S & s, D & d)const{ \
280 FUNCTION_NAME(s, d, sharedOpt_); \
281 } \
282 template<class S, class D,class SHAPE> \
283 void operator()(const S & s, D & d, const SHAPE & roiBegin, const SHAPE & roiEnd){ \
284 ConvOpt localOpt(sharedOpt_); \
285 localOpt.subarray(roiBegin, roiEnd); \
286 FUNCTION_NAME(s, d, localOpt); \
287 } \
288 private: \
289 ConvOpt sharedOpt_; \
290 };
291
292
293 CONVOLUTION_FUNCTOR(GaussianSmoothFunctor, vigra::gaussianSmoothMultiArray);
294 CONVOLUTION_FUNCTOR(GaussianGradientFunctor, vigra::gaussianGradientMultiArray);
295 CONVOLUTION_FUNCTOR(SymmetricGradientFunctor, vigra::symmetricGradientMultiArray);
296 CONVOLUTION_FUNCTOR(GaussianDivergenceFunctor, vigra::gaussianDivergenceMultiArray);
297 CONVOLUTION_FUNCTOR(HessianOfGaussianFunctor, vigra::hessianOfGaussianMultiArray);
298 CONVOLUTION_FUNCTOR(LaplacianOfGaussianFunctor, vigra::laplacianOfGaussianMultiArray);
299 CONVOLUTION_FUNCTOR(GaussianGradientMagnitudeFunctor, vigra::gaussianGradientMagnitude);
300 CONVOLUTION_FUNCTOR(StructureTensorFunctor, vigra::structureTensorMultiArray);
301
302 #undef CONVOLUTION_FUNCTOR
303
304 template<unsigned int DIM>
305 class HessianOfGaussianEigenvaluesFunctor{
306 public:
307 typedef ConvolutionOptions<DIM> ConvOpt;
308 HessianOfGaussianEigenvaluesFunctor(const ConvOpt & convOpt)
309 : sharedOpt_(convOpt){}
310 template<class S, class D>
311 void operator()(const S & s, D & d)const{
312 typedef typename vigra::NumericTraits<typename S::value_type>::RealPromote RealType;
313 vigra::MultiArray<DIM, TinyVector<RealType, int(DIM*(DIM+1)/2)> > hessianOfGaussianRes(d.shape());
314 vigra::hessianOfGaussianMultiArray(s, hessianOfGaussianRes, sharedOpt_);
315 vigra::tensorEigenvaluesMultiArray(hessianOfGaussianRes, d);
316 }
317 template<class S, class D,class SHAPE>
318 void operator()(const S & s, D & d, const SHAPE & roiBegin, const SHAPE & roiEnd){
319 typedef typename vigra::NumericTraits<typename S::value_type>::RealPromote RealType;
320 vigra::MultiArray<DIM, TinyVector<RealType, int(DIM*(DIM+1)/2)> > hessianOfGaussianRes(roiEnd-roiBegin);
321 ConvOpt localOpt(sharedOpt_);
322 localOpt.subarray(roiBegin, roiEnd);
323 vigra::hessianOfGaussianMultiArray(s, hessianOfGaussianRes, localOpt);
324 vigra::tensorEigenvaluesMultiArray(hessianOfGaussianRes, d);
325 }
326 private:
327 ConvOpt sharedOpt_;
328 };
329
330 template<unsigned int DIM, unsigned int EV>
331 class HessianOfGaussianSelectedEigenvalueFunctor{
332 public:
333 typedef ConvolutionOptions<DIM> ConvOpt;
334 HessianOfGaussianSelectedEigenvalueFunctor(const ConvOpt & convOpt)
335 : sharedOpt_(convOpt){}
336 template<class S, class D>
337 void operator()(const S & s, D & d)const{
338 typedef typename vigra::NumericTraits<typename S::value_type>::RealPromote RealType;
339
340 // compute the hessian of gaussian and extract eigenvalue
341 vigra::MultiArray<DIM, TinyVector<RealType, int(DIM*(DIM+1)/2)> > hessianOfGaussianRes(s.shape());
342 vigra::hessianOfGaussianMultiArray(s, hessianOfGaussianRes, sharedOpt_);
343
344 vigra::MultiArray<DIM, TinyVector<RealType, DIM > > allEigenvalues(s.shape());
345 vigra::tensorEigenvaluesMultiArray(hessianOfGaussianRes, allEigenvalues);
346
347 d = allEigenvalues.bindElementChannel(EV);
348 }
349 template<class S, class D,class SHAPE>
350 void operator()(const S & s, D & d, const SHAPE & roiBegin, const SHAPE & roiEnd){
351
352 typedef typename vigra::NumericTraits<typename S::value_type>::RealPromote RealType;
353
354 // compute the hessian of gaussian and extract eigenvalue
355 vigra::MultiArray<DIM, TinyVector<RealType, int(DIM*(DIM+1)/2)> > hessianOfGaussianRes(roiEnd-roiBegin);
356 ConvOpt localOpt(sharedOpt_);
357 localOpt.subarray(roiBegin, roiEnd);
358 vigra::hessianOfGaussianMultiArray(s, hessianOfGaussianRes, localOpt);
359
360 vigra::MultiArray<DIM, TinyVector<RealType, DIM > > allEigenvalues(roiEnd-roiBegin);
361 vigra::tensorEigenvaluesMultiArray(hessianOfGaussianRes, allEigenvalues);
362
363 d = allEigenvalues.bindElementChannel(EV);
364 }
365 private:
366 ConvOpt sharedOpt_;
367 };
368
369
370 template<unsigned int DIM>
371 class HessianOfGaussianFirstEigenvalueFunctor
372 : public HessianOfGaussianSelectedEigenvalueFunctor<DIM, 0>{
373 public:
374 typedef ConvolutionOptions<DIM> ConvOpt;
375 HessianOfGaussianFirstEigenvalueFunctor(const ConvOpt & convOpt)
376 : HessianOfGaussianSelectedEigenvalueFunctor<DIM, 0>(convOpt){}
377 };
378
379 template<unsigned int DIM>
380 class HessianOfGaussianLastEigenvalueFunctor
381 : public HessianOfGaussianSelectedEigenvalueFunctor<DIM, DIM-1>{
382 public:
383 typedef ConvolutionOptions<DIM> ConvOpt;
384 HessianOfGaussianLastEigenvalueFunctor(const ConvOpt & convOpt)
385 : HessianOfGaussianSelectedEigenvalueFunctor<DIM, DIM-1>(convOpt){}
386 };
387
388
389
390
391
392
393 /// \warning this functions is deprecated
394 /// and should not be used from end users
395 template<unsigned int N>
396 vigra::TinyVector< vigra::MultiArrayIndex, N > getBorder(
397 const BlockwiseConvolutionOptions<N> & opt,
398 const size_t order,
399 const bool usesOuterScale = false
400 ){
401 vigra::TinyVector< vigra::MultiArrayIndex, N > res(vigra::SkipInitialization);
402
403 if(opt.getFilterWindowSize()<=0.00001){
404 for(size_t d=0; d<N; ++d){
405 double stdDev = opt.getStdDev()[d];
406 if(usesOuterScale)
407 stdDev += opt.getOuterScale()[d];
408 res[d] = static_cast<MultiArrayIndex>(3.0 * stdDev + 0.5*static_cast<double>(order)+0.5);
409 }
410 }
411 else{
412 throw std::runtime_error("blockwise filters do not allow a user defined FilterWindowSize");
413 }
414 return res;
415 }
416
417} // end namespace blockwise
418
419#define VIGRA_BLOCKWISE(FUNCTOR, FUNCTION, ORDER, USES_OUTER_SCALE) \
420 template <unsigned int N, class T1, class S1, class T2, class S2> \
421 void FUNCTION( \
422 MultiArrayView<N, T1, S1> const & source, \
423 MultiArrayView<N, T2, S2> dest, \
424 BlockwiseConvolutionOptions<N> const & options \
425) \
426{ \
427 typedef MultiBlocking<N, vigra::MultiArrayIndex> Blocking; \
428 typedef typename Blocking::Shape Shape; \
429 const Shape border = blockwise::getBorder(options, ORDER, USES_OUTER_SCALE); \
430 BlockwiseConvolutionOptions<N> subOptions(options); \
431 subOptions.subarray(Shape(0), Shape(0)); \
432 const Blocking blocking(source.shape(), options.template getBlockShapeN<N>()); \
433 blockwise::FUNCTOR<N> f(subOptions); \
434 blockwise::blockwiseCaller(source, dest, f, blocking, border, options); \
435}
436
437VIGRA_BLOCKWISE(GaussianSmoothFunctor, gaussianSmoothMultiArray, 0, false );
438VIGRA_BLOCKWISE(GaussianGradientFunctor, gaussianGradientMultiArray, 1, false );
439VIGRA_BLOCKWISE(SymmetricGradientFunctor, symmetricGradientMultiArray, 1, false );
440VIGRA_BLOCKWISE(GaussianDivergenceFunctor, gaussianDivergenceMultiArray, 1, false );
441VIGRA_BLOCKWISE(HessianOfGaussianFunctor, hessianOfGaussianMultiArray, 2, false );
442VIGRA_BLOCKWISE(HessianOfGaussianEigenvaluesFunctor, hessianOfGaussianEigenvaluesMultiArray, 2, false );
443VIGRA_BLOCKWISE(HessianOfGaussianFirstEigenvalueFunctor, hessianOfGaussianFirstEigenvalueMultiArray, 2, false );
444VIGRA_BLOCKWISE(HessianOfGaussianLastEigenvalueFunctor, hessianOfGaussianLastEigenvalueMultiArray, 2, false );
445VIGRA_BLOCKWISE(LaplacianOfGaussianFunctor, laplacianOfGaussianMultiArray, 2, false );
446VIGRA_BLOCKWISE(GaussianGradientMagnitudeFunctor, gaussianGradientMagnitudeMultiArray, 1, false );
447VIGRA_BLOCKWISE(StructureTensorFunctor, structureTensorMultiArray, 1, true );
448
449#undef VIGRA_BLOCKWISE
450
451 // alternative name for backward compatibility
452 template <unsigned int N, class T1, class S1, class T2, class S2>
453 inline void
454 gaussianGradientMagnitude(
455 MultiArrayView<N, T1, S1> const & source,
456 MultiArrayView<N, T2, S2> dest,
457 BlockwiseConvolutionOptions<N> const & options)
458{
459 gaussianGradientMagnitudeMultiArray(source, dest, options);
460}
461
462
463} // end namespace vigra
464
465#endif // VIGRA_MULTI_BLOCKWISE_HXX
vigra::ArrayVector
Definition array_vector.hxx:514
vigra::BlockwiseConvolutionOptions
Definition multi_blockwise.hxx:162
vigra::BlockwiseOptions
Definition multi_blockwise.hxx:56
vigra::BlockwiseOptions::blockShape
BlockwiseOptions & blockShape(MultiArrayIndex blockShape)
Definition multi_blockwise.hxx:134
vigra::BlockwiseOptions::blockShape
BlockwiseOptions & blockShape(const TinyVector< T, N > &blockShape)
Definition multi_blockwise.hxx:127
vigra::BlockwiseOptions::getBlockShape
Shape const & getBlockShape() const
Definition multi_blockwise.hxx:71
vigra::BlockwiseOptions::getBlockShapeN
TinyVector< MultiArrayIndex, N > getBlockShapeN() const
Definition multi_blockwise.hxx:89
vigra::BlockwiseOptions::blockShape
BlockwiseOptions & blockShape(const Shape &blockShape)
Definition multi_blockwise.hxx:114
vigra::MultiArrayView
Base class for, and view to, MultiArray.
Definition multi_array.hxx:705
vigra::MultiBlocking
Definition multi_blocking.hxx:157
vigra::MultiBlocking::numBlocks
size_t numBlocks() const
total number of blocks
Definition multi_blocking.hxx:225
vigra::ParallelOptions::numThreads
ParallelOptions & numThreads(const int n)
Set the number of threads or one of the constants Auto, Nice and NoThreads.
Definition threadpool.hxx:111
vigra::ParallelOptions::getNumThreads
int getNumThreads() const
Get desired number of threads.
Definition threadpool.hxx:85
vigra::TinyVector
Class for fixed size vectors.
Definition tinyvector.hxx:1008
vigra
vigra::gaussianGradientMagnitude
void gaussianGradientMagnitude(...)
Calculate the gradient magnitude by means of a 1st derivatives of Gaussian filter.
vigra::structureTensorMultiArray
void structureTensorMultiArray(...)
Calculate the structure tensor of a multi-dimensional arrays.
vigra::laplacianOfGaussianMultiArray
void laplacianOfGaussianMultiArray(...)
Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters.
vigra::gaussianDivergenceMultiArray
void gaussianDivergenceMultiArray(...)
Calculate the divergence of a vector field using Gaussian derivative filters.
vigra::gaussianSmoothMultiArray
void gaussianSmoothMultiArray(...)
Isotropic Gaussian smoothing of a multi-dimensional arrays.
vigra::tensorEigenvaluesMultiArray
void tensorEigenvaluesMultiArray(...)
Calculate the tensor eigenvalues for every element of a N-D tensor array.
vigra::MultiArrayIndex
std::ptrdiff_t MultiArrayIndex
Definition multi_fwd.hxx:60
vigra::hessianOfGaussianMultiArray
void hessianOfGaussianMultiArray(...)
Calculate Hessian matrix of a N-dimensional arrays using Gaussian derivative filters.
vigra::gaussianGradientMultiArray
void gaussianGradientMultiArray(...)
Calculate Gaussian gradient of a multi-dimensional arrays.
vigra::symmetricGradientMultiArray
void symmetricGradientMultiArray(...)
Calculate gradient of a multi-dimensional arrays using symmetric difference filters.
vigra::parallel_foreach
void parallel_foreach(...)
Apply a functor to all items in a range in parallel.

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.12.4