[ 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 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::rieszTransformOfLOG
void rieszTransformOfLOG(...)
Calculate Riesz transforms of the Laplacian of Gaussian.
vigra::convolveImage
void convolveImage(...)
Convolve an image with the given kernel(s).
vigra::linalg::exp
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
vigra::odd
bool odd(int t)
Check if an integer is odd.
vigra::sq
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:382
vigra::combineTwoImages
void combineTwoImages(...)
Combine two source images into destination image.
vigra::doxygen_overloaded_function
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
vigra::even
bool even(int t)
Check if an integer is even.
vigra::boundaryTensor1
void boundaryTensor1(...)
Boundary tensor variant.
vigra::boundaryTensor
void boundaryTensor(...)
Calculate the boundary tensor for a scalar valued image.
vigra::sqrt
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.11.1 (Fri May 19 2017)