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