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

resampling_convolution.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2004 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 #ifndef VIGRA_RESAMPLING_CONVOLUTION_HXX
37 #define VIGRA_RESAMPLING_CONVOLUTION_HXX
38 
39 #include <cmath>
40 #include "stdimage.hxx"
41 #include "array_vector.hxx"
42 #include "rational.hxx"
43 #include "functortraits.hxx"
44 #include "functorexpression.hxx"
45 #include "transformimage.hxx"
46 #include "imagecontainer.hxx"
47 #include "multi_shape.hxx"
48 
49 namespace vigra {
50 
51 namespace resampling_detail
52 {
53 
54 struct MapTargetToSourceCoordinate
55 {
56  MapTargetToSourceCoordinate(Rational<int> const & samplingRatio,
57  Rational<int> const & offset)
58  : a(samplingRatio.denominator()*offset.denominator()),
59  b(samplingRatio.numerator()*offset.numerator()),
60  c(samplingRatio.numerator()*offset.denominator())
61  {}
62 
63 // the following functions are more efficient realizations of:
64 // rational_cast<T>(i / samplingRatio + offset);
65 // we need efficiency because this may be called in the inner loop
66 
67  int operator()(int i) const
68  {
69  return (i * a + b) / c;
70  }
71 
72  double toDouble(int i) const
73  {
74  return double(i * a + b) / c;
75  }
76 
77  Rational<int> toRational(int i) const
78  {
79  return Rational<int>(i * a + b, c);
80  }
81 
82  bool isExpand2() const
83  {
84  return a == 1 && b == 0 && c == 2;
85  }
86 
87  bool isReduce2() const
88  {
89  return a == 2 && b == 0 && c == 1;
90  }
91 
92  int a, b, c;
93 };
94 
95 } // namespace resampling_detail
96 
97 template <>
98 class FunctorTraits<resampling_detail::MapTargetToSourceCoordinate>
99 : public FunctorTraitsBase<resampling_detail::MapTargetToSourceCoordinate>
100 {
101  public:
102  typedef VigraTrueType isUnaryFunctor;
103 };
104 
105 template <class SrcIter, class SrcAcc,
106  class DestIter, class DestAcc,
107  class KernelArray>
108 void
109 resamplingExpandLine2(SrcIter s, SrcIter send, SrcAcc src,
110  DestIter d, DestIter dend, DestAcc dest,
111  KernelArray const & kernels)
112 {
113  typedef typename KernelArray::value_type Kernel;
114  typedef typename KernelArray::const_reference KernelRef;
115  typedef typename Kernel::const_iterator KernelIter;
116 
117  typedef typename
118  PromoteTraits<typename SrcAcc::value_type, typename Kernel::value_type>::Promote
119  TmpType;
120 
121  int wo = send - s;
122  int wn = dend - d;
123  int wo2 = 2*wo - 2;
124 
125  int ileft = std::max(kernels[0].right(), kernels[1].right());
126  int iright = wo + std::min(kernels[0].left(), kernels[1].left()) - 1;
127  for(int i = 0; i < wn; ++i, ++d)
128  {
129  int is = i / 2;
130  KernelRef kernel = kernels[i & 1];
131  KernelIter k = kernel.center() + kernel.right();
132  TmpType sum = NumericTraits<TmpType>::zero();
133  if(is < ileft)
134  {
135  for(int m=is-kernel.right(); m <= is-kernel.left(); ++m, --k)
136  {
137  int mm = (m < 0)
138  ? -m
139  : m;
140  sum += *k * src(s, mm);
141  }
142  }
143  else if(is > iright)
144  {
145  for(int m=is-kernel.right(); m <= is-kernel.left(); ++m, --k)
146  {
147  int mm = (m >= wo)
148  ? wo2 - m
149  : m;
150  sum += *k * src(s, mm);
151  }
152  }
153  else
154  {
155  SrcIter ss = s + is - kernel.right();
156  for(int m = 0; m < kernel.size(); ++m, --k, ++ss)
157  {
158  sum += *k * src(ss);
159  }
160  }
161  dest.set(sum, d);
162  }
163 }
164 
165 template <class SrcIter, class SrcAcc,
166  class DestIter, class DestAcc,
167  class KernelArray>
168 void
169 resamplingReduceLine2(SrcIter s, SrcIter send, SrcAcc src,
170  DestIter d, DestIter dend, DestAcc dest,
171  KernelArray const & kernels)
172 {
173  typedef typename KernelArray::value_type Kernel;
174  typedef typename KernelArray::const_reference KernelRef;
175  typedef typename Kernel::const_iterator KernelIter;
176 
177  KernelRef kernel = kernels[0];
178  KernelIter kbegin = kernel.center() + kernel.right();
179 
180  typedef typename
181  PromoteTraits<typename SrcAcc::value_type, typename Kernel::value_type>::Promote
182  TmpType;
183 
184  int wo = send - s;
185  int wn = dend - d;
186  int wo2 = 2*wo - 2;
187 
188  int ileft = kernel.right();
189  int iright = wo + kernel.left() - 1;
190  for(int i = 0; i < wn; ++i, ++d)
191  {
192  int is = 2 * i;
193  KernelIter k = kbegin;
194  TmpType sum = NumericTraits<TmpType>::zero();
195  if(is < ileft)
196  {
197  for(int m=is-kernel.right(); m <= is-kernel.left(); ++m, --k)
198  {
199  int mm = (m < 0)
200  ? -m
201  : m;
202  sum += *k * src(s, mm);
203  }
204  }
205  else if(is > iright)
206  {
207  for(int m=is-kernel.right(); m <= is-kernel.left(); ++m, --k)
208  {
209  int mm = (m >= wo)
210  ? wo2 - m
211  : m;
212  sum += *k * src(s, mm);
213  }
214  }
215  else
216  {
217  SrcIter ss = s + is - kernel.right();
218  for(int m = 0; m < kernel.size(); ++m, --k, ++ss)
219  {
220  sum += *k * src(ss);
221  }
222  }
223  dest.set(sum, d);
224  }
225 }
226 
227 /** \addtogroup ResamplingConvolutionFilters Resampling Convolution Filters
228 
229  These functions implement the convolution operation when the source and target images
230  have different sizes. This is realized by accessing a continuous kernel at the
231  appropriate non-integer positions. The technique is, for example, described in
232  D. Schumacher: <i>General Filtered Image Rescaling</i>, in: Graphics Gems III,
233  Academic Press, 1992.
234 */
235 //@{
236 
237 /********************************************************/
238 /* */
239 /* resamplingConvolveLine */
240 /* */
241 /********************************************************/
242 
243 /** \brief Performs a 1-dimensional resampling convolution of the source signal using the given
244  set of kernels.
245 
246  This function is mainly used internally: It is called for each dimension of a
247  higher dimensional array in order to perform a separable resize operation.
248 
249  <b> Declaration:</b>
250 
251  <b>\#include</b> <vigra/resampling_convolution.hxx><br/>
252  Namespace: vigra
253 
254  \code
255  namespace vigra {
256  template <class SrcIter, class SrcAcc,
257  class DestIter, class DestAcc,
258  class KernelArray,
259  class Functor>
260  void
261  resamplingConvolveLine(SrcIter s, SrcIter send, SrcAcc src,
262  DestIter d, DestIter dend, DestAcc dest,
263  KernelArray const & kernels,
264  Functor mapTargetToSourceCoordinate)
265  }
266  \endcode
267 
268 */
269 doxygen_overloaded_function(template <...> void resamplingConvolveLine)
270 
271 template <class SrcIter, class SrcAcc,
272  class DestIter, class DestAcc,
273  class KernelArray,
274  class Functor>
275 void
276 resamplingConvolveLine(SrcIter s, SrcIter send, SrcAcc src,
277  DestIter d, DestIter dend, DestAcc dest,
278  KernelArray const & kernels,
279  Functor mapTargetToSourceCoordinate)
280 {
281  if(mapTargetToSourceCoordinate.isExpand2())
282  {
283  resamplingExpandLine2(s, send, src, d, dend, dest, kernels);
284  return;
285  }
286  if(mapTargetToSourceCoordinate.isReduce2())
287  {
288  resamplingReduceLine2(s, send, src, d, dend, dest, kernels);
289  return;
290  }
291 
292  typedef typename
293  NumericTraits<typename SrcAcc::value_type>::RealPromote
294  TmpType;
295  typedef typename KernelArray::value_type Kernel;
296  typedef typename Kernel::const_iterator KernelIter;
297 
298  int wo = send - s;
299  int wn = dend - d;
300  int wo2 = 2*wo - 2;
301 
302  int i;
303  typename KernelArray::const_iterator kernel = kernels.begin();
304  for(i=0; i<wn; ++i, ++d, ++kernel)
305  {
306  // use the kernels periodically
307  if(kernel == kernels.end())
308  kernel = kernels.begin();
309 
310  // calculate current target point into source location
311  int is = mapTargetToSourceCoordinate(i);
312 
313  TmpType sum = NumericTraits<TmpType>::zero();
314 
315  int lbound = is - kernel->right(),
316  hbound = is - kernel->left();
317 
318  KernelIter k = kernel->center() + kernel->right();
319  if(lbound < 0 || hbound >= wo)
320  {
321  vigra_precondition(-lbound < wo && wo2 - hbound >= 0,
322  "resamplingConvolveLine(): kernel or offset larger than image.");
323  for(int m=lbound; m <= hbound; ++m, --k)
324  {
325  int mm = (m < 0) ?
326  -m :
327  (m >= wo) ?
328  wo2 - m :
329  m;
330  sum = TmpType(sum + *k * src(s, mm));
331  }
332  }
333  else
334  {
335  SrcIter ss = s + lbound;
336  SrcIter ssend = s + hbound;
337 
338  for(; ss <= ssend; ++ss, --k)
339  {
340  sum = TmpType(sum + *k * src(ss));
341  }
342  }
343 
344  dest.set(sum, d);
345  }
346 }
347 
348 template <class Kernel, class MapCoordinate, class KernelArray>
349 void
350 createResamplingKernels(Kernel const & kernel,
351  MapCoordinate const & mapCoordinate, KernelArray & kernels)
352 {
353  for(unsigned int idest = 0; idest < kernels.size(); ++idest)
354  {
355  int isrc = mapCoordinate(idest);
356  double idsrc = mapCoordinate.toDouble(idest);
357  double offset = idsrc - isrc;
358  double radius = kernel.radius();
359  int left = std::min(0, int(ceil(-radius - offset)));
360  int right = std::max(0, int(floor(radius - offset)));
361  kernels[idest].initExplicitly(left, right);
362 
363  double x = left + offset;
364  for(int i = left; i <= right; ++i, ++x)
365  kernels[idest][i] = kernel(x);
366  kernels[idest].normalize(1.0, kernel.derivativeOrder(), offset);
367  }
368 }
369 
370 /** \brief Apply a resampling filter in the x-direction.
371 
372  This function implements a convolution operation in x-direction
373  (i.e. applies a 1D filter to every row) where the width of the source
374  and destination images differ. This is typically used to avoid aliasing if
375  the image is scaled down, or to interpolate smoothly if the image is scaled up.
376  The target coordinates are transformed into source coordinates by
377 
378  \code
379  xsource = (xtarget - offset) / samplingRatio
380  \endcode
381 
382  The <tt>samplingRatio</tt> and <tt>offset</tt> must be given as \ref vigra::Rational
383  in order to avoid rounding errors in this transformation. It is required that for all
384  pixels of the target image, <tt>xsource</tt> remains within the range of the source
385  image (i.e. <tt>0 <= xsource <= sourceWidth-1</tt>. Since <tt>xsource</tt> is
386  in general not an integer, the <tt>kernel</tt> must be a functor that can be accessed at
387  arbitrary (<tt>double</tt>) coordinates. It must also provide a member function <tt>radius()</tt>
388  which specifies the support (non-zero interval) of the kernel. VIGRA already
389  provides a number of suitable functors, e.g. \ref vigra::Gaussian, \ref vigra::BSpline
390  \ref vigra::CatmullRomSpline, and \ref vigra::CoscotFunction. The function
391  \ref resizeImageSplineInterpolation() is implemented by means of resamplingConvolveX() and
392  resamplingConvolveY().
393 
394  <b> Declarations:</b>
395 
396  pass 2D array views:
397  \code
398  namespace vigra {
399  template <class T1, class S1,
400  class T2, class S2,
401  class Kernel>
402  void
403  resamplingConvolveX(MultiArrayView<2, T1, S1> const & src,
404  MultiArrayView<2, T2, S2> dest,
405  Kernel const & kernel,
406  Rational<int> const & samplingRatio, Rational<int> const & offset);
407  }
408  \endcode
409 
410  \deprecatedAPI{resamplingConvolveX}
411  pass \ref ImageIterators and \ref DataAccessors :
412  \code
413  namespace vigra {
414  template <class SrcIter, class SrcAcc,
415  class DestIter, class DestAcc,
416  class Kernel>
417  void
418  resamplingConvolveX(SrcIter sul, SrcIter slr, SrcAcc src,
419  DestIter dul, DestIter dlr, DestAcc dest,
420  Kernel const & kernel,
421  Rational<int> const & samplingRatio, Rational<int> const & offset);
422  }
423  \endcode
424  use argument objects in conjunction with \ref ArgumentObjectFactories :
425  \code
426  namespace vigra {
427  template <class SrcIter, class SrcAcc,
428  class DestIter, class DestAcc,
429  class Kernel>
430  void
431  resamplingConvolveX(triple<SrcIter, SrcIter, SrcAcc> src,
432  triple<DestIter, DestIter, DestAcc> dest,
433  Kernel const & kernel,
434  Rational<int> const & samplingRatio, Rational<int> const & offset);
435  }
436  \endcode
437  \deprecatedEnd
438 
439  <b> Usage:</b>
440 
441  <b>\#include</b> <vigra/resampling_convolution.hxx><br/>
442  Namespace: vigra
443 
444  \code
445  Rational<int> ratio(2), offset(0);
446 
447  MultiArray<2, float> src(w,h),
448  dest(rational_cast<int>(ratio*w), h);
449 
450  float sigma = 2.0;
451  Gaussian<float> smooth(sigma);
452  ...
453 
454  // simultaneously enlarge and smooth source image
455  resamplingConvolveX(src, dest,
456  smooth, ratio, offset);
457  \endcode
458 
459  \deprecatedUsage{resamplingConvolveX}
460  \code
461  Rational<int> ratio(2), offset(0);
462 
463  FImage src(w,h),
464  dest(rational_cast<int>(ratio*w), h);
465 
466  float sigma = 2.0;
467  Gaussian<float> smooth(sigma);
468  ...
469 
470  // simultaneously enlarge and smooth source image
471  resamplingConvolveX(srcImageRange(src), destImageRange(dest),
472  smooth, ratio, offset);
473  \endcode
474  <b> Required Interface:</b>
475  \code
476  Kernel kernel;
477  int kernelRadius = kernel.radius();
478  double x = ...; // must be <= radius()
479  double value = kernel(x);
480  \endcode
481  \deprecatedEnd
482 */
483 doxygen_overloaded_function(template <...> void resamplingConvolveX)
484 
485 template <class SrcIter, class SrcAcc,
486  class DestIter, class DestAcc,
487  class Kernel>
488 void
489 resamplingConvolveX(SrcIter sul, SrcIter slr, SrcAcc src,
490  DestIter dul, DestIter dlr, DestAcc dest,
491  Kernel const & kernel,
492  Rational<int> const & samplingRatio, Rational<int> const & offset)
493 {
494  int wold = slr.x - sul.x;
495  int wnew = dlr.x - dul.x;
496 
497  vigra_precondition(!samplingRatio.is_inf() && samplingRatio > 0,
498  "resamplingConvolveX(): sampling ratio must be > 0 and < infinity");
499  vigra_precondition(!offset.is_inf(),
500  "resamplingConvolveX(): offset must be < infinity");
501 
502  int period = lcm(samplingRatio.numerator(), samplingRatio.denominator());
503  resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);
504 
505  ArrayVector<Kernel1D<double> > kernels(period);
506 
507  createResamplingKernels(kernel, mapCoordinate, kernels);
508 
509  for(; sul.y < slr.y; ++sul.y, ++dul.y)
510  {
511  typename SrcIter::row_iterator sr = sul.rowIterator();
512  typename DestIter::row_iterator dr = dul.rowIterator();
513  resamplingConvolveLine(sr, sr+wold, src, dr, dr+wnew, dest,
514  kernels, mapCoordinate);
515  }
516 }
517 
518 template <class SrcIter, class SrcAcc,
519  class DestIter, class DestAcc,
520  class Kernel>
521 inline void
522 resamplingConvolveX(triple<SrcIter, SrcIter, SrcAcc> src,
523  triple<DestIter, DestIter, DestAcc> dest,
524  Kernel const & kernel,
525  Rational<int> const & samplingRatio, Rational<int> const & offset)
526 {
527  resamplingConvolveX(src.first, src.second, src.third,
528  dest.first, dest.second, dest.third,
529  kernel, samplingRatio, offset);
530 }
531 
532 template <class T1, class S1,
533  class T2, class S2,
534  class Kernel>
535 inline void
536 resamplingConvolveX(MultiArrayView<2, T1, S1> const & src,
537  MultiArrayView<2, T2, S2> dest,
538  Kernel const & kernel,
539  Rational<int> const & samplingRatio, Rational<int> const & offset)
540 {
541  resamplingConvolveX(srcImageRange(src),
542  destImageRange(dest),
543  kernel, samplingRatio, offset);
544 }
545 
546 /********************************************************/
547 /* */
548 /* resamplingConvolveY */
549 /* */
550 /********************************************************/
551 
552 /** \brief Apply a resampling filter in the y-direction.
553 
554  This function implements a convolution operation in y-direction
555  (i.e. applies a 1D filter to every column) where the height of the source
556  and destination images differ. This is typically used to avoid aliasing if
557  the image is scaled down, or to interpolate smoothly if the image is scaled up.
558  The target coordinates are transformed into source coordinates by
559 
560  \code
561  ysource = (ytarget - offset) / samplingRatio
562  \endcode
563 
564  The <tt>samplingRatio</tt> and <tt>offset</tt> must be given as \ref vigra::Rational
565  in order to avoid rounding errors in this transformation. It is required that for all
566  pixels of the target image, <tt>ysource</tt> remains within the range of the source
567  image (i.e. <tt>0 <= ysource <= sourceHeight-1</tt>. Since <tt>ysource</tt> is
568  in general not an integer, the <tt>kernel</tt> must be a functor that can be accessed at
569  arbitrary (<tt>double</tt>) coordinates. It must also provide a member function <tt>radius()</tt>
570  which specifies the support (non-zero interval) of the kernel. VIGRA already
571  provides a number of suitable functors, e.g. \ref vigra::Gaussian, \ref vigra::BSpline
572  \ref vigra::CatmullRomSpline, and \ref vigra::CoscotFunction. The function
573  \ref resizeImageSplineInterpolation() is implemented by means of resamplingConvolveX() and
574  resamplingConvolveY().
575 
576  <b> Declarations:</b>
577 
578  pass 2D array views:
579  \code
580  namespace vigra {
581  template <class T1, class S1,
582  class T2, class S2,
583  class Kernel>
584  void
585  resamplingConvolveY(MultiArrayView<2, T1, S1> const & src,
586  MultiArrayView<2, T2, S2> dest,
587  Kernel const & kernel,
588  Rational<int> const & samplingRatio, Rational<int> const & offset);
589  }
590  \endcode
591 
592  \deprecatedAPI{resamplingConvolveY}
593  pass \ref ImageIterators and \ref DataAccessors :
594  \code
595  namespace vigra {
596  template <class SrcIter, class SrcAcc,
597  class DestIter, class DestAcc,
598  class Kernel>
599  void
600  resamplingConvolveY(SrcIter sul, SrcIter slr, SrcAcc src,
601  DestIter dul, DestIter dlr, DestAcc dest,
602  Kernel const & kernel,
603  Rational<int> const & samplingRatio, Rational<int> const & offset);
604  }
605  \endcode
606  use argument objects in conjunction with \ref ArgumentObjectFactories :
607  \code
608  namespace vigra {
609  template <class SrcIter, class SrcAcc,
610  class DestIter, class DestAcc,
611  class Kernel>
612  void
613  resamplingConvolveY(triple<SrcIter, SrcIter, SrcAcc> src,
614  triple<DestIter, DestIter, DestAcc> dest,
615  Kernel const & kernel,
616  Rational<int> const & samplingRatio, Rational<int> const & offset);
617  }
618  \endcode
619  \deprecatedEnd
620 
621  <b> Usage:</b>
622 
623  <b>\#include</b> <vigra/resampling_convolution.hxx><br/>
624  Namespace: vigra
625 
626  \code
627  Rational<int> ratio(2), offset(0);
628 
629  MultiArray<2, float> src(w,h),
630  dest(w, rational_cast<int>(ratio*h));
631 
632  float sigma = 2.0;
633  Gaussian<float> smooth(sigma);
634  ...
635 
636  // simultaneously enlarge and smooth source image
637  resamplingConvolveY(src, dest,
638  smooth, ratio, offset);
639  \endcode
640 
641  \deprecatedUsage{resamplingConvolveY}
642  \code
643  Rational<int> ratio(2), offset(0);
644 
645  FImage src(w,h),
646  dest(w, rational_cast<int>(ratio*h));
647 
648  float sigma = 2.0;
649  Gaussian<float> smooth(sigma);
650  ...
651 
652  // simultaneously enlarge and smooth source image
653  resamplingConvolveY(srcImageRange(src), destImageRange(dest),
654  smooth, ratio, offset);
655  \endcode
656  <b> Required Interface:</b>
657  \code
658  Kernel kernel;
659  int kernelRadius = kernel.radius();
660  double y = ...; // must be <= radius()
661  double value = kernel(y);
662  \endcode
663  \deprecatedEnd
664 */
665 doxygen_overloaded_function(template <...> void resamplingConvolveY)
666 
667 template <class SrcIter, class SrcAcc,
668  class DestIter, class DestAcc,
669  class Kernel>
670 void
671 resamplingConvolveY(SrcIter sul, SrcIter slr, SrcAcc src,
672  DestIter dul, DestIter dlr, DestAcc dest,
673  Kernel const & kernel,
674  Rational<int> const & samplingRatio, Rational<int> const & offset)
675 {
676  int hold = slr.y - sul.y;
677  int hnew = dlr.y - dul.y;
678 
679  vigra_precondition(!samplingRatio.is_inf() && samplingRatio > 0,
680  "resamplingConvolveY(): sampling ratio must be > 0 and < infinity");
681  vigra_precondition(!offset.is_inf(),
682  "resamplingConvolveY(): offset must be < infinity");
683 
684  int period = lcm(samplingRatio.numerator(), samplingRatio.denominator());
685 
686  resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);
687 
688  ArrayVector<Kernel1D<double> > kernels(period);
689 
690  createResamplingKernels(kernel, mapCoordinate, kernels);
691 
692  for(; sul.x < slr.x; ++sul.x, ++dul.x)
693  {
694  typename SrcIter::column_iterator sc = sul.columnIterator();
695  typename DestIter::column_iterator dc = dul.columnIterator();
696  resamplingConvolveLine(sc, sc+hold, src, dc, dc+hnew, dest,
697  kernels, mapCoordinate);
698  }
699 }
700 
701 template <class SrcIter, class SrcAccessor,
702  class DestIter, class DestAccessor,
703  class Kernel>
704 inline void
705 resamplingConvolveY(triple<SrcIter, SrcIter, SrcAccessor> src,
706  triple<DestIter, DestIter, DestAccessor> dest,
707  Kernel const & kernel,
708  Rational<int> const & samplingRatio, Rational<int> const & offset)
709 {
710  resamplingConvolveY(src.first, src.second, src.third,
711  dest.first, dest.second, dest.third,
712  kernel, samplingRatio, offset);
713 }
714 
715 template <class T1, class S1,
716  class T2, class S2,
717  class Kernel>
718 inline void
719 resamplingConvolveY(MultiArrayView<2, T1, S1> const & src,
720  MultiArrayView<2, T2, S2> dest,
721  Kernel const & kernel,
722  Rational<int> const & samplingRatio, Rational<int> const & offset)
723 {
724  resamplingConvolveY(srcImageRange(src),
725  destImageRange(dest),
726  kernel, samplingRatio, offset);
727 }
728 
729 /********************************************************/
730 /* */
731 /* resamplingConvolveImage */
732 /* */
733 /********************************************************/
734 
735 /** \brief Apply two separable resampling filters successively, the first in x-direction,
736  the second in y-direction.
737 
738  This function is a shorthand for the concatenation of a call to
739  \ref resamplingConvolveX() and \ref resamplingConvolveY()
740  with the given kernels. See there for detailed documentation.
741 
742  <b> Declarations:</b>
743 
744  pass 2D array views:
745  \code
746  namespace vigra {
747  template <class T1, class S1,
748  class T2, class S2,
749  class KernelX, class KernelY>
750  void
751  resamplingConvolveImage(MultiArrayView<2, T1, S1> const & src,
752  MultiArrayView<2, T2, S2> dest,
753  KernelX const & kx,
754  Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
755  KernelY const & ky,
756  Rational<int> const & samplingRatioY, Rational<int> const & offsetY);
757  }
758  \endcode
759 
760  \deprecatedAPI{resamplingConvolveImage}
761  pass \ref ImageIterators and \ref DataAccessors :
762  \code
763  namespace vigra {
764  template <class SrcIterator, class SrcAccessor,
765  class DestIterator, class DestAccessor,
766  class KernelX, class KernelY>
767  void resamplingConvolveImage(SrcIterator sul,SrcIterator slr, SrcAccessor src,
768  DestIterator dul, DestIterator dlr, DestAccessor dest,
769  KernelX const & kx,
770  Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
771  KernelY const & ky,
772  Rational<int> const & samplingRatioY, Rational<int> const & offsetY);
773  }
774  \endcode
775  use argument objects in conjunction with \ref ArgumentObjectFactories :
776  \code
777  namespace vigra {
778  template <class SrcIterator, class SrcAccessor,
779  class DestIterator, class DestAccessor,
780  class KernelX, class KernelY>
781  void
782  resamplingConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
783  triple<DestIterator, DestIterator, DestAccessor> dest,
784  KernelX const & kx,
785  Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
786  KernelY const & ky,
787  Rational<int> const & samplingRatioY, Rational<int> const & offsetY);
788  }
789  \endcode
790  \deprecatedEnd
791 
792  <b> Usage:</b>
793 
794  <b>\#include</b> <vigra/resampling_convolution.hxx><br/>
795  Namespace: vigra
796 
797  \code
798  Rational<int> xratio(2), yratio(3), offset(0);
799 
800  MultiArray<2, float> src(w,h),
801  dest(rational_cast<int>(xratio*w), rational_cast<int>(yratio*h));
802 
803  float sigma = 2.0;
804  Gaussian<float> smooth(sigma);
805  ...
806 
807  // simultaneously enlarge and smooth source image
808  resamplingConvolveImage(src, dest,
809  smooth, xratio, offset,
810  smooth, yratio, offset);
811  \endcode
812 */
813 doxygen_overloaded_function(template <...> void resamplingConvolveImage)
814 
815 template <class SrcIterator, class SrcAccessor,
816  class DestIterator, class DestAccessor,
817  class KernelX, class KernelY>
818 void resamplingConvolveImage(SrcIterator sul,SrcIterator slr, SrcAccessor src,
819  DestIterator dul, DestIterator dlr, DestAccessor dest,
820  KernelX const & kx,
821  Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
822  KernelY const & ky,
823  Rational<int> const & samplingRatioY, Rational<int> const & offsetY)
824 {
825  typedef typename
826  NumericTraits<typename SrcAccessor::value_type>::RealPromote
827  TmpType;
828 
829  BasicImage<TmpType> tmp(dlr.x - dul.x, slr.y - sul.y);
830 
831  resamplingConvolveX(srcIterRange(sul, slr, src),
832  destImageRange(tmp),
833  kx, samplingRatioX, offsetX);
834  resamplingConvolveY(srcImageRange(tmp),
835  destIterRange(dul, dlr, dest),
836  ky, samplingRatioY, offsetY);
837 }
838 
839 template <class SrcIterator, class SrcAccessor,
840  class DestIterator, class DestAccessor,
841  class KernelX, class KernelY>
842 inline void
843 resamplingConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
844  triple<DestIterator, DestIterator, DestAccessor> dest,
845  KernelX const & kx,
846  Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
847  KernelY const & ky,
848  Rational<int> const & samplingRatioY, Rational<int> const & offsetY)
849 {
850  resamplingConvolveImage(src.first, src.second, src.third,
851  dest.first, dest.second, dest.third,
852  kx, samplingRatioX, offsetX,
853  ky, samplingRatioY, offsetY);
854 }
855 
856 template <class T1, class S1,
857  class T2, class S2,
858  class KernelX, class KernelY>
859 inline void
860 resamplingConvolveImage(MultiArrayView<2, T1, S1> const & src,
861  MultiArrayView<2, T2, S2> dest,
862  KernelX const & kx,
863  Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
864  KernelY const & ky,
865  Rational<int> const & samplingRatioY, Rational<int> const & offsetY)
866 {
867  resamplingConvolveImage(srcImageRange(src),
868  destImageRange(dest),
869  kx, samplingRatioX, offsetX,
870  ky, samplingRatioY, offsetY);
871 }
872 
873 /** \brief Two-fold down-sampling for image pyramid construction.
874 
875  Sorry, no \ref detailedDocumentation() available yet.
876 
877  <b> Declarations:</b>
878 
879  <b>\#include</b> <vigra/resampling_convolution.hxx><br>
880  Namespace: vigra
881 
882  pass 2D array views:
883  \code
884  namespace vigra {
885  template <class SrcIterator, class SrcAccessor,
886  class DestIterator, class DestAccessor>
887  void pyramidReduceBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
888  DestIterator dul, DestIterator dlr, DestAccessor dest,
889  double centerValue = 0.4);
890  }
891  \endcode
892 
893  \deprecatedAPI{pyramidReduceBurtFilter}
894  pass \ref ImageIterators and \ref DataAccessors :
895  \code
896  namespace vigra {
897  template <class SrcIterator, class SrcAccessor,
898  class DestIterator, class DestAccessor>
899  void pyramidReduceBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
900  DestIterator dul, DestIterator dlr, DestAccessor dest,
901  double centerValue = 0.4);
902  }
903  \endcode
904  use argument objects in conjunction with \ref ArgumentObjectFactories :
905  \code
906  namespace vigra {
907  template <class SrcIterator, class SrcAccessor,
908  class DestIterator, class DestAccessor>
909  void pyramidReduceBurtFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
910  triple<DestIterator, DestIterator, DestAccessor> dest,
911  double centerValue = 0.4);
912  }
913  \endcode
914  \deprecatedEnd
915 
916  use a \ref vigra::ImagePyramid :
917  \code
918  namespace vigra {
919  template <class Image, class Alloc>
920  void pyramidReduceBurtFilter(ImagePyramid<Image, Alloc> & pyramid, int fromLevel, int toLevel,
921  double centerValue = 0.4);
922  }
923  \endcode
924 */
925 doxygen_overloaded_function(template <...> void pyramidReduceBurtFilter)
926 
927 template <class SrcIterator, class SrcAccessor,
928  class DestIterator, class DestAccessor>
929 void pyramidReduceBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
930  DestIterator dul, DestIterator dlr, DestAccessor dest,
931  double centerValue = 0.4)
932 {
933  vigra_precondition(0.25 <= centerValue && centerValue <= 0.5,
934  "pyramidReduceBurtFilter(): centerValue must be between 0.25 and 0.5.");
935 
936  int wold = slr.x - sul.x;
937  int wnew = dlr.x - dul.x;
938  int hold = slr.y - sul.y;
939  int hnew = dlr.y - dul.y;
940 
941  vigra_precondition(wnew == (wold + 1) / 2 && hnew == (hold + 1) / 2,
942  "pyramidReduceBurtFilter(): oldSize = ceil(newSize / 2) required.");
943 
944  vigra_precondition(wnew == (wold + 1) / 2 && hnew == (hold + 1) / 2,
945  "pyramidReduceBurtFilter(): oldSize = ceil(newSize / 2) required.");
946 
947  Rational<int> samplingRatio(1,2), offset(0);
948  resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);
949 
950  ArrayVector<Kernel1D<double> > kernels(1);
951  kernels[0].initExplicitly(-2, 2) = 0.25 - centerValue / 2.0, 0.25, centerValue, 0.25, 0.25 - centerValue / 2.0;
952 
953  typedef typename
954  NumericTraits<typename SrcAccessor::value_type>::RealPromote
955  TmpType;
956  typedef BasicImage<TmpType> TmpImage;
957  typedef typename TmpImage::traverser TmpIterator;
958 
959  BasicImage<TmpType> tmp(wnew, hold);
960 
961  TmpIterator tul = tmp.upperLeft();
962 
963  for(; sul.y < slr.y; ++sul.y, ++tul.y)
964  {
965  typename SrcIterator::row_iterator sr = sul.rowIterator();
966  typename TmpIterator::row_iterator tr = tul.rowIterator();
967  // FIXME: replace with reduceLineBurtFilter()
968  resamplingConvolveLine(sr, sr+wold, src, tr, tr+wnew, tmp.accessor(),
969  kernels, mapCoordinate);
970  }
971 
972  tul = tmp.upperLeft();
973 
974  for(; dul.x < dlr.x; ++dul.x, ++tul.x)
975  {
976  typename DestIterator::column_iterator dc = dul.columnIterator();
977  typename TmpIterator::column_iterator tc = tul.columnIterator();
978  resamplingConvolveLine(tc, tc+hold, tmp.accessor(), dc, dc+hnew, dest,
979  kernels, mapCoordinate);
980  }
981 }
982 
983 template <class SrcIterator, class SrcAccessor,
984  class DestIterator, class DestAccessor>
985 inline
986 void pyramidReduceBurtFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
987  triple<DestIterator, DestIterator, DestAccessor> dest,
988  double centerValue = 0.4)
989 {
990  pyramidReduceBurtFilter(src.first, src.second, src.third,
991  dest.first, dest.second, dest.third, centerValue);
992 }
993 
994 template <class Image, class Alloc>
995 inline
996 void pyramidReduceBurtFilter(ImagePyramid<Image, Alloc> & pyramid, int fromLevel, int toLevel,
997  double centerValue = 0.4)
998 {
999  vigra_precondition(fromLevel < toLevel,
1000  "pyramidReduceBurtFilter(): fromLevel must be smaller than toLevel.");
1001  vigra_precondition(pyramid.lowestLevel() <= fromLevel && toLevel <= pyramid.highestLevel(),
1002  "pyramidReduceBurtFilter(): fromLevel and toLevel must be between the lowest and highest pyramid levels (inclusive).");
1003 
1004  for(int i=fromLevel+1; i <= toLevel; ++i)
1005  pyramidReduceBurtFilter(srcImageRange(pyramid[i-1]), destImageRange(pyramid[i]), centerValue);
1006 }
1007 
1008 /** \brief Two-fold up-sampling for image pyramid reconstruction.
1009 
1010  Sorry, no \ref detailedDocumentation() available yet.
1011 
1012  <b> Declarations:</b>
1013 
1014  <b>\#include</b> <vigra/resampling_convolution.hxx><br>
1015  Namespace: vigra
1016 
1017  pass 2D array views:
1018  \code
1019  namespace vigra {
1020  template <class SrcIterator, class SrcAccessor,
1021  class DestIterator, class DestAccessor>
1022  void pyramidExpandBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1023  DestIterator dul, DestIterator dlr, DestAccessor dest,
1024  double centerValue = 0.4);
1025  }
1026  \endcode
1027 
1028  \deprecatedAPI{pyramidExpandBurtFilter}
1029  pass \ref ImageIterators and \ref DataAccessors :
1030  \code
1031  namespace vigra {
1032  template <class SrcIterator, class SrcAccessor,
1033  class DestIterator, class DestAccessor>
1034  void pyramidExpandBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1035  DestIterator dul, DestIterator dlr, DestAccessor dest,
1036  double centerValue = 0.4);
1037  }
1038  \endcode
1039  use argument objects in conjunction with \ref ArgumentObjectFactories :
1040  \code
1041  namespace vigra {
1042  template <class SrcIterator, class SrcAccessor,
1043  class DestIterator, class DestAccessor>
1044  void pyramidExpandBurtFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1045  triple<DestIterator, DestIterator, DestAccessor> dest,
1046  double centerValue = 0.4);
1047  }
1048  \endcode
1049  \deprecatedEnd
1050 
1051  use a \ref vigra::ImagePyramid :
1052  \code
1053  namespace vigra {
1054  template <class Image, class Alloc>
1055  void pyramidExpandBurtFilter(ImagePyramid<Image, Alloc> & pyramid, int fromLevel, int toLevel,
1056  double centerValue = 0.4);
1057  }
1058  \endcode
1059 */
1060 doxygen_overloaded_function(template <...> void pyramidExpandBurtFilter)
1061 
1062 template <class SrcIterator, class SrcAccessor,
1063  class DestIterator, class DestAccessor>
1064 void pyramidExpandBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1065  DestIterator dul, DestIterator dlr, DestAccessor dest,
1066  double centerValue = 0.4)
1067 {
1068  vigra_precondition(0.25 <= centerValue && centerValue <= 0.5,
1069  "pyramidExpandBurtFilter(): centerValue must be between 0.25 and 0.5.");
1070 
1071  int wold = slr.x - sul.x;
1072  int wnew = dlr.x - dul.x;
1073  int hold = slr.y - sul.y;
1074  int hnew = dlr.y - dul.y;
1075 
1076  vigra_precondition(wold == (wnew + 1) / 2 && hold == (hnew + 1) / 2,
1077  "pyramidExpandBurtFilter(): oldSize = ceil(newSize / 2) required.");
1078 
1079  vigra_precondition(wold == (wnew + 1) / 2 && hold == (hnew + 1) / 2,
1080  "pyramidExpandBurtFilter(): oldSize = ceil(newSize / 2) required.");
1081 
1082  Rational<int> samplingRatio(2), offset(0);
1083  resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);
1084 
1085  ArrayVector<Kernel1D<double> > kernels(2);
1086  kernels[0].initExplicitly(-1, 1) = 0.5 - centerValue, 2.0*centerValue, 0.5 - centerValue;
1087  kernels[1].initExplicitly(-1, 0) = 0.5, 0.5;
1088 
1089  typedef typename
1090  NumericTraits<typename SrcAccessor::value_type>::RealPromote
1091  TmpType;
1092  typedef BasicImage<TmpType> TmpImage;
1093  typedef typename TmpImage::traverser TmpIterator;
1094 
1095  BasicImage<TmpType> tmp(wnew, hold);
1096 
1097  TmpIterator tul = tmp.upperLeft();
1098 
1099  for(; sul.y < slr.y; ++sul.y, ++tul.y)
1100  {
1101  typename SrcIterator::row_iterator sr = sul.rowIterator();
1102  typename TmpIterator::row_iterator tr = tul.rowIterator();
1103  // FIXME: replace with expandLineBurtFilter()
1104  resamplingConvolveLine(sr, sr+wold, src, tr, tr+wnew, tmp.accessor(),
1105  kernels, mapCoordinate);
1106  }
1107 
1108  tul = tmp.upperLeft();
1109 
1110  for(; dul.x < dlr.x; ++dul.x, ++tul.x)
1111  {
1112  typename DestIterator::column_iterator dc = dul.columnIterator();
1113  typename TmpIterator::column_iterator tc = tul.columnIterator();
1114  resamplingConvolveLine(tc, tc+hold, tmp.accessor(), dc, dc+hnew, dest,
1115  kernels, mapCoordinate);
1116  }
1117 }
1118 
1119 template <class SrcIterator, class SrcAccessor,
1120  class DestIterator, class DestAccessor>
1121 inline
1122 void pyramidExpandBurtFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1123  triple<DestIterator, DestIterator, DestAccessor> dest,
1124  double centerValue = 0.4)
1125 {
1126  pyramidExpandBurtFilter(src.first, src.second, src.third,
1127  dest.first, dest.second, dest.third, centerValue);
1128 }
1129 
1130 
1131 template <class Image, class Alloc>
1132 inline
1133 void pyramidExpandBurtFilter(ImagePyramid<Image, Alloc> & pyramid, int fromLevel, int toLevel,
1134  double centerValue = 0.4)
1135 {
1136  vigra_precondition(fromLevel > toLevel,
1137  "pyramidExpandBurtFilter(): fromLevel must be larger than toLevel.");
1138  vigra_precondition(pyramid.lowestLevel() <= toLevel && fromLevel <= pyramid.highestLevel(),
1139  "pyramidExpandBurtFilter(): fromLevel and toLevel must be between the lowest and highest pyramid levels (inclusive).");
1140 
1141  for(int i=fromLevel-1; i >= toLevel; --i)
1142  pyramidExpandBurtFilter(srcImageRange(pyramid[i+1]), destImageRange(pyramid[i]), centerValue);
1143 }
1144 
1145 /** \brief Create a Laplacian pyramid.
1146 
1147  Sorry, no \ref detailedDocumentation() available yet.
1148 
1149  <b>\#include</b> <vigra/resampling_convolution.hxx><br>
1150  Namespace: vigra
1151 */
1152 template <class Image, class Alloc>
1153 inline void
1154 pyramidReduceBurtLaplacian(ImagePyramid<Image, Alloc> & pyramid, int fromLevel, int toLevel,
1155  double centerValue = 0.4)
1156 {
1157  using namespace functor;
1158 
1159  pyramidReduceBurtFilter(pyramid, fromLevel, toLevel, centerValue);
1160  for(int i=fromLevel; i < toLevel; ++i)
1161  {
1162  typename ImagePyramid<Image, Alloc>::value_type tmpImage(pyramid[i].size());
1163  pyramidExpandBurtFilter(srcImageRange(pyramid[i+1]), destImageRange(tmpImage), centerValue);
1164  combineTwoImages(srcImageRange(tmpImage), srcImage(pyramid[i]), destImage(pyramid[i]),
1165  Arg1() - Arg2());
1166  }
1167 }
1168 
1169 /** \brief Reconstruct a Laplacian pyramid.
1170 
1171  Sorry, no \ref detailedDocumentation() available yet.
1172 
1173  <b>\#include</b> <vigra/resampling_convolution.hxx><br>
1174  Namespace: vigra
1175 */
1176 template <class Image, class Alloc>
1177 inline void
1178 pyramidExpandBurtLaplacian(ImagePyramid<Image, Alloc> & pyramid, int fromLevel, int toLevel,
1179  double centerValue = 0.4)
1180 {
1181  using namespace functor;
1182 
1183  vigra_precondition(fromLevel > toLevel,
1184  "pyramidExpandBurtLaplacian(): fromLevel must be larger than toLevel.");
1185  vigra_precondition(pyramid.lowestLevel() <= toLevel && fromLevel <= pyramid.highestLevel(),
1186  "pyramidExpandBurtLaplacian(): fromLevel and toLevel must be between the lowest and highest pyramid levels (inclusive).");
1187 
1188  for(int i=fromLevel-1; i >= toLevel; --i)
1189  {
1190  typename ImagePyramid<Image, Alloc>::value_type tmpImage(pyramid[i].size());
1191  pyramidExpandBurtFilter(srcImageRange(pyramid[i+1]), destImageRange(tmpImage), centerValue);
1192  combineTwoImages(srcImageRange(tmpImage), srcImage(pyramid[i]), destImage(pyramid[i]),
1193  Arg1() - Arg2());
1194  }
1195 }
1196 
1197 //@}
1198 
1199 } // namespace vigra
1200 
1201 
1202 #endif /* VIGRA_RESAMPLING_CONVOLUTION_HXX */
void resamplingConvolveX(...)
Apply a resampling filter in the x-direction.
void pyramidExpandBurtFilter(...)
Two-fold up-sampling for image pyramid reconstruction.
ImageType value_type
Definition: imagecontainer.hxx:476
void pyramidExpandBurtLaplacian(ImagePyramid< Image, Alloc > &pyramid, int fromLevel, int toLevel, double centerValue=0.4)
Reconstruct a Laplacian pyramid.
Definition: resampling_convolution.hxx:1178
Definition: accessor.hxx:43
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector&#39;s elements
Definition: tinyvector.hxx:1871
void resamplingConvolveImage(...)
Apply two separable resampling filters successively, the first in x-direction, the second in y-direct...
IntType lcm(IntType n, IntType m)
Definition: rational.hxx:122
void combineTwoImages(...)
Combine two source images into destination image.
void pyramidReduceBurtFilter(...)
Two-fold down-sampling for image pyramid construction.
int highestLevel() const
Definition: imagecontainer.hxx:573
void resamplingConvolveY(...)
Apply a resampling filter in the y-direction.
int lowestLevel() const
Definition: imagecontainer.hxx:566
void pyramidReduceBurtLaplacian(ImagePyramid< Image, Alloc > &pyramid, int fromLevel, int toLevel, double centerValue=0.4)
Create a Laplacian pyramid.
Definition: resampling_convolution.hxx:1154
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition: fixedpoint.hxx:675
void resamplingConvolveLine(...)
Performs a 1-dimensional resampling convolution of the source signal using the given set of kernels...
Class template for logarithmically tapering image pyramids.
Definition: imagecontainer.hxx:465
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition: fixedpoint.hxx:667

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

html generated using doxygen and Python
vigra 1.10.0