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

edgedetection.hxx VIGRA

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 #include "multi_shape.hxx"
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 2D array views:
104  \code
105  namespace vigra {
106  template <class T1, class S1,
107  class T2, class S2,
108  class GradValue, class DestValue>
109  void
110  differenceOfExponentialEdgeImage(MultiArrayView<2, T1, S1> const & src,
111  MultiArrayView<2, T2, S2> dest,
112  double scale,
113  GradValue gradient_threshold,
114  DestValue edge_marker = NumericTraits<DestValue>::one());
115  }
116  \endcode
117 
118  \deprecatedAPI{differenceOfExponentialEdgeImage}
119  pass \ref ImageIterators and \ref DataAccessors :
120  \code
121  namespace vigra {
122  template <class SrcIterator, class SrcAccessor,
123  class DestIterator, class DestAccessor,
124  class GradValue,
125  class DestValue = DestAccessor::value_type>
126  void differenceOfExponentialEdgeImage(
127  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
128  DestIterator dul, DestAccessor da,
129  double scale, GradValue gradient_threshold,
130  DestValue edge_marker = NumericTraits<DestValue>::one())
131  }
132  \endcode
133  use argument objects in conjunction with \ref ArgumentObjectFactories :
134  \code
135  namespace vigra {
136  template <class SrcIterator, class SrcAccessor,
137  class DestIterator, class DestAccessor,
138  class GradValue,
139  class DestValue = DestAccessor::value_type>
140  void differenceOfExponentialEdgeImage(
141  triple<SrcIterator, SrcIterator, SrcAccessor> src,
142  pair<DestIterator, DestAccessor> dest,
143  double scale, GradValue gradient_threshold,
144  DestValue edge_marker = NumericTraits<DestValue>::one())
145  }
146  \endcode
147  \deprecatedEnd
148 
149  <b> Usage:</b>
150 
151  <b>\#include</b> <vigra/edgedetection.hxx><br>
152  Namespace: vigra
153 
154  \code
155  MultiArray<2, unsigned char> src(w,h), edges(w,h);
156  ...
157 
158  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
159  differenceOfExponentialEdgeImage(src, edges,
160  0.8, 4.0, 1);
161  \endcode
162 
163  \deprecatedUsage{differenceOfExponentialEdgeImage}
164  \code
165  vigra::BImage src(w,h), edges(w,h);
166 
167  // empty edge image
168  edges = 0;
169  ...
170 
171  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
172  vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges),
173  0.8, 4.0, 1);
174  \endcode
175  <b> Required Interface:</b>
176  \code
177  SrcImageIterator src_upperleft, src_lowerright;
178  DestImageIterator dest_upperleft;
179 
180  SrcAccessor src_accessor;
181  DestAccessor dest_accessor;
182 
183  SrcAccessor::value_type u = src_accessor(src_upperleft);
184  double d;
185  GradValue gradient_threshold;
186 
187  u = u + u
188  u = u - u
189  u = u * u
190  u = d * u
191  u < gradient_threshold
192 
193  DestValue edge_marker;
194  dest_accessor.set(edge_marker, dest_upperleft);
195  \endcode
196  \deprecatedEnd
197 
198  <b> Preconditions:</b>
199 
200  \code
201  scale > 0
202  gradient_threshold > 0
203  \endcode
204 */
205 doxygen_overloaded_function(template <...> void differenceOfExponentialEdgeImage)
206 
207 template <class SrcIterator, class SrcAccessor,
208  class DestIterator, class DestAccessor,
209  class GradValue, class DestValue>
211  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
212  DestIterator dul, DestAccessor da,
213  double scale, GradValue gradient_threshold, DestValue edge_marker)
214 {
215  vigra_precondition(scale > 0,
216  "differenceOfExponentialEdgeImage(): scale > 0 required.");
217 
218  vigra_precondition(gradient_threshold > 0,
219  "differenceOfExponentialEdgeImage(): "
220  "gradient_threshold > 0 required.");
221 
222  int w = slr.x - sul.x;
223  int h = slr.y - sul.y;
224  int x,y;
225 
226  typedef typename
227  NumericTraits<typename SrcAccessor::value_type>::RealPromote
228  TMPTYPE;
229  typedef BasicImage<TMPTYPE> TMPIMG;
230 
231  TMPIMG tmp(w,h);
232  TMPIMG smooth(w,h);
233 
234  recursiveSmoothX(srcIterRange(sul, slr, sa), destImage(tmp), scale / 2.0);
235  recursiveSmoothY(srcImageRange(tmp), destImage(tmp), scale / 2.0);
236 
237  recursiveSmoothX(srcImageRange(tmp), destImage(smooth), scale);
238  recursiveSmoothY(srcImageRange(smooth), destImage(smooth), scale);
239 
240  typename TMPIMG::Iterator iy = smooth.upperLeft();
241  typename TMPIMG::Iterator ty = tmp.upperLeft();
242  DestIterator dy = dul;
243 
244  const Diff2D right(1, 0);
245  const Diff2D bottom(0, 1);
246 
247 
248  TMPTYPE thresh = detail::RequiresExplicitCast<TMPTYPE>::cast((gradient_threshold * gradient_threshold) *
249  NumericTraits<TMPTYPE>::one());
250  TMPTYPE zero = NumericTraits<TMPTYPE>::zero();
251 
252  for(y=0; y<h-1; ++y, ++iy.y, ++ty.y, ++dy.y)
253  {
254  typename TMPIMG::Iterator ix = iy;
255  typename TMPIMG::Iterator tx = ty;
256  DestIterator dx = dy;
257 
258  for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, ++dx.x)
259  {
260  TMPTYPE diff = *tx - *ix;
261  TMPTYPE gx = tx[right] - *tx;
262  TMPTYPE gy = tx[bottom] - *tx;
263 
264  if((gx * gx > thresh) &&
265  (diff * (tx[right] - ix[right]) < zero))
266  {
267  if(gx < zero)
268  {
269  da.set(edge_marker, dx, right);
270  }
271  else
272  {
273  da.set(edge_marker, dx);
274  }
275  }
276  if(((gy * gy > thresh) &&
277  (diff * (tx[bottom] - ix[bottom]) < zero)))
278  {
279  if(gy < zero)
280  {
281  da.set(edge_marker, dx, bottom);
282  }
283  else
284  {
285  da.set(edge_marker, dx);
286  }
287  }
288  }
289  }
290 
291  typename TMPIMG::Iterator ix = iy;
292  typename TMPIMG::Iterator tx = ty;
293  DestIterator dx = dy;
294 
295  for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, ++dx.x)
296  {
297  TMPTYPE diff = *tx - *ix;
298  TMPTYPE gx = tx[right] - *tx;
299 
300  if((gx * gx > thresh) &&
301  (diff * (tx[right] - ix[right]) < zero))
302  {
303  if(gx < zero)
304  {
305  da.set(edge_marker, dx, right);
306  }
307  else
308  {
309  da.set(edge_marker, dx);
310  }
311  }
312  }
313 }
314 
315 template <class SrcIterator, class SrcAccessor,
316  class DestIterator, class DestAccessor,
317  class GradValue>
318 inline
320  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
321  DestIterator dul, DestAccessor da,
322  double scale, GradValue gradient_threshold)
323 {
324  differenceOfExponentialEdgeImage(sul, slr, sa, dul, da,
325  scale, gradient_threshold, 1);
326 }
327 
328 template <class SrcIterator, class SrcAccessor,
329  class DestIterator, class DestAccessor,
330  class GradValue, class DestValue>
331 inline void
332 differenceOfExponentialEdgeImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
333  pair<DestIterator, DestAccessor> dest,
334  double scale, GradValue gradient_threshold,
335  DestValue edge_marker)
336 {
337  differenceOfExponentialEdgeImage(src.first, src.second, src.third,
338  dest.first, dest.second,
339  scale, gradient_threshold, edge_marker);
340 }
341 
342 template <class SrcIterator, class SrcAccessor,
343  class DestIterator, class DestAccessor,
344  class GradValue>
345 inline void
346 differenceOfExponentialEdgeImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
347  pair<DestIterator, DestAccessor> dest,
348  double scale, GradValue gradient_threshold)
349 {
350  differenceOfExponentialEdgeImage(src.first, src.second, src.third,
351  dest.first, dest.second,
352  scale, gradient_threshold, 1);
353 }
354 
355 template <class T1, class S1,
356  class T2, class S2,
357  class GradValue, class DestValue>
358 inline void
359 differenceOfExponentialEdgeImage(MultiArrayView<2, T1, S1> const & src,
360  MultiArrayView<2, T2, S2> dest,
361  double scale,
362  GradValue gradient_threshold,
363  DestValue edge_marker)
364 {
365  vigra_precondition(src.shape() == dest.shape(),
366  "differenceOfExponentialEdgeImage(): shape mismatch between input and output.");
367  differenceOfExponentialEdgeImage(srcImageRange(src),
368  destImage(dest),
369  scale, gradient_threshold, edge_marker);
370 }
371 
372 template <class T1, class S1,
373  class T2, class S2,
374  class GradValue>
375 inline void
376 differenceOfExponentialEdgeImage(MultiArrayView<2, T1, S1> const & src,
377  MultiArrayView<2, T2, S2> dest,
378  double scale, GradValue gradient_threshold)
379 {
380  vigra_precondition(src.shape() == dest.shape(),
381  "differenceOfExponentialEdgeImage(): shape mismatch between input and output.");
382  differenceOfExponentialEdgeImage(srcImageRange(src),
383  destImage(dest),
384  scale, gradient_threshold, T2(1));
385 }
386 
387 /********************************************************/
388 /* */
389 /* differenceOfExponentialCrackEdgeImage */
390 /* */
391 /********************************************************/
392 
393 /** \brief Detect and mark edges in a crack edge image using the Shen/Castan zero-crossing detector.
394 
395  This operator applies an exponential filter to the source image
396  at the given <TT>scale</TT> and subtracts the result from the original image.
397  Zero crossings are detected in the resulting difference image. Whenever the
398  gradient at a zero crossing is greater than the given <TT>gradient_threshold</TT>,
399  an edge point is marked (using <TT>edge_marker</TT>) in the destination image
400  <i>between</i> the corresponding original pixels. Topologically, this means we
401  must insert additional pixels between the original ones to represent the
402  boundaries between the pixels (the so called zero- and one-cells, with the original
403  pixels being two-cells). Within VIGRA, such an image is called \ref CrackEdgeImage.
404  To allow insertion of the zero- and one-cells, the destination image must have twice the
405  size of the original (precisely, <TT>(2*w-1)</TT> by <TT>(2*h-1)</TT> pixels). Then the algorithm
406  proceeds as follows:
407 
408  \code
409  sign of difference image insert zero- and one-cells resulting edge points (*)
410 
411  + . - . - . * . . .
412  + - - . . . . . . * * * .
413  + + - => + . + . - => . . . * .
414  + + + . . . . . . . . * *
415  + . + . + . . . . .
416  \endcode
417 
418  Thus the edge points are marked where they actually are - in between the pixels.
419  An important property of the resulting edge image is that it conforms to the notion
420  of well-composedness as defined by Latecki et al., i.e. connected regions and edges
421  obtained by a subsequent \ref Labeling do not depend on
422  whether 4- or 8-connectivity is used.
423  The non-edge pixels (<TT>.</TT>) in the destination image remain unchanged.
424  The result conforms to the requirements of a \ref CrackEdgeImage. It can be further
425  improved by the post-processing operations \ref removeShortEdges() and
426  \ref closeGapsInCrackEdgeImage().
427 
428  The source value type (<TT>SrcAccessor::value_type</TT>) must be a linear algebra, i.e. addition,
429  subtraction and multiplication of the type with itself, and multiplication
430  with double and
431  \ref NumericTraits "NumericTraits" must
432  be defined. In addition, this type must be less-comparable.
433 
434  <b> Declarations:</b>
435 
436  pass 2D array views:
437  \code
438  namespace vigra {
439  template <class T1, class S1,
440  class T2, class S2,
441  class GradValue, class DestValue>
442  void
443  differenceOfExponentialCrackEdgeImage(MultiArrayView<2, T1, S1> const & src,
444  MultiArrayView<2, T2, S2> dest,
445  double scale,
446  GradValue gradient_threshold,
447  DestValue edge_marker = NumericTraits<DestValue>::one());
448  }
449  \endcode
450 
451  \deprecatedAPI{differenceOfExponentialCrackEdgeImage}
452  pass \ref ImageIterators and \ref DataAccessors :
453  \code
454  namespace vigra {
455  template <class SrcIterator, class SrcAccessor,
456  class DestIterator, class DestAccessor,
457  class GradValue,
458  class DestValue = DestAccessor::value_type>
459  void differenceOfExponentialCrackEdgeImage(
460  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
461  DestIterator dul, DestAccessor da,
462  double scale, GradValue gradient_threshold,
463  DestValue edge_marker = NumericTraits<DestValue>::one())
464  }
465  \endcode
466  use argument objects in conjunction with \ref ArgumentObjectFactories :
467  \code
468  namespace vigra {
469  template <class SrcIterator, class SrcAccessor,
470  class DestIterator, class DestAccessor,
471  class GradValue,
472  class DestValue = DestAccessor::value_type>
473  void differenceOfExponentialCrackEdgeImage(
474  triple<SrcIterator, SrcIterator, SrcAccessor> src,
475  pair<DestIterator, DestAccessor> dest,
476  double scale, GradValue gradient_threshold,
477  DestValue edge_marker = NumericTraits<DestValue>::one())
478  }
479  \endcode
480  \deprecatedEnd
481 
482  <b> Usage:</b>
483 
484  <b>\#include</b> <vigra/edgedetection.hxx><br>
485  Namespace: vigra
486 
487  \code
488  MultiArray<2, unsigned char> src(w,h), edges(2*w-1,2*h-1);
489  ...
490 
491  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
492  differenceOfExponentialCrackEdgeImage(src, edges,
493  0.8, 4.0, 1);
494  \endcode
495 
496  \deprecatedUsage{differenceOfExponentialCrackEdgeImage}
497  \code
498  vigra::BImage src(w,h), edges(2*w-1,2*h-1);
499 
500  // empty edge image
501  edges = 0;
502  ...
503 
504  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
505  vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges),
506  0.8, 4.0, 1);
507  \endcode
508  <b> Required Interface:</b>
509  \code
510  SrcImageIterator src_upperleft, src_lowerright;
511  DestImageIterator dest_upperleft;
512 
513  SrcAccessor src_accessor;
514  DestAccessor dest_accessor;
515 
516  SrcAccessor::value_type u = src_accessor(src_upperleft);
517  double d;
518  GradValue gradient_threshold;
519 
520  u = u + u
521  u = u - u
522  u = u * u
523  u = d * u
524  u < gradient_threshold
525 
526  DestValue edge_marker;
527  dest_accessor.set(edge_marker, dest_upperleft);
528  \endcode
529  \deprecatedEnd
530 
531  <b> Preconditions:</b>
532 
533  \code
534  scale > 0
535  gradient_threshold > 0
536  \endcode
537 
538  The destination image must have twice the size of the source:
539  \code
540  w_dest = 2 * w_src - 1
541  h_dest = 2 * h_src - 1
542  \endcode
543 */
544 doxygen_overloaded_function(template <...> void differenceOfExponentialCrackEdgeImage)
545 
546 template <class SrcIterator, class SrcAccessor,
547  class DestIterator, class DestAccessor,
548  class GradValue, class DestValue>
550  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
551  DestIterator dul, DestAccessor da,
552  double scale, GradValue gradient_threshold,
553  DestValue edge_marker)
554 {
555  vigra_precondition(scale > 0,
556  "differenceOfExponentialCrackEdgeImage(): scale > 0 required.");
557 
558  vigra_precondition(gradient_threshold > 0,
559  "differenceOfExponentialCrackEdgeImage(): "
560  "gradient_threshold > 0 required.");
561 
562  int w = slr.x - sul.x;
563  int h = slr.y - sul.y;
564  int x, y;
565 
566  typedef typename
567  NumericTraits<typename SrcAccessor::value_type>::RealPromote
568  TMPTYPE;
569  typedef BasicImage<TMPTYPE> TMPIMG;
570 
571  TMPIMG tmp(w,h);
572  TMPIMG smooth(w,h);
573 
574  TMPTYPE zero = NumericTraits<TMPTYPE>::zero();
575 
576  const Diff2D right(1,0);
577  const Diff2D bottom(0,1);
578  const Diff2D left(-1,0);
579  const Diff2D top(0,-1);
580 
581  recursiveSmoothX(srcIterRange(sul, slr, sa), destImage(tmp), scale / 2.0);
582  recursiveSmoothY(srcImageRange(tmp), destImage(tmp), scale / 2.0);
583 
584  recursiveSmoothX(srcImageRange(tmp), destImage(smooth), scale);
585  recursiveSmoothY(srcImageRange(smooth), destImage(smooth), scale);
586 
587  typename TMPIMG::Iterator iy = smooth.upperLeft();
588  typename TMPIMG::Iterator ty = tmp.upperLeft();
589  DestIterator dy = dul;
590 
591  TMPTYPE thresh = detail::RequiresExplicitCast<TMPTYPE>::cast((gradient_threshold * gradient_threshold) *
592  NumericTraits<TMPTYPE>::one());
593 
594  // find zero crossings above threshold
595  for(y=0; y<h-1; ++y, ++iy.y, ++ty.y, dy.y+=2)
596  {
597  typename TMPIMG::Iterator ix = iy;
598  typename TMPIMG::Iterator tx = ty;
599  DestIterator dx = dy;
600 
601  for(int x=0; x<w-1; ++x, ++ix.x, ++tx.x, dx.x+=2)
602  {
603  TMPTYPE diff = *tx - *ix;
604  TMPTYPE gx = tx[right] - *tx;
605  TMPTYPE gy = tx[bottom] - *tx;
606 
607  if((gx * gx > thresh) &&
608  (diff * (tx[right] - ix[right]) < zero))
609  {
610  da.set(edge_marker, dx, right);
611  }
612  if((gy * gy > thresh) &&
613  (diff * (tx[bottom] - ix[bottom]) < zero))
614  {
615  da.set(edge_marker, dx, bottom);
616  }
617  }
618 
619  TMPTYPE diff = *tx - *ix;
620  TMPTYPE gy = tx[bottom] - *tx;
621 
622  if((gy * gy > thresh) &&
623  (diff * (tx[bottom] - ix[bottom]) < zero))
624  {
625  da.set(edge_marker, dx, bottom);
626  }
627  }
628 
629  typename TMPIMG::Iterator ix = iy;
630  typename TMPIMG::Iterator tx = ty;
631  DestIterator dx = dy;
632 
633  for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, dx.x+=2)
634  {
635  TMPTYPE diff = *tx - *ix;
636  TMPTYPE gx = tx[right] - *tx;
637 
638  if((gx * gx > thresh) &&
639  (diff * (tx[right] - ix[right]) < zero))
640  {
641  da.set(edge_marker, dx, right);
642  }
643  }
644 
645  iy = smooth.upperLeft() + Diff2D(0,1);
646  ty = tmp.upperLeft() + Diff2D(0,1);
647  dy = dul + Diff2D(1,2);
648 
649  const Diff2D topleft(-1,-1);
650  const Diff2D topright(1,-1);
651  const Diff2D bottomleft(-1,1);
652  const Diff2D bottomright(1,1);
653 
654  // find missing 1-cells below threshold (x-direction)
655  for(y=0; y<h-2; ++y, ++iy.y, ++ty.y, dy.y+=2)
656  {
657  typename TMPIMG::Iterator ix = iy;
658  typename TMPIMG::Iterator tx = ty;
659  DestIterator dx = dy;
660 
661  for(int x=0; x<w-2; ++x, ++ix.x, ++tx.x, dx.x+=2)
662  {
663  if(da(dx) == edge_marker) continue;
664 
665  TMPTYPE diff = *tx - *ix;
666 
667  if((diff * (tx[right] - ix[right]) < zero) &&
668  (((da(dx, bottomright) == edge_marker) &&
669  (da(dx, topleft) == edge_marker)) ||
670  ((da(dx, bottomleft) == edge_marker) &&
671  (da(dx, topright) == edge_marker))))
672 
673  {
674  da.set(edge_marker, dx);
675  }
676  }
677  }
678 
679  iy = smooth.upperLeft() + Diff2D(1,0);
680  ty = tmp.upperLeft() + Diff2D(1,0);
681  dy = dul + Diff2D(2,1);
682 
683  // find missing 1-cells below threshold (y-direction)
684  for(y=0; y<h-2; ++y, ++iy.y, ++ty.y, dy.y+=2)
685  {
686  typename TMPIMG::Iterator ix = iy;
687  typename TMPIMG::Iterator tx = ty;
688  DestIterator dx = dy;
689 
690  for(int x=0; x<w-2; ++x, ++ix.x, ++tx.x, dx.x+=2)
691  {
692  if(da(dx) == edge_marker) continue;
693 
694  TMPTYPE diff = *tx - *ix;
695 
696  if((diff * (tx[bottom] - ix[bottom]) < zero) &&
697  (((da(dx, bottomright) == edge_marker) &&
698  (da(dx, topleft) == edge_marker)) ||
699  ((da(dx, bottomleft) == edge_marker) &&
700  (da(dx, topright) == edge_marker))))
701 
702  {
703  da.set(edge_marker, dx);
704  }
705  }
706  }
707 
708  dy = dul + Diff2D(1,1);
709 
710  // find missing 0-cells
711  for(y=0; y<h-1; ++y, dy.y+=2)
712  {
713  DestIterator dx = dy;
714 
715  for(int x=0; x<w-1; ++x, dx.x+=2)
716  {
717  const Diff2D dist[] = {right, top, left, bottom };
718 
719  int i;
720  for(i=0; i<4; ++i)
721  {
722  if(da(dx, dist[i]) == edge_marker) break;
723  }
724 
725  if(i < 4) da.set(edge_marker, dx);
726  }
727  }
728 }
729 
730 template <class SrcIterator, class SrcAccessor,
731  class DestIterator, class DestAccessor,
732  class GradValue, class DestValue>
733 inline void
734 differenceOfExponentialCrackEdgeImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
735  pair<DestIterator, DestAccessor> dest,
736  double scale, GradValue gradient_threshold,
737  DestValue edge_marker)
738 {
739  differenceOfExponentialCrackEdgeImage(src.first, src.second, src.third,
740  dest.first, dest.second,
741  scale, gradient_threshold, edge_marker);
742 }
743 
744 template <class T1, class S1,
745  class T2, class S2,
746  class GradValue, class DestValue>
747 inline void
748 differenceOfExponentialCrackEdgeImage(MultiArrayView<2, T1, S1> const & src,
749  MultiArrayView<2, T2, S2> dest,
750  double scale,
751  GradValue gradient_threshold,
752  DestValue edge_marker)
753 {
754  vigra_precondition(2*src.shape() - Shape2(1) == dest.shape(),
755  "differenceOfExponentialCrackEdgeImage(): shape mismatch between input and output.");
756  differenceOfExponentialCrackEdgeImage(srcImageRange(src),
757  destImage(dest),
758  scale, gradient_threshold, edge_marker);
759 }
760 
761 /********************************************************/
762 /* */
763 /* removeShortEdges */
764 /* */
765 /********************************************************/
766 
767 /** \brief Remove short edges from an edge image.
768 
769  This algorithm can be applied as a post-processing operation of
770  \ref differenceOfExponentialEdgeImage() and \ref differenceOfExponentialCrackEdgeImage().
771  It removes all edges that are shorter than <TT>min_edge_length</TT>. The corresponding
772  pixels are set to the <TT>non_edge_marker</TT>. The idea behind this algorithms is
773  that very short edges are probably caused by noise and don't represent interesting
774  image structure. Technically, the algorithms executes a connected components labeling,
775  so the image's value type must be equality comparable.
776 
777  If the source image fulfills the requirements of a \ref CrackEdgeImage,
778  it will still do so after application of this algorithm.
779 
780  Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place,
781  i.e. on only one image. Also, the algorithm assumes that all non-edges pixels are already
782  marked with the given <TT>non_edge_marker</TT> value.
783 
784  <b> Declarations:</b>
785 
786  pass 2D array views:
787  \code
788  namespace vigra {
789  template <class T, class S, class Value>
790  void
791  removeShortEdges(MultiArrayView<2, T, S> image,
792  unsigned int min_edge_length, Value non_edge_marker);
793  }
794  \endcode
795 
796  \deprecatedAPI{removeShortEdges}
797  pass \ref ImageIterators and \ref DataAccessors :
798  \code
799  namespace vigra {
800  template <class Iterator, class Accessor, class SrcValue>
801  void removeShortEdges(
802  Iterator sul, Iterator slr, Accessor sa,
803  int min_edge_length, SrcValue non_edge_marker)
804  }
805  \endcode
806  use argument objects in conjunction with \ref ArgumentObjectFactories :
807  \code
808  namespace vigra {
809  template <class Iterator, class Accessor, class SrcValue>
810  void removeShortEdges(
811  triple<Iterator, Iterator, Accessor> src,
812  int min_edge_length, SrcValue non_edge_marker)
813  }
814  \endcode
815  \deprecatedEnd
816 
817  <b> Usage:</b>
818 
819  <b>\#include</b> <vigra/edgedetection.hxx><br>
820  Namespace: vigra
821 
822  \code
823  MultiArray<2, unsigned char> src(w,h), edges(w,h);
824  ...
825 
826  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
827  differenceOfExponentialEdgeImage(src, edges,
828  0.8, 4.0, 1);
829 
830  // zero edges shorter than 10 pixels
831  removeShortEdges(edges, 10, 0);
832  \endcode
833 
834  \deprecatedUsage{removeShortEdges}
835  \code
836  vigra::BImage src(w,h), edges(w,h);
837 
838  // empty edge image
839  edges = 0;
840  ...
841 
842  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
843  vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges),
844  0.8, 4.0, 1);
845 
846  // zero edges shorter than 10 pixels
847  vigra::removeShortEdges(srcImageRange(edges), 10, 0);
848  \endcode
849  <b> Required Interface:</b>
850  \code
851  SrcImageIterator src_upperleft, src_lowerright;
852  DestImageIterator dest_upperleft;
853 
854  SrcAccessor src_accessor;
855  DestAccessor dest_accessor;
856 
857  SrcAccessor::value_type u = src_accessor(src_upperleft);
858 
859  u == u
860 
861  SrcValue non_edge_marker;
862  src_accessor.set(non_edge_marker, src_upperleft);
863  \endcode
864  \deprecatedEnd
865 */
866 doxygen_overloaded_function(template <...> void removeShortEdges)
867 
868 template <class Iterator, class Accessor, class Value>
869 void removeShortEdges(
870  Iterator sul, Iterator slr, Accessor sa,
871  unsigned int min_edge_length, Value non_edge_marker)
872 {
873  int w = slr.x - sul.x;
874  int h = slr.y - sul.y;
875  int x,y;
876 
877  IImage labels(w, h);
878  labels = 0;
879 
880  int number_of_regions =
881  labelImageWithBackground(srcIterRange(sul,slr,sa),
882  destImage(labels), true, non_edge_marker);
883 
884  ArrayOfRegionStatistics<FindROISize<int> >
885  region_stats(number_of_regions);
886 
887  inspectTwoImages(srcImageRange(labels), srcImage(labels), region_stats);
888 
889  IImage::Iterator ly = labels.upperLeft();
890  Iterator oy = sul;
891 
892  for(y=0; y<h; ++y, ++oy.y, ++ly.y)
893  {
894  Iterator ox(oy);
895  IImage::Iterator lx(ly);
896 
897  for(x=0; x<w; ++x, ++ox.x, ++lx.x)
898  {
899  if(sa(ox) == non_edge_marker) continue;
900  if((region_stats[*lx].count) < min_edge_length)
901  {
902  sa.set(non_edge_marker, ox);
903  }
904  }
905  }
906 }
907 
908 template <class Iterator, class Accessor, class Value>
909 inline void
910 removeShortEdges(triple<Iterator, Iterator, Accessor> src,
911  unsigned int min_edge_length, Value non_edge_marker)
912 {
913  removeShortEdges(src.first, src.second, src.third,
914  min_edge_length, non_edge_marker);
915 }
916 
917 template <class T, class S, class Value>
918 inline void
919 removeShortEdges(MultiArrayView<2, T, S> image,
920  unsigned int min_edge_length, Value non_edge_marker)
921 {
922  removeShortEdges(destImageRange(image),
923  min_edge_length, non_edge_marker);
924 }
925 
926 /********************************************************/
927 /* */
928 /* closeGapsInCrackEdgeImage */
929 /* */
930 /********************************************************/
931 
932 /** \brief Close one-pixel wide gaps in a cell grid edge image.
933 
934  This algorithm is typically applied as a post-processing operation of
935  \ref differenceOfExponentialCrackEdgeImage(). The source image must fulfill
936  the requirements of a \ref CrackEdgeImage, and will still do so after
937  application of this algorithm.
938 
939  It closes one pixel wide gaps in the edges resulting from this algorithm.
940  Since these gaps are usually caused by zero crossing slightly below the gradient
941  threshold used in edge detection, this algorithms acts like a weak hysteresis
942  thresholding. The newly found edge pixels are marked with the given <TT>edge_marker</TT>.
943  The image's value type must be equality comparable.
944 
945  Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place,
946  i.e. on only one image.
947 
948  <b> Declarations:</b>
949 
950  pass 2D array views:
951  \code
952  namespace vigra {
953  template <class T, class S, class Value>
954  void
955  closeGapsInCrackEdgeImage(MultiArrayView<2, T, S> image, Value edge_marker);
956  }
957  \endcode
958 
959  \deprecatedAPI{closeGapsInCrackEdgeImage}
960  pass \ref ImageIterators and \ref DataAccessors :
961  \code
962  namespace vigra {
963  template <class SrcIterator, class SrcAccessor, class SrcValue>
964  void closeGapsInCrackEdgeImage(
965  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
966  SrcValue edge_marker)
967  }
968  \endcode
969  use argument objects in conjunction with \ref ArgumentObjectFactories :
970  \code
971  namespace vigra {
972  template <class SrcIterator, class SrcAccessor, class SrcValue>
973  void closeGapsInCrackEdgeImage(
974  triple<SrcIterator, SrcIterator, SrcAccessor> src,
975  SrcValue edge_marker)
976  }
977  \endcode
978  \deprecatedEnd
979 
980  <b> Usage:</b>
981 
982  <b>\#include</b> <vigra/edgedetection.hxx><br>
983  Namespace: vigra
984 
985  \code
986  MultiArray<2, unsigned char> src(w,h), edges(2*w-1, 2*h-1);
987  ...
988 
989  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
990  differenceOfExponentialCrackEdgeImage(src, edges,
991  0.8, 4.0, 1);
992 
993  // close gaps, mark with 1
994  closeGapsInCrackEdgeImage(edges, 1);
995 
996  // zero edges shorter than 20 pixels
997  removeShortEdges(edges, 10, 0);
998  \endcode
999 
1000  \deprecatedUsage{closeGapsInCrackEdgeImage}
1001  \code
1002  vigra::BImage src(w,h), edges(2*w-1, 2*h-1);
1003 
1004  // empty edge image
1005  edges = 0;
1006  ...
1007 
1008  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
1009  vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges),
1010  0.8, 4.0, 1);
1011 
1012  // close gaps, mark with 1
1013  vigra::closeGapsInCrackEdgeImage(srcImageRange(edges), 1);
1014 
1015  // zero edges shorter than 20 pixels
1016  vigra::removeShortEdges(srcImageRange(edges), 10, 0);
1017  \endcode
1018  <b> Required Interface:</b>
1019  \code
1020  SrcImageIterator src_upperleft, src_lowerright;
1021 
1022  SrcAccessor src_accessor;
1023  DestAccessor dest_accessor;
1024 
1025  SrcAccessor::value_type u = src_accessor(src_upperleft);
1026 
1027  u == u
1028  u != u
1029 
1030  SrcValue edge_marker;
1031  src_accessor.set(edge_marker, src_upperleft);
1032  \endcode
1033  \deprecatedEnd
1034 */
1035 doxygen_overloaded_function(template <...> void closeGapsInCrackEdgeImage)
1036 
1037 template <class SrcIterator, class SrcAccessor, class SrcValue>
1039  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1040  SrcValue edge_marker)
1041 {
1042  int w = slr.x - sul.x;
1043  int h = slr.y - sul.y;
1044 
1045  vigra_precondition(w % 2 == 1 && h % 2 == 1,
1046  "closeGapsInCrackEdgeImage(): Input is not a crack edge image (must have odd-numbered shape).");
1047 
1048  int w2 = w / 2, h2 = h / 2, x, y;
1049 
1050  int count1, count2, count3;
1051 
1052  const Diff2D right(1,0);
1053  const Diff2D bottom(0,1);
1054  const Diff2D left(-1,0);
1055  const Diff2D top(0,-1);
1056 
1057  const Diff2D leftdist[] = { Diff2D(0, 0), Diff2D(-1, 1), Diff2D(-2, 0), Diff2D(-1, -1)};
1058  const Diff2D rightdist[] = { Diff2D(2, 0), Diff2D(1, 1), Diff2D(0, 0), Diff2D(1, -1)};
1059  const Diff2D topdist[] = { Diff2D(1, -1), Diff2D(0, 0), Diff2D(-1, -1), Diff2D(0, -2)};
1060  const Diff2D bottomdist[] = { Diff2D(1, 1), Diff2D(0, 2), Diff2D(-1, 1), Diff2D(0, 0)};
1061 
1062  int i;
1063 
1064  SrcIterator sy = sul + Diff2D(0,1);
1065  SrcIterator sx;
1066 
1067  // close 1-pixel wide gaps (x-direction)
1068  for(y=0; y<h2; ++y, sy.y+=2)
1069  {
1070  sx = sy + Diff2D(2,0);
1071 
1072  for(x=2; x<w2; ++x, sx.x+=2)
1073  {
1074  if(sa(sx) == edge_marker) continue;
1075 
1076  if(sa(sx, left) != edge_marker) continue;
1077  if(sa(sx, right) != edge_marker) continue;
1078 
1079  count1 = 0;
1080  count2 = 0;
1081  count3 = 0;
1082 
1083  for(i=0; i<4; ++i)
1084  {
1085  if(sa(sx, leftdist[i]) == edge_marker)
1086  {
1087  ++count1;
1088  count3 ^= 1 << i;
1089  }
1090  if(sa(sx, rightdist[i]) == edge_marker)
1091  {
1092  ++count2;
1093  count3 ^= 1 << i;
1094  }
1095  }
1096 
1097  if(count1 <= 1 || count2 <= 1 || count3 == 15)
1098  {
1099  sa.set(edge_marker, sx);
1100  }
1101  }
1102  }
1103 
1104  sy = sul + Diff2D(1,2);
1105 
1106  // close 1-pixel wide gaps (y-direction)
1107  for(y=2; y<h2; ++y, sy.y+=2)
1108  {
1109  sx = sy;
1110 
1111  for(x=0; x<w2; ++x, sx.x+=2)
1112  {
1113  if(sa(sx) == edge_marker) continue;
1114 
1115  if(sa(sx, top) != edge_marker) continue;
1116  if(sa(sx, bottom) != edge_marker) continue;
1117 
1118  count1 = 0;
1119  count2 = 0;
1120  count3 = 0;
1121 
1122  for(i=0; i<4; ++i)
1123  {
1124  if(sa(sx, topdist[i]) == edge_marker)
1125  {
1126  ++count1;
1127  count3 ^= 1 << i;
1128  }
1129  if(sa(sx, bottomdist[i]) == edge_marker)
1130  {
1131  ++count2;
1132  count3 ^= 1 << i;
1133  }
1134  }
1135 
1136  if(count1 <= 1 || count2 <= 1 || count3 == 15)
1137  {
1138  sa.set(edge_marker, sx);
1139  }
1140  }
1141  }
1142 }
1143 
1144 template <class SrcIterator, class SrcAccessor, class SrcValue>
1145 inline void
1146 closeGapsInCrackEdgeImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1147  SrcValue edge_marker)
1148 {
1149  closeGapsInCrackEdgeImage(src.first, src.second, src.third,
1150  edge_marker);
1151 }
1152 
1153 template <class T, class S, class Value>
1154 inline void
1155 closeGapsInCrackEdgeImage(MultiArrayView<2, T, S> image, Value edge_marker)
1156 {
1157  closeGapsInCrackEdgeImage(destImageRange(image), edge_marker);
1158 }
1159 
1160 /********************************************************/
1161 /* */
1162 /* beautifyCrackEdgeImage */
1163 /* */
1164 /********************************************************/
1165 
1166 /** \brief Beautify crack edge image for visualization.
1167 
1168  This algorithm is applied as a post-processing operation of
1169  \ref differenceOfExponentialCrackEdgeImage(). The source image must fulfill
1170  the requirements of a \ref CrackEdgeImage, but will <b> not</b> do so after
1171  application of this algorithm. In particular, the algorithm removes zero-cells
1172  marked as edges to avoid staircase effects on diagonal lines like this:
1173 
1174  \code
1175  original edge points (*) resulting edge points
1176 
1177  . * . . . . * . . .
1178  . * * * . . . * . .
1179  . . . * . => . . . * .
1180  . . . * * . . . . *
1181  . . . . . . . . . .
1182  \endcode
1183 
1184  Therefore, this algorithm should only be applied as a visualization aid, i.e.
1185  for human inspection. The algorithm assumes that edges are marked with <TT>edge_marker</TT>,
1186  and background pixels with <TT>background_marker</TT>. The image's value type must be
1187  equality comparable.
1188 
1189  Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place,
1190  i.e. on only one image.
1191 
1192  <b> Declarations:</b>
1193 
1194  pass 2D array views:
1195  \code
1196  namespace vigra {
1197  template <class T, class S, class Value>
1198  void
1199  beautifyCrackEdgeImage(MultiArrayView<2, T, S> image,
1200  Value edge_marker, Value background_marker);
1201  }
1202  \endcode
1203 
1204  \deprecatedAPI{beautifyCrackEdgeImage}
1205  pass \ref ImageIterators and \ref DataAccessors :
1206  \code
1207  namespace vigra {
1208  template <class SrcIterator, class SrcAccessor, class SrcValue>
1209  void beautifyCrackEdgeImage(
1210  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1211  SrcValue edge_marker, SrcValue background_marker)
1212  }
1213  \endcode
1214  use argument objects in conjunction with \ref ArgumentObjectFactories :
1215  \code
1216  namespace vigra {
1217  template <class SrcIterator, class SrcAccessor, class SrcValue>
1218  void beautifyCrackEdgeImage(
1219  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1220  SrcValue edge_marker, SrcValue background_marker)
1221  }
1222  \endcode
1223  \deprecatedEnd
1224 
1225  <b> Usage:</b>
1226 
1227  <b>\#include</b> <vigra/edgedetection.hxx><br>
1228  Namespace: vigra
1229 
1230  \code
1231  MultiArray<2, unsigned char> src(w,h), edges(2*w-1, 2*h-1);
1232  ...
1233 
1234  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
1235  differenceOfExponentialCrackEdgeImage(src, edges,
1236  0.8, 4.0, 1);
1237 
1238  // beautify edge image for visualization
1239  beautifyCrackEdgeImage(edges, 1, 0);
1240 
1241  // show to the user ('window' is an unspecified GUI widget to display VIGRA images)
1242  window.open(edges);
1243  \endcode
1244 
1245  \deprecatedUsage{beautifyCrackEdgeImage}
1246  \code
1247  vigra::BImage src(w,h), edges(2*w-1, 2*h-1);
1248 
1249  // empty edge image
1250  edges = 0;
1251  ...
1252 
1253  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
1254  vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges),
1255  0.8, 4.0, 1);
1256 
1257  // beautify edge image for visualization
1258  vigra::beautifyCrackEdgeImage(destImageRange(edges), 1, 0);
1259 
1260  // show to the user
1261  window.open(edges);
1262  \endcode
1263  <b> Required Interface:</b>
1264  \code
1265  SrcImageIterator src_upperleft, src_lowerright;
1266 
1267  SrcAccessor src_accessor;
1268  DestAccessor dest_accessor;
1269 
1270  SrcAccessor::value_type u = src_accessor(src_upperleft);
1271 
1272  u == u
1273  u != u
1274 
1275  SrcValue background_marker;
1276  src_accessor.set(background_marker, src_upperleft);
1277  \endcode
1278  \deprecatedEnd
1279 */
1280 doxygen_overloaded_function(template <...> void beautifyCrackEdgeImage)
1281 
1282 template <class SrcIterator, class SrcAccessor, class SrcValue>
1284  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1285  SrcValue edge_marker, SrcValue background_marker)
1286 {
1287  int w = slr.x - sul.x;
1288  int h = slr.y - sul.y;
1289 
1290  vigra_precondition(w % 2 == 1 && h % 2 == 1,
1291  "beautifyCrackEdgeImage(): Input is not a crack edge image (must have odd-numbered shape).");
1292 
1293  int w2 = w / 2, h2 = h / 2, x, y;
1294 
1295  SrcIterator sy = sul + Diff2D(1,1);
1296  SrcIterator sx;
1297 
1298  const Diff2D right(1,0);
1299  const Diff2D bottom(0,1);
1300  const Diff2D left(-1,0);
1301  const Diff2D top(0,-1);
1302 
1303  // delete 0-cells at corners
1304  for(y=0; y<h2; ++y, sy.y+=2)
1305  {
1306  sx = sy;
1307 
1308  for(x=0; x<w2; ++x, sx.x+=2)
1309  {
1310  if(sa(sx) != edge_marker) continue;
1311 
1312  if(sa(sx, right) == edge_marker && sa(sx, left) == edge_marker) continue;
1313  if(sa(sx, bottom) == edge_marker && sa(sx, top) == edge_marker) continue;
1314 
1315  sa.set(background_marker, sx);
1316  }
1317  }
1318 }
1319 
1320 template <class SrcIterator, class SrcAccessor, class SrcValue>
1321 inline void
1322 beautifyCrackEdgeImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1323  SrcValue edge_marker, SrcValue background_marker)
1324 {
1325  beautifyCrackEdgeImage(src.first, src.second, src.third,
1326  edge_marker, background_marker);
1327 }
1328 
1329 template <class T, class S, class Value>
1330 inline void
1331 beautifyCrackEdgeImage(MultiArrayView<2, T, S> image,
1332  Value edge_marker, Value background_marker)
1333 {
1334  beautifyCrackEdgeImage(destImageRange(image),
1335  edge_marker, background_marker);
1336 }
1337 
1338 
1339 /** Helper class that stores edgel attributes.
1340 */
1341 class Edgel
1342 {
1343  public:
1344 
1345  /** The type of an Edgel's members.
1346  */
1347  typedef float value_type;
1348 
1349  /** The edgel's sub-pixel x coordinate.
1350  */
1352 
1353  /** The edgel's sub-pixel y coordinate.
1354  */
1356 
1357  /** The edgel's strength (magnitude of the gradient vector).
1358  */
1360 
1361  /** \brief The edgel's orientation.
1362 
1363  This is the clockwise angle in radians
1364  between the x-axis and the edge, so that the bright side of the
1365  edge is on the left when one looks along the orientation vector.
1366  The angle is measured clockwise because the y-axis increases
1367  downwards (left-handed coordinate system):
1368 
1369  \code
1370 
1371  edgel axis
1372  \
1373  (dark \ (bright side)
1374  side) \
1375  \
1376  +------------> x-axis
1377  |\ |
1378  | \ /_/ orientation angle
1379  | \\
1380  | \
1381  |
1382  y-axis V
1383  \endcode
1384 
1385  So, for example a vertical edge with its dark side on the left
1386  has orientation PI/2, and a horizontal edge with dark side on top
1387  has orientation PI. Obviously, the edge's orientation changes
1388  by PI if the contrast is reversed.
1389 
1390  Note that this convention changed as of VIGRA version 1.7.0.
1391 
1392  */
1394 
1395  Edgel()
1396  : x(0.0), y(0.0), strength(0.0), orientation(0.0)
1397  {}
1398 
1400  : x(ix), y(iy), strength(is), orientation(io)
1401  {}
1402 };
1403 
1404 template <class SrcIterator, class SrcAccessor,
1405  class MagnitudeImage, class BackInsertable, class GradValue>
1406 void internalCannyFindEdgels(SrcIterator ul, SrcAccessor grad,
1407  MagnitudeImage const & magnitude,
1408  BackInsertable & edgels, GradValue grad_thresh)
1409 {
1410  typedef typename SrcAccessor::value_type PixelType;
1411  typedef typename PixelType::value_type ValueType;
1412 
1413  vigra_precondition(grad_thresh >= NumericTraits<GradValue>::zero(),
1414  "cannyFindEdgels(): gradient threshold must not be negative.");
1415 
1416  double t = 0.5 / VIGRA_CSTD::sin(M_PI/8.0);
1417 
1418  ul += Diff2D(1,1);
1419  for(int y=1; y<magnitude.height()-1; ++y, ++ul.y)
1420  {
1421  SrcIterator ix = ul;
1422  for(int x=1; x<magnitude.width()-1; ++x, ++ix.x)
1423  {
1424  double mag = magnitude(x, y);
1425  if(mag <= grad_thresh)
1426  continue;
1427  ValueType gradx = grad.getComponent(ix, 0);
1428  ValueType grady = grad.getComponent(ix, 1);
1429 
1430  int dx = (int)VIGRA_CSTD::floor(gradx*t/mag + 0.5);
1431  int dy = (int)VIGRA_CSTD::floor(grady*t/mag + 0.5);
1432 
1433  int x1 = x - dx,
1434  x2 = x + dx,
1435  y1 = y - dy,
1436  y2 = y + dy;
1437 
1438  double m1 = magnitude(x1, y1);
1439  double m3 = magnitude(x2, y2);
1440 
1441  if(m1 < mag && m3 <= mag)
1442  {
1443  Edgel edgel;
1444 
1445  // local maximum => quadratic interpolation of sub-pixel location
1446  double del = 0.5 * (m1 - m3) / (m1 + m3 - 2.0*mag);
1447  edgel.x = Edgel::value_type(x + dx*del);
1448  edgel.y = Edgel::value_type(y + dy*del);
1449  edgel.strength = Edgel::value_type(mag);
1450  double orientation = VIGRA_CSTD::atan2(grady, gradx) + 0.5*M_PI;
1451  if(orientation < 0.0)
1452  orientation += 2.0*M_PI;
1453  edgel.orientation = Edgel::value_type(orientation);
1454  edgels.push_back(edgel);
1455  }
1456  }
1457  }
1458 }
1459 
1460 /********************************************************/
1461 /* */
1462 /* cannyEdgelList */
1463 /* */
1464 /********************************************************/
1465 
1466 /** \brief Simple implementation of Canny's edge detector.
1467 
1468  The function can be called in two modes: If you pass a 'scale', it is assumed that the
1469  original image is scalar, and the Gaussian gradient is internally computed at the
1470  given 'scale'. If the function is called without scale parameter, it is assumed that
1471  the given image already contains the gradient (i.e. its value_type must be
1472  a vector of length 2).
1473 
1474  On the basis of the gradient image, a simple non-maxima suppression is performed:
1475  for each 3x3 neighborhood, it is determined whether the center pixel has
1476  larger gradient magnitude than its two neighbors in gradient direction
1477  (where the direction is rounded into octants). If this is the case,
1478  a new \ref Edgel is appended to the given vector of <TT>edgels</TT>. The subpixel
1479  edgel position is determined by fitting a parabola to the three gradient
1480  magnitude values mentioned above. The sub-pixel location of the parabola's tip
1481  and the gradient magnitude and direction (from the pixel center)
1482  are written in the newly created edgel.
1483 
1484  <b> Declarations:</b>
1485 
1486  pass 2D array views:
1487  \code
1488  namespace vigra {
1489  // compute edgels from a scalar image (determine gradient internally at 'scale')
1490  template <class T, class S, class BackInsertable>
1491  void
1492  cannyEdgelList(MultiArrayView<2, T, S> const & src,
1493  BackInsertable & edgels,
1494  double scale);
1495 
1496  // compute edgels from a pre-computed gradient image
1497  template <class T, class S, class BackInsertable>
1498  void
1499  cannyEdgelList(MultiArrayView<2, TinyVector<T, 2>, S> const & src,
1500  BackInsertable & edgels);
1501  }
1502  \endcode
1503 
1504  \deprecatedAPI{cannyEdgelList}
1505  pass \ref ImageIterators and \ref DataAccessors :
1506  \code
1507  namespace vigra {
1508  // compute edgels from a gradient image
1509  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1510  void
1511  cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1512  BackInsertable & edgels);
1513 
1514  // compute edgels from a scalar image (determine gradient internally at 'scale')
1515  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1516  void
1517  cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1518  BackInsertable & edgels, double scale);
1519  }
1520  \endcode
1521  use argument objects in conjunction with \ref ArgumentObjectFactories :
1522  \code
1523  namespace vigra {
1524  // compute edgels from a gradient image
1525  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1526  void
1527  cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1528  BackInsertable & edgels);
1529 
1530  // compute edgels from a scalar image (determine gradient internally at 'scale')
1531  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1532  void
1533  cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1534  BackInsertable & edgels, double scale);
1535  }
1536  \endcode
1537  \deprecatedEnd
1538 
1539  <b> Usage:</b>
1540 
1541  <b>\#include</b> <vigra/edgedetection.hxx><br>
1542  Namespace: vigra
1543 
1544  \code
1545  MultiArray<2, unsigned char> src(w,h);
1546 
1547  // create empty edgel list
1548  std::vector<vigra::Edgel> edgels;
1549  ...
1550 
1551  // find edgels at scale 0.8
1552  cannyEdgelList(src, edgels, 0.8);
1553  \endcode
1554 
1555  \deprecatedUsage{cannyEdgelList}
1556  \code
1557  vigra::BImage src(w,h);
1558 
1559  // empty edgel list
1560  std::vector<vigra::Edgel> edgels;
1561  ...
1562 
1563  // find edgels at scale 0.8
1564  vigra::cannyEdgelList(srcImageRange(src), edgels, 0.8);
1565  \endcode
1566  <b> Required Interface:</b>
1567  \code
1568  SrcImageIterator src_upperleft;
1569  SrcAccessor src_accessor;
1570 
1571  src_accessor(src_upperleft);
1572 
1573  BackInsertable edgels;
1574  edgels.push_back(Edgel());
1575  \endcode
1576  SrcAccessor::value_type must be a type convertible to float
1577  \deprecatedEnd
1578 
1579  <b> Preconditions:</b>
1580 
1581  \code
1582  scale > 0
1583  \endcode
1584 */
1585 doxygen_overloaded_function(template <...> void cannyEdgelList)
1586 
1587 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1588 void
1589 cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1590  BackInsertable & edgels, double scale)
1591 {
1592  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
1593  BasicImage<TinyVector<TmpType, 2> > grad(lr-ul);
1594  gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale);
1595 
1596  cannyEdgelList(srcImageRange(grad), edgels);
1597 }
1598 
1599 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1600 void
1601 cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1602  BackInsertable & edgels)
1603 {
1604  using namespace functor;
1605 
1606  typedef typename SrcAccessor::value_type SrcType;
1607  typedef typename NumericTraits<typename SrcType::value_type>::RealPromote TmpType;
1608  BasicImage<TmpType> magnitude(lr-ul);
1609  transformImage(srcIterRange(ul, lr, src), destImage(magnitude), norm(Arg1()));
1610 
1611  // find edgels
1612  internalCannyFindEdgels(ul, src, magnitude, edgels, NumericTraits<TmpType>::zero());
1613 }
1614 
1615 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1616 inline void
1617 cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1618  BackInsertable & edgels, double scale)
1619 {
1620  cannyEdgelList(src.first, src.second, src.third, edgels, scale);
1621 }
1622 
1623 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1624 inline void
1625 cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1626  BackInsertable & edgels)
1627 {
1628  cannyEdgelList(src.first, src.second, src.third, edgels);
1629 }
1630 
1631 template <class T, class S, class BackInsertable>
1632 inline void
1633 cannyEdgelList(MultiArrayView<2, T, S> const & src,
1634  BackInsertable & edgels, double scale)
1635 {
1636  cannyEdgelList(srcImageRange(src), edgels, scale);
1637 }
1638 
1639 template <class T, class S, class BackInsertable>
1640 inline void
1641 cannyEdgelList(MultiArrayView<2, TinyVector<T, 2>, S> const & src,
1642  BackInsertable & edgels)
1643 {
1644  cannyEdgelList(srcImageRange(src), edgels);
1645 }
1646 
1647 /********************************************************/
1648 /* */
1649 /* cannyEdgelListThreshold */
1650 /* */
1651 /********************************************************/
1652 
1653 /** \brief Canny's edge detector with thresholding.
1654 
1655  This function works exactly like \ref cannyEdgelList(), but
1656  you also pass a threshold for the minimal gradient magnitude,
1657  so that edgels whose strength is below the threshold are not
1658  inserted into the edgel list.
1659 
1660  <b> Declarations:</b>
1661 
1662  pass 2D array views:
1663  \code
1664  namespace vigra {
1665  // compute edgels from a scalar image (determine gradient internally at 'scale')
1666  template <class T, class S,
1667  class BackInsertable, class GradValue>
1668  void
1669  cannyEdgelListThreshold(MultiArrayView<2, T, S> const & src,
1670  BackInsertable & edgels,
1671  double scale,
1672  GradValue grad_threshold);
1673 
1674  // compute edgels from a pre-computed gradient image
1675  template <class T, class S,
1676  class BackInsertable, class GradValue>
1677  void
1678  cannyEdgelListThreshold(MultiArrayView<2, TinyVector<T, 2>, S> const & src,
1679  BackInsertable & edgels,
1680  GradValue grad_threshold);
1681  }
1682  \endcode
1683 
1684  \deprecatedAPI{cannyEdgelListThreshold}
1685  pass \ref ImageIterators and \ref DataAccessors :
1686  \code
1687  namespace vigra {
1688  // compute edgels from a gradient image
1689  template <class SrcIterator, class SrcAccessor,
1690  class BackInsertable, class GradValue>
1691  void
1692  cannyEdgelListThreshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1693  BackInsertable & edgels, GradValue grad_threshold);
1694 
1695  // compute edgels from a scalar image (determine gradient internally at 'scale')
1696  template <class SrcIterator, class SrcAccessor,
1697  class BackInsertable, class GradValue>
1698  void
1699  cannyEdgelListThreshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1700  BackInsertable & edgels, double scale, GradValue grad_threshold);
1701  }
1702  \endcode
1703  use argument objects in conjunction with \ref ArgumentObjectFactories :
1704  \code
1705  namespace vigra {
1706  // compute edgels from a gradient image
1707  template <class SrcIterator, class SrcAccessor,
1708  class BackInsertable, class GradValue>
1709  void
1710  cannyEdgelListThreshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1711  BackInsertable & edgels, GradValue grad_threshold);
1712 
1713  // compute edgels from a scalar image (determine gradient internally at 'scale')
1714  template <class SrcIterator, class SrcAccessor,
1715  class BackInsertable, class GradValue>
1716  void
1717  cannyEdgelListThreshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1718  BackInsertable & edgels, double scale, GradValue grad_threshold);
1719  }
1720  \endcode
1721  \deprecatedEnd
1722 
1723  <b> Usage:</b>
1724 
1725  <b>\#include</b> <vigra/edgedetection.hxx><br>
1726  Namespace: vigra
1727 
1728  \code
1729  MultiArray<2, unsigned char> src(w,h);
1730 
1731  // create empty edgel list
1732  std::vector<vigra::Edgel> edgels;
1733  ...
1734 
1735  // find edgels at scale 0.8, only considering gradient magnitudes above 2.0
1736  cannyEdgelListThreshold(src, edgels, 0.8, 2.0);
1737  \endcode
1738 
1739  \deprecatedUsage{cannyEdgelListThreshold}
1740  \code
1741  vigra::BImage src(w,h);
1742 
1743  // empty edgel list
1744  std::vector<vigra::Edgel> edgels;
1745  ...
1746 
1747  // find edgels at scale 0.8, only considering gradient above 2.0
1748  vigra::cannyEdgelListThreshold(srcImageRange(src), edgels, 0.8, 2.0);
1749  \endcode
1750  <b> Required Interface:</b>
1751  \code
1752  SrcImageIterator src_upperleft;
1753  SrcAccessor src_accessor;
1754 
1755  src_accessor(src_upperleft);
1756 
1757  BackInsertable edgels;
1758  edgels.push_back(Edgel());
1759  \endcode
1760  SrcAccessor::value_type must be a type convertible to float
1761  \deprecatedEnd
1762 
1763  <b> Preconditions:</b>
1764 
1765  \code
1766  scale > 0
1767  \endcode
1768 */
1769 doxygen_overloaded_function(template <...> void cannyEdgelListThreshold)
1770 
1771 template <class SrcIterator, class SrcAccessor,
1772  class BackInsertable, class GradValue>
1773 void
1774 cannyEdgelListThreshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1775  BackInsertable & edgels, double scale, GradValue grad_threshold)
1776 {
1777  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
1778  BasicImage<TinyVector<TmpType, 2> > grad(lr-ul);
1779  gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale);
1780 
1781  cannyEdgelListThreshold(srcImageRange(grad), edgels, grad_threshold);
1782 }
1783 
1784 template <class SrcIterator, class SrcAccessor,
1785  class BackInsertable, class GradValue>
1786 void
1787 cannyEdgelListThreshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
1788  BackInsertable & edgels, GradValue grad_threshold)
1789 {
1790  using namespace functor;
1791 
1792  typedef typename SrcAccessor::value_type SrcType;
1793  typedef typename NumericTraits<typename SrcType::value_type>::RealPromote TmpType;
1794  BasicImage<TmpType> magnitude(lr-ul);
1795  transformImage(srcIterRange(ul, lr, src), destImage(magnitude), norm(Arg1()));
1796 
1797  // find edgels
1798  internalCannyFindEdgels(ul, src, magnitude, edgels, grad_threshold);
1799 }
1800 
1801 template <class SrcIterator, class SrcAccessor,
1802  class BackInsertable, class GradValue>
1803 inline void
1804 cannyEdgelListThreshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1805  BackInsertable & edgels, double scale, GradValue grad_threshold)
1806 {
1807  cannyEdgelListThreshold(src.first, src.second, src.third, edgels, scale, grad_threshold);
1808 }
1809 
1810 template <class SrcIterator, class SrcAccessor,
1811  class BackInsertable, class GradValue>
1812 inline void
1813 cannyEdgelListThreshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1814  BackInsertable & edgels, GradValue grad_threshold)
1815 {
1816  cannyEdgelListThreshold(src.first, src.second, src.third, edgels, grad_threshold);
1817 }
1818 
1819 template <class T, class S,
1820  class BackInsertable, class GradValue>
1821 inline void
1822 cannyEdgelListThreshold(MultiArrayView<2, T, S> const & src,
1823  BackInsertable & edgels,
1824  double scale,
1825  GradValue grad_threshold)
1826 {
1827  cannyEdgelListThreshold(srcImageRange(src), edgels, scale, grad_threshold);
1828 }
1829 
1830 template <class T, class S,
1831  class BackInsertable, class GradValue>
1832 inline void
1833 cannyEdgelListThreshold(MultiArrayView<2, TinyVector<T, 2>, S> const & src,
1834  BackInsertable & edgels,
1835  GradValue grad_threshold)
1836 {
1837  cannyEdgelListThreshold(srcImageRange(src), edgels, grad_threshold);
1838 }
1839 
1840 
1841 /********************************************************/
1842 /* */
1843 /* cannyEdgeImage */
1844 /* */
1845 /********************************************************/
1846 
1847 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
1848 
1849  This operator first calls \ref cannyEdgelList() with the given scale to generate an
1850  edgel list for the given image. Then it scans this list and selects edgels
1851  whose strength is above the given <TT>gradient_threshold</TT>. For each of these
1852  edgels, the edgel's location is rounded to the nearest pixel, and that
1853  pixel marked with the given <TT>edge_marker</TT>.
1854 
1855  <b> Declarations:</b>
1856 
1857  pass 2D array views:
1858  \code
1859  namespace vigra {
1860  template <class T1, class S1,
1861  class T2, class S2,
1862  class GradValue, class DestValue>
1863  void
1864  cannyEdgeImage(MultiArrayView<2, T1, S1> const & src,
1865  MultiArrayView<2, T2, S2> dest,
1866  double scale,
1867  GradValue gradient_threshold,
1868  DestValue edge_marker);
1869  }
1870  \endcode
1871 
1872  \deprecatedAPI{cannyEdgeImage}
1873  pass \ref ImageIterators and \ref DataAccessors :
1874  \code
1875  namespace vigra {
1876  template <class SrcIterator, class SrcAccessor,
1877  class DestIterator, class DestAccessor,
1878  class GradValue, class DestValue>
1879  void cannyEdgeImage(
1880  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1881  DestIterator dul, DestAccessor da,
1882  double scale, GradValue gradient_threshold, DestValue edge_marker);
1883  }
1884  \endcode
1885  use argument objects in conjunction with \ref ArgumentObjectFactories :
1886  \code
1887  namespace vigra {
1888  template <class SrcIterator, class SrcAccessor,
1889  class DestIterator, class DestAccessor,
1890  class GradValue, class DestValue>
1891  void cannyEdgeImage(
1892  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1893  pair<DestIterator, DestAccessor> dest,
1894  double scale, GradValue gradient_threshold, DestValue edge_marker);
1895  }
1896  \endcode
1897  \deprecatedEnd
1898 
1899  <b> Usage:</b>
1900 
1901  <b>\#include</b> <vigra/edgedetection.hxx><br>
1902  Namespace: vigra
1903 
1904  \code
1905  MultiArray<2, unsigned char> src(w,h), edges(w,h);
1906  ...
1907 
1908  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
1909  cannyEdgeImage(src, edges, 0.8, 4.0, 1);
1910  \endcode
1911 
1912  \deprecatedUsage{cannyEdgeImage}
1913  \code
1914  vigra::BImage src(w,h), edges(w,h);
1915 
1916  // empty edge image
1917  edges = 0;
1918  ...
1919 
1920  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
1921  vigra::cannyEdgeImage(srcImageRange(src), destImage(edges),
1922  0.8, 4.0, 1);
1923  \endcode
1924  <b> Required Interface:</b>
1925  see also: \ref cannyEdgelList().
1926  \code
1927  DestImageIterator dest_upperleft;
1928  DestAccessor dest_accessor;
1929  DestValue edge_marker;
1930 
1931  dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
1932  \endcode
1933  \deprecatedEnd
1934 
1935  <b> Preconditions:</b>
1936 
1937  \code
1938  scale > 0
1939  gradient_threshold > 0
1940  \endcode
1941 */
1942 doxygen_overloaded_function(template <...> void cannyEdgeImage)
1943 
1944 template <class SrcIterator, class SrcAccessor,
1945  class DestIterator, class DestAccessor,
1946  class GradValue, class DestValue>
1947 void cannyEdgeImage(
1948  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
1949  DestIterator dul, DestAccessor da,
1950  double scale, GradValue gradient_threshold, DestValue edge_marker)
1951 {
1952  std::vector<Edgel> edgels;
1953 
1954  cannyEdgelListThreshold(sul, slr, sa, edgels, scale, gradient_threshold);
1955 
1956  int w = slr.x - sul.x;
1957  int h = slr.y - sul.y;
1958 
1959  for(unsigned int i=0; i<edgels.size(); ++i)
1960  {
1961  Diff2D pix((int)(edgels[i].x + 0.5), (int)(edgels[i].y + 0.5));
1962 
1963  if(pix.x < 0 || pix.x >= w || pix.y < 0 || pix.y >= h)
1964  continue;
1965 
1966  da.set(edge_marker, dul, pix);
1967  }
1968 }
1969 
1970 template <class SrcIterator, class SrcAccessor,
1971  class DestIterator, class DestAccessor,
1972  class GradValue, class DestValue>
1973 inline void
1974 cannyEdgeImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1975  pair<DestIterator, DestAccessor> dest,
1976  double scale, GradValue gradient_threshold, DestValue edge_marker)
1977 {
1978  cannyEdgeImage(src.first, src.second, src.third,
1979  dest.first, dest.second,
1980  scale, gradient_threshold, edge_marker);
1981 }
1982 
1983 template <class T1, class S1,
1984  class T2, class S2,
1985  class GradValue, class DestValue>
1986 inline void
1987 cannyEdgeImage(MultiArrayView<2, T1, S1> const & src,
1988  MultiArrayView<2, T2, S2> dest,
1989  double scale, GradValue gradient_threshold, DestValue edge_marker)
1990 {
1991  vigra_precondition(src.shape() == dest.shape(),
1992  "cannyEdgeImage(): shape mismatch between input and output.");
1993  cannyEdgeImage(srcImageRange(src),
1994  destImage(dest),
1995  scale, gradient_threshold, edge_marker);
1996 }
1997 
1998 /********************************************************/
1999 
2000 namespace detail {
2001 
2002 template <class DestIterator>
2003 int neighborhoodConfiguration(DestIterator dul)
2004 {
2005  int v = 0;
2006  NeighborhoodCirculator<DestIterator, EightNeighborCode> c(dul, EightNeighborCode::SouthEast);
2007  for(int i=0; i<8; ++i, --c)
2008  {
2009  v = (v << 1) | ((*c != 0) ? 1 : 0);
2010  }
2011 
2012  return v;
2013 }
2014 
2015 template <class GradValue>
2016 struct SimplePoint
2017 {
2018  Diff2D point;
2019  GradValue grad;
2020 
2021  SimplePoint(Diff2D const & p, GradValue g)
2022  : point(p), grad(g)
2023  {}
2024 
2025  bool operator<(SimplePoint const & o) const
2026  {
2027  return grad < o.grad;
2028  }
2029 
2030  bool operator>(SimplePoint const & o) const
2031  {
2032  return grad > o.grad;
2033  }
2034 };
2035 
2036 template <class SrcIterator, class SrcAccessor,
2037  class DestIterator, class DestAccessor,
2038  class GradValue, class DestValue>
2039 void cannyEdgeImageFromGrad(
2040  SrcIterator sul, SrcIterator slr, SrcAccessor grad,
2041  DestIterator dul, DestAccessor da,
2042  GradValue gradient_threshold, DestValue edge_marker)
2043 {
2044  typedef typename SrcAccessor::value_type PixelType;
2045  typedef typename NormTraits<PixelType>::SquaredNormType NormType;
2046 
2047  NormType zero = NumericTraits<NormType>::zero();
2048  double tan22_5 = M_SQRT2 - 1.0;
2049  typename NormTraits<GradValue>::SquaredNormType g2thresh = squaredNorm(gradient_threshold);
2050 
2051  int w = slr.x - sul.x;
2052  int h = slr.y - sul.y;
2053 
2054  sul += Diff2D(1,1);
2055  dul += Diff2D(1,1);
2056  Diff2D p(0,0);
2057 
2058  for(int y = 1; y < h-1; ++y, ++sul.y, ++dul.y)
2059  {
2060  SrcIterator sx = sul;
2061  DestIterator dx = dul;
2062  for(int x = 1; x < w-1; ++x, ++sx.x, ++dx.x)
2063  {
2064  PixelType g = grad(sx);
2065  NormType g2n = squaredNorm(g);
2066  if(g2n < g2thresh)
2067  continue;
2068 
2069  NormType g2n1, g2n3;
2070  // find out quadrant
2071  if(abs(g[1]) < tan22_5*abs(g[0]))
2072  {
2073  // north-south edge
2074  g2n1 = squaredNorm(grad(sx, Diff2D(-1, 0)));
2075  g2n3 = squaredNorm(grad(sx, Diff2D(1, 0)));
2076  }
2077  else if(abs(g[0]) < tan22_5*abs(g[1]))
2078  {
2079  // west-east edge
2080  g2n1 = squaredNorm(grad(sx, Diff2D(0, -1)));
2081  g2n3 = squaredNorm(grad(sx, Diff2D(0, 1)));
2082  }
2083  else if(g[0]*g[1] < zero)
2084  {
2085  // north-west-south-east edge
2086  g2n1 = squaredNorm(grad(sx, Diff2D(1, -1)));
2087  g2n3 = squaredNorm(grad(sx, Diff2D(-1, 1)));
2088  }
2089  else
2090  {
2091  // north-east-south-west edge
2092  g2n1 = squaredNorm(grad(sx, Diff2D(-1, -1)));
2093  g2n3 = squaredNorm(grad(sx, Diff2D(1, 1)));
2094  }
2095 
2096  if(g2n1 < g2n && g2n3 <= g2n)
2097  {
2098  da.set(edge_marker, dx);
2099  }
2100  }
2101  }
2102 }
2103 
2104 } // namespace detail
2105 
2106 /********************************************************/
2107 /* */
2108 /* cannyEdgeImageFromGradWithThinning */
2109 /* */
2110 /********************************************************/
2111 
2112 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
2113 
2114  The input pixels of this algorithm must be vectors of length 2 (see Required Interface below).
2115  It first searches for all pixels whose gradient magnitude is larger
2116  than the given <tt>gradient_threshold</tt> and larger than the magnitude of its two neighbors
2117  in gradient direction (where these neighbors are determined by nearest neighbor
2118  interpolation, i.e. according to the octant where the gradient points into).
2119  The resulting edge pixel candidates are then subjected to topological thinning
2120  so that the remaining edge pixels can be linked into edgel chains with a provable,
2121  non-heuristic algorithm. Thinning is performed so that the pixels with highest gradient
2122  magnitude survive. Optionally, the outermost pixels are marked as edge pixels
2123  as well when <tt>addBorder</tt> is true. The remaining pixels will be marked in the destination
2124  image with the value of <tt>edge_marker</tt> (all non-edge pixels remain untouched).
2125 
2126  <b> Declarations:</b>
2127 
2128  pass 2D array views:
2129  \code
2130  namespace vigra {
2131  template <class T1, class S1,
2132  class T2, class S2,
2133  class GradValue, class DestValue>
2134  void
2135  cannyEdgeImageFromGradWithThinning(MultiArrayView<2, TinyVector<T1, 2>, S1> const & src,
2136  MultiArrayView<2, T2, S2> dest,
2137  GradValue gradient_threshold,
2138  DestValue edge_marker,
2139  bool addBorder = true);
2140  }
2141  \endcode
2142 
2143  \deprecatedAPI{cannyEdgeImageFromGradWithThinning}
2144  pass \ref ImageIterators and \ref DataAccessors :
2145  \code
2146  namespace vigra {
2147  template <class SrcIterator, class SrcAccessor,
2148  class DestIterator, class DestAccessor,
2149  class GradValue, class DestValue>
2150  void cannyEdgeImageFromGradWithThinning(
2151  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
2152  DestIterator dul, DestAccessor da,
2153  GradValue gradient_threshold,
2154  DestValue edge_marker, bool addBorder = true);
2155  }
2156  \endcode
2157  use argument objects in conjunction with \ref ArgumentObjectFactories :
2158  \code
2159  namespace vigra {
2160  template <class SrcIterator, class SrcAccessor,
2161  class DestIterator, class DestAccessor,
2162  class GradValue, class DestValue>
2163  void cannyEdgeImageFromGradWithThinning(
2164  triple<SrcIterator, SrcIterator, SrcAccessor> src,
2165  pair<DestIterator, DestAccessor> dest,
2166  GradValue gradient_threshold,
2167  DestValue edge_marker, bool addBorder = true);
2168  }
2169  \endcode
2170  \deprecatedEnd
2171 
2172  <b> Usage:</b>
2173 
2174  <b>\#include</b> <vigra/edgedetection.hxx><br>
2175  Namespace: vigra
2176 
2177  \code
2178  MultiArray<2, unsigned char> src(w,h), edges(w,h);
2179 
2180  MultiArray<2, TinyVector<float, 2> > grad(w,h);
2181  // compute the image gradient at scale 0.8
2182  gaussianGradient(src, grad, 0.8);
2183 
2184  // find edges with gradient larger than 4.0, mark with 1, and add border
2185  cannyEdgeImageFromGradWithThinning(grad, edges, 4.0, 1, true);
2186  \endcode
2187 
2188  \deprecatedUsage{cannyEdgeImageFromGradWithThinning}
2189  \code
2190  vigra::BImage src(w,h), edges(w,h);
2191 
2192  vigra::FVector2Image grad(w,h);
2193  // compute the image gradient at scale 0.8
2194  vigra::gaussianGradient(srcImageRange(src), destImage(grad), 0.8);
2195 
2196  // empty edge image
2197  edges = 0;
2198  // find edges gradient larger than 4.0, mark with 1, and add border
2199  vigra::cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destImage(edges),
2200  4.0, 1, true);
2201  \endcode
2202  <b> Required Interface:</b>
2203  \code
2204  // the input pixel type must be a vector with two elements
2205  SrcImageIterator src_upperleft;
2206  SrcAccessor src_accessor;
2207  typedef SrcAccessor::value_type SrcPixel;
2208  typedef NormTraits<SrcPixel>::SquaredNormType SrcSquaredNormType;
2209 
2210  SrcPixel g = src_accessor(src_upperleft);
2211  SrcPixel::value_type g0 = g[0];
2212  SrcSquaredNormType gn = squaredNorm(g);
2213 
2214  DestImageIterator dest_upperleft;
2215  DestAccessor dest_accessor;
2216  DestValue edge_marker;
2217 
2218  dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
2219  \endcode
2220  \deprecatedEnd
2221 
2222  <b> Preconditions:</b>
2223 
2224  \code
2225  gradient_threshold > 0
2226  \endcode
2227 */
2228 doxygen_overloaded_function(template <...> void cannyEdgeImageFromGradWithThinning)
2229 
2230 template <class SrcIterator, class SrcAccessor,
2231  class DestIterator, class DestAccessor,
2232  class GradValue, class DestValue>
2234  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
2235  DestIterator dul, DestAccessor da,
2236  GradValue gradient_threshold,
2237  DestValue edge_marker, bool addBorder = true)
2238 {
2239  int w = slr.x - sul.x;
2240  int h = slr.y - sul.y;
2241 
2242  BImage edgeImage(w, h, BImage::value_type(0));
2243  BImage::traverser eul = edgeImage.upperLeft();
2244  BImage::Accessor ea = edgeImage.accessor();
2245  if(addBorder)
2246  initImageBorder(destImageRange(edgeImage), 1, 1);
2247  detail::cannyEdgeImageFromGrad(sul, slr, sa, eul, ea, gradient_threshold, 1);
2248 
2249  bool isSimplePoint[256] = {
2250  0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0,
2251  0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0,
2252  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1,
2253  1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1,
2254  0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1,
2255  0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0,
2256  0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0,
2257  1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1,
2258  0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
2259  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
2260  0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,
2261  0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0,
2262  1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0,
2263  0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1,
2264  1, 0, 1, 0 };
2265 
2266  eul += Diff2D(1,1);
2267  sul += Diff2D(1,1);
2268  int w2 = w-2;
2269  int h2 = h-2;
2270 
2271  typedef detail::SimplePoint<GradValue> SP;
2272  // use std::greater because we need the smallest gradients at the top of the queue
2273  std::priority_queue<SP, std::vector<SP>, std::greater<SP> > pqueue;
2274 
2275  Diff2D p(0,0);
2276  for(; p.y < h2; ++p.y)
2277  {
2278  for(p.x = 0; p.x < w2; ++p.x)
2279  {
2280  BImage::traverser e = eul + p;
2281  if(*e == 0)
2282  continue;
2283  int v = detail::neighborhoodConfiguration(e);
2284  if(isSimplePoint[v])
2285  {
2286  pqueue.push(SP(p, norm(sa(sul+p))));
2287  *e = 2; // remember that it is already in queue
2288  }
2289  }
2290  }
2291 
2292  const Diff2D dist[] = { Diff2D(-1,0), Diff2D(0,-1), Diff2D(1,0), Diff2D(0,1) };
2293 
2294  while(pqueue.size())
2295  {
2296  p = pqueue.top().point;
2297  pqueue.pop();
2298 
2299  BImage::traverser e = eul + p;
2300  int v = detail::neighborhoodConfiguration(e);
2301  if(!isSimplePoint[v])
2302  continue; // point may no longer be simple because its neighbors changed
2303 
2304  *e = 0; // delete simple point
2305 
2306  for(int i=0; i<4; ++i)
2307  {
2308  Diff2D pneu = p + dist[i];
2309  if(pneu.x == -1 || pneu.y == -1 || pneu.x == w2 || pneu.y == h2)
2310  continue; // do not remove points at the border
2311 
2312  BImage::traverser eneu = eul + pneu;
2313  if(*eneu == 1) // point is boundary and not yet in the queue
2314  {
2315  int v = detail::neighborhoodConfiguration(eneu);
2316  if(isSimplePoint[v])
2317  {
2318  pqueue.push(SP(pneu, norm(sa(sul+pneu))));
2319  *eneu = 2; // remember that it is already in queue
2320  }
2321  }
2322  }
2323  }
2324 
2325  initImageIf(destIterRange(dul, dul+Diff2D(w,h), da),
2326  maskImage(edgeImage), edge_marker);
2327 }
2328 
2329 template <class SrcIterator, class SrcAccessor,
2330  class DestIterator, class DestAccessor,
2331  class GradValue, class DestValue>
2333  triple<SrcIterator, SrcIterator, SrcAccessor> src,
2334  pair<DestIterator, DestAccessor> dest,
2335  GradValue gradient_threshold,
2336  DestValue edge_marker, bool addBorder = true)
2337 {
2338  cannyEdgeImageFromGradWithThinning(src.first, src.second, src.third,
2339  dest.first, dest.second,
2340  gradient_threshold, edge_marker, addBorder);
2341 }
2342 
2343 template <class T1, class S1,
2344  class T2, class S2,
2345  class GradValue, class DestValue>
2346 inline void
2347 cannyEdgeImageFromGradWithThinning(MultiArrayView<2, TinyVector<T1, 2>, S1> const & src,
2348  MultiArrayView<2, T2, S2> dest,
2349  GradValue gradient_threshold,
2350  DestValue edge_marker, bool addBorder = true)
2351 {
2352  vigra_precondition(src.shape() == dest.shape(),
2353  "cannyEdgeImageFromGradWithThinning(): shape mismatch between input and output.");
2354  cannyEdgeImageFromGradWithThinning(srcImageRange(src),
2355  destImage(dest),
2356  gradient_threshold, edge_marker, addBorder);
2357 }
2358 
2359 /********************************************************/
2360 /* */
2361 /* cannyEdgeImageWithThinning */
2362 /* */
2363 /********************************************************/
2364 
2365 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
2366 
2367  This operator first calls \ref gaussianGradient() to compute the gradient of the input
2368  image, ad then \ref cannyEdgeImageFromGradWithThinning() to generate an
2369  edge image. See there for more detailed documentation.
2370 
2371  <b> Declarations:</b>
2372 
2373  pass 2D array views:
2374  \code
2375  namespace vigra {
2376  template <class T1, class S1,
2377  class T2, class S2,
2378  class GradValue, class DestValue>
2379  void
2380  cannyEdgeImageWithThinning(MultiArrayView<2, T1, S1> const & src,
2381  MultiArrayView<2, T2, S2> dest,
2382  double scale,
2383  GradValue gradient_threshold,
2384  DestValue edge_marker,
2385  bool addBorder = true);
2386  }
2387  \endcode
2388 
2389  \deprecatedAPI{cannyEdgeImageWithThinning}
2390  pass \ref ImageIterators and \ref DataAccessors :
2391  \code
2392  namespace vigra {
2393  template <class SrcIterator, class SrcAccessor,
2394  class DestIterator, class DestAccessor,
2395  class GradValue, class DestValue>
2396  void cannyEdgeImageWithThinning(
2397  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
2398  DestIterator dul, DestAccessor da,
2399  double scale, GradValue gradient_threshold,
2400  DestValue edge_marker, bool addBorder = true);
2401  }
2402  \endcode
2403  use argument objects in conjunction with \ref ArgumentObjectFactories :
2404  \code
2405  namespace vigra {
2406  template <class SrcIterator, class SrcAccessor,
2407  class DestIterator, class DestAccessor,
2408  class GradValue, class DestValue>
2409  void cannyEdgeImageWithThinning(
2410  triple<SrcIterator, SrcIterator, SrcAccessor> src,
2411  pair<DestIterator, DestAccessor> dest,
2412  double scale, GradValue gradient_threshold,
2413  DestValue edge_marker, bool addBorder = true);
2414  }
2415  \endcode
2416  \deprecatedEnd
2417 
2418  <b> Usage:</b>
2419 
2420  <b>\#include</b> <vigra/edgedetection.hxx><br>
2421  Namespace: vigra
2422 
2423  \code
2424  MultiArray<2, unsigned char> src(w,h), edges(w,h);
2425  ...
2426 
2427  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1, and add border
2428  cannyEdgeImageWithThinning(src, edges, 0.8, 4.0, 1, true);
2429  \endcode
2430 
2431  \deprecatedUsage{cannyEdgeImageWithThinning}
2432  \code
2433  vigra::BImage src(w,h), edges(w,h);
2434 
2435  // empty edge image
2436  edges = 0;
2437  ...
2438 
2439  // find edges at scale 0.8 with gradient larger than 4.0, mark with 1, annd add border
2440  vigra::cannyEdgeImageWithThinning(srcImageRange(src), destImage(edges),
2441  0.8, 4.0, 1, true);
2442  \endcode
2443  <b> Required Interface:</b>
2444  see also: \ref cannyEdgelList().
2445  \code
2446  DestImageIterator dest_upperleft;
2447  DestAccessor dest_accessor;
2448  DestValue edge_marker;
2449 
2450  dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
2451  \endcode
2452  \deprecatedEnd
2453 
2454  <b> Preconditions:</b>
2455 
2456  \code
2457  scale > 0
2458  gradient_threshold > 0
2459  \endcode
2460 */
2461 doxygen_overloaded_function(template <...> void cannyEdgeImageWithThinning)
2462 
2463 template <class SrcIterator, class SrcAccessor,
2464  class DestIterator, class DestAccessor,
2465  class GradValue, class DestValue>
2467  SrcIterator sul, SrcIterator slr, SrcAccessor sa,
2468  DestIterator dul, DestAccessor da,
2469  double scale, GradValue gradient_threshold,
2470  DestValue edge_marker, bool addBorder = true)
2471 {
2472  // mark pixels that are higher than their neighbors in gradient direction
2473  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
2474  BasicImage<TinyVector<TmpType, 2> > grad(slr-sul);
2475  gaussianGradient(srcIterRange(sul, slr, sa), destImage(grad), scale);
2476  cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destIter(dul, da),
2477  gradient_threshold, edge_marker, addBorder);
2478 }
2479 
2480 template <class SrcIterator, class SrcAccessor,
2481  class DestIterator, class DestAccessor,
2482  class GradValue, class DestValue>
2483 inline void
2484 cannyEdgeImageWithThinning(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2485  pair<DestIterator, DestAccessor> dest,
2486  double scale, GradValue gradient_threshold,
2487  DestValue edge_marker, bool addBorder = true)
2488 {
2489  cannyEdgeImageWithThinning(src.first, src.second, src.third,
2490  dest.first, dest.second,
2491  scale, gradient_threshold, edge_marker, addBorder);
2492 }
2493 
2494 template <class T1, class S1,
2495  class T2, class S2,
2496  class GradValue, class DestValue>
2497 inline void
2498 cannyEdgeImageWithThinning(MultiArrayView<2, T1, S1> const & src,
2499  MultiArrayView<2, T2, S2> dest,
2500  double scale, GradValue gradient_threshold,
2501  DestValue edge_marker, bool addBorder = true)
2502 {
2503  vigra_precondition(src.shape() == dest.shape(),
2504  "cannyEdgeImageWithThinning(): shape mismatch between input and output.");
2505  cannyEdgeImageWithThinning(srcImageRange(src),
2506  destImage(dest),
2507  scale, gradient_threshold, edge_marker, addBorder);
2508 }
2509 
2510 /********************************************************/
2511 
2512 template <class SrcIterator, class SrcAccessor,
2513  class MaskImage, class BackInsertable, class GradValue>
2514 void internalCannyFindEdgels3x3(SrcIterator ul, SrcAccessor grad,
2515  MaskImage const & mask,
2516  BackInsertable & edgels,
2517  GradValue grad_thresh)
2518 {
2519  typedef typename SrcAccessor::value_type PixelType;
2520  typedef typename PixelType::value_type ValueType;
2521 
2522  vigra_precondition(grad_thresh >= NumericTraits<GradValue>::zero(),
2523  "cannyFindEdgels3x3(): gradient threshold must not be negative.");
2524 
2525  ul += Diff2D(1,1);
2526  for(int y=1; y<mask.height()-1; ++y, ++ul.y)
2527  {
2528  SrcIterator ix = ul;
2529  for(int x=1; x<mask.width()-1; ++x, ++ix.x)
2530  {
2531  if(!mask(x,y))
2532  continue;
2533 
2534  ValueType gradx = grad.getComponent(ix, 0);
2535  ValueType grady = grad.getComponent(ix, 1);
2536  double mag = hypot(gradx, grady);
2537  if(mag <= grad_thresh)
2538  continue;
2539  double c = gradx / mag,
2540  s = grady / mag;
2541 
2542  Matrix<double> ml(3,3), mr(3,1), l(3,1), r(3,1);
2543  l(0,0) = 1.0;
2544 
2545  for(int yy = -1; yy <= 1; ++yy)
2546  {
2547  for(int xx = -1; xx <= 1; ++xx)
2548  {
2549  double u = c*xx + s*yy;
2550  double v = norm(grad(ix, Diff2D(xx, yy)));
2551  l(1,0) = u;
2552  l(2,0) = u*u;
2553  ml += outer(l);
2554  mr += v*l;
2555  }
2556  }
2557 
2558  linearSolve(ml, mr, r);
2559 
2560  Edgel edgel;
2561 
2562  // local maximum => quadratic interpolation of sub-pixel location
2563  double del = -r(1,0) / 2.0 / r(2,0);
2564  if(std::fabs(del) > 1.5) // don't move by more than about a pixel diameter
2565  del = 0.0;
2566  edgel.x = Edgel::value_type(x + c*del);
2567  edgel.y = Edgel::value_type(y + s*del);
2568  edgel.strength = Edgel::value_type(mag);
2569  double orientation = VIGRA_CSTD::atan2(grady, gradx) + 0.5*M_PI;
2570  if(orientation < 0.0)
2571  orientation += 2.0*M_PI;
2572  edgel.orientation = Edgel::value_type(orientation);
2573  edgels.push_back(edgel);
2574  }
2575  }
2576 }
2577 
2578 /********************************************************/
2579 /* */
2580 /* cannyEdgelList3x3 */
2581 /* */
2582 /********************************************************/
2583 
2584 /** \brief Improved implementation of Canny's edge detector.
2585 
2586  This operator first computes pixels which are crossed by the edge using
2587  cannyEdgeImageWithThinning(). The gradient magnitudes in the 3x3 neighborhood of these
2588  pixels are then projected onto the normal of the edge (as determined
2589  by the gradient direction). The edgel's subpixel location is found by fitting a
2590  parabola through the 9 gradient values and determining the parabola's tip.
2591  A new \ref Edgel is appended to the given vector of <TT>edgels</TT>. Since the parabola
2592  is fitted to 9 points rather than 3 points as in cannyEdgelList(), the accuracy is higher.
2593 
2594  The function can be called in two modes: If you pass a 'scale', it is assumed that the
2595  original image is scalar, and the Gaussian gradient is internally computed at the
2596  given 'scale'. If the function is called without scale parameter, it is assumed that
2597  the given image already contains the gradient (i.e. its value_type must be
2598  a vector of length 2).
2599 
2600  <b> Declarations:</b>
2601 
2602  pass 2D array views:
2603  \code
2604  namespace vigra {
2605  // compute edgels from a scalar image (determine gradient internally at 'scale')
2606  template <class T, class S, class BackInsertable>
2607  void
2608  cannyEdgelList3x3(MultiArrayView<2, T, S> const & src,
2609  BackInsertable & edgels,
2610  double scale);
2611 
2612  // compute edgels from a pre-computed gradient image
2613  template <class T, class S, class BackInsertable>
2614  void
2615  cannyEdgelList3x3(MultiArrayView<2, TinyVector<T, 2>, S> const & src,
2616  BackInsertable & edgels);
2617  }
2618  \endcode
2619 
2620  \deprecatedAPI{cannyEdgelList3x3}
2621  pass \ref ImageIterators and \ref DataAccessors :
2622  \code
2623  namespace vigra {
2624  // compute edgels from a gradient image
2625  template <class SrcIterator, class SrcAccessor, class BackInsertable>
2626  void cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2627  BackInsertable & edgels);
2628 
2629  // compute edgels from a scalar image (determine gradient internally at 'scale')
2630  template <class SrcIterator, class SrcAccessor, class BackInsertable>
2631  void cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2632  BackInsertable & edgels, double scale);
2633  }
2634  \endcode
2635  use argument objects in conjunction with \ref ArgumentObjectFactories :
2636  \code
2637  namespace vigra {
2638  // compute edgels from a gradient image
2639  template <class SrcIterator, class SrcAccessor, class BackInsertable>
2640  void
2641  cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2642  BackInsertable & edgels);
2643 
2644  // compute edgels from a scalar image (determine gradient internally at 'scale')
2645  template <class SrcIterator, class SrcAccessor, class BackInsertable>
2646  void
2647  cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2648  BackInsertable & edgels, double scale);
2649  }
2650  \endcode
2651  \deprecatedEnd
2652 
2653  <b> Usage:</b>
2654 
2655  <b>\#include</b> <vigra/edgedetection.hxx><br>
2656  Namespace: vigra
2657 
2658  \code
2659  MultiArray<2, unsigned char> src(w,h);
2660 
2661  // create empty edgel list
2662  std::vector<vigra::Edgel> edgels;
2663  ...
2664 
2665  // find edgels at scale 0.8
2666  cannyEdgelList3x3(src, edgels, 0.8);
2667  \endcode
2668 
2669  \deprecatedUsage{cannyEdgelList3x3}
2670  \code
2671  vigra::BImage src(w,h);
2672 
2673  // empty edgel list
2674  std::vector<vigra::Edgel> edgels;
2675  ...
2676 
2677  // find edgels at scale 0.8
2678  vigra::cannyEdgelList3x3(srcImageRange(src), edgels, 0.8);
2679  \endcode
2680  <b> Required Interface:</b>
2681  \code
2682  SrcImageIterator src_upperleft;
2683  SrcAccessor src_accessor;
2684 
2685  src_accessor(src_upperleft);
2686 
2687  BackInsertable edgels;
2688  edgels.push_back(Edgel());
2689  \endcode
2690  SrcAccessor::value_type must be a type convertible to float
2691  \deprecatedEnd
2692 
2693  <b> Preconditions:</b>
2694 
2695  \code
2696  scale > 0
2697  \endcode
2698 */
2699 doxygen_overloaded_function(template <...> void cannyEdgelList3x3)
2700 
2701 template <class SrcIterator, class SrcAccessor, class BackInsertable>
2702 void
2703 cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2704  BackInsertable & edgels, double scale)
2705 {
2706  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
2707  BasicImage<TinyVector<TmpType, 2> > grad(lr-ul);
2708  gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale);
2709 
2710  cannyEdgelList3x3(srcImageRange(grad), edgels);
2711 }
2712 
2713 template <class SrcIterator, class SrcAccessor, class BackInsertable>
2714 void
2715 cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2716  BackInsertable & edgels)
2717 {
2718  typedef typename NormTraits<typename SrcAccessor::value_type>::NormType NormType;
2719 
2720  UInt8Image edges(lr-ul);
2721  cannyEdgeImageFromGradWithThinning(srcIterRange(ul, lr, src), destImage(edges),
2722  0.0, 1, false);
2723 
2724  // find edgels
2725  internalCannyFindEdgels3x3(ul, src, edges, edgels, NumericTraits<NormType>::zero());
2726 }
2727 
2728 template <class SrcIterator, class SrcAccessor, class BackInsertable>
2729 inline void
2730 cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2731  BackInsertable & edgels, double scale)
2732 {
2733  cannyEdgelList3x3(src.first, src.second, src.third, edgels, scale);
2734 }
2735 
2736 template <class SrcIterator, class SrcAccessor, class BackInsertable>
2737 inline void
2738 cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2739  BackInsertable & edgels)
2740 {
2741  cannyEdgelList3x3(src.first, src.second, src.third, edgels);
2742 }
2743 
2744 template <class T, class S, class BackInsertable>
2745 inline void
2746 cannyEdgelList3x3(MultiArrayView<2, T, S> const & src,
2747  BackInsertable & edgels, double scale)
2748 {
2749  cannyEdgelList3x3(srcImageRange(src), edgels, scale);
2750 }
2751 
2752 template <class T, class S, class BackInsertable>
2753 inline void
2754 cannyEdgelList3x3(MultiArrayView<2, TinyVector<T, 2>, S> const & src,
2755  BackInsertable & edgels)
2756 {
2757  cannyEdgelList3x3(srcImageRange(src), edgels);
2758 }
2759 
2760 /********************************************************/
2761 /* */
2762 /* cannyEdgelList3x3Threshold */
2763 /* */
2764 /********************************************************/
2765 
2766 /** \brief Improved implementation of Canny's edge detector with thresholding.
2767 
2768  This function works exactly like \ref cannyEdgelList3x3(), but
2769  you also pass a threshold for the minimal gradient magnitude,
2770  so that edgels whose strength is below the threshold are not
2771  inserted into the edgel list.
2772 
2773  <b> Declarations:</b>
2774 
2775  pass 2D array views:
2776  \code
2777  namespace vigra {
2778  // compute edgels from a gradient image
2779  template <class SrcIterator, class SrcAccessor,
2780  class BackInsertable, class GradValue>
2781  void
2782  cannyEdgelList3x3Threshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2783  BackInsertable & edgels, GradValue grad_thresh);
2784 
2785  // compute edgels from a scalar image (determine gradient internally at 'scale')
2786  template <class SrcIterator, class SrcAccessor,
2787  class BackInsertable, class GradValue>
2788  void
2789  cannyEdgelList3x3Threshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2790  BackInsertable & edgels, double scale, GradValue grad_thresh);
2791  }
2792  \endcode
2793 
2794  \deprecatedAPI{cannyEdgelList3x3Threshold}
2795  pass \ref ImageIterators and \ref DataAccessors :
2796  \code
2797  namespace vigra {
2798  // compute edgels from a gradient image
2799  template <class SrcIterator, class SrcAccessor,
2800  class BackInsertable, class GradValue>
2801  void
2802  cannyEdgelList3x3Threshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2803  BackInsertable & edgels, GradValue grad_thresh);
2804 
2805  // compute edgels from a scalar image (determine gradient internally at 'scale')
2806  template <class SrcIterator, class SrcAccessor,
2807  class BackInsertable, class GradValue>
2808  void
2809  cannyEdgelList3x3Threshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2810  BackInsertable & edgels, double scale, GradValue grad_thresh);
2811  }
2812  \endcode
2813  use argument objects in conjunction with \ref ArgumentObjectFactories :
2814  \code
2815  namespace vigra {
2816  // compute edgels from a gradient image
2817  template <class SrcIterator, class SrcAccessor,
2818  class BackInsertable, class GradValue>
2819  void
2820  cannyEdgelList3x3Threshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2821  BackInsertable & edgels, GradValue grad_thresh);
2822 
2823  // compute edgels from a scalar image (determine gradient internally at 'scale')
2824  template <class SrcIterator, class SrcAccessor,
2825  class BackInsertable, class GradValue>
2826  void
2827  cannyEdgelList3x3Threshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2828  BackInsertable & edgels, double scale, GradValue grad_thresh);
2829  }
2830  \endcode
2831  \deprecatedEnd
2832 
2833  <b> Usage:</b>
2834 
2835  <b>\#include</b> <vigra/edgedetection.hxx><br>
2836  Namespace: vigra
2837 
2838  \code
2839  MultiArray<2, unsigned char> src(w,h);
2840 
2841  // create empty edgel list
2842  std::vector<vigra::Edgel> edgels;
2843  ...
2844 
2845  // find edgels at scale 0.8 whose gradient is at least 4.0
2846  cannyEdgelList3x3Threshold(src, edgels, 0.8, 4.0);
2847  \endcode
2848 
2849  \deprecatedUsage{cannyEdgelList3x3Threshold}
2850  \code
2851  vigra::BImage src(w,h);
2852 
2853  // empty edgel list
2854  std::vector<vigra::Edgel> edgels;
2855  ...
2856 
2857  // find edgels at scale 0.8 whose gradient is at least 4.0
2858  vigra::cannyEdgelList3x3Threshold(srcImageRange(src), edgels, 0.8, 4.0);
2859  \endcode
2860  <b> Required Interface:</b>
2861  \code
2862  SrcImageIterator src_upperleft;
2863  SrcAccessor src_accessor;
2864 
2865  src_accessor(src_upperleft);
2866 
2867  BackInsertable edgels;
2868  edgels.push_back(Edgel());
2869  \endcode
2870  SrcAccessor::value_type must be a type convertible to float
2871  \deprecatedEnd
2872 
2873  <b> Preconditions:</b>
2874 
2875  \code
2876  scale > 0
2877  \endcode
2878 */
2879 doxygen_overloaded_function(template <...> void cannyEdgelList3x3Threshold)
2880 
2881 template <class SrcIterator, class SrcAccessor,
2882  class BackInsertable, class GradValue>
2883 void
2884 cannyEdgelList3x3Threshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2885  BackInsertable & edgels, double scale, GradValue grad_thresh)
2886 {
2887  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
2888  BasicImage<TinyVector<TmpType, 2> > grad(lr-ul);
2889  gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale);
2890 
2891  cannyEdgelList3x3Threshold(srcImageRange(grad), edgels, grad_thresh);
2892 }
2893 
2894 template <class SrcIterator, class SrcAccessor,
2895  class BackInsertable, class GradValue>
2896 void
2897 cannyEdgelList3x3Threshold(SrcIterator ul, SrcIterator lr, SrcAccessor src,
2898  BackInsertable & edgels, GradValue grad_thresh)
2899 {
2900  UInt8Image edges(lr-ul);
2901  cannyEdgeImageFromGradWithThinning(srcIterRange(ul, lr, src), destImage(edges),
2902  0.0, 1, false);
2903 
2904  // find edgels
2905  internalCannyFindEdgels3x3(ul, src, edges, edgels, grad_thresh);
2906 }
2907 
2908 template <class SrcIterator, class SrcAccessor,
2909  class BackInsertable, class GradValue>
2910 inline void
2911 cannyEdgelList3x3Threshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2912  BackInsertable & edgels, double scale, GradValue grad_thresh)
2913 {
2914  cannyEdgelList3x3Threshold(src.first, src.second, src.third, edgels, scale, grad_thresh);
2915 }
2916 
2917 template <class SrcIterator, class SrcAccessor,
2918  class BackInsertable, class GradValue>
2919 inline void
2920 cannyEdgelList3x3Threshold(triple<SrcIterator, SrcIterator, SrcAccessor> src,
2921  BackInsertable & edgels, GradValue grad_thresh)
2922 {
2923  cannyEdgelList3x3Threshold(src.first, src.second, src.third, edgels, grad_thresh);
2924 }
2925 
2926 template <class T, class S,
2927  class BackInsertable, class GradValue>
2928 inline void
2929 cannyEdgelList3x3Threshold(MultiArrayView<2, T, S> const & src,
2930  BackInsertable & edgels, double scale, GradValue grad_thresh)
2931 {
2932  cannyEdgelList3x3Threshold(srcImageRange(src), edgels, scale, grad_thresh);
2933 }
2934 
2935 template <class T, class S,
2936  class BackInsertable, class GradValue>
2937 inline void
2938 cannyEdgelList3x3Threshold(MultiArrayView<2, TinyVector<T, 2>, S> const & src,
2939  BackInsertable & edgels,
2940  GradValue grad_thresh)
2941 {
2942  cannyEdgelList3x3Threshold(srcImageRange(src), edgels, grad_thresh);
2943 }
2944 
2945 //@}
2946 
2947 /** \page CrackEdgeImage Crack Edge Image
2948 
2949 Crack edges are marked <i>between</i> the pixels of an image.
2950 A Crack Edge Image is an image that represents these edges. In order
2951 to accommodate the cracks, the Crack Edge Image must be twice as large
2952 as the original image (precisely (2*w - 1) by (2*h - 1)). A Crack Edge Image
2953 can easily be derived from a binary image or from the signs of the
2954 response of a Laplacian filter. Consider the following sketch, where
2955 <TT>+</TT> encodes the foreground, <TT>-</TT> the background, and
2956 <TT>*</TT> the resulting crack edges.
2957 
2958  \code
2959 sign of difference image insert cracks resulting CrackEdgeImage
2960 
2961  + . - . - . * . . .
2962  + - - . . . . . . * * * .
2963  + + - => + . + . - => . . . * .
2964  + + + . . . . . . . . * *
2965  + . + . + . . . . .
2966  \endcode
2967 
2968 Starting from the original binary image (left), we insert crack pixels
2969 to get to the double-sized image (center). Finally, we mark all
2970 crack pixels whose non-crack neighbors have different signs as
2971 crack edge points, while all other pixels (crack and non-crack) become
2972 region pixels.
2973 
2974 <b>Requirements on a Crack Edge Image:</b>
2975 
2976 <ul>
2977  <li>Crack Edge Images have odd width and height.
2978  <li>Crack pixels have at least one odd coordinate.
2979  <li>Only crack pixels may be marked as edge points.
2980  <li>Crack pixels with two odd coordinates must be marked as edge points
2981  whenever any of their neighboring crack pixels was marked.
2982 </ul>
2983 
2984 The last two requirements ensure that both edges and regions are 4-connected.
2985 Thus, 4-connectivity and 8-connectivity yield identical connected
2986 components in a Crack Edge Image (so called <i>well-composedness</i>).
2987 This ensures that Crack Edge Images have nice topological properties
2988 (cf. L. J. Latecki: "Well-Composed Sets", Academic Press, 2000).
2989 */
2990 
2991 
2992 } // namespace vigra
2993 
2994 #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.10.0 (Mon Nov 18 2013)