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);
862 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
863 class DestIterator,
class DestAccessor,
class KernelIterator>
866 DestIterator d, DestAccessor dest,
867 KernelIterator kernels,
868 SrcShape start = SrcShape(),
869 SrcShape stop = SrcShape())
871 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
874 if(stop != SrcShape())
877 enum { N = 1 + SrcIterator::level };
878 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, start);
879 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, stop);
881 for(
int k=0; k<N; ++k)
882 vigra_precondition(0 <= start[k] && start[k] < stop[k] && stop[k] <= shape[k],
883 "separableConvolveMultiArray(): invalid subarray shape.");
885 detail::internalSeparableConvolveSubarray(s, shape, src, d, dest, kernels, start, stop);
887 else if(!IsSameType<TmpType, typename DestAccessor::value_type>::boolResult)
891 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src,
898 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, d, dest, kernels );
902 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
903 class DestIterator,
class DestAccessor,
class T>
906 DestIterator d, DestAccessor dest,
908 SrcShape
const & start = SrcShape(),
909 SrcShape
const & stop = SrcShape())
916 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
917 class DestIterator,
class DestAccessor,
class KernelIterator>
920 pair<DestIterator, DestAccessor>
const & dest,
922 SrcShape
const & start = SrcShape(),
923 SrcShape
const & stop = SrcShape())
926 dest.first, dest.second, kit, start, stop );
929 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
930 class DestIterator,
class DestAccessor,
class T>
933 pair<DestIterator, DestAccessor>
const & dest,
935 SrcShape
const & start = SrcShape(),
936 SrcShape
const & stop = SrcShape())
941 dest.first, dest.second, kernels.begin(), start, stop);
944 template <
unsigned int N,
class T1,
class S1,
946 class KernelIterator>
956 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), start);
957 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), stop);
958 vigra_precondition(dest.
shape() == (stop - start),
959 "separableConvolveMultiArray(): shape mismatch between ROI and output.");
963 vigra_precondition(source.
shape() == dest.
shape(),
964 "separableConvolveMultiArray(): shape mismatch between input and output.");
967 destMultiArray(dest), kit, start, stop );
970 template <
unsigned int N,
class T1,
class S1,
1073 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1074 class DestIterator,
class DestAccessor,
class T>
1077 DestIterator d, DestAccessor dest,
1079 SrcShape
const & start = SrcShape(),
1080 SrcShape
const & stop = SrcShape())
1082 enum { N = 1 + SrcIterator::level };
1083 vigra_precondition( dim < N,
1084 "convolveMultiArrayOneDimension(): The dimension number to convolve must be smaller " 1085 "than the data dimensionality" );
1087 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
1094 SrcShape sstart, sstop(shape), dstart, dstop(shape);
1096 if(stop != SrcShape())
1101 sstop[dim] = shape[dim];
1102 dstop = stop - start;
1105 SNavigator snav( s, sstart, sstop, dim );
1106 DNavigator dnav( d, dstart, dstop, dim );
1108 for( ; snav.hasMore(); snav++, dnav++ )
1111 copyLine(snav.begin(), snav.end(), src,
1115 destIter( dnav.begin(), dest ),
1116 kernel1d( kernel), start[dim], stop[dim]);
1120 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1121 class DestIterator,
class DestAccessor,
class T>
1124 pair<DestIterator, DestAccessor>
const & dest,
1127 SrcShape
const & start = SrcShape(),
1128 SrcShape
const & stop = SrcShape())
1131 dest.first, dest.second, dim, kernel, start, stop);
1134 template <
unsigned int N,
class T1,
class S1,
1147 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), start);
1148 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), stop);
1149 vigra_precondition(dest.
shape() == (stop - start),
1150 "convolveMultiArrayOneDimension(): shape mismatch between ROI and output.");
1154 vigra_precondition(source.
shape() == dest.
shape(),
1155 "convolveMultiArrayOneDimension(): shape mismatch between input and output.");
1158 destMultiArray(dest), dim, kernel, start, stop);
1299 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1300 class DestIterator,
class DestAccessor>
1303 DestIterator d, DestAccessor dest,
1305 const char *
const function_name =
"gaussianSmoothMultiArray" )
1307 static const int N = SrcShape::static_size;
1312 for (
int dim = 0; dim < N; ++dim, ++params)
1313 kernels[dim].initGaussian(params.sigma_scaled(function_name,
true),
1314 1.0, opt.window_ratio);
1319 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1320 class DestIterator,
class DestAccessor>
1323 DestIterator d, DestAccessor dest,
double sigma,
1330 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1331 class DestIterator,
class DestAccessor>
1334 pair<DestIterator, DestAccessor>
const & dest,
1338 dest.first, dest.second, opt );
1341 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1342 class DestIterator,
class DestAccessor>
1345 pair<DestIterator, DestAccessor>
const & dest,
double sigma,
1349 dest.first, dest.second, sigma, opt );
1352 template <
unsigned int N,
class T1,
class S1,
1361 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
1362 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
1363 vigra_precondition(dest.
shape() == (opt.to_point - opt.from_point),
1364 "gaussianSmoothMultiArray(): shape mismatch between ROI and output.");
1368 vigra_precondition(source.
shape() == dest.
shape(),
1369 "gaussianSmoothMultiArray(): shape mismatch between input and output.");
1373 destMultiArray(dest), opt );
1376 template <
unsigned int N,
class T1,
class S1,
1511 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1512 class DestIterator,
class DestAccessor>
1515 DestIterator di, DestAccessor dest,
1517 const char *
const function_name =
"gaussianGradientMultiArray")
1519 typedef typename DestAccessor::value_type DestType;
1520 typedef typename DestType::value_type DestValueType;
1521 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1523 static const int N = SrcShape::static_size;
1526 for(
int k=0; k<N; ++k)
1530 vigra_precondition(N == (
int)dest.size(di),
1531 "gaussianGradientMultiArray(): Wrong number of channels in output array.");
1533 ParamType params = opt.scaleParams();
1534 ParamType params2(params);
1537 for (
int dim = 0; dim < N; ++dim, ++params)
1539 double sigma = params.sigma_scaled(function_name);
1540 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
1546 for (
int dim = 0; dim < N; ++dim, ++params2)
1549 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 1, 1.0, opt.window_ratio);
1550 detail::scaleKernel(kernels[dim], 1.0 / params2.step_size());
1552 opt.from_point, opt.to_point);
1556 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1557 class DestIterator,
class DestAccessor>
1560 DestIterator di, DestAccessor dest,
double sigma,
1566 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1567 class DestIterator,
class DestAccessor>
1570 pair<DestIterator, DestAccessor>
const & dest,
1574 dest.first, dest.second, opt );
1577 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1578 class DestIterator,
class DestAccessor>
1581 pair<DestIterator, DestAccessor>
const & dest,
1586 dest.first, dest.second, sigma, opt );
1589 template <
unsigned int N,
class T1,
class S1,
1598 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
1599 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
1600 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1601 "gaussianGradientMultiArray(): shape mismatch between ROI and output.");
1605 vigra_precondition(source.
shape() == dest.shape(),
1606 "gaussianGradientMultiArray(): shape mismatch between input and output.");
1610 destMultiArray(dest), opt );
1613 template <
unsigned int N,
class T1,
class S1,
1632 template <
unsigned int N,
class T1,
class S1,
1642 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.from_point);
1643 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.to_point);
1644 vigra_precondition(dest.
shape() == (opt.to_point - opt.from_point),
1645 "gaussianGradientMagnitude(): shape mismatch between ROI and output.");
1649 vigra_precondition(shape == dest.
shape(),
1650 "gaussianGradientMagnitude(): shape mismatch between input and output.");
1655 typedef typename NumericTraits<T1>::RealPromote TmpType;
1658 using namespace multi_math;
1660 for(
int k=0; k<src.
shape(N); ++k)
1672 template <
unsigned int N,
class T1,
class S1,
1679 detail::gaussianGradientMagnitudeImpl<N, T1>(src, dest, opt);
1682 template <
unsigned int N,
class T1,
class S1,
1692 template <
unsigned int N,
class T1,
int M,
class S1,
1699 detail::gaussianGradientMagnitudeImpl<N, T1>(src.
expandElements(N), dest, opt);
1702 template <
unsigned int N,
class T1,
unsigned int R,
unsigned int G,
unsigned int B,
class S1,
1709 detail::gaussianGradientMagnitudeImpl<N, T1>(src.
expandElements(N), dest, opt);
1712 template <
unsigned int N,
class T1,
class S1,
1723 template <
unsigned int N,
class T1,
class S1,
1731 gaussianGradientMagnitude<N>(src, dest, opt.
stdDev(sigma));
1844 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1845 class DestIterator,
class DestAccessor>
1848 DestIterator di, DestAccessor dest,
1851 typedef typename DestAccessor::value_type DestType;
1852 typedef typename DestType::value_type DestValueType;
1853 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1855 static const int N = SrcShape::static_size;
1858 for(
int k=0; k<N; ++k)
1862 vigra_precondition(N == (
int)dest.size(di),
1863 "symmetricGradientMultiArray(): Wrong number of channels in output array.");
1866 filter.initSymmetricDifference();
1868 StepType step_size_it = opt.stepParams();
1873 for (
int d = 0; d < N; ++d, ++step_size_it)
1876 detail::scaleKernel(symmetric, 1 / *step_size_it);
1878 di, ElementAccessor(d, dest),
1879 d, symmetric, opt.from_point, opt.to_point);
1883 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1884 class DestIterator,
class DestAccessor>
1887 pair<DestIterator, DestAccessor>
const & dest,
1891 dest.first, dest.second, opt);
1894 template <
unsigned int N,
class T1,
class S1,
1903 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
1904 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
1905 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
1906 "symmetricGradientMultiArray(): shape mismatch between ROI and output.");
1910 vigra_precondition(source.
shape() == dest.shape(),
1911 "symmetricGradientMultiArray(): shape mismatch between input and output.");
1915 destMultiArray(dest), opt);
2038 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2039 class DestIterator,
class DestAccessor>
2042 DestIterator di, DestAccessor dest,
2045 using namespace functor;
2047 typedef typename DestAccessor::value_type DestType;
2048 typedef typename NumericTraits<DestType>::RealPromote KernelType;
2051 static const int N = SrcShape::static_size;
2054 ParamType params = opt.scaleParams();
2055 ParamType params2(params);
2058 for (
int dim = 0; dim < N; ++dim, ++params)
2060 double sigma = params.sigma_scaled(
"laplacianOfGaussianMultiArray");
2061 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
2064 SrcShape dshape(shape);
2065 if(opt.to_point != SrcShape())
2066 dshape = opt.to_point - opt.from_point;
2071 for (
int dim = 0; dim < N; ++dim, ++params2)
2074 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 2, 1.0, opt.window_ratio);
2075 detail::scaleKernel(kernels[dim], 1.0 /
sq(params2.step_size()));
2080 di, dest, kernels.
begin(), opt.from_point, opt.to_point);
2085 derivative.traverser_begin(), DerivativeAccessor(),
2086 kernels.
begin(), opt.from_point, opt.to_point);
2088 di, dest, Arg1() + Arg2() );
2093 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2094 class DestIterator,
class DestAccessor>
2097 DestIterator di, DestAccessor dest,
double sigma,
2103 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2104 class DestIterator,
class DestAccessor>
2107 pair<DestIterator, DestAccessor>
const & dest,
2111 dest.first, dest.second, opt );
2114 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2115 class DestIterator,
class DestAccessor>
2118 pair<DestIterator, DestAccessor>
const & dest,
2123 dest.first, dest.second, sigma, opt );
2126 template <
unsigned int N,
class T1,
class S1,
2135 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
2136 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
2137 vigra_precondition(dest.
shape() == (opt.to_point - opt.from_point),
2138 "laplacianOfGaussianMultiArray(): shape mismatch between ROI and output.");
2142 vigra_precondition(source.
shape() == dest.
shape(),
2143 "laplacianOfGaussianMultiArray(): shape mismatch between input and output.");
2147 destMultiArray(dest), opt );
2150 template <
unsigned int N,
class T1,
class S1,
2270 template <
class Iterator,
2271 unsigned int N,
class T,
class S>
2277 typedef typename std::iterator_traits<Iterator>::value_type ArrayType;
2278 typedef typename ArrayType::value_type SrcType;
2279 typedef typename NumericTraits<SrcType>::RealPromote TmpType;
2282 vigra_precondition(std::distance(vectorField, vectorFieldEnd) == N,
2283 "gaussianDivergenceMultiArray(): wrong number of input arrays.");
2289 for(
unsigned int k = 0; k < N; ++k, ++params)
2291 sigmas[k] = params.sigma_scaled(
"gaussianDivergenceMultiArray");
2292 kernels[k].initGaussian(sigmas[k], 1.0, opt.window_ratio);
2297 for(
unsigned int k=0; k < N; ++k, ++vectorField)
2299 kernels[k].initGaussianDerivative(sigmas[k], 1, 1.0, opt.window_ratio);
2307 divergence += tmpDeriv;
2309 kernels[k].initGaussian(sigmas[k], 1.0, opt.window_ratio);
2313 template <
class Iterator,
2314 unsigned int N,
class T,
class S>
2324 template <
unsigned int N,
class T1,
class S1,
2332 for(
unsigned int k=0; k<N; ++k)
2333 field.push_back(vectorField.bindElementChannel(k));
2338 template <
unsigned int N,
class T1,
class S1,
2471 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2472 class DestIterator,
class DestAccessor>
2475 DestIterator di, DestAccessor dest,
2478 typedef typename DestAccessor::value_type DestType;
2479 typedef typename DestType::value_type DestValueType;
2480 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
2482 static const int N = SrcShape::static_size;
2483 static const int M = N*(N+1)/2;
2486 for(
int k=0; k<N; ++k)
2490 vigra_precondition(M == (
int)dest.size(di),
2491 "hessianOfGaussianMultiArray(): Wrong number of channels in output array.");
2493 ParamType params_init = opt.scaleParams();
2496 ParamType params(params_init);
2497 for (
int dim = 0; dim < N; ++dim, ++params)
2499 double sigma = params.sigma_scaled(
"hessianOfGaussianMultiArray");
2500 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio);
2506 ParamType params_i(params_init);
2507 for (
int b=0, i=0; i<N; ++i, ++params_i)
2509 ParamType params_j(params_i);
2510 for (
int j=i; j<N; ++j, ++b, ++params_j)
2515 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 2, 1.0, opt.window_ratio);
2519 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 1, 1.0, opt.window_ratio);
2520 kernels[j].initGaussianDerivative(params_j.sigma_scaled(), 1, 1.0, opt.window_ratio);
2522 detail::scaleKernel(kernels[i], 1 / params_i.step_size());
2523 detail::scaleKernel(kernels[j], 1 / params_j.step_size());
2525 kernels.
begin(), opt.from_point, opt.to_point);
2530 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2531 class DestIterator,
class DestAccessor>
2534 DestIterator di, DestAccessor dest,
double sigma,
2540 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2541 class DestIterator,
class DestAccessor>
2544 pair<DestIterator, DestAccessor>
const & dest,
2548 dest.first, dest.second, opt );
2551 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2552 class DestIterator,
class DestAccessor>
2555 pair<DestIterator, DestAccessor>
const & dest,
2560 dest.first, dest.second, sigma, opt );
2563 template <
unsigned int N,
class T1,
class S1,
2572 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
2573 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
2574 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2575 "hessianOfGaussianMultiArray(): shape mismatch between ROI and output.");
2579 vigra_precondition(source.
shape() == dest.shape(),
2580 "hessianOfGaussianMultiArray(): shape mismatch between input and output.");
2584 destMultiArray(dest), opt );
2587 template <
unsigned int N,
class T1,
class S1,
2600 template<
int N,
class VectorType>
2601 struct StructurTensorFunctor
2603 typedef VectorType result_type;
2604 typedef typename VectorType::value_type ValueType;
2607 VectorType operator()(T
const & in)
const 2610 for(
int b=0, i=0; i<N; ++i)
2612 for(
int j=i; j<N; ++j, ++b)
2614 res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
2747 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2748 class DestIterator,
class DestAccessor>
2751 DestIterator di, DestAccessor dest,
2754 static const int N = SrcShape::static_size;
2755 static const int M = N*(N+1)/2;
2757 typedef typename DestAccessor::value_type DestType;
2758 typedef typename DestType::value_type DestValueType;
2759 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
2764 for(
int k=0; k<N; ++k)
2768 vigra_precondition(M == (
int)dest.size(di),
2769 "structureTensorMultiArray(): Wrong number of channels in output array.");
2775 SrcShape gradientShape(shape);
2776 if(opt.to_point != SrcShape())
2778 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.from_point);
2779 detail::RelativeToAbsoluteCoordinate<N-1>::exec(shape, opt.to_point);
2781 for(
int k=0; k<N; ++k, ++params)
2784 gauss.
initGaussian(params.sigma_scaled(
"structureTensorMultiArray"), 1.0, opt.window_ratio);
2785 int dilation = gauss.
right();
2786 innerOptions.from_point[k] = std::max<MultiArrayIndex>(0, opt.from_point[k] - dilation);
2787 innerOptions.to_point[k] = std::min<MultiArrayIndex>(shape[k], opt.to_point[k] + dilation);
2789 outerOptions.from_point -= innerOptions.from_point;
2790 outerOptions.to_point -= innerOptions.from_point;
2791 gradientShape = innerOptions.to_point - innerOptions.from_point;
2797 gradient.traverser_begin(), GradientAccessor(),
2799 "structureTensorMultiArray");
2802 gradientTensor.traverser_begin(), GradientTensorAccessor(),
2803 detail::StructurTensorFunctor<N, DestType>());
2806 di, dest, outerOptions,
2807 "structureTensorMultiArray");
2810 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2811 class DestIterator,
class DestAccessor>
2814 DestIterator di, DestAccessor dest,
2815 double innerScale,
double outerScale,
2819 opt.
stdDev(innerScale).outerScale(outerScale));
2822 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2823 class DestIterator,
class DestAccessor>
2826 pair<DestIterator, DestAccessor>
const & dest,
2830 dest.first, dest.second, opt );
2834 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
2835 class DestIterator,
class DestAccessor>
2838 pair<DestIterator, DestAccessor>
const & dest,
2839 double innerScale,
double outerScale,
2843 dest.first, dest.second,
2844 innerScale, outerScale, opt);
2847 template <
unsigned int N,
class T1,
class S1,
2856 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.from_point);
2857 detail::RelativeToAbsoluteCoordinate<N-1>::exec(source.
shape(), opt.to_point);
2858 vigra_precondition(dest.shape() == (opt.to_point - opt.from_point),
2859 "structureTensorMultiArray(): shape mismatch between ROI and output.");
2863 vigra_precondition(source.
shape() == dest.shape(),
2864 "structureTensorMultiArray(): shape mismatch between input and output.");
2868 destMultiArray(dest), opt );
2872 template <
unsigned int N,
class T1,
class S1,
2877 double innerScale,
double outerScale,
2888 #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 symmetricGradientMultiArray(...)
Calculate gradient of a multi-dimensional arrays using symmetric difference filters.
int right() const
Definition: separableconvolution.hxx:2157
void initGaussian(double std_dev, value_type norm, double windowRatio=0.0)
Definition: separableconvolution.hxx:2253
MultiArrayView< N-M, T, StrideTag > bindOuter(const TinyVector< Index, M > &d) const
Definition: multi_array.hxx:2184
void gaussianDivergenceMultiArray(...)
Calculate the divergence of a vector field using Gaussian derivative filters.
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2474
void separableConvolveMultiArray(...)
Separated convolution on multi-dimensional arrays.
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. ...
MultiArrayView< N+1, T, StrideTag > insertSingletonDimension(difference_type_1 i) const
Definition: multi_array.hxx:2354
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
void gaussianGradientMultiArray(...)
Calculate Gaussian gradient of a multi-dimensional arrays.
Definition: array_vector.hxx:58
Definition: multi_fwd.hxx:63
void laplacianOfGaussianMultiArray(...)
Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters.
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:382
A navigator that provides access to the 1D subranges of an n-dimensional range given by a vigra::Mult...
Definition: navigator.hxx:97
const difference_type & shape() const
Definition: multi_array.hxx:1648
iterator begin()
Definition: tinyvector.hxx:861
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
void gaussianGradientMagnitude(...)
Calculate the gradient magnitude by means of a 1st derivatives of Gaussian filter.
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:1459
void combineTwoMultiArrays(...)
Combine two multi-dimensional arrays into one using a binary function or functor. ...
const_iterator begin() const
Definition: array_vector.hxx:223
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:704
void transformMultiArray(...)
Transform a multi-dimensional array with a unary function or functor.
ConvolutionOptions< dim > & innerScale(...)
MultiArrayView & init(const U &init)
Definition: multi_array.hxx:1206
void convolveMultiArrayOneDimension(...)
Convolution along a single dimension of a multi-dimensional arrays.
void gaussianSmoothMultiArray(...)
Isotropic Gaussian smoothing of a multi-dimensional arrays.
Class for a single RGB value.
Definition: accessor.hxx:938
Encapsulate access to the values an iterator points to.
Definition: accessor.hxx:133
void hessianOfGaussianMultiArray(...)
Calculate Hessian matrix of a N-dimensional arrays using Gaussian derivative filters.
const_iterator end() const
Definition: array_vector.hxx:237
MultiArrayView< N+1, typename ExpandElementResult< T >::type, StridedArrayTag > expandElements(difference_type_1 d) const
Definition: multi_array.hxx:2325
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
ConvolutionOptions< dim > & stdDev(...)
void structureTensorMultiArray(...)
Calculate the structure tensor of a multi-dimensional arrays.