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"
56typedef ArrayVector<Kernel1D<double> > KernelArray;
58template <
class KernelArray>
60initGaussianPolarFilters1(
double std_dev, KernelArray & k)
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);
116template <
class KernelArray>
118initGaussianPolarFilters2(
double std_dev, KernelArray & k)
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);
165template <
class KernelArray>
167initGaussianPolarFilters3(
double std_dev, KernelArray & k)
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);
223template <
class SrcIterator,
class SrcAccessor,
224 class DestIterator,
class DestAccessor>
226evenPolarFilters(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);
286template <
class SrcIterator,
class SrcAccessor,
287 class DestIterator,
class DestAccessor>
289oddPolarFilters(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);
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.");
452 typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote
TmpType;
460 detail::initGaussianPolarFilters2(scale,
k2);
465 destImage(
t1),
k2[2],
k2[0]);
467 destImage(
t2),
k2[0],
k2[2]);
475 detail::initGaussianPolarFilters1(scale,
k1);
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]);
500 detail::initGaussianPolarFilters2(scale,
k2);
510 detail::initGaussianPolarFilters3(scale,
k3);
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]);
535template <
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);
546template <
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);
646 "boundaryTensor(): image for even output must have 3 bands.");
647 vigra_precondition(scale > 0.0,
648 "boundaryTensor(): scale must be positive.");
656template <
class SrcIterator,
class SrcAccessor,
657 class DestIterator,
class DestAccessor>
659void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src,
660 pair<DestIterator, DestAccessor> dest,
664 dest.first, dest.second, scale);
667template <
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);
736 "boundaryTensor1(): image for even output must have 3 bands.");
737 vigra_precondition(scale > 0.0,
738 "boundaryTensor1(): scale must be positive.");
746template <
class SrcIterator,
class SrcAccessor,
747 class DestIterator,
class DestAccessor>
750 pair<DestIterator, DestAccessor> dest,
754 dest.first, dest.second, scale);
757template <
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);
779template <
class SrcIterator,
class SrcAccessor,
780 class DestIteratorEven,
class DestAccessorEven,
781 class DestIteratorOdd,
class DestAccessorOdd>
782void 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);
865template <
class SrcIterator,
class SrcAccessor,
866 class DestIteratorEven,
class DestAccessorEven,
867 class DestIteratorOdd,
class DestAccessorOdd>
869void boundaryTensor3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
870 pair<DestIteratorEven, DestAccessorEven>
even,
871 pair<DestIteratorOdd, DestAccessorOdd>
odd,
874 boundaryTensor3(src.first, src.second, src.third,
878template <
class T1,
class S1,
879 class T2E,
class S2Even,
880 class T2O,
class S2Odd>
882void 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);
Class for a single RGB value.
Definition rgbvalue.hxx:128
Base::value_type value_type
Definition rgbvalue.hxx:141
size_type size() const
Definition tinyvector.hxx:913
bool even(int t)
Check if an integer is even.
void combineTwoImages(...)
Combine two source images into destination image.
void boundaryTensor(...)
Calculate the boundary tensor for a scalar valued image.
void rieszTransformOfLOG(...)
Calculate Riesz transforms of the Laplacian of Gaussian.
void convolveImage(...)
Convolve an image with the given kernel(s).
void boundaryTensor1(...)
Boundary tensor variant.
bool odd(int t)
Check if an integer is odd.