[ 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 
49 namespace vigra {
50 
51 namespace detail {
52 
53 /***********************************************************************/
54 
55 typedef ArrayVector<Kernel1D<double> > KernelArray;
56 
57 template <class KernelArray>
58 void
59 initGaussianPolarFilters1(double std_dev, KernelArray & k)
60 {
61  typedef typename KernelArray::value_type Kernel;
62  typedef typename Kernel::iterator iterator;
63 
64  vigra_precondition(std_dev >= 0.0,
65  "initGaussianPolarFilter1(): "
66  "Standard deviation must be >= 0.");
67 
68  k.resize(4);
69 
70  int radius = (int)(4.0*std_dev + 0.5);
71  std_dev *= 1.08179074376;
72  double f = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / std_dev; // norm
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;
76 
77 
78  for(unsigned int i=0; i<k.size(); ++i)
79  {
80  k[i].initExplicitly(-radius, radius);
81  k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
82  }
83 
84  int ix;
85  iterator c = k[0].center();
86  for(ix=-radius; ix<=radius; ++ix)
87  {
88  double x = (double)ix;
89  c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x);
90  }
91 
92  c = k[1].center();
93  for(ix=-radius; ix<=radius; ++ix)
94  {
95  double x = (double)ix;
96  c[ix] = f * x * VIGRA_CSTD::exp(sigma22 * x * x);
97  }
98 
99  c = k[2].center();
100  double b2 = b / 3.0;
101  for(ix=-radius; ix<=radius; ++ix)
102  {
103  double x = (double)ix;
104  c[ix] = f * (b2 + a * x * x) * VIGRA_CSTD::exp(sigma22 * x * x);
105  }
106 
107  c = k[3].center();
108  for(ix=-radius; ix<=radius; ++ix)
109  {
110  double x = (double)ix;
111  c[ix] = f * x * (b + a * x * x) * VIGRA_CSTD::exp(sigma22 * x * x);
112  }
113 }
114 
115 template <class KernelArray>
116 void
117 initGaussianPolarFilters2(double std_dev, KernelArray & k)
118 {
119  typedef typename KernelArray::value_type Kernel;
120  typedef typename Kernel::iterator iterator;
121 
122  vigra_precondition(std_dev >= 0.0,
123  "initGaussianPolarFilter2(): "
124  "Standard deviation must be >= 0.");
125 
126  k.resize(3);
127 
128  int radius = (int)(4.0*std_dev + 0.5);
129  double f = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / std_dev; // norm
130  double sigma2 = std_dev*std_dev;
131  double sigma22 = -0.5 / sigma2;
132 
133  for(unsigned int i=0; i<k.size(); ++i)
134  {
135  k[i].initExplicitly(-radius, radius);
136  k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
137  }
138 
139  int ix;
140  iterator c = k[0].center();
141  for(ix=-radius; ix<=radius; ++ix)
142  {
143  double x = (double)ix;
144  c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x);
145  }
146 
147  c = k[1].center();
148  double f1 = f / sigma2;
149  for(ix=-radius; ix<=radius; ++ix)
150  {
151  double x = (double)ix;
152  c[ix] = f1 * x * VIGRA_CSTD::exp(sigma22 * x * x);
153  }
154 
155  c = k[2].center();
156  double f2 = f / (sigma2 * sigma2);
157  for(ix=-radius; ix<=radius; ++ix)
158  {
159  double x = (double)ix;
160  c[ix] = f2 * (x * x - sigma2) * VIGRA_CSTD::exp(sigma22 * x * x);
161  }
162 }
163 
164 template <class KernelArray>
165 void
166 initGaussianPolarFilters3(double std_dev, KernelArray & k)
167 {
168  typedef typename KernelArray::value_type Kernel;
169  typedef typename Kernel::iterator iterator;
170 
171  vigra_precondition(std_dev >= 0.0,
172  "initGaussianPolarFilter3(): "
173  "Standard deviation must be >= 0.");
174 
175  k.resize(4);
176 
177  int radius = (int)(4.0*std_dev + 0.5);
178  std_dev *= 1.15470053838;
179  double sigma22 = -0.5 / std_dev / std_dev;
180  double f = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / std_dev; // norm
181  double a = 0.883887052922 / VIGRA_CSTD::pow(std_dev, 5);
182 
183  for(unsigned int i=0; i<k.size(); ++i)
184  {
185  k[i].initExplicitly(-radius, radius);
186  k[i].setBorderTreatment(BORDER_TREATMENT_REFLECT);
187  }
188 
189  //double b = -1.3786348292 / VIGRA_CSTD::pow(std_dev, 3);
190 
191  int ix;
192  iterator c = k[0].center();
193  for(ix=-radius; ix<=radius; ++ix)
194  {
195  double x = (double)ix;
196  c[ix] = f * VIGRA_CSTD::exp(sigma22 * x * x);
197  }
198 
199  c = k[1].center();
200  for(ix=-radius; ix<=radius; ++ix)
201  {
202  double x = (double)ix;
203  c[ix] = f * x * VIGRA_CSTD::exp(sigma22 * x * x);
204  }
205 
206  c = k[2].center();
207  double a2 = 3.0 * a;
208  for(ix=-radius; ix<=radius; ++ix)
209  {
210  double x = (double)ix;
211  c[ix] = f * a2 * x * x * VIGRA_CSTD::exp(sigma22 * x * x);
212  }
213 
214  c = k[3].center();
215  for(ix=-radius; ix<=radius; ++ix)
216  {
217  double x = (double)ix;
218  c[ix] = f * a * x * x * x * VIGRA_CSTD::exp(sigma22 * x * x);
219  }
220 }
221 
222 template <class SrcIterator, class SrcAccessor,
223  class DestIterator, class DestAccessor>
224 void
225 evenPolarFilters(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
226  DestIterator dupperleft, DestAccessor dest,
227  double scale, bool noLaplacian)
228 {
229  vigra_precondition(dest.size(dupperleft) == 3,
230  "evenPolarFilters(): image for even output must have 3 bands.");
231 
232  int w = slowerright.x - supperleft.x;
233  int h = slowerright.y - supperleft.y;
234 
235  typedef typename
236  NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
237  typedef BasicImage<TinyVector<TmpType, 3> > TmpImage;
238  typedef typename TmpImage::traverser TmpTraverser;
239  TmpImage t(w, h);
240 
241  KernelArray k2;
242  initGaussianPolarFilters2(scale, k2);
243 
244  // calculate filter responses for even filters
245  VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t.accessor());
246  convolveImage(srcIterRange(supperleft, slowerright, src),
247  destImage(t, tmpBand), k2[2], k2[0]);
248  tmpBand.setIndex(1);
249  convolveImage(srcIterRange(supperleft, slowerright, src),
250  destImage(t, tmpBand), k2[1], k2[1]);
251  tmpBand.setIndex(2);
252  convolveImage(srcIterRange(supperleft, slowerright, src),
253  destImage(t, tmpBand), k2[0], k2[2]);
254 
255  // create even tensor from filter responses
256  TmpTraverser tul(t.upperLeft());
257  TmpTraverser tlr(t.lowerRight());
258  for(; tul.y != tlr.y; ++tul.y, ++dupperleft.y)
259  {
260  typename TmpTraverser::row_iterator tr = tul.rowIterator();
261  typename TmpTraverser::row_iterator trend = tr + w;
262  typename DestIterator::row_iterator d = dupperleft.rowIterator();
263  if(noLaplacian)
264  {
265  for(; tr != trend; ++tr, ++d)
266  {
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);
271  }
272  }
273  else
274  {
275  for(; tr != trend; ++tr, ++d)
276  {
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);
280  }
281  }
282  }
283 }
284 
285 template <class SrcIterator, class SrcAccessor,
286  class DestIterator, class DestAccessor>
287 void
288 oddPolarFilters(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
289  DestIterator dupperleft, DestAccessor dest,
290  double scale, bool addResult)
291 {
292  vigra_precondition(dest.size(dupperleft) == 3,
293  "oddPolarFilters(): image for odd output must have 3 bands.");
294 
295  int w = slowerright.x - supperleft.x;
296  int h = slowerright.y - supperleft.y;
297 
298  typedef typename
299  NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
300  typedef BasicImage<TinyVector<TmpType, 4> > TmpImage;
301  typedef typename TmpImage::traverser TmpTraverser;
302  TmpImage t(w, h);
303 
304  detail::KernelArray k1;
305  detail::initGaussianPolarFilters1(scale, k1);
306 
307  // calculate filter responses for odd filters
308  VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t.accessor());
309  convolveImage(srcIterRange(supperleft, slowerright, src),
310  destImage(t, tmpBand), k1[3], k1[0]);
311  tmpBand.setIndex(1);
312  convolveImage(srcIterRange(supperleft, slowerright, src),
313  destImage(t, tmpBand), k1[2], k1[1]);
314  tmpBand.setIndex(2);
315  convolveImage(srcIterRange(supperleft, slowerright, src),
316  destImage(t, tmpBand), k1[1], k1[2]);
317  tmpBand.setIndex(3);
318  convolveImage(srcIterRange(supperleft, slowerright, src),
319  destImage(t, tmpBand), k1[0], k1[3]);
320 
321  // create odd tensor from filter responses
322  TmpTraverser tul(t.upperLeft());
323  TmpTraverser tlr(t.lowerRight());
324  for(; tul.y != tlr.y; ++tul.y, ++dupperleft.y)
325  {
326  typename TmpTraverser::row_iterator tr = tul.rowIterator();
327  typename TmpTraverser::row_iterator trend = tr + w;
328  typename DestIterator::row_iterator d = dupperleft.rowIterator();
329  if(addResult)
330  {
331  for(; tr != trend; ++tr, ++d)
332  {
333  TmpType d0 = (*tr)[0] + (*tr)[2];
334  TmpType d1 = -(*tr)[1] - (*tr)[3];
335 
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);
339  }
340  }
341  else
342  {
343  for(; tr != trend; ++tr, ++d)
344  {
345  TmpType d0 = (*tr)[0] + (*tr)[2];
346  TmpType d1 = -(*tr)[1] - (*tr)[3];
347 
348  dest.setComponent(sq(d0), d, 0);
349  dest.setComponent(d0 * d1, d, 1);
350  dest.setComponent(sq(d1), d, 2);
351  }
352  }
353  }
354 }
355 
356 } // namespace detail
357 
358 /** \addtogroup CommonConvolutionFilters Common Filters
359 */
360 //@{
361 
362 /********************************************************/
363 /* */
364 /* rieszTransformOfLOG */
365 /* */
366 /********************************************************/
367 
368 /** \brief Calculate Riesz transforms of the Laplacian of Gaussian.
369 
370  The Riesz transforms of the Laplacian of Gaussian have the following transfer
371  functions (defined in a polar coordinate representation of the frequency domain):
372 
373  \f[
374  F_{\sigma}(r, \phi)=(i \cos \phi)^n (i \sin \phi)^m r^2 e^{-r^2 \sigma^2 / 2}
375  \f]
376 
377  where <i>n</i> = <tt>xorder</tt> and <i>m</i> = <tt>yorder</tt> determine th e
378  order of the transform, and <tt>sigma > 0</tt> is the scale of the Laplacian
379  of Gaussian. This function computes a good spatial domain approximation of
380  these transforms for <tt>xorder + yorder <= 2</tt>. The filter responses may be used
381  to calculate the monogenic signal or the boundary tensor.
382 
383  <b> Declarations:</b>
384 
385  pass arguments explicitly:
386  \code
387  namespace vigra {
388  template <class SrcIterator, class SrcAccessor,
389  class DestIterator, class DestAccessor>
390  void rieszTransformOfLOG(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
391  DestIterator dupperleft, DestAccessor dest,
392  double scale, unsigned int xorder, unsigned int yorder);
393  }
394  \endcode
395 
396 
397  use argument objects in conjunction with \ref ArgumentObjectFactories :
398  \code
399  namespace vigra {
400  template <class SrcIterator, class SrcAccessor,
401  class DestIterator, class DestAccessor>
402  void rieszTransformOfLOG(triple<SrcIterator, SrcIterator, SrcAccessor> src,
403  pair<DestIterator, DestAccessor> dest,
404  double scale, unsigned int xorder, unsigned int yorder);
405  }
406  \endcode
407 
408  <b> Usage:</b>
409 
410  <b>\#include</b> <vigra/boundarytensor.hxx>
411 
412  \code
413  FImage impulse(17,17), res(17, 17);
414  impulse(8,8) = 1.0;
415 
416  // calculate the impulse response of the first order Riesz transform in x-direction
417  rieszTransformOfLOG(srcImageRange(impulse), destImage(res), 2.0, 1, 0);
418  \endcode
419 
420 */
421 doxygen_overloaded_function(template <...> void rieszTransformOfLOG)
422 
423 template <class SrcIterator, class SrcAccessor,
424  class DestIterator, class DestAccessor>
425 void rieszTransformOfLOG(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
426  DestIterator dupperleft, DestAccessor dest,
427  double scale, unsigned int xorder, unsigned int yorder)
428 {
429  unsigned int order = xorder + yorder;
430 
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.");
435 
436  int w = slowerright.x - supperleft.x;
437  int h = slowerright.y - supperleft.y;
438 
439  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
440  typedef BasicImage<TmpType> TmpImage;
441 
442  switch(order)
443  {
444  case 0:
445  {
446  detail::KernelArray k2;
447  detail::initGaussianPolarFilters2(scale, k2);
448 
449  TmpImage t1(w, h), t2(w, h);
450 
451  convolveImage(srcIterRange(supperleft, slowerright, src),
452  destImage(t1), k2[2], k2[0]);
453  convolveImage(srcIterRange(supperleft, slowerright, src),
454  destImage(t2), k2[0], k2[2]);
455  combineTwoImages(srcImageRange(t1), srcImage(t2),
456  destIter(dupperleft, dest), std::plus<TmpType>());
457  break;
458  }
459  case 1:
460  {
461  detail::KernelArray k1;
462  detail::initGaussianPolarFilters1(scale, k1);
463 
464  TmpImage t1(w, h), t2(w, h);
465 
466  if(xorder == 1)
467  {
468  convolveImage(srcIterRange(supperleft, slowerright, src),
469  destImage(t1), k1[3], k1[0]);
470  convolveImage(srcIterRange(supperleft, slowerright, src),
471  destImage(t2), k1[1], k1[2]);
472  }
473  else
474  {
475  convolveImage(srcIterRange(supperleft, slowerright, src),
476  destImage(t1), k1[0], k1[3]);
477  convolveImage(srcIterRange(supperleft, slowerright, src),
478  destImage(t2), k1[2], k1[1]);
479  }
480  combineTwoImages(srcImageRange(t1), srcImage(t2),
481  destIter(dupperleft, dest), std::plus<TmpType>());
482  break;
483  }
484  case 2:
485  {
486  detail::KernelArray k2;
487  detail::initGaussianPolarFilters2(scale, k2);
488 
489  convolveImage(srcIterRange(supperleft, slowerright, src),
490  destIter(dupperleft, dest), k2[xorder], k2[yorder]);
491  break;
492  }
493  /* for test purposes only: compute 3rd order polar filters */
494  case 3:
495  {
496  detail::KernelArray k3;
497  detail::initGaussianPolarFilters3(scale, k3);
498 
499  TmpImage t1(w, h), t2(w, h);
500 
501  if(xorder == 3)
502  {
503  convolveImage(srcIterRange(supperleft, slowerright, src),
504  destImage(t1), k3[3], k3[0]);
505  convolveImage(srcIterRange(supperleft, slowerright, src),
506  destImage(t2), k3[1], k3[2]);
507  }
508  else
509  {
510  convolveImage(srcIterRange(supperleft, slowerright, src),
511  destImage(t1), k3[0], k3[3]);
512  convolveImage(srcIterRange(supperleft, slowerright, src),
513  destImage(t2), k3[2], k3[1]);
514  }
515  combineTwoImages(srcImageRange(t1), srcImage(t2),
516  destIter(dupperleft, dest), std::minus<TmpType>());
517  break;
518  }
519  }
520 }
521 
522 template <class SrcIterator, class SrcAccessor,
523  class DestIterator, class DestAccessor>
524 inline
525 void rieszTransformOfLOG(triple<SrcIterator, SrcIterator, SrcAccessor> src,
526  pair<DestIterator, DestAccessor> dest,
527  double scale, unsigned int xorder, unsigned int yorder)
528 {
529  rieszTransformOfLOG(src.first, src.second, src.third, dest.first, dest.second,
530  scale, xorder, yorder);
531 }
532 //@}
533 
534 /** \addtogroup TensorImaging Tensor Image Processing
535 */
536 //@{
537 
538 /********************************************************/
539 /* */
540 /* boundaryTensor */
541 /* */
542 /********************************************************/
543 
544 /** \brief Calculate the boundary tensor for a scalar valued image.
545 
546  These functions calculate a spatial domain approximation of
547  the boundary tensor as described in
548 
549  U. K&ouml;the: <a href="http://hci.iwr.uni-heidelberg.de/people/ukoethe/papers/index.php#cite_polarfilters">
550  <i>"Integrated Edge and Junction Detection with the Boundary Tensor"</i></a>,
551  in: ICCV 03, Proc. of 9th Intl. Conf. on Computer Vision, Nice 2003, vol. 1,
552  pp. 424-431, Los Alamitos: IEEE Computer Society, 2003
553 
554  with the Laplacian of Gaussian as the underlying bandpass filter (see
555  \ref rieszTransformOfLOG()). The output image must have 3 bands which will hold the
556  tensor components in the order t11, t12 (== t21), t22. The function
557  \ref boundaryTensor1() with the same interface implements a variant of the
558  boundary tensor where the 0th-order Riesz transform has been dropped, so that the
559  tensor is no longer sensitive to blobs.
560 
561  <b> Declarations:</b>
562 
563  pass arguments explicitly:
564  \code
565  namespace vigra {
566  template <class SrcIterator, class SrcAccessor,
567  class DestIterator, class DestAccessor>
568  void boundaryTensor(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
569  DestIterator dupperleft, DestAccessor dest,
570  double scale);
571  }
572  \endcode
573 
574  use argument objects in conjunction with \ref ArgumentObjectFactories :
575  \code
576  namespace vigra {
577  template <class SrcIterator, class SrcAccessor,
578  class DestIterator, class DestAccessor>
579  void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src,
580  pair<DestIterator, DestAccessor> dest,
581  double scale);
582  }
583  \endcode
584 
585  <b> Usage:</b>
586 
587  <b>\#include</b> <vigra/boundarytensor.hxx>
588 
589  \code
590  FImage img(w,h);
591  FVector3Image bt(w,h);
592  ...
593  boundaryTensor(srcImageRange(img), destImage(bt), 2.0);
594  \endcode
595 
596 */
597 doxygen_overloaded_function(template <...> void boundaryTensor)
598 
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,
603  double scale)
604 {
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.");
609 
610  detail::evenPolarFilters(supperleft, slowerright, src,
611  dupperleft, dest, scale, false);
612  detail::oddPolarFilters(supperleft, slowerright, src,
613  dupperleft, dest, scale, true);
614 }
615 
616 template <class SrcIterator, class SrcAccessor,
617  class DestIterator, class DestAccessor>
618 inline
619 void boundaryTensor(triple<SrcIterator, SrcIterator, SrcAccessor> src,
620  pair<DestIterator, DestAccessor> dest,
621  double scale)
622 {
623  boundaryTensor(src.first, src.second, src.third,
624  dest.first, dest.second, scale);
625 }
626 
627 /** \brief Boundary tensor variant.
628 
629  This function implements a variant of the boundary tensor where the
630  0th-order Riesz transform has been dropped, so that the tensor is no
631  longer sensitive to blobs. See \ref boundaryTensor() for more detailed
632  documentation.
633 
634  <b> Declarations:</b>
635 
636  <b>\#include</b> <vigra/boundarytensor.hxx>
637 
638  pass arguments explicitly:
639  \code
640  namespace vigra {
641  template <class SrcIterator, class SrcAccessor,
642  class DestIterator, class DestAccessor>
643  void boundaryTensor1(SrcIterator supperleft, SrcIterator slowerright, SrcAccessor src,
644  DestIterator dupperleft, DestAccessor dest,
645  double scale);
646  }
647  \endcode
648 
649  use argument objects in conjunction with \ref ArgumentObjectFactories :
650  \code
651  namespace vigra {
652  template <class SrcIterator, class SrcAccessor,
653  class DestIterator, class DestAccessor>
654  void boundaryTensor1(triple<SrcIterator, SrcIterator, SrcAccessor> src,
655  pair<DestIterator, DestAccessor> dest,
656  double scale);
657  }
658  \endcode
659 */
660 doxygen_overloaded_function(template <...> void boundaryTensor1)
661 
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,
666  double scale)
667 {
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.");
672 
673  detail::evenPolarFilters(supperleft, slowerright, src,
674  dupperleft, dest, scale, true);
675  detail::oddPolarFilters(supperleft, slowerright, src,
676  dupperleft, dest, scale, true);
677 }
678 
679 template <class SrcIterator, class SrcAccessor,
680  class DestIterator, class DestAccessor>
681 inline
682 void boundaryTensor1(triple<SrcIterator, SrcIterator, SrcAccessor> src,
683  pair<DestIterator, DestAccessor> dest,
684  double scale)
685 {
686  boundaryTensor1(src.first, src.second, src.third,
687  dest.first, dest.second, scale);
688 }
689 
690 /********************************************************/
691 /* */
692 /* boundaryTensor3 */
693 /* */
694 /********************************************************/
695 
696 /* Add 3rd order Riesz transform to boundary tensor
697  ??? Does not work -- bug or too coarse approximation for 3rd order ???
698 */
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,
705  double scale)
706 {
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.");
711 
712  detail::evenPolarFilters(supperleft, slowerright, sa,
713  dupperleft_even, even, scale, false);
714 
715  int w = slowerright.x - supperleft.x;
716  int h = slowerright.y - supperleft.y;
717 
718  typedef typename
719  NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
720  typedef BasicImage<TinyVector<TmpType, 4> > TmpImage;
721  TmpImage t1(w, h), t2(w, h);
722 
723  detail::KernelArray k1, k3;
724  detail::initGaussianPolarFilters1(scale, k1);
725  detail::initGaussianPolarFilters3(scale, k3);
726 
727  // calculate filter responses for odd filters
728  VectorElementAccessor<typename TmpImage::Accessor> tmpBand(0, t1.accessor());
729  convolveImage(srcIterRange(supperleft, slowerright, sa),
730  destImage(t1, tmpBand), k1[3], k1[0]);
731  tmpBand.setIndex(1);
732  convolveImage(srcIterRange(supperleft, slowerright, sa),
733  destImage(t1, tmpBand), k1[1], k1[2]);
734  tmpBand.setIndex(2);
735  convolveImage(srcIterRange(supperleft, slowerright, sa),
736  destImage(t1, tmpBand), k3[3], k3[0]);
737  tmpBand.setIndex(3);
738  convolveImage(srcIterRange(supperleft, slowerright, sa),
739  destImage(t1, tmpBand), k3[1], k3[2]);
740 
741  tmpBand.setIndex(0);
742  convolveImage(srcIterRange(supperleft, slowerright, sa),
743  destImage(t2, tmpBand), k1[0], k1[3]);
744  tmpBand.setIndex(1);
745  convolveImage(srcIterRange(supperleft, slowerright, sa),
746  destImage(t2, tmpBand), k1[2], k1[1]);
747  tmpBand.setIndex(2);
748  convolveImage(srcIterRange(supperleft, slowerright, sa),
749  destImage(t2, tmpBand), k3[0], k3[3]);
750  tmpBand.setIndex(3);
751  convolveImage(srcIterRange(supperleft, slowerright, sa),
752  destImage(t2, tmpBand), k3[2], k3[1]);
753 
754  // create odd tensor from filter responses
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)
760  {
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)
766  {
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);
777 
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);
781  }
782  }
783 }
784 
785 template <class SrcIterator, class SrcAccessor,
786  class DestIteratorEven, class DestAccessorEven,
787  class DestIteratorOdd, class DestAccessorOdd>
788 inline
789 void boundaryTensor3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
790  pair<DestIteratorEven, DestAccessorEven> even,
791  pair<DestIteratorOdd, DestAccessorOdd> odd,
792  double scale)
793 {
794  boundaryTensor3(src.first, src.second, src.third,
795  even.first, even.second, odd.first, odd.second, scale);
796 }
797 
798 //@}
799 
800 } // namespace vigra
801 
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.
bool odd(int t)
NumericTraits< T >::Promote sq(T t)
Definition: mathutil.hxx:341
void combineTwoImages(...)
Combine two source images into destination image.
bool even(int t)
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

© 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.9.0 (Sun Aug 10 2014)