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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.11.2