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

multi_blockwise.hxx VIGRA

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  */
55 : public ParallelOptions
56 {
57 public:
59 
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>
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  */
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>
128  Shape(blockShape.begin(), blockShape.end()).swap(blockShape_);
129  return *this;
130  }
131 
132  /** Specify square block shape by its side length.
133  */
135  Shape(1, blockShape).swap(blockShape_);
136  return *this;
137  }
138 
139  BlockwiseOptions & numThreads(const int n)
140  {
142  return *this;
143  }
144 
145  void setNumThreads(const int n)
146  {
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>
161 : public BlockwiseOptions
162 , public ConvolutionOptions<N>{
163 public:
165  : BlockwiseOptions(),
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(
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 
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(
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  : convOpt_(convOpt){} \
278  template<class S, class D> \
279  void operator()(const S & s, D & d)const{ \
280  FUNCTION_NAME(s, d, convOpt_); \
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 convOpt(convOpt_); \
285  convOpt.subarray(roiBegin, roiEnd); \
286  FUNCTION_NAME(s, d, convOpt); \
287  } \
288  private: \
289  ConvOpt convOpt_; \
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  : convOpt_(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, convOpt_);
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_.subarray(roiBegin, roiEnd);
322  vigra::hessianOfGaussianMultiArray(s, hessianOfGaussianRes, convOpt_);
323  vigra::tensorEigenvaluesMultiArray(hessianOfGaussianRes, d);
324  }
325  private:
326  ConvOpt convOpt_;
327  };
328 
329  template<unsigned int DIM, unsigned int EV>
330  class HessianOfGaussianSelectedEigenvalueFunctor{
331  public:
332  typedef ConvolutionOptions<DIM> ConvOpt;
333  HessianOfGaussianSelectedEigenvalueFunctor(const ConvOpt & convOpt)
334  : convOpt_(convOpt){}
335  template<class S, class D>
336  void operator()(const S & s, D & d)const{
337  typedef typename vigra::NumericTraits<typename S::value_type>::RealPromote RealType;
338 
339  // compute the hessian of gaussian and extract eigenvalue
340  vigra::MultiArray<DIM, TinyVector<RealType, int(DIM*(DIM+1)/2)> > hessianOfGaussianRes(s.shape());
341  vigra::hessianOfGaussianMultiArray(s, hessianOfGaussianRes, convOpt_);
342 
343  vigra::MultiArray<DIM, TinyVector<RealType, DIM > > allEigenvalues(s.shape());
344  vigra::tensorEigenvaluesMultiArray(hessianOfGaussianRes, allEigenvalues);
345 
346  d = allEigenvalues.bindElementChannel(EV);
347  }
348  template<class S, class D,class SHAPE>
349  void operator()(const S & s, D & d, const SHAPE & roiBegin, const SHAPE & roiEnd){
350 
351  typedef typename vigra::NumericTraits<typename S::value_type>::RealPromote RealType;
352 
353  // compute the hessian of gaussian and extract eigenvalue
354  vigra::MultiArray<DIM, TinyVector<RealType, int(DIM*(DIM+1)/2)> > hessianOfGaussianRes(roiEnd-roiBegin);
355  convOpt_.subarray(roiBegin, roiEnd);
356  vigra::hessianOfGaussianMultiArray(s, hessianOfGaussianRes, convOpt_);
357 
358  vigra::MultiArray<DIM, TinyVector<RealType, DIM > > allEigenvalues(roiEnd-roiBegin);
359  vigra::tensorEigenvaluesMultiArray(hessianOfGaussianRes, allEigenvalues);
360 
361  d = allEigenvalues.bindElementChannel(EV);
362  }
363  private:
364  ConvOpt convOpt_;
365  };
366 
367 
368  template<unsigned int DIM>
369  class HessianOfGaussianFirstEigenvalueFunctor
370  : public HessianOfGaussianSelectedEigenvalueFunctor<DIM, 0>{
371  public:
372  typedef ConvolutionOptions<DIM> ConvOpt;
373  HessianOfGaussianFirstEigenvalueFunctor(const ConvOpt & convOpt)
374  : HessianOfGaussianSelectedEigenvalueFunctor<DIM, 0>(convOpt){}
375  };
376 
377  template<unsigned int DIM>
378  class HessianOfGaussianLastEigenvalueFunctor
379  : public HessianOfGaussianSelectedEigenvalueFunctor<DIM, DIM-1>{
380  public:
381  typedef ConvolutionOptions<DIM> ConvOpt;
382  HessianOfGaussianLastEigenvalueFunctor(const ConvOpt & convOpt)
383  : HessianOfGaussianSelectedEigenvalueFunctor<DIM, DIM-1>(convOpt){}
384  };
385 
386 
387 
388 
389 
390 
391  /// \warning this functions is deprecated
392  /// and should not be used from end users
393  template<unsigned int N>
395  const BlockwiseConvolutionOptions<N> & opt,
396  const size_t order,
397  const bool usesOuterScale = false
398  ){
399  vigra::TinyVector< vigra::MultiArrayIndex, N > res(vigra::SkipInitialization);
400 
401  if(opt.getFilterWindowSize()<=0.00001){
402  for(size_t d=0; d<N; ++d){
403  double stdDev = opt.getStdDev()[d];
404  if(usesOuterScale)
405  stdDev += opt.getOuterScale()[d];
406  res[d] = static_cast<MultiArrayIndex>(3.0 * stdDev + 0.5*static_cast<double>(order)+0.5);
407  }
408  }
409  else{
410  throw std::runtime_error("blockwise filters do not allow a user defined FilterWindowSize");
411  }
412  return res;
413  }
414 
415 } // end namespace blockwise
416 
417 #define VIGRA_BLOCKWISE(FUNCTOR, FUNCTION, ORDER, USES_OUTER_SCALE) \
418 template <unsigned int N, class T1, class S1, class T2, class S2> \
419 void FUNCTION( \
420  MultiArrayView<N, T1, S1> const & source, \
421  MultiArrayView<N, T2, S2> dest, \
422  BlockwiseConvolutionOptions<N> const & options \
423 ) \
424 { \
425  typedef MultiBlocking<N, vigra::MultiArrayIndex> Blocking; \
426  typedef typename Blocking::Shape Shape; \
427  const Shape border = blockwise::getBorder(options, ORDER, USES_OUTER_SCALE); \
428  BlockwiseConvolutionOptions<N> subOptions(options); \
429  subOptions.subarray(Shape(0), Shape(0)); \
430  const Blocking blocking(source.shape(), options.template getBlockShapeN<N>()); \
431  blockwise::FUNCTOR<N> f(subOptions); \
432  blockwise::blockwiseCaller(source, dest, f, blocking, border, options); \
433 }
434 
435 VIGRA_BLOCKWISE(GaussianSmoothFunctor, gaussianSmoothMultiArray, 0, false );
436 VIGRA_BLOCKWISE(GaussianGradientFunctor, gaussianGradientMultiArray, 1, false );
437 VIGRA_BLOCKWISE(SymmetricGradientFunctor, symmetricGradientMultiArray, 1, false );
438 VIGRA_BLOCKWISE(GaussianDivergenceFunctor, gaussianDivergenceMultiArray, 1, false );
439 VIGRA_BLOCKWISE(HessianOfGaussianFunctor, hessianOfGaussianMultiArray, 2, false );
440 VIGRA_BLOCKWISE(HessianOfGaussianEigenvaluesFunctor, hessianOfGaussianEigenvaluesMultiArray, 2, false );
441 VIGRA_BLOCKWISE(HessianOfGaussianFirstEigenvalueFunctor, hessianOfGaussianFirstEigenvalueMultiArray, 2, false );
442 VIGRA_BLOCKWISE(HessianOfGaussianLastEigenvalueFunctor, hessianOfGaussianLastEigenvalueMultiArray, 2, false );
443 VIGRA_BLOCKWISE(LaplacianOfGaussianFunctor, laplacianOfGaussianMultiArray, 2, false );
444 VIGRA_BLOCKWISE(GaussianGradientMagnitudeFunctor, gaussianGradientMagnitudeMultiArray, 1, false );
445 VIGRA_BLOCKWISE(StructureTensorFunctor, structureTensorMultiArray, 1, true );
446 
447 #undef VIGRA_BLOCKWISE
448 
449  // alternative name for backward compatibility
450 template <unsigned int N, class T1, class S1, class T2, class S2>
451 inline void
453  MultiArrayView<N, T1, S1> const & source,
454  MultiArrayView<N, T2, S2> dest,
455  BlockwiseConvolutionOptions<N> const & options)
456 {
457  gaussianGradientMagnitudeMultiArray(source, dest, options);
458 }
459 
460 
461 } // end namespace vigra
462 
463 #endif // VIGRA_MULTI_BLOCKWISE_HXX
void symmetricGradientMultiArray(...)
Calculate gradient of a multi-dimensional arrays using symmetric difference filters.
size_t numBlocks() const
total number of blocks
Definition: multi_blocking.hxx:225
const difference_type & shape() const
Definition: multi_array.hxx:1594
void gaussianDivergenceMultiArray(...)
Calculate the divergence of a vector field using Gaussian derivative filters.
int getNumThreads() const
Get desired number of threads.
Definition: threadpool.hxx:85
iterator begin()
Definition: multi_array.hxx:1867
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2420
iterator end()
Definition: tinyvector.hxx:864
ParallelOptions & numThreads(const int n)
Set the number of threads or one of the constants Auto, Nice and NoThreads.
Definition: threadpool.hxx:111
std::ptrdiff_t MultiArrayIndex
Definition: multi_fwd.hxx:60
BlockwiseOptions & blockShape(MultiArrayIndex blockShape)
Definition: multi_blockwise.hxx:134
Definition: multi_blocking.hxx:49
void gaussianGradientMultiArray(...)
Calculate Gaussian gradient of a multi-dimensional arrays.
void laplacianOfGaussianMultiArray(...)
Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters.
TinyVector< MultiArrayIndex, N > getBlockShapeN() const
Definition: multi_blockwise.hxx:89
iterator begin()
Definition: tinyvector.hxx:861
vigra::GridGraph< N, DirectedTag >::vertex_descriptor source(typename vigra::GridGraph< N, DirectedTag >::edge_descriptor const &e, vigra::GridGraph< N, DirectedTag > const &g)
Get a vertex descriptor for the start vertex of edge e in graph g (API: boost).
Definition: multi_gridgraph.hxx:2943
void gaussianGradientMagnitude(...)
Calculate the gradient magnitude by means of a 1st derivatives of Gaussian filter.
Definition: multi_blockwise.hxx:54
Options class template for convolutions.
Definition: multi_convolution.hxx:335
void parallel_foreach(...)
Apply a functor to all items in a range in parallel.
Class for fixed size vectors.This class contains an array of size SIZE of the specified VALUETYPE...
Definition: accessor.hxx:940
Option base class for parallel algorithms.
Definition: threadpool.hxx:61
Shape const & getBlockShape() const
Definition: multi_blockwise.hxx:71
BlockwiseOptions & blockShape(const TinyVector< T, N > &blockShape)
Definition: multi_blockwise.hxx:127
Definition: multi_blockwise.hxx:160
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:650
const_pointer data() const
Definition: array_vector.hxx:209
size_type size() const
Definition: array_vector.hxx:358
void gaussianSmoothMultiArray(...)
Isotropic Gaussian smoothing of a multi-dimensional arrays.
MultiArrayView subarray(difference_type p, difference_type q) const
Definition: multi_array.hxx:1474
void hessianOfGaussianMultiArray(...)
Calculate Hessian matrix of a N-dimensional arrays using Gaussian derivative filters.
void tensorEigenvaluesMultiArray(...)
Calculate the tensor eigenvalues for every element of a N-D tensor array.
BlockwiseOptions & blockShape(const Shape &blockShape)
Definition: multi_blockwise.hxx:114
void structureTensorMultiArray(...)
Calculate the structure tensor of a multi-dimensional arrays.

© 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.11.0 (Thu Mar 17 2016)