37 #ifndef VIGRA_RECURSIVECONVOLUTION_HXX
38 #define VIGRA_RECURSIVECONVOLUTION_HXX
42 #include "utilities.hxx"
43 #include "numerictraits.hxx"
44 #include "imageiteratoradapter.hxx"
45 #include "bordertreatment.hxx"
46 #include "array_vector.hxx"
164 template <
class SrcIterator,
class SrcAccessor,
165 class DestIterator,
class DestAccessor>
167 DestIterator
id, DestAccessor ad,
double b, BorderTreatmentMode border)
170 SrcIterator istart = is;
174 vigra_precondition(-1.0 < b && b < 1.0,
175 "recursiveFilterLine(): -1 < factor < 1 required.\n");
180 for(; is != isend; ++is, ++id)
187 double eps = 0.00001;
191 NumericTraits<typename SrcAccessor::value_type>::RealPromote TempType;
192 typedef NumericTraits<typename DestAccessor::value_type> DestTraits;
193 typedef typename DestTraits::RealPromote RealPromote;
196 std::vector<TempType> vline(w);
197 typename std::vector<TempType>::iterator line = vline.begin();
199 double norm = (1.0 - b) / (1.0 + b);
203 if(border == BORDER_TREATMENT_REPEAT ||
204 border == BORDER_TREATMENT_AVOID)
206 old = TempType((1.0 / (1.0 - b)) * as(is));
208 else if(border == BORDER_TREATMENT_REFLECT)
211 old = TempType((1.0 / (1.0 - b)) * as(is));
212 for(x = 0; x < kernelw; ++x, --is)
213 old = TempType(as(is) + b * old);
215 else if(border == BORDER_TREATMENT_WRAP)
217 is = isend - kernelw;
218 old = TempType((1.0 / (1.0 - b)) * as(is));
219 for(x = 0; x < kernelw; ++x, ++is)
220 old = TempType(as(is) + b * old);
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();
234 for(x=0, is = istart; x < w; ++x, ++is)
236 old = TempType(as(is) + b * old);
241 if(border == BORDER_TREATMENT_REPEAT ||
242 border == BORDER_TREATMENT_AVOID)
245 old = TempType((1.0 / (1.0 - b)) * as(is));
247 else if(border == BORDER_TREATMENT_REFLECT)
251 else if(border == BORDER_TREATMENT_WRAP)
253 is = istart + kernelw - 1;
254 old = TempType((1.0 / (1.0 - b)) * as(is));
255 for(x = 0; x < kernelw; ++x, --is)
256 old = TempType(as(is) + b * old);
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)
274 TempType f = TempType(b * old);
276 double norm = (1.0 - b) / (1.0 + b - bleft - bright);
279 ad.set(norm * (line[x] + f),
id);
282 else if(border == BORDER_TREATMENT_AVOID)
284 for(x=w-1; x >= kernelw; --x, --is, --id)
286 TempType f = TempType(b * old);
289 ad.set(DestTraits::fromRealPromote(RealPromote(norm * (line[x] + f))),
id);
294 for(x=w-1; x>=0; --x, --is, --id)
296 TempType f = TempType(b * old);
298 ad.set(DestTraits::fromRealPromote(RealPromote(norm * (line[x] + f))),
id);
309 template <
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;
319 typedef NumericTraits<typename DestAccessor::value_type> DestTraits;
322 std::vector<TempType> vline(w+1);
323 typename std::vector<TempType>::iterator line = vline.begin();
325 double norm = 1.0 - b1 - b2;
326 double norm1 = (1.0 - b1 - b2) / (1.0 + b1 + b2);
327 double norm2 = norm *
norm;
331 int kernelw = std::min(w-1, std::max(8, (
int)(1.0 / norm + 0.5)));
333 line[kernelw] = as(is);
334 line[kernelw-1] = as(is);
335 for(x = kernelw - 2; x > 0; --x, --is)
337 line[x] = detail::RequiresExplicitCast<TempType>::cast(as(is) + b1 * line[x+1] + b2 * line[x+2]);
339 line[0] = detail::RequiresExplicitCast<TempType>::cast(as(is) + b1 * line[1] + b2 * line[2]);
341 line[1] = detail::RequiresExplicitCast<TempType>::cast(as(is) + b1 * line[0] + b2 * line[1]);
343 for(x=2; x < w; ++x, ++is)
345 line[x] = detail::RequiresExplicitCast<TempType>::cast(as(is) + b1 * line[x-1] + b2 * line[x-2]);
349 line[w-1] = detail::RequiresExplicitCast<TempType>::cast(norm1 * (line[w-1] + b1 * line[w-2] + b2 * line[w-3]));
350 line[w-2] = detail::RequiresExplicitCast<TempType>::cast(norm1 * (line[w-2] + b1 * line[w] + b2 * line[w-2]));
352 ad.set(line[w-1],
id);
354 ad.set(line[w-2],
id);
356 for(x=w-3; x>=0; --x, --id, --is)
358 line[x] = detail::RequiresExplicitCast<TempType>::cast(norm2 * line[x] + b1 * line[x+1] + b2 * line[x+2]);
446 template <
class SrcIterator,
class SrcAccessor,
447 class DestIterator,
class DestAccessor>
450 DestIterator
id, DestAccessor ad,
454 double q = 1.31564 * (
std::sqrt(1.0 + 0.490811 * sigma*sigma) - 1.0);
457 double b0 = 1.0/(1.57825 + 2.44413*q + 1.4281*qq + 0.422205*qqq);
458 double b1 = (2.44413*q + 2.85619*qq + 1.26661*qqq)*b0;
459 double b2 = (-1.4281*qq - 1.26661*qqq)*b0;
460 double b3 = 0.422205*qqq*b0;
461 double B = 1.0 - (b1 + b2 + b3);
464 vigra_precondition(w >= 4,
465 "recursiveGaussianFilterLine(): line must have at least length 4.");
467 int kernelw = std::min(w-4, (
int)(4.0*sigma));
472 NumericTraits<typename SrcAccessor::value_type>::RealPromote TempType;
473 typedef NumericTraits<typename DestAccessor::value_type> DestTraits;
476 std::vector<TempType> yforward(w);
478 std::vector<TempType> ybackward(w, 0.0);
481 for(x=kernelw; x>=0; --x)
483 ybackward[x] = detail::RequiresExplicitCast<TempType>::cast(B*as(is, x) + (b1*ybackward[x+1]+b2*ybackward[x+2]+b3*ybackward[x+3]));
487 yforward[0] = detail::RequiresExplicitCast<TempType>::cast(B*as(is) + (b1*ybackward[1]+b2*ybackward[2]+b3*ybackward[3]));
490 yforward[1] = detail::RequiresExplicitCast<TempType>::cast(B*as(is) + (b1*yforward[0]+b2*ybackward[1]+b3*ybackward[2]));
493 yforward[2] = detail::RequiresExplicitCast<TempType>::cast(B*as(is) + (b1*yforward[1]+b2*yforward[0]+b3*ybackward[1]));
496 for(x=3; x < w; ++x, ++is)
498 yforward[x] = detail::RequiresExplicitCast<TempType>::cast(B*as(is) + (b1*yforward[x-1]+b2*yforward[x-2]+b3*yforward[x-3]));
502 ybackward[w-1] = detail::RequiresExplicitCast<TempType>::cast(B*yforward[w-1] + (b1*yforward[w-2]+b2*yforward[w-3]+b3*yforward[w-4]));
504 ybackward[w-2] = detail::RequiresExplicitCast<TempType>::cast(B*yforward[w-2] + (b1*ybackward[w-1]+b2*yforward[w-2]+b3*yforward[w-3]));
506 ybackward[w-3] = detail::RequiresExplicitCast<TempType>::cast(B*yforward[w-3] + (b1*ybackward[w-2]+b2*ybackward[w-1]+b3*yforward[w-2]));
508 for(x=w-4; x>=0; --x)
510 ybackward[x] = detail::RequiresExplicitCast<TempType>::cast(B*yforward[x]+(b1*ybackward[x+1]+b2*ybackward[x+2]+b3*ybackward[x+3]));
514 for(x=0; x < w; ++x, ++id)
516 ad.set(ybackward[x],
id);
590 template <
class SrcIterator,
class SrcAccessor,
591 class DestIterator,
class DestAccessor>
594 DestIterator
id, DestAccessor ad,
double scale)
596 vigra_precondition(scale >= 0,
597 "recursiveSmoothLine(): scale must be >= 0.\n");
599 double b = (scale == 0.0) ?
679 template <
class SrcIterator,
class SrcAccessor,
680 class DestIterator,
class DestAccessor>
682 DestIterator
id, DestAccessor ad,
double scale)
684 vigra_precondition(scale > 0,
685 "recursiveFirstDerivativeLine(): scale must be > 0.\n");
692 NumericTraits<typename SrcAccessor::value_type>::RealPromote
694 typedef NumericTraits<typename DestAccessor::value_type> DestTraits;
696 std::vector<TempType> vline(w);
697 typename std::vector<TempType>::iterator line = vline.begin();
700 double norm = (1.0 - b) * (1.0 - b) / 2.0 / b;
701 TempType old = (1.0 / (1.0 - b)) * as(is);
704 for(x=0; x<w; ++x, ++is)
706 old = as(is) + b * old;
712 old = (1.0 / (1.0 - b)) * as(is);
716 for(x=w-1; x>=0; --x)
721 old = as(is) + b * old;
723 ad.set(DestTraits::fromRealPromote(norm * (line[x] + old)),
id);
800 template <
class SrcIterator,
class SrcAccessor,
801 class DestIterator,
class DestAccessor>
803 DestIterator
id, DestAccessor ad,
double scale)
805 vigra_precondition(scale > 0,
806 "recursiveSecondDerivativeLine(): scale must be > 0.\n");
813 NumericTraits<typename SrcAccessor::value_type>::RealPromote
815 typedef NumericTraits<typename DestAccessor::value_type> DestTraits;
817 std::vector<TempType> vline(w);
818 typename std::vector<TempType>::iterator line = vline.begin();
821 double a = -2.0 / (1.0 - b);
822 double norm = (1.0 - b) * (1.0 - b) * (1.0 - b) / (1.0 + b);
823 TempType old = detail::RequiresExplicitCast<TempType>::cast((1.0 / (1.0 - b)) * as(is));
826 for(x=0; x<w; ++x, ++is)
829 old = detail::RequiresExplicitCast<TempType>::cast(as(is) + b * old);
834 old = detail::RequiresExplicitCast<TempType>::cast((1.0 / (1.0 - b)) * as(is));
838 for(x=w-1; x>=0; --x)
843 TempType f = detail::RequiresExplicitCast<TempType>::cast(old + a * as(is));
844 old = detail::RequiresExplicitCast<TempType>::cast(as(is) + b * old);
845 ad.set(DestTraits::fromRealPromote(detail::RequiresExplicitCast<TempType>::cast(norm * (line[x] + f))),
id);
923 template <
class SrcImageIterator,
class SrcAccessor,
924 class DestImageIterator,
class DestAccessor>
926 SrcImageIterator slowerright, SrcAccessor as,
927 DestImageIterator dupperleft, DestAccessor ad,
928 double b, BorderTreatmentMode border)
930 int w = slowerright.x - supperleft.x;
931 int h = slowerright.y - supperleft.y;
935 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
937 typename SrcImageIterator::row_iterator rs = supperleft.rowIterator();
938 typename DestImageIterator::row_iterator rd = dupperleft.rowIterator();
946 template <
class SrcImageIterator,
class SrcAccessor,
947 class DestImageIterator,
class DestAccessor>
949 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
950 pair<DestImageIterator, DestAccessor> dest,
951 double b, BorderTreatmentMode border)
954 dest.first, dest.second, b, border);
963 template <
class SrcImageIterator,
class SrcAccessor,
964 class DestImageIterator,
class DestAccessor>
966 SrcImageIterator slowerright, SrcAccessor as,
967 DestImageIterator dupperleft, DestAccessor ad,
968 double b1,
double b2)
970 int w = slowerright.x - supperleft.x;
971 int h = slowerright.y - supperleft.y;
975 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
977 typename SrcImageIterator::row_iterator rs = supperleft.rowIterator();
978 typename DestImageIterator::row_iterator rd = dupperleft.rowIterator();
986 template <
class SrcImageIterator,
class SrcAccessor,
987 class DestImageIterator,
class DestAccessor>
989 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
990 pair<DestImageIterator, DestAccessor> dest,
991 double b1,
double b2)
994 dest.first, dest.second, b1, b2);
1056 template <
class SrcImageIterator,
class SrcAccessor,
1057 class DestImageIterator,
class DestAccessor>
1060 DestImageIterator dupperleft, DestAccessor ad,
1063 int w = slowerright.x - supperleft.x;
1064 int h = slowerright.y - supperleft.y;
1068 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
1070 typename SrcImageIterator::row_iterator rs = supperleft.rowIterator();
1071 typename DestImageIterator::row_iterator rd = dupperleft.rowIterator();
1079 template <
class SrcImageIterator,
class SrcAccessor,
1080 class DestImageIterator,
class DestAccessor>
1083 pair<DestImageIterator, DestAccessor> dest,
1087 dest.first, dest.second, sigma);
1146 template <
class SrcImageIterator,
class SrcAccessor,
1147 class DestImageIterator,
class DestAccessor>
1149 SrcImageIterator slowerright, SrcAccessor as,
1150 DestImageIterator dupperleft, DestAccessor ad,
1153 int w = slowerright.x - supperleft.x;
1154 int h = slowerright.y - supperleft.y;
1158 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
1160 typename SrcImageIterator::row_iterator rs = supperleft.rowIterator();
1161 typename DestImageIterator::row_iterator rd = dupperleft.rowIterator();
1169 template <
class SrcImageIterator,
class SrcAccessor,
1170 class DestImageIterator,
class DestAccessor>
1172 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1173 pair<DestImageIterator, DestAccessor> dest,
1177 dest. first, dest.second, scale);
1253 template <
class SrcImageIterator,
class SrcAccessor,
1254 class DestImageIterator,
class DestAccessor>
1256 SrcImageIterator slowerright, SrcAccessor as,
1257 DestImageIterator dupperleft, DestAccessor ad,
1258 double b, BorderTreatmentMode border)
1260 int w = slowerright.x - supperleft.x;
1261 int h = slowerright.y - supperleft.y;
1265 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1267 typename SrcImageIterator::column_iterator cs = supperleft.columnIterator();
1268 typename DestImageIterator::column_iterator cd = dupperleft.columnIterator();
1276 template <
class SrcImageIterator,
class SrcAccessor,
1277 class DestImageIterator,
class DestAccessor>
1279 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1280 pair<DestImageIterator, DestAccessor> dest,
1281 double b, BorderTreatmentMode border)
1284 dest.first, dest.second, b, border);
1293 template <
class SrcImageIterator,
class SrcAccessor,
1294 class DestImageIterator,
class DestAccessor>
1296 SrcImageIterator slowerright, SrcAccessor as,
1297 DestImageIterator dupperleft, DestAccessor ad,
1298 double b1,
double b2)
1300 int w = slowerright.x - supperleft.x;
1301 int h = slowerright.y - supperleft.y;
1305 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1307 typename SrcImageIterator::column_iterator cs = supperleft.columnIterator();
1308 typename DestImageIterator::column_iterator cd = dupperleft.columnIterator();
1316 template <
class SrcImageIterator,
class SrcAccessor,
1317 class DestImageIterator,
class DestAccessor>
1319 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1320 pair<DestImageIterator, DestAccessor> dest,
1321 double b1,
double b2)
1324 dest.first, dest.second, b1, b2);
1385 template <
class SrcImageIterator,
class SrcAccessor,
1386 class DestImageIterator,
class DestAccessor>
1389 DestImageIterator dupperleft, DestAccessor ad,
1392 int w = slowerright.x - supperleft.x;
1393 int h = slowerright.y - supperleft.y;
1397 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1399 typename SrcImageIterator::column_iterator cs = supperleft.columnIterator();
1400 typename DestImageIterator::column_iterator cd = dupperleft.columnIterator();
1408 template <
class SrcImageIterator,
class SrcAccessor,
1409 class DestImageIterator,
class DestAccessor>
1412 pair<DestImageIterator, DestAccessor> dest,
1416 dest.first, dest.second, sigma);
1475 template <
class SrcImageIterator,
class SrcAccessor,
1476 class DestImageIterator,
class DestAccessor>
1478 SrcImageIterator slowerright, SrcAccessor as,
1479 DestImageIterator dupperleft, DestAccessor ad,
1482 int w = slowerright.x - supperleft.x;
1483 int h = slowerright.y - supperleft.y;
1487 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1489 typename SrcImageIterator::column_iterator cs = supperleft.columnIterator();
1490 typename DestImageIterator::column_iterator cd = dupperleft.columnIterator();
1498 template <
class SrcImageIterator,
class SrcAccessor,
1499 class DestImageIterator,
class DestAccessor>
1501 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1502 pair<DestImageIterator, DestAccessor> dest,
1506 dest. first, dest.second, scale);
1565 template <
class SrcImageIterator,
class SrcAccessor,
1566 class DestImageIterator,
class DestAccessor>
1568 SrcImageIterator slowerright, SrcAccessor as,
1569 DestImageIterator dupperleft, DestAccessor ad,
1572 int w = slowerright.x - supperleft.x;
1573 int h = slowerright.y - supperleft.y;
1577 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
1579 typename SrcImageIterator::row_iterator rs = supperleft.rowIterator();
1580 typename DestImageIterator::row_iterator rd = dupperleft.rowIterator();
1588 template <
class SrcImageIterator,
class SrcAccessor,
1589 class DestImageIterator,
class DestAccessor>
1591 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1592 pair<DestImageIterator, DestAccessor> dest,
1596 dest. first, dest.second, scale);
1655 template <
class SrcImageIterator,
class SrcAccessor,
1656 class DestImageIterator,
class DestAccessor>
1658 SrcImageIterator slowerright, SrcAccessor as,
1659 DestImageIterator dupperleft, DestAccessor ad,
1662 int w = slowerright.x - supperleft.x;
1663 int h = slowerright.y - supperleft.y;
1667 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1669 typename SrcImageIterator::column_iterator cs = supperleft.columnIterator();
1670 typename DestImageIterator::column_iterator cd = dupperleft.columnIterator();
1678 template <
class SrcImageIterator,
class SrcAccessor,
1679 class DestImageIterator,
class DestAccessor>
1681 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1682 pair<DestImageIterator, DestAccessor> dest,
1686 dest. first, dest.second, scale);
1745 template <
class SrcImageIterator,
class SrcAccessor,
1746 class DestImageIterator,
class DestAccessor>
1748 SrcImageIterator slowerright, SrcAccessor as,
1749 DestImageIterator dupperleft, DestAccessor ad,
1752 int w = slowerright.x - supperleft.x;
1753 int h = slowerright.y - supperleft.y;
1757 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
1759 typename SrcImageIterator::row_iterator rs = supperleft.rowIterator();
1760 typename DestImageIterator::row_iterator rd = dupperleft.rowIterator();
1768 template <
class SrcImageIterator,
class SrcAccessor,
1769 class DestImageIterator,
class DestAccessor>
1771 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1772 pair<DestImageIterator, DestAccessor> dest,
1776 dest. first, dest.second, scale);
1835 template <
class SrcImageIterator,
class SrcAccessor,
1836 class DestImageIterator,
class DestAccessor>
1838 SrcImageIterator slowerright, SrcAccessor as,
1839 DestImageIterator dupperleft, DestAccessor ad,
1842 int w = slowerright.x - supperleft.x;
1843 int h = slowerright.y - supperleft.y;
1847 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1849 typename SrcImageIterator::column_iterator cs = supperleft.columnIterator();
1850 typename DestImageIterator::column_iterator cd = dupperleft.columnIterator();
1858 template <
class SrcImageIterator,
class SrcAccessor,
1859 class DestImageIterator,
class DestAccessor>
1861 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1862 pair<DestImageIterator, DestAccessor> dest,
1866 dest. first, dest.second, scale);
1874 #endif // VIGRA_RECURSIVECONVOLUTION_HXX