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

separableconvolution.hxx VIGRA

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

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

html generated using doxygen and Python
vigra 1.11.1