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