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

stdconvolution.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_STDCONVOLUTION_HXX
38 #define VIGRA_STDCONVOLUTION_HXX
39 
40 #include <cmath>
41 #include "stdimage.hxx"
42 #include "bordertreatment.hxx"
43 #include "separableconvolution.hxx"
44 #include "utilities.hxx"
45 #include "sized_int.hxx"
46 
47 namespace vigra {
48 
49 template <class SrcIterator, class SrcAccessor,
50  class DestIterator, class DestAccessor,
51  class KernelIterator, class KernelAccessor,
52  class KSumType>
53 void internalPixelEvaluationByClip(int x, int y, int w, int h, SrcIterator xs,
54  SrcAccessor src_acc, DestIterator xd, DestAccessor dest_acc,
55  KernelIterator ki, Diff2D kul, Diff2D klr, KernelAccessor ak,
56  KSumType norm)
57 {
58  typedef typename
59  PromoteTraits<typename SrcAccessor::value_type,
60  typename KernelAccessor::value_type>::Promote SumType;
61  typedef typename DestAccessor::value_type DestType;
62 
63  // calculate width and height of the kernel
64  int kernel_width = klr.x - kul.x + 1;
65  int kernel_height = klr.y - kul.y + 1;
66 
67  SumType sum = NumericTraits<SumType>::zero();
68  int xx, yy;
69  int x0, y0, x1, y1;
70 
71  y0 = (y<klr.y) ? -y : -klr.y;
72  y1 = (h-y-1<-kul.y) ? h-y-1 : -kul.y;
73 
74  x0 = (x<klr.x) ? -x : -klr.x;
75  x1 = (w-x-1<-kul.x) ? w-x-1 : -kul.x;
76 
77  SrcIterator yys = xs + Diff2D(x0, y0);
78  KernelIterator yk = ki - Diff2D(x0, y0);
79 
80  KSumType ksum = NumericTraits<KSumType>::zero();
81  kernel_width = x1 - x0 + 1;
82  kernel_height = y1 - y0 + 1;
83 
84  //es wird zuerst abgeschnitten und dann gespigelt!
85 
86  for(yy=0; yy<kernel_height; ++yy, ++yys.y, --yk.y)
87  {
88  SrcIterator xxs = yys;
89  KernelIterator xk = yk;
90 
91  for(xx=0; xx<kernel_width; ++xx, ++xxs.x, --xk.x)
92  {
93  sum += ak(xk) * src_acc(xxs);
94  ksum += ak(xk);
95  }
96  }
97 
98  // store average in destination pixel
99  dest_acc.set(detail::RequiresExplicitCast<DestType>::cast((norm / ksum) * sum), xd);
100 }
101 
102 
103 #if 0
104 
105 template <class SrcIterator, class SrcAccessor,
106  class DestIterator, class DestAccessor,
107  class KernelIterator, class KernelAccessor>
108 void internalPixelEvaluationByWrapReflectRepeat(int x, int y, int src_width, int src_height, SrcIterator xs,
109  SrcAccessor src_acc, DestIterator xd, DestAccessor dest_acc,
110  KernelIterator ki, Diff2D kul, Diff2D klr, KernelAccessor ak,
111  BorderTreatmentMode border)
112 {
113 
114  typedef typename
115  NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType;
116  typedef
117  NumericTraits<typename DestAccessor::value_type> DestTraits;
118 
119  SumType sum = NumericTraits<SumType>::zero();
120 
121  SrcIterator src_ul = xs - Diff2D(x, y);
122  SrcIterator src_lr = src_ul + Diff2D(src_width, src_height);
123 
124  SrcIterator yys = xs;
125  KernelIterator yk = ki;
126 
127  // calculate width and height of the kernel
128  int kernel_width = klr.x - kul.x + 1;
129  int kernel_height = klr.y - kul.y + 1;
130 
131  // where the kernel is beyond the borders:
132  bool top_to_much = (y<klr.y) ? true : false;
133  bool down_to_much = (src_height-y-1<-kul.y)? true : false;
134  bool left_to_much = (x<klr.x)? true : false;
135  bool right_to_much = (src_width-x-1<-kul.x)? true : false;
136 
137  // direction of iteration,
138  // e.g. (-1, +1) for ll->ur or (-1, -1) for lr->ul
139  Diff2D way_increment;
140 
141  /* Iteration is always done from valid to invalid range.
142  The following tuple is composed as such:
143  - If an invalid range is reached while iterating in X,
144  a jump of border_increment.first is performed and
145  border_increment.third is used for further iterating.
146  - If an invalid range is reached while iterating in Y,
147  a jump of border_increment.second is performed and
148  border_increment.fourth is used for further iterating.
149  */
150  tuple4<int, int, int, int> border_increment;
151  if (border == BORDER_TREATMENT_REPEAT){
152  border_increment = tuple4<int, int, int, int>(1, 1, 0, 0);
153  }else if (border == BORDER_TREATMENT_REFLECT){
154  border_increment = tuple4<int, int, int, int>(2, 2, -1, -1);
155  }else{ // BORDER_TREATMENT_WRAP
156  border_increment = tuple4<int, int, int, int>(src_width, src_height, 1, 1);
157  }
158 
159  pair<int, int> valid_step_count;
160 
161  if(left_to_much && !top_to_much && !down_to_much)
162  {
163  yys += klr;
164  yk += kul;
165  way_increment = Diff2D(-1, -1);
166  border_increment.third = -border_increment.third;
167  border_increment.fourth = -border_increment.fourth;
168  valid_step_count = std::make_pair((yys - src_ul).x + 1, kernel_height);
169  }
170  else if(top_to_much && !left_to_much && !right_to_much)
171  {
172  yys += klr;
173  yk += kul;
174  way_increment = Diff2D(-1, -1);
175  border_increment.third = -border_increment.third;
176  border_increment.fourth = -border_increment.fourth;
177  valid_step_count = std::make_pair(kernel_width, (yys - src_ul).y + 1);
178  }
179  else if(right_to_much && !top_to_much && !down_to_much)
180  {
181  yys += kul;
182  yk += klr;
183  way_increment = Diff2D(1, 1);
184  border_increment.first = -border_increment.first;
185  border_increment.second = -border_increment.second;
186  valid_step_count = std::make_pair((src_lr - yys).x, kernel_height);
187  }
188  else if(down_to_much && !left_to_much && !right_to_much)
189  {
190  yys += kul;
191  yk += klr;
192  way_increment = Diff2D(1, 1);
193  border_increment.first = -border_increment.first;
194  border_increment.second = -border_increment.second;
195  valid_step_count = std::make_pair(kernel_width, (src_lr - yys).y);
196  }
197  else if(down_to_much && left_to_much)
198  {
199  yys += kul + Diff2D(kernel_width - 1, 0);
200  yk += kul + Diff2D(0, kernel_height - 1);
201  way_increment = Diff2D(-1, 1);
202  border_increment.second = -border_increment.second;
203  border_increment.third = -border_increment.third;
204  valid_step_count = std::make_pair((yys - src_ul).x + 1, (src_lr - yys).y);
205  }
206  else if(down_to_much && right_to_much)
207  {
208  yys += kul;
209  yk += klr;
210  way_increment = Diff2D(1, 1);
211  border_increment.first = -border_increment.first;
212  border_increment.second = -border_increment.second;
213  valid_step_count = std::make_pair((src_lr - yys).x, (src_lr - yys).y);
214  }
215  else if(top_to_much && left_to_much)
216  {
217  yys += klr;
218  yk += kul;
219  way_increment = Diff2D(-1, -1);
220  border_increment.third = -border_increment.third;
221  border_increment.fourth = -border_increment.fourth;
222  valid_step_count = std::make_pair((yys - src_ul).x + 1, (yys - src_ul).y + 1);
223  }
224  else
225  { //top_to_much && right_to_much
226  yys += kul + Diff2D(0, kernel_height - 1);
227  yk += kul + Diff2D(kernel_width - 1, 0);
228  way_increment = Diff2D(1, -1);
229  border_increment.first = -border_increment.first;
230  border_increment.fourth = -border_increment.fourth;
231  valid_step_count = std::make_pair((src_lr - yys).x, (yys - src_ul).y + 1);
232  }
233 
234  int yy = 0, xx;
235 
236  //laeuft den zulässigen Bereich in y-Richtung durch
237  for(; yy < valid_step_count.second; ++yy, yys.y += way_increment.y, yk.y -= way_increment.y )
238  {
239  SrcIterator xxs = yys;
240  KernelIterator xk = yk;
241 
242  //laeuft den zulässigen Bereich in x-Richtung durch
243  for(xx = 0; xx < valid_step_count.first; ++xx, xxs.x += way_increment.x, xk.x -= way_increment.x)
244  {
245  sum += ak(xk) * src_acc(xxs);
246  }
247 
248  //Nächstes ++xxs.x wuerde in unzulässigen Bereich
249  //bringen => Sprung in zulaessigen Bereich
250  xxs.x += border_increment.first;
251 
252  for( ; xx < kernel_width; ++xx, xxs.x += border_increment.third, xk.x -= way_increment.x )
253  {
254  sum += ak(xk) * src_acc(xxs);
255  }
256  }
257 
258  //Nächstes ++yys.y wuerde in unzulässigen Bereich
259  //bringen => Sprung in zulaessigen Bereich
260  yys.y += border_increment.second;
261 
262  for( ; yy < kernel_height; ++yy, yys.y += border_increment.third, yk.y -= way_increment.y)
263  {
264  SrcIterator xxs = yys;
265  KernelIterator xk = yk;
266 
267  for(xx=0; xx < valid_step_count.first; ++xx, xxs.x += way_increment.x, xk.x -= way_increment.x)
268  {
269  sum += ak(xk) * src_acc(xxs);
270  }
271 
272  //Sprung in den zulaessigen Bereich
273  xxs.x += border_increment.first;
274 
275  for( ; xx < kernel_width; ++xx, xxs.x += border_increment.third, xk.x -= way_increment.x )
276  {
277  sum += ak(xk) * src_acc(xxs);
278  }
279  }
280 
281  // store average in destination pixel
282  dest_acc.set(DestTraits::fromRealPromote(sum), xd);
283 
284 }// end of internalPixelEvaluationByWrapReflectRepeat
285 #endif /* #if 0 */
286 
287 
288 template <class SrcIterator, class SrcAccessor,
289  class KernelIterator, class KernelAccessor,
290  class SumType>
291 void
292 internalPixelEvaluationByWrapReflectRepeat(SrcIterator xs, SrcAccessor src_acc,
293  KernelIterator xk, KernelAccessor ak,
294  int left, int right, int kleft, int kright,
295  int borderskipx, int borderinc, SumType & sum)
296 {
297  SrcIterator xxs = xs + left;
298  KernelIterator xxk = xk - left;
299 
300  for(int xx = left; xx <= right; ++xx, ++xxs, --xxk)
301  {
302  sum += ak(xxk) * src_acc(xxs);
303  }
304 
305  xxs = xs + left - borderskipx;
306  xxk = xk - left + 1;
307  for(int xx = left - 1; xx >= -kright; --xx, xxs -= borderinc, ++xxk)
308  {
309  sum += ak(xxk) * src_acc(xxs);
310  }
311 
312  xxs = xs + right + borderskipx;
313  xxk = xk - right - 1;
314  for(int xx = right + 1; xx <= -kleft; ++xx, xxs += borderinc, --xxk)
315  {
316  sum += ak(xxk) * src_acc(xxs);
317  }
318 }
319 
320 
321 /** \addtogroup StandardConvolution Two-dimensional convolution functions
322 
323 Perform 2D non-separable convolution, with and without ROI mask.
324 
325 These generic convolution functions implement
326 the standard 2D convolution operation for images that fit
327 into the required interface. Arbitrary ROI's are supported
328 by the mask version of the algorithm.
329 The functions need a suitable 2D kernel to operate.
330 */
331 //@{
332 
333 /** \brief Performs a 2 dimensional convolution of the source image using the given
334  kernel.
335 
336  The KernelIterator must point to the center of the kernel, and
337  the kernel's size is given by its upper left (x and y of distance <= 0) and
338  lower right (distance >= 0) corners. The image must always be larger than the
339  kernel. At those positions where the kernel does not completely fit
340  into the image, the specified \ref BorderTreatmentMode is
341  applied. You can choice between following BorderTreatmentModes:
342  <ul>
343  <li>BORDER_TREATMENT_CLIP</li>
344  <li>BORDER_TREATMENT_AVOID</li>
345  <li>BORDER_TREATMENT_WRAP</li>
346  <li>BORDER_TREATMENT_REFLECT</li>
347  <li>BORDER_TREATMENT_REPEAT</li>
348  </ul><br>
349  The images's pixel type (SrcAccessor::value_type) must be a
350  linear space over the kernel's value_type (KernelAccessor::value_type),
351  i.e. addition of source values, multiplication with kernel values,
352  and NumericTraits must be defined.
353  The kernel's value_type must be an algebraic field,
354  i.e. the arithmetic operations (+, -, *, /) and NumericTraits must
355  be defined.
356 
357  <b> Declarations:</b>
358 
359  pass arguments explicitly:
360  \code
361  namespace vigra {
362  template <class SrcIterator, class SrcAccessor,
363  class DestIterator, class DestAccessor,
364  class KernelIterator, class KernelAccessor>
365  void convolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
366  DestIterator dest_ul, DestAccessor dest_acc,
367  KernelIterator ki, KernelAccessor ak,
368  Diff2D kul, Diff2D klr, BorderTreatmentMode border);
369  }
370  \endcode
371 
372 
373  use argument objects in conjunction with \ref ArgumentObjectFactories :
374  \code
375  namespace vigra {
376  template <class SrcIterator, class SrcAccessor,
377  class DestIterator, class DestAccessor,
378  class KernelIterator, class KernelAccessor>
379  void convolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
380  pair<DestIterator, DestAccessor> dest,
381  tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
382  BorderTreatmentMode> kernel);
383  }
384  \endcode
385 
386  <b> Usage:</b>
387 
388  <b>\#include</b> <<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>><br>
389  Namespace: vigra
390 
391 
392  \code
393  vigra::FImage src(w,h), dest(w,h);
394  ...
395 
396  // define horizontal Sobel filter
397  vigra::Kernel2D<float> sobel;
398 
399  sobel.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) = // upper left and lower right
400  0.125, 0.0, -0.125,
401  0.25, 0.0, -0.25,
402  0.125, 0.0, -0.125;
403  sobel.setBorderTreatment(vigra::BORDER_TREATMENT_REFLECT);
404 
405  vigra::convolveImage(srcImageRange(src), destImage(dest), kernel2d(sobel));
406  \endcode
407 
408  <b> Required Interface:</b>
409 
410  \code
411  ImageIterator src_ul, src_lr;
412  ImageIterator dest_ul;
413  ImageIterator ik;
414 
415  SrcAccessor src_accessor;
416  DestAccessor dest_accessor;
417  KernelAccessor kernel_accessor;
418 
419  NumericTraits<SrcAccessor::value_type>::RealPromote s = src_accessor(src_ul);
420 
421  s = s + s;
422  s = kernel_accessor(ik) * s;
423  s -= s;
424 
425  dest_accessor.set(
426  NumericTraits<DestAccessor::value_type>::fromRealPromote(s), dest_ul);
427 
428  NumericTraits<KernelAccessor::value_type>::RealPromote k = kernel_accessor(ik);
429 
430  k += k;
431  k -= k;
432  k = k / k;
433 
434  \endcode
435 
436  <b> Preconditions:</b>
437 
438  \code
439  kul.x <= 0
440  kul.y <= 0
441  klr.x >= 0
442  klr.y >= 0
443  src_lr.x - src_ul.x >= klr.x + kul.x + 1
444  src_lr.y - src_ul.y >= klr.y + kul.y + 1
445  \endcode
446 
447  If border == BORDER_TREATMENT_CLIP: Sum of kernel elements must be
448  != 0.
449 
450 */
451 template <class SrcIterator, class SrcAccessor,
452  class DestIterator, class DestAccessor,
453  class KernelIterator, class KernelAccessor>
454 void convolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
455  DestIterator dest_ul, DestAccessor dest_acc,
456  KernelIterator ki, KernelAccessor ak,
457  Diff2D kul, Diff2D klr, BorderTreatmentMode border)
458 {
459  vigra_precondition((border == BORDER_TREATMENT_CLIP ||
460  border == BORDER_TREATMENT_AVOID ||
461  border == BORDER_TREATMENT_REFLECT ||
462  border == BORDER_TREATMENT_REPEAT ||
463  border == BORDER_TREATMENT_WRAP),
464  "convolveImage():\n"
465  " Border treatment must be one of follow treatments:\n"
466  " - BORDER_TREATMENT_CLIP\n"
467  " - BORDER_TREATMENT_AVOID\n"
468  " - BORDER_TREATMENT_REFLECT\n"
469  " - BORDER_TREATMENT_REPEAT\n"
470  " - BORDER_TREATMENT_WRAP\n");
471 
472  vigra_precondition(kul.x <= 0 && kul.y <= 0,
473  "convolveImage(): coordinates of "
474  "kernel's upper left must be <= 0.");
475  vigra_precondition(klr.x >= 0 && klr.y >= 0,
476  "convolveImage(): coordinates of "
477  "kernel's lower right must be >= 0.");
478 
479  // use traits to determine SumType as to prevent possible overflow
480  typedef typename
481  PromoteTraits<typename SrcAccessor::value_type,
482  typename KernelAccessor::value_type>::Promote SumType;
483  typedef typename
484  NumericTraits<typename KernelAccessor::value_type>::RealPromote KernelSumType;
485  typedef typename DestAccessor::value_type DestType;
486 
487  // calculate width and height of the image
488  int w = src_lr.x - src_ul.x;
489  int h = src_lr.y - src_ul.y;
490 
491  // calculate width and height of the kernel
492  int kernel_width = klr.x - kul.x + 1;
493  int kernel_height = klr.y - kul.y + 1;
494 
495  vigra_precondition(w >= kernel_width && h >= kernel_height,
496  "convolveImage(): kernel larger than image.");
497 
498  KernelSumType norm = NumericTraits<KernelSumType>::zero();
499  if(border == BORDER_TREATMENT_CLIP)
500  {
501  // calculate the sum of the kernel elements for renormalization
502  KernelIterator yk = ki + klr;
503 
504  // determine sum within kernel (= norm)
505  for(int y = 0; y < kernel_height; ++y, --yk.y)
506  {
507  KernelIterator xk = yk;
508  for(int x = 0; x < kernel_width; ++x, --xk.x)
509  {
510  norm += ak(xk);
511  }
512  }
513  vigra_precondition(norm != NumericTraits<KernelSumType>::zero(),
514  "convolveImage(): Cannot use BORDER_TREATMENT_CLIP with a DC-free kernel");
515  }
516 
517  // create iterators for the interior part of the image (where the kernel always fits into the image)
518  DestIterator yd = dest_ul + Diff2D(klr.x, klr.y);
519  SrcIterator ys = src_ul + Diff2D(klr.x, klr.y);
520  SrcIterator send = src_lr + Diff2D(kul.x, kul.y);
521 
522  // iterate over the interior part
523  for(; ys.y < send.y; ++ys.y, ++yd.y)
524  {
525  // create x iterators
526  DestIterator xd(yd);
527  SrcIterator xs(ys);
528 
529  for(; xs.x < send.x; ++xs.x, ++xd.x)
530  {
531  // init the sum
532  SumType sum = NumericTraits<SumType>::zero();
533 
534  SrcIterator yys = xs - klr;
535  SrcIterator yyend = xs - kul;
536  KernelIterator yk = ki + klr;
537 
538  for(; yys.y <= yyend.y; ++yys.y, --yk.y)
539  {
540  typename SrcIterator::row_iterator xxs = yys.rowIterator();
541  typename SrcIterator::row_iterator xxe = xxs + kernel_width;
542  typename KernelIterator::row_iterator xk = yk.rowIterator();
543 
544  for(; xxs < xxe; ++xxs, --xk)
545  {
546  sum += ak(xk) * src_acc(xxs);
547  }
548  }
549 
550  // store convolution result in destination pixel
551  dest_acc.set(detail::RequiresExplicitCast<DestType>::cast(sum), xd);
552  }
553  }
554 
555  if(border == BORDER_TREATMENT_AVOID)
556  return; // skip processing near the border
557 
558  int interiorskip = w + kul.x - klr.x - 1;
559  int borderskipx = 0;
560  int borderskipy = 0;
561  int borderinc = 0;
562  if(border == BORDER_TREATMENT_REPEAT)
563  {
564  borderskipx = 0;
565  borderskipy = 0;
566  borderinc = 0;
567  }
568  else if(border == BORDER_TREATMENT_REFLECT)
569  {
570  borderskipx = -1;
571  borderskipy = -1;
572  borderinc = -1;
573  }
574  else if(border == BORDER_TREATMENT_WRAP)
575  {
576  borderskipx = -w+1;
577  borderskipy = -h+1;
578  borderinc = 1;
579  }
580 
581  // create iterators for the entire image
582  yd = dest_ul;
583  ys = src_ul;
584 
585  // work on entire image (but skip the already computed points in the loop)
586  for(int y = 0; y < h; ++y, ++ys.y, ++yd.y)
587  {
588  int top = int(std::max(static_cast<IntBiggest>(-klr.y),
589  static_cast<IntBiggest>(src_ul.y - ys.y)));
590  int bottom = int(std::min(static_cast<IntBiggest>(-kul.y),
591  static_cast<IntBiggest>(src_lr.y - ys.y - 1)));
592 
593  // create x iterators
594  DestIterator xd(yd);
595  SrcIterator xs(ys);
596 
597  for(int x = 0; x < w; ++x, ++xs.x, ++xd.x)
598  {
599  // check if we are away from the border
600  if(y >= klr.y && y < h+kul.y && x == klr.x)
601  {
602  // yes => skip the already computed points
603  x += interiorskip;
604  xs.x += interiorskip;
605  xd.x += interiorskip;
606  continue;
607  }
608  if (border == BORDER_TREATMENT_CLIP)
609  {
610  internalPixelEvaluationByClip(x, y, w, h, xs, src_acc, xd, dest_acc, ki, kul, klr, ak, norm);
611  }
612  else
613  {
614  int left = std::max(-klr.x, src_ul.x - xs.x);
615  int right = std::min(-kul.x, src_lr.x - xs.x - 1);
616 
617  // init the sum
618  SumType sum = NumericTraits<SumType>::zero();
619 
620  // create iterators for the part of the kernel that fits into the image
621  SrcIterator yys = xs + Size2D(0, top);
622  KernelIterator yk = ki - Size2D(0, top);
623 
624  int yy;
625  for(yy = top; yy <= bottom; ++yy, ++yys.y, --yk.y)
626  {
627  internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak,
628  left, right, kul.x, klr.x, borderskipx, borderinc, sum);
629  }
630  yys = xs + Size2D(0, top - borderskipy);
631  yk = ki - Size2D(0, top - 1);
632  for(yy = top - 1; yy >= -klr.y; --yy, yys.y -= borderinc, ++yk.y)
633  {
634  internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak,
635  left, right, kul.x, klr.x, borderskipx, borderinc, sum);
636  }
637  yys = xs + Size2D(0, bottom + borderskipy);
638  yk = ki - Size2D(0, bottom + 1);
639  for(yy = bottom + 1; yy <= -kul.y; ++yy, yys.y += borderinc, --yk.y)
640  {
641  internalPixelEvaluationByWrapReflectRepeat(yys.rowIterator(), src_acc, yk.rowIterator(), ak,
642  left, right, kul.x, klr.x, borderskipx, borderinc, sum);
643  }
644 
645  // store convolution result in destination pixel
646  dest_acc.set(detail::RequiresExplicitCast<DestType>::cast(sum), xd);
647 
648 // internalPixelEvaluationByWrapReflectRepeat(x, y, w, h, xs, src_acc, xd, dest_acc, ki, kul, klr, ak, border);
649  }
650  }
651  }
652 }
653 
654 
655 template <class SrcIterator, class SrcAccessor,
656  class DestIterator, class DestAccessor,
657  class KernelIterator, class KernelAccessor>
658 inline
659 void convolveImage(
660  triple<SrcIterator, SrcIterator, SrcAccessor> src,
661  pair<DestIterator, DestAccessor> dest,
662  tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
663  BorderTreatmentMode> kernel)
664 {
665  convolveImage(src.first, src.second, src.third,
666  dest.first, dest.second,
667  kernel.first, kernel.second, kernel.third,
668  kernel.fourth, kernel.fifth);
669 }
670 
671 
672 /** \brief Performs a 2-dimensional normalized convolution, i.e. convolution with a mask image.
673 
674  This functions computes
675  <a href ="http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/PIRODDI1/NormConv/NormConv.html">normalized
676  convolution</a> as defined in
677  Knutsson, H. and Westin, C-F.: <i>Normalized and differential convolution:
678  Methods for Interpolation and Filtering of incomplete and uncertain data</i>.
679  Proc. of the IEEE Conf. on Computer Vision and Pattern Recognition, 1993, 515-523.
680 
681  The mask image must be binary and encodes which pixels of the original image
682  are valid. It is used as follows:
683  Only pixel under the mask are used in the calculations. Whenever a part of the
684  kernel lies outside the mask, it is ignored, and the kernel is renormalized to its
685  original norm (analogous to the CLIP \ref BorderTreatmentMode). Thus, a useful convolution
686  result is computed whenever <i>at least one valid pixel is within the current window</i>
687  Thus, destination pixels not under the mask still receive a value if they are <i>near</i>
688  the mask. Therefore, this algorithm is useful as an interpolator of sparse input data.
689  If you are only interested in the destination values under the mask, you can perform
690  a subsequent \ref copyImageIf().
691 
692  The KernelIterator must point to the center of the kernel, and
693  the kernel's size is given by its upper left (x and y of distance <= 0) and
694  lower right (distance >= 0) corners. The image must always be larger than the
695  kernel. At those positions where the kernel does not completely fit
696  into the image, the specified \ref BorderTreatmentMode is
697  applied. Only BORDER_TREATMENT_CLIP and BORDER_TREATMENT_AVOID are currently
698  supported.
699 
700  The images's pixel type (SrcAccessor::value_type) must be a
701  linear space over the kernel's value_type (KernelAccessor::value_type),
702  i.e. addition of source values, multiplication with kernel values,
703  and NumericTraits must be defined.
704  The kernel's value_type must be an algebraic field,
705  i.e. the arithmetic operations (+, -, *, /) and NumericTraits must
706  be defined.
707 
708  <b> Declarations:</b>
709 
710  pass arguments explicitly:
711  \code
712  namespace vigra {
713  template <class SrcIterator, class SrcAccessor,
714  class MaskIterator, class MaskAccessor,
715  class DestIterator, class DestAccessor,
716  class KernelIterator, class KernelAccessor>
717  void
718  normalizedConvolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
719  MaskIterator mul, MaskAccessor am,
720  DestIterator dest_ul, DestAccessor dest_acc,
721  KernelIterator ki, KernelAccessor ak,
722  Diff2D kul, Diff2D klr, BorderTreatmentMode border);
723  }
724  \endcode
725 
726 
727  use argument objects in conjunction with \ref ArgumentObjectFactories :
728  \code
729  namespace vigra {
730  template <class SrcIterator, class SrcAccessor,
731  class MaskIterator, class MaskAccessor,
732  class DestIterator, class DestAccessor,
733  class KernelIterator, class KernelAccessor>
734  void normalizedConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
735  pair<MaskIterator, MaskAccessor> mask,
736  pair<DestIterator, DestAccessor> dest,
737  tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
738  BorderTreatmentMode> kernel);
739  }
740  \endcode
741 
742  <b> Usage:</b>
743 
744  <b>\#include</b> <<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>><br>
745  Namespace: vigra
746 
747 
748  \code
749  vigra::FImage src(w,h), dest(w,h);
750  vigra::CImage mask(w,h);
751  ...
752 
753  // define 3x3 binomial filter
754  vigra::Kernel2D<float> binom;
755 
756  binom.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) = // upper left and lower right
757  0.0625, 0.125, 0.0625,
758  0.125, 0.25, 0.125,
759  0.0625, 0.125, 0.0625;
760 
761  vigra::normalizedConvolveImage(srcImageRange(src), maskImage(mask), destImage(dest), kernel2d(binom));
762  \endcode
763 
764  <b> Required Interface:</b>
765 
766  \code
767  ImageIterator src_ul, src_lr;
768  ImageIterator mul;
769  ImageIterator dest_ul;
770  ImageIterator ik;
771 
772  SrcAccessor src_accessor;
773  MaskAccessor mask_accessor;
774  DestAccessor dest_accessor;
775  KernelAccessor kernel_accessor;
776 
777  NumericTraits<SrcAccessor::value_type>::RealPromote s = src_accessor(src_ul);
778 
779  s = s + s;
780  s = kernel_accessor(ik) * s;
781  s -= s;
782 
783  if(mask_accessor(mul)) ...;
784 
785  dest_accessor.set(
786  NumericTraits<DestAccessor::value_type>::fromRealPromote(s), dest_ul);
787 
788  NumericTraits<KernelAccessor::value_type>::RealPromote k = kernel_accessor(ik);
789 
790  k += k;
791  k -= k;
792  k = k / k;
793 
794  \endcode
795 
796  <b> Preconditions:</b>
797 
798  \code
799  kul.x <= 0
800  kul.y <= 0
801  klr.x >= 0
802  klr.y >= 0
803  src_lr.x - src_ul.x >= klr.x + kul.x + 1
804  src_lr.y - src_ul.y >= klr.y + kul.y + 1
805  border == BORDER_TREATMENT_CLIP || border == BORDER_TREATMENT_AVOID
806  \endcode
807 
808  Sum of kernel elements must be != 0.
809 
810 */
812 
813 template <class SrcIterator, class SrcAccessor,
814  class DestIterator, class DestAccessor,
815  class MaskIterator, class MaskAccessor,
816  class KernelIterator, class KernelAccessor>
817 void
818 normalizedConvolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
819  MaskIterator mul, MaskAccessor am,
820  DestIterator dest_ul, DestAccessor dest_acc,
821  KernelIterator ki, KernelAccessor ak,
822  Diff2D kul, Diff2D klr, BorderTreatmentMode border)
823 {
824  vigra_precondition((border == BORDER_TREATMENT_CLIP ||
825  border == BORDER_TREATMENT_AVOID),
826  "normalizedConvolveImage(): "
827  "Border treatment must be BORDER_TREATMENT_CLIP or BORDER_TREATMENT_AVOID.");
828 
829  vigra_precondition(kul.x <= 0 && kul.y <= 0,
830  "normalizedConvolveImage(): left borders must be <= 0.");
831  vigra_precondition(klr.x >= 0 && klr.y >= 0,
832  "normalizedConvolveImage(): right borders must be >= 0.");
833 
834  // use traits to determine SumType as to prevent possible overflow
835  typedef typename
836  NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType;
837  typedef typename
838  NumericTraits<typename KernelAccessor::value_type>::RealPromote KSumType;
839  typedef
840  NumericTraits<typename DestAccessor::value_type> DestTraits;
841 
842  // calculate width and height of the image
843  int w = src_lr.x - src_ul.x;
844  int h = src_lr.y - src_ul.y;
845  int kernel_width = klr.x - kul.x + 1;
846  int kernel_height = klr.y - kul.y + 1;
847 
848  int x,y;
849  int ystart = (border == BORDER_TREATMENT_AVOID) ? klr.y : 0;
850  int yend = (border == BORDER_TREATMENT_AVOID) ? h+kul.y : h;
851  int xstart = (border == BORDER_TREATMENT_AVOID) ? klr.x : 0;
852  int xend = (border == BORDER_TREATMENT_AVOID) ? w+kul.x : w;
853 
854  // create y iterators
855  DestIterator yd = dest_ul + Diff2D(xstart, ystart);
856  SrcIterator ys = src_ul + Diff2D(xstart, ystart);
857  MaskIterator ym = mul + Diff2D(xstart, ystart);
858 
859  KSumType norm = ak(ki);
860  int xx, yy;
861  KernelIterator yk = ki + klr;
862  for(yy=0; yy<kernel_height; ++yy, --yk.y)
863  {
864  KernelIterator xk = yk;
865 
866  for(xx=0; xx<kernel_width; ++xx, --xk.x)
867  {
868  norm += ak(xk);
869  }
870  }
871  norm -= ak(ki);
872 
873 
874  for(y=ystart; y < yend; ++y, ++ys.y, ++yd.y, ++ym.y)
875  {
876  // create x iterators
877  DestIterator xd(yd);
878  SrcIterator xs(ys);
879  MaskIterator xm(ym);
880 
881  for(x=xstart; x < xend; ++x, ++xs.x, ++xd.x, ++xm.x)
882  {
883  // how much of the kernel fits into the image ?
884  int x0, y0, x1, y1;
885 
886  y0 = (y<klr.y) ? -y : -klr.y;
887  y1 = (h-y-1<-kul.y) ? h-y-1 : -kul.y;
888  x0 = (x<klr.x) ? -x : -klr.x;
889  x1 = (w-x-1<-kul.x) ? w-x-1 : -kul.x;
890 
891  bool first = true;
892  // init the sum
893  SumType sum;
894  KSumType ksum;
895 
896  SrcIterator yys = xs + Diff2D(x0, y0);
897  MaskIterator yym = xm + Diff2D(x0, y0);
898  KernelIterator yk = ki - Diff2D(x0, y0);
899 
900  int xx, kernel_width, kernel_height;
901  kernel_width = x1 - x0 + 1;
902  kernel_height = y1 - y0 + 1;
903  for(yy=0; yy<kernel_height; ++yy, ++yys.y, --yk.y, ++yym.y)
904  {
905  typename SrcIterator::row_iterator xxs = yys.rowIterator();
906  typename SrcIterator::row_iterator xxend = xxs + kernel_width;
907  typename MaskIterator::row_iterator xxm = yym.rowIterator();
908  typename KernelIterator::row_iterator xk = yk.rowIterator();
909 
910  for(xx=0; xxs < xxend; ++xxs, --xk, ++xxm)
911  {
912  if(!am(xxm)) continue;
913 
914  if(first)
915  {
916  sum = detail::RequiresExplicitCast<SumType>::cast(ak(xk) * src_acc(xxs));
917  ksum = ak(xk);
918  first = false;
919  }
920  else
921  {
922  sum = detail::RequiresExplicitCast<SumType>::cast(sum + ak(xk) * src_acc(xxs));
923  ksum += ak(xk);
924  }
925  }
926  }
927  // store average in destination pixel
928  if(!first &&
929  ksum != NumericTraits<KSumType>::zero())
930  {
931  dest_acc.set(DestTraits::fromRealPromote(
932  detail::RequiresExplicitCast<SumType>::cast((norm / ksum) * sum)), xd);
933  }
934  }
935  }
936 }
937 
938 
939 template <class SrcIterator, class SrcAccessor,
940  class DestIterator, class DestAccessor,
941  class MaskIterator, class MaskAccessor,
942  class KernelIterator, class KernelAccessor>
943 inline
945  triple<SrcIterator, SrcIterator, SrcAccessor> src,
946  pair<MaskIterator, MaskAccessor> mask,
947  pair<DestIterator, DestAccessor> dest,
948  tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
949  BorderTreatmentMode> kernel)
950 {
951  normalizedConvolveImage(src.first, src.second, src.third,
952  mask.first, mask.second,
953  dest.first, dest.second,
954  kernel.first, kernel.second, kernel.third,
955  kernel.fourth, kernel.fifth);
956 }
957 
958 /** \brief Deprecated name of 2-dimensional normalized convolution, i.e. convolution with a mask image.
959 
960  See \ref normalizedConvolveImage() for documentation.
961 
962  <b> Declarations:</b>
963 
964  pass arguments explicitly:
965  \code
966  namespace vigra {
967  template <class SrcIterator, class SrcAccessor,
968  class MaskIterator, class MaskAccessor,
969  class DestIterator, class DestAccessor,
970  class KernelIterator, class KernelAccessor>
971  void
972  convolveImageWithMask(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
973  MaskIterator mul, MaskAccessor am,
974  DestIterator dest_ul, DestAccessor dest_acc,
975  KernelIterator ki, KernelAccessor ak,
976  Diff2D kul, Diff2D klr, BorderTreatmentMode border);
977  }
978  \endcode
979 
980 
981  use argument objects in conjunction with \ref ArgumentObjectFactories :
982  \code
983  namespace vigra {
984  template <class SrcIterator, class SrcAccessor,
985  class MaskIterator, class MaskAccessor,
986  class DestIterator, class DestAccessor,
987  class KernelIterator, class KernelAccessor>
988  void convolveImageWithMask(triple<SrcIterator, SrcIterator, SrcAccessor> src,
989  pair<MaskIterator, MaskAccessor> mask,
990  pair<DestIterator, DestAccessor> dest,
991  tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
992  BorderTreatmentMode> kernel);
993  }
994  \endcode
995 */
997 
998 template <class SrcIterator, class SrcAccessor,
999  class DestIterator, class DestAccessor,
1000  class MaskIterator, class MaskAccessor,
1001  class KernelIterator, class KernelAccessor>
1002 inline void
1003 convolveImageWithMask(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
1004  MaskIterator mul, MaskAccessor am,
1005  DestIterator dest_ul, DestAccessor dest_acc,
1006  KernelIterator ki, KernelAccessor ak,
1007  Diff2D kul, Diff2D klr, BorderTreatmentMode border)
1008 {
1009  normalizedConvolveImage(src_ul, src_lr, src_acc,
1010  mul, am,
1011  dest_ul, dest_acc,
1012  ki, ak, kul, klr, border);
1013 }
1014 
1015 template <class SrcIterator, class SrcAccessor,
1016  class DestIterator, class DestAccessor,
1017  class MaskIterator, class MaskAccessor,
1018  class KernelIterator, class KernelAccessor>
1019 inline
1021  triple<SrcIterator, SrcIterator, SrcAccessor> src,
1022  pair<MaskIterator, MaskAccessor> mask,
1023  pair<DestIterator, DestAccessor> dest,
1024  tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
1025  BorderTreatmentMode> kernel)
1026 {
1027  normalizedConvolveImage(src.first, src.second, src.third,
1028  mask.first, mask.second,
1029  dest.first, dest.second,
1030  kernel.first, kernel.second, kernel.third,
1031  kernel.fourth, kernel.fifth);
1032 }
1033 
1034 //@}
1035 
1036 /********************************************************/
1037 /* */
1038 /* Kernel2D */
1039 /* */
1040 /********************************************************/
1041 
1042 /** \brief Generic 2 dimensional convolution kernel.
1043 
1044  This kernel may be used for convolution of 2 dimensional signals.
1045 
1046  Convolution functions access the kernel via an ImageIterator
1047  which they get by calling \ref center(). This iterator
1048  points to the center of the kernel. The kernel's size is given by its upperLeft()
1049  (upperLeft().x <= 0, upperLeft().y <= 0)
1050  and lowerRight() (lowerRight().x >= 0, lowerRight().y >= 0) methods.
1051  The desired border treatment mode is returned by borderTreatment().
1052  (Note that the \ref StandardConvolution "2D convolution functions" don't currently
1053  support all modes.)
1054 
1055  The different init functions create a kernel with the specified
1056  properties. The requirements for the kernel's value_type depend
1057  on the init function used. At least NumericTraits must be defined.
1058 
1059  The kernel defines a factory function kernel2d() to create an argument object
1060  (see \ref KernelArgumentObjectFactories).
1061 
1062  <b> Usage:</b>
1063 
1064  <b>\#include</b> <<a href="stdconvolution_8hxx-source.html">vigra/stdconvolution.hxx</a>><br>
1065  Namespace: vigra
1066 
1067  \code
1068  vigra::FImage src(w,h), dest(w,h);
1069  ...
1070 
1071  // define horizontal Sobel filter
1072  vigra::Kernel2D<float> sobel;
1073 
1074  sobel.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) = // upper left and lower right
1075  0.125, 0.0, -0.125,
1076  0.25, 0.0, -0.25,
1077  0.125, 0.0, -0.125;
1078 
1079  vigra::convolveImage(srcImageRange(src), destImage(dest), kernel2d(sobel));
1080  \endcode
1081 
1082  <b> Required Interface:</b>
1083 
1084  \code
1085  value_type v = NumericTraits<value_type>::one();
1086  \endcode
1087 
1088  See also the init functions.
1089 */
1090 template <class ARITHTYPE>
1092 {
1093 public:
1094  /** the kernel's value type
1095  */
1096  typedef ARITHTYPE value_type;
1097 
1098  /** 2D random access iterator over the kernel's values
1099  */
1101 
1102  /** const 2D random access iterator over the kernel's values
1103  */
1105 
1106  /** the kernel's accessor
1107  */
1109 
1110  /** the kernel's const accessor
1111  */
1113 
1114  struct InitProxy
1115  {
1116  typedef typename
1118 
1119  InitProxy(Iterator i, int count, value_type & norm)
1120  : iter_(i), base_(i),
1121  count_(count), sum_(count),
1122  norm_(norm)
1123  {}
1124 
1125  ~InitProxy()
1126  {
1127  vigra_precondition(count_ == 1 || count_ == sum_,
1128  "Kernel2D::initExplicitly(): "
1129  "Too few init values.");
1130  }
1131 
1132  InitProxy & operator,(value_type const & v)
1133  {
1134  if(count_ == sum_) norm_ = *iter_;
1135 
1136  --count_;
1137  vigra_precondition(count_ > 0,
1138  "Kernel2D::initExplicitly(): "
1139  "Too many init values.");
1140 
1141  norm_ += v;
1142 
1143  ++iter_;
1144  *iter_ = v;
1145 
1146  return *this;
1147  }
1148 
1149  Iterator iter_, base_;
1150  int count_, sum_;
1151  value_type & norm_;
1152  };
1153 
1154  static value_type one() { return NumericTraits<value_type>::one(); }
1155 
1156  /** Default constructor.
1157  Creates a kernel of size 1x1 which would copy the signal
1158  unchanged.
1159  */
1161  : kernel_(1, 1, one()),
1162  left_(0, 0),
1163  right_(0, 0),
1164  norm_(one()),
1165  border_treatment_(BORDER_TREATMENT_CLIP)
1166  {}
1167 
1168  /** Copy constructor.
1169  */
1170  Kernel2D(Kernel2D const & k)
1171  : kernel_(k.kernel_),
1172  left_(k.left_),
1173  right_(k.right_),
1174  norm_(k.norm_),
1175  border_treatment_(k.border_treatment_)
1176  {}
1177 
1178  /** Copy assignment.
1179  */
1181  {
1182  if(this != &k)
1183  {
1184  kernel_ = k.kernel_;
1185  left_ = k.left_;
1186  right_ = k.right_;
1187  norm_ = k.norm_;
1188  border_treatment_ = k.border_treatment_;
1189  }
1190  return *this;
1191  }
1192 
1193  /** Initialization.
1194  This initializes the kernel with the given constant. The norm becomes
1195  v*width()*height().
1196 
1197  Instead of a single value an initializer list of length width()*height()
1198  can be used like this:
1199 
1200  \code
1201  vigra::Kernel2D<float> binom;
1202 
1203  binom.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) =
1204  0.0625, 0.125, 0.0625,
1205  0.125, 0.25, 0.125,
1206  0.0625, 0.125, 0.0625;
1207  \endcode
1208 
1209  In this case, the norm will be set to the sum of the init values.
1210  An initializer list of wrong length will result in a run-time error.
1211  */
1212  InitProxy operator=(value_type const & v)
1213  {
1214  int size = (right_.x - left_.x + 1) *
1215  (right_.y - left_.y + 1);
1216  kernel_ = v;
1217  norm_ = (double)size*v;
1218 
1219  return InitProxy(kernel_.begin(), size, norm_);
1220  }
1221 
1222  /** Destructor.
1223  */
1225  {}
1226 
1227  /** Init the 2D kernel as the cartesian product of two 1D kernels
1228  of type \ref Kernel1D. The norm becomes the product of the two original
1229  norms.
1230 
1231  <b> Required Interface:</b>
1232 
1233  The kernel's value_type must be a linear algebra.
1234 
1235  \code
1236  vigra::Kernel2D<...>::value_type v;
1237  v = v * v;
1238  \endcode
1239  */
1241  Kernel1D<value_type> const & ky)
1242  {
1243  left_ = Diff2D(kx.left(), ky.left());
1244  right_ = Diff2D(kx.right(), ky.right());
1245  int w = right_.x - left_.x + 1;
1246  int h = right_.y - left_.y + 1;
1247  kernel_.resize(w, h);
1248 
1249  norm_ = kx.norm() * ky.norm();
1250 
1251  typedef typename Kernel1D<value_type>::const_iterator KIter;
1252  typename Kernel1D<value_type>::Accessor ka;
1253 
1254  KIter kiy = ky.center() + left_.y;
1255  Iterator iy = center() + left_;
1256 
1257  for(int y=left_.y; y<=right_.y; ++y, ++kiy, ++iy.y)
1258  {
1259  KIter kix = kx.center() + left_.x;
1260  Iterator ix = iy;
1261  for(int x=left_.x; x<=right_.x; ++x, ++kix, ++ix.x)
1262  {
1263  *ix = ka(kix) * ka(kiy);
1264  }
1265  }
1266  }
1267 
1268  /** Init the 2D kernel as the cartesian product of two 1D kernels
1269  given explicitly by iterators and sizes. The norm becomes the
1270  sum of the resulting kernel values.
1271 
1272  <b> Required Interface:</b>
1273 
1274  The kernel's value_type must be a linear algebra.
1275 
1276  \code
1277  vigra::Kernel2D<...>::value_type v;
1278  v = v * v;
1279  v += v;
1280  \endcode
1281 
1282  <b> Preconditions:</b>
1283 
1284  \code
1285  xleft <= 0;
1286  xright >= 0;
1287  yleft <= 0;
1288  yright >= 0;
1289  \endcode
1290  */
1291  template <class KernelIterator>
1292  void initSeparable(KernelIterator kxcenter, int xleft, int xright,
1293  KernelIterator kycenter, int yleft, int yright)
1294  {
1295  vigra_precondition(xleft <= 0 && yleft <= 0,
1296  "Kernel2D::initSeparable(): left borders must be <= 0.");
1297  vigra_precondition(xright >= 0 && yright >= 0,
1298  "Kernel2D::initSeparable(): right borders must be >= 0.");
1299 
1300  left_ = Point2D(xleft, yleft);
1301  right_ = Point2D(xright, yright);
1302 
1303  int w = right_.x - left_.x + 1;
1304  int h = right_.y - left_.y + 1;
1305  kernel_.resize(w, h);
1306 
1307  KernelIterator kiy = kycenter + left_.y;
1308  Iterator iy = center() + left_;
1309 
1310  for(int y=left_.y; y<=right_.y; ++y, ++kiy, ++iy.y)
1311  {
1312  KernelIterator kix = kxcenter + left_.x;
1313  Iterator ix = iy;
1314  for(int x=left_.x; x<=right_.x; ++x, ++kix, ++ix.x)
1315  {
1316  *ix = *kix * *kiy;
1317  }
1318  }
1319 
1320  typename BasicImage<value_type>::iterator i = kernel_.begin();
1321  typename BasicImage<value_type>::iterator iend = kernel_.end();
1322  norm_ = *i;
1323  ++i;
1324 
1325  for(; i!= iend; ++i)
1326  {
1327  norm_ += *i;
1328  }
1329  }
1330 
1331  /** Init as a 2D Gaussian function with given standard deviation and norm.
1332  */
1333  void initGaussian(double std_dev, value_type norm)
1334  {
1335  Kernel1D<value_type> gauss;
1336  gauss.initGaussian(std_dev, norm);
1337  initSeparable(gauss, gauss);
1338  }
1339 
1340  /** Init as a 2D Gaussian function with given standard deviation and unit norm.
1341  */
1342  void initGaussian(double std_dev)
1343  {
1344  initGaussian(std_dev, NumericTraits<value_type>::one());
1345  }
1346 
1347  /** Init the 2D kernel as a circular averaging filter. The norm will be
1348  calculated as
1349  <TT>NumericTraits<value_type>::one() / (number of non-zero kernel values)</TT>.
1350  The kernel's value_type must be a linear space.
1351 
1352  <b> Required Interface:</b>
1353 
1354  \code
1355  value_type v = vigra::NumericTraits<value_type>::one();
1356 
1357  double d;
1358  v = d * v;
1359  \endcode
1360 
1361  <b> Precondition:</b>
1362 
1363  \code
1364  radius > 0;
1365  \endcode
1366  */
1367  void initDisk(int radius)
1368  {
1369  vigra_precondition(radius > 0,
1370  "Kernel2D::initDisk(): radius must be > 0.");
1371 
1372  left_ = Point2D(-radius, -radius);
1373  right_ = Point2D(radius, radius);
1374  int w = right_.x - left_.x + 1;
1375  int h = right_.y - left_.y + 1;
1376  kernel_.resize(w, h);
1377  norm_ = NumericTraits<value_type>::one();
1378 
1379  kernel_ = NumericTraits<value_type>::zero();
1380  double count = 0.0;
1381 
1382  Iterator k = center();
1383  double r2 = (double)radius*radius;
1384 
1385  int i;
1386  for(i=0; i<= radius; ++i)
1387  {
1388  double r = (double) i - 0.5;
1389  int w = (int)(VIGRA_CSTD::sqrt(r2 - r*r) + 0.5);
1390  for(int j=-w; j<=w; ++j)
1391  {
1392  k(j, i) = NumericTraits<value_type>::one();
1393  k(j, -i) = NumericTraits<value_type>::one();
1394  count += (i != 0) ? 2.0 : 1.0;
1395  }
1396  }
1397 
1398  count = 1.0 / count;
1399 
1400  for(int y=-radius; y<=radius; ++y)
1401  {
1402  for(int x=-radius; x<=radius; ++x)
1403  {
1404  k(x,y) = count * k(x,y);
1405  }
1406  }
1407  }
1408 
1409  /** Init the kernel by an explicit initializer list.
1410  The upper left and lower right corners of the kernel must be passed.
1411  A comma-separated initializer list is given after the assignment operator.
1412  This function is used like this:
1413 
1414  \code
1415  // define horizontal Sobel filter
1416  vigra::Kernel2D<float> sobel;
1417 
1418  sobel.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) =
1419  0.125, 0.0, -0.125,
1420  0.25, 0.0, -0.25,
1421  0.125, 0.0, -0.125;
1422  \endcode
1423 
1424  The norm is set to the sum of the initialzer values. If the wrong number of
1425  values is given, a run-time error results. It is, however, possible to give
1426  just one initializer. This creates an averaging filter with the given constant:
1427 
1428  \code
1429  vigra::Kernel2D<float> average3x3;
1430 
1431  average3x3.initExplicitly(Diff2D(-1,-1), Diff2D(1,1)) = 1.0/9.0;
1432  \endcode
1433 
1434  Here, the norm is set to value*width()*height().
1435 
1436  <b> Preconditions:</b>
1437 
1438  \code
1439  1. upperleft.x <= 0;
1440  2. upperleft.y <= 0;
1441  3. lowerright.x >= 0;
1442  4. lowerright.y >= 0;
1443  5. the number of values in the initializer list
1444  is 1 or equals the size of the kernel.
1445  \endcode
1446  */
1447  Kernel2D & initExplicitly(Diff2D upperleft, Diff2D lowerright)
1448  {
1449  vigra_precondition(upperleft.x <= 0 && upperleft.y <= 0,
1450  "Kernel2D::initExplicitly(): left borders must be <= 0.");
1451  vigra_precondition(lowerright.x >= 0 && lowerright.y >= 0,
1452  "Kernel2D::initExplicitly(): right borders must be >= 0.");
1453 
1454  left_ = Point2D(upperleft);
1455  right_ = Point2D(lowerright);
1456 
1457  int w = right_.x - left_.x + 1;
1458  int h = right_.y - left_.y + 1;
1459  kernel_.resize(w, h);
1460 
1461  return *this;
1462  }
1463 
1464  /** Coordinates of the upper left corner of the kernel.
1465  */
1466  Point2D upperLeft() const { return left_; }
1467 
1468  /** Coordinates of the lower right corner of the kernel.
1469  */
1470  Point2D lowerRight() const { return right_; }
1471 
1472  /** Width of the kernel.
1473  */
1474  int width() const { return right_.x - left_.x + 1; }
1475 
1476  /** Height of the kernel.
1477  */
1478  int height() const { return right_.y - left_.y + 1; }
1479 
1480  /** ImageIterator that points to the center of the kernel (coordinate (0,0)).
1481  */
1482  Iterator center() { return kernel_.upperLeft() - left_; }
1483 
1484  /** ImageIterator that points to the center of the kernel (coordinate (0,0)).
1485  */
1486  ConstIterator center() const { return kernel_.upperLeft() - left_; }
1487 
1488  /** Access kernel entry at given position.
1489  */
1490  value_type & operator()(int x, int y)
1491  { return kernel_[Diff2D(x,y) - left_]; }
1492 
1493  /** Read kernel entry at given position.
1494  */
1495  value_type operator()(int x, int y) const
1496  { return kernel_[Diff2D(x,y) - left_]; }
1497 
1498  /** Access kernel entry at given position.
1499  */
1501  { return kernel_[d - left_]; }
1502 
1503  /** Read kernel entry at given position.
1504  */
1505  value_type operator[](Diff2D const & d) const
1506  { return kernel_[d - left_]; }
1507 
1508  /** Norm of the kernel (i.e. sum of its elements).
1509  */
1510  value_type norm() const { return norm_; }
1511 
1512  /** The kernels default accessor.
1513  */
1514  Accessor accessor() { return Accessor(); }
1515 
1516  /** The kernels default const accessor.
1517  */
1518  ConstAccessor accessor() const { return ConstAccessor(); }
1519 
1520  /** Normalize the kernel to the given value. (The norm is the sum of all kernel
1521  elements.) The kernel's value_type must be a division algebra or
1522  algebraic field.
1523 
1524  <b> Required Interface:</b>
1525 
1526  \code
1527  value_type v = vigra::NumericTraits<value_type>::one(); // if norm is not
1528  // given explicitly
1529 
1530  v += v;
1531  v = v * v;
1532  v = v / v;
1533  \endcode
1534  */
1536  {
1537  typename BasicImage<value_type>::iterator i = kernel_.begin();
1538  typename BasicImage<value_type>::iterator iend = kernel_.end();
1539  typename NumericTraits<value_type>::RealPromote sum = *i;
1540  ++i;
1541 
1542  for(; i!= iend; ++i)
1543  {
1544  sum += *i;
1545  }
1546 
1547  sum = norm / sum;
1548  i = kernel_.begin();
1549  for(; i != iend; ++i)
1550  {
1551  *i = *i * sum;
1552  }
1553 
1554  norm_ = norm;
1555  }
1556 
1557  /** Normalize the kernel to norm 1.
1558  */
1559  void normalize()
1560  {
1561  normalize(one());
1562  }
1563 
1564  /** current border treatment mode
1565  */
1566  BorderTreatmentMode borderTreatment() const
1567  { return border_treatment_; }
1568 
1569  /** Set border treatment mode.
1570  Only <TT>BORDER_TREATMENT_CLIP</TT> and <TT>BORDER_TREATMENT_AVOID</TT> are currently
1571  allowed.
1572  */
1573  void setBorderTreatment( BorderTreatmentMode new_mode)
1574  {
1575  vigra_precondition((new_mode == BORDER_TREATMENT_CLIP ||
1576  new_mode == BORDER_TREATMENT_AVOID ||
1577  new_mode == BORDER_TREATMENT_REFLECT ||
1578  new_mode == BORDER_TREATMENT_REPEAT ||
1579  new_mode == BORDER_TREATMENT_WRAP),
1580  "convolveImage():\n"
1581  " Border treatment must be one of follow treatments:\n"
1582  " - BORDER_TREATMENT_CLIP\n"
1583  " - BORDER_TREATMENT_AVOID\n"
1584  " - BORDER_TREATMENT_REFLECT\n"
1585  " - BORDER_TREATMENT_REPEAT\n"
1586  " - BORDER_TREATMENT_WRAP\n");
1587 
1588  border_treatment_ = new_mode;
1589  }
1590 
1591 
1592 private:
1593  BasicImage<value_type> kernel_;
1594  Point2D left_, right_;
1595  value_type norm_;
1596  BorderTreatmentMode border_treatment_;
1597 };
1598 
1599 /**************************************************************/
1600 /* */
1601 /* Argument object factories for Kernel2D */
1602 /* */
1603 /* (documentation: see vigra/convolution.hxx) */
1604 /* */
1605 /**************************************************************/
1606 
1607 template <class KernelIterator, class KernelAccessor>
1608 inline
1609 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, BorderTreatmentMode>
1610 kernel2d(KernelIterator ik, KernelAccessor ak, Diff2D kul, Diff2D klr,
1611  BorderTreatmentMode border)
1612 
1613 {
1614  return
1615  tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, BorderTreatmentMode> (
1616  ik, ak, kul, klr, border);
1617 }
1618 
1619 template <class T>
1620 inline
1621 tuple5<typename Kernel2D<T>::ConstIterator,
1622  typename Kernel2D<T>::ConstAccessor,
1623  Diff2D, Diff2D, BorderTreatmentMode>
1624 kernel2d(Kernel2D<T> const & k)
1625 
1626 {
1627  return
1628  tuple5<typename Kernel2D<T>::ConstIterator,
1629  typename Kernel2D<T>::ConstAccessor,
1630  Diff2D, Diff2D, BorderTreatmentMode>(
1631  k.center(),
1632  k.accessor(),
1633  k.upperLeft(), k.lowerRight(),
1634  k.borderTreatment());
1635 }
1636 
1637 template <class T>
1638 inline
1639 tuple5<typename Kernel2D<T>::ConstIterator,
1640  typename Kernel2D<T>::ConstAccessor,
1641  Diff2D, Diff2D, BorderTreatmentMode>
1642 kernel2d(Kernel2D<T> const & k, BorderTreatmentMode border)
1643 
1644 {
1645  return
1646  tuple5<typename Kernel2D<T>::ConstIterator,
1647  typename Kernel2D<T>::ConstAccessor,
1648  Diff2D, Diff2D, BorderTreatmentMode>(
1649  k.center(),
1650  k.accessor(),
1651  k.upperLeft(), k.lowerRight(),
1652  border);
1653 }
1654 
1655 
1656 } // namespace vigra
1657 
1658 #endif // VIGRA_STDCONVOLUTION_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)