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; }
73 struct 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; }
83 struct 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; }
100 struct WrapDoubleIterator
103 typename IfBool< IsConvertibleTo<X, double>::value,
105 typename IfBool< IsIterator<X>::value || IsArray<X>::value,
106 IteratorDoubleYielder<X>,
107 SequenceDoubleYielder<X>
112 template <
class Param1,
class Param2,
class Param3>
113 struct 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 ")
const 139 sigma_precondition(sigma_eff(), function_name);
140 sigma_precondition(sigma_d(), function_name);
141 double sigma_squared =
sq(sigma_eff()) -
sq(sigma_d());
142 if (sigma_squared > 0.0)
144 return std::sqrt(sigma_squared) / step_size();
148 std::string msg =
"(): Scale would be imaginary or zero.";
149 vigra_precondition(
false, function_name + msg);
155 template <
unsigned dim>
156 struct multiArrayScaleParam
158 typedef TinyVector<double, dim> p_vector;
159 typedef typename p_vector::const_iterator return_type;
162 template <
class Param>
163 multiArrayScaleParam(Param val,
const char *
const function_name =
"multiArrayScaleParam")
165 typename WrapDoubleIterator<Param>::type in(val, dim, function_name);
166 for (
unsigned i = 0; i != dim; ++i, ++in)
169 return_type operator()()
const 173 static void precondition(
unsigned n_par,
const char *
const function_name =
"multiArrayScaleParam")
177 std::string msg =
"(): dimension parameter must be ";
178 vigra_precondition(dim == n_par, function_name + msg + n);
180 multiArrayScaleParam(
double v0,
double v1,
const char *
const function_name =
"multiArrayScaleParam")
182 precondition(2, function_name);
183 vec = p_vector(v0, v1);
185 multiArrayScaleParam(
double v0,
double v1,
double v2,
const char *
const function_name =
"multiArrayScaleParam")
187 precondition(3, function_name);
188 vec = p_vector(v0, v1, v2);
190 multiArrayScaleParam(
double v0,
double v1,
double v2,
double v3,
const char *
const function_name =
"multiArrayScaleParam")
192 precondition(4, function_name);
193 vec = p_vector(v0, v1, v2, v3);
195 multiArrayScaleParam(
double v0,
double v1,
double v2,
double v3,
double v4,
const char *
const function_name =
"multiArrayScaleParam")
197 precondition(5, function_name);
198 vec = p_vector(v0, v1, v2, v3, v4);
204 #define VIGRA_CONVOLUTION_OPTIONS(function_name, default_value, member_name, getter_setter_name) \ 205 template <class Param> \ 206 ConvolutionOptions & function_name(const Param & val) \ 208 member_name = ParamVec(val, "ConvolutionOptions::" #function_name); \ 211 ConvolutionOptions & function_name() \ 213 member_name = ParamVec(default_value, "ConvolutionOptions::" #function_name); \ 216 ConvolutionOptions & function_name(double v0, double v1) \ 218 member_name = ParamVec(v0, v1, "ConvolutionOptions::" #function_name); \ 221 ConvolutionOptions & function_name(double v0, double v1, double v2) \ 223 member_name = ParamVec(v0, v1, v2, "ConvolutionOptions::" #function_name); \ 226 ConvolutionOptions & function_name(double v0, double v1, double v2, double v3) \ 228 member_name = ParamVec(v0, v1, v2, v3, "ConvolutionOptions::" #function_name); \ 231 ConvolutionOptions & function_name(double v0, double v1, double v2, double v3, double v4) \ 233 member_name = ParamVec(v0, v1, v2, v3, v4, "ConvolutionOptions::" #function_name); \ 236 typename ParamVec::p_vector get##getter_setter_name()const{ \ 237 return member_name.vec; \ 239 void set##getter_setter_name(typename ParamVec::p_vector vec){ \ 240 member_name.vec = vec; \ 331 template <
unsigned dim>
336 typedef detail::multiArrayScaleParam<dim> ParamVec;
337 typedef typename ParamVec::return_type ParamIt;
342 ParamVec outer_scale;
344 Shape from_point, to_point;
354 typedef typename detail::WrapDoubleIteratorTriple<ParamIt, ParamIt, ParamIt>
356 typedef typename detail::WrapDoubleIterator<ParamIt>::type
359 ScaleIterator scaleParams()
const 361 return ScaleIterator(sigma_eff(), sigma_d(), step_size());
363 StepIterator stepParams()
const 365 return StepIterator(step_size());
372 return outer.
stdDev(outer_scale()).resolutionStdDev(0.0);
377 VIGRA_CONVOLUTION_OPTIONS(stepSize, 1.0, step_size, StepSize)
396 VIGRA_CONVOLUTION_OPTIONS(resolutionStdDev, 0.0, sigma_d, ResolutionStdDev)
412 VIGRA_CONVOLUTION_OPTIONS(stdDev, 0.0, sigma_eff, StdDev)
413 VIGRA_CONVOLUTION_OPTIONS(innerScale, 0.0, sigma_eff, InnerScale)
443 VIGRA_CONVOLUTION_OPTIONS(outerScale, 0.0, outer_scale, OuterScale)
475 vigra_precondition(ratio >= 0.0,
476 "ConvolutionOptions::filterWindowSize(): ratio must not be negative.");
477 window_ratio = ratio;
481 double getFilterWindowSize()
const {
502 std::pair<Shape, Shape> getSubarray()
const{
503 std::pair<Shape, Shape> res;
504 res.first = from_point;
505 res.second = to_point;
519 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
520 class DestIterator,
class DestAccessor,
class KernelIterator>
522 internalSeparableConvolveMultiArrayTmp(
523 SrcIterator si, SrcShape
const & shape, SrcAccessor src,
524 DestIterator di, DestAccessor dest, KernelIterator kit)
526 enum { N = 1 + SrcIterator::level };
528 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
541 SNavigator snav( si, shape, 0 );
542 DNavigator dnav( di, shape, 0 );
544 for( ; snav.hasMore(); snav++, dnav++ )
547 copyLine(snav.begin(), snav.end(), src, tmp.
begin(), acc);
550 destIter( dnav.begin(), dest ),
557 for(
int d = 1; d < N; ++d, ++kit )
559 DNavigator dnav( di, shape, d );
561 tmp.resize( shape[d] );
563 for( ; dnav.hasMore(); dnav++ )
566 copyLine(dnav.begin(), dnav.end(), dest, tmp.
begin(), acc);
569 destIter( dnav.begin(), dest ),
581 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
582 class DestIterator,
class DestAccessor,
class KernelIterator>
584 internalSeparableConvolveSubarray(
585 SrcIterator si, SrcShape
const & shape, SrcAccessor src,
586 DestIterator di, DestAccessor dest, KernelIterator kit,
587 SrcShape
const & start, SrcShape
const & stop)
589 enum { N = 1 + SrcIterator::level };
591 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
593 typedef typename TmpArray::traverser TmpIterator;
596 SrcShape sstart, sstop, axisorder, tmpshape;
598 for(
int k=0; k<N; ++k)
601 sstart[k] = start[k] - kit[k].right();
604 sstop[k] = stop[k] - kit[k].left();
605 if(sstop[k] > shape[k])
607 overhead[k] = double(sstop[k] - sstart[k]) / (stop[k] - start[k]);
610 indexSort(overhead.
begin(), overhead.
end(), axisorder.begin(), std::greater<double>());
611 SrcShape dstart, dstop(sstop - sstart);
612 dstop[axisorder[0]] = stop[axisorder[0]] - start[axisorder[0]];
624 SNavigator snav( si, sstart, sstop, axisorder[0]);
625 TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[0]);
629 int lstart = start[axisorder[0]] - sstart[axisorder[0]];
630 int lstop = lstart + (stop[axisorder[0]] - start[axisorder[0]]);
632 for( ; snav.hasMore(); snav++, tnav++ )
635 copyLine(snav.begin(), snav.end(), src, tmpline.begin(), acc);
637 convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
638 destIter(tnav.begin(), acc),
639 kernel1d( kit[axisorder[0]] ), lstart, lstop);
644 for(
int d = 1; d < N; ++d)
646 TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[d]);
650 int lstart = start[axisorder[d]] - sstart[axisorder[d]];
651 int lstop = lstart + (stop[axisorder[d]] - start[axisorder[d]]);
653 for( ; tnav.hasMore(); tnav++ )
656 copyLine(tnav.begin(), tnav.end(), acc, tmpline.begin(), acc );
658 convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
659 destIter( tnav.begin() + lstart, acc ),
660 kernel1d( kit[axisorder[d]] ), lstart, lstop);
663 dstart[axisorder[d]] = lstart;
664 dstop[axisorder[d]] = lstop;
667 copyMultiArray(tmp.traverser_begin()+dstart, stop-start, acc, di, dest);
673 scaleKernel(K & kernel,
double a)
675 for(
int i = kernel.left(); i <= kernel.right(); ++i)
676 kernel[i] = detail::RequiresExplicitCast<typename K::value_type>::cast(kernel[i] * a);
864 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
865 class DestIterator,
class DestAccessor,
class KernelIterator>
868 DestIterator d, DestAccessor dest,
869 KernelIterator kernels,
870 SrcShape start = SrcShape(),
871 SrcShape stop = SrcShape())
873 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
876 if(stop != SrcShape())
879 enum { N = 1 + SrcIterator::level };
880 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, start);
881 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, stop);
883 for(
int k=0; k<N; ++k)
884 vigra_precondition(0 <= start[k] && start[k] < stop[k] && stop[k] <= shape[k],
885 "separableConvolveMultiArray(): invalid subarray shape.");
887 detail::internalSeparableConvolveSubarray(s, shape, src, d, dest, kernels, start, stop);
889 else if(!IsSameType<TmpType, typename DestAccessor::value_type>::boolResult)
893 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src,
900 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, d, dest, kernels );
904 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
905 class DestIterator,
class DestAccessor,
class T>
908 DestIterator d, DestAccessor dest,
910 SrcShape
const & start = SrcShape(),
911 SrcShape
const & stop = SrcShape())
918 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
919 class DestIterator,
class DestAccessor,
class KernelIterator>
922 pair<DestIterator, DestAccessor>
const & dest,
924 SrcShape
const & start = SrcShape(),
925 SrcShape
const & stop = SrcShape())
928 dest.first, dest.second, kit, start, stop );
931 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
932 class DestIterator,
class DestAccessor,
class T>
935 pair<DestIterator, DestAccessor>
const & dest,
937 SrcShape
const & start = SrcShape(),
938 SrcShape
const & stop = SrcShape())
943 dest.first, dest.second, kernels.begin(), start, stop);
946 template <
unsigned int N,
class T1,
class S1,
948 class KernelIterator>
958 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), start);
959 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), stop);
960 vigra_precondition(dest.
shape() == (stop - start),
961 "separableConvolveMultiArray(): shape mismatch between ROI and output.");
965 vigra_precondition(source.
shape() == dest.
shape(),
966 "separableConvolveMultiArray(): shape mismatch between input and output.");
969 destMultiArray(dest), kit, start, stop );
972 template <
unsigned int N,
class T1,
class S1,
1075 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1076 class DestIterator,
class DestAccessor,
class T>
1079 DestIterator d, DestAccessor dest,
1081 SrcShape
const & start = SrcShape(),
1082 SrcShape
const & stop = SrcShape())
1084 enum { N = 1 + SrcIterator::level };
1085 vigra_precondition( dim < N,
1086 "convolveMultiArrayOneDimension(): The dimension number to convolve must be smaller " 1087 "than the data dimensionality" );
1089 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
1096 SrcShape sstart, sstop(shape), dstart, dstop(shape);
1098 if(stop != SrcShape())
1103 sstop[dim] = shape[dim];
1104 dstop = stop - start;
1107 SNavigator snav( s, sstart, sstop, dim );
1108 DNavigator dnav( d, dstart, dstop, dim );
1110 for( ; snav.hasMore(); snav++, dnav++ )
1113 copyLine(snav.begin(), snav.end(), src,
1117 destIter( dnav.begin(), dest ),
1118 kernel1d( kernel), start[dim], stop[dim]);
1122 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1123 class DestIterator,
class DestAccessor,
class T>
1126 pair<DestIterator, DestAccessor>
const & dest,
1129 SrcShape
const & start = SrcShape(),
1130 SrcShape
const & stop = SrcShape())
1133 dest.first, dest.second, dim, kernel, start, stop);
1136 template <
unsigned int N,
class T1,
class S1,
1149 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), start);
1150 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), stop);
1151 vigra_precondition(dest.
shape() == (stop - start),
1152 "convolveMultiArrayOneDimension(): shape mismatch between ROI and output.");
1156 vigra_precondition(source.
shape() == dest.
shape(),
1157 "convolveMultiArrayOneDimension(): shape mismatch between input and output.");
1160 destMultiArray(dest), dim, kernel, start, stop);
1297 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1298 class DestIterator,
class DestAccessor>
1301 DestIterator d, DestAccessor dest,
1303 const char *
const function_name =
"gaussianSmoothMultiArray" )
1305 static const int N = SrcShape::static_size;
1310 for (
int dim = 0; dim < N; ++dim, ++params)
1311 kernels[dim].initGaussian(params.sigma_scaled(function_name), 1.0, opt.window_ratio);
1316 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1317 class DestIterator,
class DestAccessor>
1320 DestIterator d, DestAccessor dest,
double sigma,
1327 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1328 class DestIterator,
class DestAccessor>
1331 pair<DestIterator, DestAccessor>
const & dest,
1335 dest.first, dest.second, opt );
1338 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1339 class DestIterator,
class DestAccessor>
1342 pair<DestIterator, DestAccessor>
const & dest,
double sigma,
1346 dest.first, dest.second, sigma, opt );
1349 template <
unsigned int N,
class T1,
class S1,
1358 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
1359 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
1360 vigra_precondition(dest.
shape() == (opt.to_point - opt.from_point),
1361 "gaussianSmoothMultiArray(): shape mismatch between ROI and output.");
1365 vigra_precondition(source.
shape() == dest.
shape(),
1366 "gaussianSmoothMultiArray(): shape mismatch between input and output.");
1370 destMultiArray(dest), opt );
1373 template <
unsigned int N,
class T1,
class S1,
1504 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1505 class DestIterator,
class DestAccessor>
1508 DestIterator di, DestAccessor dest,
1510 const char *
const function_name =
"gaussianGradientMultiArray")
1512 typedef typename DestAccessor::value_type DestType;
1513 typedef typename DestType::value_type DestValueType;
1514 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1516 static const int N = SrcShape::static_size;
1519 for(
int k=0; k<N; ++k)
1523 vigra_precondition(N == (
int)dest.size(di),
1524 "gaussianGradientMultiArray(): Wrong number of channels in output array.");
1526 ParamType params = opt.scaleParams();
1527 ParamType params2(params);
1530 for (
int dim = 0; dim < N; ++dim, ++params)
1532 double sigma = params.sigma_scaled(function_name);
1533 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
1539 for (
int dim = 0; dim < N; ++dim, ++params2)
1542 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 1, 1.0, opt.window_ratio);
1543 detail::scaleKernel(kernels[dim], 1.0 / params2.step_size());
1545 opt.from_point, opt.to_point);
1549 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1550 class DestIterator,
class DestAccessor>
1553 DestIterator di, DestAccessor dest,
double sigma,
1559 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1560 class DestIterator,
class DestAccessor>
1563 pair<DestIterator, DestAccessor>
const & dest,
1567 dest.first, dest.second, opt );
1570 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1571 class DestIterator,
class DestAccessor>
1574 pair<DestIterator, DestAccessor>
const & dest,
1579 dest.first, dest.second, sigma, opt );
1582 template <
unsigned int N,
class T1,
class S1,
1591 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
1592 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
1593 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1594 "gaussianGradientMultiArray(): shape mismatch between ROI and output.");
1598 vigra_precondition(source.
shape() == dest.shape(),
1599 "gaussianGradientMultiArray(): shape mismatch between input and output.");
1603 destMultiArray(dest), opt );
1606 template <
unsigned int N,
class T1,
class S1,
1625 template <
unsigned int N,
class T1,
class S1,
1635 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.from_point);
1636 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.to_point);
1637 vigra_precondition(dest.
shape() == (opt.to_point - opt.from_point),
1638 "gaussianGradientMagnitude(): shape mismatch between ROI and output.");
1642 vigra_precondition(shape == dest.
shape(),
1643 "gaussianGradientMagnitude(): shape mismatch between input and output.");
1648 typedef typename NumericTraits<T1>::RealPromote TmpType;
1651 using namespace multi_math;
1653 for(
int k=0; k<src.
shape(N); ++k)
1665 template <
unsigned int N,
class T1,
class S1,
1672 detail::gaussianGradientMagnitudeImpl<N, T1>(src, dest, opt);
1675 template <
unsigned int N,
class T1,
class S1,
1685 template <
unsigned int N,
class T1,
int M,
class S1,
1692 detail::gaussianGradientMagnitudeImpl<N, T1>(src.
expandElements(N), dest, opt);
1695 template <
unsigned int N,
class T1,
unsigned int R,
unsigned int G,
unsigned int B,
class S1,
1702 detail::gaussianGradientMagnitudeImpl<N, T1>(src.
expandElements(N), dest, opt);
1705 template <
unsigned int N,
class T1,
class S1,
1716 template <
unsigned int N,
class T1,
class S1,
1724 gaussianGradientMagnitude<N>(src, dest, opt.
stdDev(sigma));
1833 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1834 class DestIterator,
class DestAccessor>
1837 DestIterator di, DestAccessor dest,
1840 typedef typename DestAccessor::value_type DestType;
1841 typedef typename DestType::value_type DestValueType;
1842 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1844 static const int N = SrcShape::static_size;
1847 for(
int k=0; k<N; ++k)
1851 vigra_precondition(N == (
int)dest.size(di),
1852 "symmetricGradientMultiArray(): Wrong number of channels in output array.");
1855 filter.initSymmetricDifference();
1857 StepType step_size_it = opt.stepParams();
1862 for (
int d = 0; d < N; ++d, ++step_size_it)
1865 detail::scaleKernel(symmetric, 1 / *step_size_it);
1867 di, ElementAccessor(d, dest),
1868 d, symmetric, opt.from_point, opt.to_point);
1872 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1873 class DestIterator,
class DestAccessor>
1876 pair<DestIterator, DestAccessor>
const & dest,
1880 dest.first, dest.second, opt);
1883 template <
unsigned int N,
class T1,
class S1,
1892 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
1893 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
1894 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1895 "symmetricGradientMultiArray(): shape mismatch between ROI and output.");
1899 vigra_precondition(source.
shape() == dest.shape(),
1900 "symmetricGradientMultiArray(): shape mismatch between input and output.");
1904 destMultiArray(dest), opt);
2023 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2024 class DestIterator,
class DestAccessor>
2027 DestIterator di, DestAccessor dest,
2030 using namespace functor;
2032 typedef typename DestAccessor::value_type DestType;
2033 typedef typename NumericTraits<DestType>::RealPromote KernelType;
2036 static const int N = SrcShape::static_size;
2039 ParamType params = opt.scaleParams();
2040 ParamType params2(params);
2043 for (
int dim = 0; dim < N; ++dim, ++params)
2045 double sigma = params.sigma_scaled(
"laplacianOfGaussianMultiArray");
2046 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
2049 SrcShape dshape(shape);
2050 if(opt.to_point != SrcShape())
2051 dshape = opt.to_point - opt.from_point;
2056 for (
int dim = 0; dim < N; ++dim, ++params2)
2059 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 2, 1.0, opt.window_ratio);
2060 detail::scaleKernel(kernels[dim], 1.0 /
sq(params2.step_size()));
2065 di, dest, kernels.
begin(), opt.from_point, opt.to_point);
2070 derivative.traverser_begin(), DerivativeAccessor(),
2071 kernels.
begin(), opt.from_point, opt.to_point);
2073 di, dest, Arg1() + Arg2() );
2078 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2079 class DestIterator,
class DestAccessor>
2082 DestIterator di, DestAccessor dest,
double sigma,
2088 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2089 class DestIterator,
class DestAccessor>
2092 pair<DestIterator, DestAccessor>
const & dest,
2096 dest.first, dest.second, opt );
2099 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2100 class DestIterator,
class DestAccessor>
2103 pair<DestIterator, DestAccessor>
const & dest,
2108 dest.first, dest.second, sigma, opt );
2111 template <
unsigned int N,
class T1,
class S1,
2120 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
2121 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
2122 vigra_precondition(dest.
shape() == (opt.to_point - opt.from_point),
2123 "laplacianOfGaussianMultiArray(): shape mismatch between ROI and output.");
2127 vigra_precondition(source.
shape() == dest.
shape(),
2128 "laplacianOfGaussianMultiArray(): shape mismatch between input and output.");
2132 destMultiArray(dest), opt );
2135 template <
unsigned int N,
class T1,
class S1,
2251 template <
class Iterator,
2252 unsigned int N,
class T,
class S>
2258 typedef typename std::iterator_traits<Iterator>::value_type ArrayType;
2259 typedef typename ArrayType::value_type SrcType;
2260 typedef typename NumericTraits<SrcType>::RealPromote TmpType;
2263 vigra_precondition(std::distance(vectorField, vectorFieldEnd) == N,
2264 "gaussianDivergenceMultiArray(): wrong number of input arrays.");
2270 for(
unsigned int k = 0; k < N; ++k, ++params)
2272 sigmas[k] = params.sigma_scaled(
"gaussianDivergenceMultiArray");
2273 kernels[k].initGaussian(sigmas[k], 1.0, opt.window_ratio);
2278 for(
unsigned int k=0; k < N; ++k, ++vectorField)
2280 kernels[k].initGaussianDerivative(sigmas[k], 1, 1.0, opt.window_ratio);
2288 divergence += tmpDeriv;
2290 kernels[k].initGaussian(sigmas[k], 1.0, opt.window_ratio);
2294 template <
class Iterator,
2295 unsigned int N,
class T,
class S>
2305 template <
unsigned int N,
class T1,
class S1,
2313 for(
unsigned int k=0; k<N; ++k)
2314 field.push_back(vectorField.bindElementChannel(k));
2319 template <
unsigned int N,
class T1,
class S1,
2448 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2449 class DestIterator,
class DestAccessor>
2452 DestIterator di, DestAccessor dest,
2455 typedef typename DestAccessor::value_type DestType;
2456 typedef typename DestType::value_type DestValueType;
2457 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
2459 static const int N = SrcShape::static_size;
2460 static const int M = N*(N+1)/2;
2463 for(
int k=0; k<N; ++k)
2467 vigra_precondition(M == (
int)dest.size(di),
2468 "hessianOfGaussianMultiArray(): Wrong number of channels in output array.");
2470 ParamType params_init = opt.scaleParams();
2473 ParamType params(params_init);
2474 for (
int dim = 0; dim < N; ++dim, ++params)
2476 double sigma = params.sigma_scaled(
"hessianOfGaussianMultiArray");
2477 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
2483 ParamType params_i(params_init);
2484 for (
int b=0, i=0; i<N; ++i, ++params_i)
2486 ParamType params_j(params_i);
2487 for (
int j=i; j<N; ++j, ++b, ++params_j)
2492 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 2, 1.0, opt.window_ratio);
2496 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 1, 1.0, opt.window_ratio);
2497 kernels[j].initGaussianDerivative(params_j.sigma_scaled(), 1, 1.0, opt.window_ratio);
2499 detail::scaleKernel(kernels[i], 1 / params_i.step_size());
2500 detail::scaleKernel(kernels[j], 1 / params_j.step_size());
2502 kernels.
begin(), opt.from_point, opt.to_point);
2507 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2508 class DestIterator,
class DestAccessor>
2511 DestIterator di, DestAccessor dest,
double sigma,
2517 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2518 class DestIterator,
class DestAccessor>
2521 pair<DestIterator, DestAccessor>
const & dest,
2525 dest.first, dest.second, opt );
2528 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2529 class DestIterator,
class DestAccessor>
2532 pair<DestIterator, DestAccessor>
const & dest,
2537 dest.first, dest.second, sigma, opt );
2540 template <
unsigned int N,
class T1,
class S1,
2549 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
2550 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
2551 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2552 "hessianOfGaussianMultiArray(): shape mismatch between ROI and output.");
2556 vigra_precondition(source.
shape() == dest.shape(),
2557 "hessianOfGaussianMultiArray(): shape mismatch between input and output.");
2561 destMultiArray(dest), opt );
2564 template <
unsigned int N,
class T1,
class S1,
2577 template<
int N,
class VectorType>
2578 struct StructurTensorFunctor
2580 typedef VectorType result_type;
2581 typedef typename VectorType::value_type ValueType;
2584 VectorType operator()(T
const & in)
const 2587 for(
int b=0, i=0; i<N; ++i)
2589 for(
int j=i; j<N; ++j, ++b)
2591 res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
2720 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2721 class DestIterator,
class DestAccessor>
2724 DestIterator di, DestAccessor dest,
2727 static const int N = SrcShape::static_size;
2728 static const int M = N*(N+1)/2;
2730 typedef typename DestAccessor::value_type DestType;
2731 typedef typename DestType::value_type DestValueType;
2732 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
2737 for(
int k=0; k<N; ++k)
2741 vigra_precondition(M == (
int)dest.size(di),
2742 "structureTensorMultiArray(): Wrong number of channels in output array.");
2748 SrcShape gradientShape(shape);
2749 if(opt.to_point != SrcShape())
2751 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.from_point);
2752 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.to_point);
2754 for(
int k=0; k<N; ++k, ++params)
2757 gauss.
initGaussian(params.sigma_scaled(
"structureTensorMultiArray"), 1.0, opt.window_ratio);
2758 int dilation = gauss.
right();
2759 innerOptions.from_point[k] = std::max<MultiArrayIndex>(0, opt.from_point[k] - dilation);
2760 innerOptions.to_point[k] = std::min<MultiArrayIndex>(shape[k], opt.to_point[k] + dilation);
2762 outerOptions.from_point -= innerOptions.from_point;
2763 outerOptions.to_point -= innerOptions.from_point;
2764 gradientShape = innerOptions.to_point - innerOptions.from_point;
2770 gradient.traverser_begin(), GradientAccessor(),
2772 "structureTensorMultiArray");
2775 gradientTensor.traverser_begin(), GradientTensorAccessor(),
2776 detail::StructurTensorFunctor<N, DestType>());
2779 di, dest, outerOptions,
2780 "structureTensorMultiArray");
2783 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2784 class DestIterator,
class DestAccessor>
2787 DestIterator di, DestAccessor dest,
2788 double innerScale,
double outerScale,
2792 opt.
stdDev(innerScale).outerScale(outerScale));
2795 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2796 class DestIterator,
class DestAccessor>
2799 pair<DestIterator, DestAccessor>
const & dest,
2803 dest.first, dest.second, opt );
2807 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2808 class DestIterator,
class DestAccessor>
2811 pair<DestIterator, DestAccessor>
const & dest,
2812 double innerScale,
double outerScale,
2816 dest.first, dest.second,
2817 innerScale, outerScale, opt);
2820 template <
unsigned int N,
class T1,
class S1,
2829 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
2830 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
2831 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2832 "structureTensorMultiArray(): shape mismatch between ROI and output.");
2836 vigra_precondition(source.
shape() == dest.shape(),
2837 "structureTensorMultiArray(): shape mismatch between input and output.");
2841 destMultiArray(dest), opt );
2845 template <
unsigned int N,
class T1,
class S1,
2850 double innerScale,
double outerScale,
2861 #endif //-- VIGRA_MULTI_CONVOLUTION_H 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:414
void initGaussian(double std_dev, value_type norm, double windowRatio=0.0)
Definition: separableconvolution.hxx:2253
const difference_type & shape() const
Definition: multi_array.hxx:1596
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:2422
iterator end()
Definition: tinyvector.hxx:864
ConvolutionOptions< dim > & subarray(Shape const &from, Shape const &to)
Definition: multi_convolution.hxx:495
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.
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
MultiArrayView< N+1, T, StrideTag > insertSingletonDimension(difference_type_1 i) const
Definition: multi_array.hxx:2302
void gaussianGradientMultiArray(...)
Calculate Gaussian gradient of a multi-dimensional arrays.
Definition: array_vector.hxx:58
Definition: multi_fwd.hxx:63
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:365
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:861
void symmetricGradientMultiArray(...)
Calculate gradient of a multi-dimensional arrays using symmetric difference filters.
vigra::GridGraph< N, DirectedTag >::vertex_descriptor source(typename vigra::GridGraph< N, DirectedTag >::edge_descriptor const &e, vigra::GridGraph< N, DirectedTag > const &g)
Get a vertex descriptor for the start vertex of edge e in graph g (API: boost).
Definition: multi_gridgraph.hxx:2943
ConvolutionOptions< dim > & filterWindowSize(double ratio)
Definition: multi_convolution.hxx:473
Options class template for convolutions.
Definition: multi_convolution.hxx:332
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:652
void transformMultiArray(...)
Transform a multi-dimensional array with a unary function or functor.
int right() const
Definition: separableconvolution.hxx:2157
ConvolutionOptions< dim > & innerScale(...)
MultiArrayView & init(const U &init)
Definition: multi_array.hxx:1154
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:2273
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:2132
void separableConvolveMultiArray(...)
Separated convolution on multi-dimensional arrays.