37 #ifndef VIGRA_BOUNDARYTENSOR_HXX
38 #define VIGRA_BOUNDARYTENSOR_HXX
42 #include "utilities.hxx"
43 #include "array_vector.hxx"
44 #include "basicimage.hxx"
45 #include "combineimages.hxx"
46 #include "numerictraits.hxx"
47 #include "convolution.hxx"
55 typedef ArrayVector<Kernel1D<double> > KernelArray;
57 template <
class KernelArray>
59 initGaussianPolarFilters1(
double std_dev, KernelArray & k)
61 typedef typename KernelArray::value_type Kernel;
62 typedef typename Kernel::iterator iterator;
64 vigra_precondition(std_dev >= 0.0,
65 "initGaussianPolarFilter1(): "
66 "Standard deviation must be >= 0.");
70 int radius = (int)(4.0*std_dev + 0.5);
71 std_dev *= 1.08179074376;
73 double a = 0.558868151788 / VIGRA_CSTD::pow(std_dev, 5);
74 double b = -2.04251639729 / VIGRA_CSTD::pow(std_dev, 3);
75 double sigma22 = -0.5 / std_dev / std_dev;
78 for(
unsigned int i=0; i<k.size(); ++i)
80 k[i].initExplicitly(-radius, radius);
81 k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
85 iterator c = k[0].center();
86 for(ix=-radius; ix<=radius; ++ix)
88 double x = (double)ix;
93 for(ix=-radius; ix<=radius; ++ix)
95 double x = (double)ix;
101 for(ix=-radius; ix<=radius; ++ix)
103 double x = (double)ix;
108 for(ix=-radius; ix<=radius; ++ix)
110 double x = (double)ix;
115 template <
class KernelArray>
117 initGaussianPolarFilters2(
double std_dev, KernelArray & k)
119 typedef typename KernelArray::value_type Kernel;
120 typedef typename Kernel::iterator iterator;
122 vigra_precondition(std_dev >= 0.0,
123 "initGaussianPolarFilter2(): "
124 "Standard deviation must be >= 0.");
128 int radius = (int)(4.0*std_dev + 0.5);
130 double sigma2 = std_dev*std_dev;
131 double sigma22 = -0.5 / sigma2;
133 for(
unsigned int i=0; i<k.size(); ++i)
135 k[i].initExplicitly(-radius, radius);
136 k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
140 iterator c = k[0].center();
141 for(ix=-radius; ix<=radius; ++ix)
143 double x = (double)ix;
148 double f1 = f / sigma2;
149 for(ix=-radius; ix<=radius; ++ix)
151 double x = (double)ix;
156 double f2 = f / (sigma2 * sigma2);
157 for(ix=-radius; ix<=radius; ++ix)
159 double x = (double)ix;
164 template <
class KernelArray>
166 initGaussianPolarFilters3(
double std_dev, KernelArray & k)
168 typedef typename KernelArray::value_type Kernel;
169 typedef typename Kernel::iterator iterator;
171 vigra_precondition(std_dev >= 0.0,
172 "initGaussianPolarFilter3(): "
173 "Standard deviation must be >= 0.");
177 int radius = (int)(4.0*std_dev + 0.5);
178 std_dev *= 1.15470053838;
179 double sigma22 = -0.5 / std_dev / std_dev;
181 double a = 0.883887052922 / VIGRA_CSTD::pow(std_dev, 5);
183 for(
unsigned int i=0; i<k.size(); ++i)
185 k[i].initExplicitly(-radius, radius);
186 k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
192 iterator c = k[0].center();
193 for(ix=-radius; ix<=radius; ++ix)
195 double x = (double)ix;
200 for(ix=-radius; ix<=radius; ++ix)
202 double x = (double)ix;
208 for(ix=-radius; ix<=radius; ++ix)
210 double x = (double)ix;
215 for(ix=-radius; ix<=radius; ++ix)
217 double x = (double)ix;
222 template <
class SrcIterator,
class SrcAccessor,
223 class DestIterator,
class DestAccessor>
225 evenPolarFilters(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
226 DestIterator dupperleft, DestAccessor dest,
227 double scale,
bool noLaplacian)
229 vigra_precondition(dest.size(dupperleft) == 3,
230 "evenPolarFilters(): image for even output must have 3 bands.");
232 int w = slowerright.x - supperleft.x;
233 int h = slowerright.y - supperleft.y;
236 NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
237 typedef BasicImage<TinyVector<TmpType, 3> > TmpImage;
238 typedef typename TmpImage::traverser TmpTraverser;
242 initGaussianPolarFilters2(scale, k2);
245 VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t.accessor());
247 destImage(t, tmpBand), k2[2], k2[0]);
250 destImage(t, tmpBand), k2[1], k2[1]);
253 destImage(t, tmpBand), k2[0], k2[2]);
256 TmpTraverser tul(t.upperLeft());
257 TmpTraverser tlr(t.lowerRight());
258 for(; tul.y != tlr.y; ++tul.y, ++dupperleft.y)
260 typename TmpTraverser::row_iterator tr = tul.rowIterator();
261 typename TmpTraverser::row_iterator trend = tr + w;
262 typename DestIterator::row_iterator d = dupperleft.rowIterator();
265 for(; tr != trend; ++tr, ++d)
267 TmpType v = detail::RequiresExplicitCast<TmpType>::cast(0.5*
sq((*tr)[0]-(*tr)[2]) + 2.0*
sq((*tr)[1]));
268 dest.setComponent(v, d, 0);
269 dest.setComponent(0, d, 1);
270 dest.setComponent(v, d, 2);
275 for(; tr != trend; ++tr, ++d)
277 dest.setComponent(
sq((*tr)[0]) +
sq((*tr)[1]), d, 0);
278 dest.setComponent(-(*tr)[1] * ((*tr)[0] + (*tr)[2]), d, 1);
279 dest.setComponent(
sq((*tr)[1]) +
sq((*tr)[2]), d, 2);
285 template <
class SrcIterator,
class SrcAccessor,
286 class DestIterator,
class DestAccessor>
288 oddPolarFilters(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
289 DestIterator dupperleft, DestAccessor dest,
290 double scale,
bool addResult)
292 vigra_precondition(dest.size(dupperleft) == 3,
293 "oddPolarFilters(): image for odd output must have 3 bands.");
295 int w = slowerright.x - supperleft.x;
296 int h = slowerright.y - supperleft.y;
299 NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
300 typedef BasicImage<TinyVector<TmpType, 4> > TmpImage;
301 typedef typename TmpImage::traverser TmpTraverser;
304 detail::KernelArray k1;
305 detail::initGaussianPolarFilters1(scale, k1);
308 VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t.accessor());
310 destImage(t, tmpBand), k1[3], k1[0]);
313 destImage(t, tmpBand), k1[2], k1[1]);
316 destImage(t, tmpBand), k1[1], k1[2]);
319 destImage(t, tmpBand), k1[0], k1[3]);
322 TmpTraverser tul(t.upperLeft());
323 TmpTraverser tlr(t.lowerRight());
324 for(; tul.y != tlr.y; ++tul.y, ++dupperleft.y)
326 typename TmpTraverser::row_iterator tr = tul.rowIterator();
327 typename TmpTraverser::row_iterator trend = tr + w;
328 typename DestIterator::row_iterator d = dupperleft.rowIterator();
331 for(; tr != trend; ++tr, ++d)
333 TmpType d0 = (*tr)[0] + (*tr)[2];
334 TmpType d1 = -(*tr)[1] - (*tr)[3];
336 dest.setComponent(dest.getComponent(d, 0) +
sq(d0), d, 0);
337 dest.setComponent(dest.getComponent(d, 1) + d0 * d1, d, 1);
338 dest.setComponent(dest.getComponent(d, 2) +
sq(d1), d, 2);
343 for(; tr != trend; ++tr, ++d)
345 TmpType d0 = (*tr)[0] + (*tr)[2];
346 TmpType d1 = -(*tr)[1] - (*tr)[3];
348 dest.setComponent(
sq(d0), d, 0);
349 dest.setComponent(d0 * d1, d, 1);
350 dest.setComponent(
sq(d1), d, 2);
423 template <
class SrcIterator,
class SrcAccessor,
424 class DestIterator,
class DestAccessor>
426 DestIterator dupperleft, DestAccessor dest,
427 double scale,
unsigned int xorder,
unsigned int yorder)
429 unsigned int order = xorder + yorder;
431 vigra_precondition(order <= 2,
432 "rieszTransformOfLOG(): can only compute Riesz transforms up to order 2.");
433 vigra_precondition(scale > 0.0,
434 "rieszTransformOfLOG(): scale must be positive.");
436 int w = slowerright.x - supperleft.x;
437 int h = slowerright.y - supperleft.y;
439 typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
440 typedef BasicImage<TmpType> TmpImage;
446 detail::KernelArray k2;
447 detail::initGaussianPolarFilters2(scale, k2);
449 TmpImage t1(w, h), t2(w, h);
452 destImage(t1), k2[2], k2[0]);
454 destImage(t2), k2[0], k2[2]);
456 destIter(dupperleft, dest), std::plus<TmpType>());
461 detail::KernelArray k1;
462 detail::initGaussianPolarFilters1(scale, k1);
464 TmpImage t1(w, h), t2(w, h);
469 destImage(t1), k1[3], k1[0]);
471 destImage(t2), k1[1], k1[2]);
476 destImage(t1), k1[0], k1[3]);
478 destImage(t2), k1[2], k1[1]);
481 destIter(dupperleft, dest), std::plus<TmpType>());
486 detail::KernelArray k2;
487 detail::initGaussianPolarFilters2(scale, k2);
490 destIter(dupperleft, dest), k2[xorder], k2[yorder]);
496 detail::KernelArray k3;
497 detail::initGaussianPolarFilters3(scale, k3);
499 TmpImage t1(w, h), t2(w, h);
504 destImage(t1), k3[3], k3[0]);
506 destImage(t2), k3[1], k3[2]);
511 destImage(t1), k3[0], k3[3]);
513 destImage(t2), k3[2], k3[1]);
516 destIter(dupperleft, dest), std::minus<TmpType>());
522 template <
class SrcIterator,
class SrcAccessor,
523 class DestIterator,
class DestAccessor>
526 pair<DestIterator, DestAccessor> dest,
527 double scale,
unsigned int xorder,
unsigned int yorder)
530 scale, xorder, yorder);
599 template <
class SrcIterator,
class SrcAccessor,
600 class DestIterator,
class DestAccessor>
601 void boundaryTensor(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
602 DestIterator dupperleft, DestAccessor dest,
605 vigra_precondition(dest.size(dupperleft) == 3,
606 "boundaryTensor(): image for even output must have 3 bands.");
607 vigra_precondition(scale > 0.0,
608 "boundaryTensor(): scale must be positive.");
610 detail::evenPolarFilters(supperleft, slowerright, src,
611 dupperleft, dest, scale,
false);
612 detail::oddPolarFilters(supperleft, slowerright, src,
613 dupperleft, dest, scale,
true);
616 template <
class SrcIterator,
class SrcAccessor,
617 class DestIterator,
class DestAccessor>
619 void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src,
620 pair<DestIterator, DestAccessor> dest,
624 dest.first, dest.second, scale);
662 template <
class SrcIterator,
class SrcAccessor,
663 class DestIterator,
class DestAccessor>
664 void boundaryTensor1(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
665 DestIterator dupperleft, DestAccessor dest,
668 vigra_precondition(dest.size(dupperleft) == 3,
669 "boundaryTensor1(): image for even output must have 3 bands.");
670 vigra_precondition(scale > 0.0,
671 "boundaryTensor1(): scale must be positive.");
673 detail::evenPolarFilters(supperleft, slowerright, src,
674 dupperleft, dest, scale,
true);
675 detail::oddPolarFilters(supperleft, slowerright, src,
676 dupperleft, dest, scale,
true);
679 template <
class SrcIterator,
class SrcAccessor,
680 class DestIterator,
class DestAccessor>
682 void boundaryTensor1(triple<SrcIterator, SrcIterator, SrcAccessor> src,
683 pair<DestIterator, DestAccessor> dest,
687 dest.first, dest.second, scale);
699 template <
class SrcIterator,
class SrcAccessor,
700 class DestIteratorEven,
class DestAccessorEven,
701 class DestIteratorOdd,
class DestAccessorOdd>
702 void boundaryTensor3(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor sa,
703 DestIteratorEven dupperleft_even, DestAccessorEven
even,
704 DestIteratorOdd dupperleft_odd, DestAccessorOdd
odd,
707 vigra_precondition(even.size(dupperleft_even) == 3,
708 "boundaryTensor3(): image for even output must have 3 bands.");
709 vigra_precondition(odd.size(dupperleft_odd) == 3,
710 "boundaryTensor3(): image for odd output must have 3 bands.");
712 detail::evenPolarFilters(supperleft, slowerright, sa,
713 dupperleft_even, even, scale,
false);
715 int w = slowerright.x - supperleft.x;
716 int h = slowerright.y - supperleft.y;
719 NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
720 typedef BasicImage<TinyVector<TmpType, 4> > TmpImage;
721 TmpImage t1(w, h), t2(w, h);
723 detail::KernelArray k1, k3;
724 detail::initGaussianPolarFilters1(scale, k1);
725 detail::initGaussianPolarFilters3(scale, k3);
728 VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t1.accessor());
730 destImage(t1, tmpBand), k1[3], k1[0]);
733 destImage(t1, tmpBand), k1[1], k1[2]);
736 destImage(t1, tmpBand), k3[3], k3[0]);
739 destImage(t1, tmpBand), k3[1], k3[2]);
743 destImage(t2, tmpBand), k1[0], k1[3]);
746 destImage(t2, tmpBand), k1[2], k1[1]);
749 destImage(t2, tmpBand), k3[0], k3[3]);
752 destImage(t2, tmpBand), k3[2], k3[1]);
755 typedef typename TmpImage::traverser TmpTraverser;
756 TmpTraverser tul1(t1.upperLeft());
757 TmpTraverser tlr1(t1.lowerRight());
758 TmpTraverser tul2(t2.upperLeft());
759 for(; tul1.y != tlr1.y; ++tul1.y, ++tul2.y, ++dupperleft_odd.y)
761 typename TmpTraverser::row_iterator tr1 = tul1.rowIterator();
762 typename TmpTraverser::row_iterator trend1 = tr1 + w;
763 typename TmpTraverser::row_iterator tr2 = tul2.rowIterator();
764 typename DestIteratorOdd::row_iterator o = dupperleft_odd.rowIterator();
765 for(; tr1 != trend1; ++tr1, ++tr2, ++o)
767 TmpType d11 = (*tr1)[0] + (*tr1)[2];
768 TmpType d12 = -(*tr1)[1] - (*tr1)[3];
769 TmpType d31 = (*tr2)[0] - (*tr2)[2];
770 TmpType d32 = (*tr2)[1] - (*tr2)[3];
771 TmpType d111 = 0.75 * d11 + 0.25 * d31;
772 TmpType d112 = 0.25 * (d12 + d32);
773 TmpType d122 = 0.25 * (d11 - d31);
774 TmpType d222 = 0.75 * d12 - 0.25 * d32;
775 TmpType d2 =
sq(d112);
776 TmpType d3 =
sq(d122);
778 odd.setComponent(0.25 * (
sq(d111) + 2.0*d2 + d3), o, 0);
779 odd.setComponent(0.25 * (d111*d112 + 2.0*d112*d122 + d122*d222), o, 1);
780 odd.setComponent(0.25 * (d2 + 2.0*d3 +
sq(d222)), o, 2);
785 template <
class SrcIterator,
class SrcAccessor,
786 class DestIteratorEven,
class DestAccessorEven,
787 class DestIteratorOdd,
class DestAccessorOdd>
789 void boundaryTensor3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
790 pair<DestIteratorEven, DestAccessorEven> even,
791 pair<DestIteratorOdd, DestAccessorOdd> odd,
794 boundaryTensor3(src.first, src.second, src.third,
795 even.first, even.second, odd.first, odd.second, scale);
802 #endif // VIGRA_BOUNDARYTENSOR_HXX
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
void rieszTransformOfLOG(...)
Calculate Riesz transforms of the Laplacian of Gaussian.
NumericTraits< T >::Promote sq(T t)
Definition: mathutil.hxx:341
void combineTwoImages(...)
Combine two source images into destination image.
void boundaryTensor1(...)
Boundary tensor variant.
void convolveImage(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor sa, DestIterator dupperleft, DestAccessor da, Kernel1D< T > const &kx, Kernel1D< T > const &ky)
Apply two separable filters successively, the first in x-direction, the second in y-direction...
Definition: convolution.hxx:257
void boundaryTensor(...)
Calculate the boundary tensor for a scalar valued image.
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616