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" 63 DoubleYielder(
double v,
unsigned,
const char *
const) : value(v) {}
64 DoubleYielder(
double v) : value(v) {}
66 double operator*()
const {
return value; }
70 struct IteratorDoubleYielder
73 IteratorDoubleYielder(X i,
unsigned,
const char *
const) : it(i) {}
74 IteratorDoubleYielder(X i) : it(i) {}
75 void operator++() { ++it; }
76 double operator*()
const {
return *it; }
80 struct SequenceDoubleYielder
82 typename X::const_iterator it;
83 SequenceDoubleYielder(
const X & seq,
unsigned dim,
84 const char *
const function_name =
"SequenceDoubleYielder")
87 if (seq.size() == dim)
89 std::string msg =
"(): Parameter number be equal to the number of spatial dimensions.";
90 vigra_precondition(
false, function_name + msg);
92 void operator++() { ++it; }
93 double operator*()
const {
return *it; }
97 struct WrapDoubleIterator
100 typename IfBool< IsConvertibleTo<X, double>::value,
102 typename IfBool< IsIterator<X>::value || IsArray<X>::value,
103 IteratorDoubleYielder<X>,
104 SequenceDoubleYielder<X>
109 template <
class Param1,
class Param2,
class Param3>
110 struct WrapDoubleIteratorTriple
112 typename WrapDoubleIterator<Param1>::type sigma_eff_it;
113 typename WrapDoubleIterator<Param2>::type sigma_d_it;
114 typename WrapDoubleIterator<Param3>::type step_size_it;
115 WrapDoubleIteratorTriple(Param1 sigma_eff, Param2 sigma_d, Param3 step_size)
116 : sigma_eff_it(sigma_eff), sigma_d_it(sigma_d), step_size_it(step_size) {}
123 double sigma_eff()
const {
return *sigma_eff_it; }
124 double sigma_d()
const {
return *sigma_d_it; }
125 double step_size()
const {
return *step_size_it; }
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 =
sq(sigma_eff()) -
sq(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); \ 323 template <
unsigned dim>
328 typedef detail::multiArrayScaleParam<dim> ParamVec;
329 typedef typename ParamVec::return_type ParamIt;
334 ParamVec outer_scale;
336 Shape from_point, to_point;
346 typedef typename detail::WrapDoubleIteratorTriple<ParamIt, ParamIt, ParamIt>
348 typedef typename detail::WrapDoubleIterator<ParamIt>::type
351 ScaleIterator scaleParams()
const 353 return ScaleIterator(sigma_eff(), sigma_d(), step_size());
355 StepIterator stepParams()
const 357 return StepIterator(step_size());
364 return outer.
stdDev(outer_scale()).resolutionStdDev(0.0);
369 VIGRA_CONVOLUTION_OPTIONS(stepSize, 1.0, step_size)
388 VIGRA_CONVOLUTION_OPTIONS(resolutionStdDev, 0.0, sigma_d)
404 VIGRA_CONVOLUTION_OPTIONS(stdDev, 0.0, sigma_eff)
405 VIGRA_CONVOLUTION_OPTIONS(innerScale, 0.0, sigma_eff)
435 VIGRA_CONVOLUTION_OPTIONS(outerScale, 0.0, outer_scale)
467 vigra_precondition(ratio >= 0.0,
468 "ConvolutionOptions::filterWindowSize(): ratio must not be negative.");
469 window_ratio = ratio;
500 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
501 class DestIterator,
class DestAccessor,
class KernelIterator>
503 internalSeparableConvolveMultiArrayTmp(
504 SrcIterator si, SrcShape
const & shape, SrcAccessor src,
505 DestIterator di, DestAccessor dest, KernelIterator kit)
507 enum { N = 1 + SrcIterator::level };
509 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
522 SNavigator snav( si, shape, 0 );
523 DNavigator dnav( di, shape, 0 );
525 for( ; snav.hasMore(); snav++, dnav++ )
528 copyLine(snav.begin(), snav.end(), src, tmp.
begin(), acc);
531 destIter( dnav.begin(), dest ),
538 for(
int d = 1; d < N; ++d, ++kit )
540 DNavigator dnav( di, shape, d );
542 tmp.resize( shape[d] );
544 for( ; dnav.hasMore(); dnav++ )
547 copyLine(dnav.begin(), dnav.end(), dest, tmp.
begin(), acc);
550 destIter( dnav.begin(), dest ),
562 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
563 class DestIterator,
class DestAccessor,
class KernelIterator>
565 internalSeparableConvolveSubarray(
566 SrcIterator si, SrcShape
const & shape, SrcAccessor src,
567 DestIterator di, DestAccessor dest, KernelIterator kit,
568 SrcShape
const & start, SrcShape
const & stop)
570 enum { N = 1 + SrcIterator::level };
572 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
574 typedef typename TmpArray::traverser TmpIterator;
577 SrcShape sstart, sstop, axisorder, tmpshape;
579 for(
int k=0; k<N; ++k)
581 sstart[k] = start[k] - kit[k].right();
584 sstop[k] = stop[k] - kit[k].left();
585 if(sstop[k] > shape[k])
587 overhead[k] = double(sstop[k] - sstart[k]) / (stop[k] - start[k]);
590 indexSort(overhead.
begin(), overhead.
end(), axisorder.begin(), std::greater<double>());
592 SrcShape dstart, dstop(sstop - sstart);
593 dstop[axisorder[0]] = stop[axisorder[0]] - start[axisorder[0]];
605 SNavigator snav( si, sstart, sstop, axisorder[0]);
606 TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[0]);
610 int lstart = start[axisorder[0]] - sstart[axisorder[0]];
611 int lstop = lstart + (stop[axisorder[0]] - start[axisorder[0]]);
613 for( ; snav.hasMore(); snav++, tnav++ )
616 copyLine(snav.begin(), snav.end(), src, tmpline.begin(), acc);
618 convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
619 destIter(tnav.begin(), acc),
620 kernel1d( kit[axisorder[0]] ), lstart, lstop);
625 for(
int d = 1; d < N; ++d)
627 TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[d]);
631 int lstart = start[axisorder[d]] - sstart[axisorder[d]];
632 int lstop = lstart + (stop[axisorder[d]] - start[axisorder[d]]);
634 for( ; tnav.hasMore(); tnav++ )
637 copyLine(tnav.begin(), tnav.end(), acc, tmpline.begin(), acc );
639 convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
640 destIter( tnav.begin() + lstart, acc ),
641 kernel1d( kit[axisorder[d]] ), lstart, lstop);
644 dstart[axisorder[d]] = lstart;
645 dstop[axisorder[d]] = lstop;
648 copyMultiArray(tmp.traverser_begin()+dstart, stop-start, acc, di, dest);
654 scaleKernel(K & kernel,
double a)
656 for(
int i = kernel.left(); i <= kernel.right(); ++i)
657 kernel[i] = detail::RequiresExplicitCast<typename K::value_type>::cast(kernel[i] * a);
845 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
846 class DestIterator,
class DestAccessor,
class KernelIterator>
848 separableConvolveMultiArray( SrcIterator s, SrcShape
const & shape, SrcAccessor src,
849 DestIterator d, DestAccessor dest,
850 KernelIterator kernels,
851 SrcShape start = SrcShape(),
852 SrcShape stop = SrcShape())
854 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
857 if(stop != SrcShape())
860 enum { N = 1 + SrcIterator::level };
861 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, start);
862 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, stop);
864 for(
int k=0; k<N; ++k)
865 vigra_precondition(0 <= start[k] && start[k] < stop[k] && stop[k] <= shape[k],
866 "separableConvolveMultiArray(): invalid subarray shape.");
868 detail::internalSeparableConvolveSubarray(s, shape, src, d, dest, kernels, start, stop);
870 else if(!IsSameType<TmpType, typename DestAccessor::value_type>::boolResult)
874 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src,
881 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, d, dest, kernels );
885 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
886 class DestIterator,
class DestAccessor,
class T>
888 separableConvolveMultiArray( SrcIterator s, SrcShape
const & shape, SrcAccessor src,
889 DestIterator d, DestAccessor dest,
891 SrcShape
const & start = SrcShape(),
892 SrcShape
const & stop = SrcShape())
896 separableConvolveMultiArray( s, shape, src, d, dest, kernels.begin(), start, stop);
899 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
900 class DestIterator,
class DestAccessor,
class KernelIterator>
902 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor>
const &
source,
903 pair<DestIterator, DestAccessor>
const & dest,
905 SrcShape
const & start = SrcShape(),
906 SrcShape
const & stop = SrcShape())
908 separableConvolveMultiArray( source.first, source.second, source.third,
909 dest.first, dest.second, kit, start, stop );
912 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
913 class DestIterator,
class DestAccessor,
class T>
915 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
916 pair<DestIterator, DestAccessor>
const & dest,
918 SrcShape
const & start = SrcShape(),
919 SrcShape
const & stop = SrcShape())
923 separableConvolveMultiArray( source.first, source.second, source.third,
924 dest.first, dest.second, kernels.begin(), start, stop);
927 template <
unsigned int N,
class T1,
class S1,
929 class KernelIterator>
939 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), start);
940 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), stop);
941 vigra_precondition(dest.
shape() == (stop - start),
942 "separableConvolveMultiArray(): shape mismatch between ROI and output.");
946 vigra_precondition(source.
shape() == dest.
shape(),
947 "separableConvolveMultiArray(): shape mismatch between input and output.");
949 separableConvolveMultiArray( srcMultiArrayRange(source),
950 destMultiArray(dest), kit, start, stop );
953 template <
unsigned int N,
class T1,
class S1,
964 separableConvolveMultiArray(source, dest, kernels.
begin(), start, stop);
1056 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1057 class DestIterator,
class DestAccessor,
class T>
1059 convolveMultiArrayOneDimension(SrcIterator s, SrcShape
const & shape, SrcAccessor src,
1060 DestIterator d, DestAccessor dest,
1062 SrcShape
const & start = SrcShape(),
1063 SrcShape
const & stop = SrcShape())
1065 enum { N = 1 + SrcIterator::level };
1066 vigra_precondition( dim < N,
1067 "convolveMultiArrayOneDimension(): The dimension number to convolve must be smaller " 1068 "than the data dimensionality" );
1070 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
1077 SrcShape sstart, sstop(shape), dstart, dstop(shape);
1079 if(stop != SrcShape())
1084 sstop[dim] = shape[dim];
1085 dstop = stop - start;
1088 SNavigator snav( s, sstart, sstop, dim );
1089 DNavigator dnav( d, dstart, dstop, dim );
1091 for( ; snav.hasMore(); snav++, dnav++ )
1094 copyLine(snav.begin(), snav.end(), src,
1098 destIter( dnav.begin(), dest ),
1099 kernel1d( kernel), start[dim], stop[dim]);
1103 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1104 class DestIterator,
class DestAccessor,
class T>
1106 convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
1107 pair<DestIterator, DestAccessor>
const & dest,
1110 SrcShape
const & start = SrcShape(),
1111 SrcShape
const & stop = SrcShape())
1113 convolveMultiArrayOneDimension(source.first, source.second, source.third,
1114 dest.first, dest.second, dim, kernel, start, stop);
1117 template <
unsigned int N,
class T1,
class S1,
1130 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), start);
1131 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), stop);
1132 vigra_precondition(dest.
shape() == (stop - start),
1133 "convolveMultiArrayOneDimension(): shape mismatch between ROI and output.");
1137 vigra_precondition(source.
shape() == dest.
shape(),
1138 "convolveMultiArrayOneDimension(): shape mismatch between input and output.");
1140 convolveMultiArrayOneDimension(srcMultiArrayRange(source),
1141 destMultiArray(dest), dim, kernel, start, stop);
1247 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1248 class DestIterator,
class DestAccessor>
1250 gaussianSmoothMultiArray( SrcIterator s, SrcShape
const & shape, SrcAccessor src,
1251 DestIterator d, DestAccessor dest,
1253 const char *
const function_name =
"gaussianSmoothMultiArray" )
1255 static const int N = SrcShape::static_size;
1257 typename ConvolutionOptions<N>::ScaleIterator params = opt.scaleParams();
1260 for (
int dim = 0; dim < N; ++dim, ++params)
1261 kernels[dim].initGaussian(params.sigma_scaled(function_name), 1.0, opt.window_ratio);
1263 separableConvolveMultiArray(s, shape, src, d, dest, kernels.
begin(), opt.from_point, opt.to_point);
1266 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1267 class DestIterator,
class DestAccessor>
1269 gaussianSmoothMultiArray( SrcIterator s, SrcShape
const & shape, SrcAccessor src,
1270 DestIterator d, DestAccessor dest,
double sigma,
1274 gaussianSmoothMultiArray(s, shape, src, d, dest, par.
stdDev(sigma));
1277 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1278 class DestIterator,
class DestAccessor>
1280 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
1281 pair<DestIterator, DestAccessor>
const & dest,
1284 gaussianSmoothMultiArray( source.first, source.second, source.third,
1285 dest.first, dest.second, opt );
1288 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1289 class DestIterator,
class DestAccessor>
1291 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
1292 pair<DestIterator, DestAccessor>
const & dest,
double sigma,
1295 gaussianSmoothMultiArray( source.first, source.second, source.third,
1296 dest.first, dest.second, sigma, opt );
1299 template <
unsigned int N,
class T1,
class S1,
1308 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
1309 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
1310 vigra_precondition(dest.
shape() == (opt.to_point - opt.from_point),
1311 "gaussianSmoothMultiArray(): shape mismatch between ROI and output.");
1315 vigra_precondition(source.
shape() == dest.
shape(),
1316 "gaussianSmoothMultiArray(): shape mismatch between input and output.");
1319 gaussianSmoothMultiArray( srcMultiArrayRange(source),
1320 destMultiArray(dest), opt );
1323 template <
unsigned int N,
class T1,
class S1,
1331 gaussianSmoothMultiArray( source, dest, opt.
stdDev(sigma) );
1438 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1439 class DestIterator,
class DestAccessor>
1441 gaussianGradientMultiArray(SrcIterator si, SrcShape
const & shape, SrcAccessor src,
1442 DestIterator di, DestAccessor dest,
1444 const char *
const function_name =
"gaussianGradientMultiArray")
1446 typedef typename DestAccessor::value_type DestType;
1447 typedef typename DestType::value_type DestValueType;
1448 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1450 static const int N = SrcShape::static_size;
1451 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
1453 for(
int k=0; k<N; ++k)
1457 vigra_precondition(N == (
int)dest.size(di),
1458 "gaussianGradientMultiArray(): Wrong number of channels in output array.");
1460 ParamType params = opt.scaleParams();
1461 ParamType params2(params);
1464 for (
int dim = 0; dim < N; ++dim, ++params)
1466 double sigma = params.sigma_scaled(function_name);
1467 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
1473 for (
int dim = 0; dim < N; ++dim, ++params2)
1476 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 1, 1.0, opt.window_ratio);
1477 detail::scaleKernel(kernels[dim], 1.0 / params2.step_size());
1478 separableConvolveMultiArray(si, shape, src, di, ElementAccessor(dim, dest), kernels.
begin(),
1479 opt.from_point, opt.to_point);
1483 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1484 class DestIterator,
class DestAccessor>
1486 gaussianGradientMultiArray(SrcIterator si, SrcShape
const & shape, SrcAccessor src,
1487 DestIterator di, DestAccessor dest,
double sigma,
1490 gaussianGradientMultiArray(si, shape, src, di, dest, opt.
stdDev(sigma));
1493 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1494 class DestIterator,
class DestAccessor>
1496 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
1497 pair<DestIterator, DestAccessor>
const & dest,
1500 gaussianGradientMultiArray( source.first, source.second, source.third,
1501 dest.first, dest.second, opt );
1504 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1505 class DestIterator,
class DestAccessor>
1507 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
1508 pair<DestIterator, DestAccessor>
const & dest,
1512 gaussianGradientMultiArray( source.first, source.second, source.third,
1513 dest.first, dest.second, sigma, opt );
1516 template <
unsigned int N,
class T1,
class S1,
1525 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
1526 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
1527 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1528 "gaussianGradientMultiArray(): shape mismatch between ROI and output.");
1532 vigra_precondition(source.
shape() == dest.shape(),
1533 "gaussianGradientMultiArray(): shape mismatch between input and output.");
1536 gaussianGradientMultiArray( srcMultiArrayRange(source),
1537 destMultiArray(dest), opt );
1540 template <
unsigned int N,
class T1,
class S1,
1548 gaussianGradientMultiArray( source, dest, opt.
stdDev(sigma) );
1559 template <
unsigned int N,
class T1,
class S1,
1569 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.from_point);
1570 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.to_point);
1571 vigra_precondition(dest.
shape() == (opt.to_point - opt.from_point),
1572 "gaussianGradientMagnitude(): shape mismatch between ROI and output.");
1576 vigra_precondition(shape == dest.
shape(),
1577 "gaussianGradientMagnitude(): shape mismatch between input and output.");
1582 typedef typename NumericTraits<T1>::RealPromote TmpType;
1585 using namespace multi_math;
1587 for(
int k=0; k<src.
shape(N); ++k)
1589 gaussianGradientMultiArray(src.
bindOuter(k), grad, opt);
1599 template <
unsigned int N,
class T1,
class S1,
1606 detail::gaussianGradientMagnitudeImpl<N, T1>(src, dest, opt);
1609 template <
unsigned int N,
class T1,
class S1,
1619 template <
unsigned int N,
class T1,
int M,
class S1,
1626 detail::gaussianGradientMagnitudeImpl<N, T1>(src.
expandElements(N), dest, opt);
1629 template <
unsigned int N,
class T1,
unsigned int R,
unsigned int G,
unsigned int B,
class S1,
1636 detail::gaussianGradientMagnitudeImpl<N, T1>(src.
expandElements(N), dest, opt);
1639 template <
unsigned int N,
class T1,
class S1,
1650 template <
unsigned int N,
class T1,
class S1,
1658 gaussianGradientMagnitude<N>(src, dest, opt.
stdDev(sigma));
1750 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1751 class DestIterator,
class DestAccessor>
1753 symmetricGradientMultiArray(SrcIterator si, SrcShape
const & shape, SrcAccessor src,
1754 DestIterator di, DestAccessor dest,
1757 typedef typename DestAccessor::value_type DestType;
1758 typedef typename DestType::value_type DestValueType;
1759 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1761 static const int N = SrcShape::static_size;
1762 typedef typename ConvolutionOptions<N>::StepIterator StepType;
1764 for(
int k=0; k<N; ++k)
1768 vigra_precondition(N == (
int)dest.size(di),
1769 "symmetricGradientMultiArray(): Wrong number of channels in output array.");
1772 filter.initSymmetricDifference();
1774 StepType step_size_it = opt.stepParams();
1779 for (
int d = 0; d < N; ++d, ++step_size_it)
1782 detail::scaleKernel(symmetric, 1 / *step_size_it);
1783 convolveMultiArrayOneDimension(si, shape, src,
1784 di, ElementAccessor(d, dest),
1785 d, symmetric, opt.from_point, opt.to_point);
1789 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1790 class DestIterator,
class DestAccessor>
1792 symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
1793 pair<DestIterator, DestAccessor>
const & dest,
1796 symmetricGradientMultiArray(source.first, source.second, source.third,
1797 dest.first, dest.second, opt);
1800 template <
unsigned int N,
class T1,
class S1,
1809 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
1810 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
1811 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1812 "symmetricGradientMultiArray(): shape mismatch between ROI and output.");
1816 vigra_precondition(source.
shape() == dest.shape(),
1817 "symmetricGradientMultiArray(): shape mismatch between input and output.");
1820 symmetricGradientMultiArray(srcMultiArrayRange(source),
1821 destMultiArray(dest), opt);
1924 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1925 class DestIterator,
class DestAccessor>
1927 laplacianOfGaussianMultiArray(SrcIterator si, SrcShape
const & shape, SrcAccessor src,
1928 DestIterator di, DestAccessor dest,
1931 using namespace functor;
1933 typedef typename DestAccessor::value_type DestType;
1934 typedef typename NumericTraits<DestType>::RealPromote KernelType;
1937 static const int N = SrcShape::static_size;
1938 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
1940 ParamType params = opt.scaleParams();
1941 ParamType params2(params);
1944 for (
int dim = 0; dim < N; ++dim, ++params)
1946 double sigma = params.sigma_scaled(
"laplacianOfGaussianMultiArray");
1947 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
1950 SrcShape dshape(shape);
1951 if(opt.to_point != SrcShape())
1952 dshape = opt.to_point - opt.from_point;
1957 for (
int dim = 0; dim < N; ++dim, ++params2)
1960 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 2, 1.0, opt.window_ratio);
1961 detail::scaleKernel(kernels[dim], 1.0 /
sq(params2.step_size()));
1965 separableConvolveMultiArray( si, shape, src,
1966 di, dest, kernels.
begin(), opt.from_point, opt.to_point);
1970 separableConvolveMultiArray( si, shape, src,
1971 derivative.traverser_begin(), DerivativeAccessor(),
1972 kernels.
begin(), opt.from_point, opt.to_point);
1974 di, dest, Arg1() + Arg2() );
1979 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1980 class DestIterator,
class DestAccessor>
1982 laplacianOfGaussianMultiArray(SrcIterator si, SrcShape
const & shape, SrcAccessor src,
1983 DestIterator di, DestAccessor dest,
double sigma,
1986 laplacianOfGaussianMultiArray(si, shape, src, di, dest, opt.
stdDev(sigma));
1989 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1990 class DestIterator,
class DestAccessor>
1992 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
1993 pair<DestIterator, DestAccessor>
const & dest,
1996 laplacianOfGaussianMultiArray( source.first, source.second, source.third,
1997 dest.first, dest.second, opt );
2000 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2001 class DestIterator,
class DestAccessor>
2003 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
2004 pair<DestIterator, DestAccessor>
const & dest,
2008 laplacianOfGaussianMultiArray( source.first, source.second, source.third,
2009 dest.first, dest.second, sigma, opt );
2012 template <
unsigned int N,
class T1,
class S1,
2021 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
2022 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
2023 vigra_precondition(dest.
shape() == (opt.to_point - opt.from_point),
2024 "laplacianOfGaussianMultiArray(): shape mismatch between ROI and output.");
2028 vigra_precondition(source.
shape() == dest.
shape(),
2029 "laplacianOfGaussianMultiArray(): shape mismatch between input and output.");
2032 laplacianOfGaussianMultiArray( srcMultiArrayRange(source),
2033 destMultiArray(dest), opt );
2036 template <
unsigned int N,
class T1,
class S1,
2044 laplacianOfGaussianMultiArray( source, dest, opt.
stdDev(sigma) );
2135 template <
class Iterator,
2136 unsigned int N,
class T,
class S>
2138 gaussianDivergenceMultiArray(Iterator vectorField, Iterator vectorFieldEnd,
2142 typedef typename std::iterator_traits<Iterator>::value_type ArrayType;
2143 typedef typename ArrayType::value_type SrcType;
2144 typedef typename NumericTraits<SrcType>::RealPromote TmpType;
2147 vigra_precondition(std::distance(vectorField, vectorFieldEnd) == N,
2148 "gaussianDivergenceMultiArray(): wrong number of input arrays.");
2151 typename ConvolutionOptions<N>::ScaleIterator params = opt.scaleParams();
2154 for(
unsigned int k = 0; k < N; ++k, ++params)
2156 sigmas[k] = params.sigma_scaled(
"gaussianDivergenceMultiArray");
2157 kernels[k].initGaussian(sigmas[k], 1.0, opt.window_ratio);
2162 for(
unsigned int k=0; k < N; ++k, ++vectorField)
2164 kernels[k].initGaussianDerivative(sigmas[k], 1, 1.0, opt.window_ratio);
2167 separableConvolveMultiArray(*vectorField, divergence, kernels.
begin(), opt.from_point, opt.to_point);
2171 separableConvolveMultiArray(*vectorField, tmpDeriv, kernels.
begin(), opt.from_point, opt.to_point);
2172 divergence += tmpDeriv;
2174 kernels[k].initGaussian(sigmas[k], 1.0, opt.window_ratio);
2178 template <
class Iterator,
2179 unsigned int N,
class T,
class S>
2181 gaussianDivergenceMultiArray(Iterator vectorField, Iterator vectorFieldEnd,
2186 gaussianDivergenceMultiArray(vectorField, vectorFieldEnd, divergence, opt.
stdDev(sigma));
2189 template <
unsigned int N,
class T1,
class S1,
2197 for(
unsigned int k=0; k<N; ++k)
2198 field.push_back(vectorField.bindElementChannel(k));
2200 gaussianDivergenceMultiArray(field.
begin(), field.
end(), divergence, opt);
2203 template <
unsigned int N,
class T1,
class S1,
2211 gaussianDivergenceMultiArray(vectorField, divergence, opt.
stdDev(sigma));
2316 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2317 class DestIterator,
class DestAccessor>
2319 hessianOfGaussianMultiArray(SrcIterator si, SrcShape
const & shape, SrcAccessor src,
2320 DestIterator di, DestAccessor dest,
2323 typedef typename DestAccessor::value_type DestType;
2324 typedef typename DestType::value_type DestValueType;
2325 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
2327 static const int N = SrcShape::static_size;
2328 static const int M = N*(N+1)/2;
2329 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType;
2331 for(
int k=0; k<N; ++k)
2335 vigra_precondition(M == (
int)dest.size(di),
2336 "hessianOfGaussianMultiArray(): Wrong number of channels in output array.");
2338 ParamType params_init = opt.scaleParams();
2341 ParamType params(params_init);
2342 for (
int dim = 0; dim < N; ++dim, ++params)
2344 double sigma = params.sigma_scaled(
"hessianOfGaussianMultiArray");
2345 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
2351 ParamType params_i(params_init);
2352 for (
int b=0, i=0; i<N; ++i, ++params_i)
2354 ParamType params_j(params_i);
2355 for (
int j=i; j<N; ++j, ++b, ++params_j)
2360 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 2, 1.0, opt.window_ratio);
2364 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 1, 1.0, opt.window_ratio);
2365 kernels[j].initGaussianDerivative(params_j.sigma_scaled(), 1, 1.0, opt.window_ratio);
2367 detail::scaleKernel(kernels[i], 1 / params_i.step_size());
2368 detail::scaleKernel(kernels[j], 1 / params_j.step_size());
2369 separableConvolveMultiArray(si, shape, src, di, ElementAccessor(b, dest),
2370 kernels.
begin(), opt.from_point, opt.to_point);
2375 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2376 class DestIterator,
class DestAccessor>
2378 hessianOfGaussianMultiArray(SrcIterator si, SrcShape
const & shape, SrcAccessor src,
2379 DestIterator di, DestAccessor dest,
double sigma,
2382 hessianOfGaussianMultiArray(si, shape, src, di, dest, opt.
stdDev(sigma));
2385 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2386 class DestIterator,
class DestAccessor>
2388 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
2389 pair<DestIterator, DestAccessor>
const & dest,
2392 hessianOfGaussianMultiArray( source.first, source.second, source.third,
2393 dest.first, dest.second, opt );
2396 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2397 class DestIterator,
class DestAccessor>
2399 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
2400 pair<DestIterator, DestAccessor>
const & dest,
2404 hessianOfGaussianMultiArray( source.first, source.second, source.third,
2405 dest.first, dest.second, sigma, opt );
2408 template <
unsigned int N,
class T1,
class S1,
2417 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
2418 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
2419 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2420 "hessianOfGaussianMultiArray(): shape mismatch between ROI and output.");
2424 vigra_precondition(source.
shape() == dest.shape(),
2425 "hessianOfGaussianMultiArray(): shape mismatch between input and output.");
2428 hessianOfGaussianMultiArray( srcMultiArrayRange(source),
2429 destMultiArray(dest), opt );
2432 template <
unsigned int N,
class T1,
class S1,
2440 hessianOfGaussianMultiArray( source, dest, opt.
stdDev(sigma) );
2445 template<
int N,
class VectorType>
2446 struct StructurTensorFunctor
2448 typedef VectorType result_type;
2449 typedef typename VectorType::value_type ValueType;
2452 VectorType operator()(T
const & in)
const 2455 for(
int b=0, i=0; i<N; ++i)
2457 for(
int j=i; j<N; ++j, ++b)
2459 res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
2573 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2574 class DestIterator,
class DestAccessor>
2576 structureTensorMultiArray(SrcIterator si, SrcShape
const & shape, SrcAccessor src,
2577 DestIterator di, DestAccessor dest,
2580 static const int N = SrcShape::static_size;
2581 static const int M = N*(N+1)/2;
2583 typedef typename DestAccessor::value_type DestType;
2584 typedef typename DestType::value_type DestValueType;
2585 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
2590 for(
int k=0; k<N; ++k)
2594 vigra_precondition(M == (
int)dest.size(di),
2595 "structureTensorMultiArray(): Wrong number of channels in output array.");
2599 typename ConvolutionOptions<N>::ScaleIterator params = outerOptions.scaleParams();
2601 SrcShape gradientShape(shape);
2602 if(opt.to_point != SrcShape())
2604 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.from_point);
2605 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.to_point);
2607 for(
int k=0; k<N; ++k, ++params)
2610 gauss.
initGaussian(params.sigma_scaled(
"structureTensorMultiArray"), 1.0, opt.window_ratio);
2611 int dilation = gauss.
right();
2612 innerOptions.from_point[k] = std::max<MultiArrayIndex>(0, opt.from_point[k] - dilation);
2613 innerOptions.to_point[k] = std::min<MultiArrayIndex>(shape[k], opt.to_point[k] + dilation);
2615 outerOptions.from_point -= innerOptions.from_point;
2616 outerOptions.to_point -= innerOptions.from_point;
2617 gradientShape = innerOptions.to_point - innerOptions.from_point;
2622 gaussianGradientMultiArray(si, shape, src,
2623 gradient.traverser_begin(), GradientAccessor(),
2625 "structureTensorMultiArray");
2628 gradientTensor.traverser_begin(), GradientTensorAccessor(),
2629 detail::StructurTensorFunctor<N, DestType>());
2631 gaussianSmoothMultiArray(gradientTensor.traverser_begin(), gradientShape, GradientTensorAccessor(),
2632 di, dest, outerOptions,
2633 "structureTensorMultiArray");
2636 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2637 class DestIterator,
class DestAccessor>
2639 structureTensorMultiArray(SrcIterator si, SrcShape
const & shape, SrcAccessor src,
2640 DestIterator di, DestAccessor dest,
2641 double innerScale,
double outerScale,
2644 structureTensorMultiArray(si, shape, src, di, dest,
2645 opt.
stdDev(innerScale).outerScale(outerScale));
2648 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2649 class DestIterator,
class DestAccessor>
2651 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
2652 pair<DestIterator, DestAccessor>
const & dest,
2655 structureTensorMultiArray( source.first, source.second, source.third,
2656 dest.first, dest.second, opt );
2660 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2661 class DestIterator,
class DestAccessor>
2663 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
2664 pair<DestIterator, DestAccessor>
const & dest,
2665 double innerScale,
double outerScale,
2668 structureTensorMultiArray( source.first, source.second, source.third,
2669 dest.first, dest.second,
2670 innerScale, outerScale, opt);
2673 template <
unsigned int N,
class T1,
class S1,
2682 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
2683 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
2684 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2685 "structureTensorMultiArray(): shape mismatch between ROI and output.");
2689 vigra_precondition(source.
shape() == dest.shape(),
2690 "structureTensorMultiArray(): shape mismatch between input and output.");
2693 structureTensorMultiArray( srcMultiArrayRange(source),
2694 destMultiArray(dest), opt );
2698 template <
unsigned int N,
class T1,
class S1,
2703 double innerScale,
double outerScale,
2706 structureTensorMultiArray(source, dest, opt.
innerScale(innerScale).outerScale(outerScale));
2714 #endif //-- VIGRA_MULTI_CONVOLUTION_H 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:2786
Accessor for one component of a vector.
Definition: accessor.hxx:539
Generic 1 dimensional convolution kernel.
Definition: separableconvolution.hxx:52
void indexSort(Iterator first, Iterator last, IndexIterator index_first, Compare c)
Return the index permutation that would sort the input array.
Definition: algorithm.hxx:345
void initGaussian(double std_dev, value_type norm, double windowRatio=0.0)
Definition: separableconvolution.hxx:2251
const difference_type & shape() const
Definition: multi_array.hxx:1551
void gaussianSmoothMultiArray(...)
Isotropic Gaussian smoothing of a multi-dimensional arrays.
void gaussianGradientMagnitude(...)
Calculate the gradient magnitude by means of a 1st derivatives of Gaussian filter.
const_iterator begin() const
Definition: array_vector.hxx:223
void convolveMultiArrayOneDimension(...)
Convolution along a single dimension of a multi-dimensional arrays.
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
void laplacianOfGaussianMultiArray(...)
Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters.
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2357
iterator end()
Definition: tinyvector.hxx:838
ConvolutionOptions< dim > & subarray(Shape const &from, Shape const &to)
Definition: multi_convolution.hxx:483
void convolveLine(...)
Performs a 1-dimensional convolution of the source signal using the given kernel. ...
Definition: accessor.hxx:43
void gaussianDivergenceMultiArray(...)
Calculate the divergence of a vector field using Gaussian derivative filters.
MultiArrayView< N+1, T, StrideTag > insertSingletonDimension(difference_type_1 i) const
Definition: multi_array.hxx:2237
void gaussianGradientMultiArray(...)
Calculate Gaussian gradient of a multi-dimensional arrays.
Definition: array_vector.hxx:58
Definition: multi_shape.hxx:231
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:344
A navigator that provides access to the 1D subranges of an n-dimensional range given by a vigra::Mult...
Definition: navigator.hxx:97
iterator begin()
Definition: tinyvector.hxx:835
void symmetricGradientMultiArray(...)
Calculate gradient of a multi-dimensional arrays using symmetric difference filters.
ConvolutionOptions< dim > & filterWindowSize(double ratio)
Definition: multi_convolution.hxx:465
Options class template for convolutions.
Definition: multi_convolution.hxx:324
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. ...
Encapsulate read access to the values an iterator points to.
Definition: accessor.hxx:269
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:655
void transformMultiArray(...)
Transform a multi-dimensional array with a unary function or functor.
int right() const
Definition: separableconvolution.hxx:2155
ConvolutionOptions< dim > & innerScale(...)
MultiArrayView & init(const U &init)
Definition: multi_array.hxx:1150
const_iterator end() const
Definition: array_vector.hxx:237
Class for a single RGB value.
Definition: accessor.hxx:938
Encapsulate access to the values an iterator points to.
Definition: accessor.hxx:133
MultiArrayView< N+1, typename ExpandElementResult< T >::type, StridedArrayTag > expandElements(difference_type_1 d) const
Definition: multi_array.hxx:2208
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(...)
MultiArrayView< N-M, T, StrideTag > bindOuter(const TinyVector< Index, M > &d) const
Definition: multi_array.hxx:2067
void separableConvolveMultiArray(...)
Separated convolution on multi-dimensional arrays.