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

separableconvolution.hxx
1/************************************************************************/
2/* */
3/* Copyright 1998-2002 by Ullrich Koethe */
4/* */
5/* This file is part of the VIGRA computer vision library. */
6/* The VIGRA Website is */
7/* http://hci.iwr.uni-heidelberg.de/vigra/ */
8/* Please direct questions, bug reports, and contributions to */
9/* ullrich.koethe@iwr.uni-heidelberg.de or */
10/* vigra@informatik.uni-hamburg.de */
11/* */
12/* Permission is hereby granted, free of charge, to any person */
13/* obtaining a copy of this software and associated documentation */
14/* files (the "Software"), to deal in the Software without */
15/* restriction, including without limitation the rights to use, */
16/* copy, modify, merge, publish, distribute, sublicense, and/or */
17/* sell copies of the Software, and to permit persons to whom the */
18/* Software is furnished to do so, subject to the following */
19/* conditions: */
20/* */
21/* The above copyright notice and this permission notice shall be */
22/* included in all copies or substantial portions of the */
23/* Software. */
24/* */
25/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32/* OTHER DEALINGS IN THE SOFTWARE. */
33/* */
34/************************************************************************/
35
36
37#ifndef VIGRA_SEPARABLECONVOLUTION_HXX
38#define VIGRA_SEPARABLECONVOLUTION_HXX
39
40#include <cmath>
41#include "utilities.hxx"
42#include "numerictraits.hxx"
43#include "imageiteratoradapter.hxx"
44#include "bordertreatment.hxx"
45#include "gaussians.hxx"
46#include "array_vector.hxx"
47#include "multi_shape.hxx"
48
49namespace vigra {
50
51template <class ARITHTYPE>
52class Kernel1D;
53
54/********************************************************/
55/* */
56/* internalConvolveLineOptimistic */
57/* */
58/********************************************************/
59
60// This function assumes that the input array is actually larger than
61// the range [is, iend), so that it can safely access values outside
62// this range. This is useful if (1) we work on a small ROI, or
63// (2) we enlarge the input by copying with border treatment.
64template <class SrcIterator, class SrcAccessor,
65 class DestIterator, class DestAccessor,
66 class KernelIterator, class KernelAccessor>
67void internalConvolveLineOptimistic(SrcIterator is, SrcIterator iend, SrcAccessor sa,
68 DestIterator id, DestAccessor da,
69 KernelIterator kernel, KernelAccessor ka,
70 int kleft, int kright)
71{
72 typedef typename PromoteTraits<
73 typename SrcAccessor::value_type,
74 typename KernelAccessor::value_type>::Promote SumType;
75
76 int w = std::distance( is, iend );
77 int kw = kright - kleft + 1;
78 for(int x=0; x<w; ++x, ++is, ++id)
79 {
80 SrcIterator iss = is + (-kright);
81 KernelIterator ik = kernel + kright;
82 SumType sum = NumericTraits<SumType>::zero();
83
84 for(int k = 0; k < kw; ++k, --ik, ++iss)
85 {
86 sum += ka(ik) * sa(iss);
87 }
88
89 da.set(detail::RequiresExplicitCast<typename
90 DestAccessor::value_type>::cast(sum), id);
91 }
92}
93
94namespace detail {
95
96// dest array must have size = stop - start + kright - kleft
97template <class SrcIterator, class SrcAccessor,
98 class DestIterator, class DestAccessor>
99void
100copyLineWithBorderTreatment(SrcIterator is, SrcIterator iend, SrcAccessor sa,
101 DestIterator id, DestAccessor da,
102 int start, int stop,
103 int kleft, int kright,
104 BorderTreatmentMode borderTreatment)
105{
106 int w = std::distance( is, iend );
107 int leftBorder = start - kright;
108 int rightBorder = stop - kleft;
109 int copyEnd = std::min(w, rightBorder);
110
111 if(leftBorder < 0)
112 {
113 switch(borderTreatment)
114 {
115 case BORDER_TREATMENT_WRAP:
116 {
117 for(; leftBorder<0; ++leftBorder, ++id)
118 da.set(sa(iend, leftBorder), id);
119 break;
120 }
121 case BORDER_TREATMENT_AVOID:
122 {
123 // nothing to do
124 break;
125 }
126 case BORDER_TREATMENT_REFLECT:
127 {
128 for(; leftBorder<0; ++leftBorder, ++id)
129 da.set(sa(is, -leftBorder), id);
130 break;
131 }
132 case BORDER_TREATMENT_REPEAT:
133 {
134 for(; leftBorder<0; ++leftBorder, ++id)
135 da.set(sa(is), id);
136 break;
137 }
138 case BORDER_TREATMENT_CLIP:
139 {
140 vigra_precondition(false,
141 "copyLineWithBorderTreatment() internal error: not applicable to BORDER_TREATMENT_CLIP.");
142 break;
143 }
144 case BORDER_TREATMENT_ZEROPAD:
145 {
146 for(; leftBorder<0; ++leftBorder, ++id)
147 da.set(NumericTraits<typename DestAccessor::value_type>::zero(), id);
148 break;
149 }
150 default:
151 {
152 vigra_precondition(false,
153 "copyLineWithBorderTreatment(): Unknown border treatment mode.");
154 }
155 }
156 }
157
158 SrcIterator iss = is + leftBorder;
159 vigra_invariant( leftBorder < copyEnd,
160 "copyLineWithBorderTreatment(): assertion failed.");
161 for(; leftBorder<copyEnd; ++leftBorder, ++id, ++iss)
162 da.set(sa(iss), id);
163
164 if(copyEnd < rightBorder)
165 {
166 switch(borderTreatment)
167 {
168 case BORDER_TREATMENT_WRAP:
169 {
170 for(; copyEnd<rightBorder; ++copyEnd, ++id, ++is)
171 da.set(sa(is), id);
172 break;
173 }
174 case BORDER_TREATMENT_AVOID:
175 {
176 // nothing to do
177 break;
178 }
179 case BORDER_TREATMENT_REFLECT:
180 {
181 iss -= 2;
182 for(; copyEnd<rightBorder; ++copyEnd, ++id, --iss)
183 da.set(sa(iss), id);
184 break;
185 }
186 case BORDER_TREATMENT_REPEAT:
187 {
188 --iss;
189 for(; copyEnd<rightBorder; ++copyEnd, ++id)
190 da.set(sa(iss), id);
191 break;
192 }
193 case BORDER_TREATMENT_CLIP:
194 {
195 vigra_precondition(false,
196 "copyLineWithBorderTreatment() internal error: not applicable to BORDER_TREATMENT_CLIP.");
197 break;
198 }
199 case BORDER_TREATMENT_ZEROPAD:
200 {
201 for(; copyEnd<rightBorder; ++copyEnd, ++id)
202 da.set(NumericTraits<typename DestAccessor::value_type>::zero(), id);
203 break;
204 }
205 default:
206 {
207 vigra_precondition(false,
208 "copyLineWithBorderTreatment(): Unknown border treatment mode.");
209 }
210 }
211 }
212}
213
214} // namespace detail
215
216/********************************************************/
217/* */
218/* internalConvolveLineWrap */
219/* */
220/********************************************************/
221
222template <class SrcIterator, class SrcAccessor,
223 class DestIterator, class DestAccessor,
224 class KernelIterator, class KernelAccessor>
225void internalConvolveLineWrap(SrcIterator is, SrcIterator iend, SrcAccessor sa,
226 DestIterator id, DestAccessor da,
227 KernelIterator kernel, KernelAccessor ka,
228 int kleft, int kright,
229 int start = 0, int stop = 0)
230{
231 int w = std::distance( is, iend );
232
233 typedef typename PromoteTraits<
234 typename SrcAccessor::value_type,
235 typename KernelAccessor::value_type>::Promote SumType;
236
237 SrcIterator ibegin = is;
238
239 if(stop == 0)
240 stop = w;
241 is += start;
242
243 for(int x=start; x<stop; ++x, ++is, ++id)
244 {
245 KernelIterator ik = kernel + kright;
246 SumType sum = NumericTraits<SumType>::zero();
247
248 if(x < kright)
249 {
250 int x0 = x - kright;
251 SrcIterator iss = iend + x0;
252
253 for(; x0; ++x0, --ik, ++iss)
254 {
255 sum += ka(ik) * sa(iss);
256 }
257
258 iss = ibegin;
259 if(w-x <= -kleft)
260 {
261 SrcIterator isend = iend;
262 for(; iss != isend ; --ik, ++iss)
263 {
264 sum += ka(ik) * sa(iss);
265 }
266
267 int x0 = -kleft - w + x + 1;
268 iss = ibegin;
269
270 for(; x0; --x0, --ik, ++iss)
271 {
272 sum += ka(ik) * sa(iss);
273 }
274 }
275 else
276 {
277 SrcIterator isend = is + (1 - kleft);
278 for(; iss != isend ; --ik, ++iss)
279 {
280 sum += ka(ik) * sa(iss);
281 }
282 }
283 }
284 else if(w-x <= -kleft)
285 {
286 SrcIterator iss = is + (-kright);
287 SrcIterator isend = iend;
288 for(; iss != isend ; --ik, ++iss)
289 {
290 sum += ka(ik) * sa(iss);
291 }
292
293 int x0 = -kleft - w + x + 1;
294 iss = ibegin;
295
296 for(; x0; --x0, --ik, ++iss)
297 {
298 sum += ka(ik) * sa(iss);
299 }
300 }
301 else
302 {
303 SrcIterator iss = is - kright;
304 SrcIterator isend = is + (1 - kleft);
305 for(; iss != isend ; --ik, ++iss)
306 {
307 sum += ka(ik) * sa(iss);
308 }
309 }
310
311 da.set(detail::RequiresExplicitCast<typename
312 DestAccessor::value_type>::cast(sum), id);
313 }
314}
315
316/********************************************************/
317/* */
318/* internalConvolveLineClip */
319/* */
320/********************************************************/
321
322template <class SrcIterator, class SrcAccessor,
323 class DestIterator, class DestAccessor,
324 class KernelIterator, class KernelAccessor,
325 class Norm>
326void internalConvolveLineClip(SrcIterator is, SrcIterator iend, SrcAccessor sa,
327 DestIterator id, DestAccessor da,
328 KernelIterator kernel, KernelAccessor ka,
329 int kleft, int kright, Norm norm,
330 int start = 0, int stop = 0)
331{
332 int w = std::distance( is, iend );
333
334 typedef typename PromoteTraits<
335 typename SrcAccessor::value_type,
336 typename KernelAccessor::value_type>::Promote SumType;
337
338 SrcIterator ibegin = is;
339
340 if(stop == 0)
341 stop = w;
342 is += start;
343
344 for(int x=start; x<stop; ++x, ++is, ++id)
345 {
346 KernelIterator ik = kernel + kright;
347 SumType sum = NumericTraits<SumType>::zero();
348
349 if(x < kright)
350 {
351 int x0 = x - kright;
352 Norm clipped = NumericTraits<Norm>::zero();
353
354 for(; x0; ++x0, --ik)
355 {
356 clipped += ka(ik);
357 }
358
359 SrcIterator iss = ibegin;
360 if(w-x <= -kleft)
361 {
362 SrcIterator isend = iend;
363 for(; iss != isend ; --ik, ++iss)
364 {
365 sum += ka(ik) * sa(iss);
366 }
367
368 int x0 = -kleft - w + x + 1;
369
370 for(; x0; --x0, --ik)
371 {
372 clipped += ka(ik);
373 }
374 }
375 else
376 {
377 SrcIterator isend = is + (1 - kleft);
378 for(; iss != isend ; --ik, ++iss)
379 {
380 sum += ka(ik) * sa(iss);
381 }
382 }
383
384 sum = norm / (norm - clipped) * sum;
385 }
386 else if(w-x <= -kleft)
387 {
388 SrcIterator iss = is + (-kright);
389 SrcIterator isend = iend;
390 for(; iss != isend ; --ik, ++iss)
391 {
392 sum += ka(ik) * sa(iss);
393 }
394
395 Norm clipped = NumericTraits<Norm>::zero();
396
397 int x0 = -kleft - w + x + 1;
398
399 for(; x0; --x0, --ik)
400 {
401 clipped += ka(ik);
402 }
403
404 sum = norm / (norm - clipped) * sum;
405 }
406 else
407 {
408 SrcIterator iss = is + (-kright);
409 SrcIterator isend = is + (1 - kleft);
410 for(; iss != isend ; --ik, ++iss)
411 {
412 sum += ka(ik) * sa(iss);
413 }
414 }
415
416 da.set(detail::RequiresExplicitCast<typename
417 DestAccessor::value_type>::cast(sum), id);
418 }
419}
420
421/********************************************************/
422/* */
423/* internalConvolveLineZeropad */
424/* */
425/********************************************************/
426
427template <class SrcIterator, class SrcAccessor,
428 class DestIterator, class DestAccessor,
429 class KernelIterator, class KernelAccessor>
430void internalConvolveLineZeropad(SrcIterator is, SrcIterator iend, SrcAccessor sa,
431 DestIterator id, DestAccessor da,
432 KernelIterator kernel, KernelAccessor ka,
433 int kleft, int kright,
434 int start = 0, int stop = 0)
435{
436 int w = std::distance( is, iend );
437
438 typedef typename PromoteTraits<
439 typename SrcAccessor::value_type,
440 typename KernelAccessor::value_type>::Promote SumType;
441
442 SrcIterator ibegin = is;
443
444 if(stop == 0)
445 stop = w;
446 is += start;
447
448 for(int x=start; x<stop; ++x, ++is, ++id)
449 {
450 SumType sum = NumericTraits<SumType>::zero();
451
452 if(x < kright)
453 {
454 KernelIterator ik = kernel + x;
455 SrcIterator iss = ibegin;
456
457 if(w-x <= -kleft)
458 {
459 SrcIterator isend = iend;
460 for(; iss != isend ; --ik, ++iss)
461 {
462 sum += ka(ik) * sa(iss);
463 }
464 }
465 else
466 {
467 SrcIterator isend = is + (1 - kleft);
468 for(; iss != isend ; --ik, ++iss)
469 {
470 sum += ka(ik) * sa(iss);
471 }
472 }
473 }
474 else if(w-x <= -kleft)
475 {
476 KernelIterator ik = kernel + kright;
477 SrcIterator iss = is + (-kright);
478 SrcIterator isend = iend;
479 for(; iss != isend ; --ik, ++iss)
480 {
481 sum += ka(ik) * sa(iss);
482 }
483 }
484 else
485 {
486 KernelIterator ik = kernel + kright;
487 SrcIterator iss = is + (-kright);
488 SrcIterator isend = is + (1 - kleft);
489 for(; iss != isend ; --ik, ++iss)
490 {
491 sum += ka(ik) * sa(iss);
492 }
493 }
494
495 da.set(detail::RequiresExplicitCast<typename
496 DestAccessor::value_type>::cast(sum), id);
497 }
498}
499
500/********************************************************/
501/* */
502/* internalConvolveLineReflect */
503/* */
504/********************************************************/
505
506template <class SrcIterator, class SrcAccessor,
507 class DestIterator, class DestAccessor,
508 class KernelIterator, class KernelAccessor>
509void internalConvolveLineReflect(SrcIterator is, SrcIterator iend, SrcAccessor sa,
510 DestIterator id, DestAccessor da,
511 KernelIterator kernel, KernelAccessor ka,
512 int kleft, int kright,
513 int start = 0, int stop = 0)
514{
515 int w = std::distance( is, iend );
516
517 typedef typename PromoteTraits<
518 typename SrcAccessor::value_type,
519 typename KernelAccessor::value_type>::Promote SumType;
520
521 SrcIterator ibegin = is;
522
523 if(stop == 0)
524 stop = w;
525 is += start;
526
527 for(int x=start; x<stop; ++x, ++is, ++id)
528 {
529 KernelIterator ik = kernel + kright;
530 SumType sum = NumericTraits<SumType>::zero();
531
532 if(x < kright)
533 {
534 int x0 = x - kright;
535 SrcIterator iss = ibegin - x0;
536
537 for(; x0; ++x0, --ik, --iss)
538 {
539 sum += ka(ik) * sa(iss);
540 }
541
542 if(w-x <= -kleft)
543 {
544 SrcIterator isend = iend;
545 for(; iss != isend ; --ik, ++iss)
546 {
547 sum += ka(ik) * sa(iss);
548 }
549
550 int x0 = -kleft - w + x + 1;
551 iss = iend - 2;
552
553 for(; x0; --x0, --ik, --iss)
554 {
555 sum += ka(ik) * sa(iss);
556 }
557 }
558 else
559 {
560 SrcIterator isend = is + (1 - kleft);
561 for(; iss != isend ; --ik, ++iss)
562 {
563 sum += ka(ik) * sa(iss);
564 }
565 }
566 }
567 else if(w-x <= -kleft)
568 {
569 SrcIterator iss = is + (-kright);
570 SrcIterator isend = iend;
571 for(; iss != isend ; --ik, ++iss)
572 {
573 sum += ka(ik) * sa(iss);
574 }
575
576 int x0 = -kleft - w + x + 1;
577 iss = iend - 2;
578
579 for(; x0; --x0, --ik, --iss)
580 {
581 sum += ka(ik) * sa(iss);
582 }
583 }
584 else
585 {
586 SrcIterator iss = is + (-kright);
587 SrcIterator isend = is + (1 - kleft);
588 for(; iss != isend ; --ik, ++iss)
589 {
590 sum += ka(ik) * sa(iss);
591 }
592 }
593
594 da.set(detail::RequiresExplicitCast<typename
595 DestAccessor::value_type>::cast(sum), id);
596 }
597}
598
599/********************************************************/
600/* */
601/* internalConvolveLineRepeat */
602/* */
603/********************************************************/
604
605template <class SrcIterator, class SrcAccessor,
606 class DestIterator, class DestAccessor,
607 class KernelIterator, class KernelAccessor>
608void internalConvolveLineRepeat(SrcIterator is, SrcIterator iend, SrcAccessor sa,
609 DestIterator id, DestAccessor da,
610 KernelIterator kernel, KernelAccessor ka,
611 int kleft, int kright,
612 int start = 0, int stop = 0)
613{
614 int w = std::distance( is, iend );
615
616 typedef typename PromoteTraits<
617 typename SrcAccessor::value_type,
618 typename KernelAccessor::value_type>::Promote SumType;
619
620 SrcIterator ibegin = is;
621
622 if(stop == 0)
623 stop = w;
624 is += start;
625
626 for(int x=start; x<stop; ++x, ++is, ++id)
627 {
628 KernelIterator ik = kernel + kright;
629 SumType sum = NumericTraits<SumType>::zero();
630
631 if(x < kright)
632 {
633 int x0 = x - kright;
634 SrcIterator iss = ibegin;
635
636 for(; x0; ++x0, --ik)
637 {
638 sum += ka(ik) * sa(iss);
639 }
640
641 if(w-x <= -kleft)
642 {
643 SrcIterator isend = iend;
644 for(; iss != isend ; --ik, ++iss)
645 {
646 sum += ka(ik) * sa(iss);
647 }
648
649 int x0 = -kleft - w + x + 1;
650 iss = iend - 1;
651
652 for(; x0; --x0, --ik)
653 {
654 sum += ka(ik) * sa(iss);
655 }
656 }
657 else
658 {
659 SrcIterator isend = is + (1 - kleft);
660 for(; iss != isend ; --ik, ++iss)
661 {
662 sum += ka(ik) * sa(iss);
663 }
664 }
665 }
666 else if(w-x <= -kleft)
667 {
668 SrcIterator iss = is + (-kright);
669 SrcIterator isend = iend;
670 for(; iss != isend ; --ik, ++iss)
671 {
672 sum += ka(ik) * sa(iss);
673 }
674
675 int x0 = -kleft - w + x + 1;
676 iss = iend - 1;
677
678 for(; x0; --x0, --ik)
679 {
680 sum += ka(ik) * sa(iss);
681 }
682 }
683 else
684 {
685 SrcIterator iss = is + (-kright);
686 SrcIterator isend = is + (1 - kleft);
687 for(; iss != isend ; --ik, ++iss)
688 {
689 sum += ka(ik) * sa(iss);
690 }
691 }
692
693 da.set(detail::RequiresExplicitCast<typename
694 DestAccessor::value_type>::cast(sum), id);
695 }
696}
697
698/********************************************************/
699/* */
700/* internalConvolveLineAvoid */
701/* */
702/********************************************************/
703
704template <class SrcIterator, class SrcAccessor,
705 class DestIterator, class DestAccessor,
706 class KernelIterator, class KernelAccessor>
707void internalConvolveLineAvoid(SrcIterator is, SrcIterator iend, SrcAccessor sa,
708 DestIterator id, DestAccessor da,
709 KernelIterator kernel, KernelAccessor ka,
710 int kleft, int kright,
711 int start = 0, int stop = 0)
712{
713 int w = std::distance( is, iend );
714 if(start < stop) // we got a valid subrange
715 {
716 if(w + kleft < stop)
717 stop = w + kleft;
718 if(start < kright)
719 {
720 id += kright - start;
721 start = kright;
722 }
723 }
724 else
725 {
726 id += kright;
727 start = kright;
728 stop = w + kleft;
729 }
730
731 typedef typename PromoteTraits<
732 typename SrcAccessor::value_type,
733 typename KernelAccessor::value_type>::Promote SumType;
734
735 is += start;
736
737 for(int x=start; x<stop; ++x, ++is, ++id)
738 {
739 KernelIterator ik = kernel + kright;
740 SumType sum = NumericTraits<SumType>::zero();
741
742 SrcIterator iss = is + (-kright);
743 SrcIterator isend = is + (1 - kleft);
744 for(; iss != isend ; --ik, ++iss)
745 {
746 sum += ka(ik) * sa(iss);
747 }
748
749 da.set(detail::RequiresExplicitCast<typename
750 DestAccessor::value_type>::cast(sum), id);
751 }
752}
753
754/********************************************************/
755/* */
756/* Separable convolution functions */
757/* */
758/********************************************************/
759
760/** \addtogroup SeparableConvolution One-dimensional and separable convolution functions
761
762 Perform 1D convolution and separable filtering in 2 dimensions.
763
764 These generic convolution functions implement
765 the standard convolution operation for a wide range of images and
766 signals that fit into the required interface. They need a suitable
767 kernel to operate.
768*/
769//@{
770
771/** \brief Performs a 1-dimensional convolution of the source signal using the given
772 kernel.
773
774 The KernelIterator must point to the center iterator, and
775 the kernel's size is given by its left (kleft <= 0) and right
776 (kright >= 0) borders. The signal must always be larger than the kernel.
777 At those positions where the kernel does not completely fit
778 into the signal's range, the specified \ref BorderTreatmentMode is
779 applied.
780
781 The signal's value_type (SrcAccessor::value_type) must be a
782 linear space over the kernel's value_type (KernelAccessor::value_type),
783 i.e. addition of source values, multiplication with kernel values,
784 and NumericTraits must be defined.
785 The kernel's value_type must be an algebraic field,
786 i.e. the arithmetic operations (+, -, *, /) and NumericTraits must
787 be defined.
788
789 If <tt>start</tt> and <tt>stop</tt> are non-zero, the relation
790 <tt>0 <= start < stop <= width</tt> must hold (where <tt>width</tt>
791 is the length of the input array). The convolution is then restricted to that
792 subrange, and it is assumed that the output array only refers to that
793 subrange (i.e. <tt>id</tt> points to the element corresponding to
794 <tt>start</tt>). If <tt>start</tt> and <tt>stop</tt> are both zero
795 (the default), the entire array is convolved.
796
797 <b> Declarations:</b>
798
799 pass \ref ImageIterators and \ref DataAccessors :
800 \code
801 namespace vigra {
802 template <class SrcIterator, class SrcAccessor,
803 class DestIterator, class DestAccessor,
804 class KernelIterator, class KernelAccessor>
805 void convolveLine(SrcIterator is, SrcIterator isend, SrcAccessor sa,
806 DestIterator id, DestAccessor da,
807 KernelIterator ik, KernelAccessor ka,
808 int kleft, int kright, BorderTreatmentMode border,
809 int start = 0, int stop = 0 )
810 }
811 \endcode
812 use argument objects in conjunction with \ref ArgumentObjectFactories :
813 \code
814 namespace vigra {
815 template <class SrcIterator, class SrcAccessor,
816 class DestIterator, class DestAccessor,
817 class KernelIterator, class KernelAccessor>
818 void convolveLine(triple<SrcIterator, SrcIterator, SrcAccessor> src,
819 pair<DestIterator, DestAccessor> dest,
820 tuple5<KernelIterator, KernelAccessor,
821 int, int, BorderTreatmentMode> kernel,
822 int start = 0, int stop = 0)
823 }
824 \endcode
825
826 <b> Usage:</b>
827
828 <b>\#include</b> <vigra/separableconvolution.hxx><br/>
829 Namespace: vigra
830
831
832 \code
833 std::vector<float> src, dest;
834 ...
835
836 // define binomial filter of size 5
837 static float kernel[] =
838 { 1.0/16.0, 4.0/16.0, 6.0/16.0, 4.0/16.0, 1.0/16.0};
839
840 typedef vigra::StandardAccessor<float> FAccessor;
841 typedef vigra::StandardAccessor<float> KernelAccessor;
842
843
844 vigra::convolveLine(src.begin(), src.end(), FAccessor(), dest.begin(), FAccessor(),
845 kernel+2, KernelAccessor(), -2, 2, BORDER_TREATMENT_REFLECT);
846 // ^^^^^^^^ this is the center of the kernel
847
848 \endcode
849
850 <b> Required Interface:</b>
851
852 \code
853 RandomAccessIterator is, isend;
854 RandomAccessIterator id;
855 RandomAccessIterator ik;
856
857 SrcAccessor src_accessor;
858 DestAccessor dest_accessor;
859 KernelAccessor kernel_accessor;
860
861 NumericTraits<SrcAccessor::value_type>::RealPromote s = src_accessor(is);
862
863 s = s + s;
864 s = kernel_accessor(ik) * s;
865
866 dest_accessor.set(
867 NumericTraits<DestAccessor::value_type>::fromRealPromote(s), id);
868
869 \endcode
870
871 If border == BORDER_TREATMENT_CLIP:
872
873 \code
874 NumericTraits<KernelAccessor::value_type>::RealPromote k = kernel_accessor(ik);
875
876 k = k + k;
877 k = k - k;
878 k = k * k;
879 k = k / k;
880
881 \endcode
882
883 <b> Preconditions:</b>
884
885 \code
886 kleft <= 0
887 kright >= 0
888 iend - is >= kright + kleft + 1
889 \endcode
890
891 If border == BORDER_TREATMENT_CLIP: Sum of kernel elements must be
892 != 0.
893*/
895
896template <class SrcIterator, class SrcAccessor,
897 class DestIterator, class DestAccessor,
898 class KernelIterator, class KernelAccessor>
899void convolveLine(SrcIterator is, SrcIterator iend, SrcAccessor sa,
900 DestIterator id, DestAccessor da,
901 KernelIterator ik, KernelAccessor ka,
902 int kleft, int kright, BorderTreatmentMode border,
903 int start = 0, int stop = 0)
904{
905 vigra_precondition(kleft <= 0,
906 "convolveLine(): kleft must be <= 0.\n");
907 vigra_precondition(kright >= 0,
908 "convolveLine(): kright must be >= 0.\n");
909
910 // int w = iend - is;
911 int w = std::distance( is, iend );
912
913 vigra_precondition(w >= std::max(kright, -kleft) + 1,
914 "convolveLine(): kernel longer than line.\n");
915
916 if(stop != 0)
917 vigra_precondition(0 <= start && start < stop && stop <= w,
918 "convolveLine(): invalid subrange (start, stop).\n");
919
920 typedef typename PromoteTraits<
921 typename SrcAccessor::value_type,
922 typename KernelAccessor::value_type>::Promote SumType;
923 ArrayVector<SumType> a(iend - is);
924 switch(border)
925 {
926 case BORDER_TREATMENT_WRAP:
927 {
928 internalConvolveLineWrap(is, iend, sa, id, da, ik, ka, kleft, kright, start, stop);
929 break;
930 }
931 case BORDER_TREATMENT_AVOID:
932 {
933 internalConvolveLineAvoid(is, iend, sa, id, da, ik, ka, kleft, kright, start, stop);
934 break;
935 }
936 case BORDER_TREATMENT_REFLECT:
937 {
938 internalConvolveLineReflect(is, iend, sa, id, da, ik, ka, kleft, kright, start, stop);
939 break;
940 }
941 case BORDER_TREATMENT_REPEAT:
942 {
943 internalConvolveLineRepeat(is, iend, sa, id, da, ik, ka, kleft, kright, start, stop);
944 break;
945 }
946 case BORDER_TREATMENT_CLIP:
947 {
948 // find norm of kernel
949 typedef typename KernelAccessor::value_type KT;
950 KT norm = NumericTraits<KT>::zero();
951 KernelIterator iik = ik + kleft;
952 for(int i=kleft; i<=kright; ++i, ++iik)
953 norm += ka(iik);
954
955 vigra_precondition(norm != NumericTraits<KT>::zero(),
956 "convolveLine(): Norm of kernel must be != 0"
957 " in mode BORDER_TREATMENT_CLIP.\n");
958
959 internalConvolveLineClip(is, iend, sa, id, da, ik, ka, kleft, kright, norm, start, stop);
960 break;
961 }
962 case BORDER_TREATMENT_ZEROPAD:
963 {
964 internalConvolveLineZeropad(is, iend, sa, id, da, ik, ka, kleft, kright, start, stop);
965 break;
966 }
967 default:
968 {
969 vigra_precondition(0,
970 "convolveLine(): Unknown border treatment mode.\n");
971 }
972 }
973}
974
975template <class SrcIterator, class SrcAccessor,
976 class DestIterator, class DestAccessor,
977 class KernelIterator, class KernelAccessor>
978inline
979void convolveLine(triple<SrcIterator, SrcIterator, SrcAccessor> src,
980 pair<DestIterator, DestAccessor> dest,
981 tuple5<KernelIterator, KernelAccessor,
982 int, int, BorderTreatmentMode> kernel,
983 int start = 0, int stop = 0)
984{
985 convolveLine(src.first, src.second, src.third,
986 dest.first, dest.second,
987 kernel.first, kernel.second,
988 kernel.third, kernel.fourth, kernel.fifth, start, stop);
989}
990
991/********************************************************/
992/* */
993/* separableConvolveX */
994/* */
995/********************************************************/
996
997/** \brief Performs a 1 dimensional convolution in x direction.
998
999 It calls \ref convolveLine() for every row of the image. See \ref convolveLine()
1000 for more information about required interfaces and vigra_preconditions.
1001
1002 <b> Declarations:</b>
1003
1004 pass 2D array views:
1005 \code
1006 namespace vigra {
1007 template <class T1, class S1,
1008 class T2, class S2,
1009 class T3>
1010 void
1011 separableConvolveX(MultiArrayView<2, T1, S1> const & src,
1012 MultiArrayView<2, T2, S2> dest,
1013 Kernel1D<T3> const & kernel);
1014 }
1015 \endcode
1016
1017 \deprecatedAPI{separableConvolveX}
1018 pass \ref ImageIterators and \ref DataAccessors :
1019 \code
1020 namespace vigra {
1021 template <class SrcImageIterator, class SrcAccessor,
1022 class DestImageIterator, class DestAccessor,
1023 class KernelIterator, class KernelAccessor>
1024 void separableConvolveX(SrcImageIterator supperleft,
1025 SrcImageIterator slowerright, SrcAccessor sa,
1026 DestImageIterator dupperleft, DestAccessor da,
1027 KernelIterator ik, KernelAccessor ka,
1028 int kleft, int kright, BorderTreatmentMode border)
1029 }
1030 \endcode
1031 use argument objects in conjunction with \ref ArgumentObjectFactories :
1032 \code
1033 namespace vigra {
1034 template <class SrcImageIterator, class SrcAccessor,
1035 class DestImageIterator, class DestAccessor,
1036 class KernelIterator, class KernelAccessor>
1037 void separableConvolveX(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1038 pair<DestImageIterator, DestAccessor> dest,
1039 tuple5<KernelIterator, KernelAccessor,
1040 int, int, BorderTreatmentMode> kernel)
1041 }
1042 \endcode
1043 \deprecatedEnd
1044
1045 <b> Usage:</b>
1046
1047 <b>\#include</b> <vigra/separableconvolution.hxx><br/>
1048 Namespace: vigra
1049
1050 \code
1051 MultiArray<2, float> src(w,h), dest(w,h);
1052 ...
1053
1054 // define Gaussian kernel with std. deviation 3.0
1055 Kernel1D<double> kernel;
1056 kernel.initGaussian(3.0);
1057
1058 // apply 1D filter along the x-axis
1059 separableConvolveX(src, dest, kernel);
1060 \endcode
1061
1062 \deprecatedUsage{separableConvolveX}
1063 \code
1064 vigra::FImage src(w,h), dest(w,h);
1065 ...
1066
1067 // define Gaussian kernel with std. deviation 3.0
1068 vigra::Kernel1D<double> kernel;
1069 kernel.initGaussian(3.0);
1070
1071 // apply 1D filter along the x-axis
1072 vigra::separableConvolveX(srcImageRange(src), destImage(dest), kernel1d(kernel));
1073 \endcode
1074 \deprecatedEnd
1075
1076 <b>Preconditions:</b>
1077
1078 <ul>
1079 <li> The x-axis must be longer than the kernel radius: <tt>w > std::max(kernel.right(), -kernel.left())</tt>.
1080 <li> If <tt>border == BORDER_TREATMENT_CLIP</tt>: The sum of kernel elements must be != 0.
1081 </ul>
1082*/
1084
1085template <class SrcIterator, class SrcAccessor,
1086 class DestIterator, class DestAccessor,
1087 class KernelIterator, class KernelAccessor>
1088void separableConvolveX(SrcIterator supperleft,
1089 SrcIterator slowerright, SrcAccessor sa,
1090 DestIterator dupperleft, DestAccessor da,
1091 KernelIterator ik, KernelAccessor ka,
1092 int kleft, int kright, BorderTreatmentMode border)
1093{
1094 vigra_precondition(kleft <= 0,
1095 "separableConvolveX(): kleft must be <= 0.\n");
1096 vigra_precondition(kright >= 0,
1097 "separableConvolveX(): kright must be >= 0.\n");
1098
1099 int w = slowerright.x - supperleft.x;
1100 int h = slowerright.y - supperleft.y;
1101
1102 vigra_precondition(w >= std::max(kright, -kleft) + 1,
1103 "separableConvolveX(): kernel longer than line\n");
1104
1105 int y;
1106
1107 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
1108 {
1109 typename SrcIterator::row_iterator rs = supperleft.rowIterator();
1110 typename DestIterator::row_iterator rd = dupperleft.rowIterator();
1111
1112 convolveLine(rs, rs+w, sa, rd, da,
1113 ik, ka, kleft, kright, border);
1114 }
1115}
1116
1117template <class SrcIterator, class SrcAccessor,
1118 class DestIterator, class DestAccessor,
1119 class KernelIterator, class KernelAccessor>
1120inline void
1121separableConvolveX(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1122 pair<DestIterator, DestAccessor> dest,
1123 tuple5<KernelIterator, KernelAccessor,
1124 int, int, BorderTreatmentMode> kernel)
1125{
1126 separableConvolveX(src.first, src.second, src.third,
1127 dest.first, dest.second,
1128 kernel.first, kernel.second,
1129 kernel.third, kernel.fourth, kernel.fifth);
1130}
1131
1132template <class T1, class S1,
1133 class T2, class S2,
1134 class T3>
1135inline void
1136separableConvolveX(MultiArrayView<2, T1, S1> const & src,
1137 MultiArrayView<2, T2, S2> dest,
1138 Kernel1D<T3> const & kernel)
1139{
1140 vigra_precondition(src.shape() == dest.shape(),
1141 "separableConvolveX(): shape mismatch between input and output.");
1142 separableConvolveX(srcImageRange(src),
1143 destImage(dest), kernel1d(kernel));
1144}
1145
1146/********************************************************/
1147/* */
1148/* separableConvolveY */
1149/* */
1150/********************************************************/
1151
1152/** \brief Performs a 1 dimensional convolution in y direction.
1153
1154 It calls \ref convolveLine() for every column of the image. See \ref convolveLine()
1155 for more information about required interfaces and vigra_preconditions.
1156
1157 <b> Declarations:</b>
1158
1159 pass 2D array views:
1160 \code
1161 namespace vigra {
1162 template <class T1, class S1,
1163 class T2, class S2,
1164 class T3>
1165 void
1166 separableConvolveY(MultiArrayView<2, T1, S1> const & src,
1167 MultiArrayView<2, T2, S2> dest,
1168 Kernel1D<T3> const & kernel);
1169 }
1170 \endcode
1171
1172 \deprecatedAPI{separableConvolveY}
1173 pass \ref ImageIterators and \ref DataAccessors :
1174 \code
1175 namespace vigra {
1176 template <class SrcImageIterator, class SrcAccessor,
1177 class DestImageIterator, class DestAccessor,
1178 class KernelIterator, class KernelAccessor>
1179 void separableConvolveY(SrcImageIterator supperleft,
1180 SrcImageIterator slowerright, SrcAccessor sa,
1181 DestImageIterator dupperleft, DestAccessor da,
1182 KernelIterator ik, KernelAccessor ka,
1183 int kleft, int kright, BorderTreatmentMode border)
1184 }
1185 \endcode
1186 use argument objects in conjunction with \ref ArgumentObjectFactories :
1187 \code
1188 namespace vigra {
1189 template <class SrcImageIterator, class SrcAccessor,
1190 class DestImageIterator, class DestAccessor,
1191 class KernelIterator, class KernelAccessor>
1192 void separableConvolveY(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1193 pair<DestImageIterator, DestAccessor> dest,
1194 tuple5<KernelIterator, KernelAccessor,
1195 int, int, BorderTreatmentMode> kernel)
1196 }
1197 \endcode
1198 \deprecatedEnd
1199
1200 <b> Usage:</b>
1201
1202 <b>\#include</b> <vigra/separableconvolution.hxx><br/>
1203 Namespace: vigra
1204
1205
1206 \code
1207 MultiArray<2, float> src(w,h), dest(w,h);
1208 ...
1209
1210 // define Gaussian kernel with std. deviation 3.0
1211 Kernel1D kernel;
1212 kernel.initGaussian(3.0);
1213
1214 // apply 1D filter along the y-axis
1215 separableConvolveY(src, dest, kernel);
1216 \endcode
1217
1218 \deprecatedUsage{separableConvolveY}
1219 \code
1220 vigra::FImage src(w,h), dest(w,h);
1221 ...
1222
1223 // define Gaussian kernel with std. deviation 3.0
1224 vigra::Kernel1D kernel;
1225 kernel.initGaussian(3.0);
1226
1227 vigra::separableConvolveY(srcImageRange(src), destImage(dest), kernel1d(kernel));
1228 \endcode
1229 \deprecatedEnd
1230
1231 <b>Preconditions:</b>
1232
1233 <ul>
1234 <li> The y-axis must be longer than the kernel radius: <tt>h > std::max(kernel.right(), -kernel.left())</tt>.
1235 <li> If <tt>border == BORDER_TREATMENT_CLIP</tt>: The sum of kernel elements must be != 0.
1236 </ul>
1237*/
1239
1240template <class SrcIterator, class SrcAccessor,
1241 class DestIterator, class DestAccessor,
1242 class KernelIterator, class KernelAccessor>
1243void separableConvolveY(SrcIterator supperleft,
1244 SrcIterator slowerright, SrcAccessor sa,
1245 DestIterator dupperleft, DestAccessor da,
1246 KernelIterator ik, KernelAccessor ka,
1247 int kleft, int kright, BorderTreatmentMode border)
1248{
1249 vigra_precondition(kleft <= 0,
1250 "separableConvolveY(): kleft must be <= 0.\n");
1251 vigra_precondition(kright >= 0,
1252 "separableConvolveY(): kright must be >= 0.\n");
1253
1254 int w = slowerright.x - supperleft.x;
1255 int h = slowerright.y - supperleft.y;
1256
1257 vigra_precondition(h >= std::max(kright, -kleft) + 1,
1258 "separableConvolveY(): kernel longer than line\n");
1259
1260 int x;
1261
1262 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1263 {
1264 typename SrcIterator::column_iterator cs = supperleft.columnIterator();
1265 typename DestIterator::column_iterator cd = dupperleft.columnIterator();
1266
1267 convolveLine(cs, cs+h, sa, cd, da,
1268 ik, ka, kleft, kright, border);
1269 }
1270}
1271
1272template <class SrcIterator, class SrcAccessor,
1273 class DestIterator, class DestAccessor,
1274 class KernelIterator, class KernelAccessor>
1275inline void
1276separableConvolveY(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1277 pair<DestIterator, DestAccessor> dest,
1278 tuple5<KernelIterator, KernelAccessor,
1279 int, int, BorderTreatmentMode> kernel)
1280{
1281 separableConvolveY(src.first, src.second, src.third,
1282 dest.first, dest.second,
1283 kernel.first, kernel.second,
1284 kernel.third, kernel.fourth, kernel.fifth);
1285}
1286
1287template <class T1, class S1,
1288 class T2, class S2,
1289 class T3>
1290inline void
1291separableConvolveY(MultiArrayView<2, T1, S1> const & src,
1292 MultiArrayView<2, T2, S2> dest,
1293 Kernel1D<T3> const & kernel)
1294{
1295 vigra_precondition(src.shape() == dest.shape(),
1296 "separableConvolveY(): shape mismatch between input and output.");
1297 separableConvolveY(srcImageRange(src),
1298 destImage(dest), kernel1d(kernel));
1299}
1300
1301//@}
1302
1303/********************************************************/
1304/* */
1305/* Kernel1D */
1306/* */
1307/********************************************************/
1308
1309/** \brief Generic 1 dimensional convolution kernel.
1310
1311 This kernel may be used for convolution of 1 dimensional signals or for
1312 separable convolution of multidimensional signals.
1313
1314 Convolution functions access the kernel via a 1 dimensional random access
1315 iterator which they get by calling \ref center(). This iterator
1316 points to the center of the kernel. The kernel's size is given by its left() (<=0)
1317 and right() (>= 0) methods. The desired border treatment mode is
1318 returned by borderTreatment().
1319
1320 The different init functions create a kernel with the specified
1321 properties. The kernel's value_type must be a linear space, i.e. it
1322 must define multiplication with doubles and NumericTraits.
1323
1324 <b> Usage:</b>
1325
1326 <b>\#include</b> <vigra/separableconvolution.hxx><br/>
1327 Namespace: vigra
1328
1329 \code
1330 MultiArray<2, float> src(w,h), dest(w,h);
1331 ...
1332
1333 // define Gaussian kernel with std. deviation 3.0
1334 Kernel1D kernel;
1335 kernel.initGaussian(3.0);
1336
1337 // apply 1D kernel along the x-axis
1338 separableConvolveX(src, dest, kernel);
1339 \endcode
1340
1341 \deprecatedUsage{Kernel1D}
1342 The kernel defines a factory function kernel1d() to create an argument object
1343 (see \ref KernelArgumentObjectFactories).
1344 \code
1345 vigra::FImage src(w,h), dest(w,h);
1346 ...
1347
1348 // define Gaussian kernel with std. deviation 3.0
1349 vigra::Kernel1D kernel;
1350 kernel.initGaussian(3.0);
1351
1352 vigra::separableConvolveX(srcImageRange(src), destImage(dest), kernel1d(kernel));
1353 \endcode
1354 <b> Required Interface:</b>
1355 \code
1356 value_type v = vigra::NumericTraits<value_type>::one(); // if norm is not
1357 // given explicitly
1358 double d;
1359
1360 v = d * v;
1361 \endcode
1362 \deprecatedEnd
1363*/
1364
1365template <class ARITHTYPE = double>
1367{
1368 public:
1370
1371 /** the kernel's value type
1372 */
1373 typedef typename InternalVector::value_type value_type;
1374
1375 /** the kernel's reference type
1376 */
1377 typedef typename InternalVector::reference reference;
1378
1379 /** the kernel's const reference type
1380 */
1381 typedef typename InternalVector::const_reference const_reference;
1382
1383 /** deprecated -- use Kernel1D::iterator
1384 */
1385 typedef typename InternalVector::iterator Iterator;
1386
1387 /** 1D random access iterator over the kernel's values
1388 */
1389 typedef typename InternalVector::iterator iterator;
1390
1391 /** const 1D random access iterator over the kernel's values
1392 */
1393 typedef typename InternalVector::const_iterator const_iterator;
1394
1395 /** the kernel's accessor
1396 */
1398
1399 /** the kernel's const accessor
1400 */
1402
1403 struct InitProxy
1404 {
1405 InitProxy(Iterator i, int count, value_type & norm)
1406 : iter_(i), base_(i),
1407 count_(count), sum_(count),
1408 norm_(norm)
1409 {}
1410
1411 ~InitProxy()
1412 noexcept(false)
1413 {
1414 vigra_precondition(count_ == 1 || count_ == sum_,
1415 "Kernel1D::initExplicitly(): "
1416 "Wrong number of init values.");
1417 }
1418
1419 InitProxy & operator,(value_type const & v)
1420 {
1421 if(sum_ == count_)
1422 norm_ = *iter_;
1423
1424 norm_ += v;
1425
1426 --count_;
1427
1428 if(count_ > 0)
1429 {
1430 ++iter_;
1431 *iter_ = v;
1432 }
1433 return *this;
1434 }
1435
1436 Iterator iter_, base_;
1437 int count_, sum_;
1438 value_type & norm_;
1439 };
1440
1441 static value_type one() { return NumericTraits<value_type>::one(); }
1442
1443 /** Default constructor.
1444 Creates a kernel of size 1 which would copy the signal
1445 unchanged.
1446 */
1448 : kernel_(),
1449 left_(0),
1450 right_(0),
1451 border_treatment_(BORDER_TREATMENT_REFLECT),
1452 norm_(one())
1453 {
1454 kernel_.push_back(norm_);
1455 }
1456
1457 /** Copy constructor.
1458 */
1460 : kernel_(k.kernel_),
1461 left_(k.left_),
1462 right_(k.right_),
1463 border_treatment_(k.border_treatment_),
1464 norm_(k.norm_)
1465 {}
1466
1467 /** Construct from kernel with different element type, e.g. double => FixedPoint16.
1468 */
1469 template <class U>
1471 : kernel_(k.center()+k.left(), k.center()+k.right()+1),
1472 left_(k.left()),
1473 right_(k.right()),
1474 border_treatment_(k.borderTreatment()),
1475 norm_(k.norm())
1476 {}
1477
1478 /** Copy assignment.
1479 */
1481 {
1482 if(this != &k)
1483 {
1484 left_ = k.left_;
1485 right_ = k.right_;
1486 border_treatment_ = k.border_treatment_;
1487 norm_ = k.norm_;
1488 kernel_ = k.kernel_;
1489 }
1490 return *this;
1491 }
1492
1493 /** Initialization.
1494 This initializes the kernel with the given constant. The norm becomes
1495 v*size().
1496
1497 Instead of a single value an initializer list of length size()
1498 can be used like this:
1499
1500 \code
1501 vigra::Kernel1D<float> roberts_gradient_x;
1502
1503 roberts_gradient_x.initExplicitly(0, 1) = 1.0, -1.0;
1504 \endcode
1505
1506 In this case, the norm will be set to the sum of the init values.
1507 An initializer list of wrong length will result in a run-time error.
1508 */
1509 InitProxy operator=(value_type const & v)
1510 {
1511 int size = right_ - left_ + 1;
1512 for(unsigned int i=0; i<kernel_.size(); ++i) kernel_[i] = v;
1513 norm_ = (double)size*v;
1514
1515 return InitProxy(kernel_.begin(), size, norm_);
1516 }
1517
1518 /** Destructor.
1519 */
1521 {}
1522
1523 /**
1524 Init as a sampled Gaussian function. The radius of the kernel is
1525 always 3*std_dev. '<tt>norm</tt>' denotes the sum of all bins of the kernel
1526 (i.e. the kernel is corrected for the normalization error introduced
1527 by windowing the Gaussian to a finite interval). However,
1528 if <tt>norm</tt> is 0.0, the kernel is normalized to 1 by the analytic
1529 expression for the Gaussian, and <b>no</b> correction for the windowing
1530 error is performed. If <tt>windowRatio = 0.0</tt>, the radius of the filter
1531 window is <tt>radius = round(3.0 * std_dev)</tt>, otherwise it is
1532 <tt>radius = round(windowRatio * std_dev)</tt> (where <tt>windowRatio > 0.0</tt>
1533 is required).
1534
1535 Precondition:
1536 \code
1537 std_dev >= 0.0
1538 \endcode
1539
1540 Postconditions:
1541 \code
1542 1. left() == -(int)(3.0*std_dev + 0.5)
1543 2. right() == (int)(3.0*std_dev + 0.5)
1544 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1545 4. norm() == norm
1546 \endcode
1547 */
1548 void initGaussian(double std_dev, value_type norm, double windowRatio = 0.0);
1549
1550 /** Init as a Gaussian function with norm 1.
1551 */
1552 void initGaussian(double std_dev)
1553 {
1554 initGaussian(std_dev, one());
1555 }
1556
1557
1558 /**
1559 Init as Lindeberg's discrete analog of the Gaussian function. The radius of the kernel is
1560 always 3*std_dev. 'norm' denotes the sum of all bins of the kernel.
1561
1562 Precondition:
1563 \code
1564 std_dev >= 0.0
1565 \endcode
1566
1567 Postconditions:
1568 \code
1569 1. left() == -(int)(3.0*std_dev + 0.5)
1570 2. right() == (int)(3.0*std_dev + 0.5)
1571 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1572 4. norm() == norm
1573 \endcode
1574 */
1575 void initDiscreteGaussian(double std_dev, value_type norm);
1576
1577 /** Init as a Lindeberg's discrete analog of the Gaussian function
1578 with norm 1.
1579 */
1580 void initDiscreteGaussian(double std_dev)
1581 {
1582 initDiscreteGaussian(std_dev, one());
1583 }
1584
1585 /**
1586 Init as a Gaussian derivative of order '<tt>order</tt>'.
1587 The radius of the kernel is always <tt>3*std_dev + 0.5*order</tt>.
1588 '<tt>norm</tt>' denotes the norm of the kernel so that the
1589 following condition is fulfilled:
1590
1591 \f[ \sum_{i=left()}^{right()}
1592 \frac{(-i)^{order}kernel[i]}{order!} = norm
1593 \f]
1594
1595 Thus, the kernel will be corrected for the error introduced
1596 by windowing the Gaussian to a finite interval. However,
1597 if <tt>norm</tt> is 0.0, the kernel is normalized to 1 by the analytic
1598 expression for the Gaussian derivative, and <b>no</b> correction for the
1599 windowing error is performed. If <tt>windowRatio = 0.0</tt>, the radius
1600 of the filter window is <tt>radius = round(3.0 * std_dev + 0.5 * order)</tt>,
1601 otherwise it is <tt>radius = round(windowRatio * std_dev)</tt> (where
1602 <tt>windowRatio > 0.0</tt> is required).
1603
1604 Preconditions:
1605 \code
1606 1. std_dev >= 0.0
1607 2. order >= 1
1608 \endcode
1609
1610 Postconditions:
1611 \code
1612 1. left() == -(int)(3.0*std_dev + 0.5*order + 0.5)
1613 2. right() == (int)(3.0*std_dev + 0.5*order + 0.5)
1614 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1615 4. norm() == norm
1616 \endcode
1617 */
1618 void initGaussianDerivative(double std_dev, int order, value_type norm, double windowRatio = 0.0);
1619
1620 /** Init as a Gaussian derivative with norm 1.
1621 */
1622 void initGaussianDerivative(double std_dev, int order)
1623 {
1624 initGaussianDerivative(std_dev, order, one());
1625 }
1626
1627 /**
1628 Init an optimal 3-tap smoothing filter.
1629 The filter values are
1630
1631 \code
1632 [0.216, 0.568, 0.216]
1633 \endcode
1634
1635 These values are optimal in the sense that the 3x3 filter obtained by separable application
1636 of this filter is the best possible 3x3 approximation to a Gaussian filter.
1637 The equivalent Gaussian has sigma = 0.680.
1638
1639 Postconditions:
1640 \code
1641 1. left() == -1
1642 2. right() == 1
1643 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1644 4. norm() == 1.0
1645 \endcode
1646 */
1648 {
1649 this->initExplicitly(-1, 1) = 0.216, 0.568, 0.216;
1650 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1651 }
1652
1653 /**
1654 Init an optimal 3-tap smoothing filter to be used in the context of first derivative computation.
1655 This filter must be used in conjunction with the symmetric difference filter (see initSymmetricDifference()),
1656 such that the difference filter is applied along one dimension, and the smoothing filter along the other.
1657 The filter values are
1658
1659 \code
1660 [0.224365, 0.55127, 0.224365]
1661 \endcode
1662
1663 These values are optimal in the sense that the 3x3 filter obtained by combining
1664 this filter with the symmetric difference is the best possible 3x3 approximation to a
1665 Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.675.
1666
1667 Postconditions:
1668 \code
1669 1. left() == -1
1670 2. right() == 1
1671 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1672 4. norm() == 1.0
1673 \endcode
1674 */
1676 {
1677 this->initExplicitly(-1, 1) = 0.224365, 0.55127, 0.224365;
1678 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1679 }
1680
1681 /**
1682 Init an optimal 3-tap smoothing filter to be used in the context of second derivative computation.
1683 This filter must be used in conjunction with the 3-tap second difference filter (see initSecondDifference3()),
1684 such that the difference filter is applied along one dimension, and the smoothing filter along the other.
1685 The filter values are
1686
1687 \code
1688 [0.13, 0.74, 0.13]
1689 \endcode
1690
1691 These values are optimal in the sense that the 3x3 filter obtained by combining
1692 this filter with the 3-tap second difference is the best possible 3x3 approximation to a
1693 Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.433.
1694
1695 Postconditions:
1696 \code
1697 1. left() == -1
1698 2. right() == 1
1699 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1700 4. norm() == 1.0
1701 \endcode
1702 */
1704 {
1705 this->initExplicitly(-1, 1) = 0.13, 0.74, 0.13;
1706 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1707 }
1708
1709 /**
1710 Init an optimal 5-tap smoothing filter.
1711 The filter values are
1712
1713 \code
1714 [0.03134, 0.24, 0.45732, 0.24, 0.03134]
1715 \endcode
1716
1717 These values are optimal in the sense that the 5x5 filter obtained by separable application
1718 of this filter is the best possible 5x5 approximation to a Gaussian filter.
1719 The equivalent Gaussian has sigma = 0.867.
1720
1721 Postconditions:
1722 \code
1723 1. left() == -2
1724 2. right() == 2
1725 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1726 4. norm() == 1.0
1727 \endcode
1728 */
1730 {
1731 this->initExplicitly(-2, 2) = 0.03134, 0.24, 0.45732, 0.24, 0.03134;
1732 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1733 }
1734
1735 /**
1736 Init an optimal 5-tap smoothing filter to be used in the context of first derivative computation.
1737 This filter must be used in conjunction with the optimal 5-tap first derivative filter
1738 (see initOptimalFirstDerivative5()), such that the derivative filter is applied along one dimension,
1739 and the smoothing filter along the other. The filter values are
1740
1741 \code
1742 [0.04255, 0.241, 0.4329, 0.241, 0.04255]
1743 \endcode
1744
1745 These values are optimal in the sense that the 5x5 filter obtained by combining
1746 this filter with the optimal 5-tap first derivative is the best possible 5x5 approximation to a
1747 Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.906.
1748
1749 Postconditions:
1750 \code
1751 1. left() == -2
1752 2. right() == 2
1753 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1754 4. norm() == 1.0
1755 \endcode
1756 */
1758 {
1759 this->initExplicitly(-2, 2) = 0.04255, 0.241, 0.4329, 0.241, 0.04255;
1760 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1761 }
1762
1763 /**
1764 Init an optimal 5-tap smoothing filter to be used in the context of second derivative computation.
1765 This filter must be used in conjunction with the optimal 5-tap second derivative filter
1766 (see initOptimalSecondDerivative5()), such that the derivative filter is applied along one dimension,
1767 and the smoothing filter along the other. The filter values are
1768
1769 \code
1770 [0.0243, 0.23556, 0.48028, 0.23556, 0.0243]
1771 \endcode
1772
1773 These values are optimal in the sense that the 5x5 filter obtained by combining
1774 this filter with the optimal 5-tap second derivative is the best possible 5x5 approximation to a
1775 Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.817.
1776
1777 Postconditions:
1778 \code
1779 1. left() == -2
1780 2. right() == 2
1781 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1782 4. norm() == 1.0
1783 \endcode
1784 */
1786 {
1787 this->initExplicitly(-2, 2) = 0.0243, 0.23556, 0.48028, 0.23556, 0.0243;
1788 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1789 }
1790
1791 /**
1792 Init a 5-tap filter as defined by Peter Burt in the context of pyramid creation.
1793 The filter values are
1794
1795 \code
1796 [a, 0.25, 0.5-2*a, 0.25, a]
1797 \endcode
1798
1799 The default <tt>a = 0.04785</tt> is optimal in the sense that it minimizes the difference
1800 to a true Gaussian filter (which would have sigma = 0.975). For other values of <tt>a</tt>, the scale
1801 of the most similar Gaussian can be approximated by
1802
1803 \code
1804 sigma = 5.1 * a + 0.731
1805 \endcode
1806
1807 Preconditions:
1808 \code
1809 0 <= a <= 0.125
1810 \endcode
1811
1812 Postconditions:
1813 \code
1814 1. left() == -2
1815 2. right() == 2
1816 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1817 4. norm() == 1.0
1818 \endcode
1819 */
1820 void initBurtFilter(double a = 0.04785)
1821 {
1822 vigra_precondition(a >= 0.0 && a <= 0.125,
1823 "Kernel1D::initBurtFilter(): 0 <= a <= 0.125 required.");
1824 this->initExplicitly(-2, 2) = a, 0.25, 0.5 - 2.0*a, 0.25, a;
1825 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1826 }
1827
1828 /**
1829 Init as a Binomial filter. 'norm' denotes the sum of all bins
1830 of the kernel.
1831
1832 Precondition:
1833 \code
1834 radius >= 0
1835 \endcode
1836
1837 Postconditions:
1838 \code
1839 1. left() == -radius
1840 2. right() == radius
1841 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1842 4. norm() == norm
1843 \endcode
1844 */
1845 void initBinomial(int radius, value_type norm);
1846
1847 /** Init as a Binomial filter with norm 1.
1848 */
1849 void initBinomial(int radius)
1850 {
1851 initBinomial(radius, one());
1852 }
1853
1854 /**
1855 Init as an Averaging filter. 'norm' denotes the sum of all bins
1856 of the kernel. The window size is (2*radius+1) * (2*radius+1)
1857
1858 Precondition:
1859 \code
1860 radius >= 0
1861 \endcode
1862
1863 Postconditions:
1864 \code
1865 1. left() == -radius
1866 2. right() == radius
1867 3. borderTreatment() == BORDER_TREATMENT_CLIP
1868 4. norm() == norm
1869 \endcode
1870 */
1871 void initAveraging(int radius, value_type norm);
1872
1873 /** Init as an Averaging filter with norm 1.
1874 */
1875 void initAveraging(int radius)
1876 {
1877 initAveraging(radius, one());
1878 }
1879
1880 /**
1881 Init as a symmetric gradient filter of the form
1882 <TT>[ 0.5 * norm, 0.0 * norm, -0.5 * norm]</TT>
1883
1884 <b>Deprecated</b>. Use initSymmetricDifference() instead.
1885
1886 Postconditions:
1887 \code
1888 1. left() == -1
1889 2. right() == 1
1890 3. borderTreatment() == BORDER_TREATMENT_REPEAT
1891 4. norm() == norm
1892 \endcode
1893 */
1895 {
1897 setBorderTreatment(BORDER_TREATMENT_REPEAT);
1898 }
1899
1900 /** Init as a symmetric gradient filter with norm 1.
1901
1902 <b>Deprecated</b>. Use initSymmetricDifference() instead.
1903 */
1905 {
1906 initSymmetricGradient(one());
1907 }
1908
1909 /** Init as the 2-tap forward difference filter.
1910 The filter values are
1911
1912 \code
1913 [1.0, -1.0]
1914 \endcode
1915
1916 (note that filters are reflected by the convolution algorithm,
1917 and we get a forward difference after reflection).
1918
1919 Postconditions:
1920 \code
1921 1. left() == -1
1922 2. right() == 0
1923 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1924 4. norm() == 1.0
1925 \endcode
1926 */
1928 {
1929 this->initExplicitly(-1, 0) = 1.0, -1.0;
1930 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1931 }
1932
1933 /** Init as the 2-tap backward difference filter.
1934 The filter values are
1935
1936 \code
1937 [1.0, -1.0]
1938 \endcode
1939
1940 (note that filters are reflected by the convolution algorithm,
1941 and we get a forward difference after reflection).
1942
1943 Postconditions:
1944 \code
1945 1. left() == 0
1946 2. right() == 1
1947 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1948 4. norm() == 1.0
1949 \endcode
1950 */
1952 {
1953 this->initExplicitly(0, 1) = 1.0, -1.0;
1954 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1955 }
1956
1958
1959 /** Init as the 3-tap symmetric difference filter
1960 The filter values are
1961
1962 \code
1963 [0.5, 0, -0.5]
1964 \endcode
1965
1966 Postconditions:
1967 \code
1968 1. left() == -1
1969 2. right() == 1
1970 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1971 4. norm() == 1.0
1972 \endcode
1973 */
1975 {
1977 }
1978
1979 /**
1980 Init the 3-tap second difference filter.
1981 The filter values are
1982
1983 \code
1984 [1, -2, 1]
1985 \endcode
1986
1987 Postconditions:
1988 \code
1989 1. left() == -1
1990 2. right() == 1
1991 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1992 4. norm() == 1
1993 \endcode
1994 */
1996 {
1997 this->initExplicitly(-1, 1) = 1.0, -2.0, 1.0;
1998 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1999 }
2000
2001 /**
2002 Init an optimal 5-tap first derivative filter.
2003 This filter must be used in conjunction with the corresponding 5-tap smoothing filter
2004 (see initOptimalFirstDerivativeSmoothing5()), such that the derivative filter is applied along one dimension,
2005 and the smoothing filter along the other.
2006 The filter values are
2007
2008 \code
2009 [0.1, 0.3, 0.0, -0.3, -0.1]
2010 \endcode
2011
2012 These values are optimal in the sense that the 5x5 filter obtained by combining
2013 this filter with the corresponding 5-tap smoothing filter is the best possible 5x5 approximation to a
2014 Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.906.
2015
2016 If the filter is instead separably combined with itself, an almost optimal approximation of the
2017 mixed second Gaussian derivative at scale sigma = 0.899 results.
2018
2019 Postconditions:
2020 \code
2021 1. left() == -2
2022 2. right() == 2
2023 3. borderTreatment() == BORDER_TREATMENT_REFLECT
2024 4. norm() == 1.0
2025 \endcode
2026 */
2028 {
2029 this->initExplicitly(-2, 2) = 0.1, 0.3, 0.0, -0.3, -0.1;
2030 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
2031 }
2032
2033 /**
2034 Init an optimal 5-tap second derivative filter.
2035 This filter must be used in conjunction with the corresponding 5-tap smoothing filter
2036 (see initOptimalSecondDerivativeSmoothing5()), such that the derivative filter is applied along one dimension,
2037 and the smoothing filter along the other.
2038 The filter values are
2039
2040 \code
2041 [0.22075, 0.117, -0.6755, 0.117, 0.22075]
2042 \endcode
2043
2044 These values are optimal in the sense that the 5x5 filter obtained by combining
2045 this filter with the corresponding 5-tap smoothing filter is the best possible 5x5 approximation to a
2046 Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.817.
2047
2048 Postconditions:
2049 \code
2050 1. left() == -2
2051 2. right() == 2
2052 3. borderTreatment() == BORDER_TREATMENT_REFLECT
2053 4. norm() == 1.0
2054 \endcode
2055 */
2057 {
2058 this->initExplicitly(-2, 2) = 0.22075, 0.117, -0.6755, 0.117, 0.22075;
2059 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
2060 }
2061
2062 /** Init the kernel by an explicit initializer list.
2063 The left and right boundaries of the kernel must be passed.
2064 A comma-separated initializer list is given after the assignment
2065 operator. This function is used like this:
2066
2067 \code
2068 // define horizontal Roberts filter
2069 vigra::Kernel1D<float> roberts_gradient_x;
2070
2071 roberts_gradient_x.initExplicitly(0, 1) = 1.0, -1.0;
2072 \endcode
2073
2074 The norm is set to the sum of the initializer values. If the wrong number of
2075 values is given, a run-time error results. It is, however, possible to give
2076 just one initializer. This creates an averaging filter with the given constant:
2077
2078 \code
2079 vigra::Kernel1D<float> average5x1;
2080
2081 average5x1.initExplicitly(-2, 2) = 1.0/5.0;
2082 \endcode
2083
2084 Here, the norm is set to value*size().
2085
2086 <b> Preconditions:</b>
2087
2088 \code
2089
2090 1. left <= 0
2091 2. right >= 0
2092 3. the number of values in the initializer list
2093 is 1 or equals the size of the kernel.
2094 \endcode
2095 */
2097 {
2098 vigra_precondition(left <= 0,
2099 "Kernel1D::initExplicitly(): left border must be <= 0.");
2100 vigra_precondition(right >= 0,
2101 "Kernel1D::initExplicitly(): right border must be >= 0.");
2102
2103 right_ = right;
2104 left_ = left;
2105
2106 kernel_.resize(right - left + 1);
2107
2108 return *this;
2109 }
2110
2111 /** Get iterator to center of kernel
2112
2113 Postconditions:
2114 \code
2115
2116 center()[left()] ... center()[right()] are valid kernel positions
2117 \endcode
2118 */
2120 {
2121 return kernel_.begin() - left();
2122 }
2123
2124 const_iterator center() const
2125 {
2126 return kernel_.begin() - left();
2127 }
2128
2129 /** Access kernel value at specified location.
2130
2131 Preconditions:
2132 \code
2133
2134 left() <= location <= right()
2135 \endcode
2136 */
2137 reference operator[](int location)
2138 {
2139 return kernel_[location - left()];
2140 }
2141
2142 const_reference operator[](int location) const
2143 {
2144 return kernel_[location - left()];
2145 }
2146
2147 /** left border of kernel (inclusive), always <= 0
2148 */
2149 int left() const { return left_; }
2150
2151 /** right border of kernel (inclusive), always >= 0
2152 */
2153 int right() const { return right_; }
2154
2155 /** size of kernel (right() - left() + 1)
2156 */
2157 int size() const { return right_ - left_ + 1; }
2158
2159 /** current border treatment mode
2160 */
2161 BorderTreatmentMode borderTreatment() const
2162 { return border_treatment_; }
2163
2164 /** Set border treatment mode.
2165 */
2166 void setBorderTreatment( BorderTreatmentMode new_mode)
2167 { border_treatment_ = new_mode; }
2168
2169 /** norm of kernel
2170 */
2171 value_type norm() const { return norm_; }
2172
2173 /** set a new norm and normalize kernel, use the normalization formula
2174 for the given <tt>derivativeOrder</tt>.
2175 */
2176 void
2177 normalize(value_type norm, unsigned int derivativeOrder = 0, double offset = 0.0);
2178
2179 /** normalize kernel to norm 1.
2180 */
2181 void
2183 {
2184 normalize(one());
2185 }
2186
2187 /** get a const accessor
2188 */
2190
2191 /** get an accessor
2192 */
2194
2195 private:
2196 InternalVector kernel_;
2197 int left_, right_;
2198 BorderTreatmentMode border_treatment_;
2199 value_type norm_;
2200};
2201
2202template <class ARITHTYPE>
2204 unsigned int derivativeOrder,
2205 double offset)
2206{
2207 typedef typename NumericTraits<value_type>::RealPromote TmpType;
2208
2209 // find kernel sum
2210 Iterator k = kernel_.begin();
2211 TmpType sum = NumericTraits<TmpType>::zero();
2212
2213 if(derivativeOrder == 0)
2214 {
2215 for(; k < kernel_.end(); ++k)
2216 {
2217 sum += *k;
2218 }
2219 }
2220 else
2221 {
2222 unsigned int faculty = 1;
2223 for(unsigned int i = 2; i <= derivativeOrder; ++i)
2224 faculty *= i;
2225 for(double x = left() + offset; k < kernel_.end(); ++x, ++k)
2226 {
2227 sum = TmpType(sum + *k * VIGRA_CSTD::pow(-x, int(derivativeOrder)) / faculty);
2228 }
2229 }
2230
2231 vigra_precondition(sum != NumericTraits<value_type>::zero(),
2232 "Kernel1D<ARITHTYPE>::normalize(): "
2233 "Cannot normalize a kernel with sum = 0");
2234 // normalize
2235 sum = norm / sum;
2236 k = kernel_.begin();
2237 for(; k != kernel_.end(); ++k)
2238 {
2239 *k = *k * sum;
2240 }
2241
2242 norm_ = norm;
2243}
2244
2245/***********************************************************************/
2246
2247template <class ARITHTYPE>
2248void
2251 double windowRatio)
2252{
2253 vigra_precondition(std_dev >= 0.0,
2254 "Kernel1D::initGaussian(): Standard deviation must be >= 0.");
2255 vigra_precondition(windowRatio >= 0.0,
2256 "Kernel1D::initGaussian(): windowRatio must be >= 0.");
2257
2258 if(std_dev > 0.0)
2259 {
2260 Gaussian<ARITHTYPE> gauss((ARITHTYPE)std_dev);
2261
2262 // first calculate required kernel sizes
2263 int radius;
2264 if (windowRatio == 0.0)
2265 radius = (int)(3.0 * std_dev + 0.5);
2266 else
2267 radius = (int)(windowRatio * std_dev + 0.5);
2268 if(radius == 0)
2269 radius = 1;
2270
2271 // allocate the kernel
2272 kernel_.erase(kernel_.begin(), kernel_.end());
2273 kernel_.reserve(radius*2+1);
2274
2275 for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
2276 {
2277 kernel_.push_back(gauss(x));
2278 }
2279 left_ = -radius;
2280 right_ = radius;
2281 }
2282 else
2283 {
2284 kernel_.erase(kernel_.begin(), kernel_.end());
2285 kernel_.push_back(1.0);
2286 left_ = 0;
2287 right_ = 0;
2288 }
2289
2290 if(norm != 0.0)
2291 normalize(norm);
2292 else
2293 norm_ = 1.0;
2294
2295 // best border treatment for Gaussians is BORDER_TREATMENT_REFLECT
2296 border_treatment_ = BORDER_TREATMENT_REFLECT;
2297}
2298
2299/***********************************************************************/
2300
2301template <class ARITHTYPE>
2302void
2305{
2306 vigra_precondition(std_dev >= 0.0,
2307 "Kernel1D::initDiscreteGaussian(): Standard deviation must be >= 0.");
2308
2309 if(std_dev > 0.0)
2310 {
2311 // first calculate required kernel sizes
2312 int radius = (int)(3.0*std_dev + 0.5);
2313 if(radius == 0)
2314 radius = 1;
2315
2316 double f = 2.0 / std_dev / std_dev;
2317
2318 // allocate the working array
2319 int maxIndex = (int)(2.0 * (radius + 5.0 * VIGRA_CSTD::sqrt((double)radius)) + 0.5);
2320 InternalVector warray(maxIndex+1);
2321 warray[maxIndex] = 0.0;
2322 warray[maxIndex-1] = 1.0;
2323
2324 for(int i = maxIndex-2; i >= radius; --i)
2325 {
2326 warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
2327 if(warray[i] > 1.0e40)
2328 {
2329 warray[i+1] /= warray[i];
2330 warray[i] = 1.0;
2331 }
2332 }
2333
2334 // the following rescaling ensures that the numbers stay in a sensible range
2335 // during the rest of the iteration, so no other rescaling is needed
2336 double er = VIGRA_CSTD::exp(-radius*radius / (2.0*std_dev*std_dev));
2337 warray[radius+1] = er * warray[radius+1] / warray[radius];
2338 warray[radius] = er;
2339
2340 for(int i = radius-1; i >= 0; --i)
2341 {
2342 warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
2343 er += warray[i];
2344 }
2345
2346 double scale = norm / (2*er - warray[0]);
2347
2348 initExplicitly(-radius, radius);
2349 iterator c = center();
2350
2351 for(int i=0; i<=radius; ++i)
2352 {
2353 c[i] = c[-i] = warray[i] * scale;
2354 }
2355 }
2356 else
2357 {
2358 kernel_.erase(kernel_.begin(), kernel_.end());
2359 kernel_.push_back(norm);
2360 left_ = 0;
2361 right_ = 0;
2362 }
2363
2364 norm_ = norm;
2365
2366 // best border treatment for Gaussians is BORDER_TREATMENT_REFLECT
2367 border_treatment_ = BORDER_TREATMENT_REFLECT;
2368}
2369
2370/***********************************************************************/
2371
2372template <class ARITHTYPE>
2373void
2375 int order,
2377 double windowRatio)
2378{
2379 vigra_precondition(order >= 0,
2380 "Kernel1D::initGaussianDerivative(): Order must be >= 0.");
2381
2382 if(order == 0)
2383 {
2384 initGaussian(std_dev, norm, windowRatio);
2385 return;
2386 }
2387
2388 vigra_precondition(std_dev > 0.0,
2389 "Kernel1D::initGaussianDerivative(): "
2390 "Standard deviation must be > 0.");
2391 vigra_precondition(windowRatio >= 0.0,
2392 "Kernel1D::initGaussianDerivative(): windowRatio must be >= 0.");
2393
2394 Gaussian<ARITHTYPE> gauss((ARITHTYPE)std_dev, order);
2395
2396 // first calculate required kernel sizes
2397 int radius;
2398 if(windowRatio == 0.0)
2399 radius = (int)((3.0 + 0.5 * order) * std_dev + 0.5);
2400 else
2401 radius = (int)(windowRatio * std_dev + 0.5);
2402 if(radius == 0)
2403 radius = 1;
2404
2405 // allocate the kernels
2406 kernel_.clear();
2407 kernel_.reserve(radius*2+1);
2408
2409 // fill the kernel and calculate the DC component
2410 // introduced by truncation of the Gaussian
2411 ARITHTYPE dc = 0.0;
2412 for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
2413 {
2414 kernel_.push_back(gauss(x));
2415 dc += kernel_[kernel_.size()-1];
2416 }
2417 dc = ARITHTYPE(dc / (2.0*radius + 1.0));
2418
2419 // remove DC, but only if kernel correction is permitted by a non-zero
2420 // value for norm
2421 if(norm != 0.0)
2422 {
2423 for(unsigned int i=0; i < kernel_.size(); ++i)
2424 {
2425 kernel_[i] -= dc;
2426 }
2427 }
2428
2429 left_ = -radius;
2430 right_ = radius;
2431
2432 if(norm != 0.0)
2433 normalize(norm, order);
2434 else
2435 norm_ = 1.0;
2436
2437 // best border treatment for Gaussian derivatives is
2438 // BORDER_TREATMENT_REFLECT
2439 border_treatment_ = BORDER_TREATMENT_REFLECT;
2440}
2441
2442/***********************************************************************/
2443
2444template <class ARITHTYPE>
2445void
2448{
2449 vigra_precondition(radius > 0,
2450 "Kernel1D::initBinomial(): Radius must be > 0.");
2451
2452 // allocate the kernel
2453 InternalVector(radius*2+1).swap(kernel_);
2454 typename InternalVector::iterator x = kernel_.begin() + radius;
2455
2456 // fill kernel
2457 x[radius] = norm;
2458 for(int j=radius-1; j>=-radius; --j)
2459 {
2460 x[j] = 0.5 * x[j+1];
2461 for(int i=j+1; i<radius; ++i)
2462 {
2463 x[i] = 0.5 * (x[i] + x[i+1]);
2464 }
2465 x[radius] *= 0.5;
2466 }
2467
2468 left_ = -radius;
2469 right_ = radius;
2470 norm_ = norm;
2471
2472 // best border treatment for Binomial is BORDER_TREATMENT_REFLECT
2473 border_treatment_ = BORDER_TREATMENT_REFLECT;
2474}
2475
2476/***********************************************************************/
2477
2478template <class ARITHTYPE>
2479void
2482{
2483 vigra_precondition(radius > 0,
2484 "Kernel1D::initAveraging(): Radius must be > 0.");
2485
2486 // calculate scaling
2487 double scale = 1.0 / (radius * 2 + 1);
2488
2489 // normalize
2490 kernel_.erase(kernel_.begin(), kernel_.end());
2491 kernel_.reserve(radius*2+1);
2492
2493 for(int i=0; i<=radius*2+1; ++i)
2494 {
2495 kernel_.push_back(scale * norm);
2496 }
2497
2498 left_ = -radius;
2499 right_ = radius;
2500 norm_ = norm;
2501
2502 // best border treatment for Averaging is BORDER_TREATMENT_CLIP
2503 border_treatment_ = BORDER_TREATMENT_CLIP;
2504}
2505
2506/***********************************************************************/
2507
2508template <class ARITHTYPE>
2509void
2511{
2512 kernel_.erase(kernel_.begin(), kernel_.end());
2513 kernel_.reserve(3);
2514
2515 kernel_.push_back(ARITHTYPE(0.5 * norm));
2516 kernel_.push_back(ARITHTYPE(0.0 * norm));
2517 kernel_.push_back(ARITHTYPE(-0.5 * norm));
2518
2519 left_ = -1;
2520 right_ = 1;
2521 norm_ = norm;
2522
2523 // best border treatment for symmetric difference is
2524 // BORDER_TREATMENT_REFLECT
2525 border_treatment_ = BORDER_TREATMENT_REFLECT;
2526}
2527
2528/**************************************************************/
2529/* */
2530/* Argument object factories for Kernel1D */
2531/* */
2532/* (documentation: see vigra/convolution.hxx) */
2533/* */
2534/**************************************************************/
2535
2536template <class KernelIterator, class KernelAccessor>
2537inline
2538tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>
2539kernel1d(KernelIterator ik, KernelAccessor ka,
2540 int kleft, int kright, BorderTreatmentMode border)
2541{
2542 return
2543 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>(
2544 ik, ka, kleft, kright, border);
2545}
2546
2547template <class T>
2548inline
2549tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
2550 int, int, BorderTreatmentMode>
2551kernel1d(Kernel1D<T> const & k)
2552
2553{
2554 return
2555 tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
2556 int, int, BorderTreatmentMode>(
2557 k.center(),
2558 k.accessor(),
2559 k.left(), k.right(),
2560 k.borderTreatment());
2561}
2562
2563template <class T>
2564inline
2565tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
2566 int, int, BorderTreatmentMode>
2567kernel1d(Kernel1D<T> const & k, BorderTreatmentMode border)
2568
2569{
2570 return
2571 tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
2572 int, int, BorderTreatmentMode>(
2573 k.center(),
2574 k.accessor(),
2575 k.left(), k.right(),
2576 border);
2577}
2578
2579
2580} // namespace vigra
2581
2582#endif // VIGRA_SEPARABLECONVOLUTION_HXX
const_iterator begin() const
Definition: array_vector.hxx:223
size_type size() const
Definition: array_vector.hxx:358
Definition: array_vector.hxx:514
Definition: gaussians.hxx:64
Generic 1 dimensional convolution kernel.
Definition: separableconvolution.hxx:1367
void initBinomial(int radius)
Definition: separableconvolution.hxx:1849
void initOptimalFirstDerivativeSmoothing5()
Definition: separableconvolution.hxx:1757
void initSecondDifference3()
Definition: separableconvolution.hxx:1995
void initOptimalSecondDerivativeSmoothing3()
Definition: separableconvolution.hxx:1703
InternalVector::reference reference
Definition: separableconvolution.hxx:1377
void initBurtFilter(double a=0.04785)
Definition: separableconvolution.hxx:1820
void initDiscreteGaussian(double std_dev)
Definition: separableconvolution.hxx:1580
void initBackwardDifference()
Definition: separableconvolution.hxx:1951
void initOptimalFirstDerivative5()
Definition: separableconvolution.hxx:2027
int right() const
Definition: separableconvolution.hxx:2153
value_type norm() const
Definition: separableconvolution.hxx:2171
reference operator[](int location)
Definition: separableconvolution.hxx:2137
Kernel1D(Kernel1D< U > const &k)
Definition: separableconvolution.hxx:1470
void initGaussian(double std_dev, value_type norm, double windowRatio=0.0)
Definition: separableconvolution.hxx:2249
InitProxy operator=(value_type const &v)
Definition: separableconvolution.hxx:1509
void initSymmetricGradient()
Definition: separableconvolution.hxx:1904
BorderTreatmentMode borderTreatment() const
Definition: separableconvolution.hxx:2161
void setBorderTreatment(BorderTreatmentMode new_mode)
Definition: separableconvolution.hxx:2166
StandardAccessor< ARITHTYPE > Accessor
Definition: separableconvolution.hxx:1397
void initOptimalSmoothing5()
Definition: separableconvolution.hxx:1729
ConstAccessor accessor() const
Definition: separableconvolution.hxx:2189
void initGaussianDerivative(double std_dev, int order)
Definition: separableconvolution.hxx:1622
void initDiscreteGaussian(double std_dev, value_type norm)
Definition: separableconvolution.hxx:2303
InternalVector::value_type value_type
Definition: separableconvolution.hxx:1373
void initGaussianDerivative(double std_dev, int order, value_type norm, double windowRatio=0.0)
Definition: separableconvolution.hxx:2374
InternalVector::iterator iterator
Definition: separableconvolution.hxx:1389
~Kernel1D()
Definition: separableconvolution.hxx:1520
void initSymmetricDifference()
Definition: separableconvolution.hxx:1974
void initAveraging(int radius)
Definition: separableconvolution.hxx:1875
StandardConstAccessor< ARITHTYPE > ConstAccessor
Definition: separableconvolution.hxx:1401
Kernel1D(Kernel1D const &k)
Definition: separableconvolution.hxx:1459
void initOptimalSecondDerivative5()
Definition: separableconvolution.hxx:2056
InternalVector::const_iterator const_iterator
Definition: separableconvolution.hxx:1393
InternalVector::const_reference const_reference
Definition: separableconvolution.hxx:1381
void initAveraging(int radius, value_type norm)
Definition: separableconvolution.hxx:2480
void initSymmetricGradient(value_type norm)
Definition: separableconvolution.hxx:1894
void initGaussian(double std_dev)
Definition: separableconvolution.hxx:1552
void initOptimalSecondDerivativeSmoothing5()
Definition: separableconvolution.hxx:1785
int left() const
Definition: separableconvolution.hxx:2149
Accessor accessor()
Definition: separableconvolution.hxx:2193
Kernel1D()
Definition: separableconvolution.hxx:1447
void initBinomial(int radius, value_type norm)
Definition: separableconvolution.hxx:2446
void normalize()
Definition: separableconvolution.hxx:2182
void initForwardDifference()
Definition: separableconvolution.hxx:1927
InternalVector::iterator Iterator
Definition: separableconvolution.hxx:1385
Kernel1D & initExplicitly(int left, int right)
Definition: separableconvolution.hxx:2096
void initOptimalSmoothing3()
Definition: separableconvolution.hxx:1647
void initOptimalFirstDerivativeSmoothing3()
Definition: separableconvolution.hxx:1675
int size() const
Definition: separableconvolution.hxx:2157
iterator center()
Definition: separableconvolution.hxx:2119
Kernel1D & operator=(Kernel1D const &k)
Definition: separableconvolution.hxx:1480
Encapsulate access to the values an iterator points to.
Definition: accessor.hxx:134
Encapsulate read access to the values an iterator points to.
Definition: accessor.hxx:270
void convolveLine(...)
Performs a 1-dimensional convolution of the source signal using the given kernel.
void separableConvolveX(...)
Performs a 1 dimensional convolution in x direction.
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
void separableConvolveY(...)
Performs a 1 dimensional convolution in y direction.
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:2073
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.11.1