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)
224 old = NumericTraits<TempType>::zero();
227 vigra_fail(
"recursiveFilterLine(): Unknown border treatment mode.\n");
230 for(x=0, is = istart; x < w; ++x, ++is)
232 old = TempType(as(is) + b * old);
237 if(border == BORDER_TREATMENT_REPEAT ||
238 border == BORDER_TREATMENT_AVOID)
241 old = TempType((1.0 / (1.0 - b)) * as(is));
243 else if(border == BORDER_TREATMENT_REFLECT)
247 else if(border == BORDER_TREATMENT_WRAP)
249 is = istart + kernelw - 1;
250 old = TempType((1.0 / (1.0 - b)) * as(is));
251 for(x = 0; x < kernelw; ++x, --is)
252 old = TempType(as(is) + b * old);
254 else if(border == BORDER_TREATMENT_CLIP)
256 old = NumericTraits<TempType>::zero();
261 if(border == BORDER_TREATMENT_CLIP)
267 for(x=w-1; x>=0; --x, --is, --id)
269 TempType f = TempType(b * old);
271 double norm = (1.0 - b) / (1.0 + b - bleft - bright);
274 ad.set(norm * (line[x] + f),
id);
277 else if(border == BORDER_TREATMENT_AVOID)
279 for(x=w-1; x >= kernelw; --x, --is, --id)
281 TempType f = TempType(b * old);
284 ad.set(DestTraits::fromRealPromote(RealPromote(norm * (line[x] + f))),
id);
289 for(x=w-1; x>=0; --x, --is, --id)
291 TempType f = TempType(b * old);
293 ad.set(DestTraits::fromRealPromote(RealPromote(norm * (line[x] + f))),
id);
304 template <
class SrcIterator,
class SrcAccessor,
305 class DestIterator,
class DestAccessor>
307 DestIterator
id, DestAccessor ad,
double b1,
double b2)
310 SrcIterator istart = is;
315 NumericTraits<typename SrcAccessor::value_type>::RealPromote TempType;
316 typedef NumericTraits<typename DestAccessor::value_type> DestTraits;
319 std::vector<TempType> vline(w+1);
320 typename std::vector<TempType>::iterator line = vline.begin();
322 double norm = 1.0 - b1 - b2;
323 double norm1 = (1.0 - b1 - b2) / (1.0 + b1 + b2);
324 double norm2 = norm *
norm;
328 int kernelw = std::min(w-1, std::max(8, (
int)(1.0 / norm + 0.5)));
330 line[kernelw] = as(is);
331 line[kernelw-1] = as(is);
332 for(x = kernelw - 2; x > 0; --x, --is)
334 line[x] = detail::RequiresExplicitCast<TempType>::cast(as(is) + b1 * line[x+1] + b2 * line[x+2]);
336 line[0] = detail::RequiresExplicitCast<TempType>::cast(as(is) + b1 * line[1] + b2 * line[2]);
338 line[1] = detail::RequiresExplicitCast<TempType>::cast(as(is) + b1 * line[0] + b2 * line[1]);
340 for(x=2; x < w; ++x, ++is)
342 line[x] = detail::RequiresExplicitCast<TempType>::cast(as(is) + b1 * line[x-1] + b2 * line[x-2]);
346 line[w-1] = detail::RequiresExplicitCast<TempType>::cast(norm1 * (line[w-1] + b1 * line[w-2] + b2 * line[w-3]));
347 line[w-2] = detail::RequiresExplicitCast<TempType>::cast(norm1 * (line[w-2] + b1 * line[w] + b2 * line[w-2]));
349 ad.set(line[w-1],
id);
351 ad.set(line[w-2],
id);
353 for(x=w-3; x>=0; --x, --id, --is)
355 line[x] = detail::RequiresExplicitCast<TempType>::cast(norm2 * line[x] + b1 * line[x+1] + b2 * line[x+2]);
443 template <
class SrcIterator,
class SrcAccessor,
444 class DestIterator,
class DestAccessor>
447 DestIterator
id, DestAccessor ad,
451 double q = 1.31564 * (
std::sqrt(1.0 + 0.490811 * sigma*sigma) - 1.0);
454 double b0 = 1.0/(1.57825 + 2.44413*q + 1.4281*qq + 0.422205*qqq);
455 double b1 = (2.44413*q + 2.85619*qq + 1.26661*qqq)*b0;
456 double b2 = (-1.4281*qq - 1.26661*qqq)*b0;
457 double b3 = 0.422205*qqq*b0;
458 double B = 1.0 - (b1 + b2 + b3);
461 vigra_precondition(w >= 4,
462 "recursiveGaussianFilterLine(): line must have at least length 4.");
464 int kernelw = std::min(w-4, (
int)(4.0*sigma));
469 NumericTraits<typename SrcAccessor::value_type>::RealPromote TempType;
470 typedef NumericTraits<typename DestAccessor::value_type> DestTraits;
473 std::vector<TempType> yforward(w);
475 std::vector<TempType> ybackward(w, 0.0);
478 for(x=kernelw; x>=0; --x)
480 ybackward[x] = detail::RequiresExplicitCast<TempType>::cast(B*as(is, x) + (b1*ybackward[x+1]+b2*ybackward[x+2]+b3*ybackward[x+3]));
484 yforward[0] = detail::RequiresExplicitCast<TempType>::cast(B*as(is) + (b1*ybackward[1]+b2*ybackward[2]+b3*ybackward[3]));
487 yforward[1] = detail::RequiresExplicitCast<TempType>::cast(B*as(is) + (b1*yforward[0]+b2*ybackward[1]+b3*ybackward[2]));
490 yforward[2] = detail::RequiresExplicitCast<TempType>::cast(B*as(is) + (b1*yforward[1]+b2*yforward[0]+b3*ybackward[1]));
493 for(x=3; x < w; ++x, ++is)
495 yforward[x] = detail::RequiresExplicitCast<TempType>::cast(B*as(is) + (b1*yforward[x-1]+b2*yforward[x-2]+b3*yforward[x-3]));
499 ybackward[w-1] = detail::RequiresExplicitCast<TempType>::cast(B*yforward[w-1] + (b1*yforward[w-2]+b2*yforward[w-3]+b3*yforward[w-4]));
501 ybackward[w-2] = detail::RequiresExplicitCast<TempType>::cast(B*yforward[w-2] + (b1*ybackward[w-1]+b2*yforward[w-2]+b3*yforward[w-3]));
503 ybackward[w-3] = detail::RequiresExplicitCast<TempType>::cast(B*yforward[w-3] + (b1*ybackward[w-2]+b2*ybackward[w-1]+b3*yforward[w-2]));
505 for(x=w-4; x>=0; --x)
507 ybackward[x] = detail::RequiresExplicitCast<TempType>::cast(B*yforward[x]+(b1*ybackward[x+1]+b2*ybackward[x+2]+b3*ybackward[x+3]));
511 for(x=0; x < w; ++x, ++id)
513 ad.set(ybackward[x],
id);
587 template <
class SrcIterator,
class SrcAccessor,
588 class DestIterator,
class DestAccessor>
591 DestIterator
id, DestAccessor ad,
double scale)
593 vigra_precondition(scale >= 0,
594 "recursiveSmoothLine(): scale must be >= 0.\n");
596 double b = (scale == 0.0) ?
676 template <
class SrcIterator,
class SrcAccessor,
677 class DestIterator,
class DestAccessor>
679 DestIterator
id, DestAccessor ad,
double scale)
681 vigra_precondition(scale > 0,
682 "recursiveFirstDerivativeLine(): scale must be > 0.\n");
689 NumericTraits<typename SrcAccessor::value_type>::RealPromote
691 typedef NumericTraits<typename DestAccessor::value_type> DestTraits;
693 std::vector<TempType> vline(w);
694 typename std::vector<TempType>::iterator line = vline.begin();
697 double norm = (1.0 - b) * (1.0 - b) / 2.0 / b;
698 TempType old = (1.0 / (1.0 - b)) * as(is);
701 for(x=0; x<w; ++x, ++is)
703 old = as(is) + b * old;
709 old = (1.0 / (1.0 - b)) * as(is);
713 for(x=w-1; x>=0; --x)
718 old = as(is) + b * old;
720 ad.set(DestTraits::fromRealPromote(norm * (line[x] + old)),
id);
797 template <
class SrcIterator,
class SrcAccessor,
798 class DestIterator,
class DestAccessor>
800 DestIterator
id, DestAccessor ad,
double scale)
802 vigra_precondition(scale > 0,
803 "recursiveSecondDerivativeLine(): scale must be > 0.\n");
810 NumericTraits<typename SrcAccessor::value_type>::RealPromote
812 typedef NumericTraits<typename DestAccessor::value_type> DestTraits;
814 std::vector<TempType> vline(w);
815 typename std::vector<TempType>::iterator line = vline.begin();
818 double a = -2.0 / (1.0 - b);
819 double norm = (1.0 - b) * (1.0 - b) * (1.0 - b) / (1.0 + b);
820 TempType old = detail::RequiresExplicitCast<TempType>::cast((1.0 / (1.0 - b)) * as(is));
823 for(x=0; x<w; ++x, ++is)
826 old = detail::RequiresExplicitCast<TempType>::cast(as(is) + b * old);
831 old = detail::RequiresExplicitCast<TempType>::cast((1.0 / (1.0 - b)) * as(is));
835 for(x=w-1; x>=0; --x)
840 TempType f = detail::RequiresExplicitCast<TempType>::cast(old + a * as(is));
841 old = detail::RequiresExplicitCast<TempType>::cast(as(is) + b * old);
842 ad.set(DestTraits::fromRealPromote(detail::RequiresExplicitCast<TempType>::cast(norm * (line[x] + f))),
id);
920 template <
class SrcImageIterator,
class SrcAccessor,
921 class DestImageIterator,
class DestAccessor>
923 SrcImageIterator slowerright, SrcAccessor as,
924 DestImageIterator dupperleft, DestAccessor ad,
925 double b, BorderTreatmentMode border)
927 int w = slowerright.x - supperleft.x;
928 int h = slowerright.y - supperleft.y;
932 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
934 typename SrcImageIterator::row_iterator rs = supperleft.rowIterator();
935 typename DestImageIterator::row_iterator rd = dupperleft.rowIterator();
943 template <
class SrcImageIterator,
class SrcAccessor,
944 class DestImageIterator,
class DestAccessor>
946 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
947 pair<DestImageIterator, DestAccessor> dest,
948 double b, BorderTreatmentMode border)
951 dest.first, dest.second, b, border);
960 template <
class SrcImageIterator,
class SrcAccessor,
961 class DestImageIterator,
class DestAccessor>
963 SrcImageIterator slowerright, SrcAccessor as,
964 DestImageIterator dupperleft, DestAccessor ad,
965 double b1,
double b2)
967 int w = slowerright.x - supperleft.x;
968 int h = slowerright.y - supperleft.y;
972 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
974 typename SrcImageIterator::row_iterator rs = supperleft.rowIterator();
975 typename DestImageIterator::row_iterator rd = dupperleft.rowIterator();
983 template <
class SrcImageIterator,
class SrcAccessor,
984 class DestImageIterator,
class DestAccessor>
986 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
987 pair<DestImageIterator, DestAccessor> dest,
988 double b1,
double b2)
991 dest.first, dest.second, b1, b2);
1053 template <
class SrcImageIterator,
class SrcAccessor,
1054 class DestImageIterator,
class DestAccessor>
1057 DestImageIterator dupperleft, DestAccessor ad,
1060 int w = slowerright.x - supperleft.x;
1061 int h = slowerright.y - supperleft.y;
1065 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
1067 typename SrcImageIterator::row_iterator rs = supperleft.rowIterator();
1068 typename DestImageIterator::row_iterator rd = dupperleft.rowIterator();
1076 template <
class SrcImageIterator,
class SrcAccessor,
1077 class DestImageIterator,
class DestAccessor>
1080 pair<DestImageIterator, DestAccessor> dest,
1084 dest.first, dest.second, sigma);
1143 template <
class SrcImageIterator,
class SrcAccessor,
1144 class DestImageIterator,
class DestAccessor>
1146 SrcImageIterator slowerright, SrcAccessor as,
1147 DestImageIterator dupperleft, DestAccessor ad,
1150 int w = slowerright.x - supperleft.x;
1151 int h = slowerright.y - supperleft.y;
1155 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
1157 typename SrcImageIterator::row_iterator rs = supperleft.rowIterator();
1158 typename DestImageIterator::row_iterator rd = dupperleft.rowIterator();
1166 template <
class SrcImageIterator,
class SrcAccessor,
1167 class DestImageIterator,
class DestAccessor>
1169 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1170 pair<DestImageIterator, DestAccessor> dest,
1174 dest. first, dest.second, scale);
1250 template <
class SrcImageIterator,
class SrcAccessor,
1251 class DestImageIterator,
class DestAccessor>
1253 SrcImageIterator slowerright, SrcAccessor as,
1254 DestImageIterator dupperleft, DestAccessor ad,
1255 double b, BorderTreatmentMode border)
1257 int w = slowerright.x - supperleft.x;
1258 int h = slowerright.y - supperleft.y;
1262 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1264 typename SrcImageIterator::column_iterator cs = supperleft.columnIterator();
1265 typename DestImageIterator::column_iterator cd = dupperleft.columnIterator();
1273 template <
class SrcImageIterator,
class SrcAccessor,
1274 class DestImageIterator,
class DestAccessor>
1276 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1277 pair<DestImageIterator, DestAccessor> dest,
1278 double b, BorderTreatmentMode border)
1281 dest.first, dest.second, b, border);
1290 template <
class SrcImageIterator,
class SrcAccessor,
1291 class DestImageIterator,
class DestAccessor>
1293 SrcImageIterator slowerright, SrcAccessor as,
1294 DestImageIterator dupperleft, DestAccessor ad,
1295 double b1,
double b2)
1297 int w = slowerright.x - supperleft.x;
1298 int h = slowerright.y - supperleft.y;
1302 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1304 typename SrcImageIterator::column_iterator cs = supperleft.columnIterator();
1305 typename DestImageIterator::column_iterator cd = dupperleft.columnIterator();
1313 template <
class SrcImageIterator,
class SrcAccessor,
1314 class DestImageIterator,
class DestAccessor>
1316 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1317 pair<DestImageIterator, DestAccessor> dest,
1318 double b1,
double b2)
1321 dest.first, dest.second, b1, b2);
1382 template <
class SrcImageIterator,
class SrcAccessor,
1383 class DestImageIterator,
class DestAccessor>
1386 DestImageIterator dupperleft, DestAccessor ad,
1389 int w = slowerright.x - supperleft.x;
1390 int h = slowerright.y - supperleft.y;
1394 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1396 typename SrcImageIterator::column_iterator cs = supperleft.columnIterator();
1397 typename DestImageIterator::column_iterator cd = dupperleft.columnIterator();
1405 template <
class SrcImageIterator,
class SrcAccessor,
1406 class DestImageIterator,
class DestAccessor>
1409 pair<DestImageIterator, DestAccessor> dest,
1413 dest.first, dest.second, sigma);
1472 template <
class SrcImageIterator,
class SrcAccessor,
1473 class DestImageIterator,
class DestAccessor>
1475 SrcImageIterator slowerright, SrcAccessor as,
1476 DestImageIterator dupperleft, DestAccessor ad,
1479 int w = slowerright.x - supperleft.x;
1480 int h = slowerright.y - supperleft.y;
1484 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1486 typename SrcImageIterator::column_iterator cs = supperleft.columnIterator();
1487 typename DestImageIterator::column_iterator cd = dupperleft.columnIterator();
1495 template <
class SrcImageIterator,
class SrcAccessor,
1496 class DestImageIterator,
class DestAccessor>
1498 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1499 pair<DestImageIterator, DestAccessor> dest,
1503 dest. first, dest.second, scale);
1562 template <
class SrcImageIterator,
class SrcAccessor,
1563 class DestImageIterator,
class DestAccessor>
1565 SrcImageIterator slowerright, SrcAccessor as,
1566 DestImageIterator dupperleft, DestAccessor ad,
1569 int w = slowerright.x - supperleft.x;
1570 int h = slowerright.y - supperleft.y;
1574 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
1576 typename SrcImageIterator::row_iterator rs = supperleft.rowIterator();
1577 typename DestImageIterator::row_iterator rd = dupperleft.rowIterator();
1585 template <
class SrcImageIterator,
class SrcAccessor,
1586 class DestImageIterator,
class DestAccessor>
1588 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1589 pair<DestImageIterator, DestAccessor> dest,
1593 dest. first, dest.second, scale);
1652 template <
class SrcImageIterator,
class SrcAccessor,
1653 class DestImageIterator,
class DestAccessor>
1655 SrcImageIterator slowerright, SrcAccessor as,
1656 DestImageIterator dupperleft, DestAccessor ad,
1659 int w = slowerright.x - supperleft.x;
1660 int h = slowerright.y - supperleft.y;
1664 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1666 typename SrcImageIterator::column_iterator cs = supperleft.columnIterator();
1667 typename DestImageIterator::column_iterator cd = dupperleft.columnIterator();
1675 template <
class SrcImageIterator,
class SrcAccessor,
1676 class DestImageIterator,
class DestAccessor>
1678 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1679 pair<DestImageIterator, DestAccessor> dest,
1683 dest. first, dest.second, scale);
1742 template <
class SrcImageIterator,
class SrcAccessor,
1743 class DestImageIterator,
class DestAccessor>
1745 SrcImageIterator slowerright, SrcAccessor as,
1746 DestImageIterator dupperleft, DestAccessor ad,
1749 int w = slowerright.x - supperleft.x;
1750 int h = slowerright.y - supperleft.y;
1754 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
1756 typename SrcImageIterator::row_iterator rs = supperleft.rowIterator();
1757 typename DestImageIterator::row_iterator rd = dupperleft.rowIterator();
1765 template <
class SrcImageIterator,
class SrcAccessor,
1766 class DestImageIterator,
class DestAccessor>
1768 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1769 pair<DestImageIterator, DestAccessor> dest,
1773 dest. first, dest.second, scale);
1832 template <
class SrcImageIterator,
class SrcAccessor,
1833 class DestImageIterator,
class DestAccessor>
1835 SrcImageIterator slowerright, SrcAccessor as,
1836 DestImageIterator dupperleft, DestAccessor ad,
1839 int w = slowerright.x - supperleft.x;
1840 int h = slowerright.y - supperleft.y;
1844 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1846 typename SrcImageIterator::column_iterator cs = supperleft.columnIterator();
1847 typename DestImageIterator::column_iterator cd = dupperleft.columnIterator();
1855 template <
class SrcImageIterator,
class SrcAccessor,
1856 class DestImageIterator,
class DestAccessor>
1858 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1859 pair<DestImageIterator, DestAccessor> dest,
1863 dest. first, dest.second, scale);
1871 #endif // VIGRA_RECURSIVECONVOLUTION_HXX