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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.12.4