[ 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> <vigra/edgedetection.hxx><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 */
180 doxygen_overloaded_function(template <...> void differenceOfExponentialEdgeImage)
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 conforms 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> <vigra/edgedetection.hxx><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 */
465 doxygen_overloaded_function(template <...> void differenceOfExponentialCrackEdgeImage)
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> <vigra/edgedetection.hxx><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 */
749 doxygen_overloaded_function(template <...> void removeShortEdges)
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> <vigra/edgedetection.hxx><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 */
885 doxygen_overloaded_function(template <...> void closeGapsInCrackEdgeImage)
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  Therefore, this algorithm should only be applied as a visualization 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> <vigra/edgedetection.hxx><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 */
1102 doxygen_overloaded_function(template <...> void beautifyCrackEdgeImage)
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, class GradValue>
1219 void internalCannyFindEdgels(SrcIterator ul, SrcAccessor grad,
1220  MagnitudeImage const & magnitude,
1221  BackInsertable & edgels, GradValue grad_thresh)
1222 {
1223  typedef typename SrcAccessor::value_type PixelType;
1224  typedef typename PixelType::value_type ValueType;
1225 
1226  vigra_precondition(grad_thresh >= NumericTraits<GradValue>::zero(),
1227  "cannyFindEdgels(): gradient threshold must not be negative.");
1228 
1229  double t = 0.5 / VIGRA_CSTD::sin(M_PI/8.0);
1230 
1231  ul += Diff2D(1,1);
1232  for(int y=1; y<magnitude.height()-1; ++y, ++ul.y)
1233  {
1234  SrcIterator ix = ul;
1235  for(int x=1; x<magnitude.width()-1; ++x, ++ix.x)
1236  {
1237  double mag = magnitude(x, y);
1238  if(mag <= grad_thresh)
1239  continue;
1240  ValueType gradx = grad.getComponent(ix, 0);
1241  ValueType grady = grad.getComponent(ix, 1);
1242 
1243  int dx = (int)VIGRA_CSTD::floor(gradx*t/mag + 0.5);
1244  int dy = (int)VIGRA_CSTD::floor(grady*t/mag + 0.5);
1245 
1246  int x1 = x - dx,
1247  x2 = x + dx,
1248  y1 = y - dy,
1249  y2 = y + dy;
1250 
1251  double m1 = magnitude(x1, y1);
1252  double m3 = magnitude(x2, y2);
1253 
1254  if(m1 < mag && m3 <= mag)
1255  {
1256  Edgel edgel;
1257 
1258  // local maximum => quadratic interpolation of sub-pixel location
1259  double del = 0.5 * (m1 - m3) / (m1 + m3 - 2.0*mag);
1260  edgel.x = Edgel::value_type(x + dx*del);
1261  edgel.y = Edgel::value_type(y + dy*del);
1262  edgel.strength = Edgel::value_type(mag);
1263  double orientation = VIGRA_CSTD::atan2(grady, gradx) + 0.5*M_PI;
1264  if(orientation < 0.0)
1265  orientation += 2.0*M_PI;
1266  edgel.orientation = Edgel::value_type(orientation);
1267  edgels.push_back(edgel);
1268  }
1269  }
1270  }
1271 }
1272 
1273 /********************************************************/
1274 /* */
1275 /* cannyEdgelList */
1276 /* */
1277 /********************************************************/
1278 
1279 /** \brief Simple implementation of Canny's edge detector.
1280 
1281  The function can be called in two modes: If you pass a 'scale', it is assumed that the
1282  original image is scalar, and the Gaussian gradient is internally computed at the
1283  given 'scale'. If the function is called without scale parameter, it is assumed that
1284  the given image already contains the gradient (i.e. its value_type must be
1285  a vector of length 2).
1286 
1287  On the basis of the gradient image, a simple non-maxima suppression is performed:
1288  for each 3x3 neighborhood, it is determined whether the center pixel has
1289  larger gradient magnitude than its two neighbors in gradient direction
1290  (where the direction is rounded into octants). If this is the case,
1291  a new \ref Edgel is appended to the given vector of <TT>edgels</TT>. The subpixel
1292  edgel position is determined by fitting a parabola to the three gradient
1293  magnitude values mentioned above. The sub-pixel location of the parabola's tip
1294  and the gradient magnitude and direction (from the pixel center)
1295  are written in the newly created edgel.
1296 
1297  <b> Declarations:</b>
1298 
1299  pass arguments explicitly:
1300  \code
1301  namespace vigra {
1302  // compute edgels from a gradient image
1303  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1304  void
1305  cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1306  BackInsertable & edgels);
1307 
1308  // compute edgels from a scalar image (determine gradient internally at 'scale')
1309  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1310  void
1311  cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1312  BackInsertable & edgels, double scale);
1313  }
1314  \endcode
1315 
1316  use argument objects in conjunction with \ref ArgumentObjectFactories :
1317  \code
1318  namespace vigra {
1319  // compute edgels from a gradient image
1320  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1321  void
1322  cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1323  BackInsertable & edgels);
1324 
1325  // compute edgels from a scalar image (determine gradient internally at 'scale')
1326  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1327  void
1328  cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1329  BackInsertable & edgels, double scale);
1330  }
1331  \endcode
1332 
1333  <b> Usage:</b>
1334 
1335  <b>\#include</b> <vigra/edgedetection.hxx><br>
1336  Namespace: vigra
1337 
1338  \code
1339  vigra::BImage src(w,h);
1340 
1341  // empty edgel list
1342  std::vector<vigra::Edgel> edgels;
1343  ...
1344 
1345  // find edgels at scale 0.8
1346  vigra::cannyEdgelList(srcImageRange(src), edgels, 0.8);
1347  \endcode
1348 
1349  <b> Required Interface:</b>
1350 
1351  \code
1352  SrcImageIterator src_upperleft;
1353  SrcAccessor src_accessor;
1354 
1355  src_accessor(src_upperleft);
1356 
1357  BackInsertable edgels;
1358  edgels.push_back(Edgel());
1359  \endcode
1360 
1361  SrcAccessor::value_type must be a type convertible to float
1362 
1363  <b> Preconditions:</b>
1364 
1365  \code
1366  scale > 0
1367  \endcode
1368 */
1369 doxygen_overloaded_function(template <...> void cannyEdgelList)
1370 
1371 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1372 void
1373 cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1374  BackInsertable & edgels, double scale)
1375 {
1376  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
1377  BasicImage<TinyVector<TmpType, 2> > grad(lr-ul);
1378  gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale);
1379 
1380  cannyEdgelList(srcImageRange(grad), edgels);
1381 }
1382 
1383 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1384 inline void
1385 cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1386  BackInsertable & edgels, double scale)
1387 {
1388  cannyEdgelList(src.first, src.second, src.third, edgels, scale);
1389 }
1390 
1391 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1392 void
1393 cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1394  BackInsertable & edgels)
1395 {
1396  using namespace functor;
1397 
1398  typedef typename SrcAccessor::value_type SrcType;
1399  typedef typename NumericTraits<typename SrcType::value_type>::RealPromote TmpType;
1400  BasicImage<TmpType> magnitude(lr-ul);
1401  transformImage(srcIterRange(ul, lr, src), destImage(magnitude), norm(Arg1()));
1402 
1403  // find edgels
1404  internalCannyFindEdgels(ul, src, magnitude, edgels, NumericTraits<TmpType>::zero());
1405 }
1406 
1407 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1408 inline void
1409 cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1410  BackInsertable & edgels)
1411 {
1412  cannyEdgelList(src.first, src.second, src.third, edgels);
1413 }
1414 
1415 /********************************************************/
1416 /* */
1417 /* cannyEdgelListThreshold */
1418 /* */
1419 /********************************************************/
1420 
1421 /** \brief Canny's edge detector with thresholding.
1422 
1423  This function works exactly like \ref cannyEdgelList(), but
1424  you also pass a threshold for the minimal gradient magnitude,
1425  so that edgels whose strength is below the threshold are not
1426  inserted into the edgel list.
1427 
1428  <b> Declarations:</b>
1429 
1430  pass arguments explicitly:
1431  \code
1432  namespace vigra {
1433  // compute edgels from a gradient image
1434  template <class SrcIterator, class SrcAccessor,
1435  class BackInsertable, class GradValue>
1436  void
1437  cannyEdgelListThreshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1438  BackInsertable & edgels, GradValue grad_threshold);
1439 
1440  // compute edgels from a scalar image (determine gradient internally at 'scale')
1441  template <class SrcIterator, class SrcAccessor,
1442  class BackInsertable, class GradValue>
1443  void
1444  cannyEdgelListThreshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1445  BackInsertable & edgels, double scale, GradValue grad_threshold);
1446  }
1447  \endcode
1448 
1449  use argument objects in conjunction with \ref ArgumentObjectFactories :
1450  \code
1451  namespace vigra {
1452  // compute edgels from a gradient image
1453  template <class SrcIterator, class SrcAccessor,
1454  class BackInsertable, class GradValue>
1455  void
1456  cannyEdgelListThreshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1457  BackInsertable & edgels, GradValue grad_threshold);
1458 
1459  // compute edgels from a scalar image (determine gradient internally at 'scale')
1460  template <class SrcIterator, class SrcAccessor,
1461  class BackInsertable, class GradValue>
1462  void
1463  cannyEdgelListThreshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1464  BackInsertable & edgels, double scale, GradValue grad_threshold);
1465  }
1466  \endcode
1467 
1468  <b> Usage:</b>
1469 
1470  <b>\#include</b> <vigra/edgedetection.hxx><br>
1471  Namespace: vigra
1472 
1473  \code
1474  vigra::BImage src(w,h);
1475 
1476  // empty edgel list
1477  std::vector<vigra::Edgel> edgels;
1478  ...
1479 
1480  // find edgels at scale 0.8, only considering gradient above 2.0
1481  vigra::cannyEdgelListThreshold(srcImageRange(src), edgels, 0.8, 2.0);
1482  \endcode
1483 
1484  <b> Required Interface:</b>
1485 
1486  \code
1487  SrcImageIterator src_upperleft;
1488  SrcAccessor src_accessor;
1489 
1490  src_accessor(src_upperleft);
1491 
1492  BackInsertable edgels;
1493  edgels.push_back(Edgel());
1494  \endcode
1495 
1496  SrcAccessor::value_type must be a type convertible to float
1497 
1498  <b> Preconditions:</b>
1499 
1500  \code
1501  scale > 0
1502  \endcode
1503 */
1504 doxygen_overloaded_function(template <...> void cannyEdgelListThreshold)
1505 
1506 template <class SrcIterator, class SrcAccessor,
1507  class BackInsertable, class GradValue>
1508 void
1509 cannyEdgelListThreshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1510  BackInsertable & edgels, double scale, GradValue grad_threshold)
1511 {
1512  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
1513  BasicImage<TinyVector<TmpType, 2> > grad(lr-ul);
1514  gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale);
1515 
1516  cannyEdgelListThreshold(srcImageRange(grad), edgels, grad_threshold);
1517 }
1518 
1519 template <class SrcIterator, class SrcAccessor,
1520  class BackInsertable, class GradValue>
1521 inline void
1522 cannyEdgelListThreshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1523  BackInsertable & edgels, double scale, GradValue grad_threshold)
1524 {
1525  cannyEdgelListThreshold(src.first, src.second, src.third, edgels, scale, grad_threshold);
1526 }
1527 
1528 template <class SrcIterator, class SrcAccessor,
1529  class BackInsertable, class GradValue>
1530 void
1531 cannyEdgelListThreshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1532  BackInsertable & edgels, GradValue grad_threshold)
1533 {
1534  using namespace functor;
1535 
1536  typedef typename SrcAccessor::value_type SrcType;
1537  typedef typename NumericTraits<typename SrcType::value_type>::RealPromote TmpType;
1538  BasicImage<TmpType> magnitude(lr-ul);
1539  transformImage(srcIterRange(ul, lr, src), destImage(magnitude), norm(Arg1()));
1540 
1541  // find edgels
1542  internalCannyFindEdgels(ul, src, magnitude, edgels, grad_threshold);
1543 }
1544 
1545 template <class SrcIterator, class SrcAccessor,
1546  class BackInsertable, class GradValue>
1547 inline void
1548 cannyEdgelListThreshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1549  BackInsertable & edgels, GradValue grad_threshold)
1550 {
1551  cannyEdgelListThreshold(src.first, src.second, src.third, edgels, grad_threshold);
1552 }
1553 
1554 
1555 /********************************************************/
1556 /* */
1557 /* cannyEdgeImage */
1558 /* */
1559 /********************************************************/
1560 
1561 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
1562 
1563  This operator first calls \ref cannyEdgelList() to generate an
1564  edgel list for the given image. Then it scans this list and selects edgels
1565  whose strength is above the given <TT>gradient_threshold</TT>. For each of these
1566  edgels, the edgel's location is rounded to the nearest pixel, and that
1567  pixel marked with the given <TT>edge_marker</TT>.
1568 
1569  <b> Declarations:</b>
1570 
1571  pass arguments explicitly:
1572  \code
1573  namespace vigra {
1574  template <class SrcIterator, class SrcAccessor,
1575  class DestIterator, class DestAccessor,
1576  class GradValue, class DestValue>
1577  void cannyEdgeImage(
1578  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1579  DestIterator dul, DestAccessor da,
1580  double scale, GradValue gradient_threshold, DestValue edge_marker);
1581  }
1582  \endcode
1583 
1584  use argument objects in conjunction with \ref ArgumentObjectFactories :
1585  \code
1586  namespace vigra {
1587  template <class SrcIterator, class SrcAccessor,
1588  class DestIterator, class DestAccessor,
1589  class GradValue, class DestValue>
1590  void cannyEdgeImage(
1591  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1592  pair<DestIterator, DestAccessor> dest,
1593  double scale, GradValue gradient_threshold, DestValue edge_marker);
1594  }
1595  \endcode
1596 
1597  <b> Usage:</b>
1598 
1599  <b>\#include</b> <vigra/edgedetection.hxx><br>
1600  Namespace: vigra
1601 
1602  \code
1603  vigra::BImage src(w,h), edges(w,h);
1604 
1605  // empty edge image
1606  edges = 0;
1607  ...
1608 
1609  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
1610  vigra::cannyEdgeImage(srcImageRange(src), destImage(edges),
1611  0.8, 4.0, 1);
1612  \endcode
1613 
1614  <b> Required Interface:</b>
1615 
1616  see also: \ref cannyEdgelList().
1617 
1618  \code
1619  DestImageIterator dest_upperleft;
1620  DestAccessor dest_accessor;
1621  DestValue edge_marker;
1622 
1623  dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
1624  \endcode
1625 
1626  <b> Preconditions:</b>
1627 
1628  \code
1629  scale > 0
1630  gradient_threshold > 0
1631  \endcode
1632 */
1633 doxygen_overloaded_function(template <...> void cannyEdgeImage)
1634 
1635 template <class SrcIterator, class SrcAccessor,
1636  class DestIterator, class DestAccessor,
1637  class GradValue, class DestValue>
1638 void cannyEdgeImage(
1639  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1640  DestIterator dul, DestAccessor da,
1641  double scale, GradValue gradient_threshold, DestValue edge_marker)
1642 {
1643  std::vector<Edgel> edgels;
1644 
1645  cannyEdgelListThreshold(sul, slr, sa, edgels, scale, gradient_threshold);
1646 
1647  int w = slr.x - sul.x;
1648  int h = slr.y - sul.y;
1649 
1650  for(unsigned int i=0; i<edgels.size(); ++i)
1651  {
1652  Diff2D pix((int)(edgels[i].x + 0.5), (int)(edgels[i].y + 0.5));
1653 
1654  if(pix.x < 0 || pix.x >= w || pix.y < 0 || pix.y >= h)
1655  continue;
1656 
1657  da.set(edge_marker, dul, pix);
1658  }
1659 }
1660 
1661 template <class SrcIterator, class SrcAccessor,
1662  class DestIterator, class DestAccessor,
1663  class GradValue, class DestValue>
1664 inline void cannyEdgeImage(
1665  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1666  pair<DestIterator, DestAccessor> dest,
1667  double scale, GradValue gradient_threshold, DestValue edge_marker)
1668 {
1669  cannyEdgeImage(src.first, src.second, src.third,
1670  dest.first, dest.second,
1671  scale, gradient_threshold, edge_marker);
1672 }
1673 
1674 /********************************************************/
1675 
1676 namespace detail {
1677 
1678 template <class DestIterator>
1679 int neighborhoodConfiguration(DestIterator dul)
1680 {
1681  int v = 0;
1682  NeighborhoodCirculator<DestIterator, EightNeighborCode> c(dul, EightNeighborCode::SouthEast);
1683  for(int i=0; i<8; ++i, --c)
1684  {
1685  v = (v << 1) | ((*c != 0) ? 1 : 0);
1686  }
1687 
1688  return v;
1689 }
1690 
1691 template <class GradValue>
1692 struct SimplePoint
1693 {
1694  Diff2D point;
1695  GradValue grad;
1696 
1697  SimplePoint(Diff2D const & p, GradValue g)
1698  : point(p), grad(g)
1699  {}
1700 
1701  bool operator<(SimplePoint const & o) const
1702  {
1703  return grad < o.grad;
1704  }
1705 
1706  bool operator>(SimplePoint const & o) const
1707  {
1708  return grad > o.grad;
1709  }
1710 };
1711 
1712 template <class SrcIterator, class SrcAccessor,
1713  class DestIterator, class DestAccessor,
1714  class GradValue, class DestValue>
1715 void cannyEdgeImageFromGrad(
1716  SrcIterator sul, SrcIterator slr, SrcAccessor grad,
1717  DestIterator dul, DestAccessor da,
1718  GradValue gradient_threshold, DestValue edge_marker)
1719 {
1720  typedef typename SrcAccessor::value_type PixelType;
1721  typedef typename NormTraits<PixelType>::SquaredNormType NormType;
1722 
1723  NormType zero = NumericTraits<NormType>::zero();
1724  double tan22_5 = M_SQRT2 - 1.0;
1725  typename NormTraits<GradValue>::SquaredNormType g2thresh = squaredNorm(gradient_threshold);
1726 
1727  int w = slr.x - sul.x;
1728  int h = slr.y - sul.y;
1729 
1730  sul += Diff2D(1,1);
1731  dul += Diff2D(1,1);
1732  Diff2D p(0,0);
1733 
1734  for(int y = 1; y < h-1; ++y, ++sul.y, ++dul.y)
1735  {
1736  SrcIterator sx = sul;
1737  DestIterator dx = dul;
1738  for(int x = 1; x < w-1; ++x, ++sx.x, ++dx.x)
1739  {
1740  PixelType g = grad(sx);
1741  NormType g2n = squaredNorm(g);
1742  if(g2n < g2thresh)
1743  continue;
1744 
1745  NormType g2n1, g2n3;
1746  // find out quadrant
1747  if(abs(g[1]) < tan22_5*abs(g[0]))
1748  {
1749  // north-south edge
1750  g2n1 = squaredNorm(grad(sx, Diff2D(-1, 0)));
1751  g2n3 = squaredNorm(grad(sx, Diff2D(1, 0)));
1752  }
1753  else if(abs(g[0]) < tan22_5*abs(g[1]))
1754  {
1755  // west-east edge
1756  g2n1 = squaredNorm(grad(sx, Diff2D(0, -1)));
1757  g2n3 = squaredNorm(grad(sx, Diff2D(0, 1)));
1758  }
1759  else if(g[0]*g[1] < zero)
1760  {
1761  // north-west-south-east edge
1762  g2n1 = squaredNorm(grad(sx, Diff2D(1, -1)));
1763  g2n3 = squaredNorm(grad(sx, Diff2D(-1, 1)));
1764  }
1765  else
1766  {
1767  // north-east-south-west edge
1768  g2n1 = squaredNorm(grad(sx, Diff2D(-1, -1)));
1769  g2n3 = squaredNorm(grad(sx, Diff2D(1, 1)));
1770  }
1771 
1772  if(g2n1 < g2n && g2n3 <= g2n)
1773  {
1774  da.set(edge_marker, dx);
1775  }
1776  }
1777  }
1778 }
1779 
1780 } // namespace detail
1781 
1782 /********************************************************/
1783 /* */
1784 /* cannyEdgeImageWithThinning */
1785 /* */
1786 /********************************************************/
1787 
1788 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
1789 
1790  The input pixels of this algorithms must be vectors of length 2 (see Required Interface below).
1791  It first searches for all pixels whose gradient magnitude is larger
1792  than the given <tt>gradient_threshold</tt> and larger than the magnitude of its two neighbors
1793  in gradient direction (where these neighbors are determined by nearest neighbor
1794  interpolation, i.e. according to the octant where the gradient points into).
1795  The resulting edge pixel candidates are then subjected to topological thinning
1796  so that the remaining edge pixels can be linked into edgel chains with a provable,
1797  non-heuristic algorithm. Thinning is performed so that the pixels with highest gradient
1798  magnitude survive. Optionally, the outermost pixels are marked as edge pixels
1799  as well when <tt>addBorder</tt> is true. The remaining pixels will be marked in the destination
1800  image with the value of <tt>edge_marker</tt> (all non-edge pixels remain untouched).
1801 
1802  <b> Declarations:</b>
1803 
1804  pass arguments explicitly:
1805  \code
1806  namespace vigra {
1807  template <class SrcIterator, class SrcAccessor,
1808  class DestIterator, class DestAccessor,
1809  class GradValue, class DestValue>
1810  void cannyEdgeImageFromGradWithThinning(
1811  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1812  DestIterator dul, DestAccessor da,
1813  GradValue gradient_threshold,
1814  DestValue edge_marker, bool addBorder = true);
1815  }
1816  \endcode
1817 
1818  use argument objects in conjunction with \ref ArgumentObjectFactories :
1819  \code
1820  namespace vigra {
1821  template <class SrcIterator, class SrcAccessor,
1822  class DestIterator, class DestAccessor,
1823  class GradValue, class DestValue>
1824  void cannyEdgeImageFromGradWithThinning(
1825  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1826  pair<DestIterator, DestAccessor> dest,
1827  GradValue gradient_threshold,
1828  DestValue edge_marker, bool addBorder = true);
1829  }
1830  \endcode
1831 
1832  <b> Usage:</b>
1833 
1834  <b>\#include</b> <vigra/edgedetection.hxx><br>
1835  Namespace: vigra
1836 
1837  \code
1838  vigra::BImage src(w,h), edges(w,h);
1839 
1840  vigra::FVector2Image grad(w,h);
1841  // compute the image gradient at scale 0.8
1842  vigra::gaussianGradient(srcImageRange(src), destImage(grad), 0.8);
1843 
1844  // empty edge image
1845  edges = 0;
1846  // find edges gradient larger than 4.0, mark with 1, and add border
1847  vigra::cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destImage(edges),
1848  4.0, 1, true);
1849  \endcode
1850 
1851  <b> Required Interface:</b>
1852 
1853  \code
1854  // the input pixel type must be a vector with two elements
1855  SrcImageIterator src_upperleft;
1856  SrcAccessor src_accessor;
1857  typedef SrcAccessor::value_type SrcPixel;
1858  typedef NormTraits<SrcPixel>::SquaredNormType SrcSquaredNormType;
1859 
1860  SrcPixel g = src_accessor(src_upperleft);
1861  SrcPixel::value_type g0 = g[0];
1862  SrcSquaredNormType gn = squaredNorm(g);
1863 
1864  DestImageIterator dest_upperleft;
1865  DestAccessor dest_accessor;
1866  DestValue edge_marker;
1867 
1868  dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
1869  \endcode
1870 
1871  <b> Preconditions:</b>
1872 
1873  \code
1874  gradient_threshold > 0
1875  \endcode
1876 */
1877 doxygen_overloaded_function(template <...> void cannyEdgeImageFromGradWithThinning)
1878 
1879 template <class SrcIterator, class SrcAccessor,
1880  class DestIterator, class DestAccessor,
1881  class GradValue, class DestValue>
1883  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1884  DestIterator dul, DestAccessor da,
1885  GradValue gradient_threshold,
1886  DestValue edge_marker, bool addBorder)
1887 {
1888  int w = slr.x - sul.x;
1889  int h = slr.y - sul.y;
1890 
1891  BImage edgeImage(w, h, BImage::value_type(0));
1892  BImage::traverser eul = edgeImage.upperLeft();
1893  BImage::Accessor ea = edgeImage.accessor();
1894  if(addBorder)
1895  initImageBorder(destImageRange(edgeImage), 1, 1);
1896  detail::cannyEdgeImageFromGrad(sul, slr, sa, eul, ea, gradient_threshold, 1);
1897 
1898  static bool isSimplePoint[256] = {
1899  0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0,
1900  0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0,
1901  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1,
1902  1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1,
1903  0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1,
1904  0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0,
1905  0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0,
1906  1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1,
1907  0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
1908  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1909  0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,
1910  0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0,
1911  1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0,
1912  0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1,
1913  1, 0, 1, 0 };
1914 
1915  eul += Diff2D(1,1);
1916  sul += Diff2D(1,1);
1917  int w2 = w-2;
1918  int h2 = h-2;
1919 
1920  typedef detail::SimplePoint<GradValue> SP;
1921  // use std::greater because we need the smallest gradients at the top of the queue
1922  std::priority_queue<SP, std::vector<SP>, std::greater<SP> > pqueue;
1923 
1924  Diff2D p(0,0);
1925  for(; p.y < h2; ++p.y)
1926  {
1927  for(p.x = 0; p.x < w2; ++p.x)
1928  {
1929  BImage::traverser e = eul + p;
1930  if(*e == 0)
1931  continue;
1932  int v = detail::neighborhoodConfiguration(e);
1933  if(isSimplePoint[v])
1934  {
1935  pqueue.push(SP(p, norm(sa(sul+p))));
1936  *e = 2; // remember that it is already in queue
1937  }
1938  }
1939  }
1940 
1941  static const Diff2D dist[] = { Diff2D(-1,0), Diff2D(0,-1),
1942  Diff2D(1,0), Diff2D(0,1) };
1943 
1944  while(pqueue.size())
1945  {
1946  p = pqueue.top().point;
1947  pqueue.pop();
1948 
1949  BImage::traverser e = eul + p;
1950  int v = detail::neighborhoodConfiguration(e);
1951  if(!isSimplePoint[v])
1952  continue; // point may no longer be simple because its neighbors changed
1953 
1954  *e = 0; // delete simple point
1955 
1956  for(int i=0; i<4; ++i)
1957  {
1958  Diff2D pneu = p + dist[i];
1959  if(pneu.x == -1 || pneu.y == -1 || pneu.x == w2 || pneu.y == h2)
1960  continue; // do not remove points at the border
1961 
1962  BImage::traverser eneu = eul + pneu;
1963  if(*eneu == 1) // point is boundary and not yet in the queue
1964  {
1965  int v = detail::neighborhoodConfiguration(eneu);
1966  if(isSimplePoint[v])
1967  {
1968  pqueue.push(SP(pneu, norm(sa(sul+pneu))));
1969  *eneu = 2; // remember that it is already in queue
1970  }
1971  }
1972  }
1973  }
1974 
1975  initImageIf(destIterRange(dul, dul+Diff2D(w,h), da),
1976  maskImage(edgeImage), edge_marker);
1977 }
1978 
1979 template <class SrcIterator, class SrcAccessor,
1980  class DestIterator, class DestAccessor,
1981  class GradValue, class DestValue>
1983  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1984  pair<DestIterator, DestAccessor> dest,
1985  GradValue gradient_threshold,
1986  DestValue edge_marker, bool addBorder)
1987 {
1988  cannyEdgeImageFromGradWithThinning(src.first, src.second, src.third,
1989  dest.first, dest.second,
1990  gradient_threshold, edge_marker, addBorder);
1991 }
1992 
1993 template <class SrcIterator, class SrcAccessor,
1994  class DestIterator, class DestAccessor,
1995  class GradValue, class DestValue>
1997  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1998  DestIterator dul, DestAccessor da,
1999  GradValue gradient_threshold, DestValue edge_marker)
2000 {
2002  dul, da,
2003  gradient_threshold, edge_marker, true);
2004 }
2005 
2006 template <class SrcIterator, class SrcAccessor,
2007  class DestIterator, class DestAccessor,
2008  class GradValue, class DestValue>
2010  triple<SrcIterator, SrcIterator, SrcAccessor> src,
2011  pair<DestIterator, DestAccessor> dest,
2012  GradValue gradient_threshold, DestValue edge_marker)
2013 {
2014  cannyEdgeImageFromGradWithThinning(src.first, src.second, src.third,
2015  dest.first, dest.second,
2016  gradient_threshold, edge_marker, true);
2017 }
2018 
2019 /********************************************************/
2020 /* */
2021 /* cannyEdgeImageWithThinning */
2022 /* */
2023 /********************************************************/
2024 
2025 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
2026 
2027  This operator first calls \ref gaussianGradient() to compute the gradient of the input
2028  image, ad then \ref cannyEdgeImageFromGradWithThinning() to generate an
2029  edge image. See there for more detailed documentation.
2030 
2031  <b> Declarations:</b>
2032 
2033  pass arguments explicitly:
2034  \code
2035  namespace vigra {
2036  template <class SrcIterator, class SrcAccessor,
2037  class DestIterator, class DestAccessor,
2038  class GradValue, class DestValue>
2039  void cannyEdgeImageWithThinning(
2040  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
2041  DestIterator dul, DestAccessor da,
2042  double scale, GradValue gradient_threshold,
2043  DestValue edge_marker, bool addBorder = true);
2044  }
2045  \endcode
2046 
2047  use argument objects in conjunction with \ref ArgumentObjectFactories :
2048  \code
2049  namespace vigra {
2050  template <class SrcIterator, class SrcAccessor,
2051  class DestIterator, class DestAccessor,
2052  class GradValue, class DestValue>
2053  void cannyEdgeImageWithThinning(
2054  triple<SrcIterator, SrcIterator, SrcAccessor> src,
2055  pair<DestIterator, DestAccessor> dest,
2056  double scale, GradValue gradient_threshold,
2057  DestValue edge_marker, bool addBorder = true);
2058  }
2059  \endcode
2060 
2061  <b> Usage:</b>
2062 
2063  <b>\#include</b> <vigra/edgedetection.hxx><br>
2064  Namespace: vigra
2065 
2066  \code
2067  vigra::BImage src(w,h), edges(w,h);
2068 
2069  // empty edge image
2070  edges = 0;
2071  ...
2072 
2073  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1, annd add border
2074  vigra::cannyEdgeImageWithThinning(srcImageRange(src), destImage(edges),
2075  0.8, 4.0, 1, true);
2076  \endcode
2077 
2078  <b> Required Interface:</b>
2079 
2080  see also: \ref cannyEdgelList().
2081 
2082  \code
2083  DestImageIterator dest_upperleft;
2084  DestAccessor dest_accessor;
2085  DestValue edge_marker;
2086 
2087  dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
2088  \endcode
2089 
2090  <b> Preconditions:</b>
2091 
2092  \code
2093  scale > 0
2094  gradient_threshold > 0
2095  \endcode
2096 */
2097 doxygen_overloaded_function(template <...> void cannyEdgeImageWithThinning)
2098 
2099 template <class SrcIterator, class SrcAccessor,
2100  class DestIterator, class DestAccessor,
2101  class GradValue, class DestValue>
2103  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
2104  DestIterator dul, DestAccessor da,
2105  double scale, GradValue gradient_threshold,
2106  DestValue edge_marker, bool addBorder)
2107 {
2108  // mark pixels that are higher than their neighbors in gradient direction
2109  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
2110  BasicImage<TinyVector<TmpType, 2> > grad(slr-sul);
2111  gaussianGradient(srcIterRange(sul, slr, sa), destImage(grad), scale);
2112  cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destIter(dul, da),
2113  gradient_threshold, edge_marker, addBorder);
2114 }
2115 
2116 template <class SrcIterator, class SrcAccessor,
2117  class DestIterator, class DestAccessor,
2118  class GradValue, class DestValue>
2119 inline void cannyEdgeImageWithThinning(
2120  triple<SrcIterator, SrcIterator, SrcAccessor> src,
2121  pair<DestIterator, DestAccessor> dest,
2122  double scale, GradValue gradient_threshold,
2123  DestValue edge_marker, bool addBorder)
2124 {
2125  cannyEdgeImageWithThinning(src.first, src.second, src.third,
2126  dest.first, dest.second,
2127  scale, gradient_threshold, edge_marker, addBorder);
2128 }
2129 
2130 template <class SrcIterator, class SrcAccessor,
2131  class DestIterator, class DestAccessor,
2132  class GradValue, class DestValue>
2133 inline void cannyEdgeImageWithThinning(
2134  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
2135  DestIterator dul, DestAccessor da,
2136  double scale, GradValue gradient_threshold, DestValue edge_marker)
2137 {
2138  cannyEdgeImageWithThinning(sul, slr, sa,
2139  dul, da,
2140  scale, gradient_threshold, edge_marker, true);
2141 }
2142 
2143 template <class SrcIterator, class SrcAccessor,
2144  class DestIterator, class DestAccessor,
2145  class GradValue, class DestValue>
2146 inline void cannyEdgeImageWithThinning(
2147  triple<SrcIterator, SrcIterator, SrcAccessor> src,
2148  pair<DestIterator, DestAccessor> dest,
2149  double scale, GradValue gradient_threshold, DestValue edge_marker)
2150 {
2151  cannyEdgeImageWithThinning(src.first, src.second, src.third,
2152  dest.first, dest.second,
2153  scale, gradient_threshold, edge_marker, true);
2154 }
2155 
2156 /********************************************************/
2157 
2158 template <class SrcIterator, class SrcAccessor,
2159  class MaskImage, class BackInsertable, class GradValue>
2160 void internalCannyFindEdgels3x3(SrcIterator ul, SrcAccessor grad,
2161  MaskImage const & mask,
2162  BackInsertable & edgels,
2163  GradValue grad_thresh)
2164 {
2165  typedef typename SrcAccessor::value_type PixelType;
2166  typedef typename PixelType::value_type ValueType;
2167 
2168  vigra_precondition(grad_thresh >= NumericTraits<GradValue>::zero(),
2169  "cannyFindEdgels3x3(): gradient threshold must not be negative.");
2170 
2171  ul += Diff2D(1,1);
2172  for(int y=1; y<mask.height()-1; ++y, ++ul.y)
2173  {
2174  SrcIterator ix = ul;
2175  for(int x=1; x<mask.width()-1; ++x, ++ix.x)
2176  {
2177  if(!mask(x,y))
2178  continue;
2179 
2180  ValueType gradx = grad.getComponent(ix, 0);
2181  ValueType grady = grad.getComponent(ix, 1);
2182  double mag = hypot(gradx, grady);
2183  if(mag <= grad_thresh)
2184  continue;
2185  double c = gradx / mag,
2186  s = grady / mag;
2187 
2188  Matrix<double> ml(3,3), mr(3,1), l(3,1), r(3,1);
2189  l(0,0) = 1.0;
2190 
2191  for(int yy = -1; yy <= 1; ++yy)
2192  {
2193  for(int xx = -1; xx <= 1; ++xx)
2194  {
2195  double u = c*xx + s*yy;
2196  double v = norm(grad(ix, Diff2D(xx, yy)));
2197  l(1,0) = u;
2198  l(2,0) = u*u;
2199  ml += outer(l);
2200  mr += v*l;
2201  }
2202  }
2203 
2204  linearSolve(ml, mr, r);
2205 
2206  Edgel edgel;
2207 
2208  // local maximum => quadratic interpolation of sub-pixel location
2209  double del = -r(1,0) / 2.0 / r(2,0);
2210  if(std::fabs(del) > 1.5) // don't move by more than about a pixel diameter
2211  del = 0.0;
2212  edgel.x = Edgel::value_type(x + c*del);
2213  edgel.y = Edgel::value_type(y + s*del);
2214  edgel.strength = Edgel::value_type(mag);
2215  double orientation = VIGRA_CSTD::atan2(grady, gradx) + 0.5*M_PI;
2216  if(orientation < 0.0)
2217  orientation += 2.0*M_PI;
2218  edgel.orientation = Edgel::value_type(orientation);
2219  edgels.push_back(edgel);
2220  }
2221  }
2222 }
2223 
2224 
2225 /********************************************************/
2226 /* */
2227 /* cannyEdgelList3x3 */
2228 /* */
2229 /********************************************************/
2230 
2231 /** \brief Improved implementation of Canny's edge detector.
2232 
2233  This operator first computes pixels which are crossed by the edge using
2234  cannyEdgeImageWithThinning(). The gradient magnitudes in the 3x3 neighborhood of these
2235  pixels are then projected onto the normal of the edge (as determined
2236  by the gradient direction). The edgel's subpixel location is found by fitting a
2237  parabola through the 9 gradient values and determining the parabola's tip.
2238  A new \ref Edgel is appended to the given vector of <TT>edgels</TT>. Since the parabola
2239  is fitted to 9 points rather than 3 points as in cannyEdgelList(), the accuracy is higher.
2240 
2241  The function can be called in two modes: If you pass a 'scale', it is assumed that the
2242  original image is scalar, and the Gaussian gradient is internally computed at the
2243  given 'scale'. If the function is called without scale parameter, it is assumed that
2244  the given image already contains the gradient (i.e. its value_type must be
2245  a vector of length 2).
2246 
2247  <b> Declarations:</b>
2248 
2249  pass arguments explicitly:
2250  \code
2251  namespace vigra {
2252  // compute edgels from a gradient image
2253  template <class SrcIterator, class SrcAccessor, class BackInsertable>
2254  void cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2255  BackInsertable & edgels);
2256 
2257  // compute edgels from a scalar image (determine gradient internally at 'scale')
2258  template <class SrcIterator, class SrcAccessor, class BackInsertable>
2259  void cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2260  BackInsertable & edgels, double scale);
2261  }
2262  \endcode
2263 
2264  use argument objects in conjunction with \ref ArgumentObjectFactories :
2265  \code
2266  namespace vigra {
2267  // compute edgels from a gradient image
2268  template <class SrcIterator, class SrcAccessor, class BackInsertable>
2269  void
2270  cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2271  BackInsertable & edgels);
2272 
2273  // compute edgels from a scalar image (determine gradient internally at 'scale')
2274  template <class SrcIterator, class SrcAccessor, class BackInsertable>
2275  void
2276  cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2277  BackInsertable & edgels, double scale);
2278  }
2279  \endcode
2280 
2281  <b> Usage:</b>
2282 
2283  <b>\#include</b> <vigra/edgedetection.hxx><br>
2284  Namespace: vigra
2285 
2286  \code
2287  vigra::BImage src(w,h);
2288 
2289  // empty edgel list
2290  std::vector<vigra::Edgel> edgels;
2291  ...
2292 
2293  // find edgels at scale 0.8
2294  vigra::cannyEdgelList3x3(srcImageRange(src), edgels, 0.8);
2295  \endcode
2296 
2297  <b> Required Interface:</b>
2298 
2299  \code
2300  SrcImageIterator src_upperleft;
2301  SrcAccessor src_accessor;
2302 
2303  src_accessor(src_upperleft);
2304 
2305  BackInsertable edgels;
2306  edgels.push_back(Edgel());
2307  \endcode
2308 
2309  SrcAccessor::value_type must be a type convertible to float
2310 
2311  <b> Preconditions:</b>
2312 
2313  \code
2314  scale > 0
2315  \endcode
2316 */
2317 doxygen_overloaded_function(template <...> void cannyEdgelList3x3)
2318 
2319 template <class SrcIterator, class SrcAccessor, class BackInsertable>
2320 void
2321 cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2322  BackInsertable & edgels, double scale)
2323 {
2324  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
2325  BasicImage<TinyVector<TmpType, 2> > grad(lr-ul);
2326  gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale);
2327 
2328  cannyEdgelList3x3(srcImageRange(grad), edgels);
2329 }
2330 
2331 template <class SrcIterator, class SrcAccessor, class BackInsertable>
2332 inline void
2333 cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2334  BackInsertable & edgels, double scale)
2335 {
2336  cannyEdgelList3x3(src.first, src.second, src.third, edgels, scale);
2337 }
2338 
2339 template <class SrcIterator, class SrcAccessor, class BackInsertable>
2340 void
2341 cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2342  BackInsertable & edgels)
2343 {
2344  typedef typename NormTraits<typename SrcAccessor::value_type>::NormType NormType;
2345 
2346  UInt8Image edges(lr-ul);
2347  cannyEdgeImageFromGradWithThinning(srcIterRange(ul, lr, src), destImage(edges),
2348  0.0, 1, false);
2349 
2350  // find edgels
2351  internalCannyFindEdgels3x3(ul, src, edges, edgels, NumericTraits<NormType>::zero());
2352 }
2353 
2354 template <class SrcIterator, class SrcAccessor, class BackInsertable>
2355 inline void
2356 cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2357  BackInsertable & edgels)
2358 {
2359  cannyEdgelList3x3(src.first, src.second, src.third, edgels);
2360 }
2361 
2362 /********************************************************/
2363 /* */
2364 /* cannyEdgelList3x3Threshold */
2365 /* */
2366 /********************************************************/
2367 
2368 /** \brief Improved implementation of Canny's edge detector with thresholding.
2369 
2370  This function works exactly like \ref cannyEdgelList3x3(), but
2371  you also pass a threshold for the minimal gradient magnitude,
2372  so that edgels whose strength is below the threshold are not
2373  inserted into the edgel list.
2374 
2375  <b> Declarations:</b>
2376 
2377  pass arguments explicitly:
2378  \code
2379  namespace vigra {
2380  // compute edgels from a gradient image
2381  template <class SrcIterator, class SrcAccessor,
2382  class BackInsertable, class GradValue>
2383  void
2384  cannyEdgelList3x3Threshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2385  BackInsertable & edgels, GradValue grad_thresh);
2386 
2387  // compute edgels from a scalar image (determine gradient internally at 'scale')
2388  template <class SrcIterator, class SrcAccessor,
2389  class BackInsertable, class GradValue>
2390  void
2391  cannyEdgelList3x3Threshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2392  BackInsertable & edgels, double scale, GradValue grad_thresh);
2393  }
2394  \endcode
2395 
2396  use argument objects in conjunction with \ref ArgumentObjectFactories :
2397  \code
2398  namespace vigra {
2399  // compute edgels from a gradient image
2400  template <class SrcIterator, class SrcAccessor,
2401  class BackInsertable, class GradValue>
2402  void
2403  cannyEdgelList3x3Threshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2404  BackInsertable & edgels, GradValue grad_thresh);
2405 
2406  // compute edgels from a scalar image (determine gradient internally at 'scale')
2407  template <class SrcIterator, class SrcAccessor,
2408  class BackInsertable, class GradValue>
2409  void
2410  cannyEdgelList3x3Threshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2411  BackInsertable & edgels, double scale, GradValue grad_thresh);
2412  }
2413  \endcode
2414 
2415  <b> Usage:</b>
2416 
2417  <b>\#include</b> <vigra/edgedetection.hxx><br>
2418  Namespace: vigra
2419 
2420  \code
2421  vigra::BImage src(w,h);
2422 
2423  // empty edgel list
2424  std::vector<vigra::Edgel> edgels;
2425  ...
2426 
2427  // find edgels at scale 0.8
2428  vigra::cannyEdgelList3x3(srcImageRange(src), edgels, 0.8);
2429  \endcode
2430 
2431  <b> Required Interface:</b>
2432 
2433  \code
2434  SrcImageIterator src_upperleft;
2435  SrcAccessor src_accessor;
2436 
2437  src_accessor(src_upperleft);
2438 
2439  BackInsertable edgels;
2440  edgels.push_back(Edgel());
2441  \endcode
2442 
2443  SrcAccessor::value_type must be a type convertible to float
2444 
2445  <b> Preconditions:</b>
2446 
2447  \code
2448  scale > 0
2449  \endcode
2450 */
2451 doxygen_overloaded_function(template <...> void cannyEdgelList3x3Threshold)
2452 
2453 template <class SrcIterator, class SrcAccessor,
2454  class BackInsertable, class GradValue>
2455 void
2456 cannyEdgelList3x3Threshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2457  BackInsertable & edgels, double scale, GradValue grad_thresh)
2458 {
2459  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
2460  BasicImage<TinyVector<TmpType, 2> > grad(lr-ul);
2461  gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale);
2462 
2463  cannyEdgelList3x3Threshold(srcImageRange(grad), edgels, grad_thresh);
2464 }
2465 
2466 template <class SrcIterator, class SrcAccessor,
2467  class BackInsertable, class GradValue>
2468 inline void
2469 cannyEdgelList3x3Threshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2470  BackInsertable & edgels, double scale, GradValue grad_thresh)
2471 {
2472  cannyEdgelList3x3Threshold(src.first, src.second, src.third, edgels, scale, grad_thresh);
2473 }
2474 
2475 template <class SrcIterator, class SrcAccessor,
2476  class BackInsertable, class GradValue>
2477 void
2478 cannyEdgelList3x3Threshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2479  BackInsertable & edgels, GradValue grad_thresh)
2480 {
2481  UInt8Image edges(lr-ul);
2482  cannyEdgeImageFromGradWithThinning(srcIterRange(ul, lr, src), destImage(edges),
2483  0.0, 1, false);
2484 
2485  // find edgels
2486  internalCannyFindEdgels3x3(ul, src, edges, edgels, grad_thresh);
2487 }
2488 
2489 template <class SrcIterator, class SrcAccessor,
2490  class BackInsertable, class GradValue>
2491 inline void
2492 cannyEdgelList3x3Threshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2493  BackInsertable & edgels, GradValue grad_thresh)
2494 {
2495  cannyEdgelList3x3Threshold(src.first, src.second, src.third, edgels, grad_thresh);
2496 }
2497 
2498 //@}
2499 
2500 /** \page CrackEdgeImage Crack Edge Image
2501 
2502 Crack edges are marked <i>between</i> the pixels of an image.
2503 A Crack Edge Image is an image that represents these edges. In order
2504 to accommodate the cracks, the Crack Edge Image must be twice as large
2505 as the original image (precisely (2*w - 1) by (2*h - 1)). A Crack Edge Image
2506 can easily be derived from a binary image or from the signs of the
2507 response of a Laplacian filter. Consider the following sketch, where
2508 <TT>+</TT> encodes the foreground, <TT>-</TT> the background, and
2509 <TT>*</TT> the resulting crack edges.
2510 
2511  \code
2512 sign of difference image insert cracks resulting CrackEdgeImage
2513 
2514  + . - . - . * . . .
2515  + - - . . . . . . * * * .
2516  + + - => + . + . - => . . . * .
2517  + + + . . . . . . . . * *
2518  + . + . + . . . . .
2519  \endcode
2520 
2521 Starting from the original binary image (left), we insert crack pixels
2522 to get to the double-sized image (center). Finally, we mark all
2523 crack pixels whose non-crack neighbors have different signs as
2524 crack edge points, while all other pixels (crack and non-crack) become
2525 region pixels.
2526 
2527 <b>Requirements on a Crack Edge Image:</b>
2528 
2529 <ul>
2530  <li>Crack Edge Images have odd width and height.
2531  <li>Crack pixels have at least one odd coordinate.
2532  <li>Only crack pixels may be marked as edge points.
2533  <li>Crack pixels with two odd coordinates must be marked as edge points
2534  whenever any of their neighboring crack pixels was marked.
2535 </ul>
2536 
2537 The last two requirements ensure that both edges and regions are 4-connected.
2538 Thus, 4-connectivity and 8-connectivity yield identical connected
2539 components in a Crack Edge Image (so called <i>well-composedness</i>).
2540 This ensures that Crack Edge Images have nice topological properties
2541 (cf. L. J. Latecki: "Well-Composed Sets", Academic Press, 2000).
2542 */
2543 
2544 
2545 } // namespace vigra
2546 
2547 #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.9.0 (Tue Sep 24 2013)