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

multi_convolution.hxx VIGRA

1 //-- -*- c++ -*-
2 /************************************************************************/
3 /* */
4 /* Copyright 2003 by Christian-Dennis Rahn */
5 /* and Ullrich Koethe */
6 /* */
7 /* This file is part of the VIGRA computer vision library. */
8 /* The VIGRA Website is */
9 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
10 /* Please direct questions, bug reports, and contributions to */
11 /* ullrich.koethe@iwr.uni-heidelberg.de or */
12 /* vigra@informatik.uni-hamburg.de */
13 /* */
14 /* Permission is hereby granted, free of charge, to any person */
15 /* obtaining a copy of this software and associated documentation */
16 /* files (the "Software"), to deal in the Software without */
17 /* restriction, including without limitation the rights to use, */
18 /* copy, modify, merge, publish, distribute, sublicense, and/or */
19 /* sell copies of the Software, and to permit persons to whom the */
20 /* Software is furnished to do so, subject to the following */
21 /* conditions: */
22 /* */
23 /* The above copyright notice and this permission notice shall be */
24 /* included in all copies or substantial portions of the */
25 /* Software. */
26 /* */
27 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
28 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
29 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
30 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
31 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
32 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
33 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
34 /* OTHER DEALINGS IN THE SOFTWARE. */
35 /* */
36 /************************************************************************/
37 
38 #ifndef VIGRA_MULTI_CONVOLUTION_H
39 #define VIGRA_MULTI_CONVOLUTION_H
40 
41 #include "separableconvolution.hxx"
42 #include "array_vector.hxx"
43 #include "multi_array.hxx"
44 #include "accessor.hxx"
45 #include "numerictraits.hxx"
46 #include "navigator.hxx"
47 #include "metaprogramming.hxx"
48 #include "multi_pointoperators.hxx"
49 #include "multi_math.hxx"
50 #include "functorexpression.hxx"
51 #include "tinyvector.hxx"
52 #include "algorithm.hxx"
53 
54 
55 #include <iostream>
56 
57 namespace vigra
58 {
59 
60 namespace detail
61 {
62 
63 struct DoubleYielder
64 {
65  const double value;
66  DoubleYielder(double v, unsigned, const char *const) : value(v) {}
67  DoubleYielder(double v) : value(v) {}
68  void operator++() {}
69  double operator*() const { return value; }
70 };
71 
72 template <typename X>
73 struct IteratorDoubleYielder
74 {
75  X it;
76  IteratorDoubleYielder(X i, unsigned, const char *const) : it(i) {}
77  IteratorDoubleYielder(X i) : it(i) {}
78  void operator++() { ++it; }
79  double operator*() const { return *it; }
80 };
81 
82 template <typename X>
83 struct SequenceDoubleYielder
84 {
85  typename X::const_iterator it;
86  SequenceDoubleYielder(const X & seq, unsigned dim,
87  const char *const function_name = "SequenceDoubleYielder")
88  : it(seq.begin())
89  {
90  if (seq.size() == dim)
91  return;
92  std::string msg = "(): Parameter number be equal to the number of spatial dimensions.";
93  vigra_precondition(false, function_name + msg);
94  }
95  void operator++() { ++it; }
96  double operator*() const { return *it; }
97 };
98 
99 template <typename X>
100 struct WrapDoubleIterator
101 {
102  typedef
103  typename IfBool< IsConvertibleTo<X, double>::value,
104  DoubleYielder,
105  typename IfBool< IsIterator<X>::value || IsArray<X>::value,
106  IteratorDoubleYielder<X>,
107  SequenceDoubleYielder<X>
108  >::type
109  >::type type;
110 };
111 
112 template <class Param1, class Param2, class Param3>
113 struct WrapDoubleIteratorTriple
114 {
115  typename WrapDoubleIterator<Param1>::type sigma_eff_it;
116  typename WrapDoubleIterator<Param2>::type sigma_d_it;
117  typename WrapDoubleIterator<Param3>::type step_size_it;
118  WrapDoubleIteratorTriple(Param1 sigma_eff, Param2 sigma_d, Param3 step_size)
119  : sigma_eff_it(sigma_eff), sigma_d_it(sigma_d), step_size_it(step_size) {}
120  void operator++()
121  {
122  ++sigma_eff_it;
123  ++sigma_d_it;
124  ++step_size_it;
125  }
126  double sigma_eff() const { return *sigma_eff_it; }
127  double sigma_d() const { return *sigma_d_it; }
128  double step_size() const { return *step_size_it; }
129  static void sigma_precondition(double sigma, const char *const function_name)
130  {
131  if (sigma < 0.0)
132  {
133  std::string msg = "(): Scale must be positive.";
134  vigra_precondition(false, function_name + msg);
135  }
136  }
137  double sigma_scaled(const char *const function_name = "unknown function ",
138  bool allow_zero = false) const
139  {
140  sigma_precondition(sigma_eff(), function_name);
141  sigma_precondition(sigma_d(), function_name);
142  double sigma_squared = sq(sigma_eff()) - sq(sigma_d());
143  if (sigma_squared > 0.0 || (allow_zero && sigma_squared == 0.0))
144  {
145  return std::sqrt(sigma_squared) / step_size();
146  }
147  else
148  {
149  std::string msg = "(): Scale would be imaginary";
150  if(!allow_zero)
151  msg += " or zero";
152  vigra_precondition(false, function_name + msg + ".");
153  return 0;
154  }
155  }
156 };
157 
158 template <unsigned dim>
159 struct multiArrayScaleParam
160 {
161  typedef TinyVector<double, dim> p_vector;
162  typedef typename p_vector::const_iterator return_type;
163  p_vector vec;
164 
165  template <class Param>
166  multiArrayScaleParam(Param val, const char *const function_name = "multiArrayScaleParam")
167  {
168  typename WrapDoubleIterator<Param>::type in(val, dim, function_name);
169  for (unsigned i = 0; i != dim; ++i, ++in)
170  vec[i] = *in;
171  }
172  return_type operator()() const
173  {
174  return vec.begin();
175  }
176  static void precondition(unsigned n_par, const char *const function_name = "multiArrayScaleParam")
177  {
178  char n[3] = "0.";
179  n[0] += dim;
180  std::string msg = "(): dimension parameter must be ";
181  vigra_precondition(dim == n_par, function_name + msg + n);
182  }
183  multiArrayScaleParam(double v0, double v1, const char *const function_name = "multiArrayScaleParam")
184  {
185  precondition(2, function_name);
186  vec = p_vector(v0, v1);
187  }
188  multiArrayScaleParam(double v0, double v1, double v2, const char *const function_name = "multiArrayScaleParam")
189  {
190  precondition(3, function_name);
191  vec = p_vector(v0, v1, v2);
192  }
193  multiArrayScaleParam(double v0, double v1, double v2, double v3, const char *const function_name = "multiArrayScaleParam")
194  {
195  precondition(4, function_name);
196  vec = p_vector(v0, v1, v2, v3);
197  }
198  multiArrayScaleParam(double v0, double v1, double v2, double v3, double v4, const char *const function_name = "multiArrayScaleParam")
199  {
200  precondition(5, function_name);
201  vec = p_vector(v0, v1, v2, v3, v4);
202  }
203 };
204 
205 } // namespace detail
206 
207 #define VIGRA_CONVOLUTION_OPTIONS(function_name, default_value, member_name, getter_setter_name) \
208  template <class Param> \
209  ConvolutionOptions & function_name(const Param & val) \
210  { \
211  member_name = ParamVec(val, "ConvolutionOptions::" #function_name); \
212  return *this; \
213  } \
214  ConvolutionOptions & function_name() \
215  { \
216  member_name = ParamVec(default_value, "ConvolutionOptions::" #function_name); \
217  return *this; \
218  } \
219  ConvolutionOptions & function_name(double v0, double v1) \
220  { \
221  member_name = ParamVec(v0, v1, "ConvolutionOptions::" #function_name); \
222  return *this; \
223  } \
224  ConvolutionOptions & function_name(double v0, double v1, double v2) \
225  { \
226  member_name = ParamVec(v0, v1, v2, "ConvolutionOptions::" #function_name); \
227  return *this; \
228  } \
229  ConvolutionOptions & function_name(double v0, double v1, double v2, double v3) \
230  { \
231  member_name = ParamVec(v0, v1, v2, v3, "ConvolutionOptions::" #function_name); \
232  return *this; \
233  } \
234  ConvolutionOptions & function_name(double v0, double v1, double v2, double v3, double v4) \
235  { \
236  member_name = ParamVec(v0, v1, v2, v3, v4, "ConvolutionOptions::" #function_name); \
237  return *this; \
238  } \
239  typename ParamVec::p_vector get##getter_setter_name()const{ \
240  return member_name.vec; \
241  } \
242  void set##getter_setter_name(typename ParamVec::p_vector vec){ \
243  member_name.vec = vec; \
244  }
245 
246 /** \brief Options class template for convolutions.
247 
248  <b>\#include</b> <vigra/multi_convolution.hxx><br/>
249  Namespace: vigra
250 
251  This class enables the calculation of scale space convolutions
252  such as \ref gaussianGradientMultiArray() on data with anisotropic
253  discretization. For these, the result of the ordinary calculation
254  has to be multiplied by factors of \f$1/w^{n}\f$ for each dimension,
255  where \f$w\f$ is the step size of the grid in said dimension and
256  \f$n\f$ is the differential order of the convolution, e.g., 1 for
257  gaussianGradientMultiArray(), and 0 for gaussianSmoothMultiArray(),
258  respectively. Also for each dimension in turn, the convolution's scale
259  parameter \f$\sigma\f$ has to be replaced by
260  \f$\sqrt{\sigma_\mathrm{eff}^2 - \sigma_\mathrm{D}^2}\Big/w\f$,
261  where \f$\sigma_\mathrm{eff}\f$ is the resulting effective filtering
262  scale. The data is assumed to be already filtered by a
263  gaussian smoothing with the scale parameter \f$\sigma_\mathrm{D}\f$
264  (such as by measuring equipment). All of the above changes are
265  automatically employed by the convolution functions for <tt>MultiArray</tt>s
266  if a corresponding options object is provided.
267 
268  The <tt>ConvolutionOptions</tt> class must be parameterized by the dimension
269  <tt>dim</tt>
270  of the <tt>MultiArray</tt>s on which it is used. The actual per-axis
271  options are set by (overloaded) member functions explained below,
272  or else default to neutral values corresponding to the absence of the
273  particular option.
274 
275  All member functions set <tt>dim</tt> values of the respective convolution
276  option, one for each dimension. They may be set explicitly by multiple
277  arguments for up to five dimensions, or by a single argument to the same
278  value for all dimensions. For the general case, a single argument that is
279  either a C-syle array, an iterator, or a C++ standard library style
280  sequence (such as <tt>std::vector</tt>, with member functions <tt>begin()</tt>
281  and <tt>size()</tt>) supplies the option values for any number of dimensions.
282 
283  Note that the return value of all member functions is <tt>*this</tt>, which
284  provides the mechanism for concatenating member function calls as shown below.
285 
286  <b>usage with explicit parameters:</b>
287 
288  \code
289  ConvolutionOptions<2> opt = ConvolutionOptions<2>().stepSize(1, 2.3);
290  \endcode
291 
292  <b>usage with arrays:</b>
293 
294  \code
295  const double step_size[3] = { x_scale, y_scale, z_scale };
296  ConvolutionOptions<3> opt = ConvolutionOptions<3>().stepSize(step_size);
297  \endcode
298 
299  <b>usage with C++ standard library style sequences:</b>
300 
301  \code
302  TinyVector<double, 4> step_size(1, 1, 2.0, 1.5);
303  TinyVector<double, 4> r_sigmas(1, 1, 2.3, 3.2);
304  ConvolutionOptions<4> opt = ConvolutionOptions<4>().stepSize(step_size).resolutionStdDev(r_sigmas);
305  \endcode
306 
307  <b>usage with iterators:</b>
308 
309  \code
310  ArrayVector<double> step_size;
311  step_size.push_back(0);
312  step_size.push_back(3);
313  step_size.push_back(4);
314  ArrayVector<double>::iterator i = step_size.begin();
315  ++i;
316  ConvolutionOptions<2> opt = ConvolutionOptions<2>().stepSize(i);
317  \endcode
318 
319  <b>general usage in a convolution function call:</b>
320 
321  \code
322  MultiArray<3, double> test_image;
323  MultiArray<3, double> out_image;
324 
325  double scale = 5.0;
326  gaussianSmoothMultiArray(test_image, out_image, scale,
327  ConvolutionOptions<3>()
328  .stepSize (1, 1, 3.2)
329  .resolutionStdDev(1, 1, 4)
330  );
331  \endcode
332 
333 */
334 template <unsigned dim>
336 {
337  public:
338  typedef typename MultiArrayShape<dim>::type Shape;
339  typedef detail::multiArrayScaleParam<dim> ParamVec;
340  typedef typename ParamVec::return_type ParamIt;
341 
342  ParamVec sigma_eff;
343  ParamVec sigma_d;
344  ParamVec step_size;
345  ParamVec outer_scale;
346  double window_ratio;
347  Shape from_point, to_point;
348 
350  : sigma_eff(0.0),
351  sigma_d(0.0),
352  step_size(1.0),
353  outer_scale(0.0),
354  window_ratio(0.0)
355  {}
356 
357  typedef typename detail::WrapDoubleIteratorTriple<ParamIt, ParamIt, ParamIt>
358  ScaleIterator;
359  typedef typename detail::WrapDoubleIterator<ParamIt>::type
360  StepIterator;
361 
362  ScaleIterator scaleParams() const
363  {
364  return ScaleIterator(sigma_eff(), sigma_d(), step_size());
365  }
366  StepIterator stepParams() const
367  {
368  return StepIterator(step_size());
369  }
370 
371  ConvolutionOptions outerOptions() const
372  {
373  ConvolutionOptions outer = *this;
374  // backward-compatible values:
375  return outer.stdDev(outer_scale()).resolutionStdDev(0.0);
376  }
377 
378  // Step size per axis.
379  // Default: dim values of 1.0
380  VIGRA_CONVOLUTION_OPTIONS(stepSize, 1.0, step_size, StepSize)
381 #ifdef DOXYGEN
382  /** Step size(s) per axis, i.e., the distance between two
383  adjacent pixels. Required for <tt>MultiArray</tt>
384  containing anisotropic data.
385 
386  Note that a convolution containing a derivative operator
387  of order <tt>n</tt> results in a multiplication by
388  \f${\rm stepSize}^{-n}\f$ for each axis.
389  Also, the above standard deviations
390  are scaled according to the step size of each axis.
391  Default value for the options object if this member function is not
392  used: A value of 1.0 for each dimension.
393  */
394  ConvolutionOptions<dim> & stepSize(...);
395 #endif
396 
397  // Resolution standard deviation per axis.
398  // Default: dim values of 0.0
399  VIGRA_CONVOLUTION_OPTIONS(resolutionStdDev, 0.0, sigma_d, ResolutionStdDev)
400 #ifdef DOXYGEN
401  /** Resolution standard deviation(s) per axis, i.e., a supposed
402  pre-existing gaussian filtering by this value.
403 
404  The standard deviation actually used by the convolution operators
405  is \f$\sqrt{{\rm sigma}^{2} - {\rm resolutionStdDev}^{2}}\f$ for each
406  axis.
407  Default value for the options object if this member function is not
408  used: A value of 0.0 for each dimension.
409  */
410  ConvolutionOptions<dim> & resolutionStdDev(...);
411 #endif
412 
413  // Standard deviation of scale space operators.
414  // Default: dim values of 0.0
415  VIGRA_CONVOLUTION_OPTIONS(stdDev, 0.0, sigma_eff, StdDev)
416  VIGRA_CONVOLUTION_OPTIONS(innerScale, 0.0, sigma_eff, InnerScale)
417 
418 #ifdef DOXYGEN
419  /** Standard deviation(s) of scale space operators, or inner scale(s) for \ref structureTensorMultiArray().
420 
421  Usually not
422  needed, since a single value for all axes may be specified as a parameter
423  <tt>sigma</tt> to the call of
424  an convolution operator such as \ref gaussianGradientMultiArray(), and
425  anisotropic data requiring the use of the stepSize() member function.
426  Default value for the options object if this member function is not
427  used: A value of 0.0 for each dimension.
428  */
429  ConvolutionOptions<dim> & stdDev(...);
430 
431  /** Standard deviation(s) of scale space operators, or inner scale(s) for \ref structureTensorMultiArray().
432 
433  Usually not
434  needed, since a single value for all axes may be specified as a parameter
435  <tt>sigma</tt> to the call of
436  an convolution operator such as \ref gaussianGradientMultiArray(), and
437  anisotropic data requiring the use of the stepSize() member function.
438  Default value for the options object if this member function is not
439  used: A value of 0.0 for each dimension.
440  */
441  ConvolutionOptions<dim> & innerScale(...);
442 #endif
443 
444  // Outer scale, for structure tensor.
445  // Default: dim values of 0.0
446  VIGRA_CONVOLUTION_OPTIONS(outerScale, 0.0, outer_scale, OuterScale)
447 #ifdef DOXYGEN
448  /** Standard deviation(s) of the second convolution of the
449  structure tensor.
450 
451  Usually not needed, since a single value for
452  all axes may be specified as a parameter <tt>outerScale</tt> to
453  the call of \ref structureTensorMultiArray(), and
454  anisotropic data requiring the use of the stepSize() member
455  function.
456  Default value for the options object if this member function is not
457  used: A value of 0.0 for each dimension.
458  */
459  ConvolutionOptions<dim> & outerScale(...);
460 #endif
461 
462  /** Size of the filter window as a multiple of the scale parameter.
463 
464  This option is only used for Gaussian filters and their derivatives.
465  By default, the window size of a Gaussian filter is automatically
466  determined such that the error resulting from restricting the
467  infinitely large Gaussian function to a finite size is minimized.
468  In particular, the window radius is determined as
469  <tt>radius = round(3.0 * sigma + 0.5 * order)</tt>, where 'order' is the
470  desired derivative order. In some cases, it is desirable to trade off
471  accuracy for speed, and this function can be used to request a smaller
472  window radius.
473 
474  Default: <tt>0.0</tt> (i.e. determine the window size automatically)
475  */
477  {
478  vigra_precondition(ratio >= 0.0,
479  "ConvolutionOptions::filterWindowSize(): ratio must not be negative.");
480  window_ratio = ratio;
481  return *this;
482  }
483 
484  double getFilterWindowSize() const {
485  return window_ratio;
486  }
487 
488  /** Restrict the filter to a subregion of the input array.
489 
490  This is useful for speeding up computations by ignoring irrelevant
491  areas in the array. <b>Note:</b> It is assumed that the output array
492  of the convolution has the size given in this function. Negative ROI
493  boundaries are interpreted relative to the end of the respective dimension
494  (i.e. <tt>if(to[k] < 0) to[k] += source.shape(k);</tt>).
495 
496  Default: <tt>from = Shape(), to = Shape()</tt> (i.e. use entire array)
497  */
498  ConvolutionOptions<dim> & subarray(Shape const & from, Shape const & to)
499  {
500  from_point = from;
501  to_point = to;
502  return *this;
503  }
504 
505  std::pair<Shape, Shape> getSubarray()const{
506  std::pair<Shape, Shape> res;
507  res.first = from_point;
508  res.second = to_point;
509  return res;
510  }
511 };
512 
513 namespace detail
514 {
515 
516 /********************************************************/
517 /* */
518 /* internalSeparableConvolveMultiArray */
519 /* */
520 /********************************************************/
521 
522 template <class SrcIterator, class SrcShape, class SrcAccessor,
523  class DestIterator, class DestAccessor, class KernelIterator>
524 void
525 internalSeparableConvolveMultiArrayTmp(
526  SrcIterator si, SrcShape const & shape, SrcAccessor src,
527  DestIterator di, DestAccessor dest, KernelIterator kit)
528 {
529  enum { N = 1 + SrcIterator::level };
530 
531  typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
532  typedef typename AccessorTraits<TmpType>::default_accessor TmpAcessor;
533 
534  // temporary array to hold the current line to enable in-place operation
535  ArrayVector<TmpType> tmp( shape[0] );
536 
537  typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
538  typedef MultiArrayNavigator<DestIterator, N> DNavigator;
539 
540  TmpAcessor acc;
541 
542  {
543  // only operate on first dimension here
544  SNavigator snav( si, shape, 0 );
545  DNavigator dnav( di, shape, 0 );
546 
547  for( ; snav.hasMore(); snav++, dnav++ )
548  {
549  // first copy source to tmp for maximum cache efficiency
550  copyLine(snav.begin(), snav.end(), src, tmp.begin(), acc);
551 
552  convolveLine(srcIterRange(tmp.begin(), tmp.end(), acc),
553  destIter( dnav.begin(), dest ),
554  kernel1d( *kit ) );
555  }
556  ++kit;
557  }
558 
559  // operate on further dimensions
560  for( int d = 1; d < N; ++d, ++kit )
561  {
562  DNavigator dnav( di, shape, d );
563 
564  tmp.resize( shape[d] );
565 
566  for( ; dnav.hasMore(); dnav++ )
567  {
568  // first copy source to tmp since convolveLine() cannot work in-place
569  copyLine(dnav.begin(), dnav.end(), dest, tmp.begin(), acc);
570 
571  convolveLine(srcIterRange(tmp.begin(), tmp.end(), acc),
572  destIter( dnav.begin(), dest ),
573  kernel1d( *kit ) );
574  }
575  }
576 }
577 
578 /********************************************************/
579 /* */
580 /* internalSeparableConvolveSubarray */
581 /* */
582 /********************************************************/
583 
584 template <class SrcIterator, class SrcShape, class SrcAccessor,
585  class DestIterator, class DestAccessor, class KernelIterator>
586 void
587 internalSeparableConvolveSubarray(
588  SrcIterator si, SrcShape const & shape, SrcAccessor src,
589  DestIterator di, DestAccessor dest, KernelIterator kit,
590  SrcShape const & start, SrcShape const & stop)
591 {
592  enum { N = 1 + SrcIterator::level };
593 
594  typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
595  typedef MultiArray<N, TmpType> TmpArray;
596  typedef typename TmpArray::traverser TmpIterator;
597  typedef typename AccessorTraits<TmpType>::default_accessor TmpAcessor;
598 
599  SrcShape sstart, sstop, axisorder, tmpshape;
600  TinyVector<double, N> overhead;
601  for(int k=0; k<N; ++k)
602  {
603  axisorder[k] = k;
604  sstart[k] = start[k] - kit[k].right();
605  if(sstart[k] < 0)
606  sstart[k] = 0;
607  sstop[k] = stop[k] - kit[k].left();
608  if(sstop[k] > shape[k])
609  sstop[k] = shape[k];
610  overhead[k] = double(sstop[k] - sstart[k]) / (stop[k] - start[k]);
611  }
612 
613  indexSort(overhead.begin(), overhead.end(), axisorder.begin(), std::greater<double>());
614  SrcShape dstart, dstop(sstop - sstart);
615  dstop[axisorder[0]] = stop[axisorder[0]] - start[axisorder[0]];
616 
617  // temporary array to hold the current line to enable in-place operation
618  MultiArray<N, TmpType> tmp(dstop);
619 
620  typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
621  typedef MultiArrayNavigator<TmpIterator, N> TNavigator;
622 
623  TmpAcessor acc;
624 
625  {
626  // only operate on first dimension here
627  SNavigator snav( si, sstart, sstop, axisorder[0]);
628  TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[0]);
629 
630  ArrayVector<TmpType> tmpline(sstop[axisorder[0]] - sstart[axisorder[0]]);
631 
632  int lstart = start[axisorder[0]] - sstart[axisorder[0]];
633  int lstop = lstart + (stop[axisorder[0]] - start[axisorder[0]]);
634 
635  for( ; snav.hasMore(); snav++, tnav++ )
636  {
637  // first copy source to tmp for maximum cache efficiency
638  copyLine(snav.begin(), snav.end(), src, tmpline.begin(), acc);
639 
640  convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
641  destIter(tnav.begin(), acc),
642  kernel1d( kit[axisorder[0]] ), lstart, lstop);
643  }
644  }
645 
646  // operate on further dimensions
647  for( int d = 1; d < N; ++d)
648  {
649  TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[d]);
650 
651  ArrayVector<TmpType> tmpline(dstop[axisorder[d]] - dstart[axisorder[d]]);
652 
653  int lstart = start[axisorder[d]] - sstart[axisorder[d]];
654  int lstop = lstart + (stop[axisorder[d]] - start[axisorder[d]]);
655 
656  for( ; tnav.hasMore(); tnav++ )
657  {
658  // first copy source to tmp because convolveLine() cannot work in-place
659  copyLine(tnav.begin(), tnav.end(), acc, tmpline.begin(), acc );
660 
661  convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
662  destIter( tnav.begin() + lstart, acc ),
663  kernel1d( kit[axisorder[d]] ), lstart, lstop);
664  }
665 
666  dstart[axisorder[d]] = lstart;
667  dstop[axisorder[d]] = lstop;
668  }
669 
670  copyMultiArray(tmp.traverser_begin()+dstart, stop-start, acc, di, dest);
671 }
672 
673 
674 template <class K>
675 void
676 scaleKernel(K & kernel, double a)
677 {
678  for(int i = kernel.left(); i <= kernel.right(); ++i)
679  kernel[i] = detail::RequiresExplicitCast<typename K::value_type>::cast(kernel[i] * a);
680 }
681 
682 
683 } // namespace detail
684 
685 /** \addtogroup MultiArrayConvolutionFilters Convolution filters for multi-dimensional arrays.
686 
687  These functions realize a separable convolution on an arbitrary dimensional
688  array that is specified by iterators (compatible to \ref MultiIteratorPage)
689  and shape objects. It can therefore be applied to a wide range of data structures
690  (\ref vigra::MultiArrayView, \ref vigra::MultiArray etc.).
691 */
692 //@{
693 
694 /********************************************************/
695 /* */
696 /* separableConvolveMultiArray */
697 /* */
698 /********************************************************/
699 
700 /** \brief Separated convolution on multi-dimensional arrays.
701 
702  This function computes a separated convolution on all dimensions
703  of the given multi-dimensional array. Both source and destination
704  arrays are represented by iterators, shape objects and accessors.
705  The destination array is required to already have the correct size.
706 
707  There are two variants of this functions: one takes a single kernel
708  of type \ref vigra::Kernel1D which is then applied to all dimensions,
709  whereas the other requires an iterator referencing a sequence of
710  \ref vigra::Kernel1D objects, one for every dimension of the data.
711  Then the first kernel in this sequence is applied to the innermost
712  dimension (e.g. the x-axis of an image), while the last is applied to the
713  outermost dimension (e.g. the z-axis in a 3D image).
714 
715  This function may work in-place, which means that <tt>source.data() == dest.data()</tt> is allowed.
716  A full-sized internal array is only allocated if working on the destination
717  array directly would cause round-off errors (i.e. if
718  <tt>typeid(typename NumericTraits<T2>::RealPromote) != typeid(T2)</tt>).
719 
720  If <tt>start</tt> and <tt>stop</tt> have non-default values, they must represent
721  a valid subarray of the input array. The convolution is then restricted to that
722  subarray, and it is assumed that the output array only refers to the
723  subarray (i.e. <tt>dest.shape() == stop - start</tt>). Negative ROI boundaries are
724  interpreted relative to the end of the respective dimension
725  (i.e. <tt>if(stop[k] < 0) stop[k] += source.shape(k);</tt>).
726 
727  <b> Declarations:</b>
728 
729  pass arbitrary-dimensional array views:
730  \code
731  namespace vigra {
732  // apply each kernel from the sequence 'kernels' in turn
733  template <unsigned int N, class T1, class S1,
734  class T2, class S2,
735  class KernelIterator>
736  void
737  separableConvolveMultiArray(MultiArrayView<N, T1, S1> const & source,
738  MultiArrayView<N, T2, S2> dest,
739  KernelIterator kernels,
740  typename MultiArrayShape<N>::type start = typename MultiArrayShape<N>::type(),
741  typename MultiArrayShape<N>::type stop = typename MultiArrayShape<N>::type());
742 
743  // apply the same kernel to all dimensions
744  template <unsigned int N, class T1, class S1,
745  class T2, class S2,
746  class T>
747  void
748  separableConvolveMultiArray(MultiArrayView<N, T1, S1> const & source,
749  MultiArrayView<N, T2, S2> dest,
750  Kernel1D<T> const & kernel,
751  typename MultiArrayShape<N>::type const & start = typename MultiArrayShape<N>::type(),
752  typename MultiArrayShape<N>::type const & stop = typename MultiArrayShape<N>::type());
753  }
754  \endcode
755 
756  \deprecatedAPI{separableConvolveMultiArray}
757  pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
758  \code
759  namespace vigra {
760  // apply the same kernel to all dimensions
761  template <class SrcIterator, class SrcShape, class SrcAccessor,
762  class DestIterator, class DestAccessor, class T>
763  void
764  separableConvolveMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
765  DestIterator diter, DestAccessor dest,
766  Kernel1D<T> const & kernel,
767  SrcShape const & start = SrcShape(),
768  SrcShape const & stop = SrcShape());
769 
770  // apply each kernel from the sequence 'kernels' in turn
771  template <class SrcIterator, class SrcShape, class SrcAccessor,
772  class DestIterator, class DestAccessor, class KernelIterator>
773  void
774  separableConvolveMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
775  DestIterator diter, DestAccessor dest,
776  KernelIterator kernels,
777  SrcShape const & start = SrcShape(),
778  SrcShape const & stop = SrcShape());
779  }
780  \endcode
781  use argument objects in conjunction with \ref ArgumentObjectFactories :
782  \code
783  namespace vigra {
784  // apply the same kernel to all dimensions
785  template <class SrcIterator, class SrcShape, class SrcAccessor,
786  class DestIterator, class DestAccessor, class T>
787  void
788  separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
789  pair<DestIterator, DestAccessor> const & dest,
790  Kernel1D<T> const & kernel,
791  SrcShape const & start = SrcShape(),
792  SrcShape const & stop = SrcShape());
793 
794  // apply each kernel from the sequence 'kernels' in turn
795  template <class SrcIterator, class SrcShape, class SrcAccessor,
796  class DestIterator, class DestAccessor, class KernelIterator>
797  void
798  separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
799  pair<DestIterator, DestAccessor> const & dest,
800  KernelIterator kernels,
801  SrcShape const & start = SrcShape(),
802  SrcShape const & stop = SrcShape());
803  }
804  \endcode
805  \deprecatedEnd
806 
807  <b> Usage:</b>
808 
809  <b>\#include</b> <vigra/multi_convolution.hxx><br/>
810  Namespace: vigra
811 
812  \code
813  Shape3 shape(width, height, depth);
814  MultiArray<3, unsigned char> source(shape);
815  MultiArray<3, float> dest(shape);
816  ...
817  Kernel1D<float> gauss;
818  gauss.initGaussian(sigma);
819 
820  // smooth all dimensions with the same kernel
821  separableConvolveMultiArray(source, dest, gauss);
822 
823  // create 3 Gauss kernels, one for each dimension, but smooth the z-axis less
824  ArrayVector<Kernel1D<float> > kernels(3, gauss);
825  kernels[2].initGaussian(sigma / 2.0);
826 
827  // perform Gaussian smoothing on all dimensions
828  separableConvolveMultiArray(source, dest, kernels.begin());
829 
830  // create output array for a ROI
831  MultiArray<3, float> destROI(shape - Shape3(10,10,10));
832 
833  // only smooth the given ROI (ignore 5 pixels on all sides of the array)
834  separableConvolveMultiArray(source, destROI, gauss, Shape3(5,5,5), Shape3(-5,-5,-5));
835  \endcode
836 
837  \deprecatedUsage{separableConvolveMultiArray}
838  \code
839  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
840  MultiArray<3, unsigned char> source(shape);
841  MultiArray<3, float> dest(shape);
842  ...
843  Kernel1D<float> gauss;
844  gauss.initGaussian(sigma);
845  // create 3 Gauss kernels, one for each dimension
846  ArrayVector<Kernel1D<float> > kernels(3, gauss);
847 
848  // perform Gaussian smoothing on all dimensions
849  separableConvolveMultiArray(source, dest,
850  kernels.begin());
851  \endcode
852  <b> Required Interface:</b>
853  \code
854  see \ref separableConvolveImage(), in addition:
855 
856  NumericTraits<T1>::RealPromote s = src[0];
857 
858  s = s + s;
859  s = kernel(0) * s;
860  \endcode
861  \deprecatedEnd
862 
863  \see vigra::Kernel1D, convolveLine()
864 */
866 
867 template <class SrcIterator, class SrcShape, class SrcAccessor,
868  class DestIterator, class DestAccessor, class KernelIterator>
869 void
870 separableConvolveMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
871  DestIterator d, DestAccessor dest,
872  KernelIterator kernels,
873  SrcShape start = SrcShape(),
874  SrcShape stop = SrcShape())
875 {
876  typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
877 
878 
879  if(stop != SrcShape())
880  {
881 
882  enum { N = 1 + SrcIterator::level };
883  detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, start);
884  detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, stop);
885 
886  for(int k=0; k<N; ++k)
887  vigra_precondition(0 <= start[k] && start[k] < stop[k] && stop[k] <= shape[k],
888  "separableConvolveMultiArray(): invalid subarray shape.");
889 
890  detail::internalSeparableConvolveSubarray(s, shape, src, d, dest, kernels, start, stop);
891  }
892  else if(!IsSameType<TmpType, typename DestAccessor::value_type>::boolResult)
893  {
894  // need a temporary array to avoid rounding errors
896  detail::internalSeparableConvolveMultiArrayTmp( s, shape, src,
897  tmpArray.traverser_begin(), typename AccessorTraits<TmpType>::default_accessor(), kernels );
898  copyMultiArray(srcMultiArrayRange(tmpArray), destIter(d, dest));
899  }
900  else
901  {
902  // work directly on the destination array
903  detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, d, dest, kernels );
904  }
905 }
906 
907 template <class SrcIterator, class SrcShape, class SrcAccessor,
908  class DestIterator, class DestAccessor, class T>
909 inline void
910 separableConvolveMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
911  DestIterator d, DestAccessor dest,
912  Kernel1D<T> const & kernel,
913  SrcShape const & start = SrcShape(),
914  SrcShape const & stop = SrcShape())
915 {
916  ArrayVector<Kernel1D<T> > kernels(shape.size(), kernel);
917 
918  separableConvolveMultiArray( s, shape, src, d, dest, kernels.begin(), start, stop);
919 }
920 
921 template <class SrcIterator, class SrcShape, class SrcAccessor,
922  class DestIterator, class DestAccessor, class KernelIterator>
923 inline void
924 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
925  pair<DestIterator, DestAccessor> const & dest,
926  KernelIterator kit,
927  SrcShape const & start = SrcShape(),
928  SrcShape const & stop = SrcShape())
929 {
930  separableConvolveMultiArray( source.first, source.second, source.third,
931  dest.first, dest.second, kit, start, stop );
932 }
933 
934 template <class SrcIterator, class SrcShape, class SrcAccessor,
935  class DestIterator, class DestAccessor, class T>
936 inline void
937 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
938  pair<DestIterator, DestAccessor> const & dest,
939  Kernel1D<T> const & kernel,
940  SrcShape const & start = SrcShape(),
941  SrcShape const & stop = SrcShape())
942 {
943  ArrayVector<Kernel1D<T> > kernels(source.second.size(), kernel);
944 
945  separableConvolveMultiArray( source.first, source.second, source.third,
946  dest.first, dest.second, kernels.begin(), start, stop);
947 }
948 
949 template <unsigned int N, class T1, class S1,
950  class T2, class S2,
951  class KernelIterator>
952 inline void
955  KernelIterator kit,
956  typename MultiArrayShape<N>::type start = typename MultiArrayShape<N>::type(),
957  typename MultiArrayShape<N>::type stop = typename MultiArrayShape<N>::type())
958 {
959  if(stop != typename MultiArrayShape<N>::type())
960  {
961  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), start);
962  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), stop);
963  vigra_precondition(dest.shape() == (stop - start),
964  "separableConvolveMultiArray(): shape mismatch between ROI and output.");
965  }
966  else
967  {
968  vigra_precondition(source.shape() == dest.shape(),
969  "separableConvolveMultiArray(): shape mismatch between input and output.");
970  }
971  separableConvolveMultiArray( srcMultiArrayRange(source),
972  destMultiArray(dest), kit, start, stop );
973 }
974 
975 template <unsigned int N, class T1, class S1,
976  class T2, class S2,
977  class T>
978 inline void
981  Kernel1D<T> const & kernel,
982  typename MultiArrayShape<N>::type const & start = typename MultiArrayShape<N>::type(),
983  typename MultiArrayShape<N>::type const & stop = typename MultiArrayShape<N>::type())
984 {
985  ArrayVector<Kernel1D<T> > kernels(N, kernel);
986  separableConvolveMultiArray(source, dest, kernels.begin(), start, stop);
987 }
988 
989 /********************************************************/
990 /* */
991 /* convolveMultiArrayOneDimension */
992 /* */
993 /********************************************************/
994 
995 /** \brief Convolution along a single dimension of a multi-dimensional arrays.
996 
997  This function computes a convolution along one dimension (specified by
998  the parameter <tt>dim</tt> of the given multi-dimensional array with the given
999  <tt>kernel</tt>. The destination array must already have the correct size.
1000 
1001  If <tt>start</tt> and <tt>stop</tt> have non-default values, they must represent
1002  a valid subarray of the input array. The convolution is then restricted to that
1003  subarray, and it is assumed that the output array only refers to the
1004  subarray (i.e. <tt>dest.shape() == stop - start</tt>). Negative ROI boundaries are
1005  interpreted relative to the end of the respective dimension
1006  (i.e. <tt>if(stop[k] < 0) stop[k] += source.shape(k);</tt>).
1007 
1008  This function may work in-place, which means that <tt>source.data() == dest.data()</tt> is allowed.
1009 
1010  <b> Declarations:</b>
1011 
1012  pass arbitrary-dimensional array views:
1013  \code
1014  namespace vigra {
1015  template <unsigned int N, class T1, class S1,
1016  class T2, class S2,
1017  class T>
1018  void
1019  convolveMultiArrayOneDimension(MultiArrayView<N, T1, S1> const & source,
1020  MultiArrayView<N, T2, S2> dest,
1021  unsigned int dim,
1022  Kernel1D<T> const & kernel,
1023  typename MultiArrayShape<N>::type start = typename MultiArrayShape<N>::type(),
1024  typename MultiArrayShape<N>::type stop = typename MultiArrayShape<N>::type());
1025  }
1026  \endcode
1027 
1028  \deprecatedAPI{convolveMultiArrayOneDimension}
1029  pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
1030  \code
1031  namespace vigra {
1032  template <class SrcIterator, class SrcShape, class SrcAccessor,
1033  class DestIterator, class DestAccessor, class T>
1034  void
1035  convolveMultiArrayOneDimension(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1036  DestIterator diter, DestAccessor dest,
1037  unsigned int dim, vigra::Kernel1D<T> const & kernel,
1038  SrcShape const & start = SrcShape(),
1039  SrcShape const & stop = SrcShape());
1040  }
1041  \endcode
1042  use argument objects in conjunction with \ref ArgumentObjectFactories :
1043  \code
1044  namespace vigra {
1045  template <class SrcIterator, class SrcShape, class SrcAccessor,
1046  class DestIterator, class DestAccessor, class T>
1047  void
1048  convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1049  pair<DestIterator, DestAccessor> const & dest,
1050  unsigned int dim, vigra::Kernel1D<T> const & kernel,
1051  SrcShape const & start = SrcShape(),
1052  SrcShape const & stop = SrcShape());
1053  }
1054  \endcode
1055  \deprecatedEnd
1056 
1057  <b> Usage:</b>
1058 
1059  <b>\#include</b> <vigra/multi_convolution.hxx><br/>
1060  Namespace: vigra
1061 
1062  \code
1063  Shape3 shape(width, height, depth);
1064  MultiArray<3, unsigned char> source(shape);
1065  MultiArray<3, float> dest(shape);
1066  ...
1067  Kernel1D<float> gauss;
1068  gauss.initGaussian(sigma);
1069 
1070  // perform Gaussian smoothing along dimension 1 (height)
1071  convolveMultiArrayOneDimension(source, dest, 1, gauss);
1072  \endcode
1073 
1074  \see separableConvolveMultiArray()
1075 */
1077 
1078 template <class SrcIterator, class SrcShape, class SrcAccessor,
1079  class DestIterator, class DestAccessor, class T>
1080 void
1081 convolveMultiArrayOneDimension(SrcIterator s, SrcShape const & shape, SrcAccessor src,
1082  DestIterator d, DestAccessor dest,
1083  unsigned int dim, vigra::Kernel1D<T> const & kernel,
1084  SrcShape const & start = SrcShape(),
1085  SrcShape const & stop = SrcShape())
1086 {
1087  enum { N = 1 + SrcIterator::level };
1088  vigra_precondition( dim < N,
1089  "convolveMultiArrayOneDimension(): The dimension number to convolve must be smaller "
1090  "than the data dimensionality" );
1091 
1092  typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
1093  typedef typename AccessorTraits<TmpType>::default_const_accessor TmpAccessor;
1094  ArrayVector<TmpType> tmp( shape[dim] );
1095 
1096  typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
1097  typedef MultiArrayNavigator<DestIterator, N> DNavigator;
1098 
1099  SrcShape sstart, sstop(shape), dstart, dstop(shape);
1100 
1101  if(stop != SrcShape())
1102  {
1103  sstart = start;
1104  sstop = stop;
1105  sstart[dim] = 0;
1106  sstop[dim] = shape[dim];
1107  dstop = stop - start;
1108  }
1109 
1110  SNavigator snav( s, sstart, sstop, dim );
1111  DNavigator dnav( d, dstart, dstop, dim );
1112 
1113  for( ; snav.hasMore(); snav++, dnav++ )
1114  {
1115  // first copy source to temp for maximum cache efficiency
1116  copyLine(snav.begin(), snav.end(), src,
1118 
1119  convolveLine(srcIterRange( tmp.begin(), tmp.end(), TmpAccessor()),
1120  destIter( dnav.begin(), dest ),
1121  kernel1d( kernel), start[dim], stop[dim]);
1122  }
1123 }
1124 
1125 template <class SrcIterator, class SrcShape, class SrcAccessor,
1126  class DestIterator, class DestAccessor, class T>
1127 inline void
1128 convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1129  pair<DestIterator, DestAccessor> const & dest,
1130  unsigned int dim,
1131  Kernel1D<T> const & kernel,
1132  SrcShape const & start = SrcShape(),
1133  SrcShape const & stop = SrcShape())
1134 {
1135  convolveMultiArrayOneDimension(source.first, source.second, source.third,
1136  dest.first, dest.second, dim, kernel, start, stop);
1137 }
1138 
1139 template <unsigned int N, class T1, class S1,
1140  class T2, class S2,
1141  class T>
1142 inline void
1145  unsigned int dim,
1146  Kernel1D<T> const & kernel,
1147  typename MultiArrayShape<N>::type start = typename MultiArrayShape<N>::type(),
1148  typename MultiArrayShape<N>::type stop = typename MultiArrayShape<N>::type())
1149 {
1150  if(stop != typename MultiArrayShape<N>::type())
1151  {
1152  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), start);
1153  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), stop);
1154  vigra_precondition(dest.shape() == (stop - start),
1155  "convolveMultiArrayOneDimension(): shape mismatch between ROI and output.");
1156  }
1157  else
1158  {
1159  vigra_precondition(source.shape() == dest.shape(),
1160  "convolveMultiArrayOneDimension(): shape mismatch between input and output.");
1161  }
1162  convolveMultiArrayOneDimension(srcMultiArrayRange(source),
1163  destMultiArray(dest), dim, kernel, start, stop);
1164 }
1165 
1166 /********************************************************/
1167 /* */
1168 /* gaussianSmoothMultiArray */
1169 /* */
1170 /********************************************************/
1171 
1172 /** \brief Isotropic Gaussian smoothing of a multi-dimensional arrays.
1173 
1174  This function computes an isotropic convolution of the given N-dimensional
1175  array with a Gaussian filter at the given standard deviation <tt>sigma</tt>.
1176  Both source and destination arrays are represented by
1177  iterators, shape objects and accessors. The destination array is required to
1178  already have the correct size. This function may work in-place, which means
1179  that <tt>source.data() == dest.data()</tt> is allowed. It is implemented by a call to
1180  \ref separableConvolveMultiArray() with the appropriate kernel.
1181 
1182  Anisotropic data should be provided with appropriate \ref vigra::ConvolutionOptions
1183  to adjust the filter sizes for the resolution of each axis.
1184  Otherwise, the parameter <tt>opt</tt> is optional unless the parameter
1185  <tt>sigma</tt> is omitted.
1186 
1187  If you pass \ref vigra::BlockwiseConvolutionOptions instead, the algorithm will
1188  be executed in parallel on data blocks of a certain size. The block size can be
1189  customized via <tt>BlockwiseConvolutionOptions::blockShape()</tt>, but the defaults
1190  usually work reasonably. By default, the number of threads equals the capabilities
1191  of your hardware, but you can change this via <tt>BlockwiseConvolutionOptions::numThreads()</tt>.
1192 
1193  <b> Declarations:</b>
1194 
1195  pass arbitrary-dimensional array views:
1196  \code
1197  namespace vigra {
1198  // pass filter scale explicitly
1199  template <unsigned int N, class T1, class S1,
1200  class T2, class S2>
1201  void
1202  gaussianSmoothMultiArray(MultiArrayView<N, T1, S1> const & source,
1203  MultiArrayView<N, T2, S2> dest,
1204  double sigma,
1205  ConvolutionOptions<N> opt = ConvolutionOptions<N>());
1206 
1207  // pass filer scale(s) in the option object
1208  template <unsigned int N, class T1, class S1,
1209  class T2, class S2>
1210  void
1211  gaussianSmoothMultiArray(MultiArrayView<N, T1, S1> const & source,
1212  MultiArrayView<N, T2, S2> dest,
1213  ConvolutionOptions<N> opt);
1214 
1215  // as above, but execute algorirhm in parallel
1216  template <unsigned int N, class T1, class S1,
1217  class T2, class S2>
1218  void
1219  gaussianSmoothMultiArray(MultiArrayView<N, T1, S1> const & source,
1220  MultiArrayView<N, T2, S2> dest,
1221  BlockwiseConvolutionOptions<N> opt);
1222  }
1223  \endcode
1224 
1225  \deprecatedAPI{gaussianSmoothMultiArray}
1226  pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
1227  \code
1228  namespace vigra {
1229  template <class SrcIterator, class SrcShape, class SrcAccessor,
1230  class DestIterator, class DestAccessor>
1231  void
1232  gaussianSmoothMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1233  DestIterator diter, DestAccessor dest,
1234  double sigma,
1235  const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1236  }
1237  \endcode
1238  use argument objects in conjunction with \ref ArgumentObjectFactories :
1239  \code
1240  namespace vigra {
1241  template <class SrcIterator, class SrcShape, class SrcAccessor,
1242  class DestIterator, class DestAccessor>
1243  void
1244  gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1245  pair<DestIterator, DestAccessor> const & dest,
1246  double sigma,
1247  const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1248  }
1249  \endcode
1250  \deprecatedEnd
1251 
1252  <b> Usage:</b>
1253 
1254  <b>\#include</b> <vigra/multi_convolution.hxx> (sequential version)<br/>
1255  <b>\#include</b> <vigra/multi_blockwise.hxx> (parallel version)<br/>
1256  Namespace: vigra
1257 
1258  \code
1259  Shape3 shape(width, height, depth);
1260  MultiArray<3, unsigned char> source(shape);
1261  MultiArray<3, float> dest(shape);
1262  ...
1263  // perform isotropic Gaussian smoothing at scale 'sigma'
1264  gaussianSmoothMultiArray(source, dest, sigma);
1265  \endcode
1266 
1267  <b> Multi-threaded execution:</b>
1268 
1269  \code
1270  Shape3 shape(width, height, depth);
1271  MultiArray<3, unsigned char> source(shape);
1272  MultiArray<3, float> dest(shape);
1273  ...
1274  BlockwiseConvolutionOptions<3> opt;
1275  opt.numThreads(4); // use 4 threads (uses hardware default if not given)
1276  opt.innerScale(sigma);
1277 
1278  // perform isotropic Gaussian smoothing at scale 'sigma' in parallel
1279  gaussianSmoothMultiArray(source, dest, sigma, opt);
1280  \endcode
1281 
1282  <b> Usage with anisotropic data:</b>
1283 
1284  \code
1285  Shape3 shape(width, height, depth);
1286  MultiArray<3, unsigned char> source(shape);
1287  MultiArray<3, float> dest(shape);
1288  TinyVector<float, 3> step_size;
1289  TinyVector<float, 3> resolution_sigmas;
1290  ...
1291  // perform anisotropic Gaussian smoothing at scale 'sigma'
1292  gaussianSmoothMultiArray(source, dest, sigma,
1293  ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
1294  \endcode
1295 
1296  \see separableConvolveMultiArray()
1297 */
1299 
1300 template <class SrcIterator, class SrcShape, class SrcAccessor,
1301  class DestIterator, class DestAccessor>
1302 void
1303 gaussianSmoothMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
1304  DestIterator d, DestAccessor dest,
1306  const char *const function_name = "gaussianSmoothMultiArray" )
1307 {
1308  static const int N = SrcShape::static_size;
1309 
1310  typename ConvolutionOptions<N>::ScaleIterator params = opt.scaleParams();
1311  ArrayVector<Kernel1D<double> > kernels(N);
1312 
1313  for (int dim = 0; dim < N; ++dim, ++params)
1314  kernels[dim].initGaussian(params.sigma_scaled(function_name, true),
1315  1.0, opt.window_ratio);
1316 
1317  separableConvolveMultiArray(s, shape, src, d, dest, kernels.begin(), opt.from_point, opt.to_point);
1318 }
1319 
1320 template <class SrcIterator, class SrcShape, class SrcAccessor,
1321  class DestIterator, class DestAccessor>
1322 inline void
1323 gaussianSmoothMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
1324  DestIterator d, DestAccessor dest, double sigma,
1326 {
1328  gaussianSmoothMultiArray(s, shape, src, d, dest, par.stdDev(sigma));
1329 }
1330 
1331 template <class SrcIterator, class SrcShape, class SrcAccessor,
1332  class DestIterator, class DestAccessor>
1333 inline void
1334 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1335  pair<DestIterator, DestAccessor> const & dest,
1337 {
1338  gaussianSmoothMultiArray( source.first, source.second, source.third,
1339  dest.first, dest.second, opt );
1340 }
1341 
1342 template <class SrcIterator, class SrcShape, class SrcAccessor,
1343  class DestIterator, class DestAccessor>
1344 inline void
1345 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1346  pair<DestIterator, DestAccessor> const & dest, double sigma,
1348 {
1349  gaussianSmoothMultiArray( source.first, source.second, source.third,
1350  dest.first, dest.second, sigma, opt );
1351 }
1352 
1353 template <unsigned int N, class T1, class S1,
1354  class T2, class S2>
1355 inline void
1359 {
1360  if(opt.to_point != typename MultiArrayShape<N>::type())
1361  {
1362  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
1363  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
1364  vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1365  "gaussianSmoothMultiArray(): shape mismatch between ROI and output.");
1366  }
1367  else
1368  {
1369  vigra_precondition(source.shape() == dest.shape(),
1370  "gaussianSmoothMultiArray(): shape mismatch between input and output.");
1371  }
1372 
1373  gaussianSmoothMultiArray( srcMultiArrayRange(source),
1374  destMultiArray(dest), opt );
1375 }
1376 
1377 template <unsigned int N, class T1, class S1,
1378  class T2, class S2>
1379 inline void
1382  double sigma,
1384 {
1385  gaussianSmoothMultiArray( source, dest, opt.stdDev(sigma) );
1386 }
1387 
1388 
1389 /********************************************************/
1390 /* */
1391 /* gaussianGradientMultiArray */
1392 /* */
1393 /********************************************************/
1394 
1395 /** \brief Calculate Gaussian gradient of a multi-dimensional arrays.
1396 
1397  This function computes the Gaussian gradient of the given N-dimensional
1398  array with a sequence of first-derivative-of-Gaussian filters at the given
1399  standard deviation <tt>sigma</tt> (differentiation is applied to each dimension
1400  in turn, starting with the innermost dimension). The destination array is
1401  required to have a vector valued pixel type with as many elements as the number of
1402  dimensions. This function is implemented by calls to
1403  \ref separableConvolveMultiArray() with the appropriate kernels.
1404 
1405  Anisotropic data should be provided with appropriate \ref vigra::ConvolutionOptions
1406  to adjust the filter sizes for the resolution of each axis.
1407  Otherwise, the parameter <tt>opt</tt> is optional unless the parameter
1408  <tt>sigma</tt> is omitted.
1409 
1410  If you pass \ref vigra::BlockwiseConvolutionOptions instead, the algorithm will
1411  be executed in parallel on data blocks of a certain size. The block size can be
1412  customized via <tt>BlockwiseConvolutionOptions::blockShape()</tt>, but the defaults
1413  usually work reasonably. By default, the number of threads equals the capabilities
1414  of your hardware, but you can change this via <tt>BlockwiseConvolutionOptions::numThreads()</tt>.
1415 
1416  <b> Declarations:</b>
1417 
1418  pass arbitrary-dimensional array views:
1419  \code
1420  namespace vigra {
1421  // pass filter scale explicitly
1422  template <unsigned int N, class T1, class S1,
1423  class T2, class S2>
1424  void
1425  gaussianGradientMultiArray(MultiArrayView<N, T1, S1> const & source,
1426  MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1427  double sigma,
1428  ConvolutionOptions<N> opt = ConvolutionOptions<N>());
1429 
1430  // pass filter scale(s) in option object
1431  template <unsigned int N, class T1, class S1,
1432  class T2, class S2>
1433  void
1434  gaussianGradientMultiArray(MultiArrayView<N, T1, S1> const & source,
1435  MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1436  ConvolutionOptions<N> opt);
1437 
1438  // likewise, but execute algorithm in parallel
1439  template <unsigned int N, class T1, class S1,
1440  class T2, class S2>
1441  void
1442  gaussianGradientMultiArray(MultiArrayView<N, T1, S1> const & source,
1443  MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1444  BlockwiseConvolutionOptions<N> opt);
1445  }
1446  \endcode
1447 
1448  \deprecatedAPI{gaussianGradientMultiArray}
1449  pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
1450  \code
1451  namespace vigra {
1452  template <class SrcIterator, class SrcShape, class SrcAccessor,
1453  class DestIterator, class DestAccessor>
1454  void
1455  gaussianGradientMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1456  DestIterator diter, DestAccessor dest,
1457  double sigma,
1458  const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1459  }
1460  \endcode
1461  use argument objects in conjunction with \ref ArgumentObjectFactories :
1462  \code
1463  namespace vigra {
1464  template <class SrcIterator, class SrcShape, class SrcAccessor,
1465  class DestIterator, class DestAccessor>
1466  void
1467  gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1468  pair<DestIterator, DestAccessor> const & dest,
1469  double sigma,
1470  const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1471  }
1472  \endcode
1473  \deprecatedEnd
1474 
1475  <b> Usage:</b>
1476 
1477  <b>\#include</b> <vigra/multi_convolution.hxx> (sequential version)<br/>
1478  <b>\#include</b> <vigra/multi_blockwise.hxx> (parallel version)<br/>
1479  Namespace: vigra
1480 
1481  \code
1482  Shape3 shape(width, height, depth);
1483  MultiArray<3, unsigned char> source(shape);
1484  MultiArray<3, TinyVector<float, 3> > dest(shape);
1485  ...
1486  // compute Gaussian gradient at scale sigma
1487  gaussianGradientMultiArray(source, dest, sigma);
1488  \endcode
1489 
1490  <b> Usage with anisotropic data:</b>
1491 
1492  \code
1493  Shape3 shape(width, height, depth);
1494  MultiArray<3, unsigned char> source(shape);
1495  MultiArray<3, TinyVector<float, 3> > dest(shape);
1496  TinyVector<float, 3> step_size;
1497  TinyVector<float, 3> resolution_sigmas;
1498  ...
1499  // compute Gaussian gradient at scale sigma
1500  gaussianGradientMultiArray(source, dest, sigma,
1501  ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
1502  \endcode
1503 
1504  \see separableConvolveMultiArray()
1505 */
1507 
1508 template <class SrcIterator, class SrcShape, class SrcAccessor,
1509  class DestIterator, class DestAccessor>
1510 void
1511 gaussianGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1512  DestIterator di, DestAccessor dest,
1514  const char *const function_name = "gaussianGradientMultiArray")
1515 {
1516  typedef typename DestAccessor::value_type DestType;
1517  typedef typename DestType::value_type DestValueType;
1518  typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1519 
1520  static const int N = SrcShape::static_size;
1521  typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
1522 
1523  for(int k=0; k<N; ++k)
1524  if(shape[k] <=0)
1525  return;
1526 
1527  vigra_precondition(N == (int)dest.size(di),
1528  "gaussianGradientMultiArray(): Wrong number of channels in output array.");
1529 
1530  ParamType params = opt.scaleParams();
1531  ParamType params2(params);
1532 
1533  ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
1534  for (int dim = 0; dim < N; ++dim, ++params)
1535  {
1536  double sigma = params.sigma_scaled(function_name);
1537  plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
1538  }
1539 
1540  typedef VectorElementAccessor<DestAccessor> ElementAccessor;
1541 
1542  // compute gradient components
1543  for (int dim = 0; dim < N; ++dim, ++params2)
1544  {
1545  ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
1546  kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 1, 1.0, opt.window_ratio);
1547  detail::scaleKernel(kernels[dim], 1.0 / params2.step_size());
1548  separableConvolveMultiArray(si, shape, src, di, ElementAccessor(dim, dest), kernels.begin(),
1549  opt.from_point, opt.to_point);
1550  }
1551 }
1552 
1553 template <class SrcIterator, class SrcShape, class SrcAccessor,
1554  class DestIterator, class DestAccessor>
1555 void
1556 gaussianGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1557  DestIterator di, DestAccessor dest, double sigma,
1559 {
1560  gaussianGradientMultiArray(si, shape, src, di, dest, opt.stdDev(sigma));
1561 }
1562 
1563 template <class SrcIterator, class SrcShape, class SrcAccessor,
1564  class DestIterator, class DestAccessor>
1565 inline void
1566 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1567  pair<DestIterator, DestAccessor> const & dest,
1569 {
1570  gaussianGradientMultiArray( source.first, source.second, source.third,
1571  dest.first, dest.second, opt );
1572 }
1573 
1574 template <class SrcIterator, class SrcShape, class SrcAccessor,
1575  class DestIterator, class DestAccessor>
1576 inline void
1577 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1578  pair<DestIterator, DestAccessor> const & dest,
1579  double sigma,
1581 {
1582  gaussianGradientMultiArray( source.first, source.second, source.third,
1583  dest.first, dest.second, sigma, opt );
1584 }
1585 
1586 template <unsigned int N, class T1, class S1,
1587  class T2, class S2>
1588 inline void
1590  MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1591  ConvolutionOptions<N> opt )
1592 {
1593  if(opt.to_point != typename MultiArrayShape<N>::type())
1594  {
1595  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
1596  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
1597  vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1598  "gaussianGradientMultiArray(): shape mismatch between ROI and output.");
1599  }
1600  else
1601  {
1602  vigra_precondition(source.shape() == dest.shape(),
1603  "gaussianGradientMultiArray(): shape mismatch between input and output.");
1604  }
1605 
1606  gaussianGradientMultiArray( srcMultiArrayRange(source),
1607  destMultiArray(dest), opt );
1608 }
1609 
1610 template <unsigned int N, class T1, class S1,
1611  class T2, class S2>
1612 inline void
1614  MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1615  double sigma,
1617 {
1618  gaussianGradientMultiArray( source, dest, opt.stdDev(sigma) );
1619 }
1620 
1621 /********************************************************/
1622 /* */
1623 /* gaussianGradientMagnitude */
1624 /* */
1625 /********************************************************/
1626 
1627 namespace detail {
1628 
1629 template <unsigned int N, class T1, class S1,
1630  class T2, class S2>
1631 void
1632 gaussianGradientMagnitudeImpl(MultiArrayView<N+1, T1, S1> const & src,
1635 {
1636  typename MultiArrayShape<N>::type shape(src.shape().template subarray<0,N>());
1637  if(opt.to_point != typename MultiArrayShape<N>::type())
1638  {
1639  detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.from_point);
1640  detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.to_point);
1641  vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1642  "gaussianGradientMagnitude(): shape mismatch between ROI and output.");
1643  }
1644  else
1645  {
1646  vigra_precondition(shape == dest.shape(),
1647  "gaussianGradientMagnitude(): shape mismatch between input and output.");
1648  }
1649 
1650  dest.init(0.0);
1651 
1652  typedef typename NumericTraits<T1>::RealPromote TmpType;
1654 
1655  using namespace multi_math;
1656 
1657  for(int k=0; k<src.shape(N); ++k)
1658  {
1659  gaussianGradientMultiArray(src.bindOuter(k), grad, opt);
1660 
1661  dest += squaredNorm(grad);
1662  }
1663  dest = sqrt(dest);
1664 }
1665 
1666 } // namespace detail
1667 
1668  // documentation is in convolution.hxx
1669 template <unsigned int N, class T1, class S1,
1670  class T2, class S2>
1671 inline void
1672 gaussianGradientMagnitude(MultiArrayView<N+1, Multiband<T1>, S1> const & src,
1674  ConvolutionOptions<N> const & opt)
1675 {
1676  detail::gaussianGradientMagnitudeImpl<N, T1>(src, dest, opt);
1677 }
1678 
1679 template <unsigned int N, class T1, class S1,
1680  class T2, class S2>
1681 inline void
1684  ConvolutionOptions<N> const & opt)
1685 {
1686  detail::gaussianGradientMagnitudeImpl<N, T1>(src.insertSingletonDimension(N), dest, opt);
1687 }
1688 
1689 template <unsigned int N, class T1, int M, class S1,
1690  class T2, class S2>
1691 inline void
1694  ConvolutionOptions<N> const & opt)
1695 {
1696  detail::gaussianGradientMagnitudeImpl<N, T1>(src.expandElements(N), dest, opt);
1697 }
1698 
1699 template <unsigned int N, class T1, unsigned int R, unsigned int G, unsigned int B, class S1,
1700  class T2, class S2>
1701 inline void
1704  ConvolutionOptions<N> const & opt)
1705 {
1706  detail::gaussianGradientMagnitudeImpl<N, T1>(src.expandElements(N), dest, opt);
1707 }
1708 
1709 template <unsigned int N, class T1, class S1,
1710  class T2, class S2>
1711 inline void
1714  double sigma,
1716 {
1717  gaussianGradientMagnitude(src, dest, opt.stdDev(sigma));
1718 }
1719 
1720 template <unsigned int N, class T1, class S1,
1721  class T2, class S2>
1722 inline void
1723 gaussianGradientMagnitude(MultiArrayView<N+1, Multiband<T1>, S1> const & src,
1725  double sigma,
1727 {
1728  gaussianGradientMagnitude<N>(src, dest, opt.stdDev(sigma));
1729 }
1730 
1731 /********************************************************/
1732 /* */
1733 /* symmetricGradientMultiArray */
1734 /* */
1735 /********************************************************/
1736 
1737 /** \brief Calculate gradient of a multi-dimensional arrays using symmetric difference filters.
1738 
1739  This function computes the gradient of the given N-dimensional
1740  array with a sequence of symmetric difference filters a (differentiation is applied
1741  to each dimension in turn, starting with the innermost dimension).
1742  The destination array is required to have a vector valued pixel type with as many
1743  elements as the number of dimensions. This function is implemented by calls to
1744  \ref convolveMultiArrayOneDimension() with the symmetric difference kernel.
1745 
1746  Anisotropic data should be provided with appropriate \ref vigra::ConvolutionOptions
1747  to adjust the filter sizes for the resolution of each axis.
1748  Otherwise, the parameter <tt>opt</tt> is optional unless the parameter
1749  <tt>sigma</tt> is omitted.
1750 
1751  If you pass \ref vigra::BlockwiseConvolutionOptions instead, the algorithm will
1752  be executed in parallel on data blocks of a certain size. The block size can be
1753  customized via <tt>BlockwiseConvolutionOptions::blockShape()</tt>, but the defaults
1754  usually work reasonably. By default, the number of threads equals the capabilities
1755  of your hardware, but you can change this via <tt>BlockwiseConvolutionOptions::numThreads()</tt>.
1756 
1757  <b> Declarations:</b>
1758 
1759  pass arbitrary-dimensional array views:
1760  \code
1761  namespace vigra {
1762  // execute algorithm sequentially
1763  template <unsigned int N, class T1, class S1,
1764  class T2, class S2>
1765  void
1766  symmetricGradientMultiArray(MultiArrayView<N, T1, S1> const & source,
1767  MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1768  ConvolutionOptions<N> opt = ConvolutionOptions<N>());
1769 
1770  // execute algorithm in parallel
1771  template <unsigned int N, class T1, class S1,
1772  class T2, class S2>
1773  void
1774  symmetricGradientMultiArray(MultiArrayView<N, T1, S1> const & source,
1775  MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1776  BlockwiseConvolutionOptions<N> opt);
1777  }
1778  \endcode
1779 
1780  \deprecatedAPI{symmetricGradientMultiArray}
1781  pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
1782  \code
1783  namespace vigra {
1784  template <class SrcIterator, class SrcShape, class SrcAccessor,
1785  class DestIterator, class DestAccessor>
1786  void
1787  symmetricGradientMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1788  DestIterator diter, DestAccessor dest,
1789  const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1790  }
1791  \endcode
1792  use argument objects in conjunction with \ref ArgumentObjectFactories :
1793  \code
1794  namespace vigra {
1795  template <class SrcIterator, class SrcShape, class SrcAccessor,
1796  class DestIterator, class DestAccessor>
1797  void
1798  symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1799  pair<DestIterator, DestAccessor> const & dest,
1800  const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1801  }
1802  \endcode
1803  \deprecatedEnd
1804 
1805  <b> Usage:</b>
1806 
1807  <b>\#include</b> <vigra/multi_convolution.hxx> (sequential version)<br/>
1808  <b>\#include</b> <vigra/multi_blockwise.hxx> (parallel version)<br/>
1809  Namespace: vigra
1810 
1811  \code
1812  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
1813  MultiArray<3, unsigned char> source(shape);
1814  MultiArray<3, TinyVector<float, 3> > dest(shape);
1815  ...
1816  // compute gradient
1817  symmetricGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest));
1818  \endcode
1819 
1820  <b> Usage with anisotropic data:</b>
1821 
1822  \code
1823  Shape3 shape(width, height, depth);
1824  MultiArray<3, unsigned char> source(shape);
1825  MultiArray<3, TinyVector<float, 3> > dest(shape);
1826  TinyVector<float, 3> step_size;
1827  ...
1828  // compute gradient
1829  symmetricGradientMultiArray(source, dest,
1830  ConvolutionOptions<3>().stepSize(step_size));
1831  \endcode
1832 
1833  \see convolveMultiArrayOneDimension()
1834 */
1836 
1837 template <class SrcIterator, class SrcShape, class SrcAccessor,
1838  class DestIterator, class DestAccessor>
1839 void
1840 symmetricGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1841  DestIterator di, DestAccessor dest,
1843 {
1844  typedef typename DestAccessor::value_type DestType;
1845  typedef typename DestType::value_type DestValueType;
1846  typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1847 
1848  static const int N = SrcShape::static_size;
1849  typedef typename ConvolutionOptions<N>::StepIterator StepType;
1850 
1851  for(int k=0; k<N; ++k)
1852  if(shape[k] <=0)
1853  return;
1854 
1855  vigra_precondition(N == (int)dest.size(di),
1856  "symmetricGradientMultiArray(): Wrong number of channels in output array.");
1857 
1858  Kernel1D<KernelType> filter;
1859  filter.initSymmetricDifference();
1860 
1861  StepType step_size_it = opt.stepParams();
1862 
1863  typedef VectorElementAccessor<DestAccessor> ElementAccessor;
1864 
1865  // compute gradient components
1866  for (int d = 0; d < N; ++d, ++step_size_it)
1867  {
1868  Kernel1D<KernelType> symmetric(filter);
1869  detail::scaleKernel(symmetric, 1 / *step_size_it);
1870  convolveMultiArrayOneDimension(si, shape, src,
1871  di, ElementAccessor(d, dest),
1872  d, symmetric, opt.from_point, opt.to_point);
1873  }
1874 }
1875 
1876 template <class SrcIterator, class SrcShape, class SrcAccessor,
1877  class DestIterator, class DestAccessor>
1878 inline void
1879 symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1880  pair<DestIterator, DestAccessor> const & dest,
1882 {
1883  symmetricGradientMultiArray(source.first, source.second, source.third,
1884  dest.first, dest.second, opt);
1885 }
1886 
1887 template <unsigned int N, class T1, class S1,
1888  class T2, class S2>
1889 inline void
1891  MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1893 {
1894  if(opt.to_point != typename MultiArrayShape<N>::type())
1895  {
1896  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
1897  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
1898  vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1899  "symmetricGradientMultiArray(): shape mismatch between ROI and output.");
1900  }
1901  else
1902  {
1903  vigra_precondition(source.shape() == dest.shape(),
1904  "symmetricGradientMultiArray(): shape mismatch between input and output.");
1905  }
1906 
1907  symmetricGradientMultiArray(srcMultiArrayRange(source),
1908  destMultiArray(dest), opt);
1909 }
1910 
1911 /********************************************************/
1912 /* */
1913 /* laplacianOfGaussianMultiArray */
1914 /* */
1915 /********************************************************/
1916 
1917 /** \brief Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters.
1918 
1919  This function computes the Laplacian of the given N-dimensional
1920  array with a sequence of second-derivative-of-Gaussian filters at the given
1921  standard deviation <tt>sigma</tt>. Both source and destination
1922  arrays must have scalar value_type. This function is implemented by calls to
1923  \ref separableConvolveMultiArray() with the appropriate kernels, followed by summation.
1924 
1925  Anisotropic data should be provided with appropriate \ref vigra::ConvolutionOptions
1926  to adjust the filter sizes for the resolution of each axis.
1927  Otherwise, the parameter <tt>opt</tt> is optional unless the parameter
1928  <tt>sigma</tt> is omitted.
1929 
1930  If you pass \ref vigra::BlockwiseConvolutionOptions instead, the algorithm will
1931  be executed in parallel on data blocks of a certain size. The block size can be
1932  customized via <tt>BlockwiseConvolutionOptions::blockShape()</tt>, but the defaults
1933  usually work reasonably. By default, the number of threads equals the capabilities
1934  of your hardware, but you can change this via <tt>BlockwiseConvolutionOptions::numThreads()</tt>.
1935 
1936  <b> Declarations:</b>
1937 
1938  pass arbitrary-dimensional array views:
1939  \code
1940  namespace vigra {
1941  // pass scale explicitly
1942  template <unsigned int N, class T1, class S1,
1943  class T2, class S2>
1944  void
1945  laplacianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
1946  MultiArrayView<N, T2, S2> dest,
1947  double sigma,
1948  ConvolutionOptions<N> opt = ConvolutionOptions<N>());
1949 
1950  // pass scale(s) in option object
1951  template <unsigned int N, class T1, class S1,
1952  class T2, class S2>
1953  void
1954  laplacianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
1955  MultiArrayView<N, T2, S2> dest,
1956  ConvolutionOptions<N> opt );
1957 
1958  // likewise, but execute algorithm in parallel
1959  template <unsigned int N, class T1, class S1,
1960  class T2, class S2>
1961  void
1962  laplacianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
1963  MultiArrayView<N, T2, S2> dest,
1964  BlockwiseConvolutionOptions<N> opt );
1965  }
1966  \endcode
1967 
1968  \deprecatedAPI{laplacianOfGaussianMultiArray}
1969  pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
1970  \code
1971  namespace vigra {
1972  template <class SrcIterator, class SrcShape, class SrcAccessor,
1973  class DestIterator, class DestAccessor>
1974  void
1975  laplacianOfGaussianMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1976  DestIterator diter, DestAccessor dest,
1977  double sigma,
1978  const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1979  }
1980  \endcode
1981  use argument objects in conjunction with \ref ArgumentObjectFactories :
1982  \code
1983  namespace vigra {
1984  template <class SrcIterator, class SrcShape, class SrcAccessor,
1985  class DestIterator, class DestAccessor>
1986  void
1987  laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1988  pair<DestIterator, DestAccessor> const & dest,
1989  double sigma,
1990  const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1991  }
1992  \endcode
1993  \deprecatedEnd
1994 
1995  <b> Usage:</b>
1996 
1997  <b>\#include</b> <vigra/multi_convolution.hxx> (sequential version)<br/>
1998  <b>\#include</b> <vigra/multi_blockwise.hxx> (parallel version)<br/>
1999  Namespace: vigra
2000 
2001  \code
2002  Shape3 shape(width, height, depth);
2003  MultiArray<3, float> source(shape);
2004  MultiArray<3, float> laplacian(shape);
2005  ...
2006  // compute Laplacian at scale sigma
2007  laplacianOfGaussianMultiArray(source, laplacian, sigma);
2008  \endcode
2009 
2010  <b> Usage with anisotropic data:</b>
2011 
2012  \code
2013  MultiArray<3, float> source(shape);
2014  MultiArray<3, float> laplacian(shape);
2015  TinyVector<float, 3> step_size;
2016  TinyVector<float, 3> resolution_sigmas;
2017  ...
2018  // compute Laplacian at scale sigma
2019  laplacianOfGaussianMultiArray(source, laplacian, sigma,
2020  ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
2021  \endcode
2022 
2023  \see separableConvolveMultiArray()
2024 */
2026 
2027 template <class SrcIterator, class SrcShape, class SrcAccessor,
2028  class DestIterator, class DestAccessor>
2029 void
2030 laplacianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
2031  DestIterator di, DestAccessor dest,
2033 {
2034  using namespace functor;
2035 
2036  typedef typename DestAccessor::value_type DestType;
2037  typedef typename NumericTraits<DestType>::RealPromote KernelType;
2038  typedef typename AccessorTraits<KernelType>::default_accessor DerivativeAccessor;
2039 
2040  static const int N = SrcShape::static_size;
2041  typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
2042 
2043  ParamType params = opt.scaleParams();
2044  ParamType params2(params);
2045 
2046  ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
2047  for (int dim = 0; dim < N; ++dim, ++params)
2048  {
2049  double sigma = params.sigma_scaled("laplacianOfGaussianMultiArray");
2050  plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
2051  }
2052 
2053  SrcShape dshape(shape);
2054  if(opt.to_point != SrcShape())
2055  dshape = opt.to_point - opt.from_point;
2056 
2057  MultiArray<N, KernelType> derivative(dshape);
2058 
2059  // compute 2nd derivatives and sum them up
2060  for (int dim = 0; dim < N; ++dim, ++params2)
2061  {
2062  ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
2063  kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 2, 1.0, opt.window_ratio);
2064  detail::scaleKernel(kernels[dim], 1.0 / sq(params2.step_size()));
2065 
2066  if (dim == 0)
2067  {
2068  separableConvolveMultiArray( si, shape, src,
2069  di, dest, kernels.begin(), opt.from_point, opt.to_point);
2070  }
2071  else
2072  {
2073  separableConvolveMultiArray( si, shape, src,
2074  derivative.traverser_begin(), DerivativeAccessor(),
2075  kernels.begin(), opt.from_point, opt.to_point);
2076  combineTwoMultiArrays(di, dshape, dest, derivative.traverser_begin(), DerivativeAccessor(),
2077  di, dest, Arg1() + Arg2() );
2078  }
2079  }
2080 }
2081 
2082 template <class SrcIterator, class SrcShape, class SrcAccessor,
2083  class DestIterator, class DestAccessor>
2084 void
2085 laplacianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
2086  DestIterator di, DestAccessor dest, double sigma,
2088 {
2089  laplacianOfGaussianMultiArray(si, shape, src, di, dest, opt.stdDev(sigma));
2090 }
2091 
2092 template <class SrcIterator, class SrcShape, class SrcAccessor,
2093  class DestIterator, class DestAccessor>
2094 inline void
2095 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2096  pair<DestIterator, DestAccessor> const & dest,
2098 {
2099  laplacianOfGaussianMultiArray( source.first, source.second, source.third,
2100  dest.first, dest.second, opt );
2101 }
2102 
2103 template <class SrcIterator, class SrcShape, class SrcAccessor,
2104  class DestIterator, class DestAccessor>
2105 inline void
2106 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2107  pair<DestIterator, DestAccessor> const & dest,
2108  double sigma,
2110 {
2111  laplacianOfGaussianMultiArray( source.first, source.second, source.third,
2112  dest.first, dest.second, sigma, opt );
2113 }
2114 
2115 template <unsigned int N, class T1, class S1,
2116  class T2, class S2>
2117 inline void
2120  ConvolutionOptions<N> opt )
2121 {
2122  if(opt.to_point != typename MultiArrayShape<N>::type())
2123  {
2124  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
2125  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
2126  vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2127  "laplacianOfGaussianMultiArray(): shape mismatch between ROI and output.");
2128  }
2129  else
2130  {
2131  vigra_precondition(source.shape() == dest.shape(),
2132  "laplacianOfGaussianMultiArray(): shape mismatch between input and output.");
2133  }
2134 
2135  laplacianOfGaussianMultiArray( srcMultiArrayRange(source),
2136  destMultiArray(dest), opt );
2137 }
2138 
2139 template <unsigned int N, class T1, class S1,
2140  class T2, class S2>
2141 inline void
2144  double sigma,
2146 {
2147  laplacianOfGaussianMultiArray( source, dest, opt.stdDev(sigma) );
2148 }
2149 
2150 /********************************************************/
2151 /* */
2152 /* gaussianDivergenceMultiArray */
2153 /* */
2154 /********************************************************/
2155 
2156 /** \brief Calculate the divergence of a vector field using Gaussian derivative filters.
2157 
2158  This function computes the divergence of the given N-dimensional vector field
2159  with a sequence of first-derivative-of-Gaussian filters at the given
2160  standard deviation <tt>sigma</tt>. The input vector field can either be given as a sequence
2161  of scalar array views (one for each vector field component), represented by an
2162  iterator range, or by a single vector array with the appropriate shape.
2163  This function is implemented by calls to
2164  \ref separableConvolveMultiArray() with the suitable kernels, followed by summation.
2165 
2166  Anisotropic data should be provided with appropriate \ref vigra::ConvolutionOptions
2167  to adjust the filter sizes for the resolution of each axis.
2168  Otherwise, the parameter <tt>opt</tt> is optional unless the parameter
2169  <tt>sigma</tt> is omitted.
2170 
2171  If you pass \ref vigra::BlockwiseConvolutionOptions instead, the algorithm will
2172  be executed in parallel on data blocks of a certain size. The block size can be
2173  customized via <tt>BlockwiseConvolutionOptions::blockShape()</tt>, but the defaults
2174  usually work reasonably. By default, the number of threads equals the capabilities
2175  of your hardware, but you can change this via <tt>BlockwiseConvolutionOptions::numThreads()</tt>.
2176 
2177  <b> Declarations:</b>
2178 
2179  pass arbitrary-dimensional array views:
2180  \code
2181  namespace vigra {
2182  // specify input vector field as a sequence of scalar arrays
2183  template <class Iterator,
2184  unsigned int N, class T, class S>
2185  void
2186  gaussianDivergenceMultiArray(Iterator vectorField, Iterator vectorFieldEnd,
2187  MultiArrayView<N, T, S> divergence,
2188  ConvolutionOptions<N> const & opt);
2189 
2190  template <class Iterator,
2191  unsigned int N, class T, class S>
2192  void
2193  gaussianDivergenceMultiArray(Iterator vectorField, Iterator vectorFieldEnd,
2194  MultiArrayView<N, T, S> divergence,
2195  double sigma,
2196  ConvolutionOptions<N> opt = ConvolutionOptions<N>());
2197 
2198  // pass input vector field as an array of vectors
2199  template <unsigned int N, class T1, class S1,
2200  class T2, class S2>
2201  void
2202  gaussianDivergenceMultiArray(MultiArrayView<N, TinyVector<T1, N>, S1> const & vectorField,
2203  MultiArrayView<N, T2, S2> divergence,
2204  ConvolutionOptions<N> const & opt);
2205 
2206  template <unsigned int N, class T1, class S1,
2207  class T2, class S2>
2208  void
2209  gaussianDivergenceMultiArray(MultiArrayView<N, TinyVector<T1, N>, S1> const & vectorField,
2210  MultiArrayView<N, T2, S2> divergence,
2211  double sigma,
2212  ConvolutionOptions<N> opt = ConvolutionOptions<N>());
2213 
2214  // pass input vector field as an array of vectors and
2215  // execute algorithm in parallel
2216  template <unsigned int N, class T1, class S1,
2217  class T2, class S2>
2218  void
2219  gaussianDivergenceMultiArray(MultiArrayView<N, TinyVector<T1, N>, S1> const & vectorField,
2220  MultiArrayView<N, T2, S2> divergence,
2221  BlockwiseConvolutionOptions<N> const & opt);
2222  }
2223  \endcode
2224 
2225  <b> Usage:</b>
2226 
2227  <b>\#include</b> <vigra/multi_convolution.hxx> (sequential version)<br/>
2228  <b>\#include</b> <vigra/multi_blockwise.hxx> (parallel version)<br/>
2229  Namespace: vigra
2230 
2231  \code
2232  Shape3 shape(width, height, depth);
2233  MultiArray<3, TinyVector<float, 3> > source(shape);
2234  MultiArray<3, float> laplacian(shape);
2235  ...
2236  // compute divergence at scale sigma
2237  gaussianDivergenceMultiArray(source, laplacian, sigma);
2238  \endcode
2239 
2240  <b> Usage with anisotropic data:</b>
2241 
2242  \code
2243  MultiArray<3, TinyVector<float, 3> > source(shape);
2244  MultiArray<3, float> laplacian(shape);
2245  TinyVector<float, 3> step_size;
2246  TinyVector<float, 3> resolution_sigmas;
2247  ...
2248  // compute divergence at scale sigma
2249  gaussianDivergenceMultiArray(source, laplacian, sigma,
2250  ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
2251  \endcode
2252 */
2254 
2255 template <class Iterator,
2256  unsigned int N, class T, class S>
2257 void
2258 gaussianDivergenceMultiArray(Iterator vectorField, Iterator vectorFieldEnd,
2259  MultiArrayView<N, T, S> divergence,
2261 {
2262  typedef typename std::iterator_traits<Iterator>::value_type ArrayType;
2263  typedef typename ArrayType::value_type SrcType;
2264  typedef typename NumericTraits<SrcType>::RealPromote TmpType;
2265  typedef Kernel1D<double> Kernel;
2266 
2267  vigra_precondition(std::distance(vectorField, vectorFieldEnd) == N,
2268  "gaussianDivergenceMultiArray(): wrong number of input arrays.");
2269  // more checks are performed in separableConvolveMultiArray()
2270 
2271  typename ConvolutionOptions<N>::ScaleIterator params = opt.scaleParams();
2272  ArrayVector<double> sigmas(N);
2273  ArrayVector<Kernel> kernels(N);
2274  for(unsigned int k = 0; k < N; ++k, ++params)
2275  {
2276  sigmas[k] = params.sigma_scaled("gaussianDivergenceMultiArray");
2277  kernels[k].initGaussian(sigmas[k], 1.0, opt.window_ratio);
2278  }
2279 
2280  MultiArray<N, TmpType> tmpDeriv(divergence.shape());
2281 
2282  for(unsigned int k=0; k < N; ++k, ++vectorField)
2283  {
2284  kernels[k].initGaussianDerivative(sigmas[k], 1, 1.0, opt.window_ratio);
2285  if(k == 0)
2286  {
2287  separableConvolveMultiArray(*vectorField, divergence, kernels.begin(), opt.from_point, opt.to_point);
2288  }
2289  else
2290  {
2291  separableConvolveMultiArray(*vectorField, tmpDeriv, kernels.begin(), opt.from_point, opt.to_point);
2292  divergence += tmpDeriv;
2293  }
2294  kernels[k].initGaussian(sigmas[k], 1.0, opt.window_ratio);
2295  }
2296 }
2297 
2298 template <class Iterator,
2299  unsigned int N, class T, class S>
2300 inline void
2301 gaussianDivergenceMultiArray(Iterator vectorField, Iterator vectorFieldEnd,
2302  MultiArrayView<N, T, S> divergence,
2303  double sigma,
2305 {
2306  gaussianDivergenceMultiArray(vectorField, vectorFieldEnd, divergence, opt.stdDev(sigma));
2307 }
2308 
2309 template <unsigned int N, class T1, class S1,
2310  class T2, class S2>
2311 inline void
2313  MultiArrayView<N, T2, S2> divergence,
2314  ConvolutionOptions<N> const & opt)
2315 {
2317  for(unsigned int k=0; k<N; ++k)
2318  field.push_back(vectorField.bindElementChannel(k));
2319 
2320  gaussianDivergenceMultiArray(field.begin(), field.end(), divergence, opt);
2321 }
2322 
2323 template <unsigned int N, class T1, class S1,
2324  class T2, class S2>
2325 inline void
2327  MultiArrayView<N, T2, S2> divergence,
2328  double sigma,
2330 {
2331  gaussianDivergenceMultiArray(vectorField, divergence, opt.stdDev(sigma));
2332 }
2333 
2334 /********************************************************/
2335 /* */
2336 /* hessianOfGaussianMultiArray */
2337 /* */
2338 /********************************************************/
2339 
2340 /** \brief Calculate Hessian matrix of a N-dimensional arrays using Gaussian derivative filters.
2341 
2342  This function computes the Hessian matrix the given scalar N-dimensional
2343  array with a sequence of second-derivative-of-Gaussian filters at the given
2344  standard deviation <tt>sigma</tt>. The destination array must
2345  have a vector valued element type with N*(N+1)/2 elements (it represents the
2346  upper triangular part of the symmetric Hessian matrix, flattened row-wise).
2347  This function is implemented by calls to
2348  \ref separableConvolveMultiArray() with the appropriate kernels.
2349 
2350  Anisotropic data should be provided with appropriate \ref vigra::ConvolutionOptions
2351  to adjust the filter sizes for the resolution of each axis.
2352  Otherwise, the parameter <tt>opt</tt> is optional unless the parameter
2353  <tt>sigma</tt> is omitted.
2354 
2355  If you pass \ref vigra::BlockwiseConvolutionOptions instead, the algorithm will
2356  be executed in parallel on data blocks of a certain size. The block size can be
2357  customized via <tt>BlockwiseConvolutionOptions::blockShape()</tt>, but the defaults
2358  usually work reasonably. By default, the number of threads equals the capabilities
2359  of your hardware, but you can change this via <tt>BlockwiseConvolutionOptions::numThreads()</tt>.
2360 
2361  <b> Declarations:</b>
2362 
2363  pass arbitrary-dimensional array views:
2364  \code
2365  namespace vigra {
2366  // pass scale explicitly
2367  template <unsigned int N, class T1, class S1,
2368  class T2, class S2>
2369  void
2370  hessianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
2371  MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2372  double sigma,
2373  ConvolutionOptions<N> opt = ConvolutionOptions<N>());
2374 
2375  // pass scale(s) in option object
2376  template <unsigned int N, class T1, class S1,
2377  class T2, class S2>
2378  void
2379  hessianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
2380  MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2381  ConvolutionOptions<N> opt);
2382 
2383  // likewise, but execute algorithm in parallel
2384  template <unsigned int N, class T1, class S1,
2385  class T2, class S2>
2386  void
2387  hessianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
2388  MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2389  BlockwiseConvolutionOptions<N> opt);
2390  }
2391  \endcode
2392 
2393  \deprecatedAPI{hessianOfGaussianMultiArray}
2394  pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
2395  \code
2396  namespace vigra {
2397  template <class SrcIterator, class SrcShape, class SrcAccessor,
2398  class DestIterator, class DestAccessor>
2399  void
2400  hessianOfGaussianMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
2401  DestIterator diter, DestAccessor dest,
2402  double sigma,
2403  const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
2404  }
2405  \endcode
2406  use argument objects in conjunction with \ref ArgumentObjectFactories :
2407  \code
2408  namespace vigra {
2409  template <class SrcIterator, class SrcShape, class SrcAccessor,
2410  class DestIterator, class DestAccessor>
2411  void
2412  hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2413  pair<DestIterator, DestAccessor> const & dest,
2414  double sigma,
2415  const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
2416  }
2417  \endcode
2418  \deprecatedEnd
2419 
2420  <b> Usage:</b>
2421 
2422  <b>\#include</b> <vigra/multi_convolution.hxx> (sequential version)<br/>
2423  <b>\#include</b> <vigra/multi_blockwise.hxx> (parallel version)<br/>
2424  Namespace: vigra
2425 
2426  \code
2427  Shape3 shape(width, height, depth);
2428  MultiArray<3, float> source(shape);
2429  MultiArray<3, TinyVector<float, 6> > dest(shape);
2430  ...
2431  // compute Hessian at scale sigma
2432  hessianOfGaussianMultiArray(source, dest, sigma);
2433  \endcode
2434 
2435  <b> Usage with anisotropic data:</b>
2436 
2437  \code
2438  MultiArray<3, float> source(shape);
2439  MultiArray<3, TinyVector<float, 6> > dest(shape);
2440  TinyVector<float, 3> step_size;
2441  TinyVector<float, 3> resolution_sigmas;
2442  ...
2443  // compute Hessian at scale sigma
2444  hessianOfGaussianMultiArray(source, dest, sigma,
2445  ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
2446  \endcode
2447 
2448  \see separableConvolveMultiArray(), vectorToTensorMultiArray()
2449 */
2451 
2452 template <class SrcIterator, class SrcShape, class SrcAccessor,
2453  class DestIterator, class DestAccessor>
2454 void
2455 hessianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
2456  DestIterator di, DestAccessor dest,
2458 {
2459  typedef typename DestAccessor::value_type DestType;
2460  typedef typename DestType::value_type DestValueType;
2461  typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
2462 
2463  static const int N = SrcShape::static_size;
2464  static const int M = N*(N+1)/2;
2465  typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
2466 
2467  for(int k=0; k<N; ++k)
2468  if(shape[k] <=0)
2469  return;
2470 
2471  vigra_precondition(M == (int)dest.size(di),
2472  "hessianOfGaussianMultiArray(): Wrong number of channels in output array.");
2473 
2474  ParamType params_init = opt.scaleParams();
2475 
2476  ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
2477  ParamType params(params_init);
2478  for (int dim = 0; dim < N; ++dim, ++params)
2479  {
2480  double sigma = params.sigma_scaled("hessianOfGaussianMultiArray");
2481  plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
2482  }
2483 
2484  typedef VectorElementAccessor<DestAccessor> ElementAccessor;
2485 
2486  // compute elements of the Hessian matrix
2487  ParamType params_i(params_init);
2488  for (int b=0, i=0; i<N; ++i, ++params_i)
2489  {
2490  ParamType params_j(params_i);
2491  for (int j=i; j<N; ++j, ++b, ++params_j)
2492  {
2493  ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
2494  if(i == j)
2495  {
2496  kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 2, 1.0, opt.window_ratio);
2497  }
2498  else
2499  {
2500  kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 1, 1.0, opt.window_ratio);
2501  kernels[j].initGaussianDerivative(params_j.sigma_scaled(), 1, 1.0, opt.window_ratio);
2502  }
2503  detail::scaleKernel(kernels[i], 1 / params_i.step_size());
2504  detail::scaleKernel(kernels[j], 1 / params_j.step_size());
2505  separableConvolveMultiArray(si, shape, src, di, ElementAccessor(b, dest),
2506  kernels.begin(), opt.from_point, opt.to_point);
2507  }
2508  }
2509 }
2510 
2511 template <class SrcIterator, class SrcShape, class SrcAccessor,
2512  class DestIterator, class DestAccessor>
2513 inline void
2514 hessianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
2515  DestIterator di, DestAccessor dest, double sigma,
2517 {
2518  hessianOfGaussianMultiArray(si, shape, src, di, dest, opt.stdDev(sigma));
2519 }
2520 
2521 template <class SrcIterator, class SrcShape, class SrcAccessor,
2522  class DestIterator, class DestAccessor>
2523 inline void
2524 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2525  pair<DestIterator, DestAccessor> const & dest,
2527 {
2528  hessianOfGaussianMultiArray( source.first, source.second, source.third,
2529  dest.first, dest.second, opt );
2530 }
2531 
2532 template <class SrcIterator, class SrcShape, class SrcAccessor,
2533  class DestIterator, class DestAccessor>
2534 inline void
2535 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2536  pair<DestIterator, DestAccessor> const & dest,
2537  double sigma,
2539 {
2540  hessianOfGaussianMultiArray( source.first, source.second, source.third,
2541  dest.first, dest.second, sigma, opt );
2542 }
2543 
2544 template <unsigned int N, class T1, class S1,
2545  class T2, class S2>
2546 inline void
2548  MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2549  ConvolutionOptions<N> opt )
2550 {
2551  if(opt.to_point != typename MultiArrayShape<N>::type())
2552  {
2553  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
2554  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
2555  vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2556  "hessianOfGaussianMultiArray(): shape mismatch between ROI and output.");
2557  }
2558  else
2559  {
2560  vigra_precondition(source.shape() == dest.shape(),
2561  "hessianOfGaussianMultiArray(): shape mismatch between input and output.");
2562  }
2563 
2564  hessianOfGaussianMultiArray( srcMultiArrayRange(source),
2565  destMultiArray(dest), opt );
2566 }
2567 
2568 template <unsigned int N, class T1, class S1,
2569  class T2, class S2>
2570 inline void
2572  MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2573  double sigma,
2575 {
2576  hessianOfGaussianMultiArray( source, dest, opt.stdDev(sigma) );
2577 }
2578 
2579 namespace detail {
2580 
2581 template<int N, class VectorType>
2582 struct StructurTensorFunctor
2583 {
2584  typedef VectorType result_type;
2585  typedef typename VectorType::value_type ValueType;
2586 
2587  template <class T>
2588  VectorType operator()(T const & in) const
2589  {
2590  VectorType res;
2591  for(int b=0, i=0; i<N; ++i)
2592  {
2593  for(int j=i; j<N; ++j, ++b)
2594  {
2595  res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
2596  }
2597  }
2598  return res;
2599  }
2600 };
2601 
2602 } // namespace detail
2603 
2604 /********************************************************/
2605 /* */
2606 /* structureTensorMultiArray */
2607 /* */
2608 /********************************************************/
2609 
2610 /** \brief Calculate th structure tensor of a multi-dimensional arrays.
2611 
2612  This function computes the gradient (outer product) tensor for each element
2613  of the given N-dimensional array with first-derivative-of-Gaussian filters at
2614  the given <tt>innerScale</tt>, followed by Gaussian smoothing at <tt>outerScale</tt>.
2615  The destination array must have a vector valued pixel type with
2616  N*(N+1)/2 elements (it represents the upper triangular part of the symmetric
2617  structure tensor matrix, flattened row-wise). If the source array is also vector valued, the
2618  resulting structure tensor is the sum of the individual tensors for each channel.
2619  This function is implemented by calls to
2620  \ref separableConvolveMultiArray() with the appropriate kernels.
2621 
2622  Anisotropic data should be provided with appropriate \ref vigra::ConvolutionOptions
2623  to adjust the filter sizes for the resolution of each axis.
2624  Otherwise, the parameter <tt>opt</tt> is optional unless the parameters
2625  <tt>innerScale</tt> and <tt>outerScale</tt> are both omitted.
2626 
2627  If you pass \ref vigra::BlockwiseConvolutionOptions instead, the algorithm will
2628  be executed in parallel on data blocks of a certain size. The block size can be
2629  customized via <tt>BlockwiseConvolutionOptions::blockShape()</tt>, but the defaults
2630  usually work reasonably. By default, the number of threads equals the capabilities
2631  of your hardware, but you can change this via <tt>BlockwiseConvolutionOptions::numThreads()</tt>.
2632 
2633  <b> Declarations:</b>
2634 
2635  pass arbitrary-dimensional array views:
2636  \code
2637  namespace vigra {
2638  // pass scales explicitly
2639  template <unsigned int N, class T1, class S1,
2640  class T2, class S2>
2641  void
2642  structureTensorMultiArray(MultiArrayView<N, T1, S1> const & source,
2643  MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2644  double innerScale, double outerScale,
2645  ConvolutionOptions<N> opt = ConvolutionOptions<N>());
2646 
2647  // pass scales in option object
2648  template <unsigned int N, class T1, class S1,
2649  class T2, class S2>
2650  void
2651  structureTensorMultiArray(MultiArrayView<N, T1, S1> const & source,
2652  MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2653  ConvolutionOptions<N> opt );
2654 
2655  // likewise, but execute algorithm in parallel
2656  template <unsigned int N, class T1, class S1,
2657  class T2, class S2>
2658  void
2659  structureTensorMultiArray(MultiArrayView<N, T1, S1> const & source,
2660  MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2661  BlockwiseConvolutionOptions<N> opt );
2662  }
2663  \endcode
2664 
2665  \deprecatedAPI{structureTensorMultiArray}
2666  pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
2667  \code
2668  namespace vigra {
2669  template <class SrcIterator, class SrcShape, class SrcAccessor,
2670  class DestIterator, class DestAccessor>
2671  void
2672  structureTensorMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
2673  DestIterator diter, DestAccessor dest,
2674  double innerScale, double outerScale,
2675  ConvolutionOptions<N> opt);
2676  }
2677  \endcode
2678  use argument objects in conjunction with \ref ArgumentObjectFactories :
2679  \code
2680  namespace vigra {
2681  template <class SrcIterator, class SrcShape, class SrcAccessor,
2682  class DestIterator, class DestAccessor>
2683  void
2684  structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2685  pair<DestIterator, DestAccessor> const & dest,
2686  double innerScale, double outerScale,
2687  const ConvolutionOptions<N> & opt);
2688  }
2689  \endcode
2690  \deprecatedEnd
2691 
2692  <b> Usage:</b>
2693 
2694  <b>\#include</b> <vigra/multi_convolution.hxx> (sequential version)<br/>
2695  <b>\#include</b> <vigra/multi_blockwise.hxx> (parallel version)<br/>
2696  Namespace: vigra
2697 
2698  \code
2699  Shape3 shape(width, height, depth);
2700  MultiArray<3, RGBValue<float> > source(shape);
2701  MultiArray<3, TinyVector<float, 6> > dest(shape);
2702  ...
2703  // compute structure tensor at scales innerScale and outerScale
2704  structureTensorMultiArray(source, dest, innerScale, outerScale);
2705  \endcode
2706 
2707  <b> Usage with anisotropic data:</b>
2708 
2709  \code
2710  MultiArray<3, RGBValue<float> > source(shape);
2711  MultiArray<3, TinyVector<float, 6> > dest(shape);
2712  TinyVector<float, 3> step_size;
2713  TinyVector<float, 3> resolution_sigmas;
2714  ...
2715  // compute structure tensor at scales innerScale and outerScale
2716  structureTensorMultiArray(source, dest, innerScale, outerScale,
2717  ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
2718  \endcode
2719 
2720  \see separableConvolveMultiArray(), vectorToTensorMultiArray()
2721 */
2723 
2724 template <class SrcIterator, class SrcShape, class SrcAccessor,
2725  class DestIterator, class DestAccessor>
2726 void
2727 structureTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
2728  DestIterator di, DestAccessor dest,
2730 {
2731  static const int N = SrcShape::static_size;
2732  static const int M = N*(N+1)/2;
2733 
2734  typedef typename DestAccessor::value_type DestType;
2735  typedef typename DestType::value_type DestValueType;
2736  typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
2737  typedef TinyVector<KernelType, N> GradientVector;
2738  typedef typename AccessorTraits<GradientVector>::default_accessor GradientAccessor;
2739  typedef typename AccessorTraits<DestType>::default_accessor GradientTensorAccessor;
2740 
2741  for(int k=0; k<N; ++k)
2742  if(shape[k] <=0)
2743  return;
2744 
2745  vigra_precondition(M == (int)dest.size(di),
2746  "structureTensorMultiArray(): Wrong number of channels in output array.");
2747 
2748  ConvolutionOptions<N> innerOptions = opt;
2749  ConvolutionOptions<N> outerOptions = opt.outerOptions();
2750  typename ConvolutionOptions<N>::ScaleIterator params = outerOptions.scaleParams();
2751 
2752  SrcShape gradientShape(shape);
2753  if(opt.to_point != SrcShape())
2754  {
2755  detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.from_point);
2756  detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.to_point);
2757 
2758  for(int k=0; k<N; ++k, ++params)
2759  {
2760  Kernel1D<double> gauss;
2761  gauss.initGaussian(params.sigma_scaled("structureTensorMultiArray"), 1.0, opt.window_ratio);
2762  int dilation = gauss.right();
2763  innerOptions.from_point[k] = std::max<MultiArrayIndex>(0, opt.from_point[k] - dilation);
2764  innerOptions.to_point[k] = std::min<MultiArrayIndex>(shape[k], opt.to_point[k] + dilation);
2765  }
2766  outerOptions.from_point -= innerOptions.from_point;
2767  outerOptions.to_point -= innerOptions.from_point;
2768  gradientShape = innerOptions.to_point - innerOptions.from_point;
2769  }
2770 
2771  MultiArray<N, GradientVector> gradient(gradientShape);
2772  MultiArray<N, DestType> gradientTensor(gradientShape);
2773  gaussianGradientMultiArray(si, shape, src,
2774  gradient.traverser_begin(), GradientAccessor(),
2775  innerOptions,
2776  "structureTensorMultiArray");
2777 
2778  transformMultiArray(gradient.traverser_begin(), gradientShape, GradientAccessor(),
2779  gradientTensor.traverser_begin(), GradientTensorAccessor(),
2780  detail::StructurTensorFunctor<N, DestType>());
2781 
2782  gaussianSmoothMultiArray(gradientTensor.traverser_begin(), gradientShape, GradientTensorAccessor(),
2783  di, dest, outerOptions,
2784  "structureTensorMultiArray");
2785 }
2786 
2787 template <class SrcIterator, class SrcShape, class SrcAccessor,
2788  class DestIterator, class DestAccessor>
2789 inline void
2790 structureTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
2791  DestIterator di, DestAccessor dest,
2792  double innerScale, double outerScale,
2794 {
2795  structureTensorMultiArray(si, shape, src, di, dest,
2796  opt.stdDev(innerScale).outerScale(outerScale));
2797 }
2798 
2799 template <class SrcIterator, class SrcShape, class SrcAccessor,
2800  class DestIterator, class DestAccessor>
2801 inline void
2802 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2803  pair<DestIterator, DestAccessor> const & dest,
2805 {
2806  structureTensorMultiArray( source.first, source.second, source.third,
2807  dest.first, dest.second, opt );
2808 }
2809 
2810 
2811 template <class SrcIterator, class SrcShape, class SrcAccessor,
2812  class DestIterator, class DestAccessor>
2813 inline void
2814 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2815  pair<DestIterator, DestAccessor> const & dest,
2816  double innerScale, double outerScale,
2818 {
2819  structureTensorMultiArray( source.first, source.second, source.third,
2820  dest.first, dest.second,
2821  innerScale, outerScale, opt);
2822 }
2823 
2824 template <unsigned int N, class T1, class S1,
2825  class T2, class S2>
2826 inline void
2828  MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2829  ConvolutionOptions<N> opt )
2830 {
2831  if(opt.to_point != typename MultiArrayShape<N>::type())
2832  {
2833  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
2834  detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
2835  vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2836  "structureTensorMultiArray(): shape mismatch between ROI and output.");
2837  }
2838  else
2839  {
2840  vigra_precondition(source.shape() == dest.shape(),
2841  "structureTensorMultiArray(): shape mismatch between input and output.");
2842  }
2843 
2844  structureTensorMultiArray( srcMultiArrayRange(source),
2845  destMultiArray(dest), opt );
2846 }
2847 
2848 
2849 template <unsigned int N, class T1, class S1,
2850  class T2, class S2>
2851 inline void
2853  MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2854  double innerScale, double outerScale,
2856 {
2857  structureTensorMultiArray(source, dest, opt.innerScale(innerScale).outerScale(outerScale));
2858 }
2859 
2860 //@}
2861 
2862 } //-- namespace vigra
2863 
2864 
2865 #endif //-- VIGRA_MULTI_CONVOLUTION_H
Accessor for one component of a vector.
Definition: accessor.hxx:539
Generic 1 dimensional convolution kernel.
Definition: separableconvolution.hxx:52
void indexSort(Iterator first, Iterator last, IndexIterator index_first, Compare c)
Return the index permutation that would sort the input array.
Definition: algorithm.hxx:414
int right() const
Definition: separableconvolution.hxx:2157
void initGaussian(double std_dev, value_type norm, double windowRatio=0.0)
Definition: separableconvolution.hxx:2253
MultiArrayView< N-M, T, StrideTag > bindOuter(const TinyVector< Index, M > &d) const
Definition: multi_array.hxx:2132
void gaussianSmoothMultiArray(...)
Isotropic Gaussian smoothing of a multi-dimensional arrays.
void gaussianGradientMagnitude(...)
Calculate the gradient magnitude by means of a 1st derivatives of Gaussian filter.
void convolveMultiArrayOneDimension(...)
Convolution along a single dimension of a multi-dimensional arrays.
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
void laplacianOfGaussianMultiArray(...)
Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters.
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2422
iterator end()
Definition: tinyvector.hxx:864
ConvolutionOptions< dim > & subarray(Shape const &from, Shape const &to)
Definition: multi_convolution.hxx:498
void convolveLine(...)
Performs a 1-dimensional convolution of the source signal using the given kernel. ...
Definition: accessor.hxx:43
void gaussianDivergenceMultiArray(...)
Calculate the divergence of a vector field using Gaussian derivative filters.
MultiArrayView< N+1, T, StrideTag > insertSingletonDimension(difference_type_1 i) const
Definition: multi_array.hxx:2302
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
void gaussianGradientMultiArray(...)
Calculate Gaussian gradient of a multi-dimensional arrays.
Definition: array_vector.hxx:58
Definition: multi_fwd.hxx:63
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:365
A navigator that provides access to the 1D subranges of an n-dimensional range given by a vigra::Mult...
Definition: navigator.hxx:97
const difference_type & shape() const
Definition: multi_array.hxx:1596
iterator begin()
Definition: tinyvector.hxx:861
void symmetricGradientMultiArray(...)
Calculate gradient of a multi-dimensional arrays using symmetric difference filters.
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
ConvolutionOptions< dim > & filterWindowSize(double ratio)
Definition: multi_convolution.hxx:476
Options class template for convolutions.
Definition: multi_convolution.hxx:335
void copyMultiArray(...)
Copy a multi-dimensional array.
void outer(const MultiArrayView< 2, T, C1 > &x, const MultiArrayView< 2, T, C2 > &y, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1457
void hessianOfGaussianMultiArray(...)
Calculate Hessian matrix of a N-dimensional arrays using Gaussian derivative filters.
void combineTwoMultiArrays(...)
Combine two multi-dimensional arrays into one using a binary function or functor. ...
const_iterator begin() const
Definition: array_vector.hxx:223
Encapsulate read access to the values an iterator points to.
Definition: accessor.hxx:269
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:652
void transformMultiArray(...)
Transform a multi-dimensional array with a unary function or functor.
ConvolutionOptions< dim > & innerScale(...)
MultiArrayView & init(const U &init)
Definition: multi_array.hxx:1154
Class for a single RGB value.
Definition: accessor.hxx:938
Encapsulate access to the values an iterator points to.
Definition: accessor.hxx:133
const_iterator end() const
Definition: array_vector.hxx:237
void structureTensorMultiArray(...)
Calculate th structure tensor of a multi-dimensional arrays.
MultiArrayView< N+1, typename ExpandElementResult< T >::type, StridedArrayTag > expandElements(difference_type_1 d) const
Definition: multi_array.hxx:2273
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
ConvolutionOptions< dim > & stdDev(...)
void separableConvolveMultiArray(...)
Separated convolution on 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