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

multi_convolution.hxx
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 "functorexpression.hxx"
50 #include "tinyvector.hxx"
51 #include "algorithm.hxx"
52 
53 namespace vigra
54 {
55 
56 namespace detail
57 {
58 
59 struct DoubleYielder
60 {
61  const double value;
62  DoubleYielder(double v, unsigned, const char *const) : value(v) {}
63  DoubleYielder(double v) : value(v) {}
64  void operator++() {}
65  double operator*() const { return value; }
66 };
67 
68 template <typename X>
69 struct IteratorDoubleYielder
70 {
71  X it;
72  IteratorDoubleYielder(X i, unsigned, const char *const) : it(i) {}
73  IteratorDoubleYielder(X i) : it(i) {}
74  void operator++() { ++it; }
75  double operator*() const { return *it; }
76 };
77 
78 template <typename X>
79 struct SequenceDoubleYielder
80 {
81  typename X::const_iterator it;
82  SequenceDoubleYielder(const X & seq, unsigned dim,
83  const char *const function_name = "SequenceDoubleYielder")
84  : it(seq.begin())
85  {
86  if (seq.size() == dim)
87  return;
88  std::string msg = "(): Parameter number be equal to the number of spatial dimensions.";
89  vigra_precondition(false, function_name + msg);
90  }
91  void operator++() { ++it; }
92  double operator*() const { return *it; }
93 };
94 
95 template <typename X>
96 struct WrapDoubleIterator
97 {
98  typedef
99  typename IfBool< IsConvertibleTo<X, double>::value,
100  DoubleYielder,
101  typename IfBool< IsIterator<X>::value || IsArray<X>::value,
102  IteratorDoubleYielder<X>,
103  SequenceDoubleYielder<X>
104  >::type
105  >::type type;
106 };
107 
108 template <class Param1, class Param2, class Param3>
109 struct WrapDoubleIteratorTriple
110 {
111  typename WrapDoubleIterator<Param1>::type sigma_eff_it;
112  typename WrapDoubleIterator<Param2>::type sigma_d_it;
113  typename WrapDoubleIterator<Param3>::type step_size_it;
114  WrapDoubleIteratorTriple(Param1 sigma_eff, Param2 sigma_d, Param3 step_size)
115  : sigma_eff_it(sigma_eff), sigma_d_it(sigma_d), step_size_it(step_size) {}
116  void operator++()
117  {
118  ++sigma_eff_it;
119  ++sigma_d_it;
120  ++step_size_it;
121  }
122  double sigma_eff() const { return *sigma_eff_it; }
123  double sigma_d() const { return *sigma_d_it; }
124  double step_size() const { return *step_size_it; }
125  static double sqr(double x) { return x * x; }
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 = sqr(sigma_eff()) - sqr(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>
238 
239  This class enables the calculation of scale space convolutions
240  such as \ref gaussianGradientMultiArray() on data with anisotropic
241  discretization. For these, the result of the ordinary calculation
242  has to be multiplied by factors of \f$1/w^{n}\f$ for each dimension,
243  where \f$w\f$ is the step size of the grid in said dimension and
244  \f$n\f$ is the differential order of the convolution, e.g., 1 for
245  gaussianGradientMultiArray(), and 0 for gaussianSmoothMultiArray(),
246  respectively. Also for each dimension in turn, the convolution's scale
247  parameter \f$\sigma\f$ has to be replaced by
248  \f$\sqrt{\sigma_\mathrm{eff}^2 - \sigma_\mathrm{D}^2}\Big/w\f$,
249  where \f$\sigma_\mathrm{eff}\f$ is the resulting effective filtering
250  scale. The data is assumed to be already filtered by a
251  gaussian smoothing with the scale parameter \f$\sigma_\mathrm{D}\f$
252  (such as by measuring equipment). All of the above changes are
253  automatically employed by the convolution functions for <tt>MultiArray</tt>s
254  if a corresponding options object is provided.
255 
256  The <tt>ConvolutionOptions</tt> class must be parameterized by the dimension
257  <tt>dim</tt>
258  of the <tt>MultiArray</tt>s on which it is used. The actual per-axis
259  options are set by (overloaded) member functions explained below,
260  or else default to neutral values corresponding to the absence of the
261  particular option.
262 
263  All member functions set <tt>dim</tt> values of the respective convolution
264  option, one for each dimension. They may be set explicitly by multiple
265  arguments for up to five dimensions, or by a single argument to the same
266  value for all dimensions. For the general case, a single argument that is
267  either a C-syle array, an iterator, or a C++ standard library style
268  sequence (such as <tt>std::vector</tt>, with member functions <tt>begin()</tt>
269  and <tt>size()</tt>) supplies the option values for any number of dimensions.
270 
271  Note that the return value of all member functions is <tt>*this</tt>, which
272  provides the mechanism for concatenating member function calls as shown below.
273 
274  <b>usage with explicit parameters:</b>
275 
276  \code
277  ConvolutionOptions<2> opt = ConvolutionOptions<2>().stepSize(1, 2.3);
278  \endcode
279 
280  <b>usage with arrays:</b>
281 
282  \code
283  const double step_size[3] = { x_scale, y_scale, z_scale };
284  ConvolutionOptions<3> opt = ConvolutionOptions<3>().stepSize(step_size);
285  \endcode
286 
287  <b>usage with C++ standard library style sequences:</b>
288 
289  \code
290  TinyVector<double, 4> step_size(1, 1, 2.0, 1.5);
291  TinyVector<double, 4> r_sigmas(1, 1, 2.3, 3.2);
292  ConvolutionOptions<4> opt = ConvolutionOptions<4>().stepSize(step_size).resolutionStdDev(r_sigmas);
293  \endcode
294 
295  <b>usage with iterators:</b>
296 
297  \code
298  ArrayVector<double> step_size;
299  step_size.push_back(0);
300  step_size.push_back(3);
301  step_size.push_back(4);
302  ArrayVector<double>::iterator i = step_size.begin();
303  ++i;
304  ConvolutionOptions<2> opt = ConvolutionOptions<2>().stepSize(i);
305  \endcode
306 
307  <b>general usage in a convolution function call:</b>
308 
309  \code
310  MultiArray<3, double> test_image;
311  MultiArray<3, double> out_image;
312  gaussianSmoothMultiArray(srcMultiArrayRange(test_image),
313  destMultiArray(out_image),
314  5.0,
315  ConvolutionOptions<3>()
316  .stepSize (1, 1, 3.2)
317  .resolutionStdDev(1, 1, 4)
318  );
319  \endcode
320 
321 */
322 template <unsigned dim>
324 {
325  public:
326  typedef typename MultiArrayShape<dim>::type Shape;
327  typedef detail::multiArrayScaleParam<dim> ParamVec;
328  typedef typename ParamVec::return_type ParamIt;
329 
330  ParamVec sigma_eff;
331  ParamVec sigma_d;
332  ParamVec step_size;
333  ParamVec outer_scale;
334  double window_ratio;
335  Shape from_point, to_point;
336 
338  : sigma_eff(0.0),
339  sigma_d(0.0),
340  step_size(1.0),
341  outer_scale(0.0),
342  window_ratio(0.0)
343  {}
344 
345  typedef typename detail::WrapDoubleIteratorTriple<ParamIt, ParamIt, ParamIt>
346  ScaleIterator;
347  typedef typename detail::WrapDoubleIterator<ParamIt>::type
348  StepIterator;
349 
350  ScaleIterator scaleParams() const
351  {
352  return ScaleIterator(sigma_eff(), sigma_d(), step_size());
353  }
354  StepIterator stepParams() const
355  {
356  return StepIterator(step_size());
357  }
358 
359  ConvolutionOptions outerOptions() const
360  {
361  ConvolutionOptions outer = *this;
362  // backward-compatible values:
363  return outer.stdDev(outer_scale()).resolutionStdDev(0.0);
364  }
365 
366  // Step size per axis.
367  // Default: dim values of 1.0
368  VIGRA_CONVOLUTION_OPTIONS(stepSize, 1.0, step_size)
369 #ifdef DOXYGEN
370  /** Step size(s) per axis, i.e., the distance between two
371  adjacent pixels. Required for <tt>MultiArray</tt>
372  containing anisotropic data.
373 
374  Note that a convolution containing a derivative operator
375  of order <tt>n</tt> results in a multiplication by
376  \f${\rm stepSize}^{-n}\f$ for each axis.
377  Also, the above standard deviations
378  are scaled according to the step size of each axis.
379  Default value for the options object if this member function is not
380  used: A value of 1.0 for each dimension.
381  */
383 #endif
384 
385  // Resolution standard deviation per axis.
386  // Default: dim values of 0.0
387  VIGRA_CONVOLUTION_OPTIONS(resolutionStdDev, 0.0, sigma_d)
388 #ifdef DOXYGEN
389  /** Resolution standard deviation(s) per axis, i.e., a supposed
390  pre-existing gaussian filtering by this value.
391 
392  The standard deviation actually used by the convolution operators
393  is \f$\sqrt{{\rm sigma}^{2} - {\rm resolutionStdDev}^{2}}\f$ for each
394  axis.
395  Default value for the options object if this member function is not
396  used: A value of 0.0 for each dimension.
397  */
399 #endif
400 
401  // Standard deviation of scale space operators.
402  // Default: dim values of 0.0
403  VIGRA_CONVOLUTION_OPTIONS(stdDev, 0.0, sigma_eff)
404  VIGRA_CONVOLUTION_OPTIONS(innerScale, 0.0, sigma_eff)
405 
406 #ifdef DOXYGEN
407  /** Standard deviation(s) of scale space operators, or inner scale(s) for \ref structureTensorMultiArray().
408 
409  Usually not
410  needed, since a single value for all axes may be specified as a parameter
411  <tt>sigma</tt> to the call of
412  an convolution operator such as \ref gaussianGradientMultiArray(), and
413  anisotropic data requiring the use of the stepSize() member function.
414  Default value for the options object if this member function is not
415  used: A value of 0.0 for each dimension.
416  */
418 
419  /** Standard deviation(s) of scale space operators, or inner scale(s) for \ref structureTensorMultiArray().
420 
421  Usually not
422  needed, since a single value for all axes may be specified as a parameter
423  <tt>sigma</tt> to the call of
424  an convolution operator such as \ref gaussianGradientMultiArray(), and
425  anisotropic data requiring the use of the stepSize() member function.
426  Default value for the options object if this member function is not
427  used: A value of 0.0 for each dimension.
428  */
430 #endif
431 
432  // Outer scale, for structure tensor.
433  // Default: dim values of 0.0
434  VIGRA_CONVOLUTION_OPTIONS(outerScale, 0.0, outer_scale)
435 #ifdef DOXYGEN
436  /** Standard deviation(s) of the second convolution of the
437  structure tensor.
438 
439  Usually not needed, since a single value for
440  all axes may be specified as a parameter <tt>outerScale</tt> to
441  the call of \ref structureTensorMultiArray(), and
442  anisotropic data requiring the use of the stepSize() member
443  function.
444  Default value for the options object if this member function is not
445  used: A value of 0.0 for each dimension.
446  */
448 #endif
449 
450  /** Size of the filter window as a multiple of the scale parameter.
451 
452  This option is only used for Gaussian filters and their derivatives.
453  By default, the window size of a Gaussian filter is automatically
454  determined such that the error resulting from restricting the
455  infinitely large Gaussian function to a finite size is minimized.
456  In particular, the window radius is determined as
457  <tt>radius = round(3.0 * sigma + 0.5 * order)</tt>, where 'order' is the
458  desired derivative order. In some cases, it is desirable to trade off
459  accuracy for speed, and this function can be used to request a smaller
460  window radius.
461 
462  Default: <tt>0.0</tt> (i.e. determine the window size automatically)
463  */
465  {
466  vigra_precondition(ratio >= 0.0,
467  "ConvolutionOptions::filterWindowSize(): ratio must not be negative.");
468  window_ratio = ratio;
469  return *this;
470  }
471 
472  /** Restrict the filter to a subregion of the input array.
473 
474  This is useful for speeding up computations by ignoring irrelevant
475  areas in the array. <b>Note:</b> It is assumed that the output array
476  of the convolution has the size given in this function.
477 
478  Default: <tt>from = Shape(), to = Shape()</tt> (i.e. use entire array)
479  */
480  ConvolutionOptions<dim> & subarray(Shape const & from, Shape const & to)
481  {
482  from_point = from;
483  to_point = to;
484  return *this;
485  }
486 };
487 
488 namespace detail
489 {
490 
491 /********************************************************/
492 /* */
493 /* internalSeparableConvolveMultiArray */
494 /* */
495 /********************************************************/
496 
497 template <class SrcIterator, class SrcShape, class SrcAccessor,
498  class DestIterator, class DestAccessor, class KernelIterator>
499 void
500 internalSeparableConvolveMultiArrayTmp(
501  SrcIterator si, SrcShape const & shape, SrcAccessor src,
502  DestIterator di, DestAccessor dest, KernelIterator kit)
503 {
504  enum { N = 1 + SrcIterator::level };
505 
506  typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
507  typedef typename AccessorTraits<TmpType>::default_accessor TmpAcessor;
508 
509  // temporary array to hold the current line to enable in-place operation
510  ArrayVector<TmpType> tmp( shape[0] );
511 
512  typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
513  typedef MultiArrayNavigator<DestIterator, N> DNavigator;
514 
515  TmpAcessor acc;
516 
517  {
518  // only operate on first dimension here
519  SNavigator snav( si, shape, 0 );
520  DNavigator dnav( di, shape, 0 );
521 
522  for( ; snav.hasMore(); snav++, dnav++ )
523  {
524  // first copy source to tmp for maximum cache efficiency
525  copyLine(snav.begin(), snav.end(), src, tmp.begin(), acc);
526 
527  convolveLine(srcIterRange(tmp.begin(), tmp.end(), acc),
528  destIter( dnav.begin(), dest ),
529  kernel1d( *kit ) );
530  }
531  ++kit;
532  }
533 
534  // operate on further dimensions
535  for( int d = 1; d < N; ++d, ++kit )
536  {
537  DNavigator dnav( di, shape, d );
538 
539  tmp.resize( shape[d] );
540 
541  for( ; dnav.hasMore(); dnav++ )
542  {
543  // first copy source to tmp since convolveLine() cannot work in-place
544  copyLine(dnav.begin(), dnav.end(), dest, tmp.begin(), acc);
545 
546  convolveLine(srcIterRange(tmp.begin(), tmp.end(), acc),
547  destIter( dnav.begin(), dest ),
548  kernel1d( *kit ) );
549  }
550  }
551 }
552 
553 /********************************************************/
554 /* */
555 /* internalSeparableConvolveSubarray */
556 /* */
557 /********************************************************/
558 
559 template <class SrcIterator, class SrcShape, class SrcAccessor,
560  class DestIterator, class DestAccessor, class KernelIterator>
561 void
562 internalSeparableConvolveSubarray(
563  SrcIterator si, SrcShape const & shape, SrcAccessor src,
564  DestIterator di, DestAccessor dest, KernelIterator kit,
565  SrcShape const & start, SrcShape const & stop)
566 {
567  enum { N = 1 + SrcIterator::level };
568 
569  typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
570  typedef MultiArray<N, TmpType> TmpArray;
571  typedef typename TmpArray::traverser TmpIterator;
572  typedef typename AccessorTraits<TmpType>::default_accessor TmpAcessor;
573 
574  SrcShape sstart, sstop, axisorder, tmpshape;
575  TinyVector<double, N> overhead;
576  for(int k=0; k<N; ++k)
577  {
578  sstart[k] = start[k] - kit[k].right();
579  if(sstart[k] < 0)
580  sstart[k] = 0;
581  sstop[k] = stop[k] - kit[k].left();
582  if(sstop[k] > shape[k])
583  sstop[k] = shape[k];
584  overhead[k] = double(sstop[k] - sstart[k]) / (stop[k] - start[k]);
585  }
586 
587  indexSort(overhead.begin(), overhead.end(), axisorder.begin(), std::greater<double>());
588 
589  SrcShape dstart, dstop(sstop - sstart);
590  dstop[axisorder[0]] = stop[axisorder[0]] - start[axisorder[0]];
591 
592  // temporary array to hold the current line to enable in-place operation
593  MultiArray<N, TmpType> tmp(dstop);
594 
595  typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
596  typedef MultiArrayNavigator<TmpIterator, N> TNavigator;
597  typedef MultiArrayNavigator<DestIterator, N> DNavigator;
598 
599  TmpAcessor acc;
600 
601  {
602  // only operate on first dimension here
603  SNavigator snav( si, sstart, sstop, axisorder[0]);
604  TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[0]);
605 
606  ArrayVector<TmpType> tmpline(sstop[axisorder[0]] - sstart[axisorder[0]]);
607 
608  int lstart = start[axisorder[0]] - sstart[axisorder[0]];
609  int lstop = lstart + (stop[axisorder[0]] - start[axisorder[0]]);
610 
611  for( ; snav.hasMore(); snav++, tnav++ )
612  {
613  // first copy source to tmp for maximum cache efficiency
614  copyLine(snav.begin(), snav.end(), src, tmpline.begin(), acc);
615 
616  convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
617  destIter(tnav.begin(), acc),
618  kernel1d( kit[axisorder[0]] ), lstart, lstop);
619  }
620  }
621 
622  // operate on further dimensions
623  for( int d = 1; d < N; ++d)
624  {
625  TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[d]);
626 
627  ArrayVector<TmpType> tmpline(dstop[axisorder[d]] - dstart[axisorder[d]]);
628 
629  int lstart = start[axisorder[d]] - sstart[axisorder[d]];
630  int lstop = lstart + (stop[axisorder[d]] - start[axisorder[d]]);
631 
632  for( ; tnav.hasMore(); tnav++ )
633  {
634  // first copy source to tmp because convolveLine() cannot work in-place
635  copyLine(tnav.begin(), tnav.end(), acc, tmpline.begin(), acc );
636 
637  convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
638  destIter( tnav.begin() + lstart, acc ),
639  kernel1d( kit[axisorder[d]] ), lstart, lstop);
640  }
641 
642  dstart[axisorder[d]] = lstart;
643  dstop[axisorder[d]] = lstop;
644  }
645 
646  copyMultiArray(tmp.traverser_begin()+dstart, stop-start, acc, di, dest);
647 }
648 
649 
650 template <class K>
651 void
652 scaleKernel(K & kernel, double a)
653 {
654  for(int i = kernel.left(); i <= kernel.right(); ++i)
655  kernel[i] = detail::RequiresExplicitCast<typename K::value_type>::cast(kernel[i] * a);
656 }
657 
658 
659 } // namespace detail
660 
661 /** \addtogroup MultiArrayConvolutionFilters Convolution filters for multi-dimensional arrays.
662 
663  These functions realize a separable convolution on an arbitrary dimensional
664  array that is specified by iterators (compatible to \ref MultiIteratorPage)
665  and shape objects. It can therefore be applied to a wide range of data structures
666  (\ref vigra::MultiArrayView, \ref vigra::MultiArray etc.).
667 */
668 //@{
669 
670 /********************************************************/
671 /* */
672 /* separableConvolveMultiArray */
673 /* */
674 /********************************************************/
675 
676 /** \brief Separated convolution on multi-dimensional arrays.
677 
678  This function computes a separated convolution on all dimensions
679  of the given multi-dimensional array. Both source and destination
680  arrays are represented by iterators, shape objects and accessors.
681  The destination array is required to already have the correct size.
682 
683  There are two variants of this functions: one takes a single kernel
684  of type \ref vigra::Kernel1D which is then applied to all dimensions,
685  whereas the other requires an iterator referencing a sequence of
686  \ref vigra::Kernel1D objects, one for every dimension of the data.
687  Then the first kernel in this sequence is applied to the innermost
688  dimension (e.g. the x-dimension of an image), while the last is applied to the
689  outermost dimension (e.g. the z-dimension in a 3D image).
690 
691  This function may work in-place, which means that <tt>siter == diter</tt> is allowed.
692  A full-sized internal array is only allocated if working on the destination
693  array directly would cause round-off errors (i.e. if
694  <tt>typeid(typename NumericTraits<typename DestAccessor::value_type>::RealPromote)
695  != typeid(typename DestAccessor::value_type)</tt>.
696 
697  If <tt>start</tt> and <tt>stop</tt> have non-default values, they must represent
698  a valid subarray of the input array. The convolution is then restricted to that
699  subarray, and it is assumed that the output array only refers to the
700  subarray (i.e. <tt>diter</tt> points to the element corresponding to
701  <tt>start</tt>).
702 
703  <b> Declarations:</b>
704 
705  pass arguments explicitly:
706  \code
707  namespace vigra {
708  // apply the same kernel to all dimensions
709  template <class SrcIterator, class SrcShape, class SrcAccessor,
710  class DestIterator, class DestAccessor, class T>
711  void
712  separableConvolveMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
713  DestIterator diter, DestAccessor dest,
714  Kernel1D<T> const & kernel,
715  SrcShape const & start = SrcShape(),
716  SrcShape const & stop = SrcShape());
717 
718  // apply each kernel from the sequence 'kernels' in turn
719  template <class SrcIterator, class SrcShape, class SrcAccessor,
720  class DestIterator, class DestAccessor, class KernelIterator>
721  void
722  separableConvolveMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
723  DestIterator diter, DestAccessor dest,
724  KernelIterator kernels,
725  SrcShape const & start = SrcShape(),
726  SrcShape const & stop = SrcShape());
727  }
728  \endcode
729 
730  use argument objects in conjunction with \ref ArgumentObjectFactories :
731  \code
732  namespace vigra {
733  // apply the same kernel to all dimensions
734  template <class SrcIterator, class SrcShape, class SrcAccessor,
735  class DestIterator, class DestAccessor, class T>
736  void
737  separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
738  pair<DestIterator, DestAccessor> const & dest,
739  Kernel1D<T> const & kernel,
740  SrcShape const & start = SrcShape(),
741  SrcShape const & stop = SrcShape());
742 
743  // apply each kernel from the sequence 'kernels' in turn
744  template <class SrcIterator, class SrcShape, class SrcAccessor,
745  class DestIterator, class DestAccessor, class KernelIterator>
746  void
747  separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
748  pair<DestIterator, DestAccessor> const & dest,
749  KernelIterator kernels,
750  SrcShape const & start = SrcShape(),
751  SrcShape const & stop = SrcShape());
752  }
753  \endcode
754 
755  <b> Usage:</b>
756 
757  <b>\#include</b> <vigra/multi_convolution.hxx>
758 
759  \code
760  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
761  MultiArray<3, unsigned char> source(shape);
762  MultiArray<3, float> dest(shape);
763  ...
764  Kernel1D<float> gauss;
765  gauss.initGaussian(sigma);
766  // create 3 Gauss kernels, one for each dimension
767  ArrayVector<Kernel1D<float> > kernels(3, gauss);
768 
769  // perform Gaussian smoothing on all dimensions
770  separableConvolveMultiArray(srcMultiArrayRange(source), destMultiArray(dest),
771  kernels.begin());
772  \endcode
773 
774  <b> Required Interface:</b>
775 
776  see \ref separableConvolveMultiArray(), in addition:
777 
778  \code
779  int dimension = 0;
780  VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
781  \endcode
782 
783  \see vigra::Kernel1D, convolveLine()
784 */
785 doxygen_overloaded_function(template <...> void separableConvolveMultiArray)
786 
787 template <class SrcIterator, class SrcShape, class SrcAccessor,
788  class DestIterator, class DestAccessor, class KernelIterator>
789 void
790 separableConvolveMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
791  DestIterator d, DestAccessor dest,
792  KernelIterator kernels,
793  SrcShape const & start = SrcShape(),
794  SrcShape const & stop = SrcShape())
795 {
796  typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
797 
798  if(stop != SrcShape())
799  {
800  enum { N = 1 + SrcIterator::level };
801  for(int k=0; k<N; ++k)
802  vigra_precondition(0 <= start[k] && start[k] < stop[k] && stop[k] <= shape[k],
803  "separableConvolveMultiArray(): invalid subarray shape.");
804 
805  detail::internalSeparableConvolveSubarray(s, shape, src, d, dest, kernels, start, stop);
806  }
807  else if(!IsSameType<TmpType, typename DestAccessor::value_type>::boolResult)
808  {
809  // need a temporary array to avoid rounding errors
810  MultiArray<SrcShape::static_size, TmpType> tmpArray(shape);
811  detail::internalSeparableConvolveMultiArrayTmp( s, shape, src,
812  tmpArray.traverser_begin(), typename AccessorTraits<TmpType>::default_accessor(), kernels );
813  copyMultiArray(srcMultiArrayRange(tmpArray), destIter(d, dest));
814  }
815  else
816  {
817  // work directly on the destination array
818  detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, d, dest, kernels );
819  }
820 }
821 
822 template <class SrcIterator, class SrcShape, class SrcAccessor,
823  class DestIterator, class DestAccessor, class KernelIterator>
824 inline
826  triple<SrcIterator, SrcShape, SrcAccessor> const & source,
827  pair<DestIterator, DestAccessor> const & dest,
828  KernelIterator kit,
829  SrcShape const & start = SrcShape(),
830  SrcShape const & stop = SrcShape())
831 {
832  separableConvolveMultiArray( source.first, source.second, source.third,
833  dest.first, dest.second, kit, start, stop );
834 }
835 
836 template <class SrcIterator, class SrcShape, class SrcAccessor,
837  class DestIterator, class DestAccessor, class T>
838 inline void
839 separableConvolveMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
840  DestIterator d, DestAccessor dest,
841  Kernel1D<T> const & kernel,
842  SrcShape const & start = SrcShape(),
843  SrcShape const & stop = SrcShape())
844 {
845  ArrayVector<Kernel1D<T> > kernels(shape.size(), kernel);
846 
847  separableConvolveMultiArray( s, shape, src, d, dest, kernels.begin(), start, stop);
848 }
849 
850 template <class SrcIterator, class SrcShape, class SrcAccessor,
851  class DestIterator, class DestAccessor, class T>
852 inline void
853 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
854  pair<DestIterator, DestAccessor> const & dest,
855  Kernel1D<T> const & kernel,
856  SrcShape const & start = SrcShape(),
857  SrcShape const & stop = SrcShape())
858 {
859  ArrayVector<Kernel1D<T> > kernels(source.second.size(), kernel);
860 
861  separableConvolveMultiArray( source.first, source.second, source.third,
862  dest.first, dest.second, kernels.begin(), start, stop);
863 }
864 
865 /********************************************************/
866 /* */
867 /* convolveMultiArrayOneDimension */
868 /* */
869 /********************************************************/
870 
871 /** \brief Convolution along a single dimension of a multi-dimensional arrays.
872 
873  This function computes a convolution along one dimension (specified by
874  the parameter <tt>dim</tt> of the given multi-dimensional array with the given
875  <tt>kernel</tt>. Both source and destination arrays are represented by
876  iterators, shape objects and accessors. The destination array is required to
877  already have the correct size.
878 
879  If <tt>start</tt> and <tt>stop</tt> have non-default values, they must represent
880  a valid subarray of the input array. The convolution is then restricted to that
881  subarray, and it is assumed that the output array only refers to the
882  subarray (i.e. <tt>diter</tt> points to the element corresponding to
883  <tt>start</tt>).
884 
885  This function may work in-place, which means that <tt>siter == diter</tt> is allowed.
886 
887  <b> Declarations:</b>
888 
889  pass arguments explicitly:
890  \code
891  namespace vigra {
892  template <class SrcIterator, class SrcShape, class SrcAccessor,
893  class DestIterator, class DestAccessor, class T>
894  void
895  convolveMultiArrayOneDimension(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
896  DestIterator diter, DestAccessor dest,
897  unsigned int dim, vigra::Kernel1D<T> const & kernel,
898  SrcShape const & start = SrcShape(),
899  SrcShape const & stop = SrcShape());
900  }
901  \endcode
902 
903  use argument objects in conjunction with \ref ArgumentObjectFactories :
904  \code
905  namespace vigra {
906  template <class SrcIterator, class SrcShape, class SrcAccessor,
907  class DestIterator, class DestAccessor, class T>
908  void
909  convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
910  pair<DestIterator, DestAccessor> const & dest,
911  unsigned int dim, vigra::Kernel1D<T> const & kernel,
912  SrcShape const & start = SrcShape(),
913  SrcShape const & stop = SrcShape());
914  }
915  \endcode
916 
917  <b> Usage:</b>
918 
919  <b>\#include</b> <vigra/multi_convolution.hxx>
920 
921  \code
922  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
923  MultiArray<3, unsigned char> source(shape);
924  MultiArray<3, float> dest(shape);
925  ...
926  Kernel1D<float> gauss;
927  gauss.initGaussian(sigma);
928 
929  // perform Gaussian smoothing along dimensions 1 (height)
930  convolveMultiArrayOneDimension(srcMultiArrayRange(source), destMultiArray(dest), 1, gauss);
931  \endcode
932 
933  \see separableConvolveMultiArray()
934 */
935 doxygen_overloaded_function(template <...> void convolveMultiArrayOneDimension)
936 
937 template <class SrcIterator, class SrcShape, class SrcAccessor,
938  class DestIterator, class DestAccessor, class T>
939 void
940 convolveMultiArrayOneDimension(SrcIterator s, SrcShape const & shape, SrcAccessor src,
941  DestIterator d, DestAccessor dest,
942  unsigned int dim, vigra::Kernel1D<T> const & kernel,
943  SrcShape const & start = SrcShape(),
944  SrcShape const & stop = SrcShape())
945 {
946  enum { N = 1 + SrcIterator::level };
947  vigra_precondition( dim < N,
948  "convolveMultiArrayOneDimension(): The dimension number to convolve must be smaller "
949  "than the data dimensionality" );
950 
951  typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
952  typedef typename AccessorTraits<TmpType>::default_const_accessor TmpAccessor;
953  ArrayVector<TmpType> tmp( shape[dim] );
954 
955  typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
956  typedef MultiArrayNavigator<DestIterator, N> DNavigator;
957 
958  SrcShape sstart, sstop(shape), dstart, dstop(shape);
959 
960  if(stop != SrcShape())
961  {
962  sstart = start;
963  sstop = stop;
964  sstart[dim] = 0;
965  sstop[dim] = shape[dim];
966  dstop = stop - start;
967  }
968 
969  SNavigator snav( s, sstart, sstop, dim );
970  DNavigator dnav( d, dstart, dstop, dim );
971 
972  for( ; snav.hasMore(); snav++, dnav++ )
973  {
974  // first copy source to temp for maximum cache efficiency
975  copyLine( snav.begin(), snav.end(), src,
976  tmp.begin(), typename AccessorTraits<TmpType>::default_accessor() );
977 
978  convolveLine(srcIterRange( tmp.begin(), tmp.end(), TmpAccessor()),
979  destIter( dnav.begin(), dest ),
980  kernel1d( kernel), start[dim], stop[dim]);
981  }
982 }
983 
984 template <class SrcIterator, class SrcShape, class SrcAccessor,
985  class DestIterator, class DestAccessor, class T>
986 inline void
987 convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
988  pair<DestIterator, DestAccessor> const & dest,
989  unsigned int dim, vigra::Kernel1D<T> const & kernel,
990  SrcShape const & start = SrcShape(),
991  SrcShape const & stop = SrcShape())
992 {
993  convolveMultiArrayOneDimension(source.first, source.second, source.third,
994  dest.first, dest.second, dim, kernel, start, stop);
995 }
996 
997 /********************************************************/
998 /* */
999 /* gaussianSmoothMultiArray */
1000 /* */
1001 /********************************************************/
1002 
1003 /** \brief Isotropic Gaussian smoothing of a multi-dimensional arrays.
1004 
1005  This function computes an isotropic convolution of the given N-dimensional
1006  array with a Gaussian filter at the given standard deviation <tt>sigma</tt>.
1007  Both source and destination arrays are represented by
1008  iterators, shape objects and accessors. The destination array is required to
1009  already have the correct size. This function may work in-place, which means
1010  that <tt>siter == diter</tt> is allowed. It is implemented by a call to
1011  \ref separableConvolveMultiArray() with the appropriate kernel.
1012 
1013  Anisotropic data should be passed with appropriate
1014  \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional
1015  unless the parameter <tt>sigma</tt> is left out.
1016 
1017  <b> Declarations:</b>
1018 
1019  pass arguments explicitly:
1020  \code
1021  namespace vigra {
1022  template <class SrcIterator, class SrcShape, class SrcAccessor,
1023  class DestIterator, class DestAccessor>
1024  void
1025  gaussianSmoothMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1026  DestIterator diter, DestAccessor dest,
1027  double sigma, const ConvolutionOptions<N> & opt);
1028  }
1029  \endcode
1030 
1031  use argument objects in conjunction with \ref ArgumentObjectFactories :
1032  \code
1033  namespace vigra {
1034  template <class SrcIterator, class SrcShape, class SrcAccessor,
1035  class DestIterator, class DestAccessor>
1036  void
1037  gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1038  pair<DestIterator, DestAccessor> const & dest,
1039  double sigma, const ConvolutionOptions<N> & opt);
1040  }
1041  \endcode
1042 
1043  <b> Usage:</b>
1044 
1045  <b>\#include</b> <vigra/multi_convolution.hxx>
1046 
1047  \code
1048  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
1049  MultiArray<3, unsigned char> source(shape);
1050  MultiArray<3, float> dest(shape);
1051  ...
1052  // perform isotropic Gaussian smoothing at scale 'sigma'
1053  gaussianSmoothMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma);
1054  \endcode
1055 
1056  <b> Usage with anisotropic data:</b>
1057 
1058  <b>\#include</b> <vigra/multi_convolution.hxx>
1059 
1060  \code
1061  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
1062  MultiArray<3, unsigned char> source(shape);
1063  MultiArray<3, float> dest(shape);
1064  TinyVector<float, 3> step_size;
1065  TinyVector<float, 3> resolution_sigmas;
1066  ...
1067  // perform anisotropic Gaussian smoothing at scale 'sigma'
1068  gaussianSmoothMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma,
1069  ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
1070  \endcode
1071 
1072  \see separableConvolveMultiArray()
1073 */
1074 doxygen_overloaded_function(template <...> void gaussianSmoothMultiArray)
1075 
1076 template <class SrcIterator, class SrcShape, class SrcAccessor,
1077  class DestIterator, class DestAccessor>
1078 void
1079 gaussianSmoothMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
1080  DestIterator d, DestAccessor dest,
1081  const ConvolutionOptions<SrcShape::static_size> & opt,
1082  const char *const function_name = "gaussianSmoothMultiArray" )
1083 {
1084  typedef typename DestAccessor::value_type DestType;
1085 
1086  static const int N = SrcShape::static_size;
1087 
1088  typename ConvolutionOptions<N>::ScaleIterator params = opt.scaleParams();
1089  ArrayVector<Kernel1D<double> > kernels(N);
1090 
1091  for (int dim = 0; dim < N; ++dim, ++params)
1092  kernels[dim].initGaussian(params.sigma_scaled(function_name), 1.0, opt.window_ratio);
1093 
1094  separableConvolveMultiArray(s, shape, src, d, dest, kernels.begin(), opt.from_point, opt.to_point);
1095 }
1096 
1097 template <class SrcIterator, class SrcShape, class SrcAccessor,
1098  class DestIterator, class DestAccessor>
1099 inline void
1100 gaussianSmoothMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
1101  DestIterator d, DestAccessor dest, double sigma,
1102  const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1103 {
1104  ConvolutionOptions<SrcShape::static_size> par = opt;
1105  gaussianSmoothMultiArray(s, shape, src, d, dest, par.stdDev(sigma));
1106 }
1107 
1108 template <class SrcIterator, class SrcShape, class SrcAccessor,
1109  class DestIterator, class DestAccessor>
1110 inline void
1111 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1112  pair<DestIterator, DestAccessor> const & dest,
1113  const ConvolutionOptions<SrcShape::static_size> & opt)
1114 {
1115  gaussianSmoothMultiArray( source.first, source.second, source.third,
1116  dest.first, dest.second, opt );
1117 }
1118 
1119 template <class SrcIterator, class SrcShape, class SrcAccessor,
1120  class DestIterator, class DestAccessor>
1121 inline void
1122 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1123  pair<DestIterator, DestAccessor> const & dest, double sigma,
1124  const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1125 {
1126  gaussianSmoothMultiArray( source.first, source.second, source.third,
1127  dest.first, dest.second, sigma, opt );
1128 }
1129 
1130 
1131 /********************************************************/
1132 /* */
1133 /* gaussianGradientMultiArray */
1134 /* */
1135 /********************************************************/
1136 
1137 /** \brief Calculate Gaussian gradient of a multi-dimensional arrays.
1138 
1139  This function computes the Gaussian gradient of the given N-dimensional
1140  array with a sequence of first-derivative-of-Gaussian filters at the given
1141  standard deviation <tt>sigma</tt> (differentiation is applied to each dimension
1142  in turn, starting with the innermost dimension). Both source and destination arrays
1143  are represented by iterators, shape objects and accessors. The destination array is
1144  required to have a vector valued pixel type with as many elements as the number of
1145  dimensions. This function is implemented by calls to
1146  \ref separableConvolveMultiArray() with the appropriate kernels.
1147 
1148  Anisotropic data should be passed with appropriate
1149  \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional
1150  unless the parameter <tt>sigma</tt> is left out.
1151 
1152  <b> Declarations:</b>
1153 
1154  pass arguments explicitly:
1155  \code
1156  namespace vigra {
1157  template <class SrcIterator, class SrcShape, class SrcAccessor,
1158  class DestIterator, class DestAccessor>
1159  void
1160  gaussianGradientMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1161  DestIterator diter, DestAccessor dest,
1162  double sigma, const ConvolutionOptions<N> & opt);
1163  }
1164  \endcode
1165 
1166  use argument objects in conjunction with \ref ArgumentObjectFactories :
1167  \code
1168  namespace vigra {
1169  template <class SrcIterator, class SrcShape, class SrcAccessor,
1170  class DestIterator, class DestAccessor>
1171  void
1172  gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1173  pair<DestIterator, DestAccessor> const & dest,
1174  double sigma, const ConvolutionOptions<N> & opt);
1175  }
1176  \endcode
1177 
1178  <b> Usage:</b>
1179 
1180  <b>\#include</b> <vigra/multi_convolution.hxx>
1181 
1182  \code
1183  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
1184  MultiArray<3, unsigned char> source(shape);
1185  MultiArray<3, TinyVector<float, 3> > dest(shape);
1186  ...
1187  // compute Gaussian gradient at scale sigma
1188  gaussianGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma);
1189  \endcode
1190 
1191  <b> Usage with anisotropic data:</b>
1192 
1193  <b>\#include</b> <vigra/multi_convolution.hxx>
1194 
1195  \code
1196  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
1197  MultiArray<3, unsigned char> source(shape);
1198  MultiArray<3, TinyVector<float, 3> > dest(shape);
1199  TinyVector<float, 3> step_size;
1200  TinyVector<float, 3> resolution_sigmas;
1201  ...
1202  // compute Gaussian gradient at scale sigma
1203  gaussianGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma,
1204  ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
1205  \endcode
1206 
1207  <b> Required Interface:</b>
1208 
1209  see \ref separableConvolveMultiArray(), in addition:
1210 
1211  \code
1212  int dimension = 0;
1213  VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
1214  \endcode
1215 
1216  \see separableConvolveMultiArray()
1217 */
1218 doxygen_overloaded_function(template <...> void gaussianGradientMultiArray)
1219 
1220 template <class SrcIterator, class SrcShape, class SrcAccessor,
1221  class DestIterator, class DestAccessor>
1222 void
1223 gaussianGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1224  DestIterator di, DestAccessor dest,
1225  ConvolutionOptions<SrcShape::static_size> const & opt,
1226  const char *const function_name = "gaussianGradientMultiArray")
1227 {
1228  typedef typename DestAccessor::value_type DestType;
1229  typedef typename DestType::value_type DestValueType;
1230  typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1231 
1232  static const int N = SrcShape::static_size;
1233  typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
1234 
1235  for(int k=0; k<N; ++k)
1236  if(shape[k] <=0)
1237  return;
1238 
1239  vigra_precondition(N == (int)dest.size(di),
1240  "gaussianGradientMultiArray(): Wrong number of channels in output array.");
1241 
1242  ParamType params = opt.scaleParams();
1243  ParamType params2(params);
1244 
1245  ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
1246  for (int dim = 0; dim < N; ++dim, ++params)
1247  {
1248  double sigma = params.sigma_scaled(function_name);
1249  plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
1250  }
1251 
1252  typedef VectorElementAccessor<DestAccessor> ElementAccessor;
1253 
1254  // compute gradient components
1255  for (int dim = 0; dim < N; ++dim, ++params2)
1256  {
1257  ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
1258  kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 1, 1.0, opt.window_ratio);
1259  detail::scaleKernel(kernels[dim], 1.0 / params2.step_size());
1260  separableConvolveMultiArray(si, shape, src, di, ElementAccessor(dim, dest), kernels.begin(),
1261  opt.from_point, opt.to_point);
1262  }
1263 }
1264 
1265 template <class SrcIterator, class SrcShape, class SrcAccessor,
1266  class DestIterator, class DestAccessor>
1267 void
1268 gaussianGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1269  DestIterator di, DestAccessor dest, double sigma,
1270  const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1271 {
1272  ConvolutionOptions<SrcShape::static_size> par = opt;
1273  gaussianGradientMultiArray(si, shape, src, di, dest, par.stdDev(sigma));
1274 }
1275 
1276 template <class SrcIterator, class SrcShape, class SrcAccessor,
1277  class DestIterator, class DestAccessor>
1278 inline void
1279 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1280  pair<DestIterator, DestAccessor> const & dest,
1281  ConvolutionOptions<SrcShape::static_size> const & opt )
1282 {
1283  gaussianGradientMultiArray( source.first, source.second, source.third,
1284  dest.first, dest.second, opt );
1285 }
1286 
1287 template <class SrcIterator, class SrcShape, class SrcAccessor,
1288  class DestIterator, class DestAccessor>
1289 inline void
1290 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1291  pair<DestIterator, DestAccessor> const & dest,
1292  double sigma,
1293  const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1294 {
1295  gaussianGradientMultiArray( source.first, source.second, source.third,
1296  dest.first, dest.second, sigma, opt );
1297 }
1298 
1299 /********************************************************/
1300 /* */
1301 /* symmetricGradientMultiArray */
1302 /* */
1303 /********************************************************/
1304 
1305 /** \brief Calculate gradient of a multi-dimensional arrays using symmetric difference filters.
1306 
1307  This function computes the gradient of the given N-dimensional
1308  array with a sequence of symmetric difference filters a (differentiation is applied
1309  to each dimension in turn, starting with the innermost dimension). Both source and
1310  destination arrays are represented by iterators, shape objects and accessors.
1311  The destination array is required to have a vector valued pixel type with as many
1312  elements as the number of dimensions. This function is implemented by calls to
1313  \ref convolveMultiArrayOneDimension() with the symmetric difference kernel.
1314 
1315  Anisotropic data should be passed with appropriate
1316  \ref ConvolutionOptions, the parameter <tt>opt</tt> is optional
1317  otherwise.
1318 
1319  <b> Declarations:</b>
1320 
1321  pass arguments explicitly:
1322  \code
1323  namespace vigra {
1324  template <class SrcIterator, class SrcShape, class SrcAccessor,
1325  class DestIterator, class DestAccessor>
1326  void
1327  symmetricGradientMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1328  DestIterator diter, DestAccessor dest,
1329  const ConvolutionOptions<N> & opt);
1330  }
1331  \endcode
1332 
1333  use argument objects in conjunction with \ref ArgumentObjectFactories :
1334  \code
1335  namespace vigra {
1336  template <class SrcIterator, class SrcShape, class SrcAccessor,
1337  class DestIterator, class DestAccessor>
1338  void
1339  symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1340  pair<DestIterator, DestAccessor> const & dest,
1341  const ConvolutionOptions<N> & opt);
1342  }
1343  \endcode
1344 
1345  <b> Usage:</b>
1346 
1347  <b>\#include</b> <vigra/multi_convolution.hxx>
1348 
1349  \code
1350  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
1351  MultiArray<3, unsigned char> source(shape);
1352  MultiArray<3, TinyVector<float, 3> > dest(shape);
1353  ...
1354  // compute gradient
1355  symmetricGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest));
1356  \endcode
1357 
1358  <b> Usage with anisotropic data:</b>
1359 
1360  <b>\#include</b> <vigra/multi_convolution.hxx>
1361 
1362  \code
1363  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
1364  MultiArray<3, unsigned char> source(shape);
1365  MultiArray<3, TinyVector<float, 3> > dest(shape);
1366  TinyVector<float, 3> step_size;
1367  ...
1368  // compute gradient
1369  symmetricGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest),
1370  ConvolutionOptions<3>().stepSize(step_size));
1371  \endcode
1372 
1373  <b> Required Interface:</b>
1374 
1375  see \ref convolveMultiArrayOneDimension(), in addition:
1376 
1377  \code
1378  int dimension = 0;
1379  VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
1380  \endcode
1381 
1382  \see convolveMultiArrayOneDimension()
1383 */
1384 doxygen_overloaded_function(template <...> void symmetricGradientMultiArray)
1385 
1386 template <class SrcIterator, class SrcShape, class SrcAccessor,
1387  class DestIterator, class DestAccessor>
1388 void
1389 symmetricGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1390  DestIterator di, DestAccessor dest,
1391  const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1392 {
1393  typedef typename DestAccessor::value_type DestType;
1394  typedef typename DestType::value_type DestValueType;
1395  typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1396 
1397  static const int N = SrcShape::static_size;
1398  typedef typename ConvolutionOptions<N>::StepIterator StepType;
1399 
1400  for(int k=0; k<N; ++k)
1401  if(shape[k] <=0)
1402  return;
1403 
1404  vigra_precondition(N == (int)dest.size(di),
1405  "symmetricGradientMultiArray(): Wrong number of channels in output array.");
1406 
1407  Kernel1D<KernelType> filter;
1408  filter.initSymmetricGradient();
1409 
1410  StepType step_size_it = opt.stepParams();
1411 
1412  typedef VectorElementAccessor<DestAccessor> ElementAccessor;
1413 
1414  // compute gradient components
1415  for (int d = 0; d < N; ++d, ++step_size_it)
1416  {
1417  Kernel1D<KernelType> symmetric(filter);
1418  detail::scaleKernel(symmetric, 1 / *step_size_it);
1419  convolveMultiArrayOneDimension(si, shape, src,
1420  di, ElementAccessor(d, dest),
1421  d, symmetric, opt.from_point, opt.to_point);
1422  }
1423 }
1424 
1425 template <class SrcIterator, class SrcShape, class SrcAccessor,
1426  class DestIterator, class DestAccessor>
1427 inline void
1428 symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1429  pair<DestIterator, DestAccessor> const & dest,
1430  const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1431 {
1432  symmetricGradientMultiArray(source.first, source.second, source.third,
1433  dest.first, dest.second, opt);
1434 }
1435 
1436 /********************************************************/
1437 /* */
1438 /* laplacianOfGaussianMultiArray */
1439 /* */
1440 /********************************************************/
1441 
1442 /** \brief Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters.
1443 
1444  This function computes the Laplacian the given N-dimensional
1445  array with a sequence of second-derivative-of-Gaussian filters at the given
1446  standard deviation <tt>sigma</tt>. Both source and destination arrays
1447  are represented by iterators, shape objects and accessors. Both source and destination
1448  arrays must have scalar value_type. This function is implemented by calls to
1449  \ref separableConvolveMultiArray() with the appropriate kernels, followed by summation.
1450 
1451  Anisotropic data should be passed with appropriate
1452  \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional
1453  unless the parameter <tt>sigma</tt> is left out.
1454 
1455  <b> Declarations:</b>
1456 
1457  pass arguments explicitly:
1458  \code
1459  namespace vigra {
1460  template <class SrcIterator, class SrcShape, class SrcAccessor,
1461  class DestIterator, class DestAccessor>
1462  void
1463  laplacianOfGaussianMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1464  DestIterator diter, DestAccessor dest,
1465  double sigma, const ConvolutionOptions<N> & opt);
1466  }
1467  \endcode
1468 
1469  use argument objects in conjunction with \ref ArgumentObjectFactories :
1470  \code
1471  namespace vigra {
1472  template <class SrcIterator, class SrcShape, class SrcAccessor,
1473  class DestIterator, class DestAccessor>
1474  void
1475  laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1476  pair<DestIterator, DestAccessor> const & dest,
1477  double sigma, const ConvolutionOptions<N> & opt);
1478  }
1479  \endcode
1480 
1481  <b> Usage:</b>
1482 
1483  <b>\#include</b> <vigra/multi_convolution.hxx>
1484 
1485  \code
1486  MultiArray<3, float> source(shape);
1487  MultiArray<3, float> laplacian(shape);
1488  ...
1489  // compute Laplacian at scale sigma
1490  laplacianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(laplacian), sigma);
1491  \endcode
1492 
1493  <b> Usage with anisotropic data:</b>
1494 
1495  <b>\#include</b> <vigra/multi_convolution.hxx>
1496 
1497  \code
1498  MultiArray<3, float> source(shape);
1499  MultiArray<3, float> laplacian(shape);
1500  TinyVector<float, 3> step_size;
1501  TinyVector<float, 3> resolution_sigmas;
1502  ...
1503  // compute Laplacian at scale sigma
1504  laplacianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(laplacian), sigma,
1505  ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
1506  \endcode
1507 
1508  <b> Required Interface:</b>
1509 
1510  see \ref separableConvolveMultiArray(), in addition:
1511 
1512  \code
1513  int dimension = 0;
1514  VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
1515  \endcode
1516 
1517  \see separableConvolveMultiArray()
1518 */
1519 doxygen_overloaded_function(template <...> void laplacianOfGaussianMultiArray)
1520 
1521 template <class SrcIterator, class SrcShape, class SrcAccessor,
1522  class DestIterator, class DestAccessor>
1523 void
1524 laplacianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1525  DestIterator di, DestAccessor dest,
1526  ConvolutionOptions<SrcShape::static_size> const & opt )
1527 {
1528  using namespace functor;
1529 
1530  typedef typename DestAccessor::value_type DestType;
1531  typedef typename NumericTraits<DestType>::RealPromote KernelType;
1532  typedef typename AccessorTraits<KernelType>::default_accessor DerivativeAccessor;
1533 
1534  static const int N = SrcShape::static_size;
1535  typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
1536 
1537  ParamType params = opt.scaleParams();
1538  ParamType params2(params);
1539 
1540  ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
1541  for (int dim = 0; dim < N; ++dim, ++params)
1542  {
1543  double sigma = params.sigma_scaled("laplacianOfGaussianMultiArray");
1544  plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
1545  }
1546 
1547  SrcShape dshape(shape);
1548  if(opt.to_point != SrcShape())
1549  dshape = opt.to_point - opt.from_point;
1550 
1551  MultiArray<N, KernelType> derivative(dshape);
1552 
1553  // compute 2nd derivatives and sum them up
1554  for (int dim = 0; dim < N; ++dim, ++params2)
1555  {
1556  ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
1557  kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 2, 1.0, opt.window_ratio);
1558  detail::scaleKernel(kernels[dim], 1.0 / sq(params2.step_size()));
1559 
1560  if (dim == 0)
1561  {
1562  separableConvolveMultiArray( si, shape, src,
1563  di, dest, kernels.begin(), opt.from_point, opt.to_point);
1564  }
1565  else
1566  {
1567  separableConvolveMultiArray( si, shape, src,
1568  derivative.traverser_begin(), DerivativeAccessor(),
1569  kernels.begin(), opt.from_point, opt.to_point);
1570  combineTwoMultiArrays(di, dshape, dest, derivative.traverser_begin(), DerivativeAccessor(),
1571  di, dest, Arg1() + Arg2() );
1572  }
1573  }
1574 }
1575 
1576 template <class SrcIterator, class SrcShape, class SrcAccessor,
1577  class DestIterator, class DestAccessor>
1578 void
1579 laplacianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1580  DestIterator di, DestAccessor dest, double sigma,
1581  const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1582 {
1583  ConvolutionOptions<SrcShape::static_size> par = opt;
1584  laplacianOfGaussianMultiArray(si, shape, src, di, dest, par.stdDev(sigma));
1585 }
1586 
1587 template <class SrcIterator, class SrcShape, class SrcAccessor,
1588  class DestIterator, class DestAccessor>
1589 inline void
1590 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1591  pair<DestIterator, DestAccessor> const & dest,
1592  ConvolutionOptions<SrcShape::static_size> const & opt )
1593 {
1594  laplacianOfGaussianMultiArray( source.first, source.second, source.third,
1595  dest.first, dest.second, opt );
1596 }
1597 
1598 template <class SrcIterator, class SrcShape, class SrcAccessor,
1599  class DestIterator, class DestAccessor>
1600 inline void
1601 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1602  pair<DestIterator, DestAccessor> const & dest,
1603  double sigma,
1604  const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1605 {
1606  laplacianOfGaussianMultiArray( source.first, source.second, source.third,
1607  dest.first, dest.second, sigma, opt );
1608 }
1609 
1610 /********************************************************/
1611 /* */
1612 /* hessianOfGaussianMultiArray */
1613 /* */
1614 /********************************************************/
1615 
1616 /** \brief Calculate Hessian matrix of a N-dimensional arrays using Gaussian derivative filters.
1617 
1618  This function computes the Hessian matrix the given scalar N-dimensional
1619  array with a sequence of second-derivative-of-Gaussian filters at the given
1620  standard deviation <tt>sigma</tt>. Both source and destination arrays
1621  are represented by iterators, shape objects and accessors. The destination array must
1622  have a vector valued element type with N*(N+1)/2 elements (it represents the
1623  upper triangular part of the symmetric Hessian matrix). This function is implemented by calls to
1624  \ref separableConvolveMultiArray() with the appropriate kernels.
1625 
1626  Anisotropic data should be passed with appropriate
1627  \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional
1628  unless the parameter <tt>sigma</tt> is left out.
1629 
1630  <b> Declarations:</b>
1631 
1632  pass arguments explicitly:
1633  \code
1634  namespace vigra {
1635  template <class SrcIterator, class SrcShape, class SrcAccessor,
1636  class DestIterator, class DestAccessor>
1637  void
1638  hessianOfGaussianMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1639  DestIterator diter, DestAccessor dest,
1640  double sigma, const ConvolutionOptions<N> & opt);
1641  }
1642  \endcode
1643 
1644  use argument objects in conjunction with \ref ArgumentObjectFactories :
1645  \code
1646  namespace vigra {
1647  template <class SrcIterator, class SrcShape, class SrcAccessor,
1648  class DestIterator, class DestAccessor>
1649  void
1650  hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1651  pair<DestIterator, DestAccessor> const & dest,
1652  double sigma, const ConvolutionOptions<N> & opt);
1653  }
1654  \endcode
1655 
1656  <b> Usage:</b>
1657 
1658  <b>\#include</b> <vigra/multi_convolution.hxx>
1659 
1660  \code
1661  MultiArray<3, float> source(shape);
1662  MultiArray<3, TinyVector<float, 6> > dest(shape);
1663  ...
1664  // compute Hessian at scale sigma
1665  hessianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma);
1666  \endcode
1667 
1668  <b> Usage with anisotropic data:</b>
1669 
1670  <b>\#include</b> <vigra/multi_convolution.hxx>
1671 
1672  \code
1673  MultiArray<3, float> source(shape);
1674  MultiArray<3, TinyVector<float, 6> > dest(shape);
1675  TinyVector<float, 3> step_size;
1676  TinyVector<float, 3> resolution_sigmas;
1677  ...
1678  // compute Hessian at scale sigma
1679  hessianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma,
1680  ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
1681  \endcode
1682 
1683  <b> Required Interface:</b>
1684 
1685  see \ref separableConvolveMultiArray(), in addition:
1686 
1687  \code
1688  int dimension = 0;
1689  VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
1690  \endcode
1691 
1692  \see separableConvolveMultiArray(), vectorToTensorMultiArray()
1693 */
1694 doxygen_overloaded_function(template <...> void hessianOfGaussianMultiArray)
1695 
1696 template <class SrcIterator, class SrcShape, class SrcAccessor,
1697  class DestIterator, class DestAccessor>
1698 void
1699 hessianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1700  DestIterator di, DestAccessor dest,
1701  ConvolutionOptions<SrcShape::static_size> const & opt )
1702 {
1703  typedef typename DestAccessor::value_type DestType;
1704  typedef typename DestType::value_type DestValueType;
1705  typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1706 
1707  static const int N = SrcShape::static_size;
1708  static const int M = N*(N+1)/2;
1709  typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
1710 
1711  for(int k=0; k<N; ++k)
1712  if(shape[k] <=0)
1713  return;
1714 
1715  vigra_precondition(M == (int)dest.size(di),
1716  "hessianOfGaussianMultiArray(): Wrong number of channels in output array.");
1717 
1718  ParamType params_init = opt.scaleParams();
1719 
1720  ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
1721  ParamType params(params_init);
1722  for (int dim = 0; dim < N; ++dim, ++params)
1723  {
1724  double sigma = params.sigma_scaled("hessianOfGaussianMultiArray");
1725  plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
1726  }
1727 
1728  typedef VectorElementAccessor<DestAccessor> ElementAccessor;
1729 
1730  // compute elements of the Hessian matrix
1731  ParamType params_i(params_init);
1732  for (int b=0, i=0; i<N; ++i, ++params_i)
1733  {
1734  ParamType params_j(params_i);
1735  for (int j=i; j<N; ++j, ++b, ++params_j)
1736  {
1737  ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
1738  if(i == j)
1739  {
1740  kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 2, 1.0, opt.window_ratio);
1741  }
1742  else
1743  {
1744  kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 1, 1.0, opt.window_ratio);
1745  kernels[j].initGaussianDerivative(params_j.sigma_scaled(), 1, 1.0, opt.window_ratio);
1746  }
1747  detail::scaleKernel(kernels[i], 1 / params_i.step_size());
1748  detail::scaleKernel(kernels[j], 1 / params_j.step_size());
1749  separableConvolveMultiArray(si, shape, src, di, ElementAccessor(b, dest),
1750  kernels.begin(), opt.from_point, opt.to_point);
1751  }
1752  }
1753 }
1754 
1755 template <class SrcIterator, class SrcShape, class SrcAccessor,
1756  class DestIterator, class DestAccessor>
1757 inline void
1758 hessianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1759  DestIterator di, DestAccessor dest, double sigma,
1760  const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1761 {
1762  ConvolutionOptions<SrcShape::static_size> par = opt;
1763  hessianOfGaussianMultiArray(si, shape, src, di, dest, par.stdDev(sigma));
1764 }
1765 
1766 template <class SrcIterator, class SrcShape, class SrcAccessor,
1767  class DestIterator, class DestAccessor>
1768 inline void
1769 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1770  pair<DestIterator, DestAccessor> const & dest,
1771  ConvolutionOptions<SrcShape::static_size> const & opt )
1772 {
1773  hessianOfGaussianMultiArray( source.first, source.second, source.third,
1774  dest.first, dest.second, opt );
1775 }
1776 
1777 template <class SrcIterator, class SrcShape, class SrcAccessor,
1778  class DestIterator, class DestAccessor>
1779 inline void
1780 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1781  pair<DestIterator, DestAccessor> const & dest,
1782  double sigma,
1783  const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1784 {
1785  hessianOfGaussianMultiArray( source.first, source.second, source.third,
1786  dest.first, dest.second, sigma, opt );
1787 }
1788 
1789 namespace detail {
1790 
1791 template<int N, class VectorType>
1792 struct StructurTensorFunctor
1793 {
1794  typedef VectorType result_type;
1795  typedef typename VectorType::value_type ValueType;
1796 
1797  template <class T>
1798  VectorType operator()(T const & in) const
1799  {
1800  VectorType res;
1801  for(int b=0, i=0; i<N; ++i)
1802  {
1803  for(int j=i; j<N; ++j, ++b)
1804  {
1805  res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
1806  }
1807  }
1808  return res;
1809  }
1810 };
1811 
1812 } // namespace detail
1813 
1814 /********************************************************/
1815 /* */
1816 /* structureTensorMultiArray */
1817 /* */
1818 /********************************************************/
1819 
1820 /** \brief Calculate th structure tensor of a multi-dimensional arrays.
1821 
1822  This function computes the gradient (outer product) tensor for each element
1823  of the given N-dimensional array with first-derivative-of-Gaussian filters at
1824  the given <tt>innerScale</tt>, followed by Gaussian smoothing at <tt>outerScale</tt>.
1825  Both source and destination arrays are represented by iterators, shape objects and
1826  accessors. The destination array must have a vector valued pixel type with
1827  N*(N+1)/2 elements (it represents the upper triangular part of the symmetric
1828  structure tensor matrix). If the source array is also vector valued, the
1829  resulting structure tensor is the sum of the individual tensors for each channel.
1830  This function is implemented by calls to
1831  \ref separableConvolveMultiArray() with the appropriate kernels.
1832 
1833  Anisotropic data should be passed with appropriate
1834  \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional
1835  unless the parameters <tt>innerScale</tt> and <tt>outerScale</tt> are
1836  both left out.
1837 
1838  <b> Declarations:</b>
1839 
1840  pass arguments explicitly:
1841  \code
1842  namespace vigra {
1843  template <class SrcIterator, class SrcShape, class SrcAccessor,
1844  class DestIterator, class DestAccessor>
1845  void
1846  structureTensorMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1847  DestIterator diter, DestAccessor dest,
1848  double innerScale, double outerScale,
1849  const ConvolutionOptions<N> & opt);
1850  }
1851  \endcode
1852 
1853  use argument objects in conjunction with \ref ArgumentObjectFactories :
1854  \code
1855  namespace vigra {
1856  template <class SrcIterator, class SrcShape, class SrcAccessor,
1857  class DestIterator, class DestAccessor>
1858  void
1859  structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1860  pair<DestIterator, DestAccessor> const & dest,
1861  double innerScale, double outerScale,
1862  const ConvolutionOptions<N> & opt);
1863  }
1864  \endcode
1865 
1866  <b> Usage:</b>
1867 
1868  <b>\#include</b> <vigra/multi_convolution.hxx>
1869 
1870  \code
1871  MultiArray<3, RGBValue<float> > source(shape);
1872  MultiArray<3, TinyVector<float, 6> > dest(shape);
1873  ...
1874  // compute structure tensor at scales innerScale and outerScale
1875  structureTensorMultiArray(srcMultiArrayRange(source), destMultiArray(dest), innerScale, outerScale);
1876  \endcode
1877 
1878  <b> Usage with anisotropic data:</b>
1879 
1880  <b>\#include</b> <vigra/multi_convolution.hxx>
1881 
1882  \code
1883  MultiArray<3, RGBValue<float> > source(shape);
1884  MultiArray<3, TinyVector<float, 6> > dest(shape);
1885  TinyVector<float, 3> step_size;
1886  TinyVector<float, 3> resolution_sigmas;
1887  ...
1888  // compute structure tensor at scales innerScale and outerScale
1889  structureTensorMultiArray(srcMultiArrayRange(source), destMultiArray(dest), innerScale, outerScale,
1890  ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
1891  \endcode
1892 
1893  <b> Required Interface:</b>
1894 
1895  see \ref separableConvolveMultiArray(), in addition:
1896 
1897  \code
1898  int dimension = 0;
1899  VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
1900  \endcode
1901 
1902  \see separableConvolveMultiArray(), vectorToTensorMultiArray()
1903 */
1904 doxygen_overloaded_function(template <...> void structureTensorMultiArray)
1905 
1906 template <class SrcIterator, class SrcShape, class SrcAccessor,
1907  class DestIterator, class DestAccessor>
1908 void
1909 structureTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1910  DestIterator di, DestAccessor dest,
1911  ConvolutionOptions<SrcShape::static_size> const & opt)
1912 {
1913  static const int N = SrcShape::static_size;
1914  static const int M = N*(N+1)/2;
1915 
1916  typedef typename DestAccessor::value_type DestType;
1917  typedef typename DestType::value_type DestValueType;
1918  typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1919  typedef TinyVector<KernelType, N> GradientVector;
1920  typedef typename AccessorTraits<GradientVector>::default_accessor GradientAccessor;
1921  typedef typename AccessorTraits<DestType>::default_accessor GradientTensorAccessor;
1922 
1923  for(int k=0; k<N; ++k)
1924  if(shape[k] <=0)
1925  return;
1926 
1927  vigra_precondition(M == (int)dest.size(di),
1928  "structureTensorMultiArray(): Wrong number of channels in output array.");
1929 
1930  ConvolutionOptions<N> innerOptions = opt;
1931  ConvolutionOptions<N> outerOptions = opt.outerOptions();
1932  typename ConvolutionOptions<N>::ScaleIterator params = outerOptions.scaleParams();
1933 
1934  SrcShape gradientShape(shape);
1935  if(opt.to_point != SrcShape())
1936  {
1937  for(int k=0; k<N; ++k, ++params)
1938  {
1939  Kernel1D<double> gauss;
1940  gauss.initGaussian(params.sigma_scaled("structureTensorMultiArray"), 1.0, opt.window_ratio);
1941  int dilation = gauss.right();
1942  innerOptions.from_point[k] = std::max<MultiArrayIndex>(0, opt.from_point[k] - dilation);
1943  innerOptions.to_point[k] = std::min<MultiArrayIndex>(shape[k], opt.to_point[k] + dilation);
1944  }
1945  outerOptions.from_point -= innerOptions.from_point;
1946  outerOptions.to_point -= innerOptions.from_point;
1947  gradientShape = innerOptions.to_point - innerOptions.from_point;
1948  }
1949 
1950  MultiArray<N, GradientVector> gradient(gradientShape);
1951  MultiArray<N, DestType> gradientTensor(gradientShape);
1952  gaussianGradientMultiArray(si, shape, src,
1953  gradient.traverser_begin(), GradientAccessor(),
1954  innerOptions,
1955  "structureTensorMultiArray");
1956 
1957  transformMultiArray(gradient.traverser_begin(), gradientShape, GradientAccessor(),
1958  gradientTensor.traverser_begin(), GradientTensorAccessor(),
1959  detail::StructurTensorFunctor<N, DestType>());
1960 
1961  gaussianSmoothMultiArray(gradientTensor.traverser_begin(), gradientShape, GradientTensorAccessor(),
1962  di, dest, outerOptions,
1963  "structureTensorMultiArray");
1964 }
1965 
1966 template <class SrcIterator, class SrcShape, class SrcAccessor,
1967  class DestIterator, class DestAccessor>
1968 inline void
1969 structureTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1970  DestIterator di, DestAccessor dest,
1971  double innerScale, double outerScale,
1972  const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1973 {
1974  ConvolutionOptions<SrcShape::static_size> par = opt;
1975  structureTensorMultiArray(si, shape, src, di, dest,
1976  par.stdDev(innerScale).outerScale(outerScale));
1977 }
1978 
1979 template <class SrcIterator, class SrcShape, class SrcAccessor,
1980  class DestIterator, class DestAccessor>
1981 inline void
1982 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1983  pair<DestIterator, DestAccessor> const & dest,
1984  ConvolutionOptions<SrcShape::static_size> const & opt )
1985 {
1986  structureTensorMultiArray( source.first, source.second, source.third,
1987  dest.first, dest.second, opt );
1988 }
1989 
1990 
1991 template <class SrcIterator, class SrcShape, class SrcAccessor,
1992  class DestIterator, class DestAccessor>
1993 inline void
1994 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1995  pair<DestIterator, DestAccessor> const & dest,
1996  double innerScale, double outerScale,
1997  const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1998 {
1999  structureTensorMultiArray( source.first, source.second, source.third,
2000  dest.first, dest.second,
2001  innerScale, outerScale, opt);
2002 }
2003 
2004 //@}
2005 
2006 } //-- namespace vigra
2007 
2008 
2009 #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.9.0 (Wed Feb 27 2013)