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

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

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

html generated using doxygen and Python
vigra 1.8.0 (Wed Sep 26 2012)