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