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 "functorexpression.hxx"
50 #include "tinyvector.hxx"
51 #include "algorithm.hxx"
62 DoubleYielder(
double v,
unsigned,
const char *
const) : value(v) {}
63 DoubleYielder(
double v) : value(v) {}
65 double operator*()
const {
return value; }
69 struct IteratorDoubleYielder
72 IteratorDoubleYielder(X i,
unsigned,
const char *
const) : it(i) {}
73 IteratorDoubleYielder(X i) : it(i) {}
74 void operator++() { ++it; }
75 double operator*()
const {
return *it; }
79 struct SequenceDoubleYielder
81 typename X::const_iterator it;
82 SequenceDoubleYielder(
const X & seq,
unsigned dim,
83 const char *
const function_name =
"SequenceDoubleYielder")
86 if (seq.size() == dim)
88 std::string msg =
"(): Parameter number be equal to the number of spatial dimensions.";
89 vigra_precondition(
false, function_name + msg);
91 void operator++() { ++it; }
92 double operator*()
const {
return *it; }
96 struct WrapDoubleIterator
99 typename IfBool< IsConvertibleTo<X, double>::value,
101 typename IfBool< IsIterator<X>::value || IsArray<X>::value,
102 IteratorDoubleYielder<X>,
103 SequenceDoubleYielder<X>
108 template <
class Param1,
class Param2,
class Param3>
109 struct WrapDoubleIteratorTriple
111 typename WrapDoubleIterator<Param1>::type sigma_eff_it;
112 typename WrapDoubleIterator<Param2>::type sigma_d_it;
113 typename WrapDoubleIterator<Param3>::type step_size_it;
114 WrapDoubleIteratorTriple(Param1 sigma_eff, Param2 sigma_d, Param3 step_size)
115 : sigma_eff_it(sigma_eff), sigma_d_it(sigma_d), step_size_it(step_size) {}
122 double sigma_eff()
const {
return *sigma_eff_it; }
123 double sigma_d()
const {
return *sigma_d_it; }
124 double step_size()
const {
return *step_size_it; }
125 static double sqr(
double x) {
return x * x; }
126 static void sigma_precondition(
double sigma,
const char *
const function_name)
130 std::string msg =
"(): Scale must be positive.";
131 vigra_precondition(
false, function_name + msg);
134 double sigma_scaled(
const char *
const function_name =
"unknown function ")
const
136 sigma_precondition(sigma_eff(), function_name);
137 sigma_precondition(sigma_d(), function_name);
138 double sigma_squared = sqr(sigma_eff()) - sqr(sigma_d());
139 if (sigma_squared > 0.0)
141 return std::sqrt(sigma_squared) / step_size();
145 std::string msg =
"(): Scale would be imaginary or zero.";
146 vigra_precondition(
false, function_name + msg);
152 template <
unsigned dim>
153 struct multiArrayScaleParam
155 typedef TinyVector<double, dim> p_vector;
156 typedef typename p_vector::const_iterator return_type;
159 template <
class Param>
160 multiArrayScaleParam(Param val,
const char *
const function_name =
"multiArrayScaleParam")
162 typename WrapDoubleIterator<Param>::type in(val, dim, function_name);
163 for (
unsigned i = 0; i != dim; ++i, ++in)
166 return_type operator()()
const
170 static void precondition(
unsigned n_par,
const char *
const function_name =
"multiArrayScaleParam")
174 std::string msg =
"(): dimension parameter must be ";
175 vigra_precondition(dim == n_par, function_name + msg + n);
177 multiArrayScaleParam(
double v0,
double v1,
const char *
const function_name =
"multiArrayScaleParam")
179 precondition(2, function_name);
180 vec = p_vector(v0, v1);
182 multiArrayScaleParam(
double v0,
double v1,
double v2,
const char *
const function_name =
"multiArrayScaleParam")
184 precondition(3, function_name);
185 vec = p_vector(v0, v1, v2);
187 multiArrayScaleParam(
double v0,
double v1,
double v2,
double v3,
const char *
const function_name =
"multiArrayScaleParam")
189 precondition(4, function_name);
190 vec = p_vector(v0, v1, v2, v3);
192 multiArrayScaleParam(
double v0,
double v1,
double v2,
double v3,
double v4,
const char *
const function_name =
"multiArrayScaleParam")
194 precondition(5, function_name);
195 vec = p_vector(v0, v1, v2, v3, v4);
201 #define VIGRA_CONVOLUTION_OPTIONS(function_name, default_value, member_name) \
202 template <class Param> \
203 ConvolutionOptions & function_name(const Param & val) \
205 member_name = ParamVec(val, "ConvolutionOptions::" #function_name); \
208 ConvolutionOptions & function_name() \
210 member_name = ParamVec(default_value, "ConvolutionOptions::" #function_name); \
213 ConvolutionOptions & function_name(double v0, double v1) \
215 member_name = ParamVec(v0, v1, "ConvolutionOptions::" #function_name); \
218 ConvolutionOptions & function_name(double v0, double v1, double v2) \
220 member_name = ParamVec(v0, v1, v2, "ConvolutionOptions::" #function_name); \
223 ConvolutionOptions & function_name(double v0, double v1, double v2, double v3) \
225 member_name = ParamVec(v0, v1, v2, v3, "ConvolutionOptions::" #function_name); \
228 ConvolutionOptions & function_name(double v0, double v1, double v2, double v3, double v4) \
230 member_name = ParamVec(v0, v1, v2, v3, v4, "ConvolutionOptions::" #function_name); \
322 template <
unsigned dim>
327 typedef detail::multiArrayScaleParam<dim> ParamVec;
328 typedef typename ParamVec::return_type ParamIt;
333 ParamVec outer_scale;
335 Shape from_point, to_point;
345 typedef typename detail::WrapDoubleIteratorTriple<ParamIt, ParamIt, ParamIt>
347 typedef typename detail::WrapDoubleIterator<ParamIt>::type
350 ScaleIterator scaleParams()
const
352 return ScaleIterator(sigma_eff(), sigma_d(), step_size());
354 StepIterator stepParams()
const
356 return StepIterator(step_size());
363 return outer.
stdDev(outer_scale()).resolutionStdDev(0.0);
368 VIGRA_CONVOLUTION_OPTIONS(
stepSize, 1.0, step_size)
403 VIGRA_CONVOLUTION_OPTIONS(
stdDev, 0.0, sigma_eff)
404 VIGRA_CONVOLUTION_OPTIONS(
innerScale, 0.0, sigma_eff)
434 VIGRA_CONVOLUTION_OPTIONS(
outerScale, 0.0, outer_scale)
466 vigra_precondition(ratio >= 0.0,
467 "ConvolutionOptions::filterWindowSize(): ratio must not be negative.");
468 window_ratio = ratio;
497 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
498 class DestIterator,
class DestAccessor,
class KernelIterator>
500 internalSeparableConvolveMultiArrayTmp(
501 SrcIterator si, SrcShape
const & shape, SrcAccessor src,
502 DestIterator di, DestAccessor dest, KernelIterator kit)
504 enum { N = 1 + SrcIterator::level };
506 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
507 typedef typename AccessorTraits<TmpType>::default_accessor TmpAcessor;
510 ArrayVector<TmpType> tmp( shape[0] );
512 typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
513 typedef MultiArrayNavigator<DestIterator, N> DNavigator;
519 SNavigator snav( si, shape, 0 );
520 DNavigator dnav( di, shape, 0 );
522 for( ; snav.hasMore(); snav++, dnav++ )
525 copyLine(snav.begin(), snav.end(), src, tmp.begin(), acc);
528 destIter( dnav.begin(), dest ),
535 for(
int d = 1; d < N; ++d, ++kit )
537 DNavigator dnav( di, shape, d );
539 tmp.resize( shape[d] );
541 for( ; dnav.hasMore(); dnav++ )
544 copyLine(dnav.begin(), dnav.end(), dest, tmp.begin(), acc);
547 destIter( dnav.begin(), dest ),
559 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
560 class DestIterator,
class DestAccessor,
class KernelIterator>
562 internalSeparableConvolveSubarray(
563 SrcIterator si, SrcShape
const & shape, SrcAccessor src,
564 DestIterator di, DestAccessor dest, KernelIterator kit,
565 SrcShape
const & start, SrcShape
const & stop)
567 enum { N = 1 + SrcIterator::level };
569 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
570 typedef MultiArray<N, TmpType> TmpArray;
571 typedef typename TmpArray::traverser TmpIterator;
572 typedef typename AccessorTraits<TmpType>::default_accessor TmpAcessor;
574 SrcShape sstart, sstop, axisorder, tmpshape;
575 TinyVector<double, N> overhead;
576 for(
int k=0; k<N; ++k)
578 sstart[k] = start[k] - kit[k].right();
581 sstop[k] = stop[k] - kit[k].left();
582 if(sstop[k] > shape[k])
584 overhead[k] = double(sstop[k] - sstart[k]) / (stop[k] - start[k]);
587 indexSort(overhead.begin(), overhead.end(), axisorder.begin(), std::greater<double>());
589 SrcShape dstart, dstop(sstop - sstart);
590 dstop[axisorder[0]] = stop[axisorder[0]] - start[axisorder[0]];
593 MultiArray<N, TmpType> tmp(dstop);
595 typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
596 typedef MultiArrayNavigator<TmpIterator, N> TNavigator;
597 typedef MultiArrayNavigator<DestIterator, N> DNavigator;
603 SNavigator snav( si, sstart, sstop, axisorder[0]);
604 TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[0]);
606 ArrayVector<TmpType> tmpline(sstop[axisorder[0]] - sstart[axisorder[0]]);
608 int lstart = start[axisorder[0]] - sstart[axisorder[0]];
609 int lstop = lstart + (stop[axisorder[0]] - start[axisorder[0]]);
611 for( ; snav.hasMore(); snav++, tnav++ )
614 copyLine(snav.begin(), snav.end(), src, tmpline.begin(), acc);
616 convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
617 destIter(tnav.begin(), acc),
618 kernel1d( kit[axisorder[0]] ), lstart, lstop);
623 for(
int d = 1; d < N; ++d)
625 TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[d]);
627 ArrayVector<TmpType> tmpline(dstop[axisorder[d]] - dstart[axisorder[d]]);
629 int lstart = start[axisorder[d]] - sstart[axisorder[d]];
630 int lstop = lstart + (stop[axisorder[d]] - start[axisorder[d]]);
632 for( ; tnav.hasMore(); tnav++ )
635 copyLine(tnav.begin(), tnav.end(), acc, tmpline.begin(), acc );
637 convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
638 destIter( tnav.begin() + lstart, acc ),
639 kernel1d( kit[axisorder[d]] ), lstart, lstop);
642 dstart[axisorder[d]] = lstart;
643 dstop[axisorder[d]] = lstop;
646 copyMultiArray(tmp.traverser_begin()+dstart, stop-start, acc, di, dest);
652 scaleKernel(K & kernel,
double a)
654 for(
int i = kernel.left(); i <= kernel.right(); ++i)
655 kernel[i] = detail::RequiresExplicitCast<typename K::value_type>::cast(kernel[i] * a);
787 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
788 class DestIterator,
class DestAccessor,
class KernelIterator>
791 DestIterator d, DestAccessor dest,
792 KernelIterator kernels,
793 SrcShape
const & start = SrcShape(),
794 SrcShape
const & stop = SrcShape())
796 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
798 if(stop != SrcShape())
800 enum { N = 1 + SrcIterator::level };
801 for(
int k=0; k<N; ++k)
802 vigra_precondition(0 <= start[k] && start[k] < stop[k] && stop[k] <= shape[k],
803 "separableConvolveMultiArray(): invalid subarray shape.");
805 detail::internalSeparableConvolveSubarray(s, shape, src, d, dest, kernels, start, stop);
807 else if(!IsSameType<TmpType, typename DestAccessor::value_type>::boolResult)
810 MultiArray<SrcShape::static_size, TmpType> tmpArray(shape);
811 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src,
812 tmpArray.traverser_begin(),
typename AccessorTraits<TmpType>::default_accessor(), kernels );
818 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, d, dest, kernels );
822 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
823 class DestIterator,
class DestAccessor,
class KernelIterator>
826 triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
827 pair<DestIterator, DestAccessor>
const & dest,
829 SrcShape
const & start = SrcShape(),
830 SrcShape
const & stop = SrcShape())
833 dest.first, dest.second, kit, start, stop );
836 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
837 class DestIterator,
class DestAccessor,
class T>
840 DestIterator d, DestAccessor dest,
841 Kernel1D<T>
const & kernel,
842 SrcShape
const & start = SrcShape(),
843 SrcShape
const & stop = SrcShape())
845 ArrayVector<Kernel1D<T> > kernels(shape.size(), kernel);
850 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
851 class DestIterator,
class DestAccessor,
class T>
854 pair<DestIterator, DestAccessor>
const & dest,
855 Kernel1D<T>
const & kernel,
856 SrcShape
const & start = SrcShape(),
857 SrcShape
const & stop = SrcShape())
859 ArrayVector<Kernel1D<T> > kernels(source.second.size(), kernel);
862 dest.first, dest.second, kernels.begin(), start, stop);
937 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
938 class DestIterator,
class DestAccessor,
class T>
941 DestIterator d, DestAccessor dest,
943 SrcShape
const & start = SrcShape(),
944 SrcShape
const & stop = SrcShape())
946 enum { N = 1 + SrcIterator::level };
947 vigra_precondition( dim < N,
948 "convolveMultiArrayOneDimension(): The dimension number to convolve must be smaller "
949 "than the data dimensionality" );
951 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
952 typedef typename AccessorTraits<TmpType>::default_const_accessor TmpAccessor;
953 ArrayVector<TmpType> tmp( shape[dim] );
955 typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
956 typedef MultiArrayNavigator<DestIterator, N> DNavigator;
958 SrcShape sstart, sstop(shape), dstart, dstop(shape);
960 if(stop != SrcShape())
965 sstop[dim] = shape[dim];
966 dstop = stop - start;
969 SNavigator snav( s, sstart, sstop, dim );
970 DNavigator dnav( d, dstart, dstop, dim );
972 for( ; snav.hasMore(); snav++, dnav++ )
975 copyLine( snav.begin(), snav.end(), src,
976 tmp.begin(),
typename AccessorTraits<TmpType>::default_accessor() );
978 convolveLine(srcIterRange( tmp.begin(), tmp.end(), TmpAccessor()),
979 destIter( dnav.begin(), dest ),
980 kernel1d( kernel), start[dim], stop[dim]);
984 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
985 class DestIterator,
class DestAccessor,
class T>
988 pair<DestIterator, DestAccessor>
const & dest,
990 SrcShape
const & start = SrcShape(),
991 SrcShape
const & stop = SrcShape())
994 dest.first, dest.second, dim, kernel, start, stop);
1076 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1077 class DestIterator,
class DestAccessor>
1080 DestIterator d, DestAccessor dest,
1081 const ConvolutionOptions<SrcShape::static_size> & opt,
1082 const char *
const function_name =
"gaussianSmoothMultiArray" )
1084 typedef typename DestAccessor::value_type DestType;
1086 static const int N = SrcShape::static_size;
1088 typename ConvolutionOptions<N>::ScaleIterator params = opt.scaleParams();
1089 ArrayVector<Kernel1D<double> > kernels(N);
1091 for (
int dim = 0; dim < N; ++dim, ++params)
1092 kernels[dim].initGaussian(params.sigma_scaled(function_name), 1.0, opt.window_ratio);
1097 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1098 class DestIterator,
class DestAccessor>
1101 DestIterator d, DestAccessor dest,
double sigma,
1102 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1104 ConvolutionOptions<SrcShape::static_size> par = opt;
1108 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1109 class DestIterator,
class DestAccessor>
1112 pair<DestIterator, DestAccessor>
const & dest,
1113 const ConvolutionOptions<SrcShape::static_size> & opt)
1116 dest.first, dest.second, opt );
1119 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1120 class DestIterator,
class DestAccessor>
1123 pair<DestIterator, DestAccessor>
const & dest,
double sigma,
1124 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1127 dest.first, dest.second, sigma, opt );
1220 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1221 class DestIterator,
class DestAccessor>
1224 DestIterator di, DestAccessor dest,
1225 ConvolutionOptions<SrcShape::static_size>
const & opt,
1226 const char *
const function_name =
"gaussianGradientMultiArray")
1228 typedef typename DestAccessor::value_type DestType;
1229 typedef typename DestType::value_type DestValueType;
1230 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1232 static const int N = SrcShape::static_size;
1233 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
1235 for(
int k=0; k<N; ++k)
1239 vigra_precondition(N == (
int)dest.size(di),
1240 "gaussianGradientMultiArray(): Wrong number of channels in output array.");
1242 ParamType params = opt.scaleParams();
1243 ParamType params2(params);
1245 ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
1246 for (
int dim = 0; dim < N; ++dim, ++params)
1248 double sigma = params.sigma_scaled(function_name);
1249 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
1252 typedef VectorElementAccessor<DestAccessor> ElementAccessor;
1255 for (
int dim = 0; dim < N; ++dim, ++params2)
1257 ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
1258 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 1, 1.0, opt.window_ratio);
1259 detail::scaleKernel(kernels[dim], 1.0 / params2.step_size());
1261 opt.from_point, opt.to_point);
1265 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1266 class DestIterator,
class DestAccessor>
1269 DestIterator di, DestAccessor dest,
double sigma,
1270 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1272 ConvolutionOptions<SrcShape::static_size> par = opt;
1276 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1277 class DestIterator,
class DestAccessor>
1280 pair<DestIterator, DestAccessor>
const & dest,
1281 ConvolutionOptions<SrcShape::static_size>
const & opt )
1284 dest.first, dest.second, opt );
1287 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1288 class DestIterator,
class DestAccessor>
1291 pair<DestIterator, DestAccessor>
const & dest,
1293 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1296 dest.first, dest.second, sigma, opt );
1386 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1387 class DestIterator,
class DestAccessor>
1390 DestIterator di, DestAccessor dest,
1391 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1393 typedef typename DestAccessor::value_type DestType;
1394 typedef typename DestType::value_type DestValueType;
1395 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1397 static const int N = SrcShape::static_size;
1398 typedef typename ConvolutionOptions<N>::StepIterator StepType;
1400 for(
int k=0; k<N; ++k)
1404 vigra_precondition(N == (
int)dest.size(di),
1405 "symmetricGradientMultiArray(): Wrong number of channels in output array.");
1407 Kernel1D<KernelType> filter;
1408 filter.initSymmetricGradient();
1410 StepType step_size_it = opt.stepParams();
1412 typedef VectorElementAccessor<DestAccessor> ElementAccessor;
1415 for (
int d = 0; d < N; ++d, ++step_size_it)
1417 Kernel1D<KernelType> symmetric(filter);
1418 detail::scaleKernel(symmetric, 1 / *step_size_it);
1420 di, ElementAccessor(d, dest),
1421 d, symmetric, opt.from_point, opt.to_point);
1425 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1426 class DestIterator,
class DestAccessor>
1429 pair<DestIterator, DestAccessor>
const & dest,
1430 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1433 dest.first, dest.second, opt);
1521 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1522 class DestIterator,
class DestAccessor>
1525 DestIterator di, DestAccessor dest,
1526 ConvolutionOptions<SrcShape::static_size>
const & opt )
1528 using namespace functor;
1530 typedef typename DestAccessor::value_type DestType;
1531 typedef typename NumericTraits<DestType>::RealPromote KernelType;
1532 typedef typename AccessorTraits<KernelType>::default_accessor DerivativeAccessor;
1534 static const int N = SrcShape::static_size;
1535 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
1537 ParamType params = opt.scaleParams();
1538 ParamType params2(params);
1540 ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
1541 for (
int dim = 0; dim < N; ++dim, ++params)
1543 double sigma = params.sigma_scaled(
"laplacianOfGaussianMultiArray");
1544 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
1547 SrcShape dshape(shape);
1548 if(opt.to_point != SrcShape())
1549 dshape = opt.to_point - opt.from_point;
1551 MultiArray<N, KernelType> derivative(dshape);
1554 for (
int dim = 0; dim < N; ++dim, ++params2)
1556 ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
1557 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 2, 1.0, opt.window_ratio);
1558 detail::scaleKernel(kernels[dim], 1.0 /
sq(params2.step_size()));
1563 di, dest, kernels.begin(), opt.from_point, opt.to_point);
1568 derivative.traverser_begin(), DerivativeAccessor(),
1569 kernels.begin(), opt.from_point, opt.to_point);
1571 di, dest, Arg1() + Arg2() );
1576 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1577 class DestIterator,
class DestAccessor>
1580 DestIterator di, DestAccessor dest,
double sigma,
1581 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1583 ConvolutionOptions<SrcShape::static_size> par = opt;
1587 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1588 class DestIterator,
class DestAccessor>
1591 pair<DestIterator, DestAccessor>
const & dest,
1592 ConvolutionOptions<SrcShape::static_size>
const & opt )
1595 dest.first, dest.second, opt );
1598 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1599 class DestIterator,
class DestAccessor>
1602 pair<DestIterator, DestAccessor>
const & dest,
1604 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1607 dest.first, dest.second, sigma, opt );
1696 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1697 class DestIterator,
class DestAccessor>
1700 DestIterator di, DestAccessor dest,
1701 ConvolutionOptions<SrcShape::static_size>
const & opt )
1703 typedef typename DestAccessor::value_type DestType;
1704 typedef typename DestType::value_type DestValueType;
1705 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1707 static const int N = SrcShape::static_size;
1708 static const int M = N*(N+1)/2;
1709 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
1711 for(
int k=0; k<N; ++k)
1715 vigra_precondition(M == (
int)dest.size(di),
1716 "hessianOfGaussianMultiArray(): Wrong number of channels in output array.");
1718 ParamType params_init = opt.scaleParams();
1720 ArrayVector<Kernel1D<KernelType> > plain_kernels(N);
1721 ParamType params(params_init);
1722 for (
int dim = 0; dim < N; ++dim, ++params)
1724 double sigma = params.sigma_scaled(
"hessianOfGaussianMultiArray");
1725 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
1728 typedef VectorElementAccessor<DestAccessor> ElementAccessor;
1731 ParamType params_i(params_init);
1732 for (
int b=0, i=0; i<N; ++i, ++params_i)
1734 ParamType params_j(params_i);
1735 for (
int j=i; j<N; ++j, ++b, ++params_j)
1737 ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels);
1740 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 2, 1.0, opt.window_ratio);
1744 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 1, 1.0, opt.window_ratio);
1745 kernels[j].initGaussianDerivative(params_j.sigma_scaled(), 1, 1.0, opt.window_ratio);
1747 detail::scaleKernel(kernels[i], 1 / params_i.step_size());
1748 detail::scaleKernel(kernels[j], 1 / params_j.step_size());
1750 kernels.begin(), opt.from_point, opt.to_point);
1755 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1756 class DestIterator,
class DestAccessor>
1759 DestIterator di, DestAccessor dest,
double sigma,
1760 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1762 ConvolutionOptions<SrcShape::static_size> par = opt;
1766 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1767 class DestIterator,
class DestAccessor>
1770 pair<DestIterator, DestAccessor>
const & dest,
1771 ConvolutionOptions<SrcShape::static_size>
const & opt )
1774 dest.first, dest.second, opt );
1777 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1778 class DestIterator,
class DestAccessor>
1781 pair<DestIterator, DestAccessor>
const & dest,
1783 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1786 dest.first, dest.second, sigma, opt );
1791 template<
int N,
class VectorType>
1792 struct StructurTensorFunctor
1794 typedef VectorType result_type;
1795 typedef typename VectorType::value_type ValueType;
1798 VectorType operator()(T
const & in)
const
1801 for(
int b=0, i=0; i<N; ++i)
1803 for(
int j=i; j<N; ++j, ++b)
1805 res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
1906 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1907 class DestIterator,
class DestAccessor>
1910 DestIterator di, DestAccessor dest,
1911 ConvolutionOptions<SrcShape::static_size>
const & opt)
1913 static const int N = SrcShape::static_size;
1914 static const int M = N*(N+1)/2;
1916 typedef typename DestAccessor::value_type DestType;
1917 typedef typename DestType::value_type DestValueType;
1918 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1919 typedef TinyVector<KernelType, N> GradientVector;
1920 typedef typename AccessorTraits<GradientVector>::default_accessor GradientAccessor;
1921 typedef typename AccessorTraits<DestType>::default_accessor GradientTensorAccessor;
1923 for(
int k=0; k<N; ++k)
1927 vigra_precondition(M == (
int)dest.size(di),
1928 "structureTensorMultiArray(): Wrong number of channels in output array.");
1930 ConvolutionOptions<N> innerOptions = opt;
1931 ConvolutionOptions<N> outerOptions = opt.outerOptions();
1932 typename ConvolutionOptions<N>::ScaleIterator params = outerOptions.scaleParams();
1934 SrcShape gradientShape(shape);
1935 if(opt.to_point != SrcShape())
1937 for(
int k=0; k<N; ++k, ++params)
1939 Kernel1D<double> gauss;
1940 gauss.initGaussian(params.sigma_scaled(
"structureTensorMultiArray"), 1.0, opt.window_ratio);
1941 int dilation = gauss.right();
1942 innerOptions.from_point[k] = std::max<MultiArrayIndex>(0, opt.from_point[k] - dilation);
1943 innerOptions.to_point[k] = std::min<MultiArrayIndex>(shape[k], opt.to_point[k] + dilation);
1945 outerOptions.from_point -= innerOptions.from_point;
1946 outerOptions.to_point -= innerOptions.from_point;
1947 gradientShape = innerOptions.to_point - innerOptions.from_point;
1950 MultiArray<N, GradientVector> gradient(gradientShape);
1951 MultiArray<N, DestType> gradientTensor(gradientShape);
1953 gradient.traverser_begin(), GradientAccessor(),
1955 "structureTensorMultiArray");
1958 gradientTensor.traverser_begin(), GradientTensorAccessor(),
1959 detail::StructurTensorFunctor<N, DestType>());
1962 di, dest, outerOptions,
1963 "structureTensorMultiArray");
1966 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1967 class DestIterator,
class DestAccessor>
1970 DestIterator di, DestAccessor dest,
1971 double innerScale,
double outerScale,
1972 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
1974 ConvolutionOptions<SrcShape::static_size> par = opt;
1976 par.stdDev(innerScale).outerScale(outerScale));
1979 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1980 class DestIterator,
class DestAccessor>
1983 pair<DestIterator, DestAccessor>
const & dest,
1984 ConvolutionOptions<SrcShape::static_size>
const & opt )
1987 dest.first, dest.second, opt );
1991 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1992 class DestIterator,
class DestAccessor>
1995 pair<DestIterator, DestAccessor>
const & dest,
1996 double innerScale,
double outerScale,
1997 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>())
2000 dest.first, dest.second,
2001 innerScale, outerScale, opt);
2009 #endif //-- VIGRA_MULTI_CONVOLUTION_H
ConvolutionOptions< dim > & outerScale(...)
Generic 1 dimensional convolution kernel.
Definition: separableconvolution.hxx:1267
void indexSort(Iterator first, Iterator last, IndexIterator index_first, Compare c)
Definition: algorithm.hxx:345
void gaussianSmoothMultiArray(...)
Isotropic Gaussian smoothing of a multi-dimensional arrays.
void convolveMultiArrayOneDimension(...)
Convolution along a single dimension of a multi-dimensional arrays.
void laplacianOfGaussianMultiArray(...)
Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters.
ConvolutionOptions< dim > & subarray(Shape const &from, Shape const &to)
Definition: multi_convolution.hxx:480
void convolveLine(...)
Performs a 1-dimensional convolution of the source signal using the given kernel. ...
void gaussianGradientMultiArray(...)
Calculate Gaussian gradient of a multi-dimensional arrays.
ConvolutionOptions< dim > & stepSize(...)
Definition: multi_iterator.hxx:354
NumericTraits< T >::Promote sq(T t)
Definition: mathutil.hxx:341
void symmetricGradientMultiArray(...)
Calculate gradient of a multi-dimensional arrays using symmetric difference filters.
ConvolutionOptions< dim > & filterWindowSize(double ratio)
Definition: multi_convolution.hxx:464
Options class template for convolutions.
Definition: multi_convolution.hxx:323
void copyMultiArray(...)
Copy a multi-dimensional array.
void outer(const MultiArrayView< 2, T, C1 > &x, const MultiArrayView< 2, T, C2 > &y, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1457
void hessianOfGaussianMultiArray(...)
Calculate Hessian matrix of a N-dimensional arrays using Gaussian derivative filters.
void combineTwoMultiArrays(...)
Combine two multi-dimensional arrays into one using a binary function or functor. ...
void transformMultiArray(...)
Transform a multi-dimensional array with a unary function or functor.
ConvolutionOptions< dim > & innerScale(...)
ConvolutionOptions< dim > & resolutionStdDev(...)
void structureTensorMultiArray(...)
Calculate th structure tensor of a multi-dimensional arrays.
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
ConvolutionOptions< dim > & stdDev(...)
void separableConvolveMultiArray(...)
Separated convolution on multi-dimensional arrays.