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

© 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.10.0