37#ifndef VIGRA_RECURSIVECONVOLUTION_HXX
38#define VIGRA_RECURSIVECONVOLUTION_HXX
43#include "numerictraits.hxx"
44#include "imageiteratoradapter.hxx"
45#include "bordertreatment.hxx"
46#include "array_vector.hxx"
47#include "multi_shape.hxx"
174 vigra_precondition(-1.0 < b && b < 1.0,
175 "recursiveFilterLine(): -1 < factor < 1 required.\n");
187 double eps = 0.00001;
188 int kernelw = std::min(w-1, (
int)(VIGRA_CSTD::log(
eps)/VIGRA_CSTD::log(VIGRA_CSTD::fabs(b))));
191 NumericTraits<typename SrcAccessor::value_type>::RealPromote
TempType;
193 typedef typename DestTraits::RealPromote RealPromote;
196 std::vector<TempType>
vline(w);
199 double norm = (1.0 - b) / (1.0 + b);
203 if(
border == BORDER_TREATMENT_REPEAT ||
204 border == BORDER_TREATMENT_AVOID)
208 else if(
border == BORDER_TREATMENT_REFLECT)
215 else if(
border == BORDER_TREATMENT_WRAP)
222 else if(
border == BORDER_TREATMENT_CLIP ||
223 border == BORDER_TREATMENT_ZEROPAD)
225 old = NumericTraits<TempType>::zero();
229 vigra_fail(
"recursiveFilterLine(): Unknown border treatment mode.\n");
230 old = NumericTraits<TempType>::zero();
241 if(
border == BORDER_TREATMENT_REPEAT ||
242 border == BORDER_TREATMENT_AVOID)
247 else if(
border == BORDER_TREATMENT_REFLECT)
251 else if(
border == BORDER_TREATMENT_WRAP)
258 else if(
border == BORDER_TREATMENT_CLIP ||
259 border == BORDER_TREATMENT_ZEROPAD)
261 old = NumericTraits<TempType>::zero();
266 if(
border == BORDER_TREATMENT_CLIP)
270 double bleft = VIGRA_CSTD::pow(b, w);
272 for(x=w-1; x>=0; --x, --
is, --id)
282 else if(
border == BORDER_TREATMENT_AVOID)
284 for(x=w-1; x >=
kernelw; --x, --
is, --id)
289 ad.set(DestTraits::fromRealPromote(RealPromote(
norm * (
line[x] + f))),
id);
294 for(x=w-1; x>=0; --x, --
is, --id)
298 ad.set(DestTraits::fromRealPromote(RealPromote(
norm * (
line[x] + f))),
id);
309template <
class SrcIterator,
class SrcAccessor,
310 class DestIterator,
class DestAccessor>
312 DestIterator
id, DestAccessor ad,
double b1,
double b2)
318 NumericTraits<typename SrcAccessor::value_type>::RealPromote TempType;
321 std::vector<TempType> vline(w+1);
322 typename std::vector<TempType>::iterator line = vline.begin();
324 double norm = 1.0 - b1 - b2;
325 double norm1 = (1.0 - b1 - b2) / (1.0 + b1 + b2);
330 int kernelw = std::min(w-1, std::max(8, (
int)(1.0 /
norm + 0.5)));
332 line[kernelw] = as(is);
333 line[kernelw-1] = as(is);
334 for(x = kernelw - 2; x > 0; --x, --is)
336 line[x] = detail::RequiresExplicitCast<TempType>::cast(as(is) + b1 * line[x+1] + b2 * line[x+2]);
338 line[0] = detail::RequiresExplicitCast<TempType>::cast(as(is) + b1 * line[1] + b2 * line[2]);
340 line[1] = detail::RequiresExplicitCast<TempType>::cast(as(is) + b1 * line[0] + b2 * line[1]);
342 for(x=2; x < w; ++x, ++is)
344 line[x] = detail::RequiresExplicitCast<TempType>::cast(as(is) + b1 * line[x-1] + b2 * line[x-2]);
348 line[w-1] = detail::RequiresExplicitCast<TempType>::cast(norm1 * (line[w-1] + b1 * line[w-2] + b2 * line[w-3]));
349 line[w-2] = detail::RequiresExplicitCast<TempType>::cast(norm1 * (line[w-2] + b1 * line[w] + b2 * line[w-2]));
351 ad.set(line[w-1],
id);
353 ad.set(line[w-2],
id);
355 for(x=w-3; x>=0; --x, --id, --is)
357 line[x] = detail::RequiresExplicitCast<TempType>::cast(norm2 * line[x] + b1 * line[x+1] + b2 * line[x+2]);
452 double q = 1.31564 * (std::sqrt(1.0 + 0.490811 * sigma*sigma) - 1.0);
455 double b0 = 1.0/(1.57825 + 2.44413*
q + 1.4281*
qq + 0.422205*
qqq);
456 double b1 = (2.44413*
q + 2.85619*
qq + 1.26661*
qqq)*
b0;
457 double b2 = (-1.4281*
qq - 1.26661*
qqq)*
b0;
459 double B = 1.0 - (
b1 +
b2 +
b3);
462 vigra_precondition(w >= 4,
463 "recursiveGaussianFilterLine(): line must have at least length 4.");
465 int kernelw = std::min(w-4, (
int)(4.0*sigma));
470 NumericTraits<typename SrcAccessor::value_type>::RealPromote
TempType;
493 for(x=3; x < w; ++x, ++
is)
505 for(x=w-4; x>=0; --x)
511 for(x=0; x < w; ++x, ++id)
592 vigra_precondition(scale >= 0,
593 "recursiveSmoothLine(): scale must be >= 0.\n");
595 double b = (scale == 0.0) ?
597 VIGRA_CSTD::exp(-1.0/scale);
679 vigra_precondition(scale > 0,
680 "recursiveFirstDerivativeLine(): scale must be > 0.\n");
687 NumericTraits<typename SrcAccessor::value_type>::RealPromote
691 std::vector<TempType>
vline(w);
694 double b = VIGRA_CSTD::exp(-1.0/scale);
695 double norm = (1.0 - b) * (1.0 - b) / 2.0 / b;
699 for(x=0; x<w; ++x, ++
is)
707 old = (1.0 / (1.0 - b)) *
as(
is);
711 for(x=w-1; x>=0; --x)
718 ad.set(DestTraits::fromRealPromote(
norm * (
line[x] +
old)),
id);
799 vigra_precondition(scale > 0,
800 "recursiveSecondDerivativeLine(): scale must be > 0.\n");
807 NumericTraits<typename SrcAccessor::value_type>::RealPromote
811 std::vector<TempType>
vline(w);
814 double b = VIGRA_CSTD::exp(-1.0/scale);
815 double a = -2.0 / (1.0 - b);
816 double norm = (1.0 - b) * (1.0 - b) * (1.0 - b) / (1.0 + b);
817 TempType old = detail::RequiresExplicitCast<TempType>::cast((1.0 / (1.0 - b)) *
as(
is));
820 for(x=0; x<w; ++x, ++
is)
823 old = detail::RequiresExplicitCast<TempType>::cast(
as(
is) + b *
old);
828 old = detail::RequiresExplicitCast<TempType>::cast((1.0 / (1.0 - b)) *
as(
is));
832 for(x=w-1; x>=0; --x)
837 TempType f = detail::RequiresExplicitCast<TempType>::cast(
old + a *
as(
is));
838 old = detail::RequiresExplicitCast<TempType>::cast(
as(
is) + b *
old);
839 ad.set(DestTraits::fromRealPromote(detail::RequiresExplicitCast<TempType>::cast(
norm * (
line[x] + f))),
id);
951 double b, BorderTreatmentMode
border)
960 typename SrcImageIterator::row_iterator
rs =
supperleft.rowIterator();
961 typename DestImageIterator::row_iterator rd =
dupperleft.rowIterator();
969template <
class SrcImageIterator,
class SrcAccessor,
970 class DestImageIterator,
class DestAccessor>
972recursiveFilterX(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
973 pair<DestImageIterator, DestAccessor> dest,
974 double b, BorderTreatmentMode border)
977 dest.first, dest.second, b, border);
980template <
class T1,
class S1,
984 MultiArrayView<2, T2, S2> dest,
985 double b, BorderTreatmentMode border)
987 vigra_precondition(src.shape() == dest.shape(),
988 "recursiveFilterX(): shape mismatch between input and output.");
990 destImage(dest), b, border);
999template <
class SrcImageIterator,
class SrcAccessor,
1000 class DestImageIterator,
class DestAccessor>
1002 SrcImageIterator slowerright, SrcAccessor as,
1003 DestImageIterator dupperleft, DestAccessor ad,
1004 double b1,
double b2)
1006 int w = slowerright.x - supperleft.x;
1007 int h = slowerright.y - supperleft.y;
1011 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
1013 typename SrcImageIterator::row_iterator rs = supperleft.rowIterator();
1014 typename DestImageIterator::row_iterator rd = dupperleft.rowIterator();
1022template <
class SrcImageIterator,
class SrcAccessor,
1023 class DestImageIterator,
class DestAccessor>
1025recursiveFilterX(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1026 pair<DestImageIterator, DestAccessor> dest,
1027 double b1,
double b2)
1030 dest.first, dest.second, b1, b2);
1033template <
class T1,
class S1,
1037 MultiArrayView<2, T2, S2> dest,
1038 double b1,
double b2)
1040 vigra_precondition(src.shape() == dest.shape(),
1041 "recursiveFilterX(): shape mismatch between input and output.");
1043 destImage(dest), b1, b2);
1136 typename SrcImageIterator::row_iterator
rs =
supperleft.rowIterator();
1137 typename DestImageIterator::row_iterator rd =
dupperleft.rowIterator();
1145template <
class SrcImageIterator,
class SrcAccessor,
1146 class DestImageIterator,
class DestAccessor>
1149 pair<DestImageIterator, DestAccessor> dest,
1153 dest.first, dest.second, sigma);
1156template <
class T1,
class S1,
1160 MultiArrayView<2, T2, S2> dest,
1163 vigra_precondition(src.shape() == dest.shape(),
1164 "recursiveGaussianFilterX(): shape mismatch between input and output.");
1166 destImage(dest), sigma);
1257 typename SrcImageIterator::row_iterator
rs =
supperleft.rowIterator();
1258 typename DestImageIterator::row_iterator rd =
dupperleft.rowIterator();
1266template <
class SrcImageIterator,
class SrcAccessor,
1267 class DestImageIterator,
class DestAccessor>
1269recursiveSmoothX(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1270 pair<DestImageIterator, DestAccessor> dest,
1274 dest.first, dest.second, scale);
1277template <
class T1,
class S1,
1281 MultiArrayView<2, T2, S2> dest,
1284 vigra_precondition(src.shape() == dest.shape(),
1285 "recursiveSmoothX(): shape mismatch between input and output.");
1287 destImage(dest), scale);
1397 double b, BorderTreatmentMode
border)
1406 typename SrcImageIterator::column_iterator
cs =
supperleft.columnIterator();
1407 typename DestImageIterator::column_iterator cd =
dupperleft.columnIterator();
1415template <
class SrcImageIterator,
class SrcAccessor,
1416 class DestImageIterator,
class DestAccessor>
1418recursiveFilterY(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1419 pair<DestImageIterator, DestAccessor> dest,
1420 double b, BorderTreatmentMode border)
1423 dest.first, dest.second, b, border);
1426template <
class T1,
class S1,
1430 MultiArrayView<2, T2, S2> dest,
1431 double b, BorderTreatmentMode border)
1433 vigra_precondition(src.shape() == dest.shape(),
1434 "recursiveFilterY(): shape mismatch between input and output.");
1436 destImage(dest), b, border);
1445template <
class SrcImageIterator,
class SrcAccessor,
1446 class DestImageIterator,
class DestAccessor>
1448 SrcImageIterator slowerright, SrcAccessor as,
1449 DestImageIterator dupperleft, DestAccessor ad,
1450 double b1,
double b2)
1452 int w = slowerright.x - supperleft.x;
1453 int h = slowerright.y - supperleft.y;
1457 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1459 typename SrcImageIterator::column_iterator cs = supperleft.columnIterator();
1460 typename DestImageIterator::column_iterator cd = dupperleft.columnIterator();
1468template <
class SrcImageIterator,
class SrcAccessor,
1469 class DestImageIterator,
class DestAccessor>
1471recursiveFilterY(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1472 pair<DestImageIterator, DestAccessor> dest,
1473 double b1,
double b2)
1476 dest.first, dest.second, b1, b2);
1479template <
class T1,
class S1,
1483 MultiArrayView<2, T2, S2> dest,
1484 double b1,
double b2)
1486 vigra_precondition(src.shape() == dest.shape(),
1487 "recursiveFilterY(): shape mismatch between input and output.");
1489 destImage(dest), b1, b2);
1582 typename SrcImageIterator::column_iterator
cs =
supperleft.columnIterator();
1583 typename DestImageIterator::column_iterator cd =
dupperleft.columnIterator();
1591template <
class SrcImageIterator,
class SrcAccessor,
1592 class DestImageIterator,
class DestAccessor>
1595 pair<DestImageIterator, DestAccessor> dest,
1599 dest.first, dest.second, sigma);
1602template <
class T1,
class S1,
1606 MultiArrayView<2, T2, S2> dest,
1609 vigra_precondition(src.shape() == dest.shape(),
1610 "recursiveGaussianFilterY(): shape mismatch between input and output.");
1612 destImage(dest), sigma);
1704 typename SrcImageIterator::column_iterator
cs =
supperleft.columnIterator();
1705 typename DestImageIterator::column_iterator cd =
dupperleft.columnIterator();
1713template <
class SrcImageIterator,
class SrcAccessor,
1714 class DestImageIterator,
class DestAccessor>
1716recursiveSmoothY(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1717 pair<DestImageIterator, DestAccessor> dest,
1721 dest.first, dest.second, scale);
1724template <
class T1,
class S1,
1728 MultiArrayView<2, T2, S2> dest,
1731 vigra_precondition(src.shape() == dest.shape(),
1732 "recursiveSmoothY(): shape mismatch between input and output.");
1734 destImage(dest), scale);
1826 typename SrcImageIterator::row_iterator
rs =
supperleft.rowIterator();
1827 typename DestImageIterator::row_iterator rd =
dupperleft.rowIterator();
1835template <
class SrcImageIterator,
class SrcAccessor,
1836 class DestImageIterator,
class DestAccessor>
1839 pair<DestImageIterator, DestAccessor> dest,
1843 dest.first, dest.second, scale);
1846template <
class T1,
class S1,
1850 MultiArrayView<2, T2, S2> dest,
1853 vigra_precondition(src.shape() == dest.shape(),
1854 "recursiveFirstDerivativeX(): shape mismatch between input and output.");
1856 destImage(dest), scale);
1948 typename SrcImageIterator::column_iterator
cs =
supperleft.columnIterator();
1949 typename DestImageIterator::column_iterator cd =
dupperleft.columnIterator();
1957template <
class SrcImageIterator,
class SrcAccessor,
1958 class DestImageIterator,
class DestAccessor>
1961 pair<DestImageIterator, DestAccessor> dest,
1965 dest.first, dest.second, scale);
1968template <
class T1,
class S1,
1972 MultiArrayView<2, T2, S2> dest,
1975 vigra_precondition(src.shape() == dest.shape(),
1976 "recursiveFirstDerivativeY(): shape mismatch between input and output.");
1978 destImage(dest), scale);
2070 typename SrcImageIterator::row_iterator
rs =
supperleft.rowIterator();
2071 typename DestImageIterator::row_iterator rd =
dupperleft.rowIterator();
2079template <
class SrcImageIterator,
class SrcAccessor,
2080 class DestImageIterator,
class DestAccessor>
2083 pair<DestImageIterator, DestAccessor> dest,
2087 dest.first, dest.second, scale);
2090template <
class T1,
class S1,
2094 MultiArrayView<2, T2, S2> dest,
2097 vigra_precondition(src.shape() == dest.shape(),
2098 "recursiveSecondDerivativeX(): shape mismatch between input and output.");
2100 destImage(dest), scale);
2192 typename SrcImageIterator::column_iterator
cs =
supperleft.columnIterator();
2193 typename DestImageIterator::column_iterator cd =
dupperleft.columnIterator();
2201template <
class SrcImageIterator,
class SrcAccessor,
2202 class DestImageIterator,
class DestAccessor>
2205 pair<DestImageIterator, DestAccessor> dest,
2209 dest.first, dest.second, scale);
2212template <
class T1,
class S1,
2216 MultiArrayView<2, T2, S2> dest,
2219 vigra_precondition(src.shape() == dest.shape(),
2220 "recursiveSecondDerivativeY(): shape mismatch between input and output.");
2222 destImage(dest), scale);
Class for a single RGB value.
Definition rgbvalue.hxx:128
iterator begin()
Definition tinyvector.hxx:861
void recursiveSmoothX(...)
Performs 1 dimensional recursive smoothing in x direction.
void recursiveFirstDerivativeY(...)
Recursively calculates the 1 dimensional first derivative in y direction.
void recursiveGaussianFilterY(...)
Compute 1 dimensional recursive approximation of Gaussian smoothing in y direction.
void recursiveSmoothY(...)
Performs 1 dimensional recursive smoothing in y direction.
void recursiveGaussianFilterX(...)
Compute 1 dimensional recursive approximation of Gaussian smoothing in y direction.
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition fftw3.hxx:1037
void recursiveGaussianFilterLine(...)
Compute a 1-dimensional recursive approximation of Gaussian smoothing.
void recursiveFilterY(...)
Performs 1 dimensional recursive filtering (1st and 2nd order) in y direction.
void recursiveSecondDerivativeLine(...)
Performs a 1 dimensional recursive convolution of the source signal.
void recursiveFirstDerivativeLine(...)
Performs a 1 dimensional recursive convolution of the source signal.
void recursiveFilterX(...)
Performs 1 dimensional recursive filtering (1st and 2nd order) in x direction.
void recursiveSmoothLine(...)
Convolves the image with a 1-dimensional exponential filter.
void recursiveFirstDerivativeX(...)
Recursively calculates the 1 dimensional first derivative in x direction.
void recursiveSecondDerivativeX(...)
Recursively calculates the 1 dimensional second derivative in x direction.
void recursiveFilterLine(...)
Performs a 1-dimensional recursive convolution of the source signal.
void recursiveSecondDerivativeY(...)
Recursively calculates the 1 dimensional second derivative in y direction.