[ 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
50 namespace vigra {
51
52 namespace detail {
53
54 /***********************************************************************/
55
56 typedef ArrayVector<Kernel1D<double> > KernelArray;
57
58 template <class KernelArray>
59 void
60 initGaussianPolarFilters1(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
116 template <class KernelArray>
117 void
118 initGaussianPolarFilters2(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
165 template <class KernelArray>
166 void
167 initGaussianPolarFilters3(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
223 template <class SrcIterator, class SrcAccessor,
224 class DestIterator, class DestAccessor>
225 void
226 evenPolarFilters(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
286 template <class SrcIterator, class SrcAccessor,
287 class DestIterator, class DestAccessor>
288 void
289 oddPolarFilters(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 */
434 doxygen_overloaded_function(template <...> void rieszTransformOfLOG)
435
436 template <class SrcIterator, class SrcAccessor,
437 class DestIterator, class DestAccessor>
438 void 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
535 template <class SrcIterator, class SrcAccessor,
536 class DestIterator, class DestAccessor>
537 inline
538 void 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
546 template <class T1, class S1,
547 class T2, class S2>
548 inline void
549 rieszTransformOfLOG(MultiArrayView<2, T1, S1> const & src,
550 MultiArrayView<2, T2, S2> dest,
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 */
637 doxygen_overloaded_function(template <...> void boundaryTensor)
638
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,
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
656 template <class SrcIterator, class SrcAccessor,
657 class DestIterator, class DestAccessor>
658 inline
659 void 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
667 template <class T1, class S1,
668 class T2, class S2>
669 inline void
670 boundaryTensor(MultiArrayView<2, T1, S1> const & src,
671 MultiArrayView<2, T2, S2> dest,
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 */
727 doxygen_overloaded_function(template <...> void boundaryTensor1)
728
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,
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
746 template <class SrcIterator, class SrcAccessor,
747 class DestIterator, class DestAccessor>
748 inline
749 void 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
757 template <class T1, class S1,
758 class T2, class S2>
759 inline void
760 boundaryTensor1(MultiArrayView<2, T1, S1> const & src,
761 MultiArrayView<2, T2, S2> dest,
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 */
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,
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
808 VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t1.accessor());
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
865 template <class SrcIterator, class SrcAccessor,
866 class DestIteratorEven, class DestAccessorEven,
867 class DestIteratorOdd, class DestAccessorOdd>
868 inline
869 void 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
878 template <class T1, class S1,
879 class T2E, class S2Even,
880 class T2O, class S2Odd>
881 inline
882 void boundaryTensor3(MultiArrayView<2, T1, S1> const & src,
883 MultiArrayView<2, T2E, S2Even> even,
884 MultiArrayView<2, T2O, S2Odd> odd,
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
vigra::sq
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:382
vigra::rieszTransformOfLOG
void rieszTransformOfLOG(...)
Calculate Riesz transforms of the Laplacian of Gaussian.
vigra
vigra::doxygen_overloaded_function
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
vigra::boundaryTensor1
void boundaryTensor1(...)
Boundary tensor variant.
vigra::convolveImage
void convolveImage(...)
Convolve an image with the given kernel(s).
utilities.hxx
vigra::even
bool even(int t)
Check if an integer is even.
vigra::combineTwoImages
void combineTwoImages(...)
Combine two source images into destination image.
vigra::boundaryTensor
void boundaryTensor(...)
Calculate the boundary tensor for a scalar valued image.
vigra::odd
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.11.1