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

orientedtensorfilters.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2002-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 #ifndef VIGRA_ORIENTEDTENSORFILTERS_HXX
37 #define VIGRA_ORIENTEDTENSORFILTERS_HXX
38 
39 #include <cmath>
40 #include "utilities.hxx"
41 #include "initimage.hxx"
42 #include "stdconvolution.hxx"
43 
44 namespace vigra {
45 
46 /** \addtogroup TensorImaging Tensor Image Processing
47 */
48 //@{
49 
50 /********************************************************/
51 /* */
52 /* hourGlassFilter */
53 /* */
54 /********************************************************/
55 
56 /** \brief Anisotropic tensor smoothing with the hourglass filter.
57 
58  This function implements anisotropic tensor smoothing by an
59  hourglass-shaped filters as described in
60 
61  U. K&ouml;the: <a href="http://hci.iwr.uni-heidelberg.de/people/ukoethe/papers/index.php#cite_structureTensor">
62  <i>"Edge and Junction Detection with an Improved Structure Tensor"</i></a>,
63  in: Proc. of 25th DAGM Symposium, Magdeburg 2003, Lecture Notes in Computer Science 2781,
64  pp. 25-32, Heidelberg: Springer, 2003
65 
66  It is closely related to the structure tensor (see \ref structureTensor()), but
67  replaces the linear tensor smoothing with a smoothing along edges only.
68  Smoothing across edges is largely suppressed. This means that the
69  image structure is preserved much better because nearby features
70  such as parallel edges are not blended into each other.
71 
72  The hourglass filter is typically applied to a gradient tensor, i.e. the
73  Euclidean product of the gradient with itself, which can be obtained by a
74  gradient operator followed with \ref vectorToTensor(), see example below.
75  The hourglass shape of the filter can be interpreted as indicating the likely
76  continuations of a local edge element. The parameter <tt>sigma</tt> determines
77  the radius of the hourglass (i.e. how far the influence of the edge element
78  reaches), and <tt>rho</tt> controls its opening angle (i.e. how narrow the
79  edge orientation os followed). Recommended values are <tt>sigma = 1.4</tt>
80  (or, more generally, two to three times the scale of the gradient operator
81  used in the first step), and <tt>rho = 0.4</tt> which corresponds to an
82  opening angle of 22.5 degrees to either side of the edge.
83 
84  <b> Declarations:</b>
85 
86  pass arguments explicitly:
87  \code
88  namespace vigra {
89  template <class SrcIterator, class SrcAccessor,
90  class DestIterator, class DestAccessor>
91  void hourGlassFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
92  DestIterator dul, DestAccessor dest,
93  double sigma, double rho);
94  }
95  \endcode
96 
97 
98  use argument objects in conjunction with \ref ArgumentObjectFactories :
99  \code
100  namespace vigra {
101  template <class SrcIterator, class SrcAccessor,
102  class DestIterator, class DestAccessor>
103  inline
104  void hourGlassFilter(triple<SrcIterator, SrcIterator, SrcAccessor> s,
105  pair<DestIterator, DestAccessor> d,
106  double sigma, double rho);
107  }
108  \endcode
109 
110  <b> Usage:</b>
111 
112  <b>\#include</b> <vigra/orientedtensorfilters.hxx>
113 
114  \code
115  FImage img(w,h);
116  FVector2Image gradient(w,h);
117  FVector3Image tensor(w,h), smoothedTensor(w,h);
118 
119  gaussianGradient(srcImageRange(img), destImage(gradient), 1.0);
120  vectorToTensor(srcImageRange(gradient), destImage(tensor));
121  hourGlassFilter(srcImageRange(tensor), destImage(smoothedTensor), 2.0, 0.4);
122  \endcode
123 
124  \see vectorToTensor()
125 */
126 doxygen_overloaded_function(template <...> void hourGlassFilter)
127 
128 template <class SrcIterator, class SrcAccessor,
129  class DestIterator, class DestAccessor>
130 void hourGlassFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
131  DestIterator dul, DestAccessor dest,
132  double sigma, double rho)
133 {
134  vigra_precondition(sigma >= 0.0 && rho >= 0.0,
135  "hourGlassFilter(): sigma and rho must be >= 0.0");
136  vigra_precondition(src.size(sul) == 3,
137  "hourGlassFilter(): input image must have 3 bands.");
138  vigra_precondition(dest.size(dul) == 3,
139  "hourGlassFilter(): output image must have 3 bands.");
140 
141  // TODO: normalization
142 
143  int w = slr.x - sul.x;
144  int h = slr.y - sul.y;
145 
146  double radius = VIGRA_CSTD::floor(3.0*sigma + 0.5);
147  double sigma2 = -0.5 / sigma / sigma;
148  double rho2 = -0.5 / rho / rho;
149  double norm = 1.0 / (2.0 * M_PI * sigma * sigma);
150 
151  initImage(dul, dul+Diff2D(w,h), dest, NumericTraits<typename DestAccessor::value_type>::zero());
152 
153  for(int y=0; y<h; ++y, ++sul.y, ++dul.y)
154  {
155  SrcIterator s = sul;
156  DestIterator d = dul;
157  for(int x=0; x<w; ++x, ++s.x, ++d.x)
158  {
159  double phi = 0.5 * VIGRA_CSTD::atan2(
160  2.0*src.getComponent(s,1),
161  (double)src.getComponent(s,0) - src.getComponent(s,2));
162  double u = VIGRA_CSTD::sin(phi);
163  double v = VIGRA_CSTD::cos(phi);
164 
165  double x0 = x - radius < 0 ? -x : -radius;
166  double y0 = y - radius < 0 ? -y : -radius;
167  double x1 = x + radius >= w ? w - x - 1 : radius;
168  double y1 = y + radius >= h ? h - y - 1 : radius;
169 
170  DestIterator dwul = d + Diff2D((int)x0, (int)y0);
171 
172  for(double yy=y0; yy <= y1; ++yy, ++dwul.y)
173  {
174  typename DestIterator::row_iterator dw = dwul.rowIterator();
175  for(double xx=x0; xx <= x1; ++xx, ++dw)
176  {
177  double r2 = xx*xx + yy*yy;
178  double p = u*xx - v*yy;
179  double q = v*xx + u*yy;
180  double kernel = (p == 0.0) ?
181  (q == 0.0) ?
182  norm :
183  0.0 :
184  norm * VIGRA_CSTD::exp(sigma2*r2 + rho2*q*q/p/p);
185  dest.set(dest(dw) + kernel*src(s), dw);
186  }
187  }
188  }
189  }
190 }
191 
192 template <class SrcIterator, class SrcAccessor,
193  class DestIterator, class DestAccessor>
194 inline
195 void hourGlassFilter(triple<SrcIterator, SrcIterator, SrcAccessor> s,
196  pair<DestIterator, DestAccessor> d,
197  double sigma, double rho)
198 {
199  hourGlassFilter(s.first, s.second, s.third, d.first, d.second, sigma, rho);
200 }
201 
202 /********************************************************/
203 /* */
204 /* ellipticGaussian */
205 /* */
206 /********************************************************/
207 
208 template <class SrcIterator, class SrcAccessor,
209  class DestIterator, class DestAccessor>
210 void ellipticGaussian(SrcIterator sul, SrcIterator slr, SrcAccessor src,
211  DestIterator dul, DestAccessor dest,
212  double sigmax, double sigmin)
213 {
214  vigra_precondition(sigmax >= sigmin && sigmin >= 0.0,
215  "ellipticGaussian(): "
216  "sigmax >= sigmin and sigmin >= 0.0 required");
217 
218  int w = slr.x - sul.x;
219  int h = slr.y - sul.y;
220 
221  double radius = VIGRA_CSTD::floor(3.0*sigmax + 0.5);
222  double sigmin2 = -0.5 / sigmin / sigmin;
223  double norm = 1.0 / (2.0 * M_PI * sigmin * sigmax);
224 
225  initImage(dul, dul+Diff2D(w,h), dest, NumericTraits<typename DestAccessor::value_type>::zero());
226 
227  for(int y=0; y<h; ++y, ++sul.y, ++dul.y)
228  {
229  SrcIterator s = sul;
230  DestIterator d = dul;
231  for(int x=0; x<w; ++x, ++s.x, ++d.x)
232  {
233  typedef typename
234  NumericTraits<typename SrcAccessor::component_type>::RealPromote TmpType;
235  TmpType d1 = src.getComponent(s,0) + src.getComponent(s,2);
236  TmpType d2 = src.getComponent(s,0) - src.getComponent(s,2);
237  TmpType d3 = 2.0 * src.getComponent(s,1);
238  TmpType d4 = VIGRA_CSTD::sqrt(sq(d2) + sq(d3));
239  TmpType excentricity = 1.0 - (d1 - d4) / (d1 + d4);
240  double sigmax2 = -0.5 / sq((sigmax - sigmin)*excentricity + sigmin);
241 
242  double phi = 0.5 * VIGRA_CSTD::atan2(d3, d2);
243  double u = VIGRA_CSTD::sin(phi);
244  double v = VIGRA_CSTD::cos(phi);
245 
246  double x0 = x - radius < 0 ? -x : -radius;
247  double y0 = y - radius < 0 ? -y : -radius;
248  double x1 = x + radius >= w ? w - x - 1 : radius;
249  double y1 = y + radius >= h ? h - y - 1 : radius;
250 
251  DestIterator dwul = d + Diff2D((int)x0, (int)y0);
252 
253  for(double yy=y0; yy <= y1; ++yy, ++dwul.y)
254  {
255  typename DestIterator::row_iterator dw = dwul.rowIterator();
256  for(double xx=x0; xx <= x1; ++xx, ++dw)
257  {
258  double p = u*xx - v*yy;
259  double q = v*xx + u*yy;
260  double kernel = norm * VIGRA_CSTD::exp(sigmax2*p*p + sigmin2*q*q);
261  dest.set(dest(dw) + kernel*src(s), dw);
262  }
263  }
264  }
265  }
266 }
267 
268 template <class SrcIterator, class SrcAccessor,
269  class DestIterator, class DestAccessor>
270 inline
271 void ellipticGaussian(triple<SrcIterator, SrcIterator, SrcAccessor> s,
272  pair<DestIterator, DestAccessor> d,
273  double sigmax, double sigmin)
274 {
275  ellipticGaussian(s.first, s.second, s.third, d.first, d.second, sigmax, sigmin);
276 }
277 
278 /********************************************************/
279 /* */
280 /* kernels for orientedTrigonometricFilter */
281 /* */
282 /********************************************************/
283 
284 class FoerstnerKernelBase
285 {
286  public:
287  typedef double ResultType;
288  typedef double WeightType;
289  typedef DVector2Image::value_type VectorType;
290  typedef VectorType::value_type ValueType;
291 
292  FoerstnerKernelBase(double scale, bool ringShaped = false)
293  : radius_((int)(3.0*scale+0.5)),
294  weights_(2*radius_+1, 2*radius_+1),
295  vectors_(2*radius_+1, 2*radius_+1)
296  {
297  double norm = 1.0 / (2.0 * M_PI * scale * scale);
298  double s2 = -0.5 / scale / scale;
299 
300  for(int y = -radius_; y <= radius_; ++y)
301  {
302  for(int x = -radius_; x <= radius_; ++x)
303  {
304  double d2 = x*x + y*y;
305  double d = VIGRA_CSTD::sqrt(d2);
306  vectors_(x+radius_,y+radius_) = d != 0.0 ?
307  VectorType(x/d, -y/d) :
308  VectorType(ValueType(0), ValueType(0));
309  weights_(x+radius_,y+radius_) = ringShaped ?
310  norm * d2 * VIGRA_CSTD::exp(d2 * s2) :
311  norm * VIGRA_CSTD::exp(d2 * s2);
312  }
313  }
314  }
315 
316  ResultType operator()(int x, int y, VectorType const &) const
317  {
318  // isotropic filtering
319  return weights_(radius_, radius_);
320  }
321 
322  int radius_;
323  DImage weights_;
324  DVector2Image vectors_;
325 };
326 
327 class FoerstnerRingKernelBase
328 : public FoerstnerKernelBase
329 {
330  public:
331  FoerstnerRingKernelBase(double scale)
332  : FoerstnerKernelBase(scale, true)
333  {}
334 };
335 
336 class Cos2RingKernel
337 : public FoerstnerKernelBase
338 {
339  public:
340  typedef double ResultType;
341  typedef double WeightType;
342  typedef DVector2Image::value_type VectorType;
343 
344  Cos2RingKernel(double scale)
345  : FoerstnerKernelBase(scale, true)
346  {}
347 
348  ResultType operator()(int x, int y, VectorType const & v) const
349  {
350  if(x == 0 && y == 0)
351  return weights_(radius_, radius_);
352  double d = dot(vectors_(x+radius_, y+radius_), v);
353  return d * d * weights_(x+radius_, y+radius_);
354  }
355 };
356 
357 class Cos2Kernel
358 : public FoerstnerKernelBase
359 {
360  public:
361  typedef double ResultType;
362  typedef double WeightType;
363  typedef DVector2Image::value_type VectorType;
364 
365  Cos2Kernel(double scale)
366  : FoerstnerKernelBase(scale, false)
367  {}
368 
369  ResultType operator()(int x, int y, VectorType const & v) const
370  {
371  if(x == 0 && y == 0)
372  return weights_(radius_, radius_);
373  double d = dot(vectors_(x+radius_, y+radius_), v);
374  return d * d * weights_(x+radius_, y+radius_);
375  }
376 };
377 
378 class Sin2RingKernel
379 : public FoerstnerKernelBase
380 {
381  public:
382  typedef double ResultType;
383  typedef double WeightType;
384  typedef DVector2Image::value_type VectorType;
385 
386  Sin2RingKernel(double scale)
387  : FoerstnerKernelBase(scale, true)
388  {}
389 
390  ResultType operator()(int x, int y, VectorType const & v) const
391  {
392  if(x == 0 && y == 0)
393  return weights_(radius_, radius_);
394  double d = dot(vectors_(x+radius_, y+radius_), v);
395  return (1.0 - d * d) * weights_(x+radius_, y+radius_);
396  }
397 };
398 
399 class Sin2Kernel
400 : public FoerstnerKernelBase
401 {
402  public:
403  typedef double ResultType;
404  typedef double WeightType;
405  typedef DVector2Image::value_type VectorType;
406 
407  Sin2Kernel(double scale)
408  : FoerstnerKernelBase(scale, false)
409  {}
410 
411  ResultType operator()(int x, int y, VectorType const & v) const
412  {
413  if(x == 0 && y == 0)
414  return weights_(radius_, radius_);
415  double d = dot(vectors_(x+radius_, y+radius_), v);
416  return (1.0 - d * d) * weights_(x+radius_, y+radius_);
417  }
418 };
419 
420 class Sin6RingKernel
421 : public FoerstnerKernelBase
422 {
423  public:
424  typedef double ResultType;
425  typedef double WeightType;
426  typedef DVector2Image::value_type VectorType;
427 
428  Sin6RingKernel(double scale)
429  : FoerstnerKernelBase(scale, true)
430  {}
431 
432  ResultType operator()(int x, int y, VectorType const & v) const
433  {
434  if(x == 0 && y == 0)
435  return weights_(radius_, radius_);
436  double d = dot(vectors_(x+radius_, y+radius_), v);
437  return VIGRA_CSTD::pow(1.0 - d * d, 3) * weights_(x+radius_, y+radius_);
438  }
439 };
440 
441 class Sin6Kernel
442 : public FoerstnerKernelBase
443 {
444  public:
445  typedef double ResultType;
446  typedef double WeightType;
447  typedef DVector2Image::value_type VectorType;
448 
449  Sin6Kernel(double scale)
450  : FoerstnerKernelBase(scale, false)
451  {}
452 
453  ResultType operator()(int x, int y, VectorType const & v) const
454  {
455  if(x == 0 && y == 0)
456  return weights_(radius_, radius_);
457  double d = dot(vectors_(x+radius_, y+radius_), v);
458  return VIGRA_CSTD::pow(1.0 - d * d, 3) * weights_(x+radius_, y+radius_);
459  }
460 };
461 
462 class Cos6RingKernel
463 : public FoerstnerKernelBase
464 {
465  public:
466  typedef double ResultType;
467  typedef double WeightType;
468  typedef DVector2Image::value_type VectorType;
469 
470  Cos6RingKernel(double scale)
471  : FoerstnerKernelBase(scale, true)
472  {}
473 
474  ResultType operator()(int x, int y, VectorType const & v) const
475  {
476  if(x == 0 && y == 0)
477  return weights_(radius_, radius_);
478  double d = dot(vectors_(x+radius_, y+radius_), v);
479  return (1.0 - VIGRA_CSTD::pow(1.0 - d * d, 3)) * weights_(x+radius_, y+radius_);
480  }
481 };
482 
483 class Cos6Kernel
484 : public FoerstnerKernelBase
485 {
486  public:
487  typedef double ResultType;
488  typedef double WeightType;
489  typedef DVector2Image::value_type VectorType;
490 
491  Cos6Kernel(double scale)
492  : FoerstnerKernelBase(scale, false)
493  {}
494 
495  ResultType operator()(int x, int y, VectorType const & v) const
496  {
497  if(x == 0 && y == 0)
498  return weights_(radius_, radius_);
499  double d = dot(vectors_(x+radius_, y+radius_), v);
500  return (1.0 - VIGRA_CSTD::pow(1.0 - d * d, 3)) * weights_(x+radius_, y+radius_);
501  }
502 };
503 
504 /********************************************************/
505 /* */
506 /* orientedTrigonometricFilter */
507 /* */
508 /********************************************************/
509 
510 template <class SrcIterator, class SrcAccessor,
511  class DestIterator, class DestAccessor,
512  class Kernel>
513 void orientedTrigonometricFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
514  DestIterator dul, DestAccessor dest,
515  Kernel const & kernel)
516 {
517  vigra_precondition(src.size(sul) == 2,
518  "orientedTrigonometricFilter(): input image must have 2 bands.");
519  vigra_precondition(dest.size(dul) == 3,
520  "orientedTrigonometricFilter(): output image must have 3 bands.");
521 
522  int w = slr.x - sul.x;
523  int h = slr.y - sul.y;
524  int radius = kernel.radius_;
525 
526  typedef typename SrcAccessor::value_type VectorType;
527  typedef typename DestAccessor::value_type TensorType;
528 
529  initImage(dul, dul+Diff2D(w,h), dest, NumericTraits<TensorType>::zero());
530 
531  for(int y=0; y<h; ++y, ++sul.y, ++dul.y)
532  {
533  SrcIterator s = sul;
534  DestIterator d = dul;
535  for(int x=0; x<w; ++x, ++s.x, ++d.x)
536  {
537  int x0 = x - radius < 0 ? -x : -radius;
538  int y0 = y - radius < 0 ? -y : -radius;
539  int x1 = x + radius >= w ? w - x - 1 : radius;
540  int y1 = y + radius >= h ? h - y - 1 : radius;
541 
542  VectorType v(src(s));
543  TensorType t(sq(v[0]), v[0]*v[1], sq(v[1]));
544  double sqMag = t[0] + t[2];
545  double mag = VIGRA_CSTD::sqrt(sqMag);
546  if(mag != 0.0)
547  v /= mag;
548  else
549  v *= 0.0;
550  Diff2D dd;
551  for(dd.y = y0; dd.y <= y1; ++dd.y)
552  {
553  for(dd.x = x0; dd.x <= x1; ++dd.x)
554  {
555  dest.set(dest(d, dd) + kernel(dd.x, dd.y, v) * t, d, dd);
556  }
557  }
558  }
559  }
560 }
561 
562 template <class SrcIterator, class SrcAccessor,
563  class DestIterator, class DestAccessor,
564  class Kernel>
565 inline
566 void orientedTrigonometricFilter(triple<SrcIterator, SrcIterator, SrcAccessor> s,
567  pair<DestIterator, DestAccessor> d,
568  Kernel const & kernel)
569 {
570  orientedTrigonometricFilter(s.first, s.second, s.third, d.first, d.second, kernel);
571 }
572 
573 //@}
574 
575 } // namespace vigra
576 
577 #endif /* VIGRA_ORIENTEDTENSORFILTERS_HXX */

© 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 (Sat Oct 5 2013)