[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

boundarytensor.hxx
1/************************************************************************/
2/* */
3/* Copyright 1998-2004 by Ullrich Koethe */
4/* */
5/* This file is part of the VIGRA computer vision library. */
6/* The VIGRA Website is */
7/* http://hci.iwr.uni-heidelberg.de/vigra/ */
8/* Please direct questions, bug reports, and contributions to */
9/* ullrich.koethe@iwr.uni-heidelberg.de or */
10/* vigra@informatik.uni-hamburg.de */
11/* */
12/* Permission is hereby granted, free of charge, to any person */
13/* obtaining a copy of this software and associated documentation */
14/* files (the "Software"), to deal in the Software without */
15/* restriction, including without limitation the rights to use, */
16/* copy, modify, merge, publish, distribute, sublicense, and/or */
17/* sell copies of the Software, and to permit persons to whom the */
18/* Software is furnished to do so, subject to the following */
19/* conditions: */
20/* */
21/* The above copyright notice and this permission notice shall be */
22/* included in all copies or substantial portions of the */
23/* Software. */
24/* */
25/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32/* OTHER DEALINGS IN THE SOFTWARE. */
33/* */
34/************************************************************************/
35
36
37#ifndef VIGRA_BOUNDARYTENSOR_HXX
38#define VIGRA_BOUNDARYTENSOR_HXX
39
40#include <cmath>
41#include <functional>
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"
48#include "multi_shape.hxx"
49
50namespace vigra {
51
52namespace detail {
53
54/***********************************************************************/
55
56typedef ArrayVector<Kernel1D<double> > KernelArray;
57
58template <class KernelArray>
59void
60initGaussianPolarFilters1(double std_dev, KernelArray & k)
61{
62 typedef typename KernelArray::value_type Kernel;
63 typedef typename Kernel::iterator iterator;
64
65 vigra_precondition(std_dev >= 0.0,
66 "initGaussianPolarFilter1(): "
67 "Standard deviation must be >= 0.");
68
69 k.resize(4);
70
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; // norm
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;
77
78
79 for(unsigned int i=0; i<k.size(); ++i)
80 {
81 k[i].initExplicitly(-radius, radius);
82 k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
83 }
84
85 int ix;
86 iterator c = k[0].center();
87 for(ix=-radius; ix<=radius; ++ix)
88 {
89 double x = (double)ix;
90 c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x);
91 }
92
93 c = k[1].center();
94 for(ix=-radius; ix<=radius; ++ix)
95 {
96 double x = (double)ix;
97 c[ix] = f * x * VIGRA_CSTD::exp(sigma22 * x * x);
98 }
99
100 c = k[2].center();
101 double b2 = b / 3.0;
102 for(ix=-radius; ix<=radius; ++ix)
103 {
104 double x = (double)ix;
105 c[ix] = f * (b2 + a * x * x) * VIGRA_CSTD::exp(sigma22 * x * x);
106 }
107
108 c = k[3].center();
109 for(ix=-radius; ix<=radius; ++ix)
110 {
111 double x = (double)ix;
112 c[ix] = f * x * (b + a * x * x) * VIGRA_CSTD::exp(sigma22 * x * x);
113 }
114}
115
116template <class KernelArray>
117void
118initGaussianPolarFilters2(double std_dev, KernelArray & k)
119{
120 typedef typename KernelArray::value_type Kernel;
121 typedef typename Kernel::iterator iterator;
122
123 vigra_precondition(std_dev >= 0.0,
124 "initGaussianPolarFilter2(): "
125 "Standard deviation must be >= 0.");
126
127 k.resize(3);
128
129 int radius = (int)(4.0*std_dev + 0.5);
130 double f = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / std_dev; // norm
131 double sigma2 = std_dev*std_dev;
132 double sigma22 = -0.5 / sigma2;
133
134 for(unsigned int i=0; i<k.size(); ++i)
135 {
136 k[i].initExplicitly(-radius, radius);
137 k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
138 }
139
140 int ix;
141 iterator c = k[0].center();
142 for(ix=-radius; ix<=radius; ++ix)
143 {
144 double x = (double)ix;
145 c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x);
146 }
147
148 c = k[1].center();
149 double f1 = f / sigma2;
150 for(ix=-radius; ix<=radius; ++ix)
151 {
152 double x = (double)ix;
153 c[ix] = f1 * x * VIGRA_CSTD::exp(sigma22 * x * x);
154 }
155
156 c = k[2].center();
157 double f2 = f / (sigma2 * sigma2);
158 for(ix=-radius; ix<=radius; ++ix)
159 {
160 double x = (double)ix;
161 c[ix] = f2 * (x * x - sigma2) * VIGRA_CSTD::exp(sigma22 * x * x);
162 }
163}
164
165template <class KernelArray>
166void
167initGaussianPolarFilters3(double std_dev, KernelArray & k)
168{
169 typedef typename KernelArray::value_type Kernel;
170 typedef typename Kernel::iterator iterator;
171
172 vigra_precondition(std_dev >= 0.0,
173 "initGaussianPolarFilter3(): "
174 "Standard deviation must be >= 0.");
175
176 k.resize(4);
177
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; // norm
182 double a = 0.883887052922 / VIGRA_CSTD::pow(std_dev, 5);
183
184 for(unsigned int i=0; i<k.size(); ++i)
185 {
186 k[i].initExplicitly(-radius, radius);
187 k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
188 }
189
190 //double b = -1.3786348292 / VIGRA_CSTD::pow(std_dev, 3);
191
192 int ix;
193 iterator c = k[0].center();
194 for(ix=-radius; ix<=radius; ++ix)
195 {
196 double x = (double)ix;
197 c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x);
198 }
199
200 c = k[1].center();
201 for(ix=-radius; ix<=radius; ++ix)
202 {
203 double x = (double)ix;
204 c[ix] = f * x * VIGRA_CSTD::exp(sigma22 * x * x);
205 }
206
207 c = k[2].center();
208 double a2 = 3.0 * a;
209 for(ix=-radius; ix<=radius; ++ix)
210 {
211 double x = (double)ix;
212 c[ix] = f * a2 * x * x * VIGRA_CSTD::exp(sigma22 * x * x);
213 }
214
215 c = k[3].center();
216 for(ix=-radius; ix<=radius; ++ix)
217 {
218 double x = (double)ix;
219 c[ix] = f * a * x * x * x * VIGRA_CSTD::exp(sigma22 * x * x);
220 }
221}
222
223template <class SrcIterator, class SrcAccessor,
224 class DestIterator, class DestAccessor>
225void
226evenPolarFilters(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
227 DestIterator dupperleft, DestAccessor dest,
228 double scale, bool noLaplacian)
229{
230 vigra_precondition(dest.size(dupperleft) == 3,
231 "evenPolarFilters(): image for even output must have 3 bands.");
232
233 int w = slowerright.x - supperleft.x;
234 int h = slowerright.y - supperleft.y;
235
236 typedef typename
237 NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
238 typedef BasicImage<TinyVector<TmpType, 3> > TmpImage;
239 typedef typename TmpImage::traverser TmpTraverser;
240 TmpImage t(w, h);
241
242 KernelArray k2;
243 initGaussianPolarFilters2(scale, k2);
244
245 // calculate filter responses for even filters
246 VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t.accessor());
247 convolveImage(srcIterRange(supperleft, slowerright, src),
248 destImage(t, tmpBand), k2[2], k2[0]);
249 tmpBand.setIndex(1);
250 convolveImage(srcIterRange(supperleft, slowerright, src),
251 destImage(t, tmpBand), k2[1], k2[1]);
252 tmpBand.setIndex(2);
253 convolveImage(srcIterRange(supperleft, slowerright, src),
254 destImage(t, tmpBand), k2[0], k2[2]);
255
256 // create even tensor from filter responses
257 TmpTraverser tul(t.upperLeft());
258 TmpTraverser tlr(t.lowerRight());
259 for(; tul.y != tlr.y; ++tul.y, ++dupperleft.y)
260 {
261 typename TmpTraverser::row_iterator tr = tul.rowIterator();
262 typename TmpTraverser::row_iterator trend = tr + w;
263 typename DestIterator::row_iterator d = dupperleft.rowIterator();
264 if(noLaplacian)
265 {
266 for(; tr != trend; ++tr, ++d)
267 {
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);
272 }
273 }
274 else
275 {
276 for(; tr != trend; ++tr, ++d)
277 {
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);
281 }
282 }
283 }
284}
285
286template <class SrcIterator, class SrcAccessor,
287 class DestIterator, class DestAccessor>
288void
289oddPolarFilters(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
290 DestIterator dupperleft, DestAccessor dest,
291 double scale, bool addResult)
292{
293 vigra_precondition(dest.size(dupperleft) == 3,
294 "oddPolarFilters(): image for odd output must have 3 bands.");
295
296 int w = slowerright.x - supperleft.x;
297 int h = slowerright.y - supperleft.y;
298
299 typedef typename
300 NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
301 typedef BasicImage<TinyVector<TmpType, 4> > TmpImage;
302 typedef typename TmpImage::traverser TmpTraverser;
303 TmpImage t(w, h);
304
305 detail::KernelArray k1;
306 detail::initGaussianPolarFilters1(scale, k1);
307
308 // calculate filter responses for odd filters
309 VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t.accessor());
310 convolveImage(srcIterRange(supperleft, slowerright, src),
311 destImage(t, tmpBand), k1[3], k1[0]);
312 tmpBand.setIndex(1);
313 convolveImage(srcIterRange(supperleft, slowerright, src),
314 destImage(t, tmpBand), k1[2], k1[1]);
315 tmpBand.setIndex(2);
316 convolveImage(srcIterRange(supperleft, slowerright, src),
317 destImage(t, tmpBand), k1[1], k1[2]);
318 tmpBand.setIndex(3);
319 convolveImage(srcIterRange(supperleft, slowerright, src),
320 destImage(t, tmpBand), k1[0], k1[3]);
321
322 // create odd tensor from filter responses
323 TmpTraverser tul(t.upperLeft());
324 TmpTraverser tlr(t.lowerRight());
325 for(; tul.y != tlr.y; ++tul.y, ++dupperleft.y)
326 {
327 typename TmpTraverser::row_iterator tr = tul.rowIterator();
328 typename TmpTraverser::row_iterator trend = tr + w;
329 typename DestIterator::row_iterator d = dupperleft.rowIterator();
330 if(addResult)
331 {
332 for(; tr != trend; ++tr, ++d)
333 {
334 TmpType d0 = (*tr)[0] + (*tr)[2];
335 TmpType d1 = -(*tr)[1] - (*tr)[3];
336
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);
340 }
341 }
342 else
343 {
344 for(; tr != trend; ++tr, ++d)
345 {
346 TmpType d0 = (*tr)[0] + (*tr)[2];
347 TmpType d1 = -(*tr)[1] - (*tr)[3];
348
349 dest.setComponent(sq(d0), d, 0);
350 dest.setComponent(d0 * d1, d, 1);
351 dest.setComponent(sq(d1), d, 2);
352 }
353 }
354 }
355}
356
357} // namespace detail
358
359/** \addtogroup ConvolutionFilters
360*/
361//@{
362
363/********************************************************/
364/* */
365/* rieszTransformOfLOG */
366/* */
367/********************************************************/
368
369/** \brief Calculate Riesz transforms of the Laplacian of Gaussian.
370
371 The Riesz transforms of the Laplacian of Gaussian have the following transfer
372 functions (defined in a polar coordinate representation of the frequency domain):
373
374 \f[
375 F_{\sigma}(r, \phi)=(i \cos \phi)^n (i \sin \phi)^m r^2 e^{-r^2 \sigma^2 / 2}
376 \f]
377
378 where <i>n</i> = <tt>xorder</tt> and <i>m</i> = <tt>yorder</tt> determine th e
379 order of the transform, and <tt>sigma > 0</tt> is the scale of the Laplacian
380 of Gaussian. This function computes a good spatial domain approximation of
381 these transforms for <tt>xorder + yorder <= 2</tt>. The filter responses may be used
382 to calculate the monogenic signal or the boundary tensor.
383
384 <b> Declarations:</b>
385
386 pass 2D array views:
387 \code
388 namespace vigra {
389 template <class T1, class S1,
390 class T2, class S2>
391 void
392 rieszTransformOfLOG(MultiArrayView<2, T1, S1> const & src,
393 MultiArrayView<2, T2, S2> dest,
394 double scale, unsigned int xorder, unsigned int yorder);
395 }
396 \endcode
397
398 \deprecatedAPI{rieszTransformOfLOG}
399 pass \ref ImageIterators and \ref DataAccessors :
400 \code
401 namespace vigra {
402 template <class SrcIterator, class SrcAccessor,
403 class DestIterator, class DestAccessor>
404 void rieszTransformOfLOG(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
405 DestIterator dupperleft, DestAccessor dest,
406 double scale, unsigned int xorder, unsigned int yorder);
407 }
408 \endcode
409 use argument objects in conjunction with \ref ArgumentObjectFactories :
410 \code
411 namespace vigra {
412 template <class SrcIterator, class SrcAccessor,
413 class DestIterator, class DestAccessor>
414 void rieszTransformOfLOG(triple<SrcIterator, SrcIterator, SrcAccessor> src,
415 pair<DestIterator, DestAccessor> dest,
416 double scale, unsigned int xorder, unsigned int yorder);
417 }
418 \endcode
419 \deprecatedEnd
420
421 <b> Usage:</b>
422
423 <b>\#include</b> <vigra/boundarytensor.hxx><br>
424 Namespace: vigra
425
426 \code
427 MultiArrayView<2, double> impulse(17,17), res(17, 17);
428 impulse(8,8) = 1.0;
429
430 // calculate the impulse response of the first order Riesz transform in x-direction
431 rieszTransformOfLOG(impulse, res, 2.0, 1, 0);
432 \endcode
433*/
434doxygen_overloaded_function(template <...> void rieszTransformOfLOG)
435
436template <class SrcIterator, class SrcAccessor,
437 class DestIterator, class DestAccessor>
438void rieszTransformOfLOG(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
439 DestIterator dupperleft, DestAccessor dest,
440 double scale, unsigned int xorder, unsigned int yorder)
441{
442 unsigned int order = xorder + yorder;
443
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.");
448
449 int w = slowerright.x - supperleft.x;
450 int h = slowerright.y - supperleft.y;
451
452 typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
453 typedef BasicImage<TmpType> TmpImage;
454
455 switch(order)
456 {
457 case 0:
458 {
459 detail::KernelArray k2;
460 detail::initGaussianPolarFilters2(scale, k2);
461
462 TmpImage t1(w, h), t2(w, h);
463
464 convolveImage(srcIterRange(supperleft, slowerright, src),
465 destImage(t1), k2[2], k2[0]);
466 convolveImage(srcIterRange(supperleft, slowerright, src),
467 destImage(t2), k2[0], k2[2]);
468 combineTwoImages(srcImageRange(t1), srcImage(t2),
469 destIter(dupperleft, dest), std::plus<TmpType>());
470 break;
471 }
472 case 1:
473 {
474 detail::KernelArray k1;
475 detail::initGaussianPolarFilters1(scale, k1);
476
477 TmpImage t1(w, h), t2(w, h);
478
479 if(xorder == 1)
480 {
481 convolveImage(srcIterRange(supperleft, slowerright, src),
482 destImage(t1), k1[3], k1[0]);
483 convolveImage(srcIterRange(supperleft, slowerright, src),
484 destImage(t2), k1[1], k1[2]);
485 }
486 else
487 {
488 convolveImage(srcIterRange(supperleft, slowerright, src),
489 destImage(t1), k1[0], k1[3]);
490 convolveImage(srcIterRange(supperleft, slowerright, src),
491 destImage(t2), k1[2], k1[1]);
492 }
493 combineTwoImages(srcImageRange(t1), srcImage(t2),
494 destIter(dupperleft, dest), std::plus<TmpType>());
495 break;
496 }
497 case 2:
498 {
499 detail::KernelArray k2;
500 detail::initGaussianPolarFilters2(scale, k2);
501
502 convolveImage(srcIterRange(supperleft, slowerright, src),
503 destIter(dupperleft, dest), k2[xorder], k2[yorder]);
504 break;
505 }
506 /* for test purposes only: compute 3rd order polar filters */
507 case 3:
508 {
509 detail::KernelArray k3;
510 detail::initGaussianPolarFilters3(scale, k3);
511
512 TmpImage t1(w, h), t2(w, h);
513
514 if(xorder == 3)
515 {
516 convolveImage(srcIterRange(supperleft, slowerright, src),
517 destImage(t1), k3[3], k3[0]);
518 convolveImage(srcIterRange(supperleft, slowerright, src),
519 destImage(t2), k3[1], k3[2]);
520 }
521 else
522 {
523 convolveImage(srcIterRange(supperleft, slowerright, src),
524 destImage(t1), k3[0], k3[3]);
525 convolveImage(srcIterRange(supperleft, slowerright, src),
526 destImage(t2), k3[2], k3[1]);
527 }
528 combineTwoImages(srcImageRange(t1), srcImage(t2),
529 destIter(dupperleft, dest), std::minus<TmpType>());
530 break;
531 }
532 }
533}
534
535template <class SrcIterator, class SrcAccessor,
536 class DestIterator, class DestAccessor>
537inline
538void rieszTransformOfLOG(triple<SrcIterator, SrcIterator, SrcAccessor> src,
539 pair<DestIterator, DestAccessor> dest,
540 double scale, unsigned int xorder, unsigned int yorder)
541{
542 rieszTransformOfLOG(src.first, src.second, src.third, dest.first, dest.second,
543 scale, xorder, yorder);
544}
545
546template <class T1, class S1,
547 class T2, class S2>
548inline void
551 double scale, unsigned int xorder, unsigned int yorder)
552{
553 vigra_precondition(src.shape() == dest.shape(),
554 "rieszTransformOfLOG(): shape mismatch between input and output.");
555 rieszTransformOfLOG(srcImageRange(src), destImage(dest),
556 scale, xorder, yorder);
557}
558
559//@}
560
561/** \addtogroup TensorImaging Tensor Image Processing
562*/
563//@{
564
565/********************************************************/
566/* */
567/* boundaryTensor */
568/* */
569/********************************************************/
570
571/** \brief Calculate the boundary tensor for a scalar valued image.
572
573 These functions calculate a spatial domain approximation of
574 the boundary tensor as described in
575
576 U. K&ouml;the: <a href="http://hci.iwr.uni-heidelberg.de/people/ukoethe/papers/index.php#cite_polarfilters">
577 <i>"Integrated Edge and Junction Detection with the Boundary Tensor"</i></a>,
578 in: ICCV 03, Proc. of 9th Intl. Conf. on Computer Vision, Nice 2003, vol. 1,
579 pp. 424-431, Los Alamitos: IEEE Computer Society, 2003
580
581 with the Laplacian of Gaussian as the underlying bandpass filter (see
582 \ref rieszTransformOfLOG()). The output image must have 3 bands which will hold the
583 tensor components in the order t11, t12 (== t21), t22. The function
584 \ref boundaryTensor1() with the same interface implements a variant of the
585 boundary tensor where the 0th-order Riesz transform has been dropped, so that the
586 tensor is no longer sensitive to blobs.
587
588 <b> Declarations:</b>
589
590 pass 2D array views:
591 \code
592 namespace vigra {
593 template <class T1, class S1,
594 class T2, class S2>
595 void
596 boundaryTensor(MultiArrayView<2, T1, S1> const & src,
597 MultiArrayView<2, T2, S2> dest,
598 double scale);
599 }
600 \endcode
601
602 \deprecatedAPI{boundaryTensor}
603 pass \ref ImageIterators and \ref DataAccessors :
604 \code
605 namespace vigra {
606 template <class SrcIterator, class SrcAccessor,
607 class DestIterator, class DestAccessor>
608 void boundaryTensor(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
609 DestIterator dupperleft, DestAccessor dest,
610 double scale);
611 }
612 \endcode
613 use argument objects in conjunction with \ref ArgumentObjectFactories :
614 \code
615 namespace vigra {
616 template <class SrcIterator, class SrcAccessor,
617 class DestIterator, class DestAccessor>
618 void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src,
619 pair<DestIterator, DestAccessor> dest,
620 double scale);
621 }
622 \endcode
623 \deprecatedEnd
624
625 <b> Usage:</b>
626
627 <b>\#include</b> <vigra/boundarytensor.hxx><br/>
628 Namespace: vigra
629
630 \code
631 MultiArray<2, float> img(w,h);
632 MultiArray<2, TinyVector<float, 3> bt(w,h);
633 ...
634 boundaryTensor(img, bt, 2.0);
635 \endcode
636*/
637doxygen_overloaded_function(template <...> void boundaryTensor)
638
639template <class SrcIterator, class SrcAccessor,
640 class DestIterator, class DestAccessor>
641void boundaryTensor(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
642 DestIterator dupperleft, DestAccessor dest,
643 double scale)
644{
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.");
649
650 detail::evenPolarFilters(supperleft, slowerright, src,
651 dupperleft, dest, scale, false);
652 detail::oddPolarFilters(supperleft, slowerright, src,
653 dupperleft, dest, scale, true);
654}
655
656template <class SrcIterator, class SrcAccessor,
657 class DestIterator, class DestAccessor>
658inline
659void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src,
660 pair<DestIterator, DestAccessor> dest,
661 double scale)
662{
663 boundaryTensor(src.first, src.second, src.third,
664 dest.first, dest.second, scale);
665}
666
667template <class T1, class S1,
668 class T2, class S2>
669inline void
672 double scale)
673{
674 vigra_precondition(src.shape() == dest.shape(),
675 "boundaryTensor(): shape mismatch between input and output.");
676 boundaryTensor(srcImageRange(src),
677 destImage(dest), scale);
678}
679
680/** \brief Boundary tensor variant.
681
682 This function implements a variant of the boundary tensor where the
683 0th-order Riesz transform has been dropped, so that the tensor is no
684 longer sensitive to blobs. See \ref boundaryTensor() for more detailed
685 documentation.
686
687 <b> Declarations:</b>
688
689 <b>\#include</b> <vigra/boundarytensor.hxx><br/>
690 Namespace: vigra
691
692 pass 2D array views:
693 \code
694 namespace vigra {
695 template <class T1, class S1,
696 class T2, class S2>
697 void
698 boundaryTensor1(MultiArrayView<2, T1, S1> const & src,
699 MultiArrayView<2, T2, S2> dest,
700 double scale);
701 }
702 \endcode
703
704 \deprecatedAPI{boundaryTensor1}
705 pass \ref ImageIterators and \ref DataAccessors :
706 \code
707 namespace vigra {
708 template <class SrcIterator, class SrcAccessor,
709 class DestIterator, class DestAccessor>
710 void boundaryTensor1(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
711 DestIterator dupperleft, DestAccessor dest,
712 double scale);
713 }
714 \endcode
715 use argument objects in conjunction with \ref ArgumentObjectFactories :
716 \code
717 namespace vigra {
718 template <class SrcIterator, class SrcAccessor,
719 class DestIterator, class DestAccessor>
720 void boundaryTensor1(triple<SrcIterator, SrcIterator, SrcAccessor> src,
721 pair<DestIterator, DestAccessor> dest,
722 double scale);
723 }
724 \endcode
725 \deprecatedEnd
726*/
727doxygen_overloaded_function(template <...> void boundaryTensor1)
728
729template <class SrcIterator, class SrcAccessor,
730 class DestIterator, class DestAccessor>
731void boundaryTensor1(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
732 DestIterator dupperleft, DestAccessor dest,
733 double scale)
734{
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.");
739
740 detail::evenPolarFilters(supperleft, slowerright, src,
741 dupperleft, dest, scale, true);
742 detail::oddPolarFilters(supperleft, slowerright, src,
743 dupperleft, dest, scale, true);
744}
745
746template <class SrcIterator, class SrcAccessor,
747 class DestIterator, class DestAccessor>
748inline
749void boundaryTensor1(triple<SrcIterator, SrcIterator, SrcAccessor> src,
750 pair<DestIterator, DestAccessor> dest,
751 double scale)
752{
753 boundaryTensor1(src.first, src.second, src.third,
754 dest.first, dest.second, scale);
755}
756
757template <class T1, class S1,
758 class T2, class S2>
759inline void
762 double scale)
763{
764 vigra_precondition(src.shape() == dest.shape(),
765 "boundaryTensor1(): shape mismatch between input and output.");
766 boundaryTensor1(srcImageRange(src),
767 destImage(dest), scale);
768}
769
770/********************************************************/
771/* */
772/* boundaryTensor3 */
773/* */
774/********************************************************/
775
776/* Add 3rd order Riesz transform to boundary tensor
777 ??? Does not work -- bug or too coarse approximation for 3rd order ???
778*/
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,
785 double scale)
786{
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.");
791
792 detail::evenPolarFilters(supperleft, slowerright, sa,
793 dupperleft_even, even, scale, false);
794
795 int w = slowerright.x - supperleft.x;
796 int h = slowerright.y - supperleft.y;
797
798 typedef typename
799 NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
800 typedef BasicImage<TinyVector<TmpType, 4> > TmpImage;
801 TmpImage t1(w, h), t2(w, h);
802
803 detail::KernelArray k1, k3;
804 detail::initGaussianPolarFilters1(scale, k1);
805 detail::initGaussianPolarFilters3(scale, k3);
806
807 // calculate filter responses for odd filters
809 convolveImage(srcIterRange(supperleft, slowerright, sa),
810 destImage(t1, tmpBand), k1[3], k1[0]);
811 tmpBand.setIndex(1);
812 convolveImage(srcIterRange(supperleft, slowerright, sa),
813 destImage(t1, tmpBand), k1[1], k1[2]);
814 tmpBand.setIndex(2);
815 convolveImage(srcIterRange(supperleft, slowerright, sa),
816 destImage(t1, tmpBand), k3[3], k3[0]);
817 tmpBand.setIndex(3);
818 convolveImage(srcIterRange(supperleft, slowerright, sa),
819 destImage(t1, tmpBand), k3[1], k3[2]);
820
821 tmpBand.setIndex(0);
822 convolveImage(srcIterRange(supperleft, slowerright, sa),
823 destImage(t2, tmpBand), k1[0], k1[3]);
824 tmpBand.setIndex(1);
825 convolveImage(srcIterRange(supperleft, slowerright, sa),
826 destImage(t2, tmpBand), k1[2], k1[1]);
827 tmpBand.setIndex(2);
828 convolveImage(srcIterRange(supperleft, slowerright, sa),
829 destImage(t2, tmpBand), k3[0], k3[3]);
830 tmpBand.setIndex(3);
831 convolveImage(srcIterRange(supperleft, slowerright, sa),
832 destImage(t2, tmpBand), k3[2], k3[1]);
833
834 // create odd tensor from filter responses
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)
840 {
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)
846 {
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);
857
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);
861 }
862 }
863}
864
865template <class SrcIterator, class SrcAccessor,
866 class DestIteratorEven, class DestAccessorEven,
867 class DestIteratorOdd, class DestAccessorOdd>
868inline
869void boundaryTensor3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
870 pair<DestIteratorEven, DestAccessorEven> even,
871 pair<DestIteratorOdd, DestAccessorOdd> odd,
872 double scale)
873{
874 boundaryTensor3(src.first, src.second, src.third,
875 even.first, even.second, odd.first, odd.second, scale);
876}
877
878template <class T1, class S1,
879 class T2E, class S2Even,
880 class T2O, class S2Odd>
881inline
882void boundaryTensor3(MultiArrayView<2, T1, S1> const & src,
885 double scale)
886{
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);
891}
892
893//@}
894
895} // namespace vigra
896
897#endif // VIGRA_BOUNDARYTENSOR_HXX
Fundamental class template for images.
Definition basicimage.hxx:476
Base class for, and view to, MultiArray.
Definition multi_array.hxx:705
Accessor for one component of a vector.
Definition accessor.hxx:540
NumericTraits< T >::Promote sq(T t)
The square function.
Definition mathutil.hxx:382
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.

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.12.2 (Mon Apr 14 2025)