37 #ifndef VIGRA_BOUNDARYTENSOR_HXX
38 #define VIGRA_BOUNDARYTENSOR_HXX
43 #include "array_vector.hxx"
44 #include "basicimage.hxx"
45 #include "combineimages.hxx"
46 #include "numerictraits.hxx"
47 #include "convolution.hxx"
48 #include "multi_shape.hxx"
56 typedef ArrayVector<Kernel1D<double> > KernelArray;
58 template <
class KernelArray>
60 initGaussianPolarFilters1(
double std_dev, KernelArray & k)
62 typedef typename KernelArray::value_type Kernel;
63 typedef typename Kernel::iterator iterator;
65 vigra_precondition(std_dev >= 0.0,
66 "initGaussianPolarFilter1(): "
67 "Standard deviation must be >= 0.");
71 int radius = (int)(4.0*std_dev + 0.5);
72 std_dev *= 1.08179074376;
73 double f = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / std_dev;
74 double a = 0.558868151788 / VIGRA_CSTD::pow(std_dev, 5);
75 double b = -2.04251639729 / VIGRA_CSTD::pow(std_dev, 3);
76 double sigma22 = -0.5 / std_dev / std_dev;
79 for(
unsigned int i=0; i<k.size(); ++i)
81 k[i].initExplicitly(-radius, radius);
82 k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
86 iterator c = k[0].center();
87 for(ix=-radius; ix<=radius; ++ix)
89 double x = (double)ix;
90 c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x);
94 for(ix=-radius; ix<=radius; ++ix)
96 double x = (double)ix;
97 c[ix] = f * x * VIGRA_CSTD::exp(sigma22 * x * x);
102 for(ix=-radius; ix<=radius; ++ix)
104 double x = (double)ix;
105 c[ix] = f * (b2 + a * x * x) * VIGRA_CSTD::exp(sigma22 * x * x);
109 for(ix=-radius; ix<=radius; ++ix)
111 double x = (double)ix;
112 c[ix] = f * x * (b + a * x * x) * VIGRA_CSTD::exp(sigma22 * x * x);
116 template <
class KernelArray>
118 initGaussianPolarFilters2(
double std_dev, KernelArray & k)
120 typedef typename KernelArray::value_type Kernel;
121 typedef typename Kernel::iterator iterator;
123 vigra_precondition(std_dev >= 0.0,
124 "initGaussianPolarFilter2(): "
125 "Standard deviation must be >= 0.");
129 int radius = (int)(4.0*std_dev + 0.5);
130 double f = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / std_dev;
131 double sigma2 = std_dev*std_dev;
132 double sigma22 = -0.5 / sigma2;
134 for(
unsigned int i=0; i<k.size(); ++i)
136 k[i].initExplicitly(-radius, radius);
137 k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
141 iterator c = k[0].center();
142 for(ix=-radius; ix<=radius; ++ix)
144 double x = (double)ix;
145 c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x);
149 double f1 = f / sigma2;
150 for(ix=-radius; ix<=radius; ++ix)
152 double x = (double)ix;
153 c[ix] = f1 * x * VIGRA_CSTD::exp(sigma22 * x * x);
157 double f2 = f / (sigma2 * sigma2);
158 for(ix=-radius; ix<=radius; ++ix)
160 double x = (double)ix;
161 c[ix] = f2 * (x * x - sigma2) * VIGRA_CSTD::exp(sigma22 * x * x);
165 template <
class KernelArray>
167 initGaussianPolarFilters3(
double std_dev, KernelArray & k)
169 typedef typename KernelArray::value_type Kernel;
170 typedef typename Kernel::iterator iterator;
172 vigra_precondition(std_dev >= 0.0,
173 "initGaussianPolarFilter3(): "
174 "Standard deviation must be >= 0.");
178 int radius = (int)(4.0*std_dev + 0.5);
179 std_dev *= 1.15470053838;
180 double sigma22 = -0.5 / std_dev / std_dev;
181 double f = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / std_dev;
182 double a = 0.883887052922 / VIGRA_CSTD::pow(std_dev, 5);
184 for(
unsigned int i=0; i<k.size(); ++i)
186 k[i].initExplicitly(-radius, radius);
187 k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
193 iterator c = k[0].center();
194 for(ix=-radius; ix<=radius; ++ix)
196 double x = (double)ix;
197 c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x);
201 for(ix=-radius; ix<=radius; ++ix)
203 double x = (double)ix;
204 c[ix] = f * x * VIGRA_CSTD::exp(sigma22 * x * x);
209 for(ix=-radius; ix<=radius; ++ix)
211 double x = (double)ix;
212 c[ix] = f * a2 * x * x * VIGRA_CSTD::exp(sigma22 * x * x);
216 for(ix=-radius; ix<=radius; ++ix)
218 double x = (double)ix;
219 c[ix] = f * a * x * x * x * VIGRA_CSTD::exp(sigma22 * x * x);
223 template <
class SrcIterator,
class SrcAccessor,
224 class DestIterator,
class DestAccessor>
226 evenPolarFilters(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
227 DestIterator dupperleft, DestAccessor dest,
228 double scale,
bool noLaplacian)
230 vigra_precondition(dest.size(dupperleft) == 3,
231 "evenPolarFilters(): image for even output must have 3 bands.");
233 int w = slowerright.x - supperleft.x;
234 int h = slowerright.y - supperleft.y;
237 NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
238 typedef BasicImage<TinyVector<TmpType, 3> > TmpImage;
239 typedef typename TmpImage::traverser TmpTraverser;
243 initGaussianPolarFilters2(scale, k2);
246 VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t.accessor());
248 destImage(t, tmpBand), k2[2], k2[0]);
251 destImage(t, tmpBand), k2[1], k2[1]);
254 destImage(t, tmpBand), k2[0], k2[2]);
257 TmpTraverser tul(t.upperLeft());
258 TmpTraverser tlr(t.lowerRight());
259 for(; tul.y != tlr.y; ++tul.y, ++dupperleft.y)
261 typename TmpTraverser::row_iterator tr = tul.rowIterator();
262 typename TmpTraverser::row_iterator trend = tr + w;
263 typename DestIterator::row_iterator d = dupperleft.rowIterator();
266 for(; tr != trend; ++tr, ++d)
268 TmpType v = detail::RequiresExplicitCast<TmpType>::cast(0.5*
sq((*tr)[0]-(*tr)[2]) + 2.0*
sq((*tr)[1]));
269 dest.setComponent(v, d, 0);
270 dest.setComponent(0, d, 1);
271 dest.setComponent(v, d, 2);
276 for(; tr != trend; ++tr, ++d)
278 dest.setComponent(
sq((*tr)[0]) +
sq((*tr)[1]), d, 0);
279 dest.setComponent(-(*tr)[1] * ((*tr)[0] + (*tr)[2]), d, 1);
280 dest.setComponent(
sq((*tr)[1]) +
sq((*tr)[2]), d, 2);
286 template <
class SrcIterator,
class SrcAccessor,
287 class DestIterator,
class DestAccessor>
289 oddPolarFilters(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
290 DestIterator dupperleft, DestAccessor dest,
291 double scale,
bool addResult)
293 vigra_precondition(dest.size(dupperleft) == 3,
294 "oddPolarFilters(): image for odd output must have 3 bands.");
296 int w = slowerright.x - supperleft.x;
297 int h = slowerright.y - supperleft.y;
300 NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
301 typedef BasicImage<TinyVector<TmpType, 4> > TmpImage;
302 typedef typename TmpImage::traverser TmpTraverser;
305 detail::KernelArray k1;
306 detail::initGaussianPolarFilters1(scale, k1);
309 VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t.accessor());
311 destImage(t, tmpBand), k1[3], k1[0]);
314 destImage(t, tmpBand), k1[2], k1[1]);
317 destImage(t, tmpBand), k1[1], k1[2]);
320 destImage(t, tmpBand), k1[0], k1[3]);
323 TmpTraverser tul(t.upperLeft());
324 TmpTraverser tlr(t.lowerRight());
325 for(; tul.y != tlr.y; ++tul.y, ++dupperleft.y)
327 typename TmpTraverser::row_iterator tr = tul.rowIterator();
328 typename TmpTraverser::row_iterator trend = tr + w;
329 typename DestIterator::row_iterator d = dupperleft.rowIterator();
332 for(; tr != trend; ++tr, ++d)
334 TmpType d0 = (*tr)[0] + (*tr)[2];
335 TmpType d1 = -(*tr)[1] - (*tr)[3];
337 dest.setComponent(dest.getComponent(d, 0) +
sq(d0), d, 0);
338 dest.setComponent(dest.getComponent(d, 1) + d0 * d1, d, 1);
339 dest.setComponent(dest.getComponent(d, 2) +
sq(d1), d, 2);
344 for(; tr != trend; ++tr, ++d)
346 TmpType d0 = (*tr)[0] + (*tr)[2];
347 TmpType d1 = -(*tr)[1] - (*tr)[3];
349 dest.setComponent(
sq(d0), d, 0);
350 dest.setComponent(d0 * d1, d, 1);
351 dest.setComponent(
sq(d1), d, 2);
436 template <
class SrcIterator,
class SrcAccessor,
437 class DestIterator,
class DestAccessor>
439 DestIterator dupperleft, DestAccessor dest,
440 double scale,
unsigned int xorder,
unsigned int yorder)
442 unsigned int order = xorder + yorder;
444 vigra_precondition(order <= 2,
445 "rieszTransformOfLOG(): can only compute Riesz transforms up to order 2.");
446 vigra_precondition(scale > 0.0,
447 "rieszTransformOfLOG(): scale must be positive.");
449 int w = slowerright.x - supperleft.x;
450 int h = slowerright.y - supperleft.y;
452 typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
460 detail::initGaussianPolarFilters2(scale, k2);
462 TmpImage t1(w, h), t2(w, h);
465 destImage(t1), k2[2], k2[0]);
467 destImage(t2), k2[0], k2[2]);
469 destIter(dupperleft, dest), std::plus<TmpType>());
475 detail::initGaussianPolarFilters1(scale, k1);
477 TmpImage t1(w, h), t2(w, h);
482 destImage(t1), k1[3], k1[0]);
484 destImage(t2), k1[1], k1[2]);
489 destImage(t1), k1[0], k1[3]);
491 destImage(t2), k1[2], k1[1]);
494 destIter(dupperleft, dest), std::plus<TmpType>());
500 detail::initGaussianPolarFilters2(scale, k2);
503 destIter(dupperleft, dest), k2[xorder], k2[yorder]);
510 detail::initGaussianPolarFilters3(scale, k3);
512 TmpImage t1(w, h), t2(w, h);
517 destImage(t1), k3[3], k3[0]);
519 destImage(t2), k3[1], k3[2]);
524 destImage(t1), k3[0], k3[3]);
526 destImage(t2), k3[2], k3[1]);
529 destIter(dupperleft, dest), std::minus<TmpType>());
535 template <
class SrcIterator,
class SrcAccessor,
536 class DestIterator,
class DestAccessor>
539 pair<DestIterator, DestAccessor> dest,
540 double scale,
unsigned int xorder,
unsigned int yorder)
543 scale, xorder, yorder);
546 template <
class T1,
class S1,
550 MultiArrayView<2, T2, S2> dest,
551 double scale,
unsigned int xorder,
unsigned int yorder)
553 vigra_precondition(src.shape() == dest.shape(),
554 "rieszTransformOfLOG(): shape mismatch between input and output.");
556 scale, xorder, yorder);
639 template <
class SrcIterator,
class SrcAccessor,
640 class DestIterator,
class DestAccessor>
641 void boundaryTensor(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
642 DestIterator dupperleft, DestAccessor dest,
645 vigra_precondition(dest.size(dupperleft) == 3,
646 "boundaryTensor(): image for even output must have 3 bands.");
647 vigra_precondition(scale > 0.0,
648 "boundaryTensor(): scale must be positive.");
650 detail::evenPolarFilters(supperleft, slowerright, src,
651 dupperleft, dest, scale,
false);
652 detail::oddPolarFilters(supperleft, slowerright, src,
653 dupperleft, dest, scale,
true);
656 template <
class SrcIterator,
class SrcAccessor,
657 class DestIterator,
class DestAccessor>
659 void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src,
660 pair<DestIterator, DestAccessor> dest,
664 dest.first, dest.second, scale);
667 template <
class T1,
class S1,
671 MultiArrayView<2, T2, S2> dest,
674 vigra_precondition(src.shape() == dest.shape(),
675 "boundaryTensor(): shape mismatch between input and output.");
677 destImage(dest), scale);
729 template <
class SrcIterator,
class SrcAccessor,
730 class DestIterator,
class DestAccessor>
731 void boundaryTensor1(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
732 DestIterator dupperleft, DestAccessor dest,
735 vigra_precondition(dest.size(dupperleft) == 3,
736 "boundaryTensor1(): image for even output must have 3 bands.");
737 vigra_precondition(scale > 0.0,
738 "boundaryTensor1(): scale must be positive.");
740 detail::evenPolarFilters(supperleft, slowerright, src,
741 dupperleft, dest, scale,
true);
742 detail::oddPolarFilters(supperleft, slowerright, src,
743 dupperleft, dest, scale,
true);
746 template <
class SrcIterator,
class SrcAccessor,
747 class DestIterator,
class DestAccessor>
749 void boundaryTensor1(triple<SrcIterator, SrcIterator, SrcAccessor> src,
750 pair<DestIterator, DestAccessor> dest,
754 dest.first, dest.second, scale);
757 template <
class T1,
class S1,
761 MultiArrayView<2, T2, S2> dest,
764 vigra_precondition(src.shape() == dest.shape(),
765 "boundaryTensor1(): shape mismatch between input and output.");
767 destImage(dest), scale);
779 template <
class SrcIterator,
class SrcAccessor,
780 class DestIteratorEven,
class DestAccessorEven,
781 class DestIteratorOdd,
class DestAccessorOdd>
782 void boundaryTensor3(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor sa,
783 DestIteratorEven dupperleft_even, DestAccessorEven
even,
784 DestIteratorOdd dupperleft_odd, DestAccessorOdd
odd,
787 vigra_precondition(
even.size(dupperleft_even) == 3,
788 "boundaryTensor3(): image for even output must have 3 bands.");
789 vigra_precondition(
odd.size(dupperleft_odd) == 3,
790 "boundaryTensor3(): image for odd output must have 3 bands.");
792 detail::evenPolarFilters(supperleft, slowerright, sa,
793 dupperleft_even,
even, scale,
false);
795 int w = slowerright.x - supperleft.x;
796 int h = slowerright.y - supperleft.y;
799 NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
800 typedef BasicImage<TinyVector<TmpType, 4> > TmpImage;
801 TmpImage t1(w, h), t2(w, h);
803 detail::KernelArray k1, k3;
804 detail::initGaussianPolarFilters1(scale, k1);
805 detail::initGaussianPolarFilters3(scale, k3);
808 VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t1.accessor());
810 destImage(t1, tmpBand), k1[3], k1[0]);
813 destImage(t1, tmpBand), k1[1], k1[2]);
816 destImage(t1, tmpBand), k3[3], k3[0]);
819 destImage(t1, tmpBand), k3[1], k3[2]);
823 destImage(t2, tmpBand), k1[0], k1[3]);
826 destImage(t2, tmpBand), k1[2], k1[1]);
829 destImage(t2, tmpBand), k3[0], k3[3]);
832 destImage(t2, tmpBand), k3[2], k3[1]);
835 typedef typename TmpImage::traverser TmpTraverser;
836 TmpTraverser tul1(t1.upperLeft());
837 TmpTraverser tlr1(t1.lowerRight());
838 TmpTraverser tul2(t2.upperLeft());
839 for(; tul1.y != tlr1.y; ++tul1.y, ++tul2.y, ++dupperleft_odd.y)
841 typename TmpTraverser::row_iterator tr1 = tul1.rowIterator();
842 typename TmpTraverser::row_iterator trend1 = tr1 + w;
843 typename TmpTraverser::row_iterator tr2 = tul2.rowIterator();
844 typename DestIteratorOdd::row_iterator o = dupperleft_odd.rowIterator();
845 for(; tr1 != trend1; ++tr1, ++tr2, ++o)
847 TmpType d11 = (*tr1)[0] + (*tr1)[2];
848 TmpType d12 = -(*tr1)[1] - (*tr1)[3];
849 TmpType d31 = (*tr2)[0] - (*tr2)[2];
850 TmpType d32 = (*tr2)[1] - (*tr2)[3];
851 TmpType d111 = 0.75 * d11 + 0.25 * d31;
852 TmpType d112 = 0.25 * (d12 + d32);
853 TmpType d122 = 0.25 * (d11 - d31);
854 TmpType d222 = 0.75 * d12 - 0.25 * d32;
855 TmpType d2 =
sq(d112);
856 TmpType d3 =
sq(d122);
858 odd.setComponent(0.25 * (
sq(d111) + 2.0*d2 + d3), o, 0);
859 odd.setComponent(0.25 * (d111*d112 + 2.0*d112*d122 + d122*d222), o, 1);
860 odd.setComponent(0.25 * (d2 + 2.0*d3 +
sq(d222)), o, 2);
865 template <
class SrcIterator,
class SrcAccessor,
866 class DestIteratorEven,
class DestAccessorEven,
867 class DestIteratorOdd,
class DestAccessorOdd>
869 void boundaryTensor3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
870 pair<DestIteratorEven, DestAccessorEven>
even,
871 pair<DestIteratorOdd, DestAccessorOdd>
odd,
874 boundaryTensor3(src.first, src.second, src.third,
878 template <
class T1,
class S1,
879 class T2E,
class S2Even,
880 class T2O,
class S2Odd>
882 void boundaryTensor3(MultiArrayView<2, T1, S1>
const & src,
883 MultiArrayView<2, T2E, S2Even>
even,
884 MultiArrayView<2, T2O, S2Odd>
odd,
887 vigra_precondition(src.shape() ==
even.shape() && src.shape() ==
odd.shape(),
888 "boundaryTensor3(): shape mismatch between input and output.");
889 boundaryTensor3(srcImageRange(src),
890 destImage(
even), destImage(
odd), scale);
Definition: array_vector.hxx:514
Fundamental class template for images.
Definition: basicimage.hxx:476
void boundaryTensor1(...)
Boundary tensor variant.
void rieszTransformOfLOG(...)
Calculate Riesz transforms of the Laplacian of Gaussian.
bool even(int t)
Check if an integer is even.
bool odd(int t)
Check if an integer is odd.
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:382
void convolveImage(...)
Convolve an image with the given kernel(s).
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
void boundaryTensor(...)
Calculate the boundary tensor for a scalar valued image.
void combineTwoImages(...)
Combine two source images into destination image.