[ 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 
51 namespace vigra
52 {
53 
54 
55 namespace detail
56 {
57 
58 /********************************************************/
59 /* */
60 /* internalSeparableConvolveMultiArray */
61 /* */
62 /********************************************************/
63 
64 template <class SrcIterator, class SrcShape, class SrcAccessor,
65  class DestIterator, class DestAccessor, class KernelIterator>
66 void
67 internalSeparableConvolveMultiArrayTmp(
68  SrcIterator si, SrcShape const & shape, SrcAccessor src,
69  DestIterator di, DestAccessor dest, KernelIterator kit)
70 {
71  enum { N = 1 + SrcIterator::level };
72 
73  typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
74 
75  // temporay array to hold the current line to enable in-place operation
76  ArrayVector<TmpType> tmp( shape[0] );
77 
78  typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
79  typedef MultiArrayNavigator<DestIterator, N> DNavigator;
80 
81  { // only operate on first dimension here
82  SNavigator snav( si, shape, 0 );
83  DNavigator dnav( di, shape, 0 );
84 
85  for( ; snav.hasMore(); snav++, dnav++ )
86  {
87  // first copy source to temp for maximum cache efficiency
88  copyLine( snav.begin(), snav.end(), src,
89  tmp.begin(), typename AccessorTraits<TmpType>::default_accessor() );
90 
91  convolveLine( srcIterRange(tmp.begin(), tmp.end(),
92  typename AccessorTraits<TmpType>::default_const_accessor()),
93  destIter( dnav.begin(), dest ),
94  kernel1d( *kit ) );
95  }
96  ++kit;
97  }
98 
99  // operate on further dimensions
100  for( int d = 1; d < N; ++d, ++kit )
101  {
102  DNavigator dnav( di, shape, d );
103 
104  tmp.resize( shape[d] );
105 
106  for( ; dnav.hasMore(); dnav++ )
107  {
108  // first copy source to temp for maximum cache efficiency
109  copyLine( dnav.begin(), dnav.end(), dest,
110  tmp.begin(), typename AccessorTraits<TmpType>::default_accessor() );
111 
112  convolveLine( srcIterRange(tmp.begin(), tmp.end(),
113  typename AccessorTraits<TmpType>::default_const_accessor()),
114  destIter( dnav.begin(), dest ),
115  kernel1d( *kit ) );
116  }
117  }
118 }
119 
120 
121 } // namespace detail
122 
123 /** \addtogroup MultiArrayConvolutionFilters Convolution filters for multi-dimensional arrays.
124 
125  These functions realize a separable convolution on an arbitrary dimensional
126  array that is specified by iterators (compatible to \ref MultiIteratorPage)
127  and shape objects. It can therefore be applied to a wide range of data structures
128  (\ref vigra::MultiArrayView, \ref vigra::MultiArray etc.).
129 */
130 //@{
131 
132 /********************************************************/
133 /* */
134 /* separableConvolveMultiArray */
135 /* */
136 /********************************************************/
137 
138 /** \brief Separated convolution on multi-dimensional arrays.
139 
140  This function computes a separated convolution on all dimensions
141  of the given multi-dimensional array. Both source and destination
142  arrays are represented by iterators, shape objects and accessors.
143  The destination array is required to already have the correct size.
144 
145  There are two variants of this functions: one takes a single kernel
146  of type \ref vigra::Kernel1D which is then applied to all dimensions,
147  whereas the other requires an iterator referencing a sequence of
148  \ref vigra::Kernel1D objects, one for every dimension of the data.
149  Then the first kernel in this sequence is applied to the innermost
150  dimension (e.g. the x-dimension of an image), while the last is applied to the
151  outermost dimension (e.g. the z-dimension in a 3D image).
152 
153  This function may work in-place, which means that <tt>siter == diter</tt> is allowed.
154  A full-sized internal array is only allocated if working on the destination
155  array directly would cause round-off errors (i.e. if
156  <tt>typeid(typename NumericTraits<typename DestAccessor::value_type>::RealPromote)
157  != typeid(typename DestAccessor::value_type)</tt>.
158 
159  <b> Declarations:</b>
160 
161  pass arguments explicitly:
162  \code
163  namespace vigra {
164  // apply the same kernel to all dimensions
165  template <class SrcIterator, class SrcShape, class SrcAccessor,
166  class DestIterator, class DestAccessor, class T>
167  void
168  separableConvolveMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
169  DestIterator diter, DestAccessor dest,
170  Kernel1D<T> const & kernel);
171 
172  // apply each kernel from the sequence 'kernels' in turn
173  template <class SrcIterator, class SrcShape, class SrcAccessor,
174  class DestIterator, class DestAccessor, class KernelIterator>
175  void
176  separableConvolveMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
177  DestIterator diter, DestAccessor dest,
178  KernelIterator kernels);
179  }
180  \endcode
181 
182  use argument objects in conjunction with \ref ArgumentObjectFactories :
183  \code
184  namespace vigra {
185  // apply the same kernel to all dimensions
186  template <class SrcIterator, class SrcShape, class SrcAccessor,
187  class DestIterator, class DestAccessor, class T>
188  void
189  separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
190  pair<DestIterator, DestAccessor> const & dest,
191  Kernel1D<T> const & kernel);
192 
193  // apply each kernel from the sequence 'kernels' in turn
194  template <class SrcIterator, class SrcShape, class SrcAccessor,
195  class DestIterator, class DestAccessor, class KernelIterator>
196  void
197  separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
198  pair<DestIterator, DestAccessor> const & dest,
199  KernelIterator kernels);
200  }
201  \endcode
202 
203  <b> Usage:</b>
204 
205  <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
206 
207  \code
208  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
209  MultiArray<3, unsigned char> source(shape);
210  MultiArray<3, float> dest(shape);
211  ...
212  Kernel1D<float> gauss;
213  gauss.initGaussian(sigma);
214  // create 3 Gauss kernels, one for each dimension
215  ArrayVector<Kernel1D<float> > kernels(3, gauss);
216 
217  // perform Gaussian smoothing on all dimensions
218  separableConvolveMultiArray(srcMultiArrayRange(source), destMultiArray(dest),
219  kernels.begin());
220  \endcode
221 
222  \see vigra::Kernel1D, convolveLine()
223 */
225 
226 template <class SrcIterator, class SrcShape, class SrcAccessor,
227  class DestIterator, class DestAccessor, class KernelIterator>
228 void
229 separableConvolveMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
230  DestIterator d, DestAccessor dest, KernelIterator kernels )
231 {
232  typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
233 
234  if(!IsSameType<TmpType, typename DestAccessor::value_type>::boolResult)
235  {
236  // need a temporary array to avoid rounding errors
237  MultiArray<SrcShape::static_size, TmpType> tmpArray(shape);
238  detail::internalSeparableConvolveMultiArrayTmp( s, shape, src,
239  tmpArray.traverser_begin(), typename AccessorTraits<TmpType>::default_accessor(), kernels );
240  copyMultiArray(srcMultiArrayRange(tmpArray), destIter(d, dest));
241  }
242  else
243  {
244  // work directly on the destination array
245  detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, d, dest, kernels );
246  }
247 }
248 
249 template <class SrcIterator, class SrcShape, class SrcAccessor,
250  class DestIterator, class DestAccessor, class KernelIterator>
251 inline
253  triple<SrcIterator, SrcShape, SrcAccessor> const & source,
254  pair<DestIterator, DestAccessor> const & dest, KernelIterator kit )
255 {
256  separableConvolveMultiArray( source.first, source.second, source.third,
257  dest.first, dest.second, kit );
258 }
259 
260 template <class SrcIterator, class SrcShape, class SrcAccessor,
261  class DestIterator, class DestAccessor, class T>
262 inline void
263 separableConvolveMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
264  DestIterator d, DestAccessor dest,
265  Kernel1D<T> const & kernel )
266 {
267  ArrayVector<Kernel1D<T> > kernels(shape.size(), kernel);
268 
269  separableConvolveMultiArray( s, shape, src, d, dest, kernels.begin() );
270 }
271 
272 template <class SrcIterator, class SrcShape, class SrcAccessor,
273  class DestIterator, class DestAccessor, class T>
274 inline void
275 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
276  pair<DestIterator, DestAccessor> const & dest,
277  Kernel1D<T> const & kernel )
278 {
279  ArrayVector<Kernel1D<T> > kernels(source.second.size(), kernel);
280 
281  separableConvolveMultiArray( source.first, source.second, source.third,
282  dest.first, dest.second, kernels.begin() );
283 }
284 
285 /********************************************************/
286 /* */
287 /* convolveMultiArrayOneDimension */
288 /* */
289 /********************************************************/
290 
291 /** \brief Convolution along a single dimension of a multi-dimensional arrays.
292 
293  This function computes a convolution along one dimension (specified by
294  the parameter <tt>dim</tt> of the given multi-dimensional array with the given
295  <tt>kernel</tt>. Both source and destination arrays are represented by
296  iterators, shape objects and accessors. The destination array is required to
297  already have the correct size.
298 
299  This function may work in-place, which means that <tt>siter == diter</tt> is allowed.
300 
301  <b> Declarations:</b>
302 
303  pass arguments explicitly:
304  \code
305  namespace vigra {
306  template <class SrcIterator, class SrcShape, class SrcAccessor,
307  class DestIterator, class DestAccessor, class T>
308  void
309  convolveMultiArrayOneDimension(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
310  DestIterator diter, DestAccessor dest,
311  unsigned int dim, vigra::Kernel1D<T> const & kernel);
312  }
313  \endcode
314 
315  use argument objects in conjunction with \ref ArgumentObjectFactories :
316  \code
317  namespace vigra {
318  template <class SrcIterator, class SrcShape, class SrcAccessor,
319  class DestIterator, class DestAccessor, class T>
320  void
321  convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
322  pair<DestIterator, DestAccessor> const & dest,
323  unsigned int dim, vigra::Kernel1D<T> const & kernel);
324  }
325  \endcode
326 
327  <b> Usage:</b>
328 
329  <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
330 
331  \code
332  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
333  MultiArray<3, unsigned char> source(shape);
334  MultiArray<3, float> dest(shape);
335  ...
336  Kernel1D<float> gauss;
337  gauss.initGaussian(sigma);
338 
339  // perform Gaussian smoothing along dimensions 1 (height)
340  convolveMultiArrayOneDimension(srcMultiArrayRange(source), destMultiArray(dest), 1, gauss);
341  \endcode
342 
343  \see separableConvolveMultiArray()
344 */
346 
347 template <class SrcIterator, class SrcShape, class SrcAccessor,
348  class DestIterator, class DestAccessor, class T>
349 void
350 convolveMultiArrayOneDimension(SrcIterator s, SrcShape const & shape, SrcAccessor src,
351  DestIterator d, DestAccessor dest,
352  unsigned int dim, vigra::Kernel1D<T> const & kernel )
353 {
354  enum { N = 1 + SrcIterator::level };
355  vigra_precondition( dim < N,
356  "convolveMultiArrayOneDimension(): The dimension number to convolve must be smaller "
357  "than the data dimensionality" );
358 
359  typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
360  ArrayVector<TmpType> tmp( shape[dim] );
361 
362  typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
363  typedef MultiArrayNavigator<DestIterator, N> DNavigator;
364 
365  SNavigator snav( s, shape, dim );
366  DNavigator dnav( d, shape, dim );
367 
368  for( ; snav.hasMore(); snav++, dnav++ )
369  {
370  // first copy source to temp for maximum cache efficiency
371  copyLine( snav.begin(), snav.end(), src,
372  tmp.begin(), typename AccessorTraits<TmpType>::default_accessor() );
373 
374  convolveLine( srcIterRange( tmp.begin(), tmp.end(), typename AccessorTraits<TmpType>::default_const_accessor()),
375  destIter( dnav.begin(), dest ),
376  kernel1d( kernel ) );
377  }
378 }
379 
380 template <class SrcIterator, class SrcShape, class SrcAccessor,
381  class DestIterator, class DestAccessor, class T>
382 inline void
383 convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
384  pair<DestIterator, DestAccessor> const & dest,
385  unsigned int dim, vigra::Kernel1D<T> const & kernel )
386 {
387  convolveMultiArrayOneDimension( source.first, source.second, source.third,
388  dest.first, dest.second, dim, kernel );
389 }
390 
391 /********************************************************/
392 /* */
393 /* gaussianSmoothMultiArray */
394 /* */
395 /********************************************************/
396 
397 /** \brief Isotropic Gaussian smoothing of a multi-dimensional arrays.
398 
399  This function computes an isotropic convolution of the given multi-dimensional
400  array with a Gaussian filter at the given standard deviation <tt>sigma</tt>.
401  Both source and destination arrays are represented by
402  iterators, shape objects and accessors. The destination array is required to
403  already have the correct size. This function may work in-place, which means
404  that <tt>siter == diter</tt> is allowed. It is implemented by a call to
405  \ref separableConvolveMultiArray() with the appropriate kernel.
406  If the data are anisotropic (different pixel size along different dimensions)
407  you should call \ref separableConvolveMultiArray() directly with the appropriate
408  anisotropic Gaussians.
409 
410  <b> Declarations:</b>
411 
412  pass arguments explicitly:
413  \code
414  namespace vigra {
415  template <class SrcIterator, class SrcShape, class SrcAccessor,
416  class DestIterator, class DestAccessor>
417  void
418  gaussianSmoothMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
419  DestIterator diter, DestAccessor dest,
420  double sigma);
421  }
422  \endcode
423 
424  use argument objects in conjunction with \ref ArgumentObjectFactories :
425  \code
426  namespace vigra {
427  template <class SrcIterator, class SrcShape, class SrcAccessor,
428  class DestIterator, class DestAccessor>
429  void
430  gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
431  pair<DestIterator, DestAccessor> const & dest,
432  double sigma);
433  }
434  \endcode
435 
436  <b> Usage:</b>
437 
438  <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
439 
440  \code
441  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
442  MultiArray<3, unsigned char> source(shape);
443  MultiArray<3, float> dest(shape);
444  ...
445  // perform isotropic Gaussian smoothing at scale 'sigma'
446  gaussianSmoothMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma);
447  \endcode
448 
449  \see separableConvolveMultiArray()
450 */
452 
453 template <class SrcIterator, class SrcShape, class SrcAccessor,
454  class DestIterator, class DestAccessor>
455 void
456 gaussianSmoothMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
457  DestIterator d, DestAccessor dest, double sigma )
458 {
459  Kernel1D<double> gauss;
460  gauss.initGaussian( sigma );
461 
462  separableConvolveMultiArray( s, shape, src, d, dest, gauss);
463 }
464 
465 template <class SrcIterator, class SrcShape, class SrcAccessor,
466  class DestIterator, class DestAccessor>
467 inline void
468 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
469  pair<DestIterator, DestAccessor> const & dest,
470  double sigma )
471 {
472  gaussianSmoothMultiArray( source.first, source.second, source.third,
473  dest.first, dest.second, sigma );
474 }
475 
476 /********************************************************/
477 /* */
478 /* gaussianGradientMultiArray */
479 /* */
480 /********************************************************/
481 
482 /** \brief Calculate Gaussian gradient of a multi-dimensional arrays.
483 
484  This function computes the Gaussian gradient of the given multi-dimensional
485  array with a sequence of first-derivative-of-Gaussian filters at the given
486  standard deviation <tt>sigma</tt> (differentiation is applied to each dimension
487  in turn, starting with the innermost dimension). Both source and destination arrays
488  are represented by iterators, shape objects and accessors. The destination array is
489  required to have a vector valued pixel type with as many elements as the number of
490  dimensions. This function is implemented by calls to
491  \ref separableConvolveMultiArray() with the appropriate kernels.
492  If the data are anisotropic (different pixel size along different dimensions)
493  you should call \ref separableConvolveMultiArray() directly with the appropriate
494  anisotropic Gaussian derivatives.
495 
496  <b> Declarations:</b>
497 
498  pass arguments explicitly:
499  \code
500  namespace vigra {
501  template <class SrcIterator, class SrcShape, class SrcAccessor,
502  class DestIterator, class DestAccessor>
503  void
504  gaussianGradientMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
505  DestIterator diter, DestAccessor dest,
506  double sigma);
507  }
508  \endcode
509 
510  use argument objects in conjunction with \ref ArgumentObjectFactories :
511  \code
512  namespace vigra {
513  template <class SrcIterator, class SrcShape, class SrcAccessor,
514  class DestIterator, class DestAccessor>
515  void
516  gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
517  pair<DestIterator, DestAccessor> const & dest,
518  double sigma);
519  }
520  \endcode
521 
522  <b> Usage:</b>
523 
524  <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
525 
526  \code
527  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
528  MultiArray<3, unsigned char> source(shape);
529  MultiArray<3, TinyVector<float, 3> > dest(shape);
530  ...
531  // compute Gaussian gradient at scale sigma
532  gaussianGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma);
533  \endcode
534 
535  <b> Required Interface:</b>
536 
537  see \ref separableConvolveMultiArray(), in addition:
538 
539  \code
540  int dimension = 0;
541  VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
542  \endcode
543 
544  \see separableConvolveMultiArray()
545 */
547 
548 template <class SrcIterator, class SrcShape, class SrcAccessor,
549  class DestIterator, class DestAccessor>
550 void
551 gaussianGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
552  DestIterator di, DestAccessor dest, double sigma )
553 {
554  typedef typename DestAccessor::value_type DestType;
555  typedef typename DestType::value_type DestValueType;
556  typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
557 
558  static const int N = SrcShape::static_size;
559 
560  for(int k=0; k<N; ++k)
561  if(shape[k] <=0)
562  return;
563 
564  vigra_precondition(N == dest.size(di),
565  "gaussianGradientMultiArray(): Wrong number of channels in output array.");
566 
567  vigra_precondition(sigma > 0.0, "gaussianGradientMultiArray(): Scale must be positive.");
568 
569  Kernel1D<KernelType> gauss, derivative;
570  gauss.initGaussian(sigma);
571 
572  typedef VectorElementAccessor<DestAccessor> ElementAccessor;
573 
574  // compute gradient components
575  for(int d = 0; d < N; ++d )
576  {
577  ArrayVector<Kernel1D<KernelType> > kernels(N, gauss);
578  kernels[d].initGaussianDerivative(sigma, 1);
579  separableConvolveMultiArray( si, shape, src, di, ElementAccessor(d, dest), kernels.begin());
580  }
581 }
582 
583 template <class SrcIterator, class SrcShape, class SrcAccessor,
584  class DestIterator, class DestAccessor>
585 inline void
586 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
587  pair<DestIterator, DestAccessor> const & dest, double sigma )
588 {
589  gaussianGradientMultiArray( source.first, source.second, source.third,
590  dest.first, dest.second, sigma );
591 }
592 
593 /********************************************************/
594 /* */
595 /* symmetricGradientMultiArray */
596 /* */
597 /********************************************************/
598 
599 /** \brief Calculate gradient of a multi-dimensional arrays using symmetric difference filters.
600 
601  This function computes the gradient of the given multi-dimensional
602  array with a sequence of symmetric difference filters a (differentiation is applied
603  to each dimension in turn, starting with the innermost dimension). Both source and
604  destination arrays are represented by iterators, shape objects and accessors.
605  The destination array is required to have a vector valued pixel type with as many
606  elements as the number of dimensions. This function is implemented by calls to
607  \ref convolveMultiArrayOneDimension() with the symmetric difference kernel.
608 
609  <b> Declarations:</b>
610 
611  pass arguments explicitly:
612  \code
613  namespace vigra {
614  template <class SrcIterator, class SrcShape, class SrcAccessor,
615  class DestIterator, class DestAccessor>
616  void
617  symmetricGradientMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
618  DestIterator diter, DestAccessor dest);
619  }
620  \endcode
621 
622  use argument objects in conjunction with \ref ArgumentObjectFactories :
623  \code
624  namespace vigra {
625  template <class SrcIterator, class SrcShape, class SrcAccessor,
626  class DestIterator, class DestAccessor>
627  void
628  symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
629  pair<DestIterator, DestAccessor> const & dest);
630  }
631  \endcode
632 
633  <b> Usage:</b>
634 
635  <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
636 
637  \code
638  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
639  MultiArray<3, unsigned char> source(shape);
640  MultiArray<3, TinyVector<float, 3> > dest(shape);
641  ...
642  // compute gradient
643  symmetricGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest));
644  \endcode
645 
646  <b> Required Interface:</b>
647 
648  see \ref convolveMultiArrayOneDimension(), in addition:
649 
650  \code
651  int dimension = 0;
652  VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
653  \endcode
654 
655  \see convolveMultiArrayOneDimension()
656 */
658 
659 template <class SrcIterator, class SrcShape, class SrcAccessor,
660  class DestIterator, class DestAccessor>
661 void
662 symmetricGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
663  DestIterator di, DestAccessor dest)
664 {
665  typedef typename DestAccessor::value_type DestType;
666  typedef typename DestType::value_type DestValueType;
667  typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
668 
669  static const int N = SrcShape::static_size;
670 
671  for(int k=0; k<N; ++k)
672  if(shape[k] <=0)
673  return;
674 
675  vigra_precondition(N == dest.size(di),
676  "symmetricGradientMultiArray(): Wrong number of channels in output array.");
677 
678  Kernel1D<KernelType> filter;
679  filter.initSymmetricGradient();
680 
681  typedef VectorElementAccessor<DestAccessor> ElementAccessor;
682 
683  // compute gradient components
684  for(int d = 0; d < N; ++d )
685  {
686  convolveMultiArrayOneDimension(si, shape, src,
687  di, ElementAccessor(d, dest),
688  d, filter);
689  }
690 }
691 
692 template <class SrcIterator, class SrcShape, class SrcAccessor,
693  class DestIterator, class DestAccessor>
694 inline void
695 symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
696  pair<DestIterator, DestAccessor> const & dest )
697 {
698  symmetricGradientMultiArray(source.first, source.second, source.third,
699  dest.first, dest.second);
700 }
701 
702 
703 /********************************************************/
704 /* */
705 /* laplacianOfGaussianMultiArray */
706 /* */
707 /********************************************************/
708 
709 /** \brief Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters.
710 
711  This function computes the Laplacian the given N-dimensional
712  array with a sequence of second-derivative-of-Gaussian filters at the given
713  standard deviation <tt>sigma</tt>. Both source and destination arrays
714  are represented by iterators, shape objects and accessors. Both source and destination
715  arrays must have scalar value_type. This function is implemented by calls to
716  \ref separableConvolveMultiArray() with the appropriate kernels, followed by summation.
717 
718  <b> Declarations:</b>
719 
720  pass arguments explicitly:
721  \code
722  namespace vigra {
723  template <class SrcIterator, class SrcShape, class SrcAccessor,
724  class DestIterator, class DestAccessor>
725  void
726  laplacianOfGaussianMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
727  DestIterator diter, DestAccessor dest,
728  double sigma);
729  }
730  \endcode
731 
732  use argument objects in conjunction with \ref ArgumentObjectFactories :
733  \code
734  namespace vigra {
735  template <class SrcIterator, class SrcShape, class SrcAccessor,
736  class DestIterator, class DestAccessor>
737  void
738  laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
739  pair<DestIterator, DestAccessor> const & dest,
740  double sigma);
741  }
742  \endcode
743 
744  <b> Usage:</b>
745 
746  <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
747 
748  \code
749  MultiArray<3, float> source(shape);
750  MultiArray<3, float> laplacian(shape);
751  ...
752  // compute Laplacian at scale sigma
753  laplacianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(laplacian), sigma);
754  \endcode
755 
756  <b> Required Interface:</b>
757 
758  see \ref separableConvolveMultiArray(), in addition:
759 
760  \code
761  int dimension = 0;
762  VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
763  \endcode
764 
765  \see separableConvolveMultiArray()
766 */
768 
769 template <class SrcIterator, class SrcShape, class SrcAccessor,
770  class DestIterator, class DestAccessor>
771 void
772 laplacianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
773  DestIterator di, DestAccessor dest, double sigma )
774 {
775  using namespace functor;
776 
777  typedef typename DestAccessor::value_type DestType;
778  typedef typename NumericTraits<DestType>::RealPromote KernelType;
779  typedef typename AccessorTraits<KernelType>::default_accessor DerivativeAccessor;
780 
781  static const int N = SrcShape::static_size;
782 
783  vigra_precondition(sigma > 0.0, "laplacianOfGaussianMultiArray(): Scale must be positive.");
784 
785  Kernel1D<KernelType> gauss;
786  gauss.initGaussian(sigma);
787 
788  MultiArray<N, KernelType> derivative(shape);
789 
790  // compute 2nd derivatives and sum them up
791  for(int d = 0; d < N; ++d )
792  {
793  ArrayVector<Kernel1D<KernelType> > kernels(N, gauss);
794  kernels[d].initGaussianDerivative(sigma, 2);
795  if(d == 0)
796  {
797  separableConvolveMultiArray( si, shape, src,
798  di, dest, kernels.begin());
799  }
800  else
801  {
802  separableConvolveMultiArray( si, shape, src,
803  derivative.traverser_begin(), DerivativeAccessor(),
804  kernels.begin());
805  combineTwoMultiArrays(di, shape, dest, derivative.traverser_begin(), DerivativeAccessor(),
806  di, dest, Arg1() + Arg2() );
807  }
808  }
809 }
810 
811 template <class SrcIterator, class SrcShape, class SrcAccessor,
812  class DestIterator, class DestAccessor>
813 inline void
814 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
815  pair<DestIterator, DestAccessor> const & dest, double sigma )
816 {
817  laplacianOfGaussianMultiArray( source.first, source.second, source.third,
818  dest.first, dest.second, sigma );
819 }
820 
821 /********************************************************/
822 /* */
823 /* hessianOfGaussianMultiArray */
824 /* */
825 /********************************************************/
826 
827 /** \brief Calculate Hessian matrix of a N-dimensional arrays using Gaussian derivative filters.
828 
829  This function computes the Hessian matrix the given scalar N-dimensional
830  array with a sequence of second-derivative-of-Gaussian filters at the given
831  standard deviation <tt>sigma</tt>. Both source and destination arrays
832  are represented by iterators, shape objects and accessors. The destination array must
833  have a vector valued element type with N*(N+1)/2 elements (it represents the
834  upper triangular part of the symmetric Hessian matrix). This function is implemented by calls to
835  \ref separableConvolveMultiArray() with the appropriate kernels.
836 
837  <b> Declarations:</b>
838 
839  pass arguments explicitly:
840  \code
841  namespace vigra {
842  template <class SrcIterator, class SrcShape, class SrcAccessor,
843  class DestIterator, class DestAccessor>
844  void
845  hessianOfGaussianMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
846  DestIterator diter, DestAccessor dest,
847  double sigma);
848  }
849  \endcode
850 
851  use argument objects in conjunction with \ref ArgumentObjectFactories :
852  \code
853  namespace vigra {
854  template <class SrcIterator, class SrcShape, class SrcAccessor,
855  class DestIterator, class DestAccessor>
856  void
857  hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
858  pair<DestIterator, DestAccessor> const & dest,
859  double sigma);
860  }
861  \endcode
862 
863  <b> Usage:</b>
864 
865  <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
866 
867  \code
868  MultiArray<3, float> source(shape);
869  MultiArray<3, TinyVector<float, 6> > dest(shape);
870  ...
871  // compute Hessian at scale sigma
872  hessianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma);
873  \endcode
874 
875  <b> Required Interface:</b>
876 
877  see \ref separableConvolveMultiArray(), in addition:
878 
879  \code
880  int dimension = 0;
881  VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
882  \endcode
883 
884  \see separableConvolveMultiArray(), vectorToTensorMultiArray()
885 */
887 
888 template <class SrcIterator, class SrcShape, class SrcAccessor,
889  class DestIterator, class DestAccessor>
890 void
891 hessianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
892  DestIterator di, DestAccessor dest, double sigma )
893 {
894  typedef typename DestAccessor::value_type DestType;
895  typedef typename DestType::value_type DestValueType;
896  typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
897 
898  static const int N = SrcShape::static_size;
899  static const int M = N*(N+1)/2;
900 
901  for(int k=0; k<N; ++k)
902  if(shape[k] <=0)
903  return;
904 
905  vigra_precondition(M == dest.size(di),
906  "hessianOfGaussianMultiArray(): Wrong number of channels in output array.");
907 
908  vigra_precondition(sigma > 0.0, "hessianOfGaussianMultiArray(): Scale must be positive.");
909 
910  Kernel1D<KernelType> gauss;
911  gauss.initGaussian(sigma);
912 
913  typedef VectorElementAccessor<DestAccessor> ElementAccessor;
914 
915  // compute elements of the Hessian matrix
916  for(int b=0, i=0; i<N; ++i)
917  {
918  for(int j=i; j<N; ++j, ++b)
919  {
920  ArrayVector<Kernel1D<KernelType> > kernels(N, gauss);
921  if(i == j)
922  {
923  kernels[i].initGaussianDerivative(sigma, 2);
924  }
925  else
926  {
927  kernels[i].initGaussianDerivative(sigma, 1);
928  kernels[j].initGaussianDerivative(sigma, 1);
929  }
930  separableConvolveMultiArray(si, shape, src, di, ElementAccessor(b, dest),
931  kernels.begin());
932  }
933  }
934 }
935 
936 template <class SrcIterator, class SrcShape, class SrcAccessor,
937  class DestIterator, class DestAccessor>
938 inline void
939 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
940  pair<DestIterator, DestAccessor> const & dest, double sigma )
941 {
942  hessianOfGaussianMultiArray( source.first, source.second, source.third,
943  dest.first, dest.second, sigma );
944 }
945 
946 namespace detail {
947 
948 template<int N, class VectorType>
949 struct StructurTensorFunctor
950 {
951  typedef VectorType result_type;
952  typedef typename VectorType::value_type ValueType;
953 
954  template <class T>
955  VectorType operator()(T const & in) const
956  {
957  VectorType res;
958  for(int b=0, i=0; i<N; ++i)
959  {
960  for(int j=i; j<N; ++j, ++b)
961  {
962  res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
963  }
964  }
965  return res;
966  }
967 };
968 
969 } // namespace detail
970 
971 /********************************************************/
972 /* */
973 /* structureTensorMultiArray */
974 /* */
975 /********************************************************/
976 
977 /** \brief Calculate th structure tensor of a multi-dimensional arrays.
978 
979  This function computes the gradient (outer product) tensor for each element
980  of the given N-dimensional array with first-derivative-of-Gaussian filters at
981  the given <tt>innerScale</tt>, followed by Gaussian smoothing at <tt>outerScale</tt>.
982  Both source and destination arrays are represented by iterators, shape objects and
983  accessors. The destination array must have a vector valued pixel type with
984  N*(N+1)/2 elements (it represents the upper triangular part of the symmetric
985  structure tensor matrix). If the source array is also vector valued, the
986  resulting structure tensor is the sum of the individual tensors for each channel.
987  This function is implemented by calls to
988  \ref separableConvolveMultiArray() with the appropriate kernels.
989 
990  <b> Declarations:</b>
991 
992  pass arguments explicitly:
993  \code
994  namespace vigra {
995  template <class SrcIterator, class SrcShape, class SrcAccessor,
996  class DestIterator, class DestAccessor>
997  void
998  structureTensorMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
999  DestIterator diter, DestAccessor dest,
1000  double innerScale, double outerScale);
1001  }
1002  \endcode
1003 
1004  use argument objects in conjunction with \ref ArgumentObjectFactories :
1005  \code
1006  namespace vigra {
1007  template <class SrcIterator, class SrcShape, class SrcAccessor,
1008  class DestIterator, class DestAccessor>
1009  void
1010  structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1011  pair<DestIterator, DestAccessor> const & dest,
1012  double innerScale, double outerScale);
1013  }
1014  \endcode
1015 
1016  <b> Usage:</b>
1017 
1018  <b>\#include</b> <<a href="multi__convolution_8hxx-source.html">vigra/multi_convolution.hxx</a>>
1019 
1020  \code
1021  MultiArray<3, RGBValue<float> > source(shape);
1022  MultiArray<3, TinyVector<float, 6> > dest(shape);
1023  ...
1024  // compute structure tensor at scales innerScale and outerScale
1025  structureTensorMultiArray(srcMultiArrayRange(source), destMultiArray(dest), innerScale, outerScale);
1026  \endcode
1027 
1028  <b> Required Interface:</b>
1029 
1030  see \ref separableConvolveMultiArray(), in addition:
1031 
1032  \code
1033  int dimension = 0;
1034  VectorElementAccessor<DestAccessor> elementAccessor(0, dest);
1035  \endcode
1036 
1037  \see separableConvolveMultiArray(), vectorToTensorMultiArray()
1038 */
1040 
1041 template <class SrcIterator, class SrcShape, class SrcAccessor,
1042  class DestIterator, class DestAccessor>
1043 void
1044 structureTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1045  DestIterator di, DestAccessor dest,
1046  double innerScale, double outerScale)
1047 {
1048  static const int N = SrcShape::static_size;
1049  static const int M = N*(N+1)/2;
1050 
1051  typedef typename DestAccessor::value_type DestType;
1052  typedef typename DestType::value_type DestValueType;
1053  typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1054  typedef TinyVector<KernelType, N> GradientVector;
1055  typedef typename AccessorTraits<GradientVector>::default_accessor GradientAccessor;
1056 
1057  for(int k=0; k<N; ++k)
1058  if(shape[k] <=0)
1059  return;
1060 
1061  vigra_precondition(M == dest.size(di),
1062  "structureTensorMultiArray(): Wrong number of channels in output array.");
1063 
1064  vigra_precondition(innerScale > 0.0 && outerScale >= 0.0,
1065  "structureTensorMultiArray(): Scale must be positive.");
1066 
1067  MultiArray<N, GradientVector> gradient(shape);
1068  gaussianGradientMultiArray(si, shape, src,
1069  gradient.traverser_begin(), GradientAccessor(),
1070  innerScale);
1071 
1072  transformMultiArray(gradient.traverser_begin(), shape, GradientAccessor(),
1073  di, dest,
1074  detail::StructurTensorFunctor<N, DestType>());
1075 
1076  gaussianSmoothMultiArray(di, shape, dest, di, dest, outerScale);
1077 }
1078 
1079 template <class SrcIterator, class SrcShape, class SrcAccessor,
1080  class DestIterator, class DestAccessor>
1081 inline void
1082 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1083  pair<DestIterator, DestAccessor> const & dest,
1084  double innerScale, double outerScale)
1085 {
1086  structureTensorMultiArray( source.first, source.second, source.third,
1087  dest.first, dest.second, innerScale, outerScale );
1088 }
1089 
1090 //@}
1091 
1092 } //-- namespace vigra
1093 
1094 
1095 #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.7.1 (Thu Jun 14 2012)