1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2002 by Ullrich Koethe */
4 /* Cognitive Systems Group, University of Hamburg, Germany */
5 /* */
6 /* This file is part of the VIGRA computer vision library. */
7 /* ( Version 1.5.0, Dec 07 2006 ) */
8 /* The VIGRA Website is */
9 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */
10 /* Please direct questions, bug reports, and contributions to */
11 /* koethe@informatik.uni-hamburg.de or */
12 /* vigra@kogs1.informatik.uni-hamburg.de */
13 /* */
14 /* Permission is hereby granted, free of charge, to any person */
15 /* obtaining a copy of this software and associated documentation */
16 /* files (the "Software"), to deal in the Software without */
17 /* restriction, including without limitation the rights to use, */
18 /* copy, modify, merge, publish, distribute, sublicense, and/or */
19 /* sell copies of the Software, and to permit persons to whom the */
20 /* Software is furnished to do so, subject to the following */
21 /* conditions: */
22 /* */
23 /* The above copyright notice and this permission notice shall be */
24 /* included in all copies or substantial portions of the */
25 /* Software. */
26 /* */
27 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
28 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
29 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
30 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
31 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
32 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
33 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
34 /* OTHER DEALINGS IN THE SOFTWARE. */
35 /* */
36 /************************************************************************/
37
38
39 #ifndef VIGRA_FLATMORPHOLOGY_HXX
40 #define VIGRA_FLATMORPHOLOGY_HXX
41
42 #include <cmath>
43 #include <vector>
44 #include "utilities.hxx"
45
46 namespace vigra {
47
48 /** \addtogroup Morphology Basic Morphological Operations
49 Perform erosion, dilation, and median with disc structuring functions
50 */
51 //@{
52
53 /********************************************************/
54 /* */
55 /* discRankOrderFilter */
56 /* */
57 /********************************************************/
58
59 /** \brief Apply rank order filter with disc structuring function to the image.
60
61 The pixel values of the source image <b> must</b> be in the range
62 0...255. Radius must be >= 0. Rank must be in the range 0.0 <= rank
63 <= 1.0. The filter acts as a minimum filter if rank = 0.0,
64 as a median if rank = 0.5, and as a maximum filter if rank = 1.0.
65 Accessor are used to access the pixel data.
66
67 <b> Declarations:</b>
68
69 pass arguments explicitely:
70 \code
71 namespace vigra {
72 template <class SrcIterator, class SrcAccessor,
73 class DestIterator, class DestAccessor>
74 void
75 discRankOrderFilter(SrcIterator upperleft1,
76 SrcIterator lowerright1, SrcAccessor sa,
77 DestIterator upperleft2, DestAccessor da,
78 int radius, float rank)
79 }
80 \endcode
81
82
83 use argument objects in conjunction with \ref ArgumentObjectFactories:
84 \code
85 namespace vigra {
86 template <class SrcIterator, class SrcAccessor,
87 class DestIterator, class DestAccessor>
88 void
89 discRankOrderFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
90 pair<DestIterator, DestAccessor> dest,
91 int radius, float rank)
92 }
93 \endcode
94
95 <b> Usage:</b>
96
97 <b>\#include</b> "<a href="flatmorphology_8hxx-source.html">vigra/flatmorphology.hxx</a>"<br>
98 Namespace: vigra
99
100 \code
101 vigra::CImage src, dest;
102
103 // do median filtering
104 vigra::discRankOrderFilter(srcImageRange(src), destImage(dest), 10, 0.5);
105 \endcode
106
107 <b> Required Interface:</b>
108
109 \code
110 SrcIterator src_upperleft;
111 DestIterator dest_upperleft;
112 int x, y;
113 unsigned char value;
114
115 SrcAccessor src_accessor;
116 DestAccessor dest_accessor;
117
118 // value_type of accessor must be convertible to unsigned char
119 value = src_accessor(src_upperleft, x, y);
120
121 dest_accessor.set(value, dest_upperleft, x, y);
122 \endcode
123
124 <b> Preconditions:</b>
125
126 \code
127 for all source pixels: 0 <= value <= 255
128
129 (rank >= 0.0) && (rank <= 1.0)
130 radius >= 0
131 \endcode
132
133 */
134 template <class SrcIterator, class SrcAccessor,
135 class DestIterator, class DestAccessor>
136 void
discRankOrderFilter(SrcIterator upperleft1,SrcIterator lowerright1,SrcAccessor sa,DestIterator upperleft2,DestAccessor da,int radius,float rank)137 discRankOrderFilter(SrcIterator upperleft1,
138 SrcIterator lowerright1, SrcAccessor sa,
139 DestIterator upperleft2, DestAccessor da,
140 int radius, float rank)
141 {
142 vigra_precondition((rank >= 0.0) && (rank <= 1.0),
143 "discRankOrderFilter(): Rank must be between 0 and 1"
144 " (inclusive).");
145
146 vigra_precondition(radius >= 0,
147 "discRankOrderFilter(): Radius must be >= 0.");
148
149 int i, x, y, xmax, ymax, xx, yy;
150 int rankpos, winsize, leftsum;
151
152 long hist[256];
153
154 // prepare structuring function
155 std::vector<int> struct_function(radius+1);
156 struct_function[0] = radius;
157
158 double r2 = (double)radius*radius;
159 for(i=1; i<=radius; ++i)
160 {
161 double r = (double) i - 0.5;
162 struct_function[i] = (int)(VIGRA_CSTD::sqrt(r2 - r*r) + 0.5);
163 }
164
165 int w = lowerright1.x - upperleft1.x;
166 int h = lowerright1.y - upperleft1.y;
167
168 SrcIterator ys(upperleft1);
169 DestIterator yd(upperleft2);
170
171 for(y=0; y<h; ++y, ++ys.y, ++yd.y)
172 {
173 SrcIterator xs(ys);
174 DestIterator xd(yd);
175
176 // first column
177 int x0 = 0;
178 int y0 = y;
179 int x1 = w - 1;
180 int y1 = h - y - 1;
181
182 // clear histogram
183 for(i=0; i<256; ++i) hist[i] = 0;
184 winsize = 0;
185
186 // init histogram
187 ymax = (y1 < radius) ? y1 : radius;
188 for(yy=0; yy<=ymax; ++yy)
189 {
190 xmax = (x1 < struct_function[yy]) ? x1 : struct_function[yy];
191 for(xx=0; xx<=xmax; ++xx)
192 {
193 hist[sa(xs, Diff2D(xx, yy))]++;
194 winsize++;
195 }
196 }
197
198 ymax = (y0 < radius) ? y0 : radius;
199 for(yy=1; yy<=ymax; ++yy)
200 {
201 xmax = (x1 < struct_function[yy]) ? x1 : struct_function[yy];
202 for(xx=0; xx<=xmax; ++xx)
203 {
204 hist[sa(xs, Diff2D(xx, -yy))]++;
205 winsize++;
206 }
207 }
208
209 // find the desired histogramm bin
210 leftsum = 0;
211 if(rank == 0.0)
212 {
213 for(i=0; i<256; i++)
214 {
215 if(hist[i]) break;
216 }
217 rankpos = i;
218 }
219 else
220 {
221 for(i=0; i<256; i++)
222 {
223 if((float)(hist[i]+leftsum) / winsize >= rank) break;
224 leftsum += hist[i];
225 }
226 rankpos = i;
227 }
228
229 da.set(rankpos, xd);
230
231 ++xs.x;
232 ++xd.x;
233
234 // inner columns
235 for(x=1; x<w; ++x, ++xs.x, ++xd.x)
236 {
237 x0 = x;
238 y0 = y;
239 x1 = w - x - 1;
240 y1 = h - y - 1;
241
242 // update histogramm
243 // remove pixels at left border
244 yy = (y1 < radius) ? y1 : radius;
245 for(; yy>=0; yy--)
246 {
247 unsigned char cur;
248 xx = struct_function[yy]+1;
249 if(xx > x0) break;
250
251 cur = sa(xs, Diff2D(-xx, yy));
252
253 hist[cur]--;
254 if(cur < rankpos) leftsum--;
255 winsize--;
256 }
257 yy = (y0 < radius) ? y0 : radius;
258 for(; yy>=1; yy--)
259 {
260 unsigned char cur;
261 xx = struct_function[yy]+1;
262 if(xx > x0) break;
263
264 cur = sa(xs, Diff2D(-xx, -yy));
265
266 hist[cur]--;
267 if(cur < rankpos) leftsum--;
268 winsize--;
269 }
270
271 // add pixels at right border
272 yy = (y1 < radius) ? y1 : radius;
273 for(; yy>=0; yy--)
274 {
275 unsigned char cur;
276 xx = struct_function[yy];
277 if(xx > x1) break;
278
279 cur = sa(xs, Diff2D(xx, yy));
280
281 hist[cur]++;
282 if(cur < rankpos) leftsum++;
283 winsize++;
284 }
285 yy = (y0 < radius) ? y0 : radius;
286 for(; yy>=1; yy--)
287 {
288 unsigned char cur;
289 xx = struct_function[yy];
290 if(xx > x1) break;
291
292 cur = sa(xs, Diff2D(xx, -yy));
293
294 hist[cur]++;
295 if(cur < rankpos) leftsum++;
296 winsize++;
297 }
298
299 // find the desired histogramm bin
300 if(rank == 0.0)
301 {
302 if(leftsum == 0)
303 {
304 // search to the right
305 for(i=rankpos; i<256; i++)
306 {
307 if(hist[i]) break;
308 }
309 rankpos = i;
310 }
311 else
312 {
313 // search to the left
314 for(i=rankpos-1; i>=0; i--)
315 {
316 leftsum -= hist[i];
317 if(leftsum == 0) break;
318 }
319 rankpos = i;
320 }
321 }
322 else // rank > 0.0
323 {
324 if((float)leftsum / winsize < rank)
325 {
326 // search to the right
327 for(i=rankpos; i<256; i++)
328 {
329 if((float)(hist[i]+leftsum) / winsize >= rank) break;
330 leftsum+=hist[i];
331 }
332 rankpos = i;
333 }
334 else
335 {
336 /// search to the left
337 for(i=rankpos-1; i>=0; i--)
338 {
339 leftsum-=hist[i];
340 if((float)leftsum / winsize < rank) break;
341 }
342 rankpos = i;
343 }
344 }
345
346 da.set(rankpos, xd);
347 }
348 }
349 }
350
351 template <class SrcIterator, class SrcAccessor,
352 class DestIterator, class DestAccessor>
353 void
discRankOrderFilter(triple<SrcIterator,SrcIterator,SrcAccessor> src,pair<DestIterator,DestAccessor> dest,int radius,float rank)354 discRankOrderFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
355 pair<DestIterator, DestAccessor> dest,
356 int radius, float rank)
357 {
358 discRankOrderFilter(src.first, src.second, src.third,
359 dest.first, dest.second,
360 radius, rank);
361 }
362
363 /********************************************************/
364 /* */
365 /* discErosion */
366 /* */
367 /********************************************************/
368
369 /** \brief Apply erosion (minimum) filter with disc of given radius to image.
370
371 This is an abbreviation vor the rank order filter with rank = 0.0.
372 See \ref discRankOrderFilter() for more information.
373
374 <b> Declarations:</b>
375
376 pass arguments explicitely:
377 \code
378 namespace vigra {
379 template <class SrcIterator, class SrcAccessor,
380 class DestIterator, class DestAccessor>
381 inline void
382 discErosion(SrcIterator upperleft1,
383 SrcIterator lowerright1, SrcAccessor sa,
384 DestIterator upperleft2, DestAccessor da,
385 int radius)
386 }
387 \endcode
388
389
390 use argument objects in conjunction with \ref ArgumentObjectFactories:
391 \code
392 namespace vigra {
393 template <class SrcIterator, class SrcAccessor,
394 class DestIterator, class DestAccessor>
395 void
396 discErosion(triple<SrcIterator, SrcIterator, SrcAccessor> src,
397 pair<DestIterator, DestAccessor> dest,
398 int radius)
399 }
400 \endcode
401
402 */
403 template <class SrcIterator, class SrcAccessor,
404 class DestIterator, class DestAccessor>
405 inline void
discErosion(SrcIterator upperleft1,SrcIterator lowerright1,SrcAccessor sa,DestIterator upperleft2,DestAccessor da,int radius)406 discErosion(SrcIterator upperleft1,
407 SrcIterator lowerright1, SrcAccessor sa,
408 DestIterator upperleft2, DestAccessor da,
409 int radius)
410 {
411 vigra_precondition(radius >= 0, "discErosion(): Radius must be >= 0.");
412
413 discRankOrderFilter(upperleft1, lowerright1, sa,
414 upperleft2, da, radius, 0.0);
415 }
416
417 template <class SrcIterator, class SrcAccessor,
418 class DestIterator, class DestAccessor>
419 void
discErosion(triple<SrcIterator,SrcIterator,SrcAccessor> src,pair<DestIterator,DestAccessor> dest,int radius)420 discErosion(triple<SrcIterator, SrcIterator, SrcAccessor> src,
421 pair<DestIterator, DestAccessor> dest,
422 int radius)
423 {
424 vigra_precondition(radius >= 0, "discErosion(): Radius must be >= 0.");
425
426 discRankOrderFilter(src.first, src.second, src.third,
427 dest.first, dest.second,
428 radius, 0.0);
429 }
430
431 /********************************************************/
432 /* */
433 /* discDilation */
434 /* */
435 /********************************************************/
436
437 /** \brief Apply dilation (maximum) filter with disc of given radius to image.
438
439 This is an abbreviation vor the rank order filter with rank = 1.0.
440 See \ref discRankOrderFilter() for more information.
441
442 <b> Declarations:</b>
443
444 pass arguments explicitely:
445 \code
446 namespace vigra {
447 template <class SrcIterator, class SrcAccessor,
448 class DestIterator, class DestAccessor>
449 inline void
450 discDilation(SrcIterator upperleft1,
451 SrcIterator lowerright1, SrcAccessor sa,
452 DestIterator upperleft2, DestAccessor da,
453 int radius)
454 }
455 \endcode
456
457
458 use argument objects in conjunction with \ref ArgumentObjectFactories:
459 \code
460 namespace vigra {
461 template <class SrcIterator, class SrcAccessor,
462 class DestIterator, class DestAccessor>
463 void
464 discDilation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
465 pair<DestIterator, DestAccessor> dest,
466 int radius)
467 }
468 \endcode
469
470 */
471 template <class SrcIterator, class SrcAccessor,
472 class DestIterator, class DestAccessor>
473 inline void
discDilation(SrcIterator upperleft1,SrcIterator lowerright1,SrcAccessor sa,DestIterator upperleft2,DestAccessor da,int radius)474 discDilation(SrcIterator upperleft1,
475 SrcIterator lowerright1, SrcAccessor sa,
476 DestIterator upperleft2, DestAccessor da,
477 int radius)
478 {
479 vigra_precondition(radius >= 0, "discDilation(): Radius must be >= 0.");
480
481 discRankOrderFilter(upperleft1, lowerright1, sa,
482 upperleft2, da, radius, 1.0);
483 }
484
485 template <class SrcIterator, class SrcAccessor,
486 class DestIterator, class DestAccessor>
487 void
discDilation(triple<SrcIterator,SrcIterator,SrcAccessor> src,pair<DestIterator,DestAccessor> dest,int radius)488 discDilation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
489 pair<DestIterator, DestAccessor> dest,
490 int radius)
491 {
492 vigra_precondition(radius >= 0, "discDilation(): Radius must be >= 0.");
493
494 discRankOrderFilter(src.first, src.second, src.third,
495 dest.first, dest.second,
496 radius, 1.0);
497 }
498
499 /********************************************************/
500 /* */
501 /* discMedian */
502 /* */
503 /********************************************************/
504
505 /** \brief Apply median filter with disc of given radius to image.
506
507 This is an abbreviation vor the rank order filter with rank = 0.5.
508 See \ref discRankOrderFilter() for more information.
509
510 <b> Declarations:</b>
511
512 pass arguments explicitely:
513 \code
514 namespace vigra {
515 template <class SrcIterator, class SrcAccessor,
516 class DestIterator, class DestAccessor>
517 inline void
518 discMedian(SrcIterator upperleft1,
519 SrcIterator lowerright1, SrcAccessor sa,
520 DestIterator upperleft2, DestAccessor da,
521 int radius)
522 }
523 \endcode
524
525
526 use argument objects in conjunction with \ref ArgumentObjectFactories:
527 \code
528 namespace vigra {
529 template <class SrcIterator, class SrcAccessor,
530 class DestIterator, class DestAccessor>
531 void
532 discMedian(triple<SrcIterator, SrcIterator, SrcAccessor> src,
533 pair<DestIterator, DestAccessor> dest,
534 int radius)
535 }
536 \endcode
537
538 */
539 template <class SrcIterator, class SrcAccessor,
540 class DestIterator, class DestAccessor>
541 inline void
discMedian(SrcIterator upperleft1,SrcIterator lowerright1,SrcAccessor sa,DestIterator upperleft2,DestAccessor da,int radius)542 discMedian(SrcIterator upperleft1,
543 SrcIterator lowerright1, SrcAccessor sa,
544 DestIterator upperleft2, DestAccessor da,
545 int radius)
546 {
547 vigra_precondition(radius >= 0, "discMedian(): Radius must be >= 0.");
548
549 discRankOrderFilter(upperleft1, lowerright1, sa,
550 upperleft2, da, radius, 0.5);
551 }
552
553 template <class SrcIterator, class SrcAccessor,
554 class DestIterator, class DestAccessor>
555 void
discMedian(triple<SrcIterator,SrcIterator,SrcAccessor> src,pair<DestIterator,DestAccessor> dest,int radius)556 discMedian(triple<SrcIterator, SrcIterator, SrcAccessor> src,
557 pair<DestIterator, DestAccessor> dest,
558 int radius)
559 {
560 vigra_precondition(radius >= 0, "discMedian(): Radius must be >= 0.");
561
562 discRankOrderFilter(src.first, src.second, src.third,
563 dest.first, dest.second,
564 radius, 0.5);
565 }
566
567 /********************************************************/
568 /* */
569 /* discRankOrderFilterWithMask */
570 /* */
571 /********************************************************/
572
573 /** \brief Apply rank order filter with disc structuring function to the image
574 using a mask.
575
576 The pixel values of the source image <b> must</b> be in the range
577 0...255. Radius must be >= 0. Rank must be in the range 0.0 <= rank
578 <= 1.0. The filter acts as a minimum filter if rank = 0.0,
579 as a median if rank = 0.5, and as a maximum filter if rank = 1.0.
580 Accessor are used to access the pixel data.
581
582 The mask is only applied to th input image, i.e. the function
583 generates an output wherever the current disc contains at least
584 one pixel with mask value 'true'. Source pixels with mask value
585 'false' are ignored during the calculation of the rank order.
586
587 <b> Declarations:</b>
588
589 pass arguments explicitely:
590 \code
591 namespace vigra {
592 template <class SrcIterator, class SrcAccessor,
593 class MaskIterator, class MaskAccessor,
594 class DestIterator, class DestAccessor>
595 void
596 discRankOrderFilterWithMask(SrcIterator upperleft1,
597 SrcIterator lowerright1, SrcAccessor sa,
598 MaskIterator upperleftm, MaskAccessor mask,
599 DestIterator upperleft2, DestAccessor da,
600 int radius, float rank)
601 }
602 \endcode
603
604
605 group arguments (use in conjunction with \ref ArgumentObjectFactories):
606 \code
607 namespace vigra {
608 template <class SrcIterator, class SrcAccessor,
609 class MaskIterator, class MaskAccessor,
610 class DestIterator, class DestAccessor>
611 void
612 discRankOrderFilterWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
613 pair<MaskIterator, MaskAccessor> mask,
614 pair<DestIterator, DestAccessor> dest,
615 int radius, float rank)
616 }
617 \endcode
618
619 <b> Usage:</b>
620
621 <b>\#include</b> "<a href="flatmorphology_8hxx-source.html">vigra/flatmorphology.hxx</a>"<br>
622 Namespace: vigra
623
624 \code
625 vigra::CImage src, dest, mask;
626
627 // do median filtering
628 vigra::discRankOrderFilterWithMask(srcImageRange(src),
629 maskImage(mask), destImage(dest), 10, 0.5);
630 \endcode
631
632 <b> Required Interface:</b>
633
634 \code
635 SrcIterator src_upperleft;
636 DestIterator dest_upperleft;
637 MaskIterator mask_upperleft;
638 int x, y;
639 unsigned char value;
640
641 SrcAccessor src_accessor;
642 DestAccessor dest_accessor;
643 MaskAccessor mask_accessor;
644
645 mask_accessor(mask_upperleft, x, y) // convertible to bool
646
647 // value_type of accessor must be convertible to unsigned char
648 value = src_accessor(src_upperleft, x, y);
649
650 dest_accessor.set(value, dest_upperleft, x, y);
651 \endcode
652
653 <b> Preconditions:</b>
654
655 \code
656 for all source pixels: 0 <= value <= 255
657
658 (rank >= 0.0) && (rank <= 1.0)
659 radius >= 0
660 \endcode
661
662 */
663 template <class SrcIterator, class SrcAccessor,
664 class MaskIterator, class MaskAccessor,
665 class DestIterator, class DestAccessor>
666 void
discRankOrderFilterWithMask(SrcIterator upperleft1,SrcIterator lowerright1,SrcAccessor sa,MaskIterator upperleftm,MaskAccessor mask,DestIterator upperleft2,DestAccessor da,int radius,float rank)667 discRankOrderFilterWithMask(SrcIterator upperleft1,
668 SrcIterator lowerright1, SrcAccessor sa,
669 MaskIterator upperleftm, MaskAccessor mask,
670 DestIterator upperleft2, DestAccessor da,
671 int radius, float rank)
672 {
673 vigra_precondition((rank >= 0.0) && (rank <= 1.0),
674 "discRankOrderFilter(): Rank must be between 0 and 1"
675 " (inclusive).");
676
677 vigra_precondition(radius >= 0, "discRankOrderFilter(): Radius must be >= 0.");
678
679 int i, x, y, xmax, ymax, xx, yy;
680 int rankpos, winsize, leftsum;
681
682 long hist[256];
683
684 // prepare structuring function
685 std::vector<int> struct_function(radius+1);
686 struct_function[0] = radius;
687
688 double r2 = (double)radius*radius;
689 for(i=1; i<=radius; ++i)
690 {
691 double r = (double) i - 0.5;
692 struct_function[i] = (int)(VIGRA_CSTD::sqrt(r2 - r*r) + 0.5);
693 }
694
695 int w = lowerright1.x - upperleft1.x;
696 int h = lowerright1.y - upperleft1.y;
697
698 SrcIterator ys(upperleft1);
699 MaskIterator ym(upperleftm);
700 DestIterator yd(upperleft2);
701
702 for(y=0; y<h; ++y, ++ys.y, ++yd.y, ++ym.y)
703 {
704 SrcIterator xs(ys);
705 MaskIterator xm(ym);
706 DestIterator xd(yd);
707
708 // first column
709 int x0 = 0;
710 int y0 = y;
711 int x1 = w - 1;
712 int y1 = h - y - 1;
713
714 // clear histogram
715 for(i=0; i<256; ++i) hist[i] = 0;
716 winsize = 0;
717 leftsum = 0;
718 rankpos = 0;
719
720 // init histogram
721 ymax = (y1 < radius) ? y1 : radius;
722 for(yy=0; yy<=ymax; ++yy)
723 {
724 xmax = (x1 < struct_function[yy]) ? x1 : struct_function[yy];
725 for(xx=0; xx<=xmax; ++xx)
726 {
727 Diff2D pos(xx, yy);
728 if(mask(xm, pos))
729 {
730 hist[sa(xs, pos)]++;
731 winsize++;
732 }
733 }
734 }
735
736 ymax = (y0 < radius) ? y0 : radius;
737 for(yy=1; yy<=ymax; ++yy)
738 {
739 xmax = (x1 < struct_function[yy]) ? x1 : struct_function[yy];
740 for(xx=0; xx<=xmax; ++xx)
741 {
742 Diff2D pos(xx, -yy);
743 if(mask(xm, pos))
744 {
745 hist[sa(xs, pos)]++;
746 winsize++;
747 }
748 }
749 }
750
751 // find the desired histogramm bin
752 if(winsize)
753 {
754 if(rank == 0.0)
755 {
756 for(i=0; i<256; i++)
757 {
758 if(hist[i]) break;
759 }
760 rankpos = i;
761 }
762 else
763 {
764 for(i=0; i<256; i++)
765 {
766 if((float)(hist[i]+leftsum) / winsize >= rank) break;
767 leftsum += hist[i];
768 }
769 rankpos = i;
770 }
771
772 da.set(rankpos, xd);
773 }
774
775 ++xs.x;
776 ++xd.x;
777 ++xm.x;
778
779 // inner columns
780 for(x=1; x<w; ++x, ++xs.x, ++xd.x, ++xm.x)
781 {
782 x0 = x;
783 y0 = y;
784 x1 = w - x - 1;
785 y1 = h - y - 1;
786
787 // update histogramm
788 // remove pixels at left border
789 yy = (y1 < radius) ? y1 : radius;
790 for(; yy>=0; yy--)
791 {
792 unsigned char cur;
793 xx = struct_function[yy]+1;
794 if(xx > x0) break;
795
796 Diff2D pos(-xx, yy);
797 if(mask(xm, pos))
798 {
799 cur = sa(xs, pos);
800
801 hist[cur]--;
802 if(cur < rankpos) leftsum--;
803 winsize--;
804 }
805 }
806 yy = (y0 < radius) ? y0 : radius;
807 for(; yy>=1; yy--)
808 {
809 unsigned char cur;
810 xx = struct_function[yy]+1;
811 if(xx > x0) break;
812
813 Diff2D pos(-xx, -yy);
814 if(mask(xm, pos))
815 {
816 cur = sa(xs, pos);
817
818 hist[cur]--;
819 if(cur < rankpos) leftsum--;
820 winsize--;
821 }
822 }
823
824 // add pixels at right border
825 yy = (y1 < radius) ? y1 : radius;
826 for(; yy>=0; yy--)
827 {
828 unsigned char cur;
829 xx = struct_function[yy];
830 if(xx > x1) break;
831
832 Diff2D pos(xx, yy);
833 if(mask(xm, pos))
834 {
835 cur = sa(xs, pos);
836
837 hist[cur]++;
838 if(cur < rankpos) leftsum++;
839 winsize++;
840 }
841 }
842 yy = (y0 < radius) ? y0 : radius;
843 for(; yy>=1; yy--)
844 {
845 unsigned char cur;
846 xx = struct_function[yy];
847 if(xx > x1) break;
848
849 Diff2D pos(xx, -yy);
850 if(mask(xm, pos))
851 {
852 cur = sa(xs, pos);
853
854 hist[cur]++;
855 if(cur < rankpos) leftsum++;
856 winsize++;
857 }
858 }
859
860 // find the desired histogramm bin
861 if(winsize)
862 {
863 if(rank == 0.0)
864 {
865 if(leftsum == 0)
866 {
867 // search to the right
868 for(i=rankpos; i<256; i++)
869 {
870 if(hist[i]) break;
871 }
872 rankpos = i;
873 }
874 else
875 {
876 // search to the left
877 for(i=rankpos-1; i>=0; i--)
878 {
879 leftsum -= hist[i];
880 if(leftsum == 0) break;
881 }
882 rankpos = i;
883 }
884 }
885 else // rank > 0.0
886 {
887 if((float)leftsum / winsize < rank)
888 {
889 // search to the right
890 for(i=rankpos; i<256; i++)
891 {
892 if((float)(hist[i]+leftsum) / winsize >= rank) break;
893 leftsum+=hist[i];
894 }
895 rankpos = i;
896 }
897 else
898 {
899 /// search to the left
900 for(i=rankpos-1; i>=0; i--)
901 {
902 leftsum-=hist[i];
903 if((float)leftsum / winsize < rank) break;
904 }
905 rankpos = i;
906 }
907 }
908
909 da.set(rankpos, xd);
910 }
911 else
912 {
913 leftsum = 0;
914 rankpos = 0;
915 }
916 }
917 }
918 }
919
920 template <class SrcIterator, class SrcAccessor,
921 class MaskIterator, class MaskAccessor,
922 class DestIterator, class DestAccessor>
923 void
discRankOrderFilterWithMask(triple<SrcIterator,SrcIterator,SrcAccessor> src,pair<MaskIterator,MaskAccessor> mask,pair<DestIterator,DestAccessor> dest,int radius,float rank)924 discRankOrderFilterWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
925 pair<MaskIterator, MaskAccessor> mask,
926 pair<DestIterator, DestAccessor> dest,
927 int radius, float rank)
928 {
929 discRankOrderFilterWithMask(src.first, src.second, src.third,
930 mask.first, mask.second,
931 dest.first, dest.second,
932 radius, rank);
933 }
934
935 /********************************************************/
936 /* */
937 /* discErosionWithMask */
938 /* */
939 /********************************************************/
940
941 /** \brief Apply erosion (minimum) filter with disc of given radius to image
942 using a mask.
943
944 This is an abbreviation vor the masked rank order filter with
945 rank = 0.0. See \ref discRankOrderFilterWithMask() for more information.
946
947 <b> Declarations:</b>
948
949 pass arguments explicitely:
950 \code
951 namespace vigra {
952 template <class SrcIterator, class SrcAccessor,
953 class MaskIterator, class MaskAccessor,
954 class DestIterator, class DestAccessor>
955 inline void
956 discErosionWithMask(SrcIterator upperleft1,
957 SrcIterator lowerright1, SrcAccessor sa,
958 MaskIterator upperleftm, MaskAccessor mask,
959 DestIterator upperleft2, DestAccessor da,
960 int radius)
961 }
962 \endcode
963
964
965 group arguments (use in conjunction with \ref ArgumentObjectFactories):
966 \code
967 namespace vigra {
968 template <class SrcIterator, class SrcAccessor,
969 class MaskIterator, class MaskAccessor,
970 class DestIterator, class DestAccessor>
971 inline void
972 discErosionWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
973 pair<MaskIterator, MaskAccessor> mask,
974 pair<DestIterator, DestAccessor> dest,
975 int radius)
976 }
977 \endcode
978
979 */
980 template <class SrcIterator, class SrcAccessor,
981 class MaskIterator, class MaskAccessor,
982 class DestIterator, class DestAccessor>
983 inline void
discErosionWithMask(SrcIterator upperleft1,SrcIterator lowerright1,SrcAccessor sa,MaskIterator upperleftm,MaskAccessor mask,DestIterator upperleft2,DestAccessor da,int radius)984 discErosionWithMask(SrcIterator upperleft1,
985 SrcIterator lowerright1, SrcAccessor sa,
986 MaskIterator upperleftm, MaskAccessor mask,
987 DestIterator upperleft2, DestAccessor da,
988 int radius)
989 {
990 vigra_precondition(radius >= 0, "discErosionWithMask(): Radius must be >= 0.");
991
992 discRankOrderFilterWithMask(upperleft1, lowerright1, sa,
993 upperleftm, mask,
994 upperleft2, da,
995 radius, 0.0);
996 }
997
998 template <class SrcIterator, class SrcAccessor,
999 class MaskIterator, class MaskAccessor,
1000 class DestIterator, class DestAccessor>
1001 inline void
discErosionWithMask(triple<SrcIterator,SrcIterator,SrcAccessor> src,pair<MaskIterator,MaskAccessor> mask,pair<DestIterator,DestAccessor> dest,int radius)1002 discErosionWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1003 pair<MaskIterator, MaskAccessor> mask,
1004 pair<DestIterator, DestAccessor> dest,
1005 int radius)
1006 {
1007 vigra_precondition(radius >= 0, "discErosionWithMask(): Radius must be >= 0.");
1008
1009 discRankOrderFilterWithMask(src.first, src.second, src.third,
1010 mask.first, mask.second,
1011 dest.first, dest.second,
1012 radius, 0.0);
1013 }
1014
1015 /********************************************************/
1016 /* */
1017 /* discDilationWithMask */
1018 /* */
1019 /********************************************************/
1020
1021 /** \brief Apply dilation (maximum) filter with disc of given radius to image
1022 using a mask.
1023
1024 This is an abbreviation vor the masked rank order filter with
1025 rank = 1.0. See \ref discRankOrderFilterWithMask() for more information.
1026
1027 <b> Declarations:</b>
1028
1029 pass arguments explicitely:
1030 \code
1031 namespace vigra {
1032 template <class SrcIterator, class SrcAccessor,
1033 class MaskIterator, class MaskAccessor,
1034 class DestIterator, class DestAccessor>
1035 inline void
1036 discDilationWithMask(SrcIterator upperleft1,
1037 SrcIterator lowerright1, SrcAccessor sa,
1038 MaskIterator upperleftm, MaskAccessor mask,
1039 DestIterator upperleft2, DestAccessor da,
1040 int radius)
1041 }
1042 \endcode
1043
1044
1045 group arguments (use in conjunction with \ref ArgumentObjectFactories):
1046 \code
1047 namespace vigra {
1048 template <class SrcIterator, class SrcAccessor,
1049 class MaskIterator, class MaskAccessor,
1050 class DestIterator, class DestAccessor>
1051 inline void
1052 discDilationWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1053 pair<MaskIterator, MaskAccessor> mask,
1054 pair<DestIterator, DestAccessor> dest,
1055 int radius)
1056 }
1057 \endcode
1058
1059 */
1060 template <class SrcIterator, class SrcAccessor,
1061 class MaskIterator, class MaskAccessor,
1062 class DestIterator, class DestAccessor>
1063 inline void
discDilationWithMask(SrcIterator upperleft1,SrcIterator lowerright1,SrcAccessor sa,MaskIterator upperleftm,MaskAccessor mask,DestIterator upperleft2,DestAccessor da,int radius)1064 discDilationWithMask(SrcIterator upperleft1,
1065 SrcIterator lowerright1, SrcAccessor sa,
1066 MaskIterator upperleftm, MaskAccessor mask,
1067 DestIterator upperleft2, DestAccessor da,
1068 int radius)
1069 {
1070 vigra_precondition(radius >= 0, "discDilationWithMask(): Radius must be >= 0.");
1071
1072 discRankOrderFilterWithMask(upperleft1, lowerright1, sa,
1073 upperleftm, mask,
1074 upperleft2, da,
1075 radius, 1.0);
1076 }
1077
1078 template <class SrcIterator, class SrcAccessor,
1079 class MaskIterator, class MaskAccessor,
1080 class DestIterator, class DestAccessor>
1081 inline void
discDilationWithMask(triple<SrcIterator,SrcIterator,SrcAccessor> src,pair<MaskIterator,MaskAccessor> mask,pair<DestIterator,DestAccessor> dest,int radius)1082 discDilationWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1083 pair<MaskIterator, MaskAccessor> mask,
1084 pair<DestIterator, DestAccessor> dest,
1085 int radius)
1086 {
1087 vigra_precondition(radius >= 0, "discDilationWithMask(): Radius must be >= 0.");
1088
1089 discRankOrderFilterWithMask(src.first, src.second, src.third,
1090 mask.first, mask.second,
1091 dest.first, dest.second,
1092 radius, 1.0);
1093 }
1094
1095 /********************************************************/
1096 /* */
1097 /* discMedianWithMask */
1098 /* */
1099 /********************************************************/
1100
1101 /** \brief Apply median filter with disc of given radius to image
1102 using a mask.
1103
1104 This is an abbreviation vor the masked rank order filter with
1105 rank = 0.5. See \ref discRankOrderFilterWithMask() for more information.
1106
1107 <b> Declarations:</b>
1108
1109 pass arguments explicitely:
1110 \code
1111 namespace vigra {
1112 template <class SrcIterator, class SrcAccessor,
1113 class MaskIterator, class MaskAccessor,
1114 class DestIterator, class DestAccessor>
1115 inline void
1116 discMedianWithMask(SrcIterator upperleft1,
1117 SrcIterator lowerright1, SrcAccessor sa,
1118 MaskIterator upperleftm, MaskAccessor mask,
1119 DestIterator upperleft2, DestAccessor da,
1120 int radius)
1121 }
1122 \endcode
1123
1124
1125 group arguments (use in conjunction with \ref ArgumentObjectFactories):
1126 \code
1127 namespace vigra {
1128 template <class SrcIterator, class SrcAccessor,
1129 class MaskIterator, class MaskAccessor,
1130 class DestIterator, class DestAccessor>
1131 inline void
1132 discMedianWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1133 pair<MaskIterator, MaskAccessor> mask,
1134 pair<DestIterator, DestAccessor> dest,
1135 int radius)
1136 }
1137 \endcode
1138
1139 */
1140 template <class SrcIterator, class SrcAccessor,
1141 class MaskIterator, class MaskAccessor,
1142 class DestIterator, class DestAccessor>
1143 inline void
discMedianWithMask(SrcIterator upperleft1,SrcIterator lowerright1,SrcAccessor sa,MaskIterator upperleftm,MaskAccessor mask,DestIterator upperleft2,DestAccessor da,int radius)1144 discMedianWithMask(SrcIterator upperleft1,
1145 SrcIterator lowerright1, SrcAccessor sa,
1146 MaskIterator upperleftm, MaskAccessor mask,
1147 DestIterator upperleft2, DestAccessor da,
1148 int radius)
1149 {
1150 vigra_precondition(radius >= 0, "discMedianWithMask(): Radius must be >= 0.");
1151
1152 discRankOrderFilterWithMask(upperleft1, lowerright1, sa,
1153 upperleftm, mask,
1154 upperleft2, da,
1155 radius, 0.5);
1156 }
1157
1158 template <class SrcIterator, class SrcAccessor,
1159 class MaskIterator, class MaskAccessor,
1160 class DestIterator, class DestAccessor>
1161 inline void
discMedianWithMask(triple<SrcIterator,SrcIterator,SrcAccessor> src,pair<MaskIterator,MaskAccessor> mask,pair<DestIterator,DestAccessor> dest,int radius)1162 discMedianWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1163 pair<MaskIterator, MaskAccessor> mask,
1164 pair<DestIterator, DestAccessor> dest,
1165 int radius)
1166 {
1167 vigra_precondition(radius >= 0, "discMedianWithMask(): Radius must be >= 0.");
1168
1169 discRankOrderFilterWithMask(src.first, src.second, src.third,
1170 mask.first, mask.second,
1171 dest.first, dest.second,
1172 radius, 0.5);
1173 }
1174
1175 //@}
1176
1177 } // namespace vigra
1178
1179 #endif // VIGRA_FLATMORPHOLOGY_HXX
1180