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

boundarytensor.hxx VIGRA

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 CommonConvolutionFilters Common Filters
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
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
void rieszTransformOfLOG(...)
Calculate Riesz transforms of the Laplacian of Gaussian.
Definition: accessor.hxx:43
bool odd(int t)
Check if an integer is odd.
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:344
void combineTwoImages(...)
Combine two source images into destination image.
void convolveImage(...)
Convolve an image with the given kernel(s).
bool even(int t)
Check if an integer is even.
void boundaryTensor1(...)
Boundary tensor variant.
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.10.0