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

© 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 (Mon Nov 18 2013)