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

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

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

html generated using doxygen and Python
vigra 1.9.0 (Sat Oct 5 2013)