[ 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
57 namespace vigra
58 {
59
60 namespace detail
61 {
62
63 struct DoubleYielder
64 {
65 const double value;
66 DoubleYielder(double v, unsigned, const char *const) : value(v) {}
67 DoubleYielder(double v) : value(v) {}
68 void operator++() {}
69 double operator*() const { return value; }
70 };
71
72 template <typename X>
73 struct IteratorDoubleYielder
74 {
75 X it;
76 IteratorDoubleYielder(X i, unsigned, const char *const) : it(i) {}
77 IteratorDoubleYielder(X i) : it(i) {}
78 void operator++() { ++it; }
79 double operator*() const { return *it; }
80 };
81
82 template <typename X>
83 struct SequenceDoubleYielder
84 {
85 typename X::const_iterator it;
86 SequenceDoubleYielder(const X & seq, unsigned dim,
87 const char *const function_name = "SequenceDoubleYielder")
88 : it(seq.begin())
89 {
90 if (seq.size() == dim)
91 return;
92 std::string msg = "(): Parameter number be equal to the number of spatial dimensions.";
93 vigra_precondition(false, function_name + msg);
94 }
95 void operator++() { ++it; }
96 double operator*() const { return *it; }
97 };
98
99 template <typename X>
100 struct WrapDoubleIterator
101 {
102 typedef
103 typename IfBool< IsConvertibleTo<X, double>::value,
104 DoubleYielder,
105 typename IfBool< IsIterator<X>::value || IsArray<X>::value,
106 IteratorDoubleYielder<X>,
107 SequenceDoubleYielder<X>
108 >::type
109 >::type type;
110 };
111
112 template <class Param1, class Param2, class Param3>
113 struct WrapDoubleIteratorTriple
114 {
115 typename WrapDoubleIterator<Param1>::type sigma_eff_it;
116 typename WrapDoubleIterator<Param2>::type sigma_d_it;
117 typename WrapDoubleIterator<Param3>::type step_size_it;
118 WrapDoubleIteratorTriple(Param1 sigma_eff, Param2 sigma_d, Param3 step_size)
119 : sigma_eff_it(sigma_eff), sigma_d_it(sigma_d), step_size_it(step_size) {}
120 void operator++()
121 {
122 ++sigma_eff_it;
123 ++sigma_d_it;
124 ++step_size_it;
125 }
126 double sigma_eff() const { return *sigma_eff_it; }
127 double sigma_d() const { return *sigma_d_it; }
128 double step_size() const { return *step_size_it; }
129 static void sigma_precondition(double sigma, const char *const function_name)
130 {
131 if (sigma < 0.0)
132 {
133 std::string msg = "(): Scale must be positive.";
134 vigra_precondition(false, function_name + msg);
135 }
136 }
137 double sigma_scaled(const char *const function_name = "unknown function ",
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
158 template <unsigned dim>
159 struct 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 */
334 template <unsigned dim>
335 class ConvolutionOptions
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
349 ConvolutionOptions()
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 */
394 ConvolutionOptions<dim> & stepSize(...);
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 */
410 ConvolutionOptions<dim> & resolutionStdDev(...);
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 */
429 ConvolutionOptions<dim> & stdDev(...);
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 */
441 ConvolutionOptions<dim> & innerScale(...);
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 */
459 ConvolutionOptions<dim> & outerScale(...);
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 */
476 ConvolutionOptions<dim> & filterWindowSize(double ratio)
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 */
498 ConvolutionOptions<dim> & subarray(Shape const & from, Shape const & to)
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
513 namespace detail
514 {
515
516 /********************************************************/
517 /* */
518 /* internalSeparableConvolveMultiArray */
519 /* */
520 /********************************************************/
521
522 template <class SrcIterator, class SrcShape, class SrcAccessor,
523 class DestIterator, class DestAccessor, class KernelIterator>
524 void
525 internalSeparableConvolveMultiArrayTmp(
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
584 template <class SrcIterator, class SrcShape, class SrcAccessor,
585 class DestIterator, class DestAccessor, class KernelIterator>
586 void
587 internalSeparableConvolveSubarray(
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
674 template <class K>
675 void
676 scaleKernel(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 MultiArrayConvolutionFilters Convolution filters for multi-dimensional arrays.
686
687 These functions realize a separable convolution on an arbitrary dimensional
688 array that is specified by iterators (compatible to \ref MultiIteratorPage)
689 and shape objects. It can therefore be applied to a wide range of data structures
690 (\ref vigra::MultiArrayView, \ref vigra::MultiArray etc.).
691 */
692 //@{
693
694 /********************************************************/
695 /* */
696 /* separableConvolveMultiArray */
697 /* */
698 /********************************************************/
699
700 /** \brief Separated convolution on multi-dimensional arrays.
701
702 This function computes a separated convolution on all dimensions
703 of the given multi-dimensional array. Both source and destination
704 arrays are represented by iterators, shape objects and accessors.
705 The destination array is required to already have the correct size.
706
707 There are two variants of this functions: one takes a single kernel
708 of type \ref vigra::Kernel1D which is then applied to all dimensions,
709 whereas the other requires an iterator referencing a sequence of
710 \ref vigra::Kernel1D objects, one for every dimension of the data.
711 Then the first kernel in this sequence is applied to the innermost
712 dimension (e.g. the x-axis of an image), while the last is applied to the
713 outermost dimension (e.g. the z-axis in a 3D image).
714
715 This function may work in-place, which means that <tt>source.data() == dest.data()</tt> is allowed.
716 A full-sized internal array is only allocated if working on the destination
717 array directly would cause round-off errors (i.e. if
718 <tt>typeid(typename NumericTraits<T2>::RealPromote) != typeid(T2)</tt>).
719
720 If <tt>start</tt> and <tt>stop</tt> have non-default values, they must represent
721 a valid subarray of the input array. The convolution is then restricted to that
722 subarray, and it is assumed that the output array only refers to the
723 subarray (i.e. <tt>dest.shape() == stop - start</tt>). Negative ROI boundaries are
724 interpreted relative to the end of the respective dimension
725 (i.e. <tt>if(stop[k] < 0) stop[k] += source.shape(k);</tt>).
726
727 <b> Declarations:</b>
728
729 pass arbitrary-dimensional array views:
730 \code
731 namespace vigra {
732 // apply each kernel from the sequence 'kernels' in turn
733 template <unsigned int N, class T1, class S1,
734 class T2, class S2,
735 class KernelIterator>
736 void
737 separableConvolveMultiArray(MultiArrayView<N, T1, S1> const & source,
738 MultiArrayView<N, T2, S2> dest,
739 KernelIterator kernels,
740 typename MultiArrayShape<N>::type start = typename MultiArrayShape<N>::type(),
741 typename MultiArrayShape<N>::type stop = typename MultiArrayShape<N>::type());
742
743 // apply the same kernel to all dimensions
744 template <unsigned int N, class T1, class S1,
745 class T2, class S2,
746 class T>
747 void
748 separableConvolveMultiArray(MultiArrayView<N, T1, S1> const & source,
749 MultiArrayView<N, T2, S2> dest,
750 Kernel1D<T> const & kernel,
751 typename MultiArrayShape<N>::type const & start = typename MultiArrayShape<N>::type(),
752 typename MultiArrayShape<N>::type const & stop = typename MultiArrayShape<N>::type());
753 }
754 \endcode
755
756 \deprecatedAPI{separableConvolveMultiArray}
757 pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
758 \code
759 namespace vigra {
760 // apply the same kernel to all dimensions
761 template <class SrcIterator, class SrcShape, class SrcAccessor,
762 class DestIterator, class DestAccessor, class T>
763 void
764 separableConvolveMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
765 DestIterator diter, DestAccessor dest,
766 Kernel1D<T> const & kernel,
767 SrcShape const & start = SrcShape(),
768 SrcShape const & stop = SrcShape());
769
770 // apply each kernel from the sequence 'kernels' in turn
771 template <class SrcIterator, class SrcShape, class SrcAccessor,
772 class DestIterator, class DestAccessor, class KernelIterator>
773 void
774 separableConvolveMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
775 DestIterator diter, DestAccessor dest,
776 KernelIterator kernels,
777 SrcShape const & start = SrcShape(),
778 SrcShape const & stop = SrcShape());
779 }
780 \endcode
781 use argument objects in conjunction with \ref ArgumentObjectFactories :
782 \code
783 namespace vigra {
784 // apply the same kernel to all dimensions
785 template <class SrcIterator, class SrcShape, class SrcAccessor,
786 class DestIterator, class DestAccessor, class T>
787 void
788 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
789 pair<DestIterator, DestAccessor> const & dest,
790 Kernel1D<T> const & kernel,
791 SrcShape const & start = SrcShape(),
792 SrcShape const & stop = SrcShape());
793
794 // apply each kernel from the sequence 'kernels' in turn
795 template <class SrcIterator, class SrcShape, class SrcAccessor,
796 class DestIterator, class DestAccessor, class KernelIterator>
797 void
798 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
799 pair<DestIterator, DestAccessor> const & dest,
800 KernelIterator kernels,
801 SrcShape const & start = SrcShape(),
802 SrcShape const & stop = SrcShape());
803 }
804 \endcode
805 \deprecatedEnd
806
807 <b> Usage:</b>
808
809 <b>\#include</b> <vigra/multi_convolution.hxx><br/>
810 Namespace: vigra
811
812 \code
813 Shape3 shape(width, height, depth);
814 MultiArray<3, unsigned char> source(shape);
815 MultiArray<3, float> dest(shape);
816 ...
817 Kernel1D<float> gauss;
818 gauss.initGaussian(sigma);
819
820 // smooth all dimensions with the same kernel
821 separableConvolveMultiArray(source, dest, gauss);
822
823 // create 3 Gauss kernels, one for each dimension, but smooth the z-axis less
824 ArrayVector<Kernel1D<float> > kernels(3, gauss);
825 kernels[2].initGaussian(sigma / 2.0);
826
827 // perform Gaussian smoothing on all dimensions
828 separableConvolveMultiArray(source, dest, kernels.begin());
829
830 // create output array for a ROI
831 MultiArray<3, float> destROI(shape - Shape3(10,10,10));
832
833 // only smooth the given ROI (ignore 5 pixels on all sides of the array)
834 separableConvolveMultiArray(source, destROI, gauss, Shape3(5,5,5), Shape3(-5,-5,-5));
835 \endcode
836
837 \deprecatedUsage{separableConvolveMultiArray}
838 \code
839 MultiArray<3, unsigned char>::size_type shape(width, height, depth);
840 MultiArray<3, unsigned char> source(shape);
841 MultiArray<3, float> dest(shape);
842 ...
843 Kernel1D<float> gauss;
844 gauss.initGaussian(sigma);
845 // create 3 Gauss kernels, one for each dimension
846 ArrayVector<Kernel1D<float> > kernels(3, gauss);
847
848 // perform Gaussian smoothing on all dimensions
849 separableConvolveMultiArray(source, dest,
850 kernels.begin());
851 \endcode
852 <b> Required Interface:</b>
853 \code
854 see \ref separableConvolveImage(), in addition:
855
856 NumericTraits<T1>::RealPromote s = src[0];
857
858 s = s + s;
859 s = kernel(0) * s;
860 \endcode
861 \deprecatedEnd
862
863 \see vigra::Kernel1D, convolveLine()
864 */
865 doxygen_overloaded_function(template <...> void separableConvolveMultiArray)
866
867 template <class SrcIterator, class SrcShape, class SrcAccessor,
868 class DestIterator, class DestAccessor, class KernelIterator>
869 void
870 separableConvolveMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
871 DestIterator d, DestAccessor dest,
872 KernelIterator kernels,
873 SrcShape start = SrcShape(),
874 SrcShape stop = SrcShape())
875 {
876 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
877
878
879 if(stop != SrcShape())
880 {
881
882 enum { N = 1 + SrcIterator::level };
883 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, start);
884 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, stop);
885
886 for(int k=0; k<N; ++k)
887 vigra_precondition(0 <= start[k] && start[k] < stop[k] && stop[k] <= shape[k],
888 "separableConvolveMultiArray(): invalid subarray shape.");
889
890 detail::internalSeparableConvolveSubarray(s, shape, src, d, dest, kernels, start, stop);
891 }
892 else if(!IsSameType<TmpType, typename DestAccessor::value_type>::boolResult)
893 {
894 // need a temporary array to avoid rounding errors
895 MultiArray<SrcShape::static_size, TmpType> tmpArray(shape);
896 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src,
897 tmpArray.traverser_begin(), typename AccessorTraits<TmpType>::default_accessor(), kernels );
898 copyMultiArray(srcMultiArrayRange(tmpArray), destIter(d, dest));
899 }
900 else
901 {
902 // work directly on the destination array
903 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, d, dest, kernels );
904 }
905 }
906
907 template <class SrcIterator, class SrcShape, class SrcAccessor,
908 class DestIterator, class DestAccessor, class T>
909 inline void
910 separableConvolveMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
911 DestIterator d, DestAccessor dest,
912 Kernel1D<T> const & kernel,
913 SrcShape const & start = SrcShape(),
914 SrcShape const & stop = SrcShape())
915 {
916 ArrayVector<Kernel1D<T> > kernels(shape.size(), kernel);
917
918 separableConvolveMultiArray( s, shape, src, d, dest, kernels.begin(), start, stop);
919 }
920
921 template <class SrcIterator, class SrcShape, class SrcAccessor,
922 class DestIterator, class DestAccessor, class KernelIterator>
923 inline void
924 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
925 pair<DestIterator, DestAccessor> const & dest,
926 KernelIterator kit,
927 SrcShape const & start = SrcShape(),
928 SrcShape const & stop = SrcShape())
929 {
930 separableConvolveMultiArray( source.first, source.second, source.third,
931 dest.first, dest.second, kit, start, stop );
932 }
933
934 template <class SrcIterator, class SrcShape, class SrcAccessor,
935 class DestIterator, class DestAccessor, class T>
936 inline void
937 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
938 pair<DestIterator, DestAccessor> const & dest,
939 Kernel1D<T> const & kernel,
940 SrcShape const & start = SrcShape(),
941 SrcShape const & stop = SrcShape())
942 {
943 ArrayVector<Kernel1D<T> > kernels(source.second.size(), kernel);
944
945 separableConvolveMultiArray( source.first, source.second, source.third,
946 dest.first, dest.second, kernels.begin(), start, stop);
947 }
948
949 template <unsigned int N, class T1, class S1,
950 class T2, class S2,
951 class KernelIterator>
952 inline void
953 separableConvolveMultiArray(MultiArrayView<N, T1, S1> const & source,
954 MultiArrayView<N, T2, S2> dest,
955 KernelIterator kit,
956 typename MultiArrayShape<N>::type start = typename MultiArrayShape<N>::type(),
957 typename MultiArrayShape<N>::type stop = typename MultiArrayShape<N>::type())
958 {
959 if(stop != typename MultiArrayShape<N>::type())
960 {
961 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), start);
962 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), stop);
963 vigra_precondition(dest.shape() == (stop - start),
964 "separableConvolveMultiArray(): shape mismatch between ROI and output.");
965 }
966 else
967 {
968 vigra_precondition(source.shape() == dest.shape(),
969 "separableConvolveMultiArray(): shape mismatch between input and output.");
970 }
971 separableConvolveMultiArray( srcMultiArrayRange(source),
972 destMultiArray(dest), kit, start, stop );
973 }
974
975 template <unsigned int N, class T1, class S1,
976 class T2, class S2,
977 class T>
978 inline void
979 separableConvolveMultiArray(MultiArrayView<N, T1, S1> const & source,
980 MultiArrayView<N, T2, S2> dest,
981 Kernel1D<T> const & kernel,
982 typename MultiArrayShape<N>::type const & start = typename MultiArrayShape<N>::type(),
983 typename MultiArrayShape<N>::type const & stop = typename MultiArrayShape<N>::type())
984 {
985 ArrayVector<Kernel1D<T> > kernels(N, kernel);
986 separableConvolveMultiArray(source, dest, kernels.begin(), start, stop);
987 }
988
989 /********************************************************/
990 /* */
991 /* convolveMultiArrayOneDimension */
992 /* */
993 /********************************************************/
994
995 /** \brief Convolution along a single dimension of a multi-dimensional arrays.
996
997 This function computes a convolution along one dimension (specified by
998 the parameter <tt>dim</tt> of the given multi-dimensional array with the given
999 <tt>kernel</tt>. The destination array must already have the correct size.
1000
1001 If <tt>start</tt> and <tt>stop</tt> have non-default values, they must represent
1002 a valid subarray of the input array. The convolution is then restricted to that
1003 subarray, and it is assumed that the output array only refers to the
1004 subarray (i.e. <tt>dest.shape() == stop - start</tt>). Negative ROI boundaries are
1005 interpreted relative to the end of the respective dimension
1006 (i.e. <tt>if(stop[k] < 0) stop[k] += source.shape(k);</tt>).
1007
1008 This function may work in-place, which means that <tt>source.data() == dest.data()</tt> is allowed.
1009
1010 <b> Declarations:</b>
1011
1012 pass arbitrary-dimensional array views:
1013 \code
1014 namespace vigra {
1015 template <unsigned int N, class T1, class S1,
1016 class T2, class S2,
1017 class T>
1018 void
1019 convolveMultiArrayOneDimension(MultiArrayView<N, T1, S1> const & source,
1020 MultiArrayView<N, T2, S2> dest,
1021 unsigned int dim,
1022 Kernel1D<T> const & kernel,
1023 typename MultiArrayShape<N>::type start = typename MultiArrayShape<N>::type(),
1024 typename MultiArrayShape<N>::type stop = typename MultiArrayShape<N>::type());
1025 }
1026 \endcode
1027
1028 \deprecatedAPI{convolveMultiArrayOneDimension}
1029 pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
1030 \code
1031 namespace vigra {
1032 template <class SrcIterator, class SrcShape, class SrcAccessor,
1033 class DestIterator, class DestAccessor, class T>
1034 void
1035 convolveMultiArrayOneDimension(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1036 DestIterator diter, DestAccessor dest,
1037 unsigned int dim, vigra::Kernel1D<T> const & kernel,
1038 SrcShape const & start = SrcShape(),
1039 SrcShape const & stop = SrcShape());
1040 }
1041 \endcode
1042 use argument objects in conjunction with \ref ArgumentObjectFactories :
1043 \code
1044 namespace vigra {
1045 template <class SrcIterator, class SrcShape, class SrcAccessor,
1046 class DestIterator, class DestAccessor, class T>
1047 void
1048 convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1049 pair<DestIterator, DestAccessor> const & dest,
1050 unsigned int dim, vigra::Kernel1D<T> const & kernel,
1051 SrcShape const & start = SrcShape(),
1052 SrcShape const & stop = SrcShape());
1053 }
1054 \endcode
1055 \deprecatedEnd
1056
1057 <b> Usage:</b>
1058
1059 <b>\#include</b> <vigra/multi_convolution.hxx><br/>
1060 Namespace: vigra
1061
1062 \code
1063 Shape3 shape(width, height, depth);
1064 MultiArray<3, unsigned char> source(shape);
1065 MultiArray<3, float> dest(shape);
1066 ...
1067 Kernel1D<float> gauss;
1068 gauss.initGaussian(sigma);
1069
1070 // perform Gaussian smoothing along dimension 1 (height)
1071 convolveMultiArrayOneDimension(source, dest, 1, gauss);
1072 \endcode
1073
1074 \see separableConvolveMultiArray()
1075 */
1076 doxygen_overloaded_function(template <...> void convolveMultiArrayOneDimension)
1077
1078 template <class SrcIterator, class SrcShape, class SrcAccessor,
1079 class DestIterator, class DestAccessor, class T>
1080 void
1081 convolveMultiArrayOneDimension(SrcIterator s, SrcShape const & shape, SrcAccessor src,
1082 DestIterator d, DestAccessor dest,
1083 unsigned int dim, vigra::Kernel1D<T> const & kernel,
1084 SrcShape const & start = SrcShape(),
1085 SrcShape const & stop = SrcShape())
1086 {
1087 enum { N = 1 + SrcIterator::level };
1088 vigra_precondition( dim < N,
1089 "convolveMultiArrayOneDimension(): The dimension number to convolve must be smaller "
1090 "than the data dimensionality" );
1091
1092 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
1093 typedef typename AccessorTraits<TmpType>::default_const_accessor TmpAccessor;
1094 ArrayVector<TmpType> tmp( shape[dim] );
1095
1096 typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
1097 typedef MultiArrayNavigator<DestIterator, N> DNavigator;
1098
1099 SrcShape sstart, sstop(shape), dstart, dstop(shape);
1100
1101 if(stop != SrcShape())
1102 {
1103 sstart = start;
1104 sstop = stop;
1105 sstart[dim] = 0;
1106 sstop[dim] = shape[dim];
1107 dstop = stop - start;
1108 }
1109
1110 SNavigator snav( s, sstart, sstop, dim );
1111 DNavigator dnav( d, dstart, dstop, dim );
1112
1113 for( ; snav.hasMore(); snav++, dnav++ )
1114 {
1115 // first copy source to temp for maximum cache efficiency
1116 copyLine(snav.begin(), snav.end(), src,
1117 tmp.begin(), typename AccessorTraits<TmpType>::default_accessor() );
1118
1119 convolveLine(srcIterRange( tmp.begin(), tmp.end(), TmpAccessor()),
1120 destIter( dnav.begin(), dest ),
1121 kernel1d( kernel), start[dim], stop[dim]);
1122 }
1123 }
1124
1125 template <class SrcIterator, class SrcShape, class SrcAccessor,
1126 class DestIterator, class DestAccessor, class T>
1127 inline void
1128 convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1129 pair<DestIterator, DestAccessor> const & dest,
1130 unsigned int dim,
1131 Kernel1D<T> const & kernel,
1132 SrcShape const & start = SrcShape(),
1133 SrcShape const & stop = SrcShape())
1134 {
1135 convolveMultiArrayOneDimension(source.first, source.second, source.third,
1136 dest.first, dest.second, dim, kernel, start, stop);
1137 }
1138
1139 template <unsigned int N, class T1, class S1,
1140 class T2, class S2,
1141 class T>
1142 inline void
1143 convolveMultiArrayOneDimension(MultiArrayView<N, T1, S1> const & source,
1144 MultiArrayView<N, T2, S2> dest,
1145 unsigned int dim,
1146 Kernel1D<T> const & kernel,
1147 typename MultiArrayShape<N>::type start = typename MultiArrayShape<N>::type(),
1148 typename MultiArrayShape<N>::type stop = typename MultiArrayShape<N>::type())
1149 {
1150 if(stop != typename MultiArrayShape<N>::type())
1151 {
1152 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), start);
1153 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), stop);
1154 vigra_precondition(dest.shape() == (stop - start),
1155 "convolveMultiArrayOneDimension(): shape mismatch between ROI and output.");
1156 }
1157 else
1158 {
1159 vigra_precondition(source.shape() == dest.shape(),
1160 "convolveMultiArrayOneDimension(): shape mismatch between input and output.");
1161 }
1162 convolveMultiArrayOneDimension(srcMultiArrayRange(source),
1163 destMultiArray(dest), dim, kernel, start, stop);
1164 }
1165
1166 /********************************************************/
1167 /* */
1168 /* gaussianSmoothMultiArray */
1169 /* */
1170 /********************************************************/
1171
1172 /** \brief Isotropic Gaussian smoothing of a multi-dimensional arrays.
1173
1174 This function computes an isotropic convolution of the given N-dimensional
1175 array with a Gaussian filter at the given standard deviation <tt>sigma</tt>.
1176 Both source and destination arrays are represented by
1177 iterators, shape objects and accessors. The destination array is required to
1178 already have the correct size. This function may work in-place, which means
1179 that <tt>source.data() == dest.data()</tt> is allowed. It is implemented by a call to
1180 \ref separableConvolveMultiArray() with the appropriate kernel.
1181
1182 Anisotropic data should be provided with appropriate \ref vigra::ConvolutionOptions
1183 to adjust the filter sizes for the resolution of each axis.
1184 Otherwise, the parameter <tt>opt</tt> is optional unless the parameter
1185 <tt>sigma</tt> is omitted.
1186
1187 If you pass \ref vigra::BlockwiseConvolutionOptions instead, the algorithm will
1188 be executed in parallel on data blocks of a certain size. The block size can be
1189 customized via <tt>BlockwiseConvolutionOptions::blockShape()</tt>, but the defaults
1190 usually work reasonably. By default, the number of threads equals the capabilities
1191 of your hardware, but you can change this via <tt>BlockwiseConvolutionOptions::numThreads()</tt>.
1192
1193 <b> Declarations:</b>
1194
1195 pass arbitrary-dimensional array views:
1196 \code
1197 namespace vigra {
1198 // pass filter scale explicitly
1199 template <unsigned int N, class T1, class S1,
1200 class T2, class S2>
1201 void
1202 gaussianSmoothMultiArray(MultiArrayView<N, T1, S1> const & source,
1203 MultiArrayView<N, T2, S2> dest,
1204 double sigma,
1205 ConvolutionOptions<N> opt = ConvolutionOptions<N>());
1206
1207 // pass filer scale(s) in the option object
1208 template <unsigned int N, class T1, class S1,
1209 class T2, class S2>
1210 void
1211 gaussianSmoothMultiArray(MultiArrayView<N, T1, S1> const & source,
1212 MultiArrayView<N, T2, S2> dest,
1213 ConvolutionOptions<N> opt);
1214
1215 // as above, but execute algorirhm in parallel
1216 template <unsigned int N, class T1, class S1,
1217 class T2, class S2>
1218 void
1219 gaussianSmoothMultiArray(MultiArrayView<N, T1, S1> const & source,
1220 MultiArrayView<N, T2, S2> dest,
1221 BlockwiseConvolutionOptions<N> opt);
1222 }
1223 \endcode
1224
1225 \deprecatedAPI{gaussianSmoothMultiArray}
1226 pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
1227 \code
1228 namespace vigra {
1229 template <class SrcIterator, class SrcShape, class SrcAccessor,
1230 class DestIterator, class DestAccessor>
1231 void
1232 gaussianSmoothMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1233 DestIterator diter, DestAccessor dest,
1234 double sigma,
1235 const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1236 }
1237 \endcode
1238 use argument objects in conjunction with \ref ArgumentObjectFactories :
1239 \code
1240 namespace vigra {
1241 template <class SrcIterator, class SrcShape, class SrcAccessor,
1242 class DestIterator, class DestAccessor>
1243 void
1244 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1245 pair<DestIterator, DestAccessor> const & dest,
1246 double sigma,
1247 const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1248 }
1249 \endcode
1250 \deprecatedEnd
1251
1252 <b> Usage:</b>
1253
1254 <b>\#include</b> <vigra/multi_convolution.hxx> (sequential version)<br/>
1255 <b>\#include</b> <vigra/multi_blockwise.hxx> (parallel version)<br/>
1256 Namespace: vigra
1257
1258 \code
1259 Shape3 shape(width, height, depth);
1260 MultiArray<3, unsigned char> source(shape);
1261 MultiArray<3, float> dest(shape);
1262 ...
1263 // perform isotropic Gaussian smoothing at scale 'sigma'
1264 gaussianSmoothMultiArray(source, dest, sigma);
1265 \endcode
1266
1267 <b> Multi-threaded execution:</b>
1268
1269 \code
1270 Shape3 shape(width, height, depth);
1271 MultiArray<3, unsigned char> source(shape);
1272 MultiArray<3, float> dest(shape);
1273 ...
1274 BlockwiseConvolutionOptions<3> opt;
1275 opt.numThreads(4); // use 4 threads (uses hardware default if not given)
1276 opt.innerScale(sigma);
1277
1278 // perform isotropic Gaussian smoothing at scale 'sigma' in parallel
1279 gaussianSmoothMultiArray(source, dest, sigma, opt);
1280 \endcode
1281
1282 <b> Usage with anisotropic data:</b>
1283
1284 \code
1285 Shape3 shape(width, height, depth);
1286 MultiArray<3, unsigned char> source(shape);
1287 MultiArray<3, float> dest(shape);
1288 TinyVector<float, 3> step_size;
1289 TinyVector<float, 3> resolution_sigmas;
1290 ...
1291 // perform anisotropic Gaussian smoothing at scale 'sigma'
1292 gaussianSmoothMultiArray(source, dest, sigma,
1293 ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
1294 \endcode
1295
1296 \see separableConvolveMultiArray()
1297 */
1298 doxygen_overloaded_function(template <...> void gaussianSmoothMultiArray)
1299
1300 template <class SrcIterator, class SrcShape, class SrcAccessor,
1301 class DestIterator, class DestAccessor>
1302 void
1303 gaussianSmoothMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
1304 DestIterator d, DestAccessor dest,
1305 const ConvolutionOptions<SrcShape::static_size> & opt,
1306 const char *const function_name = "gaussianSmoothMultiArray" )
1307 {
1308 static const int N = SrcShape::static_size;
1309
1310 typename ConvolutionOptions<N>::ScaleIterator params = opt.scaleParams();
1311 ArrayVector<Kernel1D<double> > kernels(N);
1312
1313 for (int dim = 0; dim < N; ++dim, ++params)
1314 kernels[dim].initGaussian(params.sigma_scaled(function_name, true),
1315 1.0, opt.window_ratio);
1316
1317 separableConvolveMultiArray(s, shape, src, d, dest, kernels.begin(), opt.from_point, opt.to_point);
1318 }
1319
1320 template <class SrcIterator, class SrcShape, class SrcAccessor,
1321 class DestIterator, class DestAccessor>
1322 inline void
1323 gaussianSmoothMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src,
1324 DestIterator d, DestAccessor dest, double sigma,
1325 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1326 {
1327 ConvolutionOptions<SrcShape::static_size> par = opt;
1328 gaussianSmoothMultiArray(s, shape, src, d, dest, par.stdDev(sigma));
1329 }
1330
1331 template <class SrcIterator, class SrcShape, class SrcAccessor,
1332 class DestIterator, class DestAccessor>
1333 inline void
1334 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1335 pair<DestIterator, DestAccessor> const & dest,
1336 const ConvolutionOptions<SrcShape::static_size> & opt)
1337 {
1338 gaussianSmoothMultiArray( source.first, source.second, source.third,
1339 dest.first, dest.second, opt );
1340 }
1341
1342 template <class SrcIterator, class SrcShape, class SrcAccessor,
1343 class DestIterator, class DestAccessor>
1344 inline void
1345 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1346 pair<DestIterator, DestAccessor> const & dest, double sigma,
1347 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1348 {
1349 gaussianSmoothMultiArray( source.first, source.second, source.third,
1350 dest.first, dest.second, sigma, opt );
1351 }
1352
1353 template <unsigned int N, class T1, class S1,
1354 class T2, class S2>
1355 inline void
1356 gaussianSmoothMultiArray(MultiArrayView<N, T1, S1> const & source,
1357 MultiArrayView<N, T2, S2> dest,
1358 ConvolutionOptions<N> opt)
1359 {
1360 if(opt.to_point != typename MultiArrayShape<N>::type())
1361 {
1362 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
1363 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
1364 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1365 "gaussianSmoothMultiArray(): shape mismatch between ROI and output.");
1366 }
1367 else
1368 {
1369 vigra_precondition(source.shape() == dest.shape(),
1370 "gaussianSmoothMultiArray(): shape mismatch between input and output.");
1371 }
1372
1373 gaussianSmoothMultiArray( srcMultiArrayRange(source),
1374 destMultiArray(dest), opt );
1375 }
1376
1377 template <unsigned int N, class T1, class S1,
1378 class T2, class S2>
1379 inline void
1380 gaussianSmoothMultiArray(MultiArrayView<N, T1, S1> const & source,
1381 MultiArrayView<N, T2, S2> dest,
1382 double sigma,
1383 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
1384 {
1385 gaussianSmoothMultiArray( source, dest, opt.stdDev(sigma) );
1386 }
1387
1388
1389 /********************************************************/
1390 /* */
1391 /* gaussianGradientMultiArray */
1392 /* */
1393 /********************************************************/
1394
1395 /** \brief Calculate Gaussian gradient of a multi-dimensional arrays.
1396
1397 This function computes the Gaussian gradient of the given N-dimensional
1398 array with a sequence of first-derivative-of-Gaussian filters at the given
1399 standard deviation <tt>sigma</tt> (differentiation is applied to each dimension
1400 in turn, starting with the innermost dimension). The destination array is
1401 required to have a vector valued pixel type with as many elements as the number of
1402 dimensions. This function is implemented by calls to
1403 \ref separableConvolveMultiArray() with the appropriate kernels.
1404
1405 Anisotropic data should be provided with appropriate \ref vigra::ConvolutionOptions
1406 to adjust the filter sizes for the resolution of each axis.
1407 Otherwise, the parameter <tt>opt</tt> is optional unless the parameter
1408 <tt>sigma</tt> is omitted.
1409
1410 If you pass \ref vigra::BlockwiseConvolutionOptions instead, the algorithm will
1411 be executed in parallel on data blocks of a certain size. The block size can be
1412 customized via <tt>BlockwiseConvolutionOptions::blockShape()</tt>, but the defaults
1413 usually work reasonably. By default, the number of threads equals the capabilities
1414 of your hardware, but you can change this via <tt>BlockwiseConvolutionOptions::numThreads()</tt>.
1415
1416 <b> Declarations:</b>
1417
1418 pass arbitrary-dimensional array views:
1419 \code
1420 namespace vigra {
1421 // pass filter scale explicitly
1422 template <unsigned int N, class T1, class S1,
1423 class T2, class S2>
1424 void
1425 gaussianGradientMultiArray(MultiArrayView<N, T1, S1> const & source,
1426 MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1427 double sigma,
1428 ConvolutionOptions<N> opt = ConvolutionOptions<N>());
1429
1430 // pass filter scale(s) in option object
1431 template <unsigned int N, class T1, class S1,
1432 class T2, class S2>
1433 void
1434 gaussianGradientMultiArray(MultiArrayView<N, T1, S1> const & source,
1435 MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1436 ConvolutionOptions<N> opt);
1437
1438 // likewise, but execute algorithm in parallel
1439 template <unsigned int N, class T1, class S1,
1440 class T2, class S2>
1441 void
1442 gaussianGradientMultiArray(MultiArrayView<N, T1, S1> const & source,
1443 MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1444 BlockwiseConvolutionOptions<N> opt);
1445 }
1446 \endcode
1447
1448 \deprecatedAPI{gaussianGradientMultiArray}
1449 pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
1450 \code
1451 namespace vigra {
1452 template <class SrcIterator, class SrcShape, class SrcAccessor,
1453 class DestIterator, class DestAccessor>
1454 void
1455 gaussianGradientMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1456 DestIterator diter, DestAccessor dest,
1457 double sigma,
1458 const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1459 }
1460 \endcode
1461 use argument objects in conjunction with \ref ArgumentObjectFactories :
1462 \code
1463 namespace vigra {
1464 template <class SrcIterator, class SrcShape, class SrcAccessor,
1465 class DestIterator, class DestAccessor>
1466 void
1467 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1468 pair<DestIterator, DestAccessor> const & dest,
1469 double sigma,
1470 const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1471 }
1472 \endcode
1473 \deprecatedEnd
1474
1475 <b> Usage:</b>
1476
1477 <b>\#include</b> <vigra/multi_convolution.hxx> (sequential version)<br/>
1478 <b>\#include</b> <vigra/multi_blockwise.hxx> (parallel version)<br/>
1479 Namespace: vigra
1480
1481 \code
1482 Shape3 shape(width, height, depth);
1483 MultiArray<3, unsigned char> source(shape);
1484 MultiArray<3, TinyVector<float, 3> > dest(shape);
1485 ...
1486 // compute Gaussian gradient at scale sigma
1487 gaussianGradientMultiArray(source, dest, sigma);
1488 \endcode
1489
1490 <b> Usage with anisotropic data:</b>
1491
1492 \code
1493 Shape3 shape(width, height, depth);
1494 MultiArray<3, unsigned char> source(shape);
1495 MultiArray<3, TinyVector<float, 3> > dest(shape);
1496 TinyVector<float, 3> step_size;
1497 TinyVector<float, 3> resolution_sigmas;
1498 ...
1499 // compute Gaussian gradient at scale sigma
1500 gaussianGradientMultiArray(source, dest, sigma,
1501 ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
1502 \endcode
1503
1504 \see separableConvolveMultiArray()
1505 */
1506 doxygen_overloaded_function(template <...> void gaussianGradientMultiArray)
1507
1508 template <class SrcIterator, class SrcShape, class SrcAccessor,
1509 class DestIterator, class DestAccessor>
1510 void
1511 gaussianGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1512 DestIterator di, DestAccessor dest,
1513 ConvolutionOptions<SrcShape::static_size> const & opt,
1514 const char *const function_name = "gaussianGradientMultiArray")
1515 {
1516 typedef typename DestAccessor::value_type DestType;
1517 typedef typename DestType::value_type DestValueType;
1518 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1519
1520 static const int N = SrcShape::static_size;
1521 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
1522
1523 for(int k=0; k<N; ++k)
1524 if(shape[k] <=0)
1525 return;
1526
1527 vigra_precondition(N == (int)dest.size(di),
1528 "gaussianGradientMultiArray(): Wrong number of channels in output array.");
1529
1530 ParamType params = opt.scaleParams();
1531 ParamType params2(params);
1532
1533 ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
1534 for (int dim = 0; dim < N; ++dim, ++params)
1535 {
1536 double sigma = params.sigma_scaled(function_name);
1537 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
1538 }
1539
1540 typedef VectorElementAccessor<DestAccessor> ElementAccessor;
1541
1542 // compute gradient components
1543 for (int dim = 0; dim < N; ++dim, ++params2)
1544 {
1545 ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
1546 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 1, 1.0, opt.window_ratio);
1547 detail::scaleKernel(kernels[dim], 1.0 / params2.step_size());
1548 separableConvolveMultiArray(si, shape, src, di, ElementAccessor(dim, dest), kernels.begin(),
1549 opt.from_point, opt.to_point);
1550 }
1551 }
1552
1553 template <class SrcIterator, class SrcShape, class SrcAccessor,
1554 class DestIterator, class DestAccessor>
1555 void
1556 gaussianGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1557 DestIterator di, DestAccessor dest, double sigma,
1558 ConvolutionOptions<SrcShape::static_size> opt = ConvolutionOptions<SrcShape::static_size>())
1559 {
1560 gaussianGradientMultiArray(si, shape, src, di, dest, opt.stdDev(sigma));
1561 }
1562
1563 template <class SrcIterator, class SrcShape, class SrcAccessor,
1564 class DestIterator, class DestAccessor>
1565 inline void
1566 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1567 pair<DestIterator, DestAccessor> const & dest,
1568 ConvolutionOptions<SrcShape::static_size> const & opt )
1569 {
1570 gaussianGradientMultiArray( source.first, source.second, source.third,
1571 dest.first, dest.second, opt );
1572 }
1573
1574 template <class SrcIterator, class SrcShape, class SrcAccessor,
1575 class DestIterator, class DestAccessor>
1576 inline void
1577 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1578 pair<DestIterator, DestAccessor> const & dest,
1579 double sigma,
1580 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1581 {
1582 gaussianGradientMultiArray( source.first, source.second, source.third,
1583 dest.first, dest.second, sigma, opt );
1584 }
1585
1586 template <unsigned int N, class T1, class S1,
1587 class T2, class S2>
1588 inline void
1589 gaussianGradientMultiArray(MultiArrayView<N, T1, S1> const & source,
1590 MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1591 ConvolutionOptions<N> opt )
1592 {
1593 if(opt.to_point != typename MultiArrayShape<N>::type())
1594 {
1595 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
1596 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
1597 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1598 "gaussianGradientMultiArray(): shape mismatch between ROI and output.");
1599 }
1600 else
1601 {
1602 vigra_precondition(source.shape() == dest.shape(),
1603 "gaussianGradientMultiArray(): shape mismatch between input and output.");
1604 }
1605
1606 gaussianGradientMultiArray( srcMultiArrayRange(source),
1607 destMultiArray(dest), opt );
1608 }
1609
1610 template <unsigned int N, class T1, class S1,
1611 class T2, class S2>
1612 inline void
1613 gaussianGradientMultiArray(MultiArrayView<N, T1, S1> const & source,
1614 MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1615 double sigma,
1616 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
1617 {
1618 gaussianGradientMultiArray( source, dest, opt.stdDev(sigma) );
1619 }
1620
1621 /********************************************************/
1622 /* */
1623 /* gaussianGradientMagnitude */
1624 /* */
1625 /********************************************************/
1626
1627 namespace detail {
1628
1629 template <unsigned int N, class T1, class S1,
1630 class T2, class S2>
1631 void
1632 gaussianGradientMagnitudeImpl(MultiArrayView<N+1, T1, S1> const & src,
1633 MultiArrayView<N, T2, S2> dest,
1634 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
1635 {
1636 typename MultiArrayShape<N>::type shape(src.shape().template subarray<0,N>());
1637 if(opt.to_point != typename MultiArrayShape<N>::type())
1638 {
1639 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.from_point);
1640 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.to_point);
1641 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1642 "gaussianGradientMagnitude(): shape mismatch between ROI and output.");
1643 }
1644 else
1645 {
1646 vigra_precondition(shape == dest.shape(),
1647 "gaussianGradientMagnitude(): shape mismatch between input and output.");
1648 }
1649
1650 dest.init(0.0);
1651
1652 typedef typename NumericTraits<T1>::RealPromote TmpType;
1653 MultiArray<N, TinyVector<TmpType, N> > grad(dest.shape());
1654
1655 using namespace multi_math;
1656
1657 for(int k=0; k<src.shape(N); ++k)
1658 {
1659 gaussianGradientMultiArray(src.bindOuter(k), grad, opt);
1660
1661 dest += squaredNorm(grad);
1662 }
1663 dest = sqrt(dest);
1664 }
1665
1666 } // namespace detail
1667
1668 // documentation is in convolution.hxx
1669 template <unsigned int N, class T1, class S1,
1670 class T2, class S2>
1671 inline void
1672 gaussianGradientMagnitude(MultiArrayView<N+1, Multiband<T1>, S1> const & src,
1673 MultiArrayView<N, T2, S2> dest,
1674 ConvolutionOptions<N> const & opt)
1675 {
1676 detail::gaussianGradientMagnitudeImpl<N, T1>(src, dest, opt);
1677 }
1678
1679 template <unsigned int N, class T1, class S1,
1680 class T2, class S2>
1681 inline void
1682 gaussianGradientMagnitude(MultiArrayView<N, T1, S1> const & src,
1683 MultiArrayView<N, T2, S2> dest,
1684 ConvolutionOptions<N> const & opt)
1685 {
1686 detail::gaussianGradientMagnitudeImpl<N, T1>(src.insertSingletonDimension(N), dest, opt);
1687 }
1688
1689 template <unsigned int N, class T1, int M, class S1,
1690 class T2, class S2>
1691 inline void
1692 gaussianGradientMagnitude(MultiArrayView<N, TinyVector<T1, M>, S1> const & src,
1693 MultiArrayView<N, T2, S2> dest,
1694 ConvolutionOptions<N> const & opt)
1695 {
1696 detail::gaussianGradientMagnitudeImpl<N, T1>(src.expandElements(N), dest, opt);
1697 }
1698
1699 template <unsigned int N, class T1, unsigned int R, unsigned int G, unsigned int B, class S1,
1700 class T2, class S2>
1701 inline void
1702 gaussianGradientMagnitude(MultiArrayView<N, RGBValue<T1, R, G, B>, S1> const & src,
1703 MultiArrayView<N, T2, S2> dest,
1704 ConvolutionOptions<N> const & opt)
1705 {
1706 detail::gaussianGradientMagnitudeImpl<N, T1>(src.expandElements(N), dest, opt);
1707 }
1708
1709 template <unsigned int N, class T1, class S1,
1710 class T2, class S2>
1711 inline void
1712 gaussianGradientMagnitude(MultiArrayView<N, T1, S1> const & src,
1713 MultiArrayView<N, T2, S2> dest,
1714 double sigma,
1715 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
1716 {
1717 gaussianGradientMagnitude(src, dest, opt.stdDev(sigma));
1718 }
1719
1720 template <unsigned int N, class T1, class S1,
1721 class T2, class S2>
1722 inline void
1723 gaussianGradientMagnitude(MultiArrayView<N+1, Multiband<T1>, S1> const & src,
1724 MultiArrayView<N, T2, S2> dest,
1725 double sigma,
1726 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
1727 {
1728 gaussianGradientMagnitude<N>(src, dest, opt.stdDev(sigma));
1729 }
1730
1731 /********************************************************/
1732 /* */
1733 /* symmetricGradientMultiArray */
1734 /* */
1735 /********************************************************/
1736
1737 /** \brief Calculate gradient of a multi-dimensional arrays using symmetric difference filters.
1738
1739 This function computes the gradient of the given N-dimensional
1740 array with a sequence of symmetric difference filters a (differentiation is applied
1741 to each dimension in turn, starting with the innermost dimension).
1742 The destination array is required to have a vector valued pixel type with as many
1743 elements as the number of dimensions. This function is implemented by calls to
1744 \ref convolveMultiArrayOneDimension() with the symmetric difference kernel.
1745
1746 Anisotropic data should be provided with appropriate \ref vigra::ConvolutionOptions
1747 to adjust the filter sizes for the resolution of each axis.
1748 Otherwise, the parameter <tt>opt</tt> is optional unless the parameter
1749 <tt>sigma</tt> is omitted.
1750
1751 If you pass \ref vigra::BlockwiseConvolutionOptions instead, the algorithm will
1752 be executed in parallel on data blocks of a certain size. The block size can be
1753 customized via <tt>BlockwiseConvolutionOptions::blockShape()</tt>, but the defaults
1754 usually work reasonably. By default, the number of threads equals the capabilities
1755 of your hardware, but you can change this via <tt>BlockwiseConvolutionOptions::numThreads()</tt>.
1756
1757 <b> Declarations:</b>
1758
1759 pass arbitrary-dimensional array views:
1760 \code
1761 namespace vigra {
1762 // execute algorithm sequentially
1763 template <unsigned int N, class T1, class S1,
1764 class T2, class S2>
1765 void
1766 symmetricGradientMultiArray(MultiArrayView<N, T1, S1> const & source,
1767 MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1768 ConvolutionOptions<N> opt = ConvolutionOptions<N>());
1769
1770 // execute algorithm in parallel
1771 template <unsigned int N, class T1, class S1,
1772 class T2, class S2>
1773 void
1774 symmetricGradientMultiArray(MultiArrayView<N, T1, S1> const & source,
1775 MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1776 BlockwiseConvolutionOptions<N> opt);
1777 }
1778 \endcode
1779
1780 \deprecatedAPI{symmetricGradientMultiArray}
1781 pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
1782 \code
1783 namespace vigra {
1784 template <class SrcIterator, class SrcShape, class SrcAccessor,
1785 class DestIterator, class DestAccessor>
1786 void
1787 symmetricGradientMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1788 DestIterator diter, DestAccessor dest,
1789 const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1790 }
1791 \endcode
1792 use argument objects in conjunction with \ref ArgumentObjectFactories :
1793 \code
1794 namespace vigra {
1795 template <class SrcIterator, class SrcShape, class SrcAccessor,
1796 class DestIterator, class DestAccessor>
1797 void
1798 symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1799 pair<DestIterator, DestAccessor> const & dest,
1800 const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1801 }
1802 \endcode
1803 \deprecatedEnd
1804
1805 <b> Usage:</b>
1806
1807 <b>\#include</b> <vigra/multi_convolution.hxx> (sequential version)<br/>
1808 <b>\#include</b> <vigra/multi_blockwise.hxx> (parallel version)<br/>
1809 Namespace: vigra
1810
1811 \code
1812 MultiArray<3, unsigned char>::size_type shape(width, height, depth);
1813 MultiArray<3, unsigned char> source(shape);
1814 MultiArray<3, TinyVector<float, 3> > dest(shape);
1815 ...
1816 // compute gradient
1817 symmetricGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest));
1818 \endcode
1819
1820 <b> Usage with anisotropic data:</b>
1821
1822 \code
1823 Shape3 shape(width, height, depth);
1824 MultiArray<3, unsigned char> source(shape);
1825 MultiArray<3, TinyVector<float, 3> > dest(shape);
1826 TinyVector<float, 3> step_size;
1827 ...
1828 // compute gradient
1829 symmetricGradientMultiArray(source, dest,
1830 ConvolutionOptions<3>().stepSize(step_size));
1831 \endcode
1832
1833 \see convolveMultiArrayOneDimension()
1834 */
1835 doxygen_overloaded_function(template <...> void symmetricGradientMultiArray)
1836
1837 template <class SrcIterator, class SrcShape, class SrcAccessor,
1838 class DestIterator, class DestAccessor>
1839 void
1840 symmetricGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
1841 DestIterator di, DestAccessor dest,
1842 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1843 {
1844 typedef typename DestAccessor::value_type DestType;
1845 typedef typename DestType::value_type DestValueType;
1846 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1847
1848 static const int N = SrcShape::static_size;
1849 typedef typename ConvolutionOptions<N>::StepIterator StepType;
1850
1851 for(int k=0; k<N; ++k)
1852 if(shape[k] <=0)
1853 return;
1854
1855 vigra_precondition(N == (int)dest.size(di),
1856 "symmetricGradientMultiArray(): Wrong number of channels in output array.");
1857
1858 Kernel1D<KernelType> filter;
1859 filter.initSymmetricDifference();
1860
1861 StepType step_size_it = opt.stepParams();
1862
1863 typedef VectorElementAccessor<DestAccessor> ElementAccessor;
1864
1865 // compute gradient components
1866 for (int d = 0; d < N; ++d, ++step_size_it)
1867 {
1868 Kernel1D<KernelType> symmetric(filter);
1869 detail::scaleKernel(symmetric, 1 / *step_size_it);
1870 convolveMultiArrayOneDimension(si, shape, src,
1871 di, ElementAccessor(d, dest),
1872 d, symmetric, opt.from_point, opt.to_point);
1873 }
1874 }
1875
1876 template <class SrcIterator, class SrcShape, class SrcAccessor,
1877 class DestIterator, class DestAccessor>
1878 inline void
1879 symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1880 pair<DestIterator, DestAccessor> const & dest,
1881 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1882 {
1883 symmetricGradientMultiArray(source.first, source.second, source.third,
1884 dest.first, dest.second, opt);
1885 }
1886
1887 template <unsigned int N, class T1, class S1,
1888 class T2, class S2>
1889 inline void
1890 symmetricGradientMultiArray(MultiArrayView<N, T1, S1> const & source,
1891 MultiArrayView<N, TinyVector<T2, N>, S2> dest,
1892 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
1893 {
1894 if(opt.to_point != typename MultiArrayShape<N>::type())
1895 {
1896 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
1897 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
1898 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1899 "symmetricGradientMultiArray(): shape mismatch between ROI and output.");
1900 }
1901 else
1902 {
1903 vigra_precondition(source.shape() == dest.shape(),
1904 "symmetricGradientMultiArray(): shape mismatch between input and output.");
1905 }
1906
1907 symmetricGradientMultiArray(srcMultiArrayRange(source),
1908 destMultiArray(dest), opt);
1909 }
1910
1911 /********************************************************/
1912 /* */
1913 /* laplacianOfGaussianMultiArray */
1914 /* */
1915 /********************************************************/
1916
1917 /** \brief Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters.
1918
1919 This function computes the Laplacian of the given N-dimensional
1920 array with a sequence of second-derivative-of-Gaussian filters at the given
1921 standard deviation <tt>sigma</tt>. Both source and destination
1922 arrays must have scalar value_type. This function is implemented by calls to
1923 \ref separableConvolveMultiArray() with the appropriate kernels, followed by summation.
1924
1925 Anisotropic data should be provided with appropriate \ref vigra::ConvolutionOptions
1926 to adjust the filter sizes for the resolution of each axis.
1927 Otherwise, the parameter <tt>opt</tt> is optional unless the parameter
1928 <tt>sigma</tt> is omitted.
1929
1930 If you pass \ref vigra::BlockwiseConvolutionOptions instead, the algorithm will
1931 be executed in parallel on data blocks of a certain size. The block size can be
1932 customized via <tt>BlockwiseConvolutionOptions::blockShape()</tt>, but the defaults
1933 usually work reasonably. By default, the number of threads equals the capabilities
1934 of your hardware, but you can change this via <tt>BlockwiseConvolutionOptions::numThreads()</tt>.
1935
1936 <b> Declarations:</b>
1937
1938 pass arbitrary-dimensional array views:
1939 \code
1940 namespace vigra {
1941 // pass scale explicitly
1942 template <unsigned int N, class T1, class S1,
1943 class T2, class S2>
1944 void
1945 laplacianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
1946 MultiArrayView<N, T2, S2> dest,
1947 double sigma,
1948 ConvolutionOptions<N> opt = ConvolutionOptions<N>());
1949
1950 // pass scale(s) in option object
1951 template <unsigned int N, class T1, class S1,
1952 class T2, class S2>
1953 void
1954 laplacianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
1955 MultiArrayView<N, T2, S2> dest,
1956 ConvolutionOptions<N> opt );
1957
1958 // likewise, but execute algorithm in parallel
1959 template <unsigned int N, class T1, class S1,
1960 class T2, class S2>
1961 void
1962 laplacianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
1963 MultiArrayView<N, T2, S2> dest,
1964 BlockwiseConvolutionOptions<N> opt );
1965 }
1966 \endcode
1967
1968 \deprecatedAPI{laplacianOfGaussianMultiArray}
1969 pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
1970 \code
1971 namespace vigra {
1972 template <class SrcIterator, class SrcShape, class SrcAccessor,
1973 class DestIterator, class DestAccessor>
1974 void
1975 laplacianOfGaussianMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
1976 DestIterator diter, DestAccessor dest,
1977 double sigma,
1978 const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1979 }
1980 \endcode
1981 use argument objects in conjunction with \ref ArgumentObjectFactories :
1982 \code
1983 namespace vigra {
1984 template <class SrcIterator, class SrcShape, class SrcAccessor,
1985 class DestIterator, class DestAccessor>
1986 void
1987 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
1988 pair<DestIterator, DestAccessor> const & dest,
1989 double sigma,
1990 const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
1991 }
1992 \endcode
1993 \deprecatedEnd
1994
1995 <b> Usage:</b>
1996
1997 <b>\#include</b> <vigra/multi_convolution.hxx> (sequential version)<br/>
1998 <b>\#include</b> <vigra/multi_blockwise.hxx> (parallel version)<br/>
1999 Namespace: vigra
2000
2001 \code
2002 Shape3 shape(width, height, depth);
2003 MultiArray<3, float> source(shape);
2004 MultiArray<3, float> laplacian(shape);
2005 ...
2006 // compute Laplacian at scale sigma
2007 laplacianOfGaussianMultiArray(source, laplacian, sigma);
2008 \endcode
2009
2010 <b> Usage with anisotropic data:</b>
2011
2012 \code
2013 MultiArray<3, float> source(shape);
2014 MultiArray<3, float> laplacian(shape);
2015 TinyVector<float, 3> step_size;
2016 TinyVector<float, 3> resolution_sigmas;
2017 ...
2018 // compute Laplacian at scale sigma
2019 laplacianOfGaussianMultiArray(source, laplacian, sigma,
2020 ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
2021 \endcode
2022
2023 \see separableConvolveMultiArray()
2024 */
2025 doxygen_overloaded_function(template <...> void laplacianOfGaussianMultiArray)
2026
2027 template <class SrcIterator, class SrcShape, class SrcAccessor,
2028 class DestIterator, class DestAccessor>
2029 void
2030 laplacianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
2031 DestIterator di, DestAccessor dest,
2032 ConvolutionOptions<SrcShape::static_size> const & opt )
2033 {
2034 using namespace functor;
2035
2036 typedef typename DestAccessor::value_type DestType;
2037 typedef typename NumericTraits<DestType>::RealPromote KernelType;
2038 typedef typename AccessorTraits<KernelType>::default_accessor DerivativeAccessor;
2039
2040 static const int N = SrcShape::static_size;
2041 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
2042
2043 ParamType params = opt.scaleParams();
2044 ParamType params2(params);
2045
2046 ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
2047 for (int dim = 0; dim < N; ++dim, ++params)
2048 {
2049 double sigma = params.sigma_scaled("laplacianOfGaussianMultiArray");
2050 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
2051 }
2052
2053 SrcShape dshape(shape);
2054 if(opt.to_point != SrcShape())
2055 dshape = opt.to_point - opt.from_point;
2056
2057 MultiArray<N, KernelType> derivative(dshape);
2058
2059 // compute 2nd derivatives and sum them up
2060 for (int dim = 0; dim < N; ++dim, ++params2)
2061 {
2062 ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
2063 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 2, 1.0, opt.window_ratio);
2064 detail::scaleKernel(kernels[dim], 1.0 / sq(params2.step_size()));
2065
2066 if (dim == 0)
2067 {
2068 separableConvolveMultiArray( si, shape, src,
2069 di, dest, kernels.begin(), opt.from_point, opt.to_point);
2070 }
2071 else
2072 {
2073 separableConvolveMultiArray( si, shape, src,
2074 derivative.traverser_begin(), DerivativeAccessor(),
2075 kernels.begin(), opt.from_point, opt.to_point);
2076 combineTwoMultiArrays(di, dshape, dest, derivative.traverser_begin(), DerivativeAccessor(),
2077 di, dest, Arg1() + Arg2() );
2078 }
2079 }
2080 }
2081
2082 template <class SrcIterator, class SrcShape, class SrcAccessor,
2083 class DestIterator, class DestAccessor>
2084 void
2085 laplacianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
2086 DestIterator di, DestAccessor dest, double sigma,
2087 ConvolutionOptions<SrcShape::static_size> opt = ConvolutionOptions<SrcShape::static_size>())
2088 {
2089 laplacianOfGaussianMultiArray(si, shape, src, di, dest, opt.stdDev(sigma));
2090 }
2091
2092 template <class SrcIterator, class SrcShape, class SrcAccessor,
2093 class DestIterator, class DestAccessor>
2094 inline void
2095 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2096 pair<DestIterator, DestAccessor> const & dest,
2097 ConvolutionOptions<SrcShape::static_size> const & opt )
2098 {
2099 laplacianOfGaussianMultiArray( source.first, source.second, source.third,
2100 dest.first, dest.second, opt );
2101 }
2102
2103 template <class SrcIterator, class SrcShape, class SrcAccessor,
2104 class DestIterator, class DestAccessor>
2105 inline void
2106 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2107 pair<DestIterator, DestAccessor> const & dest,
2108 double sigma,
2109 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
2110 {
2111 laplacianOfGaussianMultiArray( source.first, source.second, source.third,
2112 dest.first, dest.second, sigma, opt );
2113 }
2114
2115 template <unsigned int N, class T1, class S1,
2116 class T2, class S2>
2117 inline void
2118 laplacianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
2119 MultiArrayView<N, T2, S2> dest,
2120 ConvolutionOptions<N> opt )
2121 {
2122 if(opt.to_point != typename MultiArrayShape<N>::type())
2123 {
2124 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
2125 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
2126 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2127 "laplacianOfGaussianMultiArray(): shape mismatch between ROI and output.");
2128 }
2129 else
2130 {
2131 vigra_precondition(source.shape() == dest.shape(),
2132 "laplacianOfGaussianMultiArray(): shape mismatch between input and output.");
2133 }
2134
2135 laplacianOfGaussianMultiArray( srcMultiArrayRange(source),
2136 destMultiArray(dest), opt );
2137 }
2138
2139 template <unsigned int N, class T1, class S1,
2140 class T2, class S2>
2141 inline void
2142 laplacianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
2143 MultiArrayView<N, T2, S2> dest,
2144 double sigma,
2145 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
2146 {
2147 laplacianOfGaussianMultiArray( source, dest, opt.stdDev(sigma) );
2148 }
2149
2150 /********************************************************/
2151 /* */
2152 /* gaussianDivergenceMultiArray */
2153 /* */
2154 /********************************************************/
2155
2156 /** \brief Calculate the divergence of a vector field using Gaussian derivative filters.
2157
2158 This function computes the divergence of the given N-dimensional vector field
2159 with a sequence of first-derivative-of-Gaussian filters at the given
2160 standard deviation <tt>sigma</tt>. The input vector field can either be given as a sequence
2161 of scalar array views (one for each vector field component), represented by an
2162 iterator range, or by a single vector array with the appropriate shape.
2163 This function is implemented by calls to
2164 \ref separableConvolveMultiArray() with the suitable kernels, followed by summation.
2165
2166 Anisotropic data should be provided with appropriate \ref vigra::ConvolutionOptions
2167 to adjust the filter sizes for the resolution of each axis.
2168 Otherwise, the parameter <tt>opt</tt> is optional unless the parameter
2169 <tt>sigma</tt> is omitted.
2170
2171 If you pass \ref vigra::BlockwiseConvolutionOptions instead, the algorithm will
2172 be executed in parallel on data blocks of a certain size. The block size can be
2173 customized via <tt>BlockwiseConvolutionOptions::blockShape()</tt>, but the defaults
2174 usually work reasonably. By default, the number of threads equals the capabilities
2175 of your hardware, but you can change this via <tt>BlockwiseConvolutionOptions::numThreads()</tt>.
2176
2177 <b> Declarations:</b>
2178
2179 pass arbitrary-dimensional array views:
2180 \code
2181 namespace vigra {
2182 // specify input vector field as a sequence of scalar arrays
2183 template <class Iterator,
2184 unsigned int N, class T, class S>
2185 void
2186 gaussianDivergenceMultiArray(Iterator vectorField, Iterator vectorFieldEnd,
2187 MultiArrayView<N, T, S> divergence,
2188 ConvolutionOptions<N> const & opt);
2189
2190 template <class Iterator,
2191 unsigned int N, class T, class S>
2192 void
2193 gaussianDivergenceMultiArray(Iterator vectorField, Iterator vectorFieldEnd,
2194 MultiArrayView<N, T, S> divergence,
2195 double sigma,
2196 ConvolutionOptions<N> opt = ConvolutionOptions<N>());
2197
2198 // pass input vector field as an array of vectors
2199 template <unsigned int N, class T1, class S1,
2200 class T2, class S2>
2201 void
2202 gaussianDivergenceMultiArray(MultiArrayView<N, TinyVector<T1, N>, S1> const & vectorField,
2203 MultiArrayView<N, T2, S2> divergence,
2204 ConvolutionOptions<N> const & opt);
2205
2206 template <unsigned int N, class T1, class S1,
2207 class T2, class S2>
2208 void
2209 gaussianDivergenceMultiArray(MultiArrayView<N, TinyVector<T1, N>, S1> const & vectorField,
2210 MultiArrayView<N, T2, S2> divergence,
2211 double sigma,
2212 ConvolutionOptions<N> opt = ConvolutionOptions<N>());
2213
2214 // pass input vector field as an array of vectors and
2215 // execute algorithm in parallel
2216 template <unsigned int N, class T1, class S1,
2217 class T2, class S2>
2218 void
2219 gaussianDivergenceMultiArray(MultiArrayView<N, TinyVector<T1, N>, S1> const & vectorField,
2220 MultiArrayView<N, T2, S2> divergence,
2221 BlockwiseConvolutionOptions<N> const & opt);
2222 }
2223 \endcode
2224
2225 <b> Usage:</b>
2226
2227 <b>\#include</b> <vigra/multi_convolution.hxx> (sequential version)<br/>
2228 <b>\#include</b> <vigra/multi_blockwise.hxx> (parallel version)<br/>
2229 Namespace: vigra
2230
2231 \code
2232 Shape3 shape(width, height, depth);
2233 MultiArray<3, TinyVector<float, 3> > source(shape);
2234 MultiArray<3, float> laplacian(shape);
2235 ...
2236 // compute divergence at scale sigma
2237 gaussianDivergenceMultiArray(source, laplacian, sigma);
2238 \endcode
2239
2240 <b> Usage with anisotropic data:</b>
2241
2242 \code
2243 MultiArray<3, TinyVector<float, 3> > source(shape);
2244 MultiArray<3, float> laplacian(shape);
2245 TinyVector<float, 3> step_size;
2246 TinyVector<float, 3> resolution_sigmas;
2247 ...
2248 // compute divergence at scale sigma
2249 gaussianDivergenceMultiArray(source, laplacian, sigma,
2250 ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
2251 \endcode
2252 */
2253 doxygen_overloaded_function(template <...> void gaussianDivergenceMultiArray)
2254
2255 template <class Iterator,
2256 unsigned int N, class T, class S>
2257 void
2258 gaussianDivergenceMultiArray(Iterator vectorField, Iterator vectorFieldEnd,
2259 MultiArrayView<N, T, S> divergence,
2260 ConvolutionOptions<N> opt)
2261 {
2262 typedef typename std::iterator_traits<Iterator>::value_type ArrayType;
2263 typedef typename ArrayType::value_type SrcType;
2264 typedef typename NumericTraits<SrcType>::RealPromote TmpType;
2265 typedef Kernel1D<double> Kernel;
2266
2267 vigra_precondition(std::distance(vectorField, vectorFieldEnd) == N,
2268 "gaussianDivergenceMultiArray(): wrong number of input arrays.");
2269 // more checks are performed in separableConvolveMultiArray()
2270
2271 typename ConvolutionOptions<N>::ScaleIterator params = opt.scaleParams();
2272 ArrayVector<double> sigmas(N);
2273 ArrayVector<Kernel> kernels(N);
2274 for(unsigned int k = 0; k < N; ++k, ++params)
2275 {
2276 sigmas[k] = params.sigma_scaled("gaussianDivergenceMultiArray");
2277 kernels[k].initGaussian(sigmas[k], 1.0, opt.window_ratio);
2278 }
2279
2280 MultiArray<N, TmpType> tmpDeriv(divergence.shape());
2281
2282 for(unsigned int k=0; k < N; ++k, ++vectorField)
2283 {
2284 kernels[k].initGaussianDerivative(sigmas[k], 1, 1.0, opt.window_ratio);
2285 if(k == 0)
2286 {
2287 separableConvolveMultiArray(*vectorField, divergence, kernels.begin(), opt.from_point, opt.to_point);
2288 }
2289 else
2290 {
2291 separableConvolveMultiArray(*vectorField, tmpDeriv, kernels.begin(), opt.from_point, opt.to_point);
2292 divergence += tmpDeriv;
2293 }
2294 kernels[k].initGaussian(sigmas[k], 1.0, opt.window_ratio);
2295 }
2296 }
2297
2298 template <class Iterator,
2299 unsigned int N, class T, class S>
2300 inline void
2301 gaussianDivergenceMultiArray(Iterator vectorField, Iterator vectorFieldEnd,
2302 MultiArrayView<N, T, S> divergence,
2303 double sigma,
2304 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
2305 {
2306 gaussianDivergenceMultiArray(vectorField, vectorFieldEnd, divergence, opt.stdDev(sigma));
2307 }
2308
2309 template <unsigned int N, class T1, class S1,
2310 class T2, class S2>
2311 inline void
2312 gaussianDivergenceMultiArray(MultiArrayView<N, TinyVector<T1, N>, S1> const & vectorField,
2313 MultiArrayView<N, T2, S2> divergence,
2314 ConvolutionOptions<N> const & opt)
2315 {
2316 ArrayVector<MultiArrayView<N, T1> > field;
2317 for(unsigned int k=0; k<N; ++k)
2318 field.push_back(vectorField.bindElementChannel(k));
2319
2320 gaussianDivergenceMultiArray(field.begin(), field.end(), divergence, opt);
2321 }
2322
2323 template <unsigned int N, class T1, class S1,
2324 class T2, class S2>
2325 inline void
2326 gaussianDivergenceMultiArray(MultiArrayView<N, TinyVector<T1, N>, S1> const & vectorField,
2327 MultiArrayView<N, T2, S2> divergence,
2328 double sigma,
2329 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
2330 {
2331 gaussianDivergenceMultiArray(vectorField, divergence, opt.stdDev(sigma));
2332 }
2333
2334 /********************************************************/
2335 /* */
2336 /* hessianOfGaussianMultiArray */
2337 /* */
2338 /********************************************************/
2339
2340 /** \brief Calculate Hessian matrix of a N-dimensional arrays using Gaussian derivative filters.
2341
2342 This function computes the Hessian matrix the given scalar N-dimensional
2343 array with a sequence of second-derivative-of-Gaussian filters at the given
2344 standard deviation <tt>sigma</tt>. The destination array must
2345 have a vector valued element type with N*(N+1)/2 elements (it represents the
2346 upper triangular part of the symmetric Hessian matrix, flattened row-wise).
2347 This function is implemented by calls to
2348 \ref separableConvolveMultiArray() with the appropriate kernels.
2349
2350 Anisotropic data should be provided with appropriate \ref vigra::ConvolutionOptions
2351 to adjust the filter sizes for the resolution of each axis.
2352 Otherwise, the parameter <tt>opt</tt> is optional unless the parameter
2353 <tt>sigma</tt> is omitted.
2354
2355 If you pass \ref vigra::BlockwiseConvolutionOptions instead, the algorithm will
2356 be executed in parallel on data blocks of a certain size. The block size can be
2357 customized via <tt>BlockwiseConvolutionOptions::blockShape()</tt>, but the defaults
2358 usually work reasonably. By default, the number of threads equals the capabilities
2359 of your hardware, but you can change this via <tt>BlockwiseConvolutionOptions::numThreads()</tt>.
2360
2361 <b> Declarations:</b>
2362
2363 pass arbitrary-dimensional array views:
2364 \code
2365 namespace vigra {
2366 // pass scale explicitly
2367 template <unsigned int N, class T1, class S1,
2368 class T2, class S2>
2369 void
2370 hessianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
2371 MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2372 double sigma,
2373 ConvolutionOptions<N> opt = ConvolutionOptions<N>());
2374
2375 // pass scale(s) in option object
2376 template <unsigned int N, class T1, class S1,
2377 class T2, class S2>
2378 void
2379 hessianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
2380 MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2381 ConvolutionOptions<N> opt);
2382
2383 // likewise, but execute algorithm in parallel
2384 template <unsigned int N, class T1, class S1,
2385 class T2, class S2>
2386 void
2387 hessianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
2388 MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2389 BlockwiseConvolutionOptions<N> opt);
2390 }
2391 \endcode
2392
2393 \deprecatedAPI{hessianOfGaussianMultiArray}
2394 pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
2395 \code
2396 namespace vigra {
2397 template <class SrcIterator, class SrcShape, class SrcAccessor,
2398 class DestIterator, class DestAccessor>
2399 void
2400 hessianOfGaussianMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
2401 DestIterator diter, DestAccessor dest,
2402 double sigma,
2403 const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
2404 }
2405 \endcode
2406 use argument objects in conjunction with \ref ArgumentObjectFactories :
2407 \code
2408 namespace vigra {
2409 template <class SrcIterator, class SrcShape, class SrcAccessor,
2410 class DestIterator, class DestAccessor>
2411 void
2412 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2413 pair<DestIterator, DestAccessor> const & dest,
2414 double sigma,
2415 const ConvolutionOptions<N> & opt = ConvolutionOptions<N>());
2416 }
2417 \endcode
2418 \deprecatedEnd
2419
2420 <b> Usage:</b>
2421
2422 <b>\#include</b> <vigra/multi_convolution.hxx> (sequential version)<br/>
2423 <b>\#include</b> <vigra/multi_blockwise.hxx> (parallel version)<br/>
2424 Namespace: vigra
2425
2426 \code
2427 Shape3 shape(width, height, depth);
2428 MultiArray<3, float> source(shape);
2429 MultiArray<3, TinyVector<float, 6> > dest(shape);
2430 ...
2431 // compute Hessian at scale sigma
2432 hessianOfGaussianMultiArray(source, dest, sigma);
2433 \endcode
2434
2435 <b> Usage with anisotropic data:</b>
2436
2437 \code
2438 MultiArray<3, float> source(shape);
2439 MultiArray<3, TinyVector<float, 6> > dest(shape);
2440 TinyVector<float, 3> step_size;
2441 TinyVector<float, 3> resolution_sigmas;
2442 ...
2443 // compute Hessian at scale sigma
2444 hessianOfGaussianMultiArray(source, dest, sigma,
2445 ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
2446 \endcode
2447
2448 \see separableConvolveMultiArray(), vectorToTensorMultiArray()
2449 */
2450 doxygen_overloaded_function(template <...> void hessianOfGaussianMultiArray)
2451
2452 template <class SrcIterator, class SrcShape, class SrcAccessor,
2453 class DestIterator, class DestAccessor>
2454 void
2455 hessianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
2456 DestIterator di, DestAccessor dest,
2457 ConvolutionOptions<SrcShape::static_size> const & opt )
2458 {
2459 typedef typename DestAccessor::value_type DestType;
2460 typedef typename DestType::value_type DestValueType;
2461 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
2462
2463 static const int N = SrcShape::static_size;
2464 static const int M = N*(N+1)/2;
2465 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
2466
2467 for(int k=0; k<N; ++k)
2468 if(shape[k] <=0)
2469 return;
2470
2471 vigra_precondition(M == (int)dest.size(di),
2472 "hessianOfGaussianMultiArray(): Wrong number of channels in output array.");
2473
2474 ParamType params_init = opt.scaleParams();
2475
2476 ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
2477 ParamType params(params_init);
2478 for (int dim = 0; dim < N; ++dim, ++params)
2479 {
2480 double sigma = params.sigma_scaled("hessianOfGaussianMultiArray");
2481 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
2482 }
2483
2484 typedef VectorElementAccessor<DestAccessor> ElementAccessor;
2485
2486 // compute elements of the Hessian matrix
2487 ParamType params_i(params_init);
2488 for (int b=0, i=0; i<N; ++i, ++params_i)
2489 {
2490 ParamType params_j(params_i);
2491 for (int j=i; j<N; ++j, ++b, ++params_j)
2492 {
2493 ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
2494 if(i == j)
2495 {
2496 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 2, 1.0, opt.window_ratio);
2497 }
2498 else
2499 {
2500 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 1, 1.0, opt.window_ratio);
2501 kernels[j].initGaussianDerivative(params_j.sigma_scaled(), 1, 1.0, opt.window_ratio);
2502 }
2503 detail::scaleKernel(kernels[i], 1 / params_i.step_size());
2504 detail::scaleKernel(kernels[j], 1 / params_j.step_size());
2505 separableConvolveMultiArray(si, shape, src, di, ElementAccessor(b, dest),
2506 kernels.begin(), opt.from_point, opt.to_point);
2507 }
2508 }
2509 }
2510
2511 template <class SrcIterator, class SrcShape, class SrcAccessor,
2512 class DestIterator, class DestAccessor>
2513 inline void
2514 hessianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
2515 DestIterator di, DestAccessor dest, double sigma,
2516 ConvolutionOptions<SrcShape::static_size> opt = ConvolutionOptions<SrcShape::static_size>())
2517 {
2518 hessianOfGaussianMultiArray(si, shape, src, di, dest, opt.stdDev(sigma));
2519 }
2520
2521 template <class SrcIterator, class SrcShape, class SrcAccessor,
2522 class DestIterator, class DestAccessor>
2523 inline void
2524 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2525 pair<DestIterator, DestAccessor> const & dest,
2526 ConvolutionOptions<SrcShape::static_size> const & opt )
2527 {
2528 hessianOfGaussianMultiArray( source.first, source.second, source.third,
2529 dest.first, dest.second, opt );
2530 }
2531
2532 template <class SrcIterator, class SrcShape, class SrcAccessor,
2533 class DestIterator, class DestAccessor>
2534 inline void
2535 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2536 pair<DestIterator, DestAccessor> const & dest,
2537 double sigma,
2538 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
2539 {
2540 hessianOfGaussianMultiArray( source.first, source.second, source.third,
2541 dest.first, dest.second, sigma, opt );
2542 }
2543
2544 template <unsigned int N, class T1, class S1,
2545 class T2, class S2>
2546 inline void
2547 hessianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
2548 MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2549 ConvolutionOptions<N> opt )
2550 {
2551 if(opt.to_point != typename MultiArrayShape<N>::type())
2552 {
2553 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
2554 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
2555 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2556 "hessianOfGaussianMultiArray(): shape mismatch between ROI and output.");
2557 }
2558 else
2559 {
2560 vigra_precondition(source.shape() == dest.shape(),
2561 "hessianOfGaussianMultiArray(): shape mismatch between input and output.");
2562 }
2563
2564 hessianOfGaussianMultiArray( srcMultiArrayRange(source),
2565 destMultiArray(dest), opt );
2566 }
2567
2568 template <unsigned int N, class T1, class S1,
2569 class T2, class S2>
2570 inline void
2571 hessianOfGaussianMultiArray(MultiArrayView<N, T1, S1> const & source,
2572 MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2573 double sigma,
2574 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
2575 {
2576 hessianOfGaussianMultiArray( source, dest, opt.stdDev(sigma) );
2577 }
2578
2579 namespace detail {
2580
2581 template<int N, class VectorType>
2582 struct StructurTensorFunctor
2583 {
2584 typedef VectorType result_type;
2585 typedef typename VectorType::value_type ValueType;
2586
2587 template <class T>
2588 VectorType operator()(T const & in) const
2589 {
2590 VectorType res;
2591 for(int b=0, i=0; i<N; ++i)
2592 {
2593 for(int j=i; j<N; ++j, ++b)
2594 {
2595 res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
2596 }
2597 }
2598 return res;
2599 }
2600 };
2601
2602 } // namespace detail
2603
2604 /********************************************************/
2605 /* */
2606 /* structureTensorMultiArray */
2607 /* */
2608 /********************************************************/
2609
2610 /** \brief Calculate th structure tensor of a multi-dimensional arrays.
2611
2612 This function computes the gradient (outer product) tensor for each element
2613 of the given N-dimensional array with first-derivative-of-Gaussian filters at
2614 the given <tt>innerScale</tt>, followed by Gaussian smoothing at <tt>outerScale</tt>.
2615 The destination array must have a vector valued pixel type with
2616 N*(N+1)/2 elements (it represents the upper triangular part of the symmetric
2617 structure tensor matrix, flattened row-wise). If the source array is also vector valued, the
2618 resulting structure tensor is the sum of the individual tensors for each channel.
2619 This function is implemented by calls to
2620 \ref separableConvolveMultiArray() with the appropriate kernels.
2621
2622 Anisotropic data should be provided with appropriate \ref vigra::ConvolutionOptions
2623 to adjust the filter sizes for the resolution of each axis.
2624 Otherwise, the parameter <tt>opt</tt> is optional unless the parameters
2625 <tt>innerScale</tt> and <tt>outerScale</tt> are both omitted.
2626
2627 If you pass \ref vigra::BlockwiseConvolutionOptions instead, the algorithm will
2628 be executed in parallel on data blocks of a certain size. The block size can be
2629 customized via <tt>BlockwiseConvolutionOptions::blockShape()</tt>, but the defaults
2630 usually work reasonably. By default, the number of threads equals the capabilities
2631 of your hardware, but you can change this via <tt>BlockwiseConvolutionOptions::numThreads()</tt>.
2632
2633 <b> Declarations:</b>
2634
2635 pass arbitrary-dimensional array views:
2636 \code
2637 namespace vigra {
2638 // pass scales explicitly
2639 template <unsigned int N, class T1, class S1,
2640 class T2, class S2>
2641 void
2642 structureTensorMultiArray(MultiArrayView<N, T1, S1> const & source,
2643 MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2644 double innerScale, double outerScale,
2645 ConvolutionOptions<N> opt = ConvolutionOptions<N>());
2646
2647 // pass scales in option object
2648 template <unsigned int N, class T1, class S1,
2649 class T2, class S2>
2650 void
2651 structureTensorMultiArray(MultiArrayView<N, T1, S1> const & source,
2652 MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2653 ConvolutionOptions<N> opt );
2654
2655 // likewise, but execute algorithm in parallel
2656 template <unsigned int N, class T1, class S1,
2657 class T2, class S2>
2658 void
2659 structureTensorMultiArray(MultiArrayView<N, T1, S1> const & source,
2660 MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2661 BlockwiseConvolutionOptions<N> opt );
2662 }
2663 \endcode
2664
2665 \deprecatedAPI{structureTensorMultiArray}
2666 pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
2667 \code
2668 namespace vigra {
2669 template <class SrcIterator, class SrcShape, class SrcAccessor,
2670 class DestIterator, class DestAccessor>
2671 void
2672 structureTensorMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
2673 DestIterator diter, DestAccessor dest,
2674 double innerScale, double outerScale,
2675 ConvolutionOptions<N> opt);
2676 }
2677 \endcode
2678 use argument objects in conjunction with \ref ArgumentObjectFactories :
2679 \code
2680 namespace vigra {
2681 template <class SrcIterator, class SrcShape, class SrcAccessor,
2682 class DestIterator, class DestAccessor>
2683 void
2684 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2685 pair<DestIterator, DestAccessor> const & dest,
2686 double innerScale, double outerScale,
2687 const ConvolutionOptions<N> & opt);
2688 }
2689 \endcode
2690 \deprecatedEnd
2691
2692 <b> Usage:</b>
2693
2694 <b>\#include</b> <vigra/multi_convolution.hxx> (sequential version)<br/>
2695 <b>\#include</b> <vigra/multi_blockwise.hxx> (parallel version)<br/>
2696 Namespace: vigra
2697
2698 \code
2699 Shape3 shape(width, height, depth);
2700 MultiArray<3, RGBValue<float> > source(shape);
2701 MultiArray<3, TinyVector<float, 6> > dest(shape);
2702 ...
2703 // compute structure tensor at scales innerScale and outerScale
2704 structureTensorMultiArray(source, dest, innerScale, outerScale);
2705 \endcode
2706
2707 <b> Usage with anisotropic data:</b>
2708
2709 \code
2710 MultiArray<3, RGBValue<float> > source(shape);
2711 MultiArray<3, TinyVector<float, 6> > dest(shape);
2712 TinyVector<float, 3> step_size;
2713 TinyVector<float, 3> resolution_sigmas;
2714 ...
2715 // compute structure tensor at scales innerScale and outerScale
2716 structureTensorMultiArray(source, dest, innerScale, outerScale,
2717 ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas));
2718 \endcode
2719
2720 \see separableConvolveMultiArray(), vectorToTensorMultiArray()
2721 */
2722 doxygen_overloaded_function(template <...> void structureTensorMultiArray)
2723
2724 template <class SrcIterator, class SrcShape, class SrcAccessor,
2725 class DestIterator, class DestAccessor>
2726 void
2727 structureTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
2728 DestIterator di, DestAccessor dest,
2729 ConvolutionOptions<SrcShape::static_size> opt)
2730 {
2731 static const int N = SrcShape::static_size;
2732 static const int M = N*(N+1)/2;
2733
2734 typedef typename DestAccessor::value_type DestType;
2735 typedef typename DestType::value_type DestValueType;
2736 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
2737 typedef TinyVector<KernelType, N> GradientVector;
2738 typedef typename AccessorTraits<GradientVector>::default_accessor GradientAccessor;
2739 typedef typename AccessorTraits<DestType>::default_accessor GradientTensorAccessor;
2740
2741 for(int k=0; k<N; ++k)
2742 if(shape[k] <=0)
2743 return;
2744
2745 vigra_precondition(M == (int)dest.size(di),
2746 "structureTensorMultiArray(): Wrong number of channels in output array.");
2747
2748 ConvolutionOptions<N> innerOptions = opt;
2749 ConvolutionOptions<N> outerOptions = opt.outerOptions();
2750 typename ConvolutionOptions<N>::ScaleIterator params = outerOptions.scaleParams();
2751
2752 SrcShape gradientShape(shape);
2753 if(opt.to_point != SrcShape())
2754 {
2755 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.from_point);
2756 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.to_point);
2757
2758 for(int k=0; k<N; ++k, ++params)
2759 {
2760 Kernel1D<double> gauss;
2761 gauss.initGaussian(params.sigma_scaled("structureTensorMultiArray"), 1.0, opt.window_ratio);
2762 int dilation = gauss.right();
2763 innerOptions.from_point[k] = std::max<MultiArrayIndex>(0, opt.from_point[k] - dilation);
2764 innerOptions.to_point[k] = std::min<MultiArrayIndex>(shape[k], opt.to_point[k] + dilation);
2765 }
2766 outerOptions.from_point -= innerOptions.from_point;
2767 outerOptions.to_point -= innerOptions.from_point;
2768 gradientShape = innerOptions.to_point - innerOptions.from_point;
2769 }
2770
2771 MultiArray<N, GradientVector> gradient(gradientShape);
2772 MultiArray<N, DestType> gradientTensor(gradientShape);
2773 gaussianGradientMultiArray(si, shape, src,
2774 gradient.traverser_begin(), GradientAccessor(),
2775 innerOptions,
2776 "structureTensorMultiArray");
2777
2778 transformMultiArray(gradient.traverser_begin(), gradientShape, GradientAccessor(),
2779 gradientTensor.traverser_begin(), GradientTensorAccessor(),
2780 detail::StructurTensorFunctor<N, DestType>());
2781
2782 gaussianSmoothMultiArray(gradientTensor.traverser_begin(), gradientShape, GradientTensorAccessor(),
2783 di, dest, outerOptions,
2784 "structureTensorMultiArray");
2785 }
2786
2787 template <class SrcIterator, class SrcShape, class SrcAccessor,
2788 class DestIterator, class DestAccessor>
2789 inline void
2790 structureTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
2791 DestIterator di, DestAccessor dest,
2792 double innerScale, double outerScale,
2793 ConvolutionOptions<SrcShape::static_size> opt = ConvolutionOptions<SrcShape::static_size>())
2794 {
2795 structureTensorMultiArray(si, shape, src, di, dest,
2796 opt.stdDev(innerScale).outerScale(outerScale));
2797 }
2798
2799 template <class SrcIterator, class SrcShape, class SrcAccessor,
2800 class DestIterator, class DestAccessor>
2801 inline void
2802 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2803 pair<DestIterator, DestAccessor> const & dest,
2804 ConvolutionOptions<SrcShape::static_size> const & opt )
2805 {
2806 structureTensorMultiArray( source.first, source.second, source.third,
2807 dest.first, dest.second, opt );
2808 }
2809
2810
2811 template <class SrcIterator, class SrcShape, class SrcAccessor,
2812 class DestIterator, class DestAccessor>
2813 inline void
2814 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
2815 pair<DestIterator, DestAccessor> const & dest,
2816 double innerScale, double outerScale,
2817 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
2818 {
2819 structureTensorMultiArray( source.first, source.second, source.third,
2820 dest.first, dest.second,
2821 innerScale, outerScale, opt);
2822 }
2823
2824 template <unsigned int N, class T1, class S1,
2825 class T2, class S2>
2826 inline void
2827 structureTensorMultiArray(MultiArrayView<N, T1, S1> const & source,
2828 MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2829 ConvolutionOptions<N> opt )
2830 {
2831 if(opt.to_point != typename MultiArrayShape<N>::type())
2832 {
2833 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
2834 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
2835 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2836 "structureTensorMultiArray(): shape mismatch between ROI and output.");
2837 }
2838 else
2839 {
2840 vigra_precondition(source.shape() == dest.shape(),
2841 "structureTensorMultiArray(): shape mismatch between input and output.");
2842 }
2843
2844 structureTensorMultiArray( srcMultiArrayRange(source),
2845 destMultiArray(dest), opt );
2846 }
2847
2848
2849 template <unsigned int N, class T1, class S1,
2850 class T2, class S2>
2851 inline void
2852 structureTensorMultiArray(MultiArrayView<N, T1, S1> const & source,
2853 MultiArrayView<N, TinyVector<T2, int(N*(N+1)/2)>, S2> dest,
2854 double innerScale, double outerScale,
2855 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
2856 {
2857 structureTensorMultiArray(source, dest, opt.innerScale(innerScale).outerScale(outerScale));
2858 }
2859
2860 //@}
2861
2862 } //-- namespace vigra
2863
2864
2865 #endif //-- VIGRA_MULTI_CONVOLUTION_H
vigra::VectorElementAccessor
Accessor for one component of a vector.
Definition: accessor.hxx:539
vigra::Kernel1D
Generic 1 dimensional convolution kernel.
Definition: separableconvolution.hxx:52
vigra::indexSort
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
vigra::Kernel1D::right
int right() const
Definition: separableconvolution.hxx:2157
vigra::Kernel1D::initGaussian
void initGaussian(double std_dev, value_type norm, double windowRatio=0.0)
Definition: separableconvolution.hxx:2253
vigra::MultiArrayView::bindOuter
MultiArrayView< N-M, T, StrideTag > bindOuter(const TinyVector< Index, M > &d) const
Definition: multi_array.hxx:2132
vigra::gaussianSmoothMultiArray
void gaussianSmoothMultiArray(...)
Isotropic Gaussian smoothing of a multi-dimensional arrays.
vigra::gaussianGradientMagnitude
void gaussianGradientMagnitude(...)
Calculate the gradient magnitude by means of a 1st derivatives of Gaussian filter.
vigra::convolveMultiArrayOneDimension
void convolveMultiArrayOneDimension(...)
Convolution along a single dimension of a multi-dimensional arrays.
vigra::squaredNorm
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
vigra::laplacianOfGaussianMultiArray
void laplacianOfGaussianMultiArray(...)
Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters.
vigra::MultiArray
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2422
vigra::TinyVectorBase::end
iterator end()
Definition: tinyvector.hxx:864
vigra::ConvolutionOptions::subarray
ConvolutionOptions< dim > & subarray(Shape const &from, Shape const &to)
Definition: multi_convolution.hxx:498
vigra::convolveLine
void convolveLine(...)
Performs a 1-dimensional convolution of the source signal using the given kernel. ...
vigra
Definition: accessor.hxx:43
vigra::gaussianDivergenceMultiArray
void gaussianDivergenceMultiArray(...)
Calculate the divergence of a vector field using Gaussian derivative filters.
vigra::MultiArrayView::insertSingletonDimension
MultiArrayView< N+1, T, StrideTag > insertSingletonDimension(difference_type_1 i) const
Definition: multi_array.hxx:2302
vigra::doxygen_overloaded_function
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
vigra::gaussianGradientMultiArray
void gaussianGradientMultiArray(...)
Calculate Gaussian gradient of a multi-dimensional arrays.
vigra::ArrayVector
Definition: array_vector.hxx:58
vigra::MultiArrayShape
Definition: multi_fwd.hxx:63
vigra::sq
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:365
vigra::MultiArrayNavigator
A navigator that provides access to the 1D subranges of an n-dimensional range given by a vigra::Mult...
Definition: navigator.hxx:97
vigra::MultiArrayView::shape
const difference_type & shape() const
Definition: multi_array.hxx:1596
vigra::TinyVectorBase::begin
iterator begin()
Definition: tinyvector.hxx:861
vigra::symmetricGradientMultiArray
void symmetricGradientMultiArray(...)
Calculate gradient of a multi-dimensional arrays using symmetric difference filters.
vigra::boost_graph::source
vigra::GridGraph< N, DirectedTag >::vertex_descriptor source(typename vigra::GridGraph< N, DirectedTag >::edge_descriptor const &e, vigra::GridGraph< N, DirectedTag > const &g)
Get a vertex descriptor for the start vertex of edge e in graph g (API: boost).
Definition: multi_gridgraph.hxx:2943
vigra::ConvolutionOptions::filterWindowSize
ConvolutionOptions< dim > & filterWindowSize(double ratio)
Definition: multi_convolution.hxx:476
vigra::ConvolutionOptions
Options class template for convolutions.
Definition: multi_convolution.hxx:335
vigra::copyMultiArray
void copyMultiArray(...)
Copy a multi-dimensional array.
vigra::TinyVector< double, N >
vigra::linalg::outer
void outer(const MultiArrayView< 2, T, C1 > &x, const MultiArrayView< 2, T, C2 > &y, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1457
vigra::hessianOfGaussianMultiArray
void hessianOfGaussianMultiArray(...)
Calculate Hessian matrix of a N-dimensional arrays using Gaussian derivative filters.
vigra::combineTwoMultiArrays
void combineTwoMultiArrays(...)
Combine two multi-dimensional arrays into one using a binary function or functor. ...
vigra::ArrayVectorView::begin
const_iterator begin() const
Definition: array_vector.hxx:223
vigra::StandardConstAccessor
Encapsulate read access to the values an iterator points to.
Definition: accessor.hxx:269
vigra::MultiArrayView
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:652
vigra::transformMultiArray
void transformMultiArray(...)
Transform a multi-dimensional array with a unary function or functor.
vigra::ConvolutionOptions::innerScale
ConvolutionOptions< dim > & innerScale(...)
vigra::MultiArrayView::init
MultiArrayView & init(const U &init)
Definition: multi_array.hxx:1154
vigra::RGBValue
Class for a single RGB value.
Definition: accessor.hxx:938
vigra::StandardAccessor
Encapsulate access to the values an iterator points to.
Definition: accessor.hxx:133
vigra::ArrayVectorView::end
const_iterator end() const
Definition: array_vector.hxx:237
vigra::structureTensorMultiArray
void structureTensorMultiArray(...)
Calculate th structure tensor of a multi-dimensional arrays.
vigra::MultiArrayView::expandElements
MultiArrayView< N+1, typename ExpandElementResult< T >::type, StridedArrayTag > expandElements(difference_type_1 d) const
Definition: multi_array.hxx:2273
vigra::sqrt
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
vigra::ConvolutionOptions::stdDev
ConvolutionOptions< dim > & stdDev(...)
vigra::separableConvolveMultiArray
void separableConvolveMultiArray(...)
Separated convolution on multi-dimensional arrays.

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.11.0