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

edgedetection.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2002 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_EDGEDETECTION_HXX
38 #define VIGRA_EDGEDETECTION_HXX
39 
40 #include <vector>
41 #include <queue>
42 #include <cmath> // sqrt(), abs()
43 #include "utilities.hxx"
44 #include "numerictraits.hxx"
45 #include "stdimage.hxx"
46 #include "stdimagefunctions.hxx"
47 #include "recursiveconvolution.hxx"
48 #include "separableconvolution.hxx"
49 #include "convolution.hxx"
50 #include "labelimage.hxx"
51 #include "mathutil.hxx"
52 #include "pixelneighborhood.hxx"
53 #include "linear_solve.hxx"
54 #include "functorexpression.hxx"
55 
56 
57 namespace vigra {
58 
59 /** \addtogroup EdgeDetection Edge Detection
60  Edge detectors based on first and second derivatives,
61  and related post-processing.
62 */
63 //@{
64 
65 /********************************************************/
66 /* */
67 /* differenceOfExponentialEdgeImage */
68 /* */
69 /********************************************************/
70 
71 /** \brief Detect and mark edges in an edge image using the Shen/Castan zero-crossing detector.
72 
73  This operator applies an exponential filter to the source image
74  at the given <TT>scale</TT> and subtracts the result from the original image.
75  Zero crossings are detected in the resulting difference image. Whenever the
76  gradient at a zero crossing is greater than the given <TT>gradient_threshold</TT>,
77  an edge point is marked (using <TT>edge_marker</TT>) in the destination image on
78  the darker side of the zero crossing (note that zero crossings occur
79  <i>between</i> pixels). For example:
80 
81  \code
82  sign of difference image resulting edge points (*)
83 
84  + - - * * .
85  + + - => . * *
86  + + + . . .
87  \endcode
88 
89  Non-edge pixels (<TT>.</TT>) remain untouched in the destination image.
90  The result can be improved by the post-processing operation \ref removeShortEdges().
91  A more accurate edge placement can be achieved with the function
92  \ref differenceOfExponentialCrackEdgeImage().
93 
94  The source value type
95  (<TT>SrcAccessor::value_type</TT>) must be a linear algebra, i.e. addition,
96  subtraction and multiplication of the type with itself, and multiplication
97  with double and
98  \ref NumericTraits "NumericTraits" must
99  be defined. In addition, this type must be less-comparable.
100 
101  <b> Declarations:</b>
102 
103  pass arguments explicitly:
104  \code
105  namespace vigra {
106  template <class SrcIterator, class SrcAccessor,
107  class DestIterator, class DestAccessor,
108  class GradValue,
109  class DestValue = DestAccessor::value_type>
110  void differenceOfExponentialEdgeImage(
111  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
112  DestIterator dul, DestAccessor da,
113  double scale, GradValue gradient_threshold,
114  DestValue edge_marker = NumericTraits<DestValue>::one())
115  }
116  \endcode
117 
118  use argument objects in conjunction with \ref ArgumentObjectFactories :
119  \code
120  namespace vigra {
121  template <class SrcIterator, class SrcAccessor,
122  class DestIterator, class DestAccessor,
123  class GradValue,
124  class DestValue = DestAccessor::value_type>
125  void differenceOfExponentialEdgeImage(
126  triple<SrcIterator, SrcIterator, SrcAccessor> src,
127  pair<DestIterator, DestAccessor> dest,
128  double scale, GradValue gradient_threshold,
129  DestValue edge_marker = NumericTraits<DestValue>::one())
130  }
131  \endcode
132 
133  <b> Usage:</b>
134 
135  <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
136  Namespace: vigra
137 
138  \code
139  vigra::BImage src(w,h), edges(w,h);
140 
141  // empty edge image
142  edges = 0;
143  ...
144 
145  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
146  vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges),
147  0.8, 4.0, 1);
148  \endcode
149 
150  <b> Required Interface:</b>
151 
152  \code
153  SrcImageIterator src_upperleft, src_lowerright;
154  DestImageIterator dest_upperleft;
155 
156  SrcAccessor src_accessor;
157  DestAccessor dest_accessor;
158 
159  SrcAccessor::value_type u = src_accessor(src_upperleft);
160  double d;
161  GradValue gradient_threshold;
162 
163  u = u + u
164  u = u - u
165  u = u * u
166  u = d * u
167  u < gradient_threshold
168 
169  DestValue edge_marker;
170  dest_accessor.set(edge_marker, dest_upperleft);
171  \endcode
172 
173  <b> Preconditions:</b>
174 
175  \code
176  scale > 0
177  gradient_threshold > 0
178  \endcode
179 */
181 
182 template <class SrcIterator, class SrcAccessor,
183  class DestIterator, class DestAccessor,
184  class GradValue, class DestValue>
186  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
187  DestIterator dul, DestAccessor da,
188  double scale, GradValue gradient_threshold, DestValue edge_marker)
189 {
190  vigra_precondition(scale > 0,
191  "differenceOfExponentialEdgeImage(): scale > 0 required.");
192 
193  vigra_precondition(gradient_threshold > 0,
194  "differenceOfExponentialEdgeImage(): "
195  "gradient_threshold > 0 required.");
196 
197  int w = slr.x - sul.x;
198  int h = slr.y - sul.y;
199  int x,y;
200 
201  typedef typename
202  NumericTraits<typename SrcAccessor::value_type>::RealPromote
203  TMPTYPE;
204  typedef BasicImage<TMPTYPE> TMPIMG;
205 
206  TMPIMG tmp(w,h);
207  TMPIMG smooth(w,h);
208 
209  recursiveSmoothX(srcIterRange(sul, slr, sa), destImage(tmp), scale / 2.0);
210  recursiveSmoothY(srcImageRange(tmp), destImage(tmp), scale / 2.0);
211 
212  recursiveSmoothX(srcImageRange(tmp), destImage(smooth), scale);
213  recursiveSmoothY(srcImageRange(smooth), destImage(smooth), scale);
214 
215  typename TMPIMG::Iterator iy = smooth.upperLeft();
216  typename TMPIMG::Iterator ty = tmp.upperLeft();
217  DestIterator dy = dul;
218 
219  static const Diff2D right(1, 0);
220  static const Diff2D bottom(0, 1);
221 
222 
223  TMPTYPE thresh = detail::RequiresExplicitCast<TMPTYPE>::cast((gradient_threshold * gradient_threshold) *
224  NumericTraits<TMPTYPE>::one());
225  TMPTYPE zero = NumericTraits<TMPTYPE>::zero();
226 
227  for(y=0; y<h-1; ++y, ++iy.y, ++ty.y, ++dy.y)
228  {
229  typename TMPIMG::Iterator ix = iy;
230  typename TMPIMG::Iterator tx = ty;
231  DestIterator dx = dy;
232 
233  for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, ++dx.x)
234  {
235  TMPTYPE diff = *tx - *ix;
236  TMPTYPE gx = tx[right] - *tx;
237  TMPTYPE gy = tx[bottom] - *tx;
238 
239  if((gx * gx > thresh) &&
240  (diff * (tx[right] - ix[right]) < zero))
241  {
242  if(gx < zero)
243  {
244  da.set(edge_marker, dx, right);
245  }
246  else
247  {
248  da.set(edge_marker, dx);
249  }
250  }
251  if(((gy * gy > thresh) &&
252  (diff * (tx[bottom] - ix[bottom]) < zero)))
253  {
254  if(gy < zero)
255  {
256  da.set(edge_marker, dx, bottom);
257  }
258  else
259  {
260  da.set(edge_marker, dx);
261  }
262  }
263  }
264  }
265 
266  typename TMPIMG::Iterator ix = iy;
267  typename TMPIMG::Iterator tx = ty;
268  DestIterator dx = dy;
269 
270  for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, ++dx.x)
271  {
272  TMPTYPE diff = *tx - *ix;
273  TMPTYPE gx = tx[right] - *tx;
274 
275  if((gx * gx > thresh) &&
276  (diff * (tx[right] - ix[right]) < zero))
277  {
278  if(gx < zero)
279  {
280  da.set(edge_marker, dx, right);
281  }
282  else
283  {
284  da.set(edge_marker, dx);
285  }
286  }
287  }
288 }
289 
290 template <class SrcIterator, class SrcAccessor,
291  class DestIterator, class DestAccessor,
292  class GradValue>
293 inline
295  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
296  DestIterator dul, DestAccessor da,
297  double scale, GradValue gradient_threshold)
298 {
299  differenceOfExponentialEdgeImage(sul, slr, sa, dul, da,
300  scale, gradient_threshold, 1);
301 }
302 
303 template <class SrcIterator, class SrcAccessor,
304  class DestIterator, class DestAccessor,
305  class GradValue, class DestValue>
306 inline
308  triple<SrcIterator, SrcIterator, SrcAccessor> src,
309  pair<DestIterator, DestAccessor> dest,
310  double scale, GradValue gradient_threshold,
311  DestValue edge_marker)
312 {
313  differenceOfExponentialEdgeImage(src.first, src.second, src.third,
314  dest.first, dest.second,
315  scale, gradient_threshold,
316  edge_marker);
317 }
318 
319 template <class SrcIterator, class SrcAccessor,
320  class DestIterator, class DestAccessor,
321  class GradValue>
322 inline
324  triple<SrcIterator, SrcIterator, SrcAccessor> src,
325  pair<DestIterator, DestAccessor> dest,
326  double scale, GradValue gradient_threshold)
327 {
328  differenceOfExponentialEdgeImage(src.first, src.second, src.third,
329  dest.first, dest.second,
330  scale, gradient_threshold, 1);
331 }
332 
333 /********************************************************/
334 /* */
335 /* differenceOfExponentialCrackEdgeImage */
336 /* */
337 /********************************************************/
338 
339 /** \brief Detect and mark edges in a crack edge image using the Shen/Castan zero-crossing detector.
340 
341  This operator applies an exponential filter to the source image
342  at the given <TT>scale</TT> and subtracts the result from the original image.
343  Zero crossings are detected in the resulting difference image. Whenever the
344  gradient at a zero crossing is greater than the given <TT>gradient_threshold</TT>,
345  an edge point is marked (using <TT>edge_marker</TT>) in the destination image
346  <i>between</i> the corresponding original pixels. Topologically, this means we
347  must insert additional pixels between the original ones to represent the
348  boundaries between the pixels (the so called zero- and one-cells, with the original
349  pixels being two-cells). Within VIGRA, such an image is called \ref CrackEdgeImage.
350  To allow insertion of the zero- and one-cells, the destination image must have twice the
351  size of the original (precisely, <TT>(2*w-1)</TT> by <TT>(2*h-1)</TT> pixels). Then the algorithm
352  proceeds as follows:
353 
354  \code
355 sign of difference image insert zero- and one-cells resulting edge points (*)
356 
357  + . - . - . * . . .
358  + - - . . . . . . * * * .
359  + + - => + . + . - => . . . * .
360  + + + . . . . . . . . * *
361  + . + . + . . . . .
362  \endcode
363 
364  Thus the edge points are marked where they actually are - in between the pixels.
365  An important property of the resulting edge image is that it conforms to the notion
366  of well-composedness as defined by Latecki et al., i.e. connected regions and edges
367  obtained by a subsequent \ref Labeling do not depend on
368  whether 4- or 8-connectivity is used.
369  The non-edge pixels (<TT>.</TT>) in the destination image remain unchanged.
370  The result conformes to the requirements of a \ref CrackEdgeImage. It can be further
371  improved by the post-processing operations \ref removeShortEdges() and
372  \ref closeGapsInCrackEdgeImage().
373 
374  The source value type (<TT>SrcAccessor::value_type</TT>) must be a linear algebra, i.e. addition,
375  subtraction and multiplication of the type with itself, and multiplication
376  with double and
377  \ref NumericTraits "NumericTraits" must
378  be defined. In addition, this type must be less-comparable.
379 
380  <b> Declarations:</b>
381 
382  pass arguments explicitly:
383  \code
384  namespace vigra {
385  template <class SrcIterator, class SrcAccessor,
386  class DestIterator, class DestAccessor,
387  class GradValue,
388  class DestValue = DestAccessor::value_type>
389  void differenceOfExponentialCrackEdgeImage(
390  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
391  DestIterator dul, DestAccessor da,
392  double scale, GradValue gradient_threshold,
393  DestValue edge_marker = NumericTraits<DestValue>::one())
394  }
395  \endcode
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  class GradValue,
403  class DestValue = DestAccessor::value_type>
404  void differenceOfExponentialCrackEdgeImage(
405  triple<SrcIterator, SrcIterator, SrcAccessor> src,
406  pair<DestIterator, DestAccessor> dest,
407  double scale, GradValue gradient_threshold,
408  DestValue edge_marker = NumericTraits<DestValue>::one())
409  }
410  \endcode
411 
412  <b> Usage:</b>
413 
414  <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
415  Namespace: vigra
416 
417  \code
418  vigra::BImage src(w,h), edges(2*w-1,2*h-1);
419 
420  // empty edge image
421  edges = 0;
422  ...
423 
424  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
425  vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges),
426  0.8, 4.0, 1);
427  \endcode
428 
429  <b> Required Interface:</b>
430 
431  \code
432  SrcImageIterator src_upperleft, src_lowerright;
433  DestImageIterator dest_upperleft;
434 
435  SrcAccessor src_accessor;
436  DestAccessor dest_accessor;
437 
438  SrcAccessor::value_type u = src_accessor(src_upperleft);
439  double d;
440  GradValue gradient_threshold;
441 
442  u = u + u
443  u = u - u
444  u = u * u
445  u = d * u
446  u < gradient_threshold
447 
448  DestValue edge_marker;
449  dest_accessor.set(edge_marker, dest_upperleft);
450  \endcode
451 
452  <b> Preconditions:</b>
453 
454  \code
455  scale > 0
456  gradient_threshold > 0
457  \endcode
458 
459  The destination image must have twice the size of the source:
460  \code
461  w_dest = 2 * w_src - 1
462  h_dest = 2 * h_src - 1
463  \endcode
464 */
466 
467 template <class SrcIterator, class SrcAccessor,
468  class DestIterator, class DestAccessor,
469  class GradValue, class DestValue>
471  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
472  DestIterator dul, DestAccessor da,
473  double scale, GradValue gradient_threshold,
474  DestValue edge_marker)
475 {
476  vigra_precondition(scale > 0,
477  "differenceOfExponentialCrackEdgeImage(): scale > 0 required.");
478 
479  vigra_precondition(gradient_threshold > 0,
480  "differenceOfExponentialCrackEdgeImage(): "
481  "gradient_threshold > 0 required.");
482 
483  int w = slr.x - sul.x;
484  int h = slr.y - sul.y;
485  int x, y;
486 
487  typedef typename
488  NumericTraits<typename SrcAccessor::value_type>::RealPromote
489  TMPTYPE;
490  typedef BasicImage<TMPTYPE> TMPIMG;
491 
492  TMPIMG tmp(w,h);
493  TMPIMG smooth(w,h);
494 
495  TMPTYPE zero = NumericTraits<TMPTYPE>::zero();
496 
497  static const Diff2D right(1,0);
498  static const Diff2D bottom(0,1);
499  static const Diff2D left(-1,0);
500  static const Diff2D top(0,-1);
501 
502  recursiveSmoothX(srcIterRange(sul, slr, sa), destImage(tmp), scale / 2.0);
503  recursiveSmoothY(srcImageRange(tmp), destImage(tmp), scale / 2.0);
504 
505  recursiveSmoothX(srcImageRange(tmp), destImage(smooth), scale);
506  recursiveSmoothY(srcImageRange(smooth), destImage(smooth), scale);
507 
508  typename TMPIMG::Iterator iy = smooth.upperLeft();
509  typename TMPIMG::Iterator ty = tmp.upperLeft();
510  DestIterator dy = dul;
511 
512  TMPTYPE thresh = detail::RequiresExplicitCast<TMPTYPE>::cast((gradient_threshold * gradient_threshold) *
513  NumericTraits<TMPTYPE>::one());
514 
515  // find zero crossings above threshold
516  for(y=0; y<h-1; ++y, ++iy.y, ++ty.y, dy.y+=2)
517  {
518  typename TMPIMG::Iterator ix = iy;
519  typename TMPIMG::Iterator tx = ty;
520  DestIterator dx = dy;
521 
522  for(int x=0; x<w-1; ++x, ++ix.x, ++tx.x, dx.x+=2)
523  {
524  TMPTYPE diff = *tx - *ix;
525  TMPTYPE gx = tx[right] - *tx;
526  TMPTYPE gy = tx[bottom] - *tx;
527 
528  if((gx * gx > thresh) &&
529  (diff * (tx[right] - ix[right]) < zero))
530  {
531  da.set(edge_marker, dx, right);
532  }
533  if((gy * gy > thresh) &&
534  (diff * (tx[bottom] - ix[bottom]) < zero))
535  {
536  da.set(edge_marker, dx, bottom);
537  }
538  }
539 
540  TMPTYPE diff = *tx - *ix;
541  TMPTYPE gy = tx[bottom] - *tx;
542 
543  if((gy * gy > thresh) &&
544  (diff * (tx[bottom] - ix[bottom]) < zero))
545  {
546  da.set(edge_marker, dx, bottom);
547  }
548  }
549 
550  typename TMPIMG::Iterator ix = iy;
551  typename TMPIMG::Iterator tx = ty;
552  DestIterator dx = dy;
553 
554  for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, dx.x+=2)
555  {
556  TMPTYPE diff = *tx - *ix;
557  TMPTYPE gx = tx[right] - *tx;
558 
559  if((gx * gx > thresh) &&
560  (diff * (tx[right] - ix[right]) < zero))
561  {
562  da.set(edge_marker, dx, right);
563  }
564  }
565 
566  iy = smooth.upperLeft() + Diff2D(0,1);
567  ty = tmp.upperLeft() + Diff2D(0,1);
568  dy = dul + Diff2D(1,2);
569 
570  static const Diff2D topleft(-1,-1);
571  static const Diff2D topright(1,-1);
572  static const Diff2D bottomleft(-1,1);
573  static const Diff2D bottomright(1,1);
574 
575  // find missing 1-cells below threshold (x-direction)
576  for(y=0; y<h-2; ++y, ++iy.y, ++ty.y, dy.y+=2)
577  {
578  typename TMPIMG::Iterator ix = iy;
579  typename TMPIMG::Iterator tx = ty;
580  DestIterator dx = dy;
581 
582  for(int x=0; x<w-2; ++x, ++ix.x, ++tx.x, dx.x+=2)
583  {
584  if(da(dx) == edge_marker) continue;
585 
586  TMPTYPE diff = *tx - *ix;
587 
588  if((diff * (tx[right] - ix[right]) < zero) &&
589  (((da(dx, bottomright) == edge_marker) &&
590  (da(dx, topleft) == edge_marker)) ||
591  ((da(dx, bottomleft) == edge_marker) &&
592  (da(dx, topright) == edge_marker))))
593 
594  {
595  da.set(edge_marker, dx);
596  }
597  }
598  }
599 
600  iy = smooth.upperLeft() + Diff2D(1,0);
601  ty = tmp.upperLeft() + Diff2D(1,0);
602  dy = dul + Diff2D(2,1);
603 
604  // find missing 1-cells below threshold (y-direction)
605  for(y=0; y<h-2; ++y, ++iy.y, ++ty.y, dy.y+=2)
606  {
607  typename TMPIMG::Iterator ix = iy;
608  typename TMPIMG::Iterator tx = ty;
609  DestIterator dx = dy;
610 
611  for(int x=0; x<w-2; ++x, ++ix.x, ++tx.x, dx.x+=2)
612  {
613  if(da(dx) == edge_marker) continue;
614 
615  TMPTYPE diff = *tx - *ix;
616 
617  if((diff * (tx[bottom] - ix[bottom]) < zero) &&
618  (((da(dx, bottomright) == edge_marker) &&
619  (da(dx, topleft) == edge_marker)) ||
620  ((da(dx, bottomleft) == edge_marker) &&
621  (da(dx, topright) == edge_marker))))
622 
623  {
624  da.set(edge_marker, dx);
625  }
626  }
627  }
628 
629  dy = dul + Diff2D(1,1);
630 
631  // find missing 0-cells
632  for(y=0; y<h-1; ++y, dy.y+=2)
633  {
634  DestIterator dx = dy;
635 
636  for(int x=0; x<w-1; ++x, dx.x+=2)
637  {
638  static const Diff2D dist[] = {right, top, left, bottom };
639 
640  int i;
641  for(i=0; i<4; ++i)
642  {
643  if(da(dx, dist[i]) == edge_marker) break;
644  }
645 
646  if(i < 4) da.set(edge_marker, dx);
647  }
648  }
649 }
650 
651 template <class SrcIterator, class SrcAccessor,
652  class DestIterator, class DestAccessor,
653  class GradValue, class DestValue>
654 inline
656  triple<SrcIterator, SrcIterator, SrcAccessor> src,
657  pair<DestIterator, DestAccessor> dest,
658  double scale, GradValue gradient_threshold,
659  DestValue edge_marker)
660 {
661  differenceOfExponentialCrackEdgeImage(src.first, src.second, src.third,
662  dest.first, dest.second,
663  scale, gradient_threshold,
664  edge_marker);
665 }
666 
667 /********************************************************/
668 /* */
669 /* removeShortEdges */
670 /* */
671 /********************************************************/
672 
673 /** \brief Remove short edges from an edge image.
674 
675  This algorithm can be applied as a post-processing operation of
676  \ref differenceOfExponentialEdgeImage() and \ref differenceOfExponentialCrackEdgeImage().
677  It removes all edges that are shorter than <TT>min_edge_length</TT>. The corresponding
678  pixels are set to the <TT>non_edge_marker</TT>. The idea behind this algorithms is
679  that very short edges are probably caused by noise and don't represent interesting
680  image structure. Technically, the algorithms executes a connected components labeling,
681  so the image's value type must be equality comparable.
682 
683  If the source image fulfills the requirements of a \ref CrackEdgeImage,
684  it will still do so after application of this algorithm.
685 
686  Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place,
687  i.e. on only one image. Also, the algorithm assumes that all non-edges pixels are already
688  marked with the given <TT>non_edge_marker</TT> value.
689 
690  <b> Declarations:</b>
691 
692  pass arguments explicitly:
693  \code
694  namespace vigra {
695  template <class Iterator, class Accessor, class SrcValue>
696  void removeShortEdges(
697  Iterator sul, Iterator slr, Accessor sa,
698  int min_edge_length, SrcValue non_edge_marker)
699  }
700  \endcode
701 
702  use argument objects in conjunction with \ref ArgumentObjectFactories :
703  \code
704  namespace vigra {
705  template <class Iterator, class Accessor, class SrcValue>
706  void removeShortEdges(
707  triple<Iterator, Iterator, Accessor> src,
708  int min_edge_length, SrcValue non_edge_marker)
709  }
710  \endcode
711 
712  <b> Usage:</b>
713 
714  <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
715  Namespace: vigra
716 
717  \code
718  vigra::BImage src(w,h), edges(w,h);
719 
720  // empty edge image
721  edges = 0;
722  ...
723 
724  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
725  vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges),
726  0.8, 4.0, 1);
727 
728  // zero edges shorter than 10 pixels
729  vigra::removeShortEdges(srcImageRange(edges), 10, 0);
730  \endcode
731 
732  <b> Required Interface:</b>
733 
734  \code
735  SrcImageIterator src_upperleft, src_lowerright;
736  DestImageIterator dest_upperleft;
737 
738  SrcAccessor src_accessor;
739  DestAccessor dest_accessor;
740 
741  SrcAccessor::value_type u = src_accessor(src_upperleft);
742 
743  u == u
744 
745  SrcValue non_edge_marker;
746  src_accessor.set(non_edge_marker, src_upperleft);
747  \endcode
748 */
750 
751 template <class Iterator, class Accessor, class Value>
752 void removeShortEdges(
753  Iterator sul, Iterator slr, Accessor sa,
754  unsigned int min_edge_length, Value non_edge_marker)
755 {
756  int w = slr.x - sul.x;
757  int h = slr.y - sul.y;
758  int x,y;
759 
760  IImage labels(w, h);
761  labels = 0;
762 
763  int number_of_regions =
764  labelImageWithBackground(srcIterRange(sul,slr,sa),
765  destImage(labels), true, non_edge_marker);
766 
767  ArrayOfRegionStatistics<FindROISize<int> >
768  region_stats(number_of_regions);
769 
770  inspectTwoImages(srcImageRange(labels), srcImage(labels), region_stats);
771 
772  IImage::Iterator ly = labels.upperLeft();
773  Iterator oy = sul;
774 
775  for(y=0; y<h; ++y, ++oy.y, ++ly.y)
776  {
777  Iterator ox(oy);
778  IImage::Iterator lx(ly);
779 
780  for(x=0; x<w; ++x, ++ox.x, ++lx.x)
781  {
782  if(sa(ox) == non_edge_marker) continue;
783  if((region_stats[*lx].count) < min_edge_length)
784  {
785  sa.set(non_edge_marker, ox);
786  }
787  }
788  }
789 }
790 
791 template <class Iterator, class Accessor, class Value>
792 inline
793 void removeShortEdges(
794  triple<Iterator, Iterator, Accessor> src,
795  unsigned int min_edge_length, Value non_edge_marker)
796 {
797  removeShortEdges(src.first, src.second, src.third,
798  min_edge_length, non_edge_marker);
799 }
800 
801 /********************************************************/
802 /* */
803 /* closeGapsInCrackEdgeImage */
804 /* */
805 /********************************************************/
806 
807 /** \brief Close one-pixel wide gaps in a cell grid edge image.
808 
809  This algorithm is typically applied as a post-processing operation of
810  \ref differenceOfExponentialCrackEdgeImage(). The source image must fulfill
811  the requirements of a \ref CrackEdgeImage, and will still do so after
812  application of this algorithm.
813 
814  It closes one pixel wide gaps in the edges resulting from this algorithm.
815  Since these gaps are usually caused by zero crossing slightly below the gradient
816  threshold used in edge detection, this algorithms acts like a weak hysteresis
817  thresholding. The newly found edge pixels are marked with the given <TT>edge_marker</TT>.
818  The image's value type must be equality comparable.
819 
820  Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place,
821  i.e. on only one image.
822 
823  <b> Declarations:</b>
824 
825  pass arguments explicitly:
826  \code
827  namespace vigra {
828  template <class SrcIterator, class SrcAccessor, class SrcValue>
829  void closeGapsInCrackEdgeImage(
830  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
831  SrcValue edge_marker)
832  }
833  \endcode
834 
835  use argument objects in conjunction with \ref ArgumentObjectFactories :
836  \code
837  namespace vigra {
838  template <class SrcIterator, class SrcAccessor, class SrcValue>
839  void closeGapsInCrackEdgeImage(
840  triple<SrcIterator, SrcIterator, SrcAccessor> src,
841  SrcValue edge_marker)
842  }
843  \endcode
844 
845  <b> Usage:</b>
846 
847  <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
848  Namespace: vigra
849 
850  \code
851  vigra::BImage src(w,h), edges(2*w-1, 2*h-1);
852 
853  // empty edge image
854  edges = 0;
855  ...
856 
857  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
858  vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges),
859  0.8, 4.0, 1);
860 
861  // close gaps, mark with 1
862  vigra::closeGapsInCrackEdgeImage(srcImageRange(edges), 1);
863 
864  // zero edges shorter than 20 pixels
865  vigra::removeShortEdges(srcImageRange(edges), 10, 0);
866  \endcode
867 
868  <b> Required Interface:</b>
869 
870  \code
871  SrcImageIterator src_upperleft, src_lowerright;
872 
873  SrcAccessor src_accessor;
874  DestAccessor dest_accessor;
875 
876  SrcAccessor::value_type u = src_accessor(src_upperleft);
877 
878  u == u
879  u != u
880 
881  SrcValue edge_marker;
882  src_accessor.set(edge_marker, src_upperleft);
883  \endcode
884 */
886 
887 template <class SrcIterator, class SrcAccessor, class SrcValue>
889  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
890  SrcValue edge_marker)
891 {
892  int w = slr.x - sul.x;
893  int h = slr.y - sul.y;
894 
895  vigra_precondition(w % 2 == 1 && h % 2 == 1,
896  "closeGapsInCrackEdgeImage(): Input is not a crack edge image (must have odd-numbered shape).");
897 
898  int w2 = w / 2, h2 = h / 2, x, y;
899 
900  int count1, count2, count3;
901 
902  static const Diff2D right(1,0);
903  static const Diff2D bottom(0,1);
904  static const Diff2D left(-1,0);
905  static const Diff2D top(0,-1);
906 
907  static const Diff2D leftdist[] = {
908  Diff2D(0, 0), Diff2D(-1, 1), Diff2D(-2, 0), Diff2D(-1, -1)};
909  static const Diff2D rightdist[] = {
910  Diff2D(2, 0), Diff2D(1, 1), Diff2D(0, 0), Diff2D(1, -1)};
911  static const Diff2D topdist[] = {
912  Diff2D(1, -1), Diff2D(0, 0), Diff2D(-1, -1), Diff2D(0, -2)};
913  static const Diff2D bottomdist[] = {
914  Diff2D(1, 1), Diff2D(0, 2), Diff2D(-1, 1), Diff2D(0, 0)};
915 
916  int i;
917 
918  SrcIterator sy = sul + Diff2D(0,1);
919  SrcIterator sx;
920 
921  // close 1-pixel wide gaps (x-direction)
922  for(y=0; y<h2; ++y, sy.y+=2)
923  {
924  sx = sy + Diff2D(2,0);
925 
926  for(x=2; x<w2; ++x, sx.x+=2)
927  {
928  if(sa(sx) == edge_marker) continue;
929 
930  if(sa(sx, left) != edge_marker) continue;
931  if(sa(sx, right) != edge_marker) continue;
932 
933  count1 = 0;
934  count2 = 0;
935  count3 = 0;
936 
937  for(i=0; i<4; ++i)
938  {
939  if(sa(sx, leftdist[i]) == edge_marker)
940  {
941  ++count1;
942  count3 ^= 1 << i;
943  }
944  if(sa(sx, rightdist[i]) == edge_marker)
945  {
946  ++count2;
947  count3 ^= 1 << i;
948  }
949  }
950 
951  if(count1 <= 1 || count2 <= 1 || count3 == 15)
952  {
953  sa.set(edge_marker, sx);
954  }
955  }
956  }
957 
958  sy = sul + Diff2D(1,2);
959 
960  // close 1-pixel wide gaps (y-direction)
961  for(y=2; y<h2; ++y, sy.y+=2)
962  {
963  sx = sy;
964 
965  for(x=0; x<w2; ++x, sx.x+=2)
966  {
967  if(sa(sx) == edge_marker) continue;
968 
969  if(sa(sx, top) != edge_marker) continue;
970  if(sa(sx, bottom) != edge_marker) continue;
971 
972  count1 = 0;
973  count2 = 0;
974  count3 = 0;
975 
976  for(i=0; i<4; ++i)
977  {
978  if(sa(sx, topdist[i]) == edge_marker)
979  {
980  ++count1;
981  count3 ^= 1 << i;
982  }
983  if(sa(sx, bottomdist[i]) == edge_marker)
984  {
985  ++count2;
986  count3 ^= 1 << i;
987  }
988  }
989 
990  if(count1 <= 1 || count2 <= 1 || count3 == 15)
991  {
992  sa.set(edge_marker, sx);
993  }
994  }
995  }
996 }
997 
998 template <class SrcIterator, class SrcAccessor, class SrcValue>
999 inline
1001  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1002  SrcValue edge_marker)
1003 {
1004  closeGapsInCrackEdgeImage(src.first, src.second, src.third,
1005  edge_marker);
1006 }
1007 
1008 /********************************************************/
1009 /* */
1010 /* beautifyCrackEdgeImage */
1011 /* */
1012 /********************************************************/
1013 
1014 /** \brief Beautify crack edge image for visualization.
1015 
1016  This algorithm is applied as a post-processing operation of
1017  \ref differenceOfExponentialCrackEdgeImage(). The source image must fulfill
1018  the requirements of a \ref CrackEdgeImage, but will <b> not</b> do so after
1019  application of this algorithm. In particular, the algorithm removes zero-cells
1020  marked as edges to avoid staircase effects on diagonal lines like this:
1021 
1022  \code
1023  original edge points (*) resulting edge points
1024 
1025  . * . . . . * . . .
1026  . * * * . . . * . .
1027  . . . * . => . . . * .
1028  . . . * * . . . . *
1029  . . . . . . . . . .
1030  \endcode
1031 
1032  Therfore, this algorithm should only be applied as a vizualization aid, i.e.
1033  for human inspection. The algorithm assumes that edges are marked with <TT>edge_marker</TT>,
1034  and background pixels with <TT>background_marker</TT>. The image's value type must be
1035  equality comparable.
1036 
1037  Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place,
1038  i.e. on only one image.
1039 
1040  <b> Declarations:</b>
1041 
1042  pass arguments explicitly:
1043  \code
1044  namespace vigra {
1045  template <class SrcIterator, class SrcAccessor, class SrcValue>
1046  void beautifyCrackEdgeImage(
1047  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1048  SrcValue edge_marker, SrcValue background_marker)
1049  }
1050  \endcode
1051 
1052  use argument objects in conjunction with \ref ArgumentObjectFactories :
1053  \code
1054  namespace vigra {
1055  template <class SrcIterator, class SrcAccessor, class SrcValue>
1056  void beautifyCrackEdgeImage(
1057  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1058  SrcValue edge_marker, SrcValue background_marker)
1059  }
1060  \endcode
1061 
1062  <b> Usage:</b>
1063 
1064  <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
1065  Namespace: vigra
1066 
1067  \code
1068  vigra::BImage src(w,h), edges(2*w-1, 2*h-1);
1069 
1070  // empty edge image
1071  edges = 0;
1072  ...
1073 
1074  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
1075  vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges),
1076  0.8, 4.0, 1);
1077 
1078  // beautify edge image for visualization
1079  vigra::beautifyCrackEdgeImage(destImageRange(edges), 1, 0);
1080 
1081  // show to the user
1082  window.open(edges);
1083  \endcode
1084 
1085  <b> Required Interface:</b>
1086 
1087  \code
1088  SrcImageIterator src_upperleft, src_lowerright;
1089 
1090  SrcAccessor src_accessor;
1091  DestAccessor dest_accessor;
1092 
1093  SrcAccessor::value_type u = src_accessor(src_upperleft);
1094 
1095  u == u
1096  u != u
1097 
1098  SrcValue background_marker;
1099  src_accessor.set(background_marker, src_upperleft);
1100  \endcode
1101 */
1103 
1104 template <class SrcIterator, class SrcAccessor, class SrcValue>
1106  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1107  SrcValue edge_marker, SrcValue background_marker)
1108 {
1109  int w = slr.x - sul.x;
1110  int h = slr.y - sul.y;
1111 
1112  vigra_precondition(w % 2 == 1 && h % 2 == 1,
1113  "beautifyCrackEdgeImage(): Input is not a crack edge image (must have odd-numbered shape).");
1114 
1115  int w2 = w / 2, h2 = h / 2, x, y;
1116 
1117  SrcIterator sy = sul + Diff2D(1,1);
1118  SrcIterator sx;
1119 
1120  static const Diff2D right(1,0);
1121  static const Diff2D bottom(0,1);
1122  static const Diff2D left(-1,0);
1123  static const Diff2D top(0,-1);
1124 
1125  // delete 0-cells at corners
1126  for(y=0; y<h2; ++y, sy.y+=2)
1127  {
1128  sx = sy;
1129 
1130  for(x=0; x<w2; ++x, sx.x+=2)
1131  {
1132  if(sa(sx) != edge_marker) continue;
1133 
1134  if(sa(sx, right) == edge_marker && sa(sx, left) == edge_marker) continue;
1135  if(sa(sx, bottom) == edge_marker && sa(sx, top) == edge_marker) continue;
1136 
1137  sa.set(background_marker, sx);
1138  }
1139  }
1140 }
1141 
1142 template <class SrcIterator, class SrcAccessor, class SrcValue>
1143 inline
1145  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1146  SrcValue edge_marker, SrcValue background_marker)
1147 {
1148  beautifyCrackEdgeImage(src.first, src.second, src.third,
1149  edge_marker, background_marker);
1150 }
1151 
1152 
1153 /** Helper class that stores edgel attributes.
1154 */
1155 class Edgel
1156 {
1157  public:
1158 
1159  /** The type of an Edgel's members.
1160  */
1161  typedef float value_type;
1162 
1163  /** The edgel's sub-pixel x coordinate.
1164  */
1166 
1167  /** The edgel's sub-pixel y coordinate.
1168  */
1170 
1171  /** The edgel's strength (magnitude of the gradient vector).
1172  */
1174 
1175  /**
1176  The edgel's orientation. This is the clockwise angle in radians
1177  between the x-axis and the edge, so that the bright side of the
1178  edge is on the left when one looks along the orientation vector.
1179  The angle is measured clockwise because the y-axis increases
1180  downwards (left-handed coordinate system):
1181 
1182  \code
1183 
1184  edgel axis
1185  \
1186  (dark \ (bright side)
1187  side) \
1188  \
1189  +------------> x-axis
1190  |\ |
1191  | \ /_/ orientation angle
1192  | \\
1193  | \
1194  |
1195  y-axis V
1196  \endcode
1197 
1198  So, for example a vertical edge with its dark side on the left
1199  has orientation PI/2, and a horizontal edge with dark side on top
1200  has orientation PI. Obviously, the edge's orientation changes
1201  by PI if the contrast is reversed.
1202 
1203  Note that this convention changed as of VIGRA version 1.7.0.
1204 
1205  */
1207 
1208  Edgel()
1209  : x(0.0), y(0.0), strength(0.0), orientation(0.0)
1210  {}
1211 
1213  : x(ix), y(iy), strength(is), orientation(io)
1214  {}
1215 };
1216 
1217 template <class SrcIterator, class SrcAccessor,
1218  class MagnitudeImage, class BackInsertable>
1219 void internalCannyFindEdgels(SrcIterator ul, SrcAccessor grad,
1220  MagnitudeImage const & magnitude,
1221  BackInsertable & edgels)
1222 {
1223  typedef typename SrcAccessor::value_type PixelType;
1224  typedef typename PixelType::value_type ValueType;
1225 
1226  double t = 0.5 / VIGRA_CSTD::sin(M_PI/8.0);
1227 
1228  ul += Diff2D(1,1);
1229  for(int y=1; y<magnitude.height()-1; ++y, ++ul.y)
1230  {
1231  SrcIterator ix = ul;
1232  for(int x=1; x<magnitude.width()-1; ++x, ++ix.x)
1233  {
1234  double mag = magnitude(x, y);
1235  if(mag == 0.0)
1236  continue;
1237  ValueType gradx = grad.getComponent(ix, 0);
1238  ValueType grady = grad.getComponent(ix, 1);
1239 
1240  int dx = (int)VIGRA_CSTD::floor(gradx*t/mag + 0.5);
1241  int dy = (int)VIGRA_CSTD::floor(grady*t/mag + 0.5);
1242 
1243  int x1 = x - dx,
1244  x2 = x + dx,
1245  y1 = y - dy,
1246  y2 = y + dy;
1247 
1248  double m1 = magnitude(x1, y1);
1249  double m3 = magnitude(x2, y2);
1250 
1251  if(m1 < mag && m3 <= mag)
1252  {
1253  Edgel edgel;
1254 
1255  // local maximum => quadratic interpolation of sub-pixel location
1256  double del = 0.5 * (m1 - m3) / (m1 + m3 - 2.0*mag);
1257  edgel.x = Edgel::value_type(x + dx*del);
1258  edgel.y = Edgel::value_type(y + dy*del);
1259  edgel.strength = Edgel::value_type(mag);
1260  double orientation = VIGRA_CSTD::atan2(grady, gradx) + 0.5*M_PI;
1261  if(orientation < 0.0)
1262  orientation += 2.0*M_PI;
1263  edgel.orientation = Edgel::value_type(orientation);
1264  edgels.push_back(edgel);
1265  }
1266  }
1267  }
1268 }
1269 
1270 /********************************************************/
1271 /* */
1272 /* cannyEdgelList */
1273 /* */
1274 /********************************************************/
1275 
1276 /** \brief Simple implementation of Canny's edge detector.
1277 
1278  The function can be called in two modes: If you pass a 'scale', it is assumed that the
1279  original image is scalar, and the Gaussian gradient is internally computed at the
1280  given 'scale'. If the function is called without scale parameter, it is assumed that
1281  the given image already contains the gradient (i.e. its value_type must be
1282  a vector of length 2).
1283 
1284  On the basis of the gradient image, a simple non-maxima supression is performed:
1285  for each 3x3 neighborhood, it is determined whether the center pixel has
1286  larger gradient magnitude than its two neighbors in gradient direction
1287  (where the direction is rounded into octands). If this is the case,
1288  a new \ref Edgel is appended to the given vector of <TT>edgels</TT>. The subpixel
1289  edgel position is determined by fitting a parabola to the three gradient
1290  magnitude values mentioned above. The sub-pixel location of the parabola's tip
1291  and the gradient magnitude and direction (from the pixel center)
1292  are written in the newly created edgel.
1293 
1294  <b> Declarations:</b>
1295 
1296  pass arguments explicitly:
1297  \code
1298  namespace vigra {
1299  // compute edgels from a gradient image
1300  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1301  void
1302  cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1303  BackInsertable & edgels);
1304 
1305  // compute edgels from a scaler image (determine gradient internally at 'scale')
1306  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1307  void
1308  cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1309  BackInsertable & edgels, double scale);
1310  }
1311  \endcode
1312 
1313  use argument objects in conjunction with \ref ArgumentObjectFactories :
1314  \code
1315  namespace vigra {
1316  // compute edgels from a gradient image
1317  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1318  void
1319  cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1320  BackInsertable & edgels);
1321 
1322  // compute edgels from a scaler image (determine gradient internally at 'scale')
1323  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1324  void
1325  cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1326  BackInsertable & edgels, double scale);
1327  }
1328  \endcode
1329 
1330  <b> Usage:</b>
1331 
1332  <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
1333  Namespace: vigra
1334 
1335  \code
1336  vigra::BImage src(w,h);
1337 
1338  // empty edgel list
1339  std::vector<vigra::Edgel> edgels;
1340  ...
1341 
1342  // find edgels at scale 0.8
1343  vigra::cannyEdgelList(srcImageRange(src), edgels, 0.8);
1344  \endcode
1345 
1346  <b> Required Interface:</b>
1347 
1348  \code
1349  SrcImageIterator src_upperleft;
1350  SrcAccessor src_accessor;
1351 
1352  src_accessor(src_upperleft);
1353 
1354  BackInsertable edgels;
1355  edgels.push_back(Edgel());
1356  \endcode
1357 
1358  SrcAccessor::value_type must be a type convertible to float
1359 
1360  <b> Preconditions:</b>
1361 
1362  \code
1363  scale > 0
1364  \endcode
1365 */
1366 doxygen_overloaded_function(template <...> void cannyEdgelList)
1367 
1368 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1369 void
1370 cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1371  BackInsertable & edgels, double scale)
1372 {
1373  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
1374  BasicImage<TinyVector<TmpType, 2> > grad(lr-ul);
1375  gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale);
1376 
1377  cannyEdgelList(srcImageRange(grad), edgels);
1378 }
1379 
1380 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1381 inline void
1382 cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1383  BackInsertable & edgels, double scale)
1384 {
1385  cannyEdgelList(src.first, src.second, src.third, edgels, scale);
1386 }
1387 
1388 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1389 void
1390 cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1391  BackInsertable & edgels)
1392 {
1393  using namespace functor;
1394 
1395  typedef typename SrcAccessor::value_type SrcType;
1396  typedef typename NumericTraits<typename SrcType::value_type>::RealPromote TmpType;
1397  BasicImage<TmpType> magnitude(lr-ul);
1398  transformImage(srcIterRange(ul, lr, src), destImage(magnitude), norm(Arg1()));
1399 
1400  // find edgels
1401  internalCannyFindEdgels(ul, src, magnitude, edgels);
1402 }
1403 
1404 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1405 inline void
1406 cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1407  BackInsertable & edgels)
1408 {
1409  cannyEdgelList(src.first, src.second, src.third, edgels);
1410 }
1411 
1412 
1413 /********************************************************/
1414 /* */
1415 /* cannyEdgeImage */
1416 /* */
1417 /********************************************************/
1418 
1419 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
1420 
1421  This operator first calls \ref cannyEdgelList() to generate an
1422  edgel list for the given image. Then it scans this list and selects edgels
1423  whose strength is above the given <TT>gradient_threshold</TT>. For each of these
1424  edgels, the edgel's location is rounded to the nearest pixel, and that
1425  pixel marked with the given <TT>edge_marker</TT>.
1426 
1427  <b> Declarations:</b>
1428 
1429  pass arguments explicitly:
1430  \code
1431  namespace vigra {
1432  template <class SrcIterator, class SrcAccessor,
1433  class DestIterator, class DestAccessor,
1434  class GradValue, class DestValue>
1435  void cannyEdgeImage(
1436  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1437  DestIterator dul, DestAccessor da,
1438  double scale, GradValue gradient_threshold, DestValue edge_marker);
1439  }
1440  \endcode
1441 
1442  use argument objects in conjunction with \ref ArgumentObjectFactories :
1443  \code
1444  namespace vigra {
1445  template <class SrcIterator, class SrcAccessor,
1446  class DestIterator, class DestAccessor,
1447  class GradValue, class DestValue>
1448  void cannyEdgeImage(
1449  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1450  pair<DestIterator, DestAccessor> dest,
1451  double scale, GradValue gradient_threshold, DestValue edge_marker);
1452  }
1453  \endcode
1454 
1455  <b> Usage:</b>
1456 
1457  <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
1458  Namespace: vigra
1459 
1460  \code
1461  vigra::BImage src(w,h), edges(w,h);
1462 
1463  // empty edge image
1464  edges = 0;
1465  ...
1466 
1467  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
1468  vigra::cannyEdgeImage(srcImageRange(src), destImage(edges),
1469  0.8, 4.0, 1);
1470  \endcode
1471 
1472  <b> Required Interface:</b>
1473 
1474  see also: \ref cannyEdgelList().
1475 
1476  \code
1477  DestImageIterator dest_upperleft;
1478  DestAccessor dest_accessor;
1479  DestValue edge_marker;
1480 
1481  dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
1482  \endcode
1483 
1484  <b> Preconditions:</b>
1485 
1486  \code
1487  scale > 0
1488  gradient_threshold > 0
1489  \endcode
1490 */
1491 doxygen_overloaded_function(template <...> void cannyEdgeImage)
1492 
1493 template <class SrcIterator, class SrcAccessor,
1494  class DestIterator, class DestAccessor,
1495  class GradValue, class DestValue>
1496 void cannyEdgeImage(
1497  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1498  DestIterator dul, DestAccessor da,
1499  double scale, GradValue gradient_threshold, DestValue edge_marker)
1500 {
1501  std::vector<Edgel> edgels;
1502 
1503  cannyEdgelList(sul, slr, sa, edgels, scale);
1504 
1505  int w = slr.x - sul.x;
1506  int h = slr.y - sul.y;
1507 
1508  for(unsigned int i=0; i<edgels.size(); ++i)
1509  {
1510  if(gradient_threshold < edgels[i].strength)
1511  {
1512  Diff2D pix((int)(edgels[i].x + 0.5), (int)(edgels[i].y + 0.5));
1513 
1514  if(pix.x < 0 || pix.x >= w || pix.y < 0 || pix.y >= h)
1515  continue;
1516 
1517  da.set(edge_marker, dul, pix);
1518  }
1519  }
1520 }
1521 
1522 template <class SrcIterator, class SrcAccessor,
1523  class DestIterator, class DestAccessor,
1524  class GradValue, class DestValue>
1525 inline void cannyEdgeImage(
1526  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1527  pair<DestIterator, DestAccessor> dest,
1528  double scale, GradValue gradient_threshold, DestValue edge_marker)
1529 {
1530  cannyEdgeImage(src.first, src.second, src.third,
1531  dest.first, dest.second,
1532  scale, gradient_threshold, edge_marker);
1533 }
1534 
1535 /********************************************************/
1536 
1537 namespace detail {
1538 
1539 template <class DestIterator>
1540 int neighborhoodConfiguration(DestIterator dul)
1541 {
1542  int v = 0;
1543  NeighborhoodCirculator<DestIterator, EightNeighborCode> c(dul, EightNeighborCode::SouthEast);
1544  for(int i=0; i<8; ++i, --c)
1545  {
1546  v = (v << 1) | ((*c != 0) ? 1 : 0);
1547  }
1548 
1549  return v;
1550 }
1551 
1552 template <class GradValue>
1553 struct SimplePoint
1554 {
1555  Diff2D point;
1556  GradValue grad;
1557 
1558  SimplePoint(Diff2D const & p, GradValue g)
1559  : point(p), grad(g)
1560  {}
1561 
1562  bool operator<(SimplePoint const & o) const
1563  {
1564  return grad < o.grad;
1565  }
1566 
1567  bool operator>(SimplePoint const & o) const
1568  {
1569  return grad > o.grad;
1570  }
1571 };
1572 
1573 template <class SrcIterator, class SrcAccessor,
1574  class DestIterator, class DestAccessor,
1575  class GradValue, class DestValue>
1576 void cannyEdgeImageFromGrad(
1577  SrcIterator sul, SrcIterator slr, SrcAccessor grad,
1578  DestIterator dul, DestAccessor da,
1579  GradValue gradient_threshold, DestValue edge_marker)
1580 {
1581  typedef typename SrcAccessor::value_type PixelType;
1582  typedef typename NormTraits<PixelType>::SquaredNormType NormType;
1583 
1584  NormType zero = NumericTraits<NormType>::zero();
1585  double tan22_5 = M_SQRT2 - 1.0;
1586  typename NormTraits<GradValue>::SquaredNormType g2thresh = squaredNorm(gradient_threshold);
1587 
1588  int w = slr.x - sul.x;
1589  int h = slr.y - sul.y;
1590 
1591  sul += Diff2D(1,1);
1592  dul += Diff2D(1,1);
1593  Diff2D p(0,0);
1594 
1595  for(int y = 1; y < h-1; ++y, ++sul.y, ++dul.y)
1596  {
1597  SrcIterator sx = sul;
1598  DestIterator dx = dul;
1599  for(int x = 1; x < w-1; ++x, ++sx.x, ++dx.x)
1600  {
1601  PixelType g = grad(sx);
1602  NormType g2n = squaredNorm(g);
1603  if(g2n < g2thresh)
1604  continue;
1605 
1606  NormType g2n1, g2n3;
1607  // find out quadrant
1608  if(abs(g[1]) < tan22_5*abs(g[0]))
1609  {
1610  // north-south edge
1611  g2n1 = squaredNorm(grad(sx, Diff2D(-1, 0)));
1612  g2n3 = squaredNorm(grad(sx, Diff2D(1, 0)));
1613  }
1614  else if(abs(g[0]) < tan22_5*abs(g[1]))
1615  {
1616  // west-east edge
1617  g2n1 = squaredNorm(grad(sx, Diff2D(0, -1)));
1618  g2n3 = squaredNorm(grad(sx, Diff2D(0, 1)));
1619  }
1620  else if(g[0]*g[1] < zero)
1621  {
1622  // north-west-south-east edge
1623  g2n1 = squaredNorm(grad(sx, Diff2D(1, -1)));
1624  g2n3 = squaredNorm(grad(sx, Diff2D(-1, 1)));
1625  }
1626  else
1627  {
1628  // north-east-south-west edge
1629  g2n1 = squaredNorm(grad(sx, Diff2D(-1, -1)));
1630  g2n3 = squaredNorm(grad(sx, Diff2D(1, 1)));
1631  }
1632 
1633  if(g2n1 < g2n && g2n3 <= g2n)
1634  {
1635  da.set(edge_marker, dx);
1636  }
1637  }
1638  }
1639 }
1640 
1641 } // namespace detail
1642 
1643 /********************************************************/
1644 /* */
1645 /* cannyEdgeImageWithThinning */
1646 /* */
1647 /********************************************************/
1648 
1649 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
1650 
1651  The input pixels of this algorithms must be vectors of length 2 (see Required Interface below).
1652  It first searches for all pixels whose gradient magnitude is larger
1653  than the given <tt>gradient_threshold</tt> and larger than the magnitude of its two neighbors
1654  in gradient direction (where these neighbors are determined by nearest neighbor
1655  interpolation, i.e. according to the octant where the gradient points into).
1656  The resulting edge pixel candidates are then subjected to topological thinning
1657  so that the remaining edge pixels can be linked into edgel chains with a provable,
1658  non-heuristic algorithm. Thinning is performed so that the pixels with highest gradient
1659  magnitude survive. Optionally, the outermost pixels are marked as edge pixels
1660  as well when <tt>addBorder</tt> is true. The remaining pixels will be marked in the destination
1661  image with the value of <tt>edge_marker</tt> (all non-edge pixels remain untouched).
1662 
1663  <b> Declarations:</b>
1664 
1665  pass arguments explicitly:
1666  \code
1667  namespace vigra {
1668  template <class SrcIterator, class SrcAccessor,
1669  class DestIterator, class DestAccessor,
1670  class GradValue, class DestValue>
1671  void cannyEdgeImageFromGradWithThinning(
1672  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1673  DestIterator dul, DestAccessor da,
1674  GradValue gradient_threshold,
1675  DestValue edge_marker, bool addBorder = true);
1676  }
1677  \endcode
1678 
1679  use argument objects in conjunction with \ref ArgumentObjectFactories :
1680  \code
1681  namespace vigra {
1682  template <class SrcIterator, class SrcAccessor,
1683  class DestIterator, class DestAccessor,
1684  class GradValue, class DestValue>
1685  void cannyEdgeImageFromGradWithThinning(
1686  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1687  pair<DestIterator, DestAccessor> dest,
1688  GradValue gradient_threshold,
1689  DestValue edge_marker, bool addBorder = true);
1690  }
1691  \endcode
1692 
1693  <b> Usage:</b>
1694 
1695  <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
1696  Namespace: vigra
1697 
1698  \code
1699  vigra::BImage src(w,h), edges(w,h);
1700 
1701  vigra::FVector2Image grad(w,h);
1702  // compute the image gradient at scale 0.8
1703  vigra::gaussianGradient(srcImageRange(src), destImage(grad), 0.8);
1704 
1705  // empty edge image
1706  edges = 0;
1707  // find edges gradient larger than 4.0, mark with 1, and add border
1708  vigra::cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destImage(edges),
1709  4.0, 1, true);
1710  \endcode
1711 
1712  <b> Required Interface:</b>
1713 
1714  \code
1715  // the input pixel type must be a vector with two elements
1716  SrcImageIterator src_upperleft;
1717  SrcAccessor src_accessor;
1718  typedef SrcAccessor::value_type SrcPixel;
1719  typedef NormTraits<SrcPixel>::SquaredNormType SrcSquaredNormType;
1720 
1721  SrcPixel g = src_accessor(src_upperleft);
1722  SrcPixel::value_type g0 = g[0];
1723  SrcSquaredNormType gn = squaredNorm(g);
1724 
1725  DestImageIterator dest_upperleft;
1726  DestAccessor dest_accessor;
1727  DestValue edge_marker;
1728 
1729  dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
1730  \endcode
1731 
1732  <b> Preconditions:</b>
1733 
1734  \code
1735  gradient_threshold > 0
1736  \endcode
1737 */
1739 
1740 template <class SrcIterator, class SrcAccessor,
1741  class DestIterator, class DestAccessor,
1742  class GradValue, class DestValue>
1744  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1745  DestIterator dul, DestAccessor da,
1746  GradValue gradient_threshold,
1747  DestValue edge_marker, bool addBorder)
1748 {
1749  int w = slr.x - sul.x;
1750  int h = slr.y - sul.y;
1751 
1752  BImage edgeImage(w, h, BImage::value_type(0));
1753  BImage::traverser eul = edgeImage.upperLeft();
1754  BImage::Accessor ea = edgeImage.accessor();
1755  if(addBorder)
1756  initImageBorder(destImageRange(edgeImage), 1, 1);
1757  detail::cannyEdgeImageFromGrad(sul, slr, sa, eul, ea, gradient_threshold, 1);
1758 
1759  static bool isSimplePoint[256] = {
1760  0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0,
1761  0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0,
1762  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1,
1763  1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1,
1764  0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1,
1765  0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0,
1766  0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0,
1767  1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1,
1768  0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
1769  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1770  0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,
1771  0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0,
1772  1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0,
1773  0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1,
1774  1, 0, 1, 0 };
1775 
1776  eul += Diff2D(1,1);
1777  sul += Diff2D(1,1);
1778  int w2 = w-2;
1779  int h2 = h-2;
1780 
1781  typedef detail::SimplePoint<GradValue> SP;
1782  // use std::greater becaus we need the smallest gradients at the top of the queue
1783  std::priority_queue<SP, std::vector<SP>, std::greater<SP> > pqueue;
1784 
1785  Diff2D p(0,0);
1786  for(; p.y < h2; ++p.y)
1787  {
1788  for(p.x = 0; p.x < w2; ++p.x)
1789  {
1790  BImage::traverser e = eul + p;
1791  if(*e == 0)
1792  continue;
1793  int v = detail::neighborhoodConfiguration(e);
1794  if(isSimplePoint[v])
1795  {
1796  pqueue.push(SP(p, norm(sa(sul+p))));
1797  *e = 2; // remember that it is already in queue
1798  }
1799  }
1800  }
1801 
1802  static const Diff2D dist[] = { Diff2D(-1,0), Diff2D(0,-1),
1803  Diff2D(1,0), Diff2D(0,1) };
1804 
1805  while(pqueue.size())
1806  {
1807  p = pqueue.top().point;
1808  pqueue.pop();
1809 
1810  BImage::traverser e = eul + p;
1811  int v = detail::neighborhoodConfiguration(e);
1812  if(!isSimplePoint[v])
1813  continue; // point may no longer be simple because its neighbors changed
1814 
1815  *e = 0; // delete simple point
1816 
1817  for(int i=0; i<4; ++i)
1818  {
1819  Diff2D pneu = p + dist[i];
1820  if(pneu.x == -1 || pneu.y == -1 || pneu.x == w2 || pneu.y == h2)
1821  continue; // do not remove points at the border
1822 
1823  BImage::traverser eneu = eul + pneu;
1824  if(*eneu == 1) // point is boundary and not yet in the queue
1825  {
1826  int v = detail::neighborhoodConfiguration(eneu);
1827  if(isSimplePoint[v])
1828  {
1829  pqueue.push(SP(pneu, norm(sa(sul+pneu))));
1830  *eneu = 2; // remember that it is already in queue
1831  }
1832  }
1833  }
1834  }
1835 
1836  initImageIf(destIterRange(dul, dul+Diff2D(w,h), da),
1837  maskImage(edgeImage), edge_marker);
1838 }
1839 
1840 template <class SrcIterator, class SrcAccessor,
1841  class DestIterator, class DestAccessor,
1842  class GradValue, class DestValue>
1844  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1845  pair<DestIterator, DestAccessor> dest,
1846  GradValue gradient_threshold,
1847  DestValue edge_marker, bool addBorder)
1848 {
1849  cannyEdgeImageFromGradWithThinning(src.first, src.second, src.third,
1850  dest.first, dest.second,
1851  gradient_threshold, edge_marker, addBorder);
1852 }
1853 
1854 template <class SrcIterator, class SrcAccessor,
1855  class DestIterator, class DestAccessor,
1856  class GradValue, class DestValue>
1858  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1859  DestIterator dul, DestAccessor da,
1860  GradValue gradient_threshold, DestValue edge_marker)
1861 {
1863  dul, da,
1864  gradient_threshold, edge_marker, true);
1865 }
1866 
1867 template <class SrcIterator, class SrcAccessor,
1868  class DestIterator, class DestAccessor,
1869  class GradValue, class DestValue>
1871  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1872  pair<DestIterator, DestAccessor> dest,
1873  GradValue gradient_threshold, DestValue edge_marker)
1874 {
1875  cannyEdgeImageFromGradWithThinning(src.first, src.second, src.third,
1876  dest.first, dest.second,
1877  gradient_threshold, edge_marker, true);
1878 }
1879 
1880 /********************************************************/
1881 /* */
1882 /* cannyEdgeImageWithThinning */
1883 /* */
1884 /********************************************************/
1885 
1886 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
1887 
1888  This operator first calls \ref gaussianGradient() to compute the gradient of the input
1889  image, ad then \ref cannyEdgeImageFromGradWithThinning() to generate an
1890  edge image. See there for more detailed documentation.
1891 
1892  <b> Declarations:</b>
1893 
1894  pass arguments explicitly:
1895  \code
1896  namespace vigra {
1897  template <class SrcIterator, class SrcAccessor,
1898  class DestIterator, class DestAccessor,
1899  class GradValue, class DestValue>
1900  void cannyEdgeImageWithThinning(
1901  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1902  DestIterator dul, DestAccessor da,
1903  double scale, GradValue gradient_threshold,
1904  DestValue edge_marker, bool addBorder = true);
1905  }
1906  \endcode
1907 
1908  use argument objects in conjunction with \ref ArgumentObjectFactories :
1909  \code
1910  namespace vigra {
1911  template <class SrcIterator, class SrcAccessor,
1912  class DestIterator, class DestAccessor,
1913  class GradValue, class DestValue>
1914  void cannyEdgeImageWithThinning(
1915  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1916  pair<DestIterator, DestAccessor> dest,
1917  double scale, GradValue gradient_threshold,
1918  DestValue edge_marker, bool addBorder = true);
1919  }
1920  \endcode
1921 
1922  <b> Usage:</b>
1923 
1924  <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
1925  Namespace: vigra
1926 
1927  \code
1928  vigra::BImage src(w,h), edges(w,h);
1929 
1930  // empty edge image
1931  edges = 0;
1932  ...
1933 
1934  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1, annd add border
1935  vigra::cannyEdgeImageWithThinning(srcImageRange(src), destImage(edges),
1936  0.8, 4.0, 1, true);
1937  \endcode
1938 
1939  <b> Required Interface:</b>
1940 
1941  see also: \ref cannyEdgelList().
1942 
1943  \code
1944  DestImageIterator dest_upperleft;
1945  DestAccessor dest_accessor;
1946  DestValue edge_marker;
1947 
1948  dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
1949  \endcode
1950 
1951  <b> Preconditions:</b>
1952 
1953  \code
1954  scale > 0
1955  gradient_threshold > 0
1956  \endcode
1957 */
1959 
1960 template <class SrcIterator, class SrcAccessor,
1961  class DestIterator, class DestAccessor,
1962  class GradValue, class DestValue>
1964  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1965  DestIterator dul, DestAccessor da,
1966  double scale, GradValue gradient_threshold,
1967  DestValue edge_marker, bool addBorder)
1968 {
1969  // mark pixels that are higher than their neighbors in gradient direction
1970  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
1971  BasicImage<TinyVector<TmpType, 2> > grad(slr-sul);
1972  gaussianGradient(srcIterRange(sul, slr, sa), destImage(grad), scale);
1973  cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destIter(dul, da),
1974  gradient_threshold, edge_marker, addBorder);
1975 }
1976 
1977 template <class SrcIterator, class SrcAccessor,
1978  class DestIterator, class DestAccessor,
1979  class GradValue, class DestValue>
1980 inline void cannyEdgeImageWithThinning(
1981  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1982  pair<DestIterator, DestAccessor> dest,
1983  double scale, GradValue gradient_threshold,
1984  DestValue edge_marker, bool addBorder)
1985 {
1986  cannyEdgeImageWithThinning(src.first, src.second, src.third,
1987  dest.first, dest.second,
1988  scale, gradient_threshold, edge_marker, addBorder);
1989 }
1990 
1991 template <class SrcIterator, class SrcAccessor,
1992  class DestIterator, class DestAccessor,
1993  class GradValue, class DestValue>
1994 inline void cannyEdgeImageWithThinning(
1995  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1996  DestIterator dul, DestAccessor da,
1997  double scale, GradValue gradient_threshold, DestValue edge_marker)
1998 {
1999  cannyEdgeImageWithThinning(sul, slr, sa,
2000  dul, da,
2001  scale, gradient_threshold, edge_marker, true);
2002 }
2003 
2004 template <class SrcIterator, class SrcAccessor,
2005  class DestIterator, class DestAccessor,
2006  class GradValue, class DestValue>
2007 inline void cannyEdgeImageWithThinning(
2008  triple<SrcIterator, SrcIterator, SrcAccessor> src,
2009  pair<DestIterator, DestAccessor> dest,
2010  double scale, GradValue gradient_threshold, DestValue edge_marker)
2011 {
2012  cannyEdgeImageWithThinning(src.first, src.second, src.third,
2013  dest.first, dest.second,
2014  scale, gradient_threshold, edge_marker, true);
2015 }
2016 
2017 /********************************************************/
2018 
2019 template <class SrcIterator, class SrcAccessor,
2020  class MaskImage, class BackInsertable>
2021 void internalCannyFindEdgels3x3(SrcIterator ul, SrcAccessor grad,
2022  MaskImage const & mask,
2023  BackInsertable & edgels)
2024 {
2025  typedef typename SrcAccessor::value_type PixelType;
2026  typedef typename PixelType::value_type ValueType;
2027 
2028  ul += Diff2D(1,1);
2029  for(int y=1; y<mask.height()-1; ++y, ++ul.y)
2030  {
2031  SrcIterator ix = ul;
2032  for(int x=1; x<mask.width()-1; ++x, ++ix.x)
2033  {
2034  if(!mask(x,y))
2035  continue;
2036 
2037  ValueType gradx = grad.getComponent(ix, 0);
2038  ValueType grady = grad.getComponent(ix, 1);
2039  double mag = hypot(gradx, grady);
2040  if(mag == 0.0)
2041  continue;
2042  double c = gradx / mag,
2043  s = grady / mag;
2044 
2045  Matrix<double> ml(3,3), mr(3,1), l(3,1), r(3,1);
2046  l(0,0) = 1.0;
2047 
2048  for(int yy = -1; yy <= 1; ++yy)
2049  {
2050  for(int xx = -1; xx <= 1; ++xx)
2051  {
2052  double u = c*xx + s*yy;
2053  double v = norm(grad(ix, Diff2D(xx, yy)));
2054  l(1,0) = u;
2055  l(2,0) = u*u;
2056  ml += outer(l);
2057  mr += v*l;
2058  }
2059  }
2060 
2061  linearSolve(ml, mr, r);
2062 
2063  Edgel edgel;
2064 
2065  // local maximum => quadratic interpolation of sub-pixel location
2066  double del = -r(1,0) / 2.0 / r(2,0);
2067  if(std::fabs(del) > 1.5) // don't move by more than about a pixel diameter
2068  del = 0.0;
2069  edgel.x = Edgel::value_type(x + c*del);
2070  edgel.y = Edgel::value_type(y + s*del);
2071  edgel.strength = Edgel::value_type(mag);
2072  double orientation = VIGRA_CSTD::atan2(grady, gradx) + 0.5*M_PI;
2073  if(orientation < 0.0)
2074  orientation += 2.0*M_PI;
2075  edgel.orientation = Edgel::value_type(orientation);
2076  edgels.push_back(edgel);
2077  }
2078  }
2079 }
2080 
2081 
2082 /********************************************************/
2083 /* */
2084 /* cannyEdgelList3x3 */
2085 /* */
2086 /********************************************************/
2087 
2088 /** \brief Improved implementation of Canny's edge detector.
2089 
2090  This operator first computes pixels which are crossed by the edge using
2091  cannyEdgeImageWithThinning(). The gradient magnitudes in the 3x3 neighborhood of these
2092  pixels are then projected onto the normal of the edge (as determined
2093  by the gradient direction). The edgel's subpixel location is found by fitting a
2094  parabola through the 9 gradient values and determining the parabola's tip.
2095  A new \ref Edgel is appended to the given vector of <TT>edgels</TT>. Since the parabola
2096  is fitted to 9 points rather than 3 points as in cannyEdgelList(), the accuracy is higher.
2097 
2098  The function can be called in two modes: If you pass a 'scale', it is assumed that the
2099  original image is scalar, and the Gaussian gradient is internally computed at the
2100  given 'scale'. If the function is called without scale parameter, it is assumed that
2101  the given image already contains the gradient (i.e. its value_type must be
2102  a vector of length 2).
2103 
2104  <b> Declarations:</b>
2105 
2106  pass arguments explicitly:
2107  \code
2108  namespace vigra {
2109  // compute edgels from a gradient image
2110  template <class SrcIterator, class SrcAccessor, class BackInsertable>
2111  void cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2112  BackInsertable & edgels);
2113 
2114  // compute edgels from a scaler image (determine gradient internally at 'scale')
2115  template <class SrcIterator, class SrcAccessor, class BackInsertable>
2116  void cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2117  BackInsertable & edgels, double scale);
2118  }
2119  \endcode
2120 
2121  use argument objects in conjunction with \ref ArgumentObjectFactories :
2122  \code
2123  namespace vigra {
2124  // compute edgels from a gradient image
2125  template <class SrcIterator, class SrcAccessor, class BackInsertable>
2126  void
2127  cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2128  BackInsertable & edgels);
2129 
2130  // compute edgels from a scaler image (determine gradient internally at 'scale')
2131  template <class SrcIterator, class SrcAccessor, class BackInsertable>
2132  void
2133  cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2134  BackInsertable & edgels, double scale);
2135  }
2136  \endcode
2137 
2138  <b> Usage:</b>
2139 
2140  <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
2141  Namespace: vigra
2142 
2143  \code
2144  vigra::BImage src(w,h);
2145 
2146  // empty edgel list
2147  std::vector<vigra::Edgel> edgels;
2148  ...
2149 
2150  // find edgels at scale 0.8
2151  vigra::cannyEdgelList3x3(srcImageRange(src), edgels, 0.8);
2152  \endcode
2153 
2154  <b> Required Interface:</b>
2155 
2156  \code
2157  SrcImageIterator src_upperleft;
2158  SrcAccessor src_accessor;
2159 
2160  src_accessor(src_upperleft);
2161 
2162  BackInsertable edgels;
2163  edgels.push_back(Edgel());
2164  \endcode
2165 
2166  SrcAccessor::value_type must be a type convertible to float
2167 
2168  <b> Preconditions:</b>
2169 
2170  \code
2171  scale > 0
2172  \endcode
2173 */
2175 
2176 template <class SrcIterator, class SrcAccessor, class BackInsertable>
2177 void
2178 cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2179  BackInsertable & edgels, double scale)
2180 {
2181  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
2182  BasicImage<TinyVector<TmpType, 2> > grad(lr-ul);
2183  gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale);
2184 
2185  cannyEdgelList3x3(srcImageRange(grad), edgels);
2186 }
2187 
2188 template <class SrcIterator, class SrcAccessor, class BackInsertable>
2189 inline void
2190 cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2191  BackInsertable & edgels, double scale)
2192 {
2193  cannyEdgelList3x3(src.first, src.second, src.third, edgels, scale);
2194 }
2195 
2196 template <class SrcIterator, class SrcAccessor, class BackInsertable>
2197 void
2198 cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2199  BackInsertable & edgels)
2200 {
2201  UInt8Image edges(lr-ul);
2202  cannyEdgeImageFromGradWithThinning(srcIterRange(ul, lr, src), destImage(edges),
2203  0.0, 1, false);
2204 
2205  // find edgels
2206  internalCannyFindEdgels3x3(ul, src, edges, edgels);
2207 }
2208 
2209 template <class SrcIterator, class SrcAccessor, class BackInsertable>
2210 inline void
2211 cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2212  BackInsertable & edgels)
2213 {
2214  cannyEdgelList3x3(src.first, src.second, src.third, edgels);
2215 }
2216 
2217 //@}
2218 
2219 /** \page CrackEdgeImage Crack Edge Image
2220 
2221 Crack edges are marked <i>between</i> the pixels of an image.
2222 A Crack Edge Image is an image that represents these edges. In order
2223 to accomodate the cracks, the Crack Edge Image must be twice as large
2224 as the original image (precisely (2*w - 1) by (2*h - 1)). A Crack Edge Image
2225 can easily be derived from a binary image or from the signs of the
2226 response of a Laplacean filter. Consider the following sketch, where
2227 <TT>+</TT> encodes the foreground, <TT>-</TT> the background, and
2228 <TT>*</TT> the resulting crack edges.
2229 
2230  \code
2231 sign of difference image insert cracks resulting CrackEdgeImage
2232 
2233  + . - . - . * . . .
2234  + - - . . . . . . * * * .
2235  + + - => + . + . - => . . . * .
2236  + + + . . . . . . . . * *
2237  + . + . + . . . . .
2238  \endcode
2239 
2240 Starting from the original binary image (left), we insert crack pixels
2241 to get to the double-sized image (center). Finally, we mark all
2242 crack pixels whose non-crack neighbors have different signs as
2243 crack edge points, while all other pixels (crack and non-crack) become
2244 region pixels.
2245 
2246 <b>Requirements on a Crack Edge Image:</b>
2247 
2248 <ul>
2249  <li>Crack Edge Images have odd width and height.
2250  <li>Crack pixels have at least one odd coordinate.
2251  <li>Only crack pixels may be marked as edge points.
2252  <li>Crack pixels with two odd coordinates must be marked as edge points
2253  whenever any of their neighboring crack pixels was marked.
2254 </ul>
2255 
2256 The last two requirements ensure that both edges and regions are 4-connected.
2257 Thus, 4-connectivity and 8-connectivity yield identical connected
2258 components in a Crack Edge Image (so called <i>well-composedness</i>).
2259 This ensures that Crack Edge Images have nice topological properties
2260 (cf. L. J. Latecki: "Well-Composed Sets", Academic Press, 2000).
2261 */
2262 
2263 
2264 } // namespace vigra
2265 
2266 #endif // VIGRA_EDGEDETECTION_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.7.1 (Wed Mar 12 2014)