[ 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 /* internalConvolveLineWrap */
53 /* */
54 /********************************************************/
55 
56 template <class SrcIterator, class SrcAccessor,
57  class DestIterator, class DestAccessor,
58  class KernelIterator, class KernelAccessor>
59 void internalConvolveLineWrap(SrcIterator is, SrcIterator iend, SrcAccessor sa,
60  DestIterator id, DestAccessor da,
61  KernelIterator kernel, KernelAccessor ka,
62  int kleft, int kright,
63  int start = 0, int stop = 0)
64 {
65  int w = std::distance( is, iend );
66 
67  typedef typename PromoteTraits<
68  typename SrcAccessor::value_type,
69  typename KernelAccessor::value_type>::Promote SumType;
70 
71  SrcIterator ibegin = is;
72 
73  if(stop == 0)
74  stop = w;
75  is += start;
76 
77  for(int x=start; x<stop; ++x, ++is, ++id)
78  {
79  KernelIterator ik = kernel + kright;
80  SumType sum = NumericTraits<SumType>::zero();
81 
82  if(x < kright)
83  {
84  int x0 = x - kright;
85  SrcIterator iss = iend + x0;
86 
87  for(; x0; ++x0, --ik, ++iss)
88  {
89  sum += ka(ik) * sa(iss);
90  }
91 
92  iss = ibegin;
93  SrcIterator isend = is + (1 - kleft);
94  for(; iss != isend ; --ik, ++iss)
95  {
96  sum += ka(ik) * sa(iss);
97  }
98  }
99  else if(w-x <= -kleft)
100  {
101  SrcIterator iss = is + (-kright);
102  SrcIterator isend = iend;
103  for(; iss != isend ; --ik, ++iss)
104  {
105  sum += ka(ik) * sa(iss);
106  }
107 
108  int x0 = -kleft - w + x + 1;
109  iss = ibegin;
110 
111  for(; x0; --x0, --ik, ++iss)
112  {
113  sum += ka(ik) * sa(iss);
114  }
115  }
116  else
117  {
118  SrcIterator iss = is - kright;
119  SrcIterator isend = is + (1 - kleft);
120  for(; iss != isend ; --ik, ++iss)
121  {
122  sum += ka(ik) * sa(iss);
123  }
124  }
125 
126  da.set(detail::RequiresExplicitCast<typename
127  DestAccessor::value_type>::cast(sum), id);
128  }
129 }
130 
131 /********************************************************/
132 /* */
133 /* internalConvolveLineClip */
134 /* */
135 /********************************************************/
136 
137 template <class SrcIterator, class SrcAccessor,
138  class DestIterator, class DestAccessor,
139  class KernelIterator, class KernelAccessor,
140  class Norm>
141 void internalConvolveLineClip(SrcIterator is, SrcIterator iend, SrcAccessor sa,
142  DestIterator id, DestAccessor da,
143  KernelIterator kernel, KernelAccessor ka,
144  int kleft, int kright, Norm norm,
145  int start = 0, int stop = 0)
146 {
147  int w = std::distance( is, iend );
148 
149  typedef typename PromoteTraits<
150  typename SrcAccessor::value_type,
151  typename KernelAccessor::value_type>::Promote SumType;
152 
153  SrcIterator ibegin = is;
154 
155  if(stop == 0)
156  stop = w;
157  is += start;
158 
159  for(int x=start; x<stop; ++x, ++is, ++id)
160  {
161  KernelIterator ik = kernel + kright;
162  SumType sum = NumericTraits<SumType>::zero();
163 
164  if(x < kright)
165  {
166  int x0 = x - kright;
167  Norm clipped = NumericTraits<Norm>::zero();
168 
169  for(; x0; ++x0, --ik)
170  {
171  clipped += ka(ik);
172  }
173 
174  SrcIterator iss = ibegin;
175  SrcIterator isend = is + (1 - kleft);
176  for(; iss != isend ; --ik, ++iss)
177  {
178  sum += ka(ik) * sa(iss);
179  }
180 
181  sum = norm / (norm - clipped) * sum;
182  }
183  else if(w-x <= -kleft)
184  {
185  SrcIterator iss = is + (-kright);
186  SrcIterator isend = iend;
187  for(; iss != isend ; --ik, ++iss)
188  {
189  sum += ka(ik) * sa(iss);
190  }
191 
192  Norm clipped = NumericTraits<Norm>::zero();
193 
194  int x0 = -kleft - w + x + 1;
195 
196  for(; x0; --x0, --ik)
197  {
198  clipped += ka(ik);
199  }
200 
201  sum = norm / (norm - clipped) * sum;
202  }
203  else
204  {
205  SrcIterator iss = is + (-kright);
206  SrcIterator isend = is + (1 - kleft);
207  for(; iss != isend ; --ik, ++iss)
208  {
209  sum += ka(ik) * sa(iss);
210  }
211  }
212 
213  da.set(detail::RequiresExplicitCast<typename
214  DestAccessor::value_type>::cast(sum), id);
215  }
216 }
217 
218 /********************************************************/
219 /* */
220 /* internalConvolveLineReflect */
221 /* */
222 /********************************************************/
223 
224 template <class SrcIterator, class SrcAccessor,
225  class DestIterator, class DestAccessor,
226  class KernelIterator, class KernelAccessor>
227 void internalConvolveLineReflect(SrcIterator is, SrcIterator iend, SrcAccessor sa,
228  DestIterator id, DestAccessor da,
229  KernelIterator kernel, KernelAccessor ka,
230  int kleft, int kright,
231  int start = 0, int stop = 0)
232 {
233  int w = std::distance( is, iend );
234 
235  typedef typename PromoteTraits<
236  typename SrcAccessor::value_type,
237  typename KernelAccessor::value_type>::Promote SumType;
238 
239  SrcIterator ibegin = is;
240 
241  if(stop == 0)
242  stop = w;
243  is += start;
244 
245  for(int x=start; x<stop; ++x, ++is, ++id)
246  {
247  KernelIterator ik = kernel + kright;
248  SumType sum = NumericTraits<SumType>::zero();
249 
250  if(x < kright)
251  {
252  int x0 = x - kright;
253  SrcIterator iss = ibegin - x0;
254 
255  for(; x0; ++x0, --ik, --iss)
256  {
257  sum += ka(ik) * sa(iss);
258  }
259 
260  SrcIterator isend = is + (1 - kleft);
261  for(; iss != isend ; --ik, ++iss)
262  {
263  sum += ka(ik) * sa(iss);
264  }
265  }
266  else if(w-x <= -kleft)
267  {
268  SrcIterator iss = is + (-kright);
269  SrcIterator isend = iend;
270  for(; iss != isend ; --ik, ++iss)
271  {
272  sum += ka(ik) * sa(iss);
273  }
274 
275  int x0 = -kleft - w + x + 1;
276  iss = iend - 2;
277 
278  for(; x0; --x0, --ik, --iss)
279  {
280  sum += ka(ik) * sa(iss);
281  }
282  }
283  else
284  {
285  SrcIterator iss = is + (-kright);
286  SrcIterator isend = is + (1 - kleft);
287  for(; iss != isend ; --ik, ++iss)
288  {
289  sum += ka(ik) * sa(iss);
290  }
291  }
292 
293  da.set(detail::RequiresExplicitCast<typename
294  DestAccessor::value_type>::cast(sum), id);
295  }
296 }
297 
298 /********************************************************/
299 /* */
300 /* internalConvolveLineRepeat */
301 /* */
302 /********************************************************/
303 
304 template <class SrcIterator, class SrcAccessor,
305  class DestIterator, class DestAccessor,
306  class KernelIterator, class KernelAccessor>
307 void internalConvolveLineRepeat(SrcIterator is, SrcIterator iend, SrcAccessor sa,
308  DestIterator id, DestAccessor da,
309  KernelIterator kernel, KernelAccessor ka,
310  int kleft, int kright,
311  int start = 0, int stop = 0)
312 {
313  int w = std::distance( is, iend );
314 
315  typedef typename PromoteTraits<
316  typename SrcAccessor::value_type,
317  typename KernelAccessor::value_type>::Promote SumType;
318 
319  SrcIterator ibegin = is;
320 
321  if(stop == 0)
322  stop = w;
323  is += start;
324 
325  for(int x=start; x<stop; ++x, ++is, ++id)
326  {
327  KernelIterator ik = kernel + kright;
328  SumType sum = NumericTraits<SumType>::zero();
329 
330  if(x < kright)
331  {
332  int x0 = x - kright;
333  SrcIterator iss = ibegin;
334 
335  for(; x0; ++x0, --ik)
336  {
337  sum += ka(ik) * sa(iss);
338  }
339 
340  SrcIterator isend = is + (1 - kleft);
341  for(; iss != isend ; --ik, ++iss)
342  {
343  sum += ka(ik) * sa(iss);
344  }
345  }
346  else if(w-x <= -kleft)
347  {
348  SrcIterator iss = is + (-kright);
349  SrcIterator isend = iend;
350  for(; iss != isend ; --ik, ++iss)
351  {
352  sum += ka(ik) * sa(iss);
353  }
354 
355  int x0 = -kleft - w + x + 1;
356  iss = iend - 1;
357 
358  for(; x0; --x0, --ik)
359  {
360  sum += ka(ik) * sa(iss);
361  }
362  }
363  else
364  {
365  SrcIterator iss = is + (-kright);
366  SrcIterator isend = is + (1 - kleft);
367  for(; iss != isend ; --ik, ++iss)
368  {
369  sum += ka(ik) * sa(iss);
370  }
371  }
372 
373  da.set(detail::RequiresExplicitCast<typename
374  DestAccessor::value_type>::cast(sum), id);
375  }
376 }
377 
378 /********************************************************/
379 /* */
380 /* internalConvolveLineAvoid */
381 /* */
382 /********************************************************/
383 
384 template <class SrcIterator, class SrcAccessor,
385  class DestIterator, class DestAccessor,
386  class KernelIterator, class KernelAccessor>
387 void internalConvolveLineAvoid(SrcIterator is, SrcIterator iend, SrcAccessor sa,
388  DestIterator id, DestAccessor da,
389  KernelIterator kernel, KernelAccessor ka,
390  int kleft, int kright,
391  int start = 0, int stop = 0)
392 {
393  int w = std::distance( is, iend );
394  if(start < stop) // we got a valid subrange
395  {
396  if(w + kleft < stop)
397  stop = w + kleft;
398  if(start < kright)
399  {
400  id += kright - start;
401  start = kright;
402  }
403  }
404  else
405  {
406  id += kright;
407  start = kright;
408  stop = w + kleft;
409  }
410 
411  typedef typename PromoteTraits<
412  typename SrcAccessor::value_type,
413  typename KernelAccessor::value_type>::Promote SumType;
414 
415  is += start;
416 
417  for(int x=start; x<stop; ++x, ++is, ++id)
418  {
419  KernelIterator ik = kernel + kright;
420  SumType sum = NumericTraits<SumType>::zero();
421 
422  SrcIterator iss = is + (-kright);
423  SrcIterator isend = is + (1 - kleft);
424  for(; iss != isend ; --ik, ++iss)
425  {
426  sum += ka(ik) * sa(iss);
427  }
428 
429  da.set(detail::RequiresExplicitCast<typename
430  DestAccessor::value_type>::cast(sum), id);
431  }
432 }
433 
434 /********************************************************/
435 /* */
436 /* Separable convolution functions */
437 /* */
438 /********************************************************/
439 
440 /** \addtogroup SeparableConvolution One-dimensional and separable convolution functions
441 
442  Perform 1D convolution and separable filtering in 2 dimensions.
443 
444  These generic convolution functions implement
445  the standard convolution operation for a wide range of images and
446  signals that fit into the required interface. They need a suitable
447  kernel to operate.
448 */
449 //@{
450 
451 /** \brief Performs a 1-dimensional convolution of the source signal using the given
452  kernel.
453 
454  The KernelIterator must point to the center iterator, and
455  the kernel's size is given by its left (kleft <= 0) and right
456  (kright >= 0) borders. The signal must always be larger than the kernel.
457  At those positions where the kernel does not completely fit
458  into the signal's range, the specified \ref BorderTreatmentMode is
459  applied.
460 
461  The signal's value_type (SrcAccessor::value_type) must be a
462  linear space over the kernel's value_type (KernelAccessor::value_type),
463  i.e. addition of source values, multiplication with kernel values,
464  and NumericTraits must be defined.
465  The kernel's value_type must be an algebraic field,
466  i.e. the arithmetic operations (+, -, *, /) and NumericTraits must
467  be defined.
468 
469  If <tt>start</tt> and <tt>stop</tt> are non-zero, the relation
470  <tt>0 <= start < stop <= width</tt> must hold (where <tt>width</tt>
471  is the length of the input array). The convolution is then restricted to that
472  subrange, and it is assumed that the output array only refers to that
473  subrange (i.e. <tt>id</tt> points to the element corresponding to
474  <tt>start</tt>). If <tt>start</tt> and <tt>stop</tt> are both zero
475  (the default), the entire array is convolved.
476 
477  <b> Declarations:</b>
478 
479  pass arguments explicitly:
480  \code
481  namespace vigra {
482  template <class SrcIterator, class SrcAccessor,
483  class DestIterator, class DestAccessor,
484  class KernelIterator, class KernelAccessor>
485  void convolveLine(SrcIterator is, SrcIterator isend, SrcAccessor sa,
486  DestIterator id, DestAccessor da,
487  KernelIterator ik, KernelAccessor ka,
488  int kleft, int kright, BorderTreatmentMode border,
489  int start = 0, int stop = 0 )
490  }
491  \endcode
492 
493 
494  use argument objects in conjunction with \ref ArgumentObjectFactories :
495  \code
496  namespace vigra {
497  template <class SrcIterator, class SrcAccessor,
498  class DestIterator, class DestAccessor,
499  class KernelIterator, class KernelAccessor>
500  void convolveLine(triple<SrcIterator, SrcIterator, SrcAccessor> src,
501  pair<DestIterator, DestAccessor> dest,
502  tuple5<KernelIterator, KernelAccessor,
503  int, int, BorderTreatmentMode> kernel,
504  int start = 0, int stop = 0)
505  }
506  \endcode
507 
508  <b> Usage:</b>
509 
510  <b>\#include</b> <vigra/separableconvolution.hxx>
511 
512 
513  \code
514  std::vector<float> src, dest;
515  ...
516 
517  // define binomial filter of size 5
518  static float kernel[] =
519  { 1.0/16.0, 4.0/16.0, 6.0/16.0, 4.0/16.0, 1.0/16.0};
520 
521  typedef vigra::StandardAccessor<float> FAccessor;
522  typedef vigra::StandardAccessor<float> KernelAccessor;
523 
524 
525  vigra::convolveLine(src.begin(), src.end(), FAccessor(), dest.begin(), FAccessor(),
526  kernel+2, KernelAccessor(), -2, 2, BORDER_TREATMENT_REFLECT);
527  // ^^^^^^^^ this is the center of the kernel
528 
529  \endcode
530 
531  <b> Required Interface:</b>
532 
533  \code
534  RandomAccessIterator is, isend;
535  RandomAccessIterator id;
536  RandomAccessIterator ik;
537 
538  SrcAccessor src_accessor;
539  DestAccessor dest_accessor;
540  KernelAccessor kernel_accessor;
541 
542  NumericTraits<SrcAccessor::value_type>::RealPromote s = src_accessor(is);
543 
544  s = s + s;
545  s = kernel_accessor(ik) * s;
546 
547  dest_accessor.set(
548  NumericTraits<DestAccessor::value_type>::fromRealPromote(s), id);
549 
550  \endcode
551 
552  If border == BORDER_TREATMENT_CLIP:
553 
554  \code
555  NumericTraits<KernelAccessor::value_type>::RealPromote k = kernel_accessor(ik);
556 
557  k = k + k;
558  k = k - k;
559  k = k * k;
560  k = k / k;
561 
562  \endcode
563 
564  <b> Preconditions:</b>
565 
566  \code
567  kleft <= 0
568  kright >= 0
569  iend - is >= kright + kleft + 1
570  \endcode
571 
572  If border == BORDER_TREATMENT_CLIP: Sum of kernel elements must be
573  != 0.
574 
575 */
576 doxygen_overloaded_function(template <...> void convolveLine)
577 
578 template <class SrcIterator, class SrcAccessor,
579  class DestIterator, class DestAccessor,
580  class KernelIterator, class KernelAccessor>
581 void convolveLine(SrcIterator is, SrcIterator iend, SrcAccessor sa,
582  DestIterator id, DestAccessor da,
583  KernelIterator ik, KernelAccessor ka,
584  int kleft, int kright, BorderTreatmentMode border,
585  int start = 0, int stop = 0)
586 {
587  typedef typename KernelAccessor::value_type KernelValue;
588 
589  vigra_precondition(kleft <= 0,
590  "convolveLine(): kleft must be <= 0.\n");
591  vigra_precondition(kright >= 0,
592  "convolveLine(): kright must be >= 0.\n");
593 
594  // int w = iend - is;
595  int w = std::distance( is, iend );
596 
597  vigra_precondition(w >= std::max(kright, -kleft) + 1,
598  "convolveLine(): kernel longer than line.\n");
599 
600  if(stop != 0)
601  vigra_precondition(0 <= start && start < stop && stop <= w,
602  "convolveLine(): invalid subrange (start, stop).\n");
603 
604  switch(border)
605  {
606  case BORDER_TREATMENT_WRAP:
607  {
608  internalConvolveLineWrap(is, iend, sa, id, da, ik, ka, kleft, kright, start, stop);
609  break;
610  }
611  case BORDER_TREATMENT_AVOID:
612  {
613  internalConvolveLineAvoid(is, iend, sa, id, da, ik, ka, kleft, kright, start, stop);
614  break;
615  }
616  case BORDER_TREATMENT_REFLECT:
617  {
618  internalConvolveLineReflect(is, iend, sa, id, da, ik, ka, kleft, kright, start, stop);
619  break;
620  }
621  case BORDER_TREATMENT_REPEAT:
622  {
623  internalConvolveLineRepeat(is, iend, sa, id, da, ik, ka, kleft, kright, start, stop);
624  break;
625  }
626  case BORDER_TREATMENT_CLIP:
627  {
628  // find norm of kernel
629  typedef typename KernelAccessor::value_type KT;
630  KT norm = NumericTraits<KT>::zero();
631  KernelIterator iik = ik + kleft;
632  for(int i=kleft; i<=kright; ++i, ++iik)
633  norm += ka(iik);
634 
635  vigra_precondition(norm != NumericTraits<KT>::zero(),
636  "convolveLine(): Norm of kernel must be != 0"
637  " in mode BORDER_TREATMENT_CLIP.\n");
638 
639  internalConvolveLineClip(is, iend, sa, id, da, ik, ka, kleft, kright, norm, start, stop);
640  break;
641  }
642  default:
643  {
644  vigra_precondition(0,
645  "convolveLine(): Unknown border treatment mode.\n");
646  }
647  }
648 }
649 
650 template <class SrcIterator, class SrcAccessor,
651  class DestIterator, class DestAccessor,
652  class KernelIterator, class KernelAccessor>
653 inline
654 void convolveLine(triple<SrcIterator, SrcIterator, SrcAccessor> src,
655  pair<DestIterator, DestAccessor> dest,
656  tuple5<KernelIterator, KernelAccessor,
657  int, int, BorderTreatmentMode> kernel,
658  int start = 0, int stop = 0)
659 {
660  convolveLine(src.first, src.second, src.third,
661  dest.first, dest.second,
662  kernel.first, kernel.second,
663  kernel.third, kernel.fourth, kernel.fifth, start, stop);
664 }
665 
666 /********************************************************/
667 /* */
668 /* separableConvolveX */
669 /* */
670 /********************************************************/
671 
672 /** \brief Performs a 1 dimensional convolution in x direction.
673 
674  It calls \ref convolveLine() for every row of the image. See \ref convolveLine()
675  for more information about required interfaces and vigra_preconditions.
676 
677  <b> Declarations:</b>
678 
679  pass arguments explicitly:
680  \code
681  namespace vigra {
682  template <class SrcImageIterator, class SrcAccessor,
683  class DestImageIterator, class DestAccessor,
684  class KernelIterator, class KernelAccessor>
685  void separableConvolveX(SrcImageIterator supperleft,
686  SrcImageIterator slowerright, SrcAccessor sa,
687  DestImageIterator dupperleft, DestAccessor da,
688  KernelIterator ik, KernelAccessor ka,
689  int kleft, int kright, BorderTreatmentMode border)
690  }
691  \endcode
692 
693 
694  use argument objects in conjunction with \ref ArgumentObjectFactories :
695  \code
696  namespace vigra {
697  template <class SrcImageIterator, class SrcAccessor,
698  class DestImageIterator, class DestAccessor,
699  class KernelIterator, class KernelAccessor>
700  void separableConvolveX(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
701  pair<DestImageIterator, DestAccessor> dest,
702  tuple5<KernelIterator, KernelAccessor,
703  int, int, BorderTreatmentMode> kernel)
704  }
705  \endcode
706 
707  <b> Usage:</b>
708 
709  <b>\#include</b> <vigra/separableconvolution.hxx>
710 
711 
712  \code
713  vigra::FImage src(w,h), dest(w,h);
714  ...
715 
716  // define Gaussian kernel with std. deviation 3.0
717  vigra::Kernel1D<double> kernel;
718  kernel.initGaussian(3.0);
719 
720  vigra::separableConvolveX(srcImageRange(src), destImage(dest), kernel1d(kernel));
721 
722  \endcode
723 
724 */
725 doxygen_overloaded_function(template <...> void separableConvolveX)
726 
727 template <class SrcIterator, class SrcAccessor,
728  class DestIterator, class DestAccessor,
729  class KernelIterator, class KernelAccessor>
730 void separableConvolveX(SrcIterator supperleft,
731  SrcIterator slowerright, SrcAccessor sa,
732  DestIterator dupperleft, DestAccessor da,
733  KernelIterator ik, KernelAccessor ka,
734  int kleft, int kright, BorderTreatmentMode border)
735 {
736  typedef typename KernelAccessor::value_type KernelValue;
737 
738  vigra_precondition(kleft <= 0,
739  "separableConvolveX(): kleft must be <= 0.\n");
740  vigra_precondition(kright >= 0,
741  "separableConvolveX(): kright must be >= 0.\n");
742 
743  int w = slowerright.x - supperleft.x;
744  int h = slowerright.y - supperleft.y;
745 
746  vigra_precondition(w >= std::max(kright, -kleft) + 1,
747  "separableConvolveX(): kernel longer than line\n");
748 
749  int y;
750 
751  for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
752  {
753  typename SrcIterator::row_iterator rs = supperleft.rowIterator();
754  typename DestIterator::row_iterator rd = dupperleft.rowIterator();
755 
756  convolveLine(rs, rs+w, sa, rd, da,
757  ik, ka, kleft, kright, border);
758  }
759 }
760 
761 template <class SrcIterator, class SrcAccessor,
762  class DestIterator, class DestAccessor,
763  class KernelIterator, class KernelAccessor>
764 inline void
765 separableConvolveX(triple<SrcIterator, SrcIterator, SrcAccessor> src,
766  pair<DestIterator, DestAccessor> dest,
767  tuple5<KernelIterator, KernelAccessor,
768  int, int, BorderTreatmentMode> kernel)
769 {
770  separableConvolveX(src.first, src.second, src.third,
771  dest.first, dest.second,
772  kernel.first, kernel.second,
773  kernel.third, kernel.fourth, kernel.fifth);
774 }
775 
776 
777 
778 /********************************************************/
779 /* */
780 /* separableConvolveY */
781 /* */
782 /********************************************************/
783 
784 /** \brief Performs a 1 dimensional convolution in y direction.
785 
786  It calls \ref convolveLine() for every column of the image. See \ref convolveLine()
787  for more information about required interfaces and vigra_preconditions.
788 
789  <b> Declarations:</b>
790 
791  pass arguments explicitly:
792  \code
793  namespace vigra {
794  template <class SrcImageIterator, class SrcAccessor,
795  class DestImageIterator, class DestAccessor,
796  class KernelIterator, class KernelAccessor>
797  void separableConvolveY(SrcImageIterator supperleft,
798  SrcImageIterator slowerright, SrcAccessor sa,
799  DestImageIterator dupperleft, DestAccessor da,
800  KernelIterator ik, KernelAccessor ka,
801  int kleft, int kright, BorderTreatmentMode border)
802  }
803  \endcode
804 
805 
806  use argument objects in conjunction with \ref ArgumentObjectFactories :
807  \code
808  namespace vigra {
809  template <class SrcImageIterator, class SrcAccessor,
810  class DestImageIterator, class DestAccessor,
811  class KernelIterator, class KernelAccessor>
812  void separableConvolveY(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
813  pair<DestImageIterator, DestAccessor> dest,
814  tuple5<KernelIterator, KernelAccessor,
815  int, int, BorderTreatmentMode> kernel)
816  }
817  \endcode
818 
819  <b> Usage:</b>
820 
821  <b>\#include</b> <vigra/separableconvolution.hxx>
822 
823 
824  \code
825  vigra::FImage src(w,h), dest(w,h);
826  ...
827 
828  // define Gaussian kernel with std. deviation 3.0
829  vigra::Kernel1D kernel;
830  kernel.initGaussian(3.0);
831 
832  vigra::separableConvolveY(srcImageRange(src), destImage(dest), kernel1d(kernel));
833 
834  \endcode
835 
836 */
837 doxygen_overloaded_function(template <...> void separableConvolveY)
838 
839 template <class SrcIterator, class SrcAccessor,
840  class DestIterator, class DestAccessor,
841  class KernelIterator, class KernelAccessor>
842 void separableConvolveY(SrcIterator supperleft,
843  SrcIterator slowerright, SrcAccessor sa,
844  DestIterator dupperleft, DestAccessor da,
845  KernelIterator ik, KernelAccessor ka,
846  int kleft, int kright, BorderTreatmentMode border)
847 {
848  typedef typename KernelAccessor::value_type KernelValue;
849 
850  vigra_precondition(kleft <= 0,
851  "separableConvolveY(): kleft must be <= 0.\n");
852  vigra_precondition(kright >= 0,
853  "separableConvolveY(): kright must be >= 0.\n");
854 
855  int w = slowerright.x - supperleft.x;
856  int h = slowerright.y - supperleft.y;
857 
858  vigra_precondition(h >= std::max(kright, -kleft) + 1,
859  "separableConvolveY(): kernel longer than line\n");
860 
861  int x;
862 
863  for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
864  {
865  typename SrcIterator::column_iterator cs = supperleft.columnIterator();
866  typename DestIterator::column_iterator cd = dupperleft.columnIterator();
867 
868  convolveLine(cs, cs+h, sa, cd, da,
869  ik, ka, kleft, kright, border);
870  }
871 }
872 
873 template <class SrcIterator, class SrcAccessor,
874  class DestIterator, class DestAccessor,
875  class KernelIterator, class KernelAccessor>
876 inline void
877 separableConvolveY(triple<SrcIterator, SrcIterator, SrcAccessor> src,
878  pair<DestIterator, DestAccessor> dest,
879  tuple5<KernelIterator, KernelAccessor,
880  int, int, BorderTreatmentMode> kernel)
881 {
882  separableConvolveY(src.first, src.second, src.third,
883  dest.first, dest.second,
884  kernel.first, kernel.second,
885  kernel.third, kernel.fourth, kernel.fifth);
886 }
887 
888 //@}
889 
890 /********************************************************/
891 /* */
892 /* Kernel1D */
893 /* */
894 /********************************************************/
895 
896 /** \brief Generic 1 dimensional convolution kernel.
897 
898  This kernel may be used for convolution of 1 dimensional signals or for
899  separable convolution of multidimensional signals.
900 
901  Convolution functions access the kernel via a 1 dimensional random access
902  iterator which they get by calling \ref center(). This iterator
903  points to the center of the kernel. The kernel's size is given by its left() (<=0)
904  and right() (>= 0) methods. The desired border treatment mode is
905  returned by borderTreatment().
906 
907  The different init functions create a kernel with the specified
908  properties. The kernel's value_type must be a linear space, i.e. it
909  must define multiplication with doubles and NumericTraits.
910 
911 
912  The kernel defines a factory function kernel1d() to create an argument object
913  (see \ref KernelArgumentObjectFactories).
914 
915  <b> Usage:</b>
916 
917  <b>\#include</b> <vigra/stdconvolution.hxx>
918 
919  \code
920  vigra::FImage src(w,h), dest(w,h);
921  ...
922 
923  // define Gaussian kernel with std. deviation 3.0
924  vigra::Kernel1D kernel;
925  kernel.initGaussian(3.0);
926 
927  vigra::separableConvolveX(srcImageRange(src), destImage(dest), kernel1d(kernel));
928  \endcode
929 
930  <b> Required Interface:</b>
931 
932  \code
933  value_type v = vigra::NumericTraits<value_type>::one(); // if norm is not
934  // given explicitly
935  double d;
936 
937  v = d * v;
938  \endcode
939 */
940 
941 template <class ARITHTYPE>
942 class Kernel1D
943 {
944  public:
946 
947  /** the kernel's value type
948  */
949  typedef typename InternalVector::value_type value_type;
950 
951  /** the kernel's reference type
952  */
953  typedef typename InternalVector::reference reference;
954 
955  /** the kernel's const reference type
956  */
957  typedef typename InternalVector::const_reference const_reference;
958 
959  /** deprecated -- use Kernel1D::iterator
960  */
961  typedef typename InternalVector::iterator Iterator;
962 
963  /** 1D random access iterator over the kernel's values
964  */
965  typedef typename InternalVector::iterator iterator;
966 
967  /** const 1D random access iterator over the kernel's values
968  */
969  typedef typename InternalVector::const_iterator const_iterator;
970 
971  /** the kernel's accessor
972  */
974 
975  /** the kernel's const accessor
976  */
978 
979  struct InitProxy
980  {
981  InitProxy(Iterator i, int count, value_type & norm)
982  : iter_(i), base_(i),
983  count_(count), sum_(count),
984  norm_(norm)
985  {}
986 
987  ~InitProxy()
988  {
989  vigra_precondition(count_ == 1 || count_ == sum_,
990  "Kernel1D::initExplicitly(): "
991  "Wrong number of init values.");
992  }
993 
994  InitProxy & operator,(value_type const & v)
995  {
996  if(sum_ == count_)
997  norm_ = *iter_;
998 
999  norm_ += v;
1000 
1001  --count_;
1002 
1003  if(count_ > 0)
1004  {
1005  ++iter_;
1006  *iter_ = v;
1007  }
1008 
1009  return *this;
1010  }
1011 
1012  Iterator iter_, base_;
1013  int count_, sum_;
1014  value_type & norm_;
1015  };
1016 
1017  static value_type one() { return NumericTraits<value_type>::one(); }
1018 
1019  /** Default constructor.
1020  Creates a kernel of size 1 which would copy the signal
1021  unchanged.
1022  */
1024  : kernel_(),
1025  left_(0),
1026  right_(0),
1027  border_treatment_(BORDER_TREATMENT_REFLECT),
1028  norm_(one())
1029  {
1030  kernel_.push_back(norm_);
1031  }
1032 
1033  /** Copy constructor.
1034  */
1035  Kernel1D(Kernel1D const & k)
1036  : kernel_(k.kernel_),
1037  left_(k.left_),
1038  right_(k.right_),
1039  border_treatment_(k.border_treatment_),
1040  norm_(k.norm_)
1041  {}
1042 
1043  /** Construct from kernel with different element type, e.g. double => FixedPoint16.
1044  */
1045  template <class U>
1047  : kernel_(k.center()+k.left(), k.center()+k.right()+1),
1048  left_(k.left()),
1049  right_(k.right()),
1050  border_treatment_(k.borderTreatment()),
1051  norm_(k.norm())
1052  {}
1053 
1054  /** Copy assignment.
1055  */
1057  {
1058  if(this != &k)
1059  {
1060  left_ = k.left_;
1061  right_ = k.right_;
1062  border_treatment_ = k.border_treatment_;
1063  norm_ = k.norm_;
1064  kernel_ = k.kernel_;
1065  }
1066  return *this;
1067  }
1068 
1069  /** Initialization.
1070  This initializes the kernel with the given constant. The norm becomes
1071  v*size().
1072 
1073  Instead of a single value an initializer list of length size()
1074  can be used like this:
1075 
1076  \code
1077  vigra::Kernel1D<float> roberts_gradient_x;
1078 
1079  roberts_gradient_x.initExplicitly(0, 1) = 1.0, -1.0;
1080  \endcode
1081 
1082  In this case, the norm will be set to the sum of the init values.
1083  An initializer list of wrong length will result in a run-time error.
1084  */
1085  InitProxy operator=(value_type const & v)
1086  {
1087  int size = right_ - left_ + 1;
1088  for(unsigned int i=0; i<kernel_.size(); ++i) kernel_[i] = v;
1089  norm_ = (double)size*v;
1090 
1091  return InitProxy(kernel_.begin(), size, norm_);
1092  }
1093 
1094  /** Destructor.
1095  */
1097  {}
1098 
1099  /**
1100  Init as a sampled Gaussian function. The radius of the kernel is
1101  always 3*std_dev. '<tt>norm</tt>' denotes the sum of all bins of the kernel
1102  (i.e. the kernel is corrected for the normalization error introduced
1103  by windowing the Gaussian to a finite interval). However,
1104  if <tt>norm</tt> is 0.0, the kernel is normalized to 1 by the analytic
1105  expression for the Gaussian, and <b>no</b> correction for the windowing
1106  error is performed. If <tt>windowRatio = 0.0</tt>, the radius of the filter
1107  window is <tt>radius = round(3.0 * std_dev)</tt>, otherwise it is
1108  <tt>radius = round(windowRatio * std_dev)</tt> (where <tt>windowRatio > 0.0</tt>
1109  is required).
1110 
1111  Precondition:
1112  \code
1113  std_dev >= 0.0
1114  \endcode
1115 
1116  Postconditions:
1117  \code
1118  1. left() == -(int)(3.0*std_dev + 0.5)
1119  2. right() == (int)(3.0*std_dev + 0.5)
1120  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1121  4. norm() == norm
1122  \endcode
1123  */
1124  void initGaussian(double std_dev, value_type norm, double windowRatio = 0.0);
1125 
1126  /** Init as a Gaussian function with norm 1.
1127  */
1128  void initGaussian(double std_dev)
1129  {
1130  initGaussian(std_dev, one());
1131  }
1132 
1133 
1134  /**
1135  Init as Lindeberg's discrete analog of the Gaussian function. The radius of the kernel is
1136  always 3*std_dev. 'norm' denotes the sum of all bins of the kernel.
1137 
1138  Precondition:
1139  \code
1140  std_dev >= 0.0
1141  \endcode
1142 
1143  Postconditions:
1144  \code
1145  1. left() == -(int)(3.0*std_dev + 0.5)
1146  2. right() == (int)(3.0*std_dev + 0.5)
1147  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1148  4. norm() == norm
1149  \endcode
1150  */
1151  void initDiscreteGaussian(double std_dev, value_type norm);
1152 
1153  /** Init as a Lindeberg's discrete analog of the Gaussian function
1154  with norm 1.
1155  */
1156  void initDiscreteGaussian(double std_dev)
1157  {
1158  initDiscreteGaussian(std_dev, one());
1159  }
1160 
1161  /**
1162  Init as a Gaussian derivative of order '<tt>order</tt>'.
1163  The radius of the kernel is always <tt>3*std_dev + 0.5*order</tt>.
1164  '<tt>norm</tt>' denotes the norm of the kernel so that the
1165  following condition is fulfilled:
1166 
1167  \f[ \sum_{i=left()}^{right()}
1168  \frac{(-i)^{order}kernel[i]}{order!} = norm
1169  \f]
1170 
1171  Thus, the kernel will be corrected for the error introduced
1172  by windowing the Gaussian to a finite interval. However,
1173  if <tt>norm</tt> is 0.0, the kernel is normalized to 1 by the analytic
1174  expression for the Gaussian derivative, and <b>no</b> correction for the
1175  windowing error is performed. If <tt>windowRatio = 0.0</tt>, the radius
1176  of the filter window is <tt>radius = round(3.0 * std_dev + 0.5 * order)</tt>,
1177  otherwise it is <tt>radius = round(windowRatio * std_dev)</tt> (where
1178  <tt>windowRatio > 0.0</tt> is required).
1179 
1180  Preconditions:
1181  \code
1182  1. std_dev >= 0.0
1183  2. order >= 1
1184  \endcode
1185 
1186  Postconditions:
1187  \code
1188  1. left() == -(int)(3.0*std_dev + 0.5*order + 0.5)
1189  2. right() == (int)(3.0*std_dev + 0.5*order + 0.5)
1190  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1191  4. norm() == norm
1192  \endcode
1193  */
1194  void initGaussianDerivative(double std_dev, int order, value_type norm, double windowRatio = 0.0);
1195 
1196  /** Init as a Gaussian derivative with norm 1.
1197  */
1198  void initGaussianDerivative(double std_dev, int order)
1199  {
1200  initGaussianDerivative(std_dev, order, one());
1201  }
1202 
1203  /**
1204  Init an optimal 3-tap smoothing filter.
1205  The filter values are
1206 
1207  \code
1208  [0.216, 0.568, 0.216]
1209  \endcode
1210 
1211  These values are optimal in the sense that the 3x3 filter obtained by separable application
1212  of this filter is the best possible 3x3 approximation to a Gaussian filter.
1213  The equivalent Gaussian has sigma = 0.680.
1214 
1215  Postconditions:
1216  \code
1217  1. left() == -1
1218  2. right() == 1
1219  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1220  4. norm() == 1.0
1221  \endcode
1222  */
1224  {
1225  this->initExplicitly(-1, 1) = 0.216, 0.568, 0.216;
1226  this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1227  }
1228 
1229  /**
1230  Init an optimal 3-tap smoothing filter to be used in the context of first derivative computation.
1231  This filter must be used in conjunction with the symmetric difference filter (see initSymmetricDifference()),
1232  such that the difference filter is applied along one dimension, and the smoothing filter along the other.
1233  The filter values are
1234 
1235  \code
1236  [0.224365, 0.55127, 0.224365]
1237  \endcode
1238 
1239  These values are optimal in the sense that the 3x3 filter obtained by combining
1240  this filter with the symmetric difference is the best possible 3x3 approximation to a
1241  Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.675.
1242 
1243  Postconditions:
1244  \code
1245  1. left() == -1
1246  2. right() == 1
1247  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1248  4. norm() == 1.0
1249  \endcode
1250  */
1252  {
1253  this->initExplicitly(-1, 1) = 0.224365, 0.55127, 0.224365;
1254  this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1255  }
1256 
1257  /**
1258  Init an optimal 3-tap smoothing filter to be used in the context of second derivative computation.
1259  This filter must be used in conjunction with the 3-tap second difference filter (see initSecondDifference3()),
1260  such that the difference filter is applied along one dimension, and the smoothing filter along the other.
1261  The filter values are
1262 
1263  \code
1264  [0.13, 0.74, 0.13]
1265  \endcode
1266 
1267  These values are optimal in the sense that the 3x3 filter obtained by combining
1268  this filter with the 3-tap second difference is the best possible 3x3 approximation to a
1269  Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.433.
1270 
1271  Postconditions:
1272  \code
1273  1. left() == -1
1274  2. right() == 1
1275  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1276  4. norm() == 1.0
1277  \endcode
1278  */
1280  {
1281  this->initExplicitly(-1, 1) = 0.13, 0.74, 0.13;
1282  this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1283  }
1284 
1285  /**
1286  Init an optimal 5-tap smoothing filter.
1287  The filter values are
1288 
1289  \code
1290  [0.03134, 0.24, 0.45732, 0.24, 0.03134]
1291  \endcode
1292 
1293  These values are optimal in the sense that the 5x5 filter obtained by separable application
1294  of this filter is the best possible 5x5 approximation to a Gaussian filter.
1295  The equivalent Gaussian has sigma = 0.867.
1296 
1297  Postconditions:
1298  \code
1299  1. left() == -2
1300  2. right() == 2
1301  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1302  4. norm() == 1.0
1303  \endcode
1304  */
1306  {
1307  this->initExplicitly(-2, 2) = 0.03134, 0.24, 0.45732, 0.24, 0.03134;
1308  this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1309  }
1310 
1311  /**
1312  Init an optimal 5-tap smoothing filter to be used in the context of first derivative computation.
1313  This filter must be used in conjunction with the optimal 5-tap first derivative filter
1314  (see initOptimalFirstDerivative5()), such that the derivative filter is applied along one dimension,
1315  and the smoothing filter along the other. The filter values are
1316 
1317  \code
1318  [0.04255, 0.241, 0.4329, 0.241, 0.04255]
1319  \endcode
1320 
1321  These values are optimal in the sense that the 5x5 filter obtained by combining
1322  this filter with the optimal 5-tap first derivative is the best possible 5x5 approximation to a
1323  Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.906.
1324 
1325  Postconditions:
1326  \code
1327  1. left() == -2
1328  2. right() == 2
1329  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1330  4. norm() == 1.0
1331  \endcode
1332  */
1334  {
1335  this->initExplicitly(-2, 2) = 0.04255, 0.241, 0.4329, 0.241, 0.04255;
1336  this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1337  }
1338 
1339  /**
1340  Init an optimal 5-tap smoothing filter to be used in the context of second derivative computation.
1341  This filter must be used in conjunction with the optimal 5-tap second derivative filter
1342  (see initOptimalSecondDerivative5()), such that the derivative filter is applied along one dimension,
1343  and the smoothing filter along the other. The filter values are
1344 
1345  \code
1346  [0.0243, 0.23556, 0.48028, 0.23556, 0.0243]
1347  \endcode
1348 
1349  These values are optimal in the sense that the 5x5 filter obtained by combining
1350  this filter with the optimal 5-tap second derivative is the best possible 5x5 approximation to a
1351  Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.817.
1352 
1353  Postconditions:
1354  \code
1355  1. left() == -2
1356  2. right() == 2
1357  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1358  4. norm() == 1.0
1359  \endcode
1360  */
1362  {
1363  this->initExplicitly(-2, 2) = 0.0243, 0.23556, 0.48028, 0.23556, 0.0243;
1364  this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1365  }
1366 
1367  /**
1368  Init a 5-tap filter as defined by Peter Burt in the context of pyramid creation.
1369  The filter values are
1370 
1371  \code
1372  [a, 0.25, 0.5-2*a, 0.25, a]
1373  \endcode
1374 
1375  The default <tt>a = 0.04785</tt> is optimal in the sense that it minimizes the difference
1376  to a true Gaussian filter (which would have sigma = 0.975). For other values of <tt>a</tt>, the scale
1377  of the most similar Gaussian can be approximated by
1378 
1379  \code
1380  sigma = 5.1 * a + 0.731
1381  \endcode
1382 
1383  Preconditions:
1384  \code
1385  0 <= a <= 0.125
1386  \endcode
1387 
1388  Postconditions:
1389  \code
1390  1. left() == -2
1391  2. right() == 2
1392  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1393  4. norm() == 1.0
1394  \endcode
1395  */
1396  void initBurtFilter(double a = 0.04785)
1397  {
1398  vigra_precondition(a >= 0.0 && a <= 0.125,
1399  "Kernel1D::initBurtFilter(): 0 <= a <= 0.125 required.");
1400  this->initExplicitly(-2, 2) = a, 0.25, 0.5 - 2.0*a, 0.25, a;
1401  this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1402  }
1403 
1404  /**
1405  Init as a Binomial filter. 'norm' denotes the sum of all bins
1406  of the kernel.
1407 
1408  Precondition:
1409  \code
1410  radius >= 0
1411  \endcode
1412 
1413  Postconditions:
1414  \code
1415  1. left() == -radius
1416  2. right() == radius
1417  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1418  4. norm() == norm
1419  \endcode
1420  */
1421  void initBinomial(int radius, value_type norm);
1422 
1423  /** Init as a Binomial filter with norm 1.
1424  */
1425  void initBinomial(int radius)
1426  {
1427  initBinomial(radius, one());
1428  }
1429 
1430  /**
1431  Init as an Averaging filter. 'norm' denotes the sum of all bins
1432  of the kernel. The window size is (2*radius+1) * (2*radius+1)
1433 
1434  Precondition:
1435  \code
1436  radius >= 0
1437  \endcode
1438 
1439  Postconditions:
1440  \code
1441  1. left() == -radius
1442  2. right() == radius
1443  3. borderTreatment() == BORDER_TREATMENT_CLIP
1444  4. norm() == norm
1445  \endcode
1446  */
1447  void initAveraging(int radius, value_type norm);
1448 
1449  /** Init as an Averaging filter with norm 1.
1450  */
1451  void initAveraging(int radius)
1452  {
1453  initAveraging(radius, one());
1454  }
1455 
1456  /**
1457  Init as a symmetric gradient filter of the form
1458  <TT>[ 0.5 * norm, 0.0 * norm, -0.5 * norm]</TT>
1459 
1460  <b>Deprecated</b>. Use initSymmetricDifference() instead.
1461 
1462  Postconditions:
1463  \code
1464  1. left() == -1
1465  2. right() == 1
1466  3. borderTreatment() == BORDER_TREATMENT_REPEAT
1467  4. norm() == norm
1468  \endcode
1469  */
1470  void
1472  {
1474  setBorderTreatment(BORDER_TREATMENT_REPEAT);
1475  }
1476 
1477  /** Init as a symmetric gradient filter with norm 1.
1478 
1479  <b>Deprecated</b>. Use initSymmetricDifference() instead.
1480  */
1482  {
1483  initSymmetricGradient(one());
1484  }
1485 
1486  void
1488 
1489  /** Init as the 3-tap symmetric difference filter
1490  The filter values are
1491 
1492  \code
1493  [0.5, 0, -0.5]
1494  \endcode
1495 
1496  Postconditions:
1497  \code
1498  1. left() == -1
1499  2. right() == 1
1500  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1501  4. norm() == 1.0
1502  \endcode
1503  */
1505  {
1506  initSymmetricDifference(one());
1507  }
1508 
1509  /**
1510  Init the 3-tap second difference filter.
1511  The filter values are
1512 
1513  \code
1514  [1, -2, 1]
1515  \endcode
1516 
1517  Postconditions:
1518  \code
1519  1. left() == -1
1520  2. right() == 1
1521  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1522  4. norm() == 1
1523  \endcode
1524  */
1525  void
1527  {
1528  this->initExplicitly(-1, 1) = 1.0, -2.0, 1.0;
1529  this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1530  }
1531 
1532  /**
1533  Init an optimal 5-tap first derivative filter.
1534  This filter must be used in conjunction with the corresponding 5-tap smoothing filter
1535  (see initOptimalFirstDerivativeSmoothing5()), such that the derivative filter is applied along one dimension,
1536  and the smoothing filter along the other.
1537  The filter values are
1538 
1539  \code
1540  [0.1, 0.3, 0.0, -0.3, -0.1]
1541  \endcode
1542 
1543  These values are optimal in the sense that the 5x5 filter obtained by combining
1544  this filter with the corresponding 5-tap smoothing filter is the best possible 5x5 approximation to a
1545  Gaussian first derivative filter. The equivalent Gaussian has sigma = 0.906.
1546 
1547  If the filter is instead separably combined with itself, an almost optimal approximation of the
1548  mixed second Gaussian derivative at scale sigma = 0.899 results.
1549 
1550  Postconditions:
1551  \code
1552  1. left() == -2
1553  2. right() == 2
1554  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1555  4. norm() == 1.0
1556  \endcode
1557  */
1559  {
1560  this->initExplicitly(-2, 2) = 0.1, 0.3, 0.0, -0.3, -0.1;
1561  this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1562  }
1563 
1564  /**
1565  Init an optimal 5-tap second derivative filter.
1566  This filter must be used in conjunction with the corresponding 5-tap smoothing filter
1567  (see initOptimalSecondDerivativeSmoothing5()), such that the derivative filter is applied along one dimension,
1568  and the smoothing filter along the other.
1569  The filter values are
1570 
1571  \code
1572  [0.22075, 0.117, -0.6755, 0.117, 0.22075]
1573  \endcode
1574 
1575  These values are optimal in the sense that the 5x5 filter obtained by combining
1576  this filter with the corresponding 5-tap smoothing filter is the best possible 5x5 approximation to a
1577  Gaussian second derivative filter. The equivalent Gaussian has sigma = 0.817.
1578 
1579  Postconditions:
1580  \code
1581  1. left() == -2
1582  2. right() == 2
1583  3. borderTreatment() == BORDER_TREATMENT_REFLECT
1584  4. norm() == 1.0
1585  \endcode
1586  */
1588  {
1589  this->initExplicitly(-2, 2) = 0.22075, 0.117, -0.6755, 0.117, 0.22075;
1590  this->setBorderTreatment(BORDER_TREATMENT_REFLECT);
1591  }
1592 
1593  /** Init the kernel by an explicit initializer list.
1594  The left and right boundaries of the kernel must be passed.
1595  A comma-separated initializer list is given after the assignment
1596  operator. This function is used like this:
1597 
1598  \code
1599  // define horizontal Roberts filter
1600  vigra::Kernel1D<float> roberts_gradient_x;
1601 
1602  roberts_gradient_x.initExplicitly(0, 1) = 1.0, -1.0;
1603  \endcode
1604 
1605  The norm is set to the sum of the initializer values. If the wrong number of
1606  values is given, a run-time error results. It is, however, possible to give
1607  just one initializer. This creates an averaging filter with the given constant:
1608 
1609  \code
1610  vigra::Kernel1D<float> average5x1;
1611 
1612  average5x1.initExplicitly(-2, 2) = 1.0/5.0;
1613  \endcode
1614 
1615  Here, the norm is set to value*size().
1616 
1617  <b> Preconditions:</b>
1618 
1619  \code
1620 
1621  1. left <= 0
1622  2. right >= 0
1623  3. the number of values in the initializer list
1624  is 1 or equals the size of the kernel.
1625  \endcode
1626  */
1628  {
1629  vigra_precondition(left <= 0,
1630  "Kernel1D::initExplicitly(): left border must be <= 0.");
1631  vigra_precondition(right >= 0,
1632  "Kernel1D::initExplicitly(): right border must be >= 0.");
1633 
1634  right_ = right;
1635  left_ = left;
1636 
1637  kernel_.resize(right - left + 1);
1638 
1639  return *this;
1640  }
1641 
1642  /** Get iterator to center of kernel
1643 
1644  Postconditions:
1645  \code
1646 
1647  center()[left()] ... center()[right()] are valid kernel positions
1648  \endcode
1649  */
1651  {
1652  return kernel_.begin() - left();
1653  }
1654 
1655  const_iterator center() const
1656  {
1657  return kernel_.begin() - left();
1658  }
1659 
1660  /** Access kernel value at specified location.
1661 
1662  Preconditions:
1663  \code
1664 
1665  left() <= location <= right()
1666  \endcode
1667  */
1668  reference operator[](int location)
1669  {
1670  return kernel_[location - left()];
1671  }
1672 
1673  const_reference operator[](int location) const
1674  {
1675  return kernel_[location - left()];
1676  }
1677 
1678  /** left border of kernel (inclusive), always <= 0
1679  */
1680  int left() const { return left_; }
1681 
1682  /** right border of kernel (inclusive), always >= 0
1683  */
1684  int right() const { return right_; }
1685 
1686  /** size of kernel (right() - left() + 1)
1687  */
1688  int size() const { return right_ - left_ + 1; }
1689 
1690  /** current border treatment mode
1691  */
1692  BorderTreatmentMode borderTreatment() const
1693  { return border_treatment_; }
1694 
1695  /** Set border treatment mode.
1696  */
1697  void setBorderTreatment( BorderTreatmentMode new_mode)
1698  { border_treatment_ = new_mode; }
1699 
1700  /** norm of kernel
1701  */
1702  value_type norm() const { return norm_; }
1703 
1704  /** set a new norm and normalize kernel, use the normalization formula
1705  for the given <tt>derivativeOrder</tt>.
1706  */
1707  void
1708  normalize(value_type norm, unsigned int derivativeOrder = 0, double offset = 0.0);
1709 
1710  /** normalize kernel to norm 1.
1711  */
1712  void
1714  {
1715  normalize(one());
1716  }
1717 
1718  /** get a const accessor
1719  */
1720  ConstAccessor accessor() const { return ConstAccessor(); }
1721 
1722  /** get an accessor
1723  */
1724  Accessor accessor() { return Accessor(); }
1725 
1726  private:
1727  InternalVector kernel_;
1728  int left_, right_;
1729  BorderTreatmentMode border_treatment_;
1730  value_type norm_;
1731 };
1732 
1733 template <class ARITHTYPE>
1735  unsigned int derivativeOrder,
1736  double offset)
1737 {
1738  typedef typename NumericTraits<value_type>::RealPromote TmpType;
1739 
1740  // find kernel sum
1741  Iterator k = kernel_.begin();
1742  TmpType sum = NumericTraits<TmpType>::zero();
1743 
1744  if(derivativeOrder == 0)
1745  {
1746  for(; k < kernel_.end(); ++k)
1747  {
1748  sum += *k;
1749  }
1750  }
1751  else
1752  {
1753  unsigned int faculty = 1;
1754  for(unsigned int i = 2; i <= derivativeOrder; ++i)
1755  faculty *= i;
1756  for(double x = left() + offset; k < kernel_.end(); ++x, ++k)
1757  {
1758  sum = TmpType(sum + *k * VIGRA_CSTD::pow(-x, int(derivativeOrder)) / faculty);
1759  }
1760  }
1761 
1762  vigra_precondition(sum != NumericTraits<value_type>::zero(),
1763  "Kernel1D<ARITHTYPE>::normalize(): "
1764  "Cannot normalize a kernel with sum = 0");
1765  // normalize
1766  sum = norm / sum;
1767  k = kernel_.begin();
1768  for(; k != kernel_.end(); ++k)
1769  {
1770  *k = *k * sum;
1771  }
1772 
1773  norm_ = norm;
1774 }
1775 
1776 /***********************************************************************/
1777 
1778 template <class ARITHTYPE>
1780  value_type norm,
1781  double windowRatio)
1782 {
1783  vigra_precondition(std_dev >= 0.0,
1784  "Kernel1D::initGaussian(): Standard deviation must be >= 0.");
1785  vigra_precondition(windowRatio >= 0.0,
1786  "Kernel1D::initGaussian(): windowRatio must be >= 0.");
1787 
1788  if(std_dev > 0.0)
1789  {
1790  Gaussian<ARITHTYPE> gauss((ARITHTYPE)std_dev);
1791 
1792  // first calculate required kernel sizes
1793  int radius;
1794  if (windowRatio == 0.0)
1795  radius = (int)(3.0 * std_dev + 0.5);
1796  else
1797  radius = (int)(windowRatio * std_dev + 0.5);
1798  if(radius == 0)
1799  radius = 1;
1800 
1801  // allocate the kernel
1802  kernel_.erase(kernel_.begin(), kernel_.end());
1803  kernel_.reserve(radius*2+1);
1804 
1805  for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
1806  {
1807  kernel_.push_back(gauss(x));
1808  }
1809  left_ = -radius;
1810  right_ = radius;
1811  }
1812  else
1813  {
1814  kernel_.erase(kernel_.begin(), kernel_.end());
1815  kernel_.push_back(1.0);
1816  left_ = 0;
1817  right_ = 0;
1818  }
1819 
1820  if(norm != 0.0)
1821  normalize(norm);
1822  else
1823  norm_ = 1.0;
1824 
1825  // best border treatment for Gaussians is BORDER_TREATMENT_REFLECT
1826  border_treatment_ = BORDER_TREATMENT_REFLECT;
1827 }
1828 
1829 /***********************************************************************/
1830 
1831 template <class ARITHTYPE>
1833  value_type norm)
1834 {
1835  vigra_precondition(std_dev >= 0.0,
1836  "Kernel1D::initDiscreteGaussian(): Standard deviation must be >= 0.");
1837 
1838  if(std_dev > 0.0)
1839  {
1840  // first calculate required kernel sizes
1841  int radius = (int)(3.0*std_dev + 0.5);
1842  if(radius == 0)
1843  radius = 1;
1844 
1845  double f = 2.0 / std_dev / std_dev;
1846 
1847  // allocate the working array
1848  int maxIndex = (int)(2.0 * (radius + 5.0 * VIGRA_CSTD::sqrt((double)radius)) + 0.5);
1849  InternalVector warray(maxIndex+1);
1850  warray[maxIndex] = 0.0;
1851  warray[maxIndex-1] = 1.0;
1852 
1853  for(int i = maxIndex-2; i >= radius; --i)
1854  {
1855  warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
1856  if(warray[i] > 1.0e40)
1857  {
1858  warray[i+1] /= warray[i];
1859  warray[i] = 1.0;
1860  }
1861  }
1862 
1863  // the following rescaling ensures that the numbers stay in a sensible range
1864  // during the rest of the iteration, so no other rescaling is needed
1865  double er = VIGRA_CSTD::exp(-radius*radius / (2.0*std_dev*std_dev));
1866  warray[radius+1] = er * warray[radius+1] / warray[radius];
1867  warray[radius] = er;
1868 
1869  for(int i = radius-1; i >= 0; --i)
1870  {
1871  warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
1872  er += warray[i];
1873  }
1874 
1875  double scale = norm / (2*er - warray[0]);
1876 
1877  initExplicitly(-radius, radius);
1878  iterator c = center();
1879 
1880  for(int i=0; i<=radius; ++i)
1881  {
1882  c[i] = c[-i] = warray[i] * scale;
1883  }
1884  }
1885  else
1886  {
1887  kernel_.erase(kernel_.begin(), kernel_.end());
1888  kernel_.push_back(norm);
1889  left_ = 0;
1890  right_ = 0;
1891  }
1892 
1893  norm_ = norm;
1894 
1895  // best border treatment for Gaussians is BORDER_TREATMENT_REFLECT
1896  border_treatment_ = BORDER_TREATMENT_REFLECT;
1897 }
1898 
1899 /***********************************************************************/
1900 
1901 template <class ARITHTYPE>
1902 void
1904  int order,
1905  value_type norm,
1906  double windowRatio)
1907 {
1908  vigra_precondition(order >= 0,
1909  "Kernel1D::initGaussianDerivative(): Order must be >= 0.");
1910 
1911  if(order == 0)
1912  {
1913  initGaussian(std_dev, norm, windowRatio);
1914  return;
1915  }
1916 
1917  vigra_precondition(std_dev > 0.0,
1918  "Kernel1D::initGaussianDerivative(): "
1919  "Standard deviation must be > 0.");
1920  vigra_precondition(windowRatio >= 0.0,
1921  "Kernel1D::initGaussianDerivative(): windowRatio must be >= 0.");
1922 
1923  Gaussian<ARITHTYPE> gauss((ARITHTYPE)std_dev, order);
1924 
1925  // first calculate required kernel sizes
1926  int radius;
1927  if(windowRatio == 0.0)
1928  radius = (int)(3.0 * std_dev + 0.5 * order + 0.5);
1929  else
1930  radius = (int)(windowRatio * std_dev + 0.5);
1931  if(radius == 0)
1932  radius = 1;
1933 
1934  // allocate the kernels
1935  kernel_.clear();
1936  kernel_.reserve(radius*2+1);
1937 
1938  // fill the kernel and calculate the DC component
1939  // introduced by truncation of the Gaussian
1940  ARITHTYPE dc = 0.0;
1941  for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
1942  {
1943  kernel_.push_back(gauss(x));
1944  dc += kernel_[kernel_.size()-1];
1945  }
1946  dc = ARITHTYPE(dc / (2.0*radius + 1.0));
1947 
1948  // remove DC, but only if kernel correction is permitted by a non-zero
1949  // value for norm
1950  if(norm != 0.0)
1951  {
1952  for(unsigned int i=0; i < kernel_.size(); ++i)
1953  {
1954  kernel_[i] -= dc;
1955  }
1956  }
1957 
1958  left_ = -radius;
1959  right_ = radius;
1960 
1961  if(norm != 0.0)
1962  normalize(norm, order);
1963  else
1964  norm_ = 1.0;
1965 
1966  // best border treatment for Gaussian derivatives is
1967  // BORDER_TREATMENT_REFLECT
1968  border_treatment_ = BORDER_TREATMENT_REFLECT;
1969 }
1970 
1971 /***********************************************************************/
1972 
1973 template <class ARITHTYPE>
1974 void
1976  value_type norm)
1977 {
1978  vigra_precondition(radius > 0,
1979  "Kernel1D::initBinomial(): Radius must be > 0.");
1980 
1981  // allocate the kernel
1982  InternalVector(radius*2+1).swap(kernel_);
1983  typename InternalVector::iterator x = kernel_.begin() + radius;
1984 
1985  // fill kernel
1986  x[radius] = norm;
1987  for(int j=radius-1; j>=-radius; --j)
1988  {
1989  x[j] = 0.5 * x[j+1];
1990  for(int i=j+1; i<radius; ++i)
1991  {
1992  x[i] = 0.5 * (x[i] + x[i+1]);
1993  }
1994  x[radius] *= 0.5;
1995  }
1996 
1997  left_ = -radius;
1998  right_ = radius;
1999  norm_ = norm;
2000 
2001  // best border treatment for Binomial is BORDER_TREATMENT_REFLECT
2002  border_treatment_ = BORDER_TREATMENT_REFLECT;
2003 }
2004 
2005 /***********************************************************************/
2006 
2007 template <class ARITHTYPE>
2009  value_type norm)
2010 {
2011  vigra_precondition(radius > 0,
2012  "Kernel1D::initAveraging(): Radius must be > 0.");
2013 
2014  // calculate scaling
2015  double scale = 1.0 / (radius * 2 + 1);
2016 
2017  // normalize
2018  kernel_.erase(kernel_.begin(), kernel_.end());
2019  kernel_.reserve(radius*2+1);
2020 
2021  for(int i=0; i<=radius*2+1; ++i)
2022  {
2023  kernel_.push_back(scale * norm);
2024  }
2025 
2026  left_ = -radius;
2027  right_ = radius;
2028  norm_ = norm;
2029 
2030  // best border treatment for Averaging is BORDER_TREATMENT_CLIP
2031  border_treatment_ = BORDER_TREATMENT_CLIP;
2032 }
2033 
2034 /***********************************************************************/
2035 
2036 template <class ARITHTYPE>
2037 void
2039 {
2040  kernel_.erase(kernel_.begin(), kernel_.end());
2041  kernel_.reserve(3);
2042 
2043  kernel_.push_back(ARITHTYPE(0.5 * norm));
2044  kernel_.push_back(ARITHTYPE(0.0 * norm));
2045  kernel_.push_back(ARITHTYPE(-0.5 * norm));
2046 
2047  left_ = -1;
2048  right_ = 1;
2049  norm_ = norm;
2050 
2051  // best border treatment for symmetric difference is
2052  // BORDER_TREATMENT_REFLECT
2053  border_treatment_ = BORDER_TREATMENT_REFLECT;
2054 }
2055 
2056 /**************************************************************/
2057 /* */
2058 /* Argument object factories for Kernel1D */
2059 /* */
2060 /* (documentation: see vigra/convolution.hxx) */
2061 /* */
2062 /**************************************************************/
2063 
2064 template <class KernelIterator, class KernelAccessor>
2065 inline
2066 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>
2067 kernel1d(KernelIterator ik, KernelAccessor ka,
2068  int kleft, int kright, BorderTreatmentMode border)
2069 {
2070  return
2071  tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>(
2072  ik, ka, kleft, kright, border);
2073 }
2074 
2075 template <class T>
2076 inline
2077 tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
2078  int, int, BorderTreatmentMode>
2079 kernel1d(Kernel1D<T> const & k)
2080 
2081 {
2082  return
2083  tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
2084  int, int, BorderTreatmentMode>(
2085  k.center(),
2086  k.accessor(),
2087  k.left(), k.right(),
2088  k.borderTreatment());
2089 }
2090 
2091 template <class T>
2092 inline
2093 tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
2094  int, int, BorderTreatmentMode>
2095 kernel1d(Kernel1D<T> const & k, BorderTreatmentMode border)
2096 
2097 {
2098  return
2099  tuple5<typename Kernel1D<T>::const_iterator, typename Kernel1D<T>::ConstAccessor,
2100  int, int, BorderTreatmentMode>(
2101  k.center(),
2102  k.accessor(),
2103  k.left(), k.right(),
2104  border);
2105 }
2106 
2107 
2108 } // namespace vigra
2109 
2110 #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.8.0 (Wed Sep 26 2012)