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 ",
138 bool allow_zero =
false)
const 140 sigma_precondition(sigma_eff(), function_name);
141 sigma_precondition(sigma_d(), function_name);
142 double sigma_squared =
sq(sigma_eff()) -
sq(sigma_d());
143 if (sigma_squared > 0.0 || (allow_zero && sigma_squared == 0.0))
145 return std::sqrt(sigma_squared) / step_size();
149 std::string msg =
"(): Scale would be imaginary";
152 vigra_precondition(
false, function_name + msg +
".");
158 template <
unsigned dim>
159 struct multiArrayScaleParam
161 typedef TinyVector<double, dim> p_vector;
162 typedef typename p_vector::const_iterator return_type;
165 template <
class Param>
166 multiArrayScaleParam(Param val,
const char *
const function_name =
"multiArrayScaleParam")
168 typename WrapDoubleIterator<Param>::type in(val, dim, function_name);
169 for (
unsigned i = 0; i != dim; ++i, ++in)
172 return_type operator()()
const 176 static void precondition(
unsigned n_par,
const char *
const function_name =
"multiArrayScaleParam")
180 std::string msg =
"(): dimension parameter must be ";
181 vigra_precondition(dim == n_par, function_name + msg + n);
183 multiArrayScaleParam(
double v0,
double v1,
const char *
const function_name =
"multiArrayScaleParam")
185 precondition(2, function_name);
186 vec = p_vector(v0, v1);
188 multiArrayScaleParam(
double v0,
double v1,
double v2,
const char *
const function_name =
"multiArrayScaleParam")
190 precondition(3, function_name);
191 vec = p_vector(v0, v1, v2);
193 multiArrayScaleParam(
double v0,
double v1,
double v2,
double v3,
const char *
const function_name =
"multiArrayScaleParam")
195 precondition(4, function_name);
196 vec = p_vector(v0, v1, v2, v3);
198 multiArrayScaleParam(
double v0,
double v1,
double v2,
double v3,
double v4,
const char *
const function_name =
"multiArrayScaleParam")
200 precondition(5, function_name);
201 vec = p_vector(v0, v1, v2, v3, v4);
207 #define VIGRA_CONVOLUTION_OPTIONS(function_name, default_value, member_name, getter_setter_name) \ 208 template <class Param> \ 209 ConvolutionOptions & function_name(const Param & val) \ 211 member_name = ParamVec(val, "ConvolutionOptions::" #function_name); \ 214 ConvolutionOptions & function_name() \ 216 member_name = ParamVec(default_value, "ConvolutionOptions::" #function_name); \ 219 ConvolutionOptions & function_name(double v0, double v1) \ 221 member_name = ParamVec(v0, v1, "ConvolutionOptions::" #function_name); \ 224 ConvolutionOptions & function_name(double v0, double v1, double v2) \ 226 member_name = ParamVec(v0, v1, v2, "ConvolutionOptions::" #function_name); \ 229 ConvolutionOptions & function_name(double v0, double v1, double v2, double v3) \ 231 member_name = ParamVec(v0, v1, v2, v3, "ConvolutionOptions::" #function_name); \ 234 ConvolutionOptions & function_name(double v0, double v1, double v2, double v3, double v4) \ 236 member_name = ParamVec(v0, v1, v2, v3, v4, "ConvolutionOptions::" #function_name); \ 239 typename ParamVec::p_vector get##getter_setter_name()const{ \ 240 return member_name.vec; \ 242 void set##getter_setter_name(typename ParamVec::p_vector vec){ \ 243 member_name.vec = vec; \ 334 template <
unsigned dim>
339 typedef detail::multiArrayScaleParam<dim> ParamVec;
340 typedef typename ParamVec::return_type ParamIt;
345 ParamVec outer_scale;
347 Shape from_point, to_point;
357 typedef typename detail::WrapDoubleIteratorTriple<ParamIt, ParamIt, ParamIt>
359 typedef typename detail::WrapDoubleIterator<ParamIt>::type
362 ScaleIterator scaleParams()
const 364 return ScaleIterator(sigma_eff(), sigma_d(), step_size());
366 StepIterator stepParams()
const 368 return StepIterator(step_size());
375 return outer.
stdDev(outer_scale()).resolutionStdDev(0.0);
380 VIGRA_CONVOLUTION_OPTIONS(stepSize, 1.0, step_size, StepSize)
399 VIGRA_CONVOLUTION_OPTIONS(resolutionStdDev, 0.0, sigma_d, ResolutionStdDev)
415 VIGRA_CONVOLUTION_OPTIONS(stdDev, 0.0, sigma_eff, StdDev)
416 VIGRA_CONVOLUTION_OPTIONS(innerScale, 0.0, sigma_eff, InnerScale)
446 VIGRA_CONVOLUTION_OPTIONS(outerScale, 0.0, outer_scale, OuterScale)
478 vigra_precondition(ratio >= 0.0,
479 "ConvolutionOptions::filterWindowSize(): ratio must not be negative.");
480 window_ratio = ratio;
484 double getFilterWindowSize()
const {
505 std::pair<Shape, Shape> getSubarray()
const{
506 std::pair<Shape, Shape> res;
507 res.first = from_point;
508 res.second = to_point;
522 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
523 class DestIterator,
class DestAccessor,
class KernelIterator>
525 internalSeparableConvolveMultiArrayTmp(
526 SrcIterator si, SrcShape
const & shape, SrcAccessor src,
527 DestIterator di, DestAccessor dest, KernelIterator kit)
529 enum { N = 1 + SrcIterator::level };
531 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
544 SNavigator snav( si, shape, 0 );
545 DNavigator dnav( di, shape, 0 );
547 for( ; snav.hasMore(); snav++, dnav++ )
550 copyLine(snav.begin(), snav.end(), src, tmp.
begin(), acc);
553 destIter( dnav.begin(), dest ),
560 for(
int d = 1; d < N; ++d, ++kit )
562 DNavigator dnav( di, shape, d );
564 tmp.resize( shape[d] );
566 for( ; dnav.hasMore(); dnav++ )
569 copyLine(dnav.begin(), dnav.end(), dest, tmp.
begin(), acc);
572 destIter( dnav.begin(), dest ),
584 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
585 class DestIterator,
class DestAccessor,
class KernelIterator>
587 internalSeparableConvolveSubarray(
588 SrcIterator si, SrcShape
const & shape, SrcAccessor src,
589 DestIterator di, DestAccessor dest, KernelIterator kit,
590 SrcShape
const & start, SrcShape
const & stop)
592 enum { N = 1 + SrcIterator::level };
594 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
596 typedef typename TmpArray::traverser TmpIterator;
599 SrcShape sstart, sstop, axisorder, tmpshape;
601 for(
int k=0; k<N; ++k)
604 sstart[k] = start[k] - kit[k].right();
607 sstop[k] = stop[k] - kit[k].left();
608 if(sstop[k] > shape[k])
610 overhead[k] = double(sstop[k] - sstart[k]) / (stop[k] - start[k]);
613 indexSort(overhead.
begin(), overhead.
end(), axisorder.begin(), std::greater<double>());
614 SrcShape dstart, dstop(sstop - sstart);
615 dstop[axisorder[0]] = stop[axisorder[0]] - start[axisorder[0]];
627 SNavigator snav( si, sstart, sstop, axisorder[0]);
628 TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[0]);
632 int lstart = start[axisorder[0]] - sstart[axisorder[0]];
633 int lstop = lstart + (stop[axisorder[0]] - start[axisorder[0]]);
635 for( ; snav.hasMore(); snav++, tnav++ )
638 copyLine(snav.begin(), snav.end(), src, tmpline.begin(), acc);
640 convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
641 destIter(tnav.begin(), acc),
642 kernel1d( kit[axisorder[0]] ), lstart, lstop);
647 for(
int d = 1; d < N; ++d)
649 TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[d]);
653 int lstart = start[axisorder[d]] - sstart[axisorder[d]];
654 int lstop = lstart + (stop[axisorder[d]] - start[axisorder[d]]);
656 for( ; tnav.hasMore(); tnav++ )
659 copyLine(tnav.begin(), tnav.end(), acc, tmpline.begin(), acc );
661 convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc),
662 destIter( tnav.begin() + lstart, acc ),
663 kernel1d( kit[axisorder[d]] ), lstart, lstop);
666 dstart[axisorder[d]] = lstart;
667 dstop[axisorder[d]] = lstop;
670 copyMultiArray(tmp.traverser_begin()+dstart, stop-start, acc, di, dest);
676 scaleKernel(K & kernel,
double a)
678 for(
int i = kernel.left(); i <= kernel.right(); ++i)
679 kernel[i] = detail::RequiresExplicitCast<typename K::value_type>::cast(kernel[i] * a);
867 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
868 class DestIterator,
class DestAccessor,
class KernelIterator>
871 DestIterator d, DestAccessor dest,
872 KernelIterator kernels,
873 SrcShape start = SrcShape(),
874 SrcShape stop = SrcShape())
876 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
879 if(stop != SrcShape())
882 enum { N = 1 + SrcIterator::level };
883 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, start);
884 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, stop);
886 for(
int k=0; k<N; ++k)
887 vigra_precondition(0 <= start[k] && start[k] < stop[k] && stop[k] <= shape[k],
888 "separableConvolveMultiArray(): invalid subarray shape.");
890 detail::internalSeparableConvolveSubarray(s, shape, src, d, dest, kernels, start, stop);
892 else if(!IsSameType<TmpType, typename DestAccessor::value_type>::boolResult)
896 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src,
903 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, d, dest, kernels );
907 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
908 class DestIterator,
class DestAccessor,
class T>
911 DestIterator d, DestAccessor dest,
913 SrcShape
const & start = SrcShape(),
914 SrcShape
const & stop = SrcShape())
921 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
922 class DestIterator,
class DestAccessor,
class KernelIterator>
925 pair<DestIterator, DestAccessor>
const & dest,
927 SrcShape
const & start = SrcShape(),
928 SrcShape
const & stop = SrcShape())
931 dest.first, dest.second, kit, start, stop );
934 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
935 class DestIterator,
class DestAccessor,
class T>
938 pair<DestIterator, DestAccessor>
const & dest,
940 SrcShape
const & start = SrcShape(),
941 SrcShape
const & stop = SrcShape())
946 dest.first, dest.second, kernels.begin(), start, stop);
949 template <
unsigned int N,
class T1,
class S1,
951 class KernelIterator>
961 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), start);
962 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), stop);
963 vigra_precondition(dest.
shape() == (stop - start),
964 "separableConvolveMultiArray(): shape mismatch between ROI and output.");
968 vigra_precondition(source.
shape() == dest.
shape(),
969 "separableConvolveMultiArray(): shape mismatch between input and output.");
972 destMultiArray(dest), kit, start, stop );
975 template <
unsigned int N,
class T1,
class S1,
1078 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1079 class DestIterator,
class DestAccessor,
class T>
1082 DestIterator d, DestAccessor dest,
1084 SrcShape
const & start = SrcShape(),
1085 SrcShape
const & stop = SrcShape())
1087 enum { N = 1 + SrcIterator::level };
1088 vigra_precondition( dim < N,
1089 "convolveMultiArrayOneDimension(): The dimension number to convolve must be smaller " 1090 "than the data dimensionality" );
1092 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
1099 SrcShape sstart, sstop(shape), dstart, dstop(shape);
1101 if(stop != SrcShape())
1106 sstop[dim] = shape[dim];
1107 dstop = stop - start;
1110 SNavigator snav( s, sstart, sstop, dim );
1111 DNavigator dnav( d, dstart, dstop, dim );
1113 for( ; snav.hasMore(); snav++, dnav++ )
1116 copyLine(snav.begin(), snav.end(), src,
1120 destIter( dnav.begin(), dest ),
1121 kernel1d( kernel), start[dim], stop[dim]);
1125 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1126 class DestIterator,
class DestAccessor,
class T>
1129 pair<DestIterator, DestAccessor>
const & dest,
1132 SrcShape
const & start = SrcShape(),
1133 SrcShape
const & stop = SrcShape())
1136 dest.first, dest.second, dim, kernel, start, stop);
1139 template <
unsigned int N,
class T1,
class S1,
1152 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), start);
1153 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), stop);
1154 vigra_precondition(dest.
shape() == (stop - start),
1155 "convolveMultiArrayOneDimension(): shape mismatch between ROI and output.");
1159 vigra_precondition(source.
shape() == dest.
shape(),
1160 "convolveMultiArrayOneDimension(): shape mismatch between input and output.");
1163 destMultiArray(dest), dim, kernel, start, stop);
1300 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1301 class DestIterator,
class DestAccessor>
1304 DestIterator d, DestAccessor dest,
1306 const char *
const function_name =
"gaussianSmoothMultiArray" )
1308 static const int N = SrcShape::static_size;
1313 for (
int dim = 0; dim < N; ++dim, ++params)
1314 kernels[dim].initGaussian(params.sigma_scaled(function_name,
true),
1315 1.0, opt.window_ratio);
1320 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1321 class DestIterator,
class DestAccessor>
1324 DestIterator d, DestAccessor dest,
double sigma,
1331 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1332 class DestIterator,
class DestAccessor>
1335 pair<DestIterator, DestAccessor>
const & dest,
1339 dest.first, dest.second, opt );
1342 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1343 class DestIterator,
class DestAccessor>
1346 pair<DestIterator, DestAccessor>
const & dest,
double sigma,
1350 dest.first, dest.second, sigma, opt );
1353 template <
unsigned int N,
class T1,
class S1,
1362 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
1363 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
1364 vigra_precondition(dest.
shape() == (opt.to_point - opt.from_point),
1365 "gaussianSmoothMultiArray(): shape mismatch between ROI and output.");
1369 vigra_precondition(source.
shape() == dest.
shape(),
1370 "gaussianSmoothMultiArray(): shape mismatch between input and output.");
1374 destMultiArray(dest), opt );
1377 template <
unsigned int N,
class T1,
class S1,
1508 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1509 class DestIterator,
class DestAccessor>
1512 DestIterator di, DestAccessor dest,
1514 const char *
const function_name =
"gaussianGradientMultiArray")
1516 typedef typename DestAccessor::value_type DestType;
1517 typedef typename DestType::value_type DestValueType;
1518 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1520 static const int N = SrcShape::static_size;
1523 for(
int k=0; k<N; ++k)
1527 vigra_precondition(N == (
int)dest.size(di),
1528 "gaussianGradientMultiArray(): Wrong number of channels in output array.");
1530 ParamType params = opt.scaleParams();
1531 ParamType params2(params);
1534 for (
int dim = 0; dim < N; ++dim, ++params)
1536 double sigma = params.sigma_scaled(function_name);
1537 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
1543 for (
int dim = 0; dim < N; ++dim, ++params2)
1546 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 1, 1.0, opt.window_ratio);
1547 detail::scaleKernel(kernels[dim], 1.0 / params2.step_size());
1549 opt.from_point, opt.to_point);
1553 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1554 class DestIterator,
class DestAccessor>
1557 DestIterator di, DestAccessor dest,
double sigma,
1563 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1564 class DestIterator,
class DestAccessor>
1567 pair<DestIterator, DestAccessor>
const & dest,
1571 dest.first, dest.second, opt );
1574 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1575 class DestIterator,
class DestAccessor>
1578 pair<DestIterator, DestAccessor>
const & dest,
1583 dest.first, dest.second, sigma, opt );
1586 template <
unsigned int N,
class T1,
class S1,
1595 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
1596 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
1597 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1598 "gaussianGradientMultiArray(): shape mismatch between ROI and output.");
1602 vigra_precondition(source.
shape() == dest.shape(),
1603 "gaussianGradientMultiArray(): shape mismatch between input and output.");
1607 destMultiArray(dest), opt );
1610 template <
unsigned int N,
class T1,
class S1,
1629 template <
unsigned int N,
class T1,
class S1,
1639 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.from_point);
1640 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.to_point);
1641 vigra_precondition(dest.
shape() == (opt.to_point - opt.from_point),
1642 "gaussianGradientMagnitude(): shape mismatch between ROI and output.");
1646 vigra_precondition(shape == dest.
shape(),
1647 "gaussianGradientMagnitude(): shape mismatch between input and output.");
1652 typedef typename NumericTraits<T1>::RealPromote TmpType;
1655 using namespace multi_math;
1657 for(
int k=0; k<src.
shape(N); ++k)
1669 template <
unsigned int N,
class T1,
class S1,
1676 detail::gaussianGradientMagnitudeImpl<N, T1>(src, dest, opt);
1679 template <
unsigned int N,
class T1,
class S1,
1689 template <
unsigned int N,
class T1,
int M,
class S1,
1696 detail::gaussianGradientMagnitudeImpl<N, T1>(src.
expandElements(N), dest, opt);
1699 template <
unsigned int N,
class T1,
unsigned int R,
unsigned int G,
unsigned int B,
class S1,
1706 detail::gaussianGradientMagnitudeImpl<N, T1>(src.
expandElements(N), dest, opt);
1709 template <
unsigned int N,
class T1,
class S1,
1720 template <
unsigned int N,
class T1,
class S1,
1728 gaussianGradientMagnitude<N>(src, dest, opt.
stdDev(sigma));
1837 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1838 class DestIterator,
class DestAccessor>
1841 DestIterator di, DestAccessor dest,
1844 typedef typename DestAccessor::value_type DestType;
1845 typedef typename DestType::value_type DestValueType;
1846 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1848 static const int N = SrcShape::static_size;
1851 for(
int k=0; k<N; ++k)
1855 vigra_precondition(N == (
int)dest.size(di),
1856 "symmetricGradientMultiArray(): Wrong number of channels in output array.");
1859 filter.initSymmetricDifference();
1861 StepType step_size_it = opt.stepParams();
1866 for (
int d = 0; d < N; ++d, ++step_size_it)
1869 detail::scaleKernel(symmetric, 1 / *step_size_it);
1871 di, ElementAccessor(d, dest),
1872 d, symmetric, opt.from_point, opt.to_point);
1876 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1877 class DestIterator,
class DestAccessor>
1880 pair<DestIterator, DestAccessor>
const & dest,
1884 dest.first, dest.second, opt);
1887 template <
unsigned int N,
class T1,
class S1,
1896 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
1897 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
1898 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1899 "symmetricGradientMultiArray(): shape mismatch between ROI and output.");
1903 vigra_precondition(source.
shape() == dest.shape(),
1904 "symmetricGradientMultiArray(): shape mismatch between input and output.");
1908 destMultiArray(dest), opt);
2027 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2028 class DestIterator,
class DestAccessor>
2031 DestIterator di, DestAccessor dest,
2034 using namespace functor;
2036 typedef typename DestAccessor::value_type DestType;
2037 typedef typename NumericTraits<DestType>::RealPromote KernelType;
2040 static const int N = SrcShape::static_size;
2043 ParamType params = opt.scaleParams();
2044 ParamType params2(params);
2047 for (
int dim = 0; dim < N; ++dim, ++params)
2049 double sigma = params.sigma_scaled(
"laplacianOfGaussianMultiArray");
2050 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
2053 SrcShape dshape(shape);
2054 if(opt.to_point != SrcShape())
2055 dshape = opt.to_point - opt.from_point;
2060 for (
int dim = 0; dim < N; ++dim, ++params2)
2063 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 2, 1.0, opt.window_ratio);
2064 detail::scaleKernel(kernels[dim], 1.0 /
sq(params2.step_size()));
2069 di, dest, kernels.
begin(), opt.from_point, opt.to_point);
2074 derivative.traverser_begin(), DerivativeAccessor(),
2075 kernels.
begin(), opt.from_point, opt.to_point);
2077 di, dest, Arg1() + Arg2() );
2082 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2083 class DestIterator,
class DestAccessor>
2086 DestIterator di, DestAccessor dest,
double sigma,
2092 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2093 class DestIterator,
class DestAccessor>
2096 pair<DestIterator, DestAccessor>
const & dest,
2100 dest.first, dest.second, opt );
2103 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2104 class DestIterator,
class DestAccessor>
2107 pair<DestIterator, DestAccessor>
const & dest,
2112 dest.first, dest.second, sigma, opt );
2115 template <
unsigned int N,
class T1,
class S1,
2124 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
2125 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
2126 vigra_precondition(dest.
shape() == (opt.to_point - opt.from_point),
2127 "laplacianOfGaussianMultiArray(): shape mismatch between ROI and output.");
2131 vigra_precondition(source.
shape() == dest.
shape(),
2132 "laplacianOfGaussianMultiArray(): shape mismatch between input and output.");
2136 destMultiArray(dest), opt );
2139 template <
unsigned int N,
class T1,
class S1,
2255 template <
class Iterator,
2256 unsigned int N,
class T,
class S>
2262 typedef typename std::iterator_traits<Iterator>::value_type ArrayType;
2263 typedef typename ArrayType::value_type SrcType;
2264 typedef typename NumericTraits<SrcType>::RealPromote TmpType;
2267 vigra_precondition(std::distance(vectorField, vectorFieldEnd) == N,
2268 "gaussianDivergenceMultiArray(): wrong number of input arrays.");
2274 for(
unsigned int k = 0; k < N; ++k, ++params)
2276 sigmas[k] = params.sigma_scaled(
"gaussianDivergenceMultiArray");
2277 kernels[k].initGaussian(sigmas[k], 1.0, opt.window_ratio);
2282 for(
unsigned int k=0; k < N; ++k, ++vectorField)
2284 kernels[k].initGaussianDerivative(sigmas[k], 1, 1.0, opt.window_ratio);
2292 divergence += tmpDeriv;
2294 kernels[k].initGaussian(sigmas[k], 1.0, opt.window_ratio);
2298 template <
class Iterator,
2299 unsigned int N,
class T,
class S>
2309 template <
unsigned int N,
class T1,
class S1,
2317 for(
unsigned int k=0; k<N; ++k)
2318 field.push_back(vectorField.bindElementChannel(k));
2323 template <
unsigned int N,
class T1,
class S1,
2452 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2453 class DestIterator,
class DestAccessor>
2456 DestIterator di, DestAccessor dest,
2459 typedef typename DestAccessor::value_type DestType;
2460 typedef typename DestType::value_type DestValueType;
2461 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
2463 static const int N = SrcShape::static_size;
2464 static const int M = N*(N+1)/2;
2467 for(
int k=0; k<N; ++k)
2471 vigra_precondition(M == (
int)dest.size(di),
2472 "hessianOfGaussianMultiArray(): Wrong number of channels in output array.");
2474 ParamType params_init = opt.scaleParams();
2477 ParamType params(params_init);
2478 for (
int dim = 0; dim < N; ++dim, ++params)
2480 double sigma = params.sigma_scaled(
"hessianOfGaussianMultiArray");
2481 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
2487 ParamType params_i(params_init);
2488 for (
int b=0, i=0; i<N; ++i, ++params_i)
2490 ParamType params_j(params_i);
2491 for (
int j=i; j<N; ++j, ++b, ++params_j)
2496 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 2, 1.0, opt.window_ratio);
2500 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 1, 1.0, opt.window_ratio);
2501 kernels[j].initGaussianDerivative(params_j.sigma_scaled(), 1, 1.0, opt.window_ratio);
2503 detail::scaleKernel(kernels[i], 1 / params_i.step_size());
2504 detail::scaleKernel(kernels[j], 1 / params_j.step_size());
2506 kernels.
begin(), opt.from_point, opt.to_point);
2511 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2512 class DestIterator,
class DestAccessor>
2515 DestIterator di, DestAccessor dest,
double sigma,
2521 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2522 class DestIterator,
class DestAccessor>
2525 pair<DestIterator, DestAccessor>
const & dest,
2529 dest.first, dest.second, opt );
2532 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2533 class DestIterator,
class DestAccessor>
2536 pair<DestIterator, DestAccessor>
const & dest,
2541 dest.first, dest.second, sigma, opt );
2544 template <
unsigned int N,
class T1,
class S1,
2553 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
2554 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
2555 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2556 "hessianOfGaussianMultiArray(): shape mismatch between ROI and output.");
2560 vigra_precondition(source.
shape() == dest.shape(),
2561 "hessianOfGaussianMultiArray(): shape mismatch between input and output.");
2565 destMultiArray(dest), opt );
2568 template <
unsigned int N,
class T1,
class S1,
2581 template<
int N,
class VectorType>
2582 struct StructurTensorFunctor
2584 typedef VectorType result_type;
2585 typedef typename VectorType::value_type ValueType;
2588 VectorType operator()(T
const & in)
const 2591 for(
int b=0, i=0; i<N; ++i)
2593 for(
int j=i; j<N; ++j, ++b)
2595 res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
2724 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2725 class DestIterator,
class DestAccessor>
2728 DestIterator di, DestAccessor dest,
2731 static const int N = SrcShape::static_size;
2732 static const int M = N*(N+1)/2;
2734 typedef typename DestAccessor::value_type DestType;
2735 typedef typename DestType::value_type DestValueType;
2736 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
2741 for(
int k=0; k<N; ++k)
2745 vigra_precondition(M == (
int)dest.size(di),
2746 "structureTensorMultiArray(): Wrong number of channels in output array.");
2752 SrcShape gradientShape(shape);
2753 if(opt.to_point != SrcShape())
2755 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.from_point);
2756 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.to_point);
2758 for(
int k=0; k<N; ++k, ++params)
2761 gauss.
initGaussian(params.sigma_scaled(
"structureTensorMultiArray"), 1.0, opt.window_ratio);
2762 int dilation = gauss.
right();
2763 innerOptions.from_point[k] = std::max<MultiArrayIndex>(0, opt.from_point[k] - dilation);
2764 innerOptions.to_point[k] = std::min<MultiArrayIndex>(shape[k], opt.to_point[k] + dilation);
2766 outerOptions.from_point -= innerOptions.from_point;
2767 outerOptions.to_point -= innerOptions.from_point;
2768 gradientShape = innerOptions.to_point - innerOptions.from_point;
2774 gradient.traverser_begin(), GradientAccessor(),
2776 "structureTensorMultiArray");
2779 gradientTensor.traverser_begin(), GradientTensorAccessor(),
2780 detail::StructurTensorFunctor<N, DestType>());
2783 di, dest, outerOptions,
2784 "structureTensorMultiArray");
2787 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2788 class DestIterator,
class DestAccessor>
2791 DestIterator di, DestAccessor dest,
2792 double innerScale,
double outerScale,
2796 opt.
stdDev(innerScale).outerScale(outerScale));
2799 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2800 class DestIterator,
class DestAccessor>
2803 pair<DestIterator, DestAccessor>
const & dest,
2807 dest.first, dest.second, opt );
2811 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2812 class DestIterator,
class DestAccessor>
2815 pair<DestIterator, DestAccessor>
const & dest,
2816 double innerScale,
double outerScale,
2820 dest.first, dest.second,
2821 innerScale, outerScale, opt);
2824 template <
unsigned int N,
class T1,
class S1,
2833 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
2834 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
2835 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2836 "structureTensorMultiArray(): shape mismatch between ROI and output.");
2840 vigra_precondition(source.
shape() == dest.shape(),
2841 "structureTensorMultiArray(): shape mismatch between input and output.");
2845 destMultiArray(dest), opt );
2849 template <
unsigned int N,
class T1,
class S1,
2854 double innerScale,
double outerScale,
2865 #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:498
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:476
Options class template for convolutions.
Definition: multi_convolution.hxx:335
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.