[ 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 // 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
222 template <class SrcIterator, class SrcAccessor,
223 class DestIterator, class DestAccessor,
224 class KernelIterator, class KernelAccessor>
225 void 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
322 template <class SrcIterator, class SrcAccessor,
323 class DestIterator, class DestAccessor,
324 class KernelIterator, class KernelAccessor,
325 class Norm>
326 void 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
427 template <class SrcIterator, class SrcAccessor,
428 class DestIterator, class DestAccessor,
429 class KernelIterator, class KernelAccessor>
430 void 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
506 template <class SrcIterator, class SrcAccessor,
507 class DestIterator, class DestAccessor,
508 class KernelIterator, class KernelAccessor>
509 void 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
605 template <class SrcIterator, class SrcAccessor,
606 class DestIterator, class DestAccessor,
607 class KernelIterator, class KernelAccessor>
608 void 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
704 template <class SrcIterator, class SrcAccessor,
705 class DestIterator, class DestAccessor,
706 class KernelIterator, class KernelAccessor>
707 void 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*/
894doxygen_overloaded_function(template <...> void convolveLine)
895
896 template <class SrcIterator, class SrcAccessor,
897 class DestIterator, class DestAccessor,
898 class KernelIterator, class KernelAccessor>
899 void 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
975 template <class SrcIterator, class SrcAccessor,
976 class DestIterator, class DestAccessor,
977 class KernelIterator, class KernelAccessor>
978 inline
979 void 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*/
1083doxygen_overloaded_function(template <...> void separableConvolveX)
1084
1085 template <class SrcIterator, class SrcAccessor,
1086 class DestIterator, class DestAccessor,
1087 class KernelIterator, class KernelAccessor>
1088 void 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
1117 template <class SrcIterator, class SrcAccessor,
1118 class DestIterator, class DestAccessor,
1119 class KernelIterator, class KernelAccessor>
1120 inline void
1121 separableConvolveX(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
1132 template <class T1, class S1,
1133 class T2, class S2,
1134 class T3>
1135 inline void
1136 separableConvolveX(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*/
1238doxygen_overloaded_function(template <...> void separableConvolveY)
1239
1240 template <class SrcIterator, class SrcAccessor,
1241 class DestIterator, class DestAccessor,
1242 class KernelIterator, class KernelAccessor>
1243 void 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
1272 template <class SrcIterator, class SrcAccessor,
1273 class DestIterator, class DestAccessor,
1274 class KernelIterator, class KernelAccessor>
1275 inline void
1276 separableConvolveY(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
1287 template <class T1, class S1,
1288 class T2, class S2,
1289 class T3>
1290 inline void
1291 separableConvolveY(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
1365 template <class ARITHTYPE = double>
1366 class Kernel1D
1367{
1368 public:
1369 typedef ArrayVector<ARITHTYPE> InternalVector;
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 */
1397 typedef StandardAccessor<ARITHTYPE> Accessor;
1398
1399 /** the kernel's const accessor
1400 */
1401 typedef StandardConstAccessor<ARITHTYPE> ConstAccessor;
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#if _MSC_VER >= 1900 || __cplusplus >= 201103L
1413 noexcept(false)
1414#else
1415 throw(PreconditionViolation)
1416#endif
1417 {
1418 vigra_precondition(count_ == 1 || count_ == sum_,
1419 "Kernel1D::initExplicitly(): "
1420 "Wrong number of init values.");
1421 }
1422
1423 InitProxy & operator,(value_type const & v)
1424 {
1425 if(sum_ == count_)
1426 norm_ = *iter_;
1427
1428 norm_ += v;
1429
1430 --count_;
1431
1432 if(count_ > 0)
1433 {
1434 ++iter_;
1435 *iter_ = v;
1436 }
1437 return *this;
1438 }
1439
1440 Iterator iter_, base_;
1441 int count_, sum_;
1442 value_type & norm_;
1443 };
1444
1445 static value_type one() { return NumericTraits<value_type>::one(); }
1446
1447 /** Default constructor.
1448 Creates a kernel of size 1 which would copy the signal
1449 unchanged.
1450 */
1451 Kernel1D()
1452 : kernel_(),
1453 left_(0),
1454 right_(0),
1455 border_treatment_(BORDER_TREATMENT_REFLECT),
1456 norm_(one())
1457 {
1458 kernel_.push_back(norm_);
1459 }
1460
1461 /** Copy constructor.
1462 */
1463 Kernel1D(Kernel1D const & k)
1464 : kernel_(k.kernel_),
1465 left_(k.left_),
1466 right_(k.right_),
1467 border_treatment_(k.border_treatment_),
1468 norm_(k.norm_)
1469 {}
1470
1471 /** Construct from kernel with different element type, e.g. double => FixedPoint16.
1472 */
1473 template <class U>
1474 Kernel1D(Kernel1D<U> const & k)
1475 : kernel_(k.center()+k.left(), k.center()+k.right()+1),
1476 left_(k.left()),
1477 right_(k.right()),
1478 border_treatment_(k.borderTreatment()),
1479 norm_(k.norm())
1480 {}
1481
1482 /** Copy assignment.
1483 */
1484 Kernel1D & operator=(Kernel1D const & k)
1485 {
1486 if(this != &k)
1487 {
1488 left_ = k.left_;
1489 right_ = k.right_;
1490 border_treatment_ = k.border_treatment_;
1491 norm_ = k.norm_;
1492 kernel_ = k.kernel_;
1493 }
1494 return *this;
1495 }
1496
1497 /** Initialization.
1498 This initializes the kernel with the given constant. The norm becomes
1499 v*size().
1500
1501 Instead of a single value an initializer list of length size()
1502 can be used like this:
1503
1504 \code
1505 vigra::Kernel1D<float> roberts_gradient_x;
1506
1507 roberts_gradient_x.initExplicitly(0, 1) = 1.0, -1.0;
1508 \endcode
1509
1510 In this case, the norm will be set to the sum of the init values.
1511 An initializer list of wrong length will result in a run-time error.
1512 */
1513 InitProxy operator=(value_type const & v)
1514 {
1515 int size = right_ - left_ + 1;
1516 for(unsigned int i=0; i<kernel_.size(); ++i) kernel_[i] = v;
1517 norm_ = (double)size*v;
1518
1519 return InitProxy(kernel_.begin(), size, norm_);
1520 }
1521
1522 /** Destructor.
1523 */
1524 ~Kernel1D()
1525 {}
1526
1527 /**
1528 Init as a sampled Gaussian function. The radius of the kernel is
1529 always 3*std_dev. '<tt>norm</tt>' denotes the sum of all bins of the kernel
1530 (i.e. the kernel is corrected for the normalization error introduced
1531 by windowing the Gaussian to a finite interval). However,
1532 if <tt>norm</tt> is 0.0, the kernel is normalized to 1 by the analytic
1533 expression for the Gaussian, and <b>no</b> correction for the windowing
1534 error is performed. If <tt>windowRatio = 0.0</tt>, the radius of the filter
1535 window is <tt>radius = round(3.0 * std_dev)</tt>, otherwise it is
1536 <tt>radius = round(windowRatio * std_dev)</tt> (where <tt>windowRatio > 0.0</tt>
1537 is required).
1538
1539 Precondition:
1540 \code
1541 std_dev >= 0.0
1542 \endcode
1543
1544 Postconditions:
1545 \code
1546 1. left() == -(int)(3.0*std_dev + 0.5)
1547 2. right() == (int)(3.0*std_dev + 0.5)
1548 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1549 4. norm() == norm
1550 \endcode
1551 */
1552 void initGaussian(double std_dev, value_type norm, double windowRatio = 0.0);
1553
1554 /** Init as a Gaussian function with norm 1.
1555 */
1556 void initGaussian(double std_dev)
1557 {
1558 initGaussian(std_dev, one());
1559 }
1560
1561
1562 /**
1563 Init as Lindeberg's discrete analog of the Gaussian function. The radius of the kernel is
1564 always 3*std_dev. 'norm' denotes the sum of all bins of the kernel.
1565
1566 Precondition:
1567 \code
1568 std_dev >= 0.0
1569 \endcode
1570
1571 Postconditions:
1572 \code
1573 1. left() == -(int)(3.0*std_dev + 0.5)
1574 2. right() == (int)(3.0*std_dev + 0.5)
1575 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1576 4. norm() == norm
1577 \endcode
1578 */
1579 void initDiscreteGaussian(double std_dev, value_type norm);
1580
1581 /** Init as a Lindeberg's discrete analog of the Gaussian function
1582 with norm 1.
1583 */
1584 void initDiscreteGaussian(double std_dev)
1585 {
1586 initDiscreteGaussian(std_dev, one());
1587 }
1588
1589 /**
1590 Init as a Gaussian derivative of order '<tt>order</tt>'.
1591 The radius of the kernel is always <tt>3*std_dev + 0.5*order</tt>.
1592 '<tt>norm</tt>' denotes the norm of the kernel so that the
1593 following condition is fulfilled:
1594
1595 \f[ \sum_{i=left()}^{right()}
1596 \frac{(-i)^{order}kernel[i]}{order!} = norm
1597 \f]
1598
1599 Thus, the kernel will be corrected for the error introduced
1600 by windowing the Gaussian to a finite interval. However,
1601 if <tt>norm</tt> is 0.0, the kernel is normalized to 1 by the analytic
1602 expression for the Gaussian derivative, and <b>no</b> correction for the
1603 windowing error is performed. If <tt>windowRatio = 0.0</tt>, the radius
1604 of the filter window is <tt>radius = round(3.0 * std_dev + 0.5 * order)</tt>,
1605 otherwise it is <tt>radius = round(windowRatio * std_dev)</tt> (where
1606 <tt>windowRatio > 0.0</tt> is required).
1607
1608 Preconditions:
1609 \code
1610 1. std_dev >= 0.0
1611 2. order >= 1
1612 \endcode
1613
1614 Postconditions:
1615 \code
1616 1. left() == -(int)(3.0*std_dev + 0.5*order + 0.5)
1617 2. right() == (int)(3.0*std_dev + 0.5*order + 0.5)
1618 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1619 4. norm() == norm
1620 \endcode
1621 */
1622 void initGaussianDerivative(double std_dev, int order, value_type norm, double windowRatio = 0.0);
1623
1624 /** Init as a Gaussian derivative with norm 1.
1625 */
1626 void initGaussianDerivative(double std_dev, int order)
1627 {
1628 initGaussianDerivative(std_dev, order, one());
1629 }
1630
1631 /**
1632 Init an optimal 3-tap smoothing filter.
1633 The filter values are
1634
1635 \code
1636 [0.216, 0.568, 0.216]
1637 \endcode
1638
1639 These values are optimal in the sense that the 3x3 filter obtained by separable application
1640 of this filter is the best possible 3x3 approximation to a Gaussian filter.
1641 The equivalent Gaussian has sigma = 0.680.
1642
1643 Postconditions:
1644 \code
1645 1. left() == -1
1646 2. right() == 1
1647 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1648 4. norm() == 1.0
1649 \endcode
1650 */
1651 void initOptimalSmoothing3()
1652 {
1653 this->initExplicitly(-1, 1) = 0.216, 0.568, 0.216;
1654 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1655 }
1656
1657 /**
1658 Init an optimal 3-tap smoothing filter to be used in the context of first derivative computation.
1659 This filter must be used in conjunction with the symmetric difference filter (see initSymmetricDifference()),
1660 such that the difference filter is applied along one dimension, and the smoothing filter along the other.
1661 The filter values are
1662
1663 \code
1664 [0.224365, 0.55127, 0.224365]
1665 \endcode
1666
1667 These values are optimal in the sense that the 3x3 filter obtained by combining
1668 this filter with the symmetric difference is the best possible 3x3 approximation to a
1669 Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.675.
1670
1671 Postconditions:
1672 \code
1673 1. left() == -1
1674 2. right() == 1
1675 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1676 4. norm() == 1.0
1677 \endcode
1678 */
1679 void initOptimalFirstDerivativeSmoothing3()
1680 {
1681 this->initExplicitly(-1, 1) = 0.224365, 0.55127, 0.224365;
1682 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1683 }
1684
1685 /**
1686 Init an optimal 3-tap smoothing filter to be used in the context of second derivative computation.
1687 This filter must be used in conjunction with the 3-tap second difference filter (see initSecondDifference3()),
1688 such that the difference filter is applied along one dimension, and the smoothing filter along the other.
1689 The filter values are
1690
1691 \code
1692 [0.13, 0.74, 0.13]
1693 \endcode
1694
1695 These values are optimal in the sense that the 3x3 filter obtained by combining
1696 this filter with the 3-tap second difference is the best possible 3x3 approximation to a
1697 Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.433.
1698
1699 Postconditions:
1700 \code
1701 1. left() == -1
1702 2. right() == 1
1703 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1704 4. norm() == 1.0
1705 \endcode
1706 */
1707 void initOptimalSecondDerivativeSmoothing3()
1708 {
1709 this->initExplicitly(-1, 1) = 0.13, 0.74, 0.13;
1710 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1711 }
1712
1713 /**
1714 Init an optimal 5-tap smoothing filter.
1715 The filter values are
1716
1717 \code
1718 [0.03134, 0.24, 0.45732, 0.24, 0.03134]
1719 \endcode
1720
1721 These values are optimal in the sense that the 5x5 filter obtained by separable application
1722 of this filter is the best possible 5x5 approximation to a Gaussian filter.
1723 The equivalent Gaussian has sigma = 0.867.
1724
1725 Postconditions:
1726 \code
1727 1. left() == -2
1728 2. right() == 2
1729 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1730 4. norm() == 1.0
1731 \endcode
1732 */
1733 void initOptimalSmoothing5()
1734 {
1735 this->initExplicitly(-2, 2) = 0.03134, 0.24, 0.45732, 0.24, 0.03134;
1736 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1737 }
1738
1739 /**
1740 Init an optimal 5-tap smoothing filter to be used in the context of first derivative computation.
1741 This filter must be used in conjunction with the optimal 5-tap first derivative filter
1742 (see initOptimalFirstDerivative5()), such that the derivative filter is applied along one dimension,
1743 and the smoothing filter along the other. The filter values are
1744
1745 \code
1746 [0.04255, 0.241, 0.4329, 0.241, 0.04255]
1747 \endcode
1748
1749 These values are optimal in the sense that the 5x5 filter obtained by combining
1750 this filter with the optimal 5-tap first derivative is the best possible 5x5 approximation to a
1751 Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.906.
1752
1753 Postconditions:
1754 \code
1755 1. left() == -2
1756 2. right() == 2
1757 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1758 4. norm() == 1.0
1759 \endcode
1760 */
1761 void initOptimalFirstDerivativeSmoothing5()
1762 {
1763 this->initExplicitly(-2, 2) = 0.04255, 0.241, 0.4329, 0.241, 0.04255;
1764 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1765 }
1766
1767 /**
1768 Init an optimal 5-tap smoothing filter to be used in the context of second derivative computation.
1769 This filter must be used in conjunction with the optimal 5-tap second derivative filter
1770 (see initOptimalSecondDerivative5()), such that the derivative filter is applied along one dimension,
1771 and the smoothing filter along the other. The filter values are
1772
1773 \code
1774 [0.0243, 0.23556, 0.48028, 0.23556, 0.0243]
1775 \endcode
1776
1777 These values are optimal in the sense that the 5x5 filter obtained by combining
1778 this filter with the optimal 5-tap second derivative is the best possible 5x5 approximation to a
1779 Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.817.
1780
1781 Postconditions:
1782 \code
1783 1. left() == -2
1784 2. right() == 2
1785 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1786 4. norm() == 1.0
1787 \endcode
1788 */
1789 void initOptimalSecondDerivativeSmoothing5()
1790 {
1791 this->initExplicitly(-2, 2) = 0.0243, 0.23556, 0.48028, 0.23556, 0.0243;
1792 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1793 }
1794
1795 /**
1796 Init a 5-tap filter as defined by Peter Burt in the context of pyramid creation.
1797 The filter values are
1798
1799 \code
1800 [a, 0.25, 0.5-2*a, 0.25, a]
1801 \endcode
1802
1803 The default <tt>a = 0.04785</tt> is optimal in the sense that it minimizes the difference
1804 to a true Gaussian filter (which would have sigma = 0.975). For other values of <tt>a</tt>, the scale
1805 of the most similar Gaussian can be approximated by
1806
1807 \code
1808 sigma = 5.1 * a + 0.731
1809 \endcode
1810
1811 Preconditions:
1812 \code
1813 0 <= a <= 0.125
1814 \endcode
1815
1816 Postconditions:
1817 \code
1818 1. left() == -2
1819 2. right() == 2
1820 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1821 4. norm() == 1.0
1822 \endcode
1823 */
1824 void initBurtFilter(double a = 0.04785)
1825 {
1826 vigra_precondition(a >= 0.0 && a <= 0.125,
1827 "Kernel1D::initBurtFilter(): 0 <= a <= 0.125 required.");
1828 this->initExplicitly(-2, 2) = a, 0.25, 0.5 - 2.0*a, 0.25, a;
1829 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1830 }
1831
1832 /**
1833 Init as a Binomial filter. 'norm' denotes the sum of all bins
1834 of the kernel.
1835
1836 Precondition:
1837 \code
1838 radius >= 0
1839 \endcode
1840
1841 Postconditions:
1842 \code
1843 1. left() == -radius
1844 2. right() == radius
1845 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1846 4. norm() == norm
1847 \endcode
1848 */
1849 void initBinomial(int radius, value_type norm);
1850
1851 /** Init as a Binomial filter with norm 1.
1852 */
1853 void initBinomial(int radius)
1854 {
1855 initBinomial(radius, one());
1856 }
1857
1858 /**
1859 Init as an Averaging filter. 'norm' denotes the sum of all bins
1860 of the kernel. The window size is (2*radius+1) * (2*radius+1)
1861
1862 Precondition:
1863 \code
1864 radius >= 0
1865 \endcode
1866
1867 Postconditions:
1868 \code
1869 1. left() == -radius
1870 2. right() == radius
1871 3. borderTreatment() == BORDER_TREATMENT_CLIP
1872 4. norm() == norm
1873 \endcode
1874 */
1875 void initAveraging(int radius, value_type norm);
1876
1877 /** Init as an Averaging filter with norm 1.
1878 */
1879 void initAveraging(int radius)
1880 {
1881 initAveraging(radius, one());
1882 }
1883
1884 /**
1885 Init as a symmetric gradient filter of the form
1886 <TT>[ 0.5 * norm, 0.0 * norm, -0.5 * norm]</TT>
1887
1888 <b>Deprecated</b>. Use initSymmetricDifference() instead.
1889
1890 Postconditions:
1891 \code
1892 1. left() == -1
1893 2. right() == 1
1894 3. borderTreatment() == BORDER_TREATMENT_REPEAT
1895 4. norm() == norm
1896 \endcode
1897 */
1898 void initSymmetricGradient(value_type norm )
1899 {
1900 initSymmetricDifference(norm);
1901 setBorderTreatment(BORDER_TREATMENT_REPEAT);
1902 }
1903
1904 /** Init as a symmetric gradient filter with norm 1.
1905
1906 <b>Deprecated</b>. Use initSymmetricDifference() instead.
1907 */
1908 void initSymmetricGradient()
1909 {
1910 initSymmetricGradient(one());
1911 }
1912
1913 /** Init as the 2-tap forward difference filter.
1914 The filter values are
1915
1916 \code
1917 [1.0, -1.0]
1918 \endcode
1919
1920 (note that filters are reflected by the convolution algorithm,
1921 and we get a forward difference after reflection).
1922
1923 Postconditions:
1924 \code
1925 1. left() == -1
1926 2. right() == 0
1927 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1928 4. norm() == 1.0
1929 \endcode
1930 */
1931 void initForwardDifference()
1932 {
1933 this->initExplicitly(-1, 0) = 1.0, -1.0;
1934 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1935 }
1936
1937 /** Init as the 2-tap backward difference filter.
1938 The filter values are
1939
1940 \code
1941 [1.0, -1.0]
1942 \endcode
1943
1944 (note that filters are reflected by the convolution algorithm,
1945 and we get a forward difference after reflection).
1946
1947 Postconditions:
1948 \code
1949 1. left() == 0
1950 2. right() == 1
1951 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1952 4. norm() == 1.0
1953 \endcode
1954 */
1955 void initBackwardDifference()
1956 {
1957 this->initExplicitly(0, 1) = 1.0, -1.0;
1958 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1959 }
1960
1961 void initSymmetricDifference(value_type norm );
1962
1963 /** Init as the 3-tap symmetric difference filter
1964 The filter values are
1965
1966 \code
1967 [0.5, 0, -0.5]
1968 \endcode
1969
1970 Postconditions:
1971 \code
1972 1. left() == -1
1973 2. right() == 1
1974 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1975 4. norm() == 1.0
1976 \endcode
1977 */
1978 void initSymmetricDifference()
1979 {
1980 initSymmetricDifference(one());
1981 }
1982
1983 /**
1984 Init the 3-tap second difference filter.
1985 The filter values are
1986
1987 \code
1988 [1, -2, 1]
1989 \endcode
1990
1991 Postconditions:
1992 \code
1993 1. left() == -1
1994 2. right() == 1
1995 3. borderTreatment() == BORDER_TREATMENT_REFLECT
1996 4. norm() == 1
1997 \endcode
1998 */
1999 void initSecondDifference3()
2000 {
2001 this->initExplicitly(-1, 1) = 1.0, -2.0, 1.0;
2002 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
2003 }
2004
2005 /**
2006 Init an optimal 5-tap first derivative filter.
2007 This filter must be used in conjunction with the corresponding 5-tap smoothing filter
2008 (see initOptimalFirstDerivativeSmoothing5()), such that the derivative filter is applied along one dimension,
2009 and the smoothing filter along the other.
2010 The filter values are
2011
2012 \code
2013 [0.1, 0.3, 0.0, -0.3, -0.1]
2014 \endcode
2015
2016 These values are optimal in the sense that the 5x5 filter obtained by combining
2017 this filter with the corresponding 5-tap smoothing filter is the best possible 5x5 approximation to a
2018 Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.906.
2019
2020 If the filter is instead separably combined with itself, an almost optimal approximation of the
2021 mixed second Gaussian derivative at scale sigma = 0.899 results.
2022
2023 Postconditions:
2024 \code
2025 1. left() == -2
2026 2. right() == 2
2027 3. borderTreatment() == BORDER_TREATMENT_REFLECT
2028 4. norm() == 1.0
2029 \endcode
2030 */
2031 void initOptimalFirstDerivative5()
2032 {
2033 this->initExplicitly(-2, 2) = 0.1, 0.3, 0.0, -0.3, -0.1;
2034 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
2035 }
2036
2037 /**
2038 Init an optimal 5-tap second derivative filter.
2039 This filter must be used in conjunction with the corresponding 5-tap smoothing filter
2040 (see initOptimalSecondDerivativeSmoothing5()), such that the derivative filter is applied along one dimension,
2041 and the smoothing filter along the other.
2042 The filter values are
2043
2044 \code
2045 [0.22075, 0.117, -0.6755, 0.117, 0.22075]
2046 \endcode
2047
2048 These values are optimal in the sense that the 5x5 filter obtained by combining
2049 this filter with the corresponding 5-tap smoothing filter is the best possible 5x5 approximation to a
2050 Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.817.
2051
2052 Postconditions:
2053 \code
2054 1. left() == -2
2055 2. right() == 2
2056 3. borderTreatment() == BORDER_TREATMENT_REFLECT
2057 4. norm() == 1.0
2058 \endcode
2059 */
2060 void initOptimalSecondDerivative5()
2061 {
2062 this->initExplicitly(-2, 2) = 0.22075, 0.117, -0.6755, 0.117, 0.22075;
2063 this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
2064 }
2065
2066 /** Init the kernel by an explicit initializer list.
2067 The left and right boundaries of the kernel must be passed.
2068 A comma-separated initializer list is given after the assignment
2069 operator. This function is used like this:
2070
2071 \code
2072 // define horizontal Roberts filter
2073 vigra::Kernel1D<float> roberts_gradient_x;
2074
2075 roberts_gradient_x.initExplicitly(0, 1) = 1.0, -1.0;
2076 \endcode
2077
2078 The norm is set to the sum of the initializer values. If the wrong number of
2079 values is given, a run-time error results. It is, however, possible to give
2080 just one initializer. This creates an averaging filter with the given constant:
2081
2082 \code
2083 vigra::Kernel1D<float> average5x1;
2084
2085 average5x1.initExplicitly(-2, 2) = 1.0/5.0;
2086 \endcode
2087
2088 Here, the norm is set to value*size().
2089
2090 <b> Preconditions:</b>
2091
2092 \code
2093
2094 1. left <= 0
2095 2. right >= 0
2096 3. the number of values in the initializer list
2097 is 1 or equals the size of the kernel.
2098 \endcode
2099 */
2100 Kernel1D & initExplicitly(int left, int right)
2101 {
2102 vigra_precondition(left <= 0,
2103 "Kernel1D::initExplicitly(): left border must be <= 0.");
2104 vigra_precondition(right >= 0,
2105 "Kernel1D::initExplicitly(): right border must be >= 0.");
2106
2107 right_ = right;
2108 left_ = left;
2109
2110 kernel_.resize(right - left + 1);
2111
2112 return *this;
2113 }
2114
2115 /** Get iterator to center of kernel
2116
2117 Postconditions:
2118 \code
2119
2120 center()[left()] ... center()[right()] are valid kernel positions
2121 \endcode
2122 */
2123 iterator center()
2124 {
2125 return kernel_.begin() - left();
2126 }
2127
2128 const_iterator center() const
2129 {
2130 return kernel_.begin() - left();
2131 }
2132
2133 /** Access kernel value at specified location.
2134
2135 Preconditions:
2136 \code
2137
2138 left() <= location <= right()
2139 \endcode
2140 */
2141 reference operator[](int location)
2142 {
2143 return kernel_[location - left()];
2144 }
2145
2146 const_reference operator[](int location) const
2147 {
2148 return kernel_[location - left()];
2149 }
2150
2151 /** left border of kernel (inclusive), always <= 0
2152 */
2153 int left() const { return left_; }
2154
2155 /** right border of kernel (inclusive), always >= 0
2156 */
2157 int right() const { return right_; }
2158
2159 /** size of kernel (right() - left() + 1)
2160 */
2161 int size() const { return right_ - left_ + 1; }
2162
2163 /** current border treatment mode
2164 */
2165 BorderTreatmentMode borderTreatment() const
2166 { return border_treatment_; }
2167
2168 /** Set border treatment mode.
2169 */
2170 void setBorderTreatment( BorderTreatmentMode new_mode)
2171 { border_treatment_ = new_mode; }
2172
2173 /** norm of kernel
2174 */
2175 value_type norm() const { return norm_; }
2176
2177 /** set a new norm and normalize kernel, use the normalization formula
2178 for the given <tt>derivativeOrder</tt>.
2179 */
2180 void
2181 normalize(value_type norm, unsigned int derivativeOrder = 0, double offset = 0.0);
2182
2183 /** normalize kernel to norm 1.
2184 */
2185 void
2186 normalize()
2187 {
2188 normalize(one());
2189 }
2190
2191 /** get a const accessor
2192 */
2193 ConstAccessor accessor() const { return ConstAccessor(); }
2194
2195 /** get an accessor
2196 */
2197 Accessor accessor() { return Accessor(); }
2198
2199 private:
2200 InternalVector kernel_;
2201 int left_, right_;
2202 BorderTreatmentMode border_treatment_;
2203 value_type norm_;
2204};
2205
2206 template <class ARITHTYPE>
2207 void Kernel1D<ARITHTYPE>::normalize(value_type norm,
2208 unsigned int derivativeOrder,
2209 double offset)
2210{
2211 typedef typename NumericTraits<value_type>::RealPromote TmpType;
2212
2213 // find kernel sum
2214 Iterator k = kernel_.begin();
2215 TmpType sum = NumericTraits<TmpType>::zero();
2216
2217 if(derivativeOrder == 0)
2218 {
2219 for(; k < kernel_.end(); ++k)
2220 {
2221 sum += *k;
2222 }
2223 }
2224 else
2225 {
2226 unsigned int faculty = 1;
2227 for(unsigned int i = 2; i <= derivativeOrder; ++i)
2228 faculty *= i;
2229 for(double x = left() + offset; k < kernel_.end(); ++x, ++k)
2230 {
2231 sum = TmpType(sum + *k * VIGRA_CSTD::pow(-x, int(derivativeOrder)) / faculty);
2232 }
2233 }
2234
2235 vigra_precondition(sum != NumericTraits<value_type>::zero(),
2236 "Kernel1D<ARITHTYPE>::normalize(): "
2237 "Cannot normalize a kernel with sum = 0");
2238 // normalize
2239 sum = norm / sum;
2240 k = kernel_.begin();
2241 for(; k != kernel_.end(); ++k)
2242 {
2243 *k = *k * sum;
2244 }
2245
2246 norm_ = norm;
2247}
2248
2249/***********************************************************************/
2250
2251 template <class ARITHTYPE>
2252 void
2253 Kernel1D<ARITHTYPE>::initGaussian(double std_dev,
2254 value_type norm,
2255 double windowRatio)
2256{
2257 vigra_precondition(std_dev >= 0.0,
2258 "Kernel1D::initGaussian(): Standard deviation must be >= 0.");
2259 vigra_precondition(windowRatio >= 0.0,
2260 "Kernel1D::initGaussian(): windowRatio must be >= 0.");
2261
2262 if(std_dev > 0.0)
2263 {
2264 Gaussian<ARITHTYPE> gauss((ARITHTYPE)std_dev);
2265
2266 // first calculate required kernel sizes
2267 int radius;
2268 if (windowRatio == 0.0)
2269 radius = (int)(3.0 * std_dev + 0.5);
2270 else
2271 radius = (int)(windowRatio * std_dev + 0.5);
2272 if(radius == 0)
2273 radius = 1;
2274
2275 // allocate the kernel
2276 kernel_.erase(kernel_.begin(), kernel_.end());
2277 kernel_.reserve(radius*2+1);
2278
2279 for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
2280 {
2281 kernel_.push_back(gauss(x));
2282 }
2283 left_ = -radius;
2284 right_ = radius;
2285 }
2286 else
2287 {
2288 kernel_.erase(kernel_.begin(), kernel_.end());
2289 kernel_.push_back(1.0);
2290 left_ = 0;
2291 right_ = 0;
2292 }
2293
2294 if(norm != 0.0)
2295 normalize(norm);
2296 else
2297 norm_ = 1.0;
2298
2299 // best border treatment for Gaussians is BORDER_TREATMENT_REFLECT
2300 border_treatment_ = BORDER_TREATMENT_REFLECT;
2301}
2302
2303/***********************************************************************/
2304
2305 template <class ARITHTYPE>
2306 void
2307 Kernel1D<ARITHTYPE>::initDiscreteGaussian(double std_dev,
2308 value_type norm)
2309{
2310 vigra_precondition(std_dev >= 0.0,
2311 "Kernel1D::initDiscreteGaussian(): Standard deviation must be >= 0.");
2312
2313 if(std_dev > 0.0)
2314 {
2315 // first calculate required kernel sizes
2316 int radius = (int)(3.0*std_dev + 0.5);
2317 if(radius == 0)
2318 radius = 1;
2319
2320 double f = 2.0 / std_dev / std_dev;
2321
2322 // allocate the working array
2323 int maxIndex = (int)(2.0 * (radius + 5.0 * VIGRA_CSTD::sqrt((double)radius)) + 0.5);
2324 InternalVector warray(maxIndex+1);
2325 warray[maxIndex] = 0.0;
2326 warray[maxIndex-1] = 1.0;
2327
2328 for(int i = maxIndex-2; i >= radius; --i)
2329 {
2330 warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
2331 if(warray[i] > 1.0e40)
2332 {
2333 warray[i+1] /= warray[i];
2334 warray[i] = 1.0;
2335 }
2336 }
2337
2338 // the following rescaling ensures that the numbers stay in a sensible range
2339 // during the rest of the iteration, so no other rescaling is needed
2340 double er = VIGRA_CSTD::exp(-radius*radius / (2.0*std_dev*std_dev));
2341 warray[radius+1] = er * warray[radius+1] / warray[radius];
2342 warray[radius] = er;
2343
2344 for(int i = radius-1; i >= 0; --i)
2345 {
2346 warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
2347 er += warray[i];
2348 }
2349
2350 double scale = norm / (2*er - warray[0]);
2351
2352 initExplicitly(-radius, radius);
2353 iterator c = center();
2354
2355 for(int i=0; i<=radius; ++i)
2356 {
2357 c[i] = c[-i] = warray[i] * scale;
2358 }
2359 }
2360 else
2361 {
2362 kernel_.erase(kernel_.begin(), kernel_.end());
2363 kernel_.push_back(norm);
2364 left_ = 0;
2365 right_ = 0;
2366 }
2367
2368 norm_ = norm;
2369
2370 // best border treatment for Gaussians is BORDER_TREATMENT_REFLECT
2371 border_treatment_ = BORDER_TREATMENT_REFLECT;
2372}
2373
2374/***********************************************************************/
2375
2376 template <class ARITHTYPE>
2377 void
2378 Kernel1D<ARITHTYPE>::initGaussianDerivative(double std_dev,
2379 int order,
2380 value_type norm,
2381 double windowRatio)
2382{
2383 vigra_precondition(order >= 0,
2384 "Kernel1D::initGaussianDerivative(): Order must be >= 0.");
2385
2386 if(order == 0)
2387 {
2388 initGaussian(std_dev, norm, windowRatio);
2389 return;
2390 }
2391
2392 vigra_precondition(std_dev > 0.0,
2393 "Kernel1D::initGaussianDerivative(): "
2394 "Standard deviation must be > 0.");
2395 vigra_precondition(windowRatio >= 0.0,
2396 "Kernel1D::initGaussianDerivative(): windowRatio must be >= 0.");
2397
2398 Gaussian<ARITHTYPE> gauss((ARITHTYPE)std_dev, order);
2399
2400 // first calculate required kernel sizes
2401 int radius;
2402 if(windowRatio == 0.0)
2403 radius = (int)((3.0 + 0.5 * order) * std_dev + 0.5);
2404 else
2405 radius = (int)(windowRatio * std_dev + 0.5);
2406 if(radius == 0)
2407 radius = 1;
2408
2409 // allocate the kernels
2410 kernel_.clear();
2411 kernel_.reserve(radius*2+1);
2412
2413 // fill the kernel and calculate the DC component
2414 // introduced by truncation of the Gaussian
2415 ARITHTYPE dc = 0.0;
2416 for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
2417 {
2418 kernel_.push_back(gauss(x));
2419 dc += kernel_[kernel_.size()-1];
2420 }
2421 dc = ARITHTYPE(dc / (2.0*radius + 1.0));
2422
2423 // remove DC, but only if kernel correction is permitted by a non-zero
2424 // value for norm
2425 if(norm != 0.0)
2426 {
2427 for(unsigned int i=0; i < kernel_.size(); ++i)
2428 {
2429 kernel_[i] -= dc;
2430 }
2431 }
2432
2433 left_ = -radius;
2434 right_ = radius;
2435
2436 if(norm != 0.0)
2437 normalize(norm, order);
2438 else
2439 norm_ = 1.0;
2440
2441 // best border treatment for Gaussian derivatives is
2442 // BORDER_TREATMENT_REFLECT
2443 border_treatment_ = BORDER_TREATMENT_REFLECT;
2444}
2445
2446/***********************************************************************/
2447
2448 template <class ARITHTYPE>
2449 void
2450 Kernel1D<ARITHTYPE>::initBinomial(int radius,
2451 value_type norm)
2452{
2453 vigra_precondition(radius > 0,
2454 "Kernel1D::initBinomial(): Radius must be > 0.");
2455
2456 // allocate the kernel
2457 InternalVector(radius*2+1).swap(kernel_);
2458 typename InternalVector::iterator x = kernel_.begin() + radius;
2459
2460 // fill kernel
2461 x[radius] = norm;
2462 for(int j=radius-1; j>=-radius; --j)
2463 {
2464 x[j] = 0.5 * x[j+1];
2465 for(int i=j+1; i<radius; ++i)
2466 {
2467 x[i] = 0.5 * (x[i] + x[i+1]);
2468 }
2469 x[radius] *= 0.5;
2470 }
2471
2472 left_ = -radius;
2473 right_ = radius;
2474 norm_ = norm;
2475
2476 // best border treatment for Binomial is BORDER_TREATMENT_REFLECT
2477 border_treatment_ = BORDER_TREATMENT_REFLECT;
2478}
2479
2480/***********************************************************************/
2481
2482 template <class ARITHTYPE>
2483 void
2484 Kernel1D<ARITHTYPE>::initAveraging(int radius,
2485 value_type norm)
2486{
2487 vigra_precondition(radius > 0,
2488 "Kernel1D::initAveraging(): Radius must be > 0.");
2489
2490 // calculate scaling
2491 double scale = 1.0 / (radius * 2 + 1);
2492
2493 // normalize
2494 kernel_.erase(kernel_.begin(), kernel_.end());
2495 kernel_.reserve(radius*2+1);
2496
2497 for(int i=0; i<=radius*2+1; ++i)
2498 {
2499 kernel_.push_back(scale * norm);
2500 }
2501
2502 left_ = -radius;
2503 right_ = radius;
2504 norm_ = norm;
2505
2506 // best border treatment for Averaging is BORDER_TREATMENT_CLIP
2507 border_treatment_ = BORDER_TREATMENT_CLIP;
2508}
2509
2510/***********************************************************************/
2511
2512 template <class ARITHTYPE>
2513 void
2514 Kernel1D<ARITHTYPE>::initSymmetricDifference(value_type norm)
2515{
2516 kernel_.erase(kernel_.begin(), kernel_.end());
2517 kernel_.reserve(3);
2518
2519 kernel_.push_back(ARITHTYPE(0.5 * norm));
2520 kernel_.push_back(ARITHTYPE(0.0 * norm));
2521 kernel_.push_back(ARITHTYPE(-0.5 * norm));
2522
2523 left_ = -1;
2524 right_ = 1;
2525 norm_ = norm;
2526
2527 // best border treatment for symmetric difference is
2528 // BORDER_TREATMENT_REFLECT
2529 border_treatment_ = BORDER_TREATMENT_REFLECT;
2530}
2531
2532/**************************************************************/
2533/* */
2534/* Argument object factories for Kernel1D */
2535/* */
2536/* (documentation: see vigra/convolution.hxx) */
2537/* */
2538/**************************************************************/
2539
2540 template <class KernelIterator, class KernelAccessor>
2541 inline
2542tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>
2543kernel1d(KernelIterator ik, KernelAccessor ka,
2544 int kleft, int kright, BorderTreatmentMode border)
2545{
2546 return
2547 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>(
2548 ik, ka, kleft, kright, border);
2549}
2550
2551 template <class T>
2552 inline
2553tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
2554 int, int, BorderTreatmentMode>
2555kernel1d(Kernel1D<T> const & k)
2556
2557{
2558 return
2559 tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
2560 int, int, BorderTreatmentMode>(
2561 k.center(),
2562 k.accessor(),
2563 k.left(), k.right(),
2564 k.borderTreatment());
2565}
2566
2567 template <class T>
2568 inline
2569tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
2570 int, int, BorderTreatmentMode>
2571kernel1d(Kernel1D<T> const & k, BorderTreatmentMode border)
2572
2573{
2574 return
2575 tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
2576 int, int, BorderTreatmentMode>(
2577 k.center(),
2578 k.accessor(),
2579 k.left(), k.right(),
2580 border);
2581}
2582
2583
2584} // namespace vigra
2585
2586#endif // VIGRA_SEPARABLECONVOLUTION_HXX
vigra::Kernel1D
Generic 1 dimensional convolution kernel.
Definition separableconvolution.hxx:1367
vigra::Kernel1D::initBinomial
void initBinomial(int radius)
Definition separableconvolution.hxx:1853
vigra::Kernel1D::initOptimalFirstDerivativeSmoothing5
void initOptimalFirstDerivativeSmoothing5()
Definition separableconvolution.hxx:1761
vigra::Kernel1D::initSecondDifference3
void initSecondDifference3()
Definition separableconvolution.hxx:1999
vigra::Kernel1D::initOptimalSecondDerivativeSmoothing3
void initOptimalSecondDerivativeSmoothing3()
Definition separableconvolution.hxx:1707
vigra::Kernel1D::reference
InternalVector::reference reference
Definition separableconvolution.hxx:1377
vigra::Kernel1D::initBurtFilter
void initBurtFilter(double a=0.04785)
Definition separableconvolution.hxx:1824
vigra::Kernel1D::initDiscreteGaussian
void initDiscreteGaussian(double std_dev)
Definition separableconvolution.hxx:1584
vigra::Kernel1D::initBackwardDifference
void initBackwardDifference()
Definition separableconvolution.hxx:1955
vigra::Kernel1D::initOptimalFirstDerivative5
void initOptimalFirstDerivative5()
Definition separableconvolution.hxx:2031
vigra::Kernel1D::right
int right() const
Definition separableconvolution.hxx:2157
vigra::Kernel1D::norm
value_type norm() const
Definition separableconvolution.hxx:2175
vigra::Kernel1D::operator[]
reference operator[](int location)
Definition separableconvolution.hxx:2141
vigra::Kernel1D::Kernel1D
Kernel1D(Kernel1D< U > const &k)
Definition separableconvolution.hxx:1474
vigra::Kernel1D::initGaussian
void initGaussian(double std_dev, value_type norm, double windowRatio=0.0)
Definition separableconvolution.hxx:2253
vigra::Kernel1D::operator=
InitProxy operator=(value_type const &v)
Definition separableconvolution.hxx:1513
vigra::Kernel1D::initSymmetricGradient
void initSymmetricGradient()
Definition separableconvolution.hxx:1908
vigra::Kernel1D::borderTreatment
BorderTreatmentMode borderTreatment() const
Definition separableconvolution.hxx:2165
vigra::Kernel1D::setBorderTreatment
void setBorderTreatment(BorderTreatmentMode new_mode)
Definition separableconvolution.hxx:2170
vigra::Kernel1D::Accessor
StandardAccessor< ARITHTYPE > Accessor
Definition separableconvolution.hxx:1397
vigra::Kernel1D::initOptimalSmoothing5
void initOptimalSmoothing5()
Definition separableconvolution.hxx:1733
vigra::Kernel1D::accessor
ConstAccessor accessor() const
Definition separableconvolution.hxx:2193
vigra::Kernel1D::initGaussianDerivative
void initGaussianDerivative(double std_dev, int order)
Definition separableconvolution.hxx:1626
vigra::Kernel1D::initDiscreteGaussian
void initDiscreteGaussian(double std_dev, value_type norm)
Definition separableconvolution.hxx:2307
vigra::Kernel1D::value_type
InternalVector::value_type value_type
Definition separableconvolution.hxx:1373
vigra::Kernel1D::initGaussianDerivative
void initGaussianDerivative(double std_dev, int order, value_type norm, double windowRatio=0.0)
Definition separableconvolution.hxx:2378
vigra::Kernel1D::iterator
InternalVector::iterator iterator
Definition separableconvolution.hxx:1389
vigra::Kernel1D::~Kernel1D
~Kernel1D()
Definition separableconvolution.hxx:1524
vigra::Kernel1D::initSymmetricDifference
void initSymmetricDifference()
Definition separableconvolution.hxx:1978
vigra::Kernel1D::initAveraging
void initAveraging(int radius)
Definition separableconvolution.hxx:1879
vigra::Kernel1D::ConstAccessor
StandardConstAccessor< ARITHTYPE > ConstAccessor
Definition separableconvolution.hxx:1401
vigra::Kernel1D::Kernel1D
Kernel1D(Kernel1D const &k)
Definition separableconvolution.hxx:1463
vigra::Kernel1D::initOptimalSecondDerivative5
void initOptimalSecondDerivative5()
Definition separableconvolution.hxx:2060
vigra::Kernel1D::const_iterator
InternalVector::const_iterator const_iterator
Definition separableconvolution.hxx:1393
vigra::Kernel1D::const_reference
InternalVector::const_reference const_reference
Definition separableconvolution.hxx:1381
vigra::Kernel1D::initAveraging
void initAveraging(int radius, value_type norm)
Definition separableconvolution.hxx:2484
vigra::Kernel1D::initSymmetricGradient
void initSymmetricGradient(value_type norm)
Definition separableconvolution.hxx:1898
vigra::Kernel1D::initGaussian
void initGaussian(double std_dev)
Definition separableconvolution.hxx:1556
vigra::Kernel1D::initOptimalSecondDerivativeSmoothing5
void initOptimalSecondDerivativeSmoothing5()
Definition separableconvolution.hxx:1789
vigra::Kernel1D::left
int left() const
Definition separableconvolution.hxx:2153
vigra::Kernel1D::accessor
Accessor accessor()
Definition separableconvolution.hxx:2197
vigra::Kernel1D::Kernel1D
Kernel1D()
Definition separableconvolution.hxx:1451
vigra::Kernel1D::initBinomial
void initBinomial(int radius, value_type norm)
Definition separableconvolution.hxx:2450
vigra::Kernel1D::normalize
void normalize()
Definition separableconvolution.hxx:2186
vigra::Kernel1D::initForwardDifference
void initForwardDifference()
Definition separableconvolution.hxx:1931
vigra::Kernel1D::Iterator
InternalVector::iterator Iterator
Definition separableconvolution.hxx:1385
vigra::Kernel1D::initExplicitly
Kernel1D & initExplicitly(int left, int right)
Definition separableconvolution.hxx:2100
vigra::Kernel1D::initOptimalSmoothing3
void initOptimalSmoothing3()
Definition separableconvolution.hxx:1651
vigra::Kernel1D::initOptimalFirstDerivativeSmoothing3
void initOptimalFirstDerivativeSmoothing3()
Definition separableconvolution.hxx:1679
vigra::Kernel1D::size
int size() const
Definition separableconvolution.hxx:2161
vigra::Kernel1D::center
iterator center()
Definition separableconvolution.hxx:2123
vigra::Kernel1D::operator=
Kernel1D & operator=(Kernel1D const &k)
Definition separableconvolution.hxx:1484
vigra::RGBValue
Class for a single RGB value.
Definition rgbvalue.hxx:128
vigra::RGBValue::RGBValue
RGBValue()
Definition rgbvalue.hxx:209
vigra::TinyVectorBase::size
size_type size() const
Definition tinyvector.hxx:913
vigra::TinyVectorBase::end
iterator end()
Definition tinyvector.hxx:864
vigra::TinyVectorBase::begin
iterator begin()
Definition tinyvector.hxx:861
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
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.2