38#ifndef VIGRA_MULTI_CONVOLUTION_H
39#define VIGRA_MULTI_CONVOLUTION_H
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"
66 DoubleYielder(
double v,
unsigned,
const char *
const) : value(v) {}
67 DoubleYielder(
double v) : value(v) {}
69 double operator*()
const {
return value; }
73struct IteratorDoubleYielder
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; }
83struct SequenceDoubleYielder
85 typename X::const_iterator it;
86 SequenceDoubleYielder(
const X & seq,
unsigned dim,
87 const char *
const function_name =
"SequenceDoubleYielder")
90 if (seq.size() == dim)
92 std::string msg =
"(): Parameter number be equal to the number of spatial dimensions.";
93 vigra_precondition(
false, function_name + msg);
95 void operator++() { ++it; }
96 double operator*()
const {
return *it; }
100struct WrapDoubleIterator
103 typename IfBool< IsConvertibleTo<X, double>::value,
105 typename IfBool< IsIterator<X>::value || IsArray<X>::value,
106 IteratorDoubleYielder<X>,
107 SequenceDoubleYielder<X>
112template <
class Param1,
class Param2,
class Param3>
113struct WrapDoubleIteratorTriple
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) {}
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)
133 std::string msg =
"(): Scale must be positive.";
134 vigra_precondition(
false, function_name + msg);
137 double sigma_scaled(
const char *
const function_name =
"unknown function ",
138 bool allow_zero =
false)
const
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))
145 return std::sqrt(sigma_squared) / step_size();
149 std::string msg =
"(): Scale would be imaginary";
152 vigra_precondition(
false, function_name + msg +
".");
158template <
unsigned dim>
159struct multiArrayScaleParam
161 typedef TinyVector<double, dim> p_vector;
162 typedef typename p_vector::const_iterator return_type;
165 template <
class Param>
166 multiArrayScaleParam(Param val,
const char *
const function_name =
"multiArrayScaleParam")
168 typename WrapDoubleIterator<Param>::type in(val, dim, function_name);
169 for (
unsigned i = 0; i != dim; ++i, ++in)
172 return_type operator()()
const
176 static void precondition(
unsigned n_par,
const char *
const function_name =
"multiArrayScaleParam")
180 std::string msg =
"(): dimension parameter must be ";
181 vigra_precondition(dim == n_par, function_name + msg + n);
183 multiArrayScaleParam(
double v0,
double v1,
const char *
const function_name =
"multiArrayScaleParam")
185 precondition(2, function_name);
186 vec = p_vector(v0, v1);
188 multiArrayScaleParam(
double v0,
double v1,
double v2,
const char *
const function_name =
"multiArrayScaleParam")
190 precondition(3, function_name);
191 vec = p_vector(v0, v1, v2);
193 multiArrayScaleParam(
double v0,
double v1,
double v2,
double v3,
const char *
const function_name =
"multiArrayScaleParam")
195 precondition(4, function_name);
196 vec = p_vector(v0, v1, v2, v3);
198 multiArrayScaleParam(
double v0,
double v1,
double v2,
double v3,
double v4,
const char *
const function_name =
"multiArrayScaleParam")
200 precondition(5, function_name);
201 vec = p_vector(v0, v1, v2, v3, v4);
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) \
211 member_name = ParamVec(val, "ConvolutionOptions::" #function_name); \
214 ConvolutionOptions & function_name() \
216 member_name = ParamVec(default_value, "ConvolutionOptions::" #function_name); \
219 ConvolutionOptions & function_name(double v0, double v1) \
221 member_name = ParamVec(v0, v1, "ConvolutionOptions::" #function_name); \
224 ConvolutionOptions & function_name(double v0, double v1, double v2) \
226 member_name = ParamVec(v0, v1, v2, "ConvolutionOptions::" #function_name); \
229 ConvolutionOptions & function_name(double v0, double v1, double v2, double v3) \
231 member_name = ParamVec(v0, v1, v2, v3, "ConvolutionOptions::" #function_name); \
234 ConvolutionOptions & function_name(double v0, double v1, double v2, double v3, double v4) \
236 member_name = ParamVec(v0, v1, v2, v3, v4, "ConvolutionOptions::" #function_name); \
239 typename ParamVec::p_vector get##getter_setter_name()const{ \
240 return member_name.vec; \
242 void set##getter_setter_name(typename ParamVec::p_vector vec){ \
243 member_name.vec = vec; \
334template <
unsigned dim>
339 typedef detail::multiArrayScaleParam<dim> ParamVec;
340 typedef typename ParamVec::return_type ParamIt;
345 ParamVec outer_scale;
347 Shape from_point, to_point;
357 typedef typename detail::WrapDoubleIteratorTriple<ParamIt, ParamIt, ParamIt>
359 typedef typename detail::WrapDoubleIterator<ParamIt>::type
362 ScaleIterator scaleParams()
const
364 return ScaleIterator(sigma_eff(), sigma_d(), step_size());
366 StepIterator stepParams()
const
368 return StepIterator(step_size());
375 return outer.
stdDev(outer_scale()).resolutionStdDev(0.0);
380 VIGRA_CONVOLUTION_OPTIONS(
stepSize, 1.0, step_size, StepSize)
415 VIGRA_CONVOLUTION_OPTIONS(
stdDev, 0.0, sigma_eff, StdDev)
416 VIGRA_CONVOLUTION_OPTIONS(
innerScale, 0.0, sigma_eff, InnerScale)
446 VIGRA_CONVOLUTION_OPTIONS(
outerScale, 0.0, outer_scale, OuterScale)
478 vigra_precondition(ratio >= 0.0,
479 "ConvolutionOptions::filterWindowSize(): ratio must not be negative.");
480 window_ratio = ratio;
484 double getFilterWindowSize()
const {
505 std::pair<Shape, Shape> getSubarray()
const{
506 std::pair<Shape, Shape> res;
507 res.first = from_point;
508 res.second = to_point;
522template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
523 class DestIterator,
class DestAccessor,
class KernelIterator>
525internalSeparableConvolveMultiArrayTmp(
526 SrcIterator si, SrcShape
const & shape, SrcAccessor src,
527 DestIterator di, DestAccessor dest, KernelIterator kit)
529 enum { N = 1 + SrcIterator::level };
531 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
532 typedef typename AccessorTraits<TmpType>::default_accessor TmpAcessor;
535 ArrayVector<TmpType> tmp( shape[0] );
537 typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
538 typedef MultiArrayNavigator<DestIterator, N> DNavigator;
544 SNavigator snav( si, shape, 0 );
545 DNavigator dnav( di, shape, 0 );
547 for( ; snav.hasMore(); snav++, dnav++ )
550 copyLine(snav.begin(), snav.end(), src, tmp.begin(), acc);
553 destIter( dnav.begin(), dest ),
560 for(
int d = 1; d < N; ++d, ++kit )
562 DNavigator dnav( di, shape, d );
564 tmp.resize( shape[d] );
566 for( ; dnav.hasMore(); dnav++ )
569 copyLine(dnav.begin(), dnav.end(), dest, tmp.begin(), acc);
572 destIter( dnav.begin(), dest ),
584template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
585 class DestIterator,
class DestAccessor,
class KernelIterator>
587internalSeparableConvolveSubarray(
588 SrcIterator si, SrcShape
const & shape, SrcAccessor src,
589 DestIterator di, DestAccessor dest, KernelIterator kit,
590 SrcShape
const & start, SrcShape
const & stop)
592 enum { N = 1 + SrcIterator::level };
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;
599 SrcShape sstart, sstop, axisorder, tmpshape;
600 TinyVector<double, N> overhead;
601 for(
int k=0; k<N; ++k)
604 sstart[k] = start[k] - kit[k].right();
607 sstop[k] = stop[k] - kit[k].left();
608 if(sstop[k] > shape[k])
610 overhead[k] = double(sstop[k] - sstart[k]) / (stop[k] - start[k]);
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]];
618 MultiArray<N, TmpType> tmp(dstop);
620 typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
621 typedef MultiArrayNavigator<TmpIterator, N> TNavigator;
627 SNavigator snav( si, sstart, sstop, axisorder[0]);
628 TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[0]);
630 ArrayVector<TmpType> tmpline(sstop[axisorder[0]] - sstart[axisorder[0]]);
632 int lstart = start[axisorder[0]] - sstart[axisorder[0]];
633 int lstop = lstart + (stop[axisorder[0]] - start[axisorder[0]]);
635 for( ; snav.hasMore(); snav++, tnav++ )
638 copyLine(snav.begin(), snav.end(), src, tmpline.begin(), acc);
640 convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
641 destIter(tnav.begin(), acc),
642 kernel1d( kit[axisorder[0]] ), lstart, lstop);
647 for(
int d = 1; d < N; ++d)
649 TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[d]);
651 ArrayVector<TmpType> tmpline(dstop[axisorder[d]] - dstart[axisorder[d]]);
653 int lstart = start[axisorder[d]] - sstart[axisorder[d]];
654 int lstop = lstart + (stop[axisorder[d]] - start[axisorder[d]]);
656 for( ; tnav.hasMore(); tnav++ )
659 copyLine(tnav.begin(), tnav.end(), acc, tmpline.begin(), acc );
661 convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
662 destIter( tnav.begin() + lstart, acc ),
663 kernel1d( kit[axisorder[d]] ), lstart, lstop);
666 dstart[axisorder[d]] = lstart;
667 dstop[axisorder[d]] = lstop;
670 copyMultiArray(tmp.traverser_begin()+dstart, stop-start, acc, di, dest);
676scaleKernel(K & kernel,
double a)
678 for(
int i = kernel.left(); i <= kernel.right(); ++i)
679 kernel[i] = detail::RequiresExplicitCast<typename K::value_type>::cast(kernel[i] * a);
862template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
863 class DestIterator,
class DestAccessor,
class KernelIterator>
866 DestIterator d, DestAccessor dest,
867 KernelIterator kernels,
868 SrcShape start = SrcShape(),
869 SrcShape stop = SrcShape())
871 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
874 if(stop != SrcShape())
877 enum { N = 1 + SrcIterator::level };
878 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, start);
879 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, stop);
881 for(
int k=0; k<N; ++k)
882 vigra_precondition(0 <= start[k] && start[k] < stop[k] && stop[k] <= shape[k],
883 "separableConvolveMultiArray(): invalid subarray shape.");
885 detail::internalSeparableConvolveSubarray(s, shape, src, d, dest, kernels, start, stop);
887 else if(!IsSameType<TmpType, typename DestAccessor::value_type>::boolResult)
891 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src,
898 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, d, dest, kernels );
902template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
903 class DestIterator,
class DestAccessor,
class T>
906 DestIterator d, DestAccessor dest,
907 Kernel1D<T>
const & kernel,
908 SrcShape
const & start = SrcShape(),
909 SrcShape
const & stop = SrcShape())
911 ArrayVector<Kernel1D<T> > kernels(shape.size(), kernel);
916template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
917 class DestIterator,
class DestAccessor,
class KernelIterator>
920 pair<DestIterator, DestAccessor>
const & dest,
922 SrcShape
const & start = SrcShape(),
923 SrcShape
const & stop = SrcShape())
926 dest.first, dest.second, kit, start, stop );
929template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
930 class DestIterator,
class DestAccessor,
class T>
933 pair<DestIterator, DestAccessor>
const & dest,
934 Kernel1D<T>
const & kernel,
935 SrcShape
const & start = SrcShape(),
936 SrcShape
const & stop = SrcShape())
938 ArrayVector<Kernel1D<T> > kernels(source.second.size(), kernel);
941 dest.first, dest.second, kernels.begin(), start, stop);
944template <
unsigned int N,
class T1,
class S1,
946 class KernelIterator,
class SHAPE>
949 MultiArrayView<N, T2, S2> dest,
951 SHAPE start, SHAPE stop)
955 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), start);
956 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), stop);
957 vigra_precondition(dest.shape() == (stop - start),
958 "separableConvolveMultiArray(): shape mismatch between ROI and output.");
962 vigra_precondition(source.shape() == dest.shape(),
963 "separableConvolveMultiArray(): shape mismatch between input and output.");
966 destMultiArray(dest), kit, start, stop );
969template <
unsigned int N,
class T1,
class S1,
971 class KernelIterator>
974 MultiArrayView<N, T2, S2> dest,
981template <
unsigned int N,
class T1,
class S1,
983 class T,
class SHAPE>
986 MultiArrayView<N, T2, S2> dest,
987 Kernel1D<T>
const & kernel,
988 SHAPE
const & start, SHAPE
const & stop)
990 ArrayVector<Kernel1D<T> > kernels(N, kernel);
994template <
unsigned int N,
class T1,
class S1,
999 MultiArrayView<N, T2, S2> dest,
1000 Kernel1D<T>
const & kernel)
1095template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1096 class DestIterator,
class DestAccessor,
class T>
1099 DestIterator d, DestAccessor dest,
1101 SrcShape
const & start = SrcShape(),
1102 SrcShape
const & stop = SrcShape())
1104 enum { N = 1 + SrcIterator::level };
1105 vigra_precondition( dim < N,
1106 "convolveMultiArrayOneDimension(): The dimension number to convolve must be smaller "
1107 "than the data dimensionality" );
1109 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
1116 SrcShape sstart, sstop(shape), dstart, dstop(shape);
1118 if(stop != SrcShape())
1123 sstop[dim] = shape[dim];
1124 dstop = stop - start;
1127 SNavigator snav( s, sstart, sstop, dim );
1128 DNavigator dnav( d, dstart, dstop, dim );
1130 for( ; snav.hasMore(); snav++, dnav++ )
1133 copyLine(snav.begin(), snav.end(), src,
1137 destIter( dnav.begin(), dest ),
1138 kernel1d( kernel), start[dim], stop[dim]);
1142template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1143 class DestIterator,
class DestAccessor,
class T>
1146 pair<DestIterator, DestAccessor>
const & dest,
1148 Kernel1D<T>
const & kernel,
1149 SrcShape
const & start = SrcShape(),
1150 SrcShape
const & stop = SrcShape())
1153 dest.first, dest.second, dim, kernel, start, stop);
1156template <
unsigned int N,
class T1,
class S1,
1158 class T,
class SHAPE>
1161 MultiArrayView<N, T2, S2> dest,
1163 Kernel1D<T>
const & kernel,
1164 SHAPE start, SHAPE stop)
1168 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), start);
1169 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), stop);
1170 vigra_precondition(dest.shape() == (stop - start),
1171 "convolveMultiArrayOneDimension(): shape mismatch between ROI and output.");
1175 vigra_precondition(source.shape() == dest.shape(),
1176 "convolveMultiArrayOneDimension(): shape mismatch between input and output.");
1179 destMultiArray(dest), dim, kernel, start, stop);
1182template <
unsigned int N,
class T1,
class S1,
1187 MultiArrayView<N, T2, S2> dest,
1189 Kernel1D<T>
const & kernel)
1333template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1334 class DestIterator,
class DestAccessor>
1337 DestIterator d, DestAccessor dest,
1339 const char *
const function_name =
"gaussianSmoothMultiArray" )
1341 static const int N = SrcShape::static_size;
1343 typename ConvolutionOptions<N>::ScaleIterator params = opt.scaleParams();
1346 for (
int dim = 0; dim < N; ++dim, ++params)
1347 kernels[dim].initGaussian(params.sigma_scaled(function_name,
true),
1348 1.0, opt.window_ratio);
1353template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1354 class DestIterator,
class DestAccessor>
1357 DestIterator d, DestAccessor dest,
double sigma,
1358 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1360 ConvolutionOptions<SrcShape::static_size> par = opt;
1364template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1365 class DestIterator,
class DestAccessor>
1368 pair<DestIterator, DestAccessor>
const & dest,
1369 const ConvolutionOptions<SrcShape::static_size> & opt)
1372 dest.first, dest.second, opt );
1375template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1376 class DestIterator,
class DestAccessor>
1379 pair<DestIterator, DestAccessor>
const & dest,
double sigma,
1380 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1383 dest.first, dest.second, sigma, opt );
1386template <
unsigned int N,
class T1,
class S1,
1390 MultiArrayView<N, T2, S2> dest,
1391 ConvolutionOptions<N> opt)
1395 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
1396 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
1397 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1398 "gaussianSmoothMultiArray(): shape mismatch between ROI and output.");
1402 vigra_precondition(source.shape() == dest.shape(),
1403 "gaussianSmoothMultiArray(): shape mismatch between input and output.");
1407 destMultiArray(dest), opt );
1410template <
unsigned int N,
class T1,
class S1,
1414 MultiArrayView<N, T2, S2> dest,
1416 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
1545template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1546 class DestIterator,
class DestAccessor>
1549 DestIterator di, DestAccessor dest,
1551 const char *
const function_name =
"gaussianGradientMultiArray")
1553 typedef typename DestAccessor::value_type DestType;
1554 typedef typename DestType::value_type DestValueType;
1555 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1557 static const int N = SrcShape::static_size;
1558 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
1560 for(
int k=0; k<N; ++k)
1564 vigra_precondition(N == (
int)dest.size(di),
1565 "gaussianGradientMultiArray(): Wrong number of channels in output array.");
1567 ParamType params = opt.scaleParams();
1568 ParamType params2(params);
1571 for (
int dim = 0; dim < N; ++dim, ++params)
1573 double sigma = params.sigma_scaled(function_name);
1574 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
1580 for (
int dim = 0; dim < N; ++dim, ++params2)
1583 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 1, 1.0, opt.window_ratio);
1584 detail::scaleKernel(kernels[dim], 1.0 / params2.step_size());
1586 opt.from_point, opt.to_point);
1590template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1591 class DestIterator,
class DestAccessor>
1594 DestIterator di, DestAccessor dest,
double sigma,
1595 ConvolutionOptions<SrcShape::static_size> opt = ConvolutionOptions<SrcShape::static_size>())
1600template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1601 class DestIterator,
class DestAccessor>
1604 pair<DestIterator, DestAccessor>
const & dest,
1605 ConvolutionOptions<SrcShape::static_size>
const & opt )
1608 dest.first, dest.second, opt );
1611template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1612 class DestIterator,
class DestAccessor>
1615 pair<DestIterator, DestAccessor>
const & dest,
1617 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1620 dest.first, dest.second, sigma, opt );
1623template <
unsigned int N,
class T1,
class S1,
1627 MultiArrayView<N, TinyVector<T2,
int(N)>, S2> dest,
1628 ConvolutionOptions<N> opt )
1632 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
1633 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
1634 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1635 "gaussianGradientMultiArray(): shape mismatch between ROI and output.");
1639 vigra_precondition(source.shape() == dest.shape(),
1640 "gaussianGradientMultiArray(): shape mismatch between input and output.");
1644 destMultiArray(dest), opt );
1647template <
unsigned int N,
class T1,
class S1,
1651 MultiArrayView<N, TinyVector<T2,
int(N)>, S2> dest,
1653 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
1666template <
unsigned int N,
class T1,
class S1,
1669gaussianGradientMagnitudeImpl(MultiArrayView<N+1, T1, S1>
const & src,
1670 MultiArrayView<N, T2, S2> dest,
1671 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
1676 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.from_point);
1677 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.to_point);
1678 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1679 "gaussianGradientMagnitude(): shape mismatch between ROI and output.");
1683 vigra_precondition(shape == dest.shape(),
1684 "gaussianGradientMagnitude(): shape mismatch between input and output.");
1689 typedef typename NumericTraits<T1>::RealPromote TmpType;
1690 MultiArray<N, TinyVector<TmpType, int(N)> > grad(dest.shape());
1692 using namespace multi_math;
1694 for(
int k=0; k<src.shape(N); ++k)
1706template <
unsigned int N,
class T1,
class S1,
1710 MultiArrayView<N, T2, S2> dest,
1711 ConvolutionOptions<N>
const & opt)
1713 detail::gaussianGradientMagnitudeImpl<N, T1>(src, dest, opt);
1716template <
unsigned int N,
class T1,
class S1,
1720 MultiArrayView<N, T2, S2> dest,
1721 ConvolutionOptions<N>
const & opt)
1723 detail::gaussianGradientMagnitudeImpl<N, T1>(src.insertSingletonDimension(N), dest, opt);
1726template <
unsigned int N,
class T1,
int M,
class S1,
1730 MultiArrayView<N, T2, S2> dest,
1731 ConvolutionOptions<N>
const & opt)
1733 detail::gaussianGradientMagnitudeImpl<N, T1>(src.expandElements(N), dest, opt);
1736template <
unsigned int N,
class T1,
unsigned int R,
unsigned int G,
unsigned int B,
class S1,
1740 MultiArrayView<N, T2, S2> dest,
1741 ConvolutionOptions<N>
const & opt)
1743 detail::gaussianGradientMagnitudeImpl<N, T1>(src.expandElements(N), dest, opt);
1746template <
unsigned int N,
class T1,
class S1,
1750 MultiArrayView<N, T2, S2> dest,
1752 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
1757template <
unsigned int N,
class T1,
class S1,
1761 MultiArrayView<N, T2, S2> dest,
1763 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
1878template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1879 class DestIterator,
class DestAccessor>
1882 DestIterator di, DestAccessor dest,
1885 typedef typename DestAccessor::value_type DestType;
1886 typedef typename DestType::value_type DestValueType;
1887 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1889 static const int N = SrcShape::static_size;
1890 typedef typename ConvolutionOptions<N>::StepIterator StepType;
1892 for(
int k=0; k<N; ++k)
1896 vigra_precondition(N == (
int)dest.size(di),
1897 "symmetricGradientMultiArray(): Wrong number of channels in output array.");
1900 filter.initSymmetricDifference();
1902 StepType step_size_it = opt.stepParams();
1907 for (
int d = 0; d < N; ++d, ++step_size_it)
1910 detail::scaleKernel(symmetric, 1 / *step_size_it);
1912 di, ElementAccessor(d, dest),
1913 d, symmetric, opt.from_point, opt.to_point);
1917template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1918 class DestIterator,
class DestAccessor>
1921 pair<DestIterator, DestAccessor>
const & dest,
1922 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1925 dest.first, dest.second, opt);
1928template <
unsigned int N,
class T1,
class S1,
1932 MultiArrayView<N, TinyVector<T2,
int(N)>, S2> dest,
1933 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
1937 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
1938 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
1939 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1940 "symmetricGradientMultiArray(): shape mismatch between ROI and output.");
1944 vigra_precondition(source.shape() == dest.shape(),
1945 "symmetricGradientMultiArray(): shape mismatch between input and output.");
1949 destMultiArray(dest), opt);
2072template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2073 class DestIterator,
class DestAccessor>
2076 DestIterator di, DestAccessor dest,
2079 using namespace functor;
2081 typedef typename DestAccessor::value_type DestType;
2082 typedef typename NumericTraits<DestType>::RealPromote KernelType;
2085 static const int N = SrcShape::static_size;
2086 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
2088 ParamType params = opt.scaleParams();
2089 ParamType params2(params);
2092 for (
int dim = 0; dim < N; ++dim, ++params)
2094 double sigma = params.sigma_scaled(
"laplacianOfGaussianMultiArray");
2095 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
2098 SrcShape dshape(shape);
2099 if(opt.to_point != SrcShape())
2100 dshape = opt.to_point - opt.from_point;
2105 for (
int dim = 0; dim < N; ++dim, ++params2)
2108 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 2, 1.0, opt.window_ratio);
2109 detail::scaleKernel(kernels[dim], 1.0 /
sq(params2.step_size()));
2114 di, dest, kernels.
begin(), opt.from_point, opt.to_point);
2119 derivative.traverser_begin(), DerivativeAccessor(),
2120 kernels.
begin(), opt.from_point, opt.to_point);
2122 di, dest, Arg1() + Arg2() );
2127template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2128 class DestIterator,
class DestAccessor>
2131 DestIterator di, DestAccessor dest,
double sigma,
2132 ConvolutionOptions<SrcShape::static_size> opt = ConvolutionOptions<SrcShape::static_size>())
2137template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2138 class DestIterator,
class DestAccessor>
2141 pair<DestIterator, DestAccessor>
const & dest,
2142 ConvolutionOptions<SrcShape::static_size>
const & opt )
2145 dest.first, dest.second, opt );
2148template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2149 class DestIterator,
class DestAccessor>
2152 pair<DestIterator, DestAccessor>
const & dest,
2154 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
2157 dest.first, dest.second, sigma, opt );
2160template <
unsigned int N,
class T1,
class S1,
2164 MultiArrayView<N, T2, S2> dest,
2165 ConvolutionOptions<N> opt )
2169 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
2170 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
2171 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2172 "laplacianOfGaussianMultiArray(): shape mismatch between ROI and output.");
2176 vigra_precondition(source.shape() == dest.shape(),
2177 "laplacianOfGaussianMultiArray(): shape mismatch between input and output.");
2181 destMultiArray(dest), opt );
2184template <
unsigned int N,
class T1,
class S1,
2188 MultiArrayView<N, T2, S2> dest,
2190 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
2304template <
class Iterator,
2305 unsigned int N,
class T,
class S>
2311 typedef typename std::iterator_traits<Iterator>::value_type ArrayType;
2312 typedef typename ArrayType::value_type SrcType;
2313 typedef typename NumericTraits<SrcType>::RealPromote TmpType;
2316 vigra_precondition(std::distance(vectorField, vectorFieldEnd) == N,
2317 "gaussianDivergenceMultiArray(): wrong number of input arrays.");
2320 typename ConvolutionOptions<N>::ScaleIterator params = opt.scaleParams();
2323 for(
unsigned int k = 0; k < N; ++k, ++params)
2325 sigmas[k] = params.sigma_scaled(
"gaussianDivergenceMultiArray");
2326 kernels[k].initGaussian(sigmas[k], 1.0, opt.window_ratio);
2331 for(
unsigned int k=0; k < N; ++k, ++vectorField)
2333 kernels[k].initGaussianDerivative(sigmas[k], 1, 1.0, opt.window_ratio);
2341 divergence += tmpDeriv;
2343 kernels[k].initGaussian(sigmas[k], 1.0, opt.window_ratio);
2347template <
class Iterator,
2348 unsigned int N,
class T,
class S>
2351 MultiArrayView<N, T, S> divergence,
2353 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
2358template <
unsigned int N,
class T1,
class S1,
2362 MultiArrayView<N, T2, S2> divergence,
2363 ConvolutionOptions<N>
const & opt)
2365 ArrayVector<MultiArrayView<N, T1> > field;
2366 for(
unsigned int k=0; k<N; ++k)
2367 field.push_back(vectorField.bindElementChannel(k));
2372template <
unsigned int N,
class T1,
class S1,
2376 MultiArrayView<N, T2, S2> divergence,
2378 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
2505template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2506 class DestIterator,
class DestAccessor>
2509 DestIterator di, DestAccessor dest,
2512 typedef typename DestAccessor::value_type DestType;
2513 typedef typename DestType::value_type DestValueType;
2514 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
2516 static const int N = SrcShape::static_size;
2517 static const int M = N*(N+1)/2;
2518 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
2520 for(
int k=0; k<N; ++k)
2524 vigra_precondition(M == (
int)dest.size(di),
2525 "hessianOfGaussianMultiArray(): Wrong number of channels in output array.");
2527 ParamType params_init = opt.scaleParams();
2530 ParamType params(params_init);
2531 for (
int dim = 0; dim < N; ++dim, ++params)
2533 double sigma = params.sigma_scaled(
"hessianOfGaussianMultiArray");
2534 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
2540 ParamType params_i(params_init);
2541 for (
int b=0, i=0; i<N; ++i, ++params_i)
2543 ParamType params_j(params_i);
2544 for (
int j=i; j<N; ++j, ++b, ++params_j)
2549 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 2, 1.0, opt.window_ratio);
2553 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 1, 1.0, opt.window_ratio);
2554 kernels[j].initGaussianDerivative(params_j.sigma_scaled(), 1, 1.0, opt.window_ratio);
2556 detail::scaleKernel(kernels[i], 1 / params_i.step_size());
2557 detail::scaleKernel(kernels[j], 1 / params_j.step_size());
2559 kernels.
begin(), opt.from_point, opt.to_point);
2564template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2565 class DestIterator,
class DestAccessor>
2568 DestIterator di, DestAccessor dest,
double sigma,
2569 ConvolutionOptions<SrcShape::static_size> opt = ConvolutionOptions<SrcShape::static_size>())
2574template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2575 class DestIterator,
class DestAccessor>
2578 pair<DestIterator, DestAccessor>
const & dest,
2579 ConvolutionOptions<SrcShape::static_size>
const & opt )
2582 dest.first, dest.second, opt );
2585template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2586 class DestIterator,
class DestAccessor>
2589 pair<DestIterator, DestAccessor>
const & dest,
2591 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
2594 dest.first, dest.second, sigma, opt );
2597template <
unsigned int N,
class T1,
class S1,
2601 MultiArrayView<N, TinyVector<T2,
int(N*(N+1)/2)>, S2> dest,
2602 ConvolutionOptions<N> opt )
2606 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
2607 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
2608 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2609 "hessianOfGaussianMultiArray(): shape mismatch between ROI and output.");
2613 vigra_precondition(source.shape() == dest.shape(),
2614 "hessianOfGaussianMultiArray(): shape mismatch between input and output.");
2618 destMultiArray(dest), opt );
2621template <
unsigned int N,
class T1,
class S1,
2625 MultiArrayView<N, TinyVector<T2,
int(N*(N+1)/2)>, S2> dest,
2627 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
2634template<
int N,
class VectorType>
2635struct StructurTensorFunctor
2637 typedef VectorType result_type;
2638 typedef typename VectorType::value_type ValueType;
2641 VectorType operator()(T
const & in)
const
2644 for(
int b=0, i=0; i<N; ++i)
2646 for(
int j=i; j<N; ++j, ++b)
2648 res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
2781template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2782 class DestIterator,
class DestAccessor>
2785 DestIterator di, DestAccessor dest,
2788 static const int N = SrcShape::static_size;
2789 static const int M = N*(N+1)/2;
2791 typedef typename DestAccessor::value_type DestType;
2792 typedef typename DestType::value_type DestValueType;
2793 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
2798 for(
int k=0; k<N; ++k)
2802 vigra_precondition(M == (
int)dest.size(di),
2803 "structureTensorMultiArray(): Wrong number of channels in output array.");
2807 typename ConvolutionOptions<N>::ScaleIterator params = outerOptions.scaleParams();
2809 SrcShape gradientShape(shape);
2810 if(opt.to_point != SrcShape())
2812 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.from_point);
2813 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.to_point);
2815 for(
int k=0; k<N; ++k, ++params)
2818 gauss.
initGaussian(params.sigma_scaled(
"structureTensorMultiArray"), 1.0, opt.window_ratio);
2819 int dilation = gauss.
right();
2820 innerOptions.from_point[k] = std::max<MultiArrayIndex>(0, opt.from_point[k] - dilation);
2821 innerOptions.to_point[k] = std::min<MultiArrayIndex>(shape[k], opt.to_point[k] + dilation);
2823 outerOptions.from_point -= innerOptions.from_point;
2824 outerOptions.to_point -= innerOptions.from_point;
2825 gradientShape = innerOptions.to_point - innerOptions.from_point;
2831 gradient.traverser_begin(), GradientAccessor(),
2833 "structureTensorMultiArray");
2836 gradientTensor.traverser_begin(), GradientTensorAccessor(),
2837 detail::StructurTensorFunctor<N, DestType>());
2840 di, dest, outerOptions,
2841 "structureTensorMultiArray");
2844template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2845 class DestIterator,
class DestAccessor>
2848 DestIterator di, DestAccessor dest,
2849 double innerScale,
double outerScale,
2850 ConvolutionOptions<SrcShape::static_size> opt = ConvolutionOptions<SrcShape::static_size>())
2853 opt.stdDev(innerScale).outerScale(outerScale));
2856template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2857 class DestIterator,
class DestAccessor>
2860 pair<DestIterator, DestAccessor>
const & dest,
2861 ConvolutionOptions<SrcShape::static_size>
const & opt )
2864 dest.first, dest.second, opt );
2868template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2869 class DestIterator,
class DestAccessor>
2872 pair<DestIterator, DestAccessor>
const & dest,
2873 double innerScale,
double outerScale,
2874 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
2877 dest.first, dest.second,
2878 innerScale, outerScale, opt);
2881template <
unsigned int N,
class T1,
class S1,
2885 MultiArrayView<N, TinyVector<T2,
int(N*(N+1)/2)>, S2> dest,
2886 ConvolutionOptions<N> opt )
2890 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.from_point);
2891 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.shape(), opt.to_point);
2892 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2893 "structureTensorMultiArray(): shape mismatch between ROI and output.");
2897 vigra_precondition(source.shape() == dest.shape(),
2898 "structureTensorMultiArray(): shape mismatch between input and output.");
2902 destMultiArray(dest), opt );
2906template <
unsigned int N,
class T1,
class S1,
2910 MultiArrayView<N, TinyVector<T2,
int(N*(N+1)/2)>, S2> dest,
2911 double innerScale,
double outerScale,
2912 ConvolutionOptions<N> opt = ConvolutionOptions<N>())
const_iterator begin() const
Definition array_vector.hxx:223
const_iterator end() const
Definition array_vector.hxx:237
Definition array_vector.hxx:514
Options class template for convolutions.
Definition multi_convolution.hxx:336
ConvolutionOptions< dim > & stdDev(...)
ConvolutionOptions< dim > & resolutionStdDev(...)
ConvolutionOptions< dim > & filterWindowSize(double ratio)
Definition multi_convolution.hxx:476
ConvolutionOptions< dim > & innerScale(...)
ConvolutionOptions< dim > & outerScale(...)
ConvolutionOptions< dim > & subarray(Shape const &from, Shape const &to)
Definition multi_convolution.hxx:498
ConvolutionOptions< dim > & stepSize(...)
Generic 1 dimensional convolution kernel.
Definition separableconvolution.hxx:1367
int right() const
Definition separableconvolution.hxx:2157
void initGaussian(double std_dev, value_type norm, double windowRatio=0.0)
Definition separableconvolution.hxx:2253
A navigator that provides access to the 1D subranges of an n-dimensional range given by a vigra::Mult...
Definition navigator.hxx:101
TinyVector< MultiArrayIndex, N > type
Definition multi_shape.hxx:272
Base class for, and view to, vigra::MultiArray.
Definition multi_fwd.hxx:127
const difference_type & shape() const
Definition multi_array.hxx:1650
Main MultiArray class containing the memory management.
Definition multi_fwd.hxx:131
Encapsulate access to the values an iterator points to.
Definition accessor.hxx:134
Encapsulate read access to the values an iterator points to.
Definition accessor.hxx:270
Class for fixed size vectors.
Definition tinyvector.hxx:1008
Accessor for one component of a vector.
Definition accessor.hxx:540
void gaussianGradientMagnitude(...)
Calculate the gradient magnitude by means of a 1st derivatives of Gaussian filter.
void convolveLine(...)
Performs a 1-dimensional convolution of the source signal using the given kernel.
void separableConvolveMultiArray(...)
Separated convolution on multi-dimensional arrays.
void structureTensorMultiArray(...)
Calculate the structure tensor of a multi-dimensional arrays.
void indexSort(Iterator first, Iterator last, IndexIterator index_first, Compare c)
Return the index permutation that would sort the input array.
Definition algorithm.hxx:414
void laplacianOfGaussianMultiArray(...)
Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters.
void gaussianDivergenceMultiArray(...)
Calculate the divergence of a vector field using Gaussian derivative filters.
void gaussianSmoothMultiArray(...)
Isotropic Gaussian smoothing of a multi-dimensional arrays.
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition fixedpoint.hxx:616
NumericTraits< T >::Promote sq(T t)
The square function.
Definition mathutil.hxx:382
void convolveMultiArrayOneDimension(...)
Convolution along a single dimension of a multi-dimensional arrays.
void transformMultiArray(...)
Transform a multi-dimensional array with a unary function or functor.
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition fftw3.hxx:1044
void combineTwoMultiArrays(...)
Combine two multi-dimensional arrays into one using a binary function or functor.
void hessianOfGaussianMultiArray(...)
Calculate Hessian matrix of a N-dimensional arrays using Gaussian derivative filters.
void copyMultiArray(...)
Copy a multi-dimensional array.
void gaussianGradientMultiArray(...)
Calculate Gaussian gradient of a multi-dimensional arrays.
void symmetricGradientMultiArray(...)
Calculate gradient of a multi-dimensional arrays using symmetric difference filters.