1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2004 by Ullrich Koethe */
4 /* Cognitive Systems Group, University of Hamburg, Germany */
5 /* */
6 /* This file is part of the VIGRA computer vision library. */
7 /* ( Version 1.5.0, Dec 07 2006 ) */
8 /* The VIGRA Website is */
9 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */
10 /* Please direct questions, bug reports, and contributions to */
11 /* koethe@informatik.uni-hamburg.de or */
12 /* vigra@kogs1.informatik.uni-hamburg.de */
13 /* */
14 /* Permission is hereby granted, free of charge, to any person */
15 /* obtaining a copy of this software and associated documentation */
16 /* files (the "Software"), to deal in the Software without */
17 /* restriction, including without limitation the rights to use, */
18 /* copy, modify, merge, publish, distribute, sublicense, and/or */
19 /* sell copies of the Software, and to permit persons to whom the */
20 /* Software is furnished to do so, subject to the following */
21 /* conditions: */
22 /* */
23 /* The above copyright notice and this permission notice shall be */
24 /* included in all copies or substantial portions of the */
25 /* Software. */
26 /* */
27 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
28 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
29 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
30 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
31 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
32 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
33 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
34 /* OTHER DEALINGS IN THE SOFTWARE. */
35 /* */
36 /************************************************************************/
37
38 #ifndef VIGRA_RESAMPLING_CONVOLUTION_HXX
39 #define VIGRA_RESAMPLING_CONVOLUTION_HXX
40
41 #include <cmath>
42 #include "stdimage.hxx"
43 #include "array_vector.hxx"
44 #include "rational.hxx"
45 #include "functortraits.hxx"
46
47 namespace vigra {
48
49 namespace resampling_detail
50 {
51
52 struct MapTargetToSourceCoordinate
53 {
54 MapTargetToSourceCoordinate(Rational<int> const & samplingRatio,
55 Rational<int> const & offset)
56 : a(samplingRatio.denominator()*offset.denominator()),
57 b(samplingRatio.numerator()*offset.numerator()),
58 c(samplingRatio.numerator()*offset.denominator())
59 {}
60
61 // the following funcions are more efficient realizations of:
62 // rational_cast<T>(i / samplingRatio + offset);
63 // we need efficiency because this may be called in the inner loop
64
65 int operator()(int i) const
66 {
67 return (i * a + b) / c;
68 }
69
70 double toDouble(int i) const
71 {
72 return double(i * a + b) / c;
73 }
74
75 Rational<int> toRational(int i) const
76 {
77 return Rational<int>(i * a + b, c);
78 }
79
EnsureNewOldDoclet()80 int a, b, c;
81 };
82
83 } // namespace resampling_detail
84
85 template <>
86 class FunctorTraits<resampling_detail::MapTargetToSourceCoordinate>
87 : public FunctorTraitsBase<resampling_detail::MapTargetToSourceCoordinate>
88 {
89 public:
90 typedef VigraTrueType isUnaryFunctor;
generateSample(File testSrc)91 };
92
93 template <class SrcIter, class SrcAcc,
94 class DestIter, class DestAcc,
95 class KernelArray,
96 class Functor>
97 void
98 resamplingConvolveLine(SrcIter s, SrcIter send, SrcAcc src,
99 DestIter d, DestIter dend, DestAcc dest,
100 KernelArray const & kernels,
main(String... args)101 Functor mapTargetToSourceCoordinate)
102 {
103 typedef typename
104 NumericTraits<typename SrcAcc::value_type>::RealPromote
105 TmpType;
106 typedef typename KernelArray::value_type Kernel;
107 typedef typename Kernel::const_iterator KernelIter;
108
109 int wo = send - s;
110 int wn = dend - d;
111 int wo2 = 2*wo - 2;
112
113 int i;
114 typename KernelArray::const_iterator kernel = kernels.begin();
115 for(i=0; i<wn; ++i, ++d, ++kernel)
116 {
117 // use the kernels periodically
118 if(kernel == kernels.end())
119 kernel = kernels.begin();
120
121 // calculate current target point into source location
122 int is = mapTargetToSourceCoordinate(i);
123
124 TmpType sum = NumericTraits<TmpType>::zero();
125
126 int lbound = is - kernel->right(),
127 hbound = is - kernel->left();
128
129 KernelIter k = kernel->center() + kernel->right();
130 if(lbound < 0 || hbound >= wo)
131 {
132 vigra_precondition(-lbound < wo && wo2 - hbound >= 0,
133 "resamplingConvolveLine(): kernel or offset larger than image.");
134 for(int m=lbound; m <= hbound; ++m, --k)
135 {
136 int mm = (m < 0) ?
137 -m :
138 (m >= wo) ?
139 wo2 - m :
140 m;
141 sum += *k * src(s, mm);
142 }
143 }
144 else
145 {
146 SrcIter ss = s + lbound;
147 SrcIter ssend = s + hbound;
148
149 for(; ss <= ssend; ++ss, --k)
150 {
151 sum += *k * src(ss);
152 }
153 }
154
155 dest.set(sum, d);
156 }
157 }
158
getAllowedLocations()159 template <class Kernel, class MapCoordinate, class KernelArray>
160 void
161 createResamplingKernels(Kernel const & kernel,
162 MapCoordinate const & mapCoordinate, KernelArray & kernels)
163 {
164 for(unsigned int idest = 0; idest < kernels.size(); ++idest)
165 {
166 int isrc = mapCoordinate(idest);
167 double idsrc = mapCoordinate.toDouble(idest);
168 double offset = idsrc - isrc;
169 double radius = kernel.radius();
170 int left = int(ceil(-radius - offset));
171 int right = int(floor(radius - offset));
172 kernels[idest].initExplicitly(left, right);
173
174 double x = left + offset;
175 for(int i = left; i <= right; ++i, ++x)
176 kernels[idest][i] = kernel(x);
177 kernels[idest].normalize(1.0, kernel.derivativeOrder(), offset);
178 }
179 }
180
181 /** \addtogroup ResamplingConvolutionFilters Resampling Convolution Filters
182
183 These functions implement the convolution operation when the source and target images
184 have different sizes. This is realized by accessing a continous kernel at the
185 appropriate non-integer positions. The technique is, for example, described in
186 D. Schumacher: <i>General Filtered Image Rescaling</i>, in: Graphics Gems III,
187 Academic Press, 1992.
188 */
189 //@{
190
191 /********************************************************/
192 /* */
193 /* resamplingConvolveX */
194 /* */
195 /********************************************************/
196
197 /** \brief Apply a resampling filter in the x-direction.
198
199 This function implements a convolution operation in x-direction
200 (i.e. applies a 1D filter to every row) where the width of the source
201 and destination images differ. This is typically used to avoid aliasing if
202 the image is scaled down, or to interpolate smoothly if the image is scaled up.
203 The target coordinates are transformed into source coordinates by
204
205 \code
206 xsource = (xtarget - offset) / samplingRatio
207 \endcode
208
209 The <tt>samplingRatio</tt> and <tt>offset</tt> must be given as \ref vigra::Rational
210 in order to avoid rounding errors in this transformation. It is required that for all
211 pixels of the target image, <tt>xsource</tt> remains within the range of the source
212 image (i.e. <tt>0 <= xsource <= sourceWidth-1</tt>. Since <tt>xsource</tt> is
213 in general not an integer, the <tt>kernel</tt> must be a functor that can be accessed at
214 arbitrary (<tt>double</tt>) coordinates. It must also provide a member function <tt>radius()</tt>
215 which specifies the support (non-zero interval) of the kernel. VIGRA already
216 provides a number of suitable functors, e.g. \ref vigra::Gaussian, \ref vigra::BSpline
217 \ref vigra::CatmullRomSpline, and \ref vigra::CoscotFunction. The function
218 \ref resizeImageSplineInterpolation() is implemented by means resamplingConvolveX() and
219 resamplingConvolveY().
220
221 <b> Declarations:</b>
222
223 pass arguments explicitly:
224 \code
225 namespace vigra {
226 template <class SrcIter, class SrcAcc,
227 class DestIter, class DestAcc,
228 class Kernel>
229 void
230 resamplingConvolveX(SrcIter sul, SrcIter slr, SrcAcc src,
231 DestIter dul, DestIter dlr, DestAcc dest,
232 Kernel const & kernel,
233 Rational<int> const & samplingRatio, Rational<int> const & offset);
234 }
235 \endcode
236
237
238 use argument objects in conjunction with \ref ArgumentObjectFactories:
239 \code
240 namespace vigra {
241 template <class SrcIter, class SrcAcc,
242 class DestIter, class DestAcc,
243 class Kernel>
244 void
245 resamplingConvolveX(triple<SrcIter, SrcIter, SrcAcc> src,
246 triple<DestIter, DestIter, DestAcc> dest,
247 Kernel const & kernel,
248 Rational<int> const & samplingRatio, Rational<int> const & offset);
249 }
250 \endcode
251
252 <b> Usage:</b>
253
254 <b>\#include</b> "<a href="resampling__convolution_8hxx-source.html">vigra/resampling_convolution.hxx</a>"
255
256
257 \code
258 Rational<int> ratio(2), offset(0);
259
260 FImage src(w,h),
261 dest(rational_cast<int>(ratio*w), h);
262
263 float sigma = 2.0;
264 Gaussian<float> smooth(sigma);
265 ...
266
267 // simpultaneously enlarge and smooth source image
268 resamplingConvolveX(srcImageRange(src), destImageRange(dest),
269 smooth, ratio, offset);
270 \endcode
271
272 <b> Required Interface:</b>
273
274 \code
275 Kernel kernel;
276 int kernelRadius = kernel.radius();
277 double x = ...; // must be <= radius()
278 double value = kernel(x);
279 \endcode
280 */
281 template <class SrcIter, class SrcAcc,
282 class DestIter, class DestAcc,
283 class Kernel>
284 void
285 resamplingConvolveX(SrcIter sul, SrcIter slr, SrcAcc src,
286 DestIter dul, DestIter dlr, DestAcc dest,
287 Kernel const & kernel,
288 Rational<int> const & samplingRatio, Rational<int> const & offset)
289 {
290 int wold = slr.x - sul.x;
291 int wnew = dlr.x - dul.x;
292
293 vigra_precondition(!samplingRatio.is_inf() && samplingRatio > 0,
294 "resamplingConvolveX(): sampling ratio must be > 0 and < infinity");
295 vigra_precondition(!offset.is_inf(),
296 "resamplingConvolveX(): offset must be < infinity");
297
298 int period = lcm(samplingRatio.numerator(), samplingRatio.denominator());
299 resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);
300
301 ArrayVector<Kernel1D<double> > kernels(period);
302
303 createResamplingKernels(kernel, mapCoordinate, kernels);
304
305 for(; sul.y < slr.y; ++sul.y, ++dul.y)
306 {
307 typename SrcIter::row_iterator sr = sul.rowIterator();
308 typename DestIter::row_iterator dr = dul.rowIterator();
309 resamplingConvolveLine(sr, sr+wold, src, dr, dr+wnew, dest,
310 kernels, mapCoordinate);
311 }
312 }
313
314 template <class SrcIter, class SrcAcc,
315 class DestIter, class DestAcc,
316 class Kernel>
317 inline void
318 resamplingConvolveX(triple<SrcIter, SrcIter, SrcAcc> src,
319 triple<DestIter, DestIter, DestAcc> dest,
320 Kernel const & kernel,
321 Rational<int> const & samplingRatio, Rational<int> const & offset)
322 {
323 resamplingConvolveX(src.first, src.second, src.third,
324 dest.first, dest.second, dest.third,
325 kernel, samplingRatio, offset);
326 }
327
328 /********************************************************/
329 /* */
330 /* resamplingConvolveY */
331 /* */
332 /********************************************************/
333
334 /** \brief Apply a resampling filter in the y-direction.
335
336 This function implements a convolution operation in y-direction
337 (i.e. applies a 1D filter to every column) where the height of the source
338 and destination images differ. This is typically used to avoid aliasing if
339 the image is scaled down, or to interpolate smoothly if the image is scaled up.
340 The target coordinates are transformed into source coordinates by
341
342 \code
343 ysource = (ytarget - offset) / samplingRatio
344 \endcode
345
346 The <tt>samplingRatio</tt> and <tt>offset</tt> must be given as \ref vigra::Rational
347 in order to avoid rounding errors in this transformation. It is required that for all
348 pixels of the target image, <tt>ysource</tt> remains within the range of the source
349 image (i.e. <tt>0 <= ysource <= sourceHeight-1</tt>. Since <tt>ysource</tt> is
350 in general not an integer, the <tt>kernel</tt> must be a functor that can be accessed at
351 arbitrary (<tt>double</tt>) coordinates. It must also provide a member function <tt>radius()</tt>
352 which specifies the support (non-zero interval) of the kernel. VIGRA already
353 provides a number of suitable functors, e.g. \ref vigra::Gaussian, \ref vigra::BSpline
354 \ref vigra::CatmullRomSpline, and \ref vigra::CoscotFunction. The function
355 \ref resizeImageSplineInterpolation() is implemented by means resamplingConvolveX() and
356 resamplingConvolveY().
357
358 <b> Declarations:</b>
359
360 pass arguments explicitly:
361 \code
362 namespace vigra {
363 template <class SrcIter, class SrcAcc,
364 class DestIter, class DestAcc,
365 class Kernel>
366 void
367 resamplingConvolveY(SrcIter sul, SrcIter slr, SrcAcc src,
368 DestIter dul, DestIter dlr, DestAcc dest,
369 Kernel const & kernel,
370 Rational<int> const & samplingRatio, Rational<int> const & offset);
371 }
372 \endcode
373
374
375 use argument objects in conjunction with \ref ArgumentObjectFactories:
376 \code
377 namespace vigra {
378 template <class SrcIter, class SrcAcc,
379 class DestIter, class DestAcc,
380 class Kernel>
381 void
382 resamplingConvolveY(triple<SrcIter, SrcIter, SrcAcc> src,
383 triple<DestIter, DestIter, DestAcc> dest,
384 Kernel const & kernel,
385 Rational<int> const & samplingRatio, Rational<int> const & offset);
386 }
387 \endcode
388
389 <b> Usage:</b>
390
391 <b>\#include</b> "<a href="resampling__convolution_8hxx-source.html">vigra/resampling_convolution.hxx</a>"
392
393
394 \code
395 Rational<int> ratio(2), offset(0);
396
397 FImage src(w,h),
398 dest(w, rational_cast<int>(ratio*h));
399
400 float sigma = 2.0;
401 Gaussian<float> smooth(sigma);
402 ...
403
404 // simpultaneously enlarge and smooth source image
405 resamplingConvolveY(srcImageRange(src), destImageRange(dest),
406 smooth, ratio, offset);
407 \endcode
408
409 <b> Required Interface:</b>
410
411 \code
412 Kernel kernel;
413 int kernelRadius = kernel.radius();
414 double y = ...; // must be <= radius()
415 double value = kernel(y);
416 \endcode
417 */
418 template <class SrcIter, class SrcAcc,
419 class DestIter, class DestAcc,
420 class Kernel>
421 void
422 resamplingConvolveY(SrcIter sul, SrcIter slr, SrcAcc src,
423 DestIter dul, DestIter dlr, DestAcc dest,
424 Kernel const & kernel,
425 Rational<int> const & samplingRatio, Rational<int> const & offset)
426 {
427 int hold = slr.y - sul.y;
428 int hnew = dlr.y - dul.y;
429
430 vigra_precondition(!samplingRatio.is_inf() && samplingRatio > 0,
431 "resamplingConvolveY(): sampling ratio must be > 0 and < infinity");
432 vigra_precondition(!offset.is_inf(),
433 "resamplingConvolveY(): offset must be < infinity");
434
435 int period = lcm(samplingRatio.numerator(), samplingRatio.denominator());
436
437 resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);
438
439 ArrayVector<Kernel1D<double> > kernels(period);
440
441 createResamplingKernels(kernel, mapCoordinate, kernels);
442
443 for(; sul.x < slr.x; ++sul.x, ++dul.x)
444 {
445 typename SrcIter::column_iterator sc = sul.columnIterator();
446 typename DestIter::column_iterator dc = dul.columnIterator();
447 resamplingConvolveLine(sc, sc+hold, src, dc, dc+hnew, dest,
448 kernels, mapCoordinate);
449 }
450 }
451
452 template <class SrcIter, class SrcAcc,
453 class DestIter, class DestAcc,
454 class Kernel>
455 inline void
456 resamplingConvolveY(triple<SrcIter, SrcIter, SrcAcc> src,
457 triple<DestIter, DestIter, DestAcc> dest,
458 Kernel const & kernel,
459 Rational<int> const & samplingRatio, Rational<int> const & offset)
460 {
461 resamplingConvolveY(src.first, src.second, src.third,
462 dest.first, dest.second, dest.third,
463 kernel, samplingRatio, offset);
464 }
465
466 /********************************************************/
467 /* */
468 /* resamplingConvolveImage */
469 /* */
470 /********************************************************/
471
472 /** \brief Apply two separable resampling filters successively, the first in x-direction,
473 the second in y-direction.
474
475 This function is a shorthand for the concatenation of a call to
476 \link ResamplingConvolutionFilters#resamplingConvolveX resamplingConvolveX\endlink()
477 and \link ResamplingConvolutionFilters#resamplingConvolveY resamplingConvolveY\endlink()
478 with the given kernels. See there for detailed documentation.
479
480 <b> Declarations:</b>
481
482 pass arguments explicitly:
483 \code
484 namespace vigra {
485 template <class SrcIterator, class SrcAccessor,
486 class DestIterator, class DestAccessor,
487 class KernelX, class KernelY>
488 void resamplingConvolveImage(SrcIterator sul,SrcIterator slr, SrcAccessor src,
489 DestIterator dul, DestIterator dlr, DestAccessor dest,
490 KernelX const & kx,
491 Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
492 KernelY const & ky,
493 Rational<int> const & samplingRatioY, Rational<int> const & offsetY);
494 }
495 \endcode
496
497
498 use argument objects in conjunction with \ref ArgumentObjectFactories:
499 \code
500 namespace vigra {
501 template <class SrcIterator, class SrcAccessor,
502 class DestIterator, class DestAccessor,
503 class KernelX, class KernelY>
504 void
505 resamplingConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
506 triple<DestIterator, DestIterator, DestAccessor> dest,
507 KernelX const & kx,
508 Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
509 KernelY const & ky,
510 Rational<int> const & samplingRatioY, Rational<int> const & offsetY);
511 }
512 \endcode
513
514 <b> Usage:</b>
515
516 <b>\#include</b> "<a href="resampling__convolution_8hxx-source.html">vigra/resampling_convolution.hxx</a>"
517
518
519 \code
520 Rational<int> xratio(2), yratio(3), offset(0);
521
522 FImage src(w,h),
523 dest(rational_cast<int>(xratio*w), rational_cast<int>(yratio*h));
524
525 float sigma = 2.0;
526 Gaussian<float> smooth(sigma);
527 ...
528
529 // simpultaneously enlarge and smooth source image
530 resamplingConvolveImage(srcImageRange(src), destImageRange(dest),
531 smooth, xratio, offset,
532 smooth, yratio, offset);
533
534 \endcode
535
536 */
537 template <class SrcIterator, class SrcAccessor,
538 class DestIterator, class DestAccessor,
539 class KernelX, class KernelY>
540 void resamplingConvolveImage(SrcIterator sul,SrcIterator slr, SrcAccessor src,
541 DestIterator dul, DestIterator dlr, DestAccessor dest,
542 KernelX const & kx,
543 Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
544 KernelY const & ky,
545 Rational<int> const & samplingRatioY, Rational<int> const & offsetY)
546 {
547 typedef typename
548 NumericTraits<typename SrcAccessor::value_type>::RealPromote
549 TmpType;
550
551 BasicImage<TmpType> tmp(dlr.x - dul.x, slr.y - sul.y);
552
553 resamplingConvolveX(srcIterRange(sul, slr, src),
554 destImageRange(tmp),
555 kx, samplingRatioX, offsetX);
556 resamplingConvolveY(srcImageRange(tmp),
557 destIterRange(dul, dlr, dest),
558 ky, samplingRatioY, offsetY);
559 }
560
561 template <class SrcIterator, class SrcAccessor,
562 class DestIterator, class DestAccessor,
563 class KernelX, class KernelY>
564 inline void
565 resamplingConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
566 triple<DestIterator, DestIterator, DestAccessor> dest,
567 KernelX const & kx,
568 Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
569 KernelY const & ky,
570 Rational<int> const & samplingRatioY, Rational<int> const & offsetY)
571 {
572 resamplingConvolveImage(src.first, src.second, src.third,
573 dest.first, dest.second, dest.third,
574 kx, samplingRatioX, offsetX,
575 ky, samplingRatioY, offsetY);
576 }
577
578 } // namespace vigra
579
580
581 #endif /* VIGRA_RESAMPLING_CONVOLUTION_HXX */
582