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 {
MapTargetToSourceCoordinatevigra::resampling_detail::MapTargetToSourceCoordinate56 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
operator ()vigra::resampling_detail::MapTargetToSourceCoordinate67 int operator()(int i) const
68 {
69 return (i * a + b) / c;
70 }
71
toDoublevigra::resampling_detail::MapTargetToSourceCoordinate72 double toDouble(int i) const
73 {
74 return double(i * a + b) / c;
75 }
76
toRationalvigra::resampling_detail::MapTargetToSourceCoordinate77 Rational<int> toRational(int i) const
78 {
79 return Rational<int>(i * a + b, c);
80 }
81
isExpand2vigra::resampling_detail::MapTargetToSourceCoordinate82 bool isExpand2() const
83 {
84 return a == 1 && b == 0 && c == 2;
85 }
86
isReduce2vigra::resampling_detail::MapTargetToSourceCoordinate87 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
resamplingExpandLine2(SrcIter s,SrcIter send,SrcAcc src,DestIter d,DestIter dend,DestAcc dest,KernelArray const & kernels)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
resamplingReduceLine2(SrcIter s,SrcIter send,SrcAcc src,DestIter d,DestIter dend,DestAcc dest,KernelArray const & kernels)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 */
doxygen_overloaded_function(template<...> void resamplingConvolveLine)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
createResamplingKernels(Kernel const & kernel,MapCoordinate const & mapCoordinate,KernelArray & kernels)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 */
doxygen_overloaded_function(template<...> void resamplingConvolveX)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
resamplingConvolveX(triple<SrcIter,SrcIter,SrcAcc> src,triple<DestIter,DestIter,DestAcc> dest,Kernel const & kernel,Rational<int> const & samplingRatio,Rational<int> const & offset)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
resamplingConvolveX(MultiArrayView<2,T1,S1> const & src,MultiArrayView<2,T2,S2> dest,Kernel const & kernel,Rational<int> const & samplingRatio,Rational<int> const & offset)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 */
doxygen_overloaded_function(template<...> void resamplingConvolveY)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
resamplingConvolveY(triple<SrcIter,SrcIter,SrcAccessor> src,triple<DestIter,DestIter,DestAccessor> dest,Kernel const & kernel,Rational<int> const & samplingRatio,Rational<int> const & offset)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
resamplingConvolveY(MultiArrayView<2,T1,S1> const & src,MultiArrayView<2,T2,S2> dest,Kernel const & kernel,Rational<int> const & samplingRatio,Rational<int> const & offset)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 */
doxygen_overloaded_function(template<...> void resamplingConvolveImage)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
resamplingConvolveImage(triple<SrcIterator,SrcIterator,SrcAccessor> src,triple<DestIterator,DestIterator,DestAccessor> dest,KernelX const & kx,Rational<int> const & samplingRatioX,Rational<int> const & offsetX,KernelY const & ky,Rational<int> const & samplingRatioY,Rational<int> const & offsetY)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
resamplingConvolveImage(MultiArrayView<2,T1,S1> const & src,MultiArrayView<2,T2,S2> dest,KernelX const & kx,Rational<int> const & samplingRatioX,Rational<int> const & offsetX,KernelY const & ky,Rational<int> const & samplingRatioY,Rational<int> const & offsetY)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 This function implements the reduction by one resolution level (first signature)
876 or across several pyramid levels (last signature) of a Gaussian pyramid as described in
877
878 P. Burt, E. Adelson: <i>"The Laplacian Pyramid as a Compact Image Code"</i>, IEEE Trans. Communications, 9(4):532-540, 1983
879
880 That is, it applies the smoothing filter
881 \code
882 [0.25 - centerValue / 2.0, 0.25, centerValue, 0.25, 0.25 - centerValue / 2.0]
883 \endcode
884 to the source image and then copies all pixels with even coordinates to the destination
885 image. The destination image shape must be <tt>dest_shape = ceil(src_shape / 2.0)</tt>.
886 <tt>centerValue</tt> must be between 0.25 and 0.5 and determines the strength of smoothing
887 (bigger values correspond to less smoothing). If <tt>toLevel - fromLevel > 1</tt> in the
888 pyramid variant of the function, this process is repeated until <tt>toLevel</tt> is
889 reached.
890
891 Typically, this functions is used in connection with a \ref vigra::ImagePyramid
892 (last signature below) to perform several levels of reduction in one go.
893
894 <b> Declarations:</b>
895
896 <b>\#include</b> \<vigra/resampling_convolution.hxx\><br>
897 Namespace: vigra
898
899 pass 2D array views:
900 \code
901 namespace vigra {
902 template <class SrcIterator, class SrcAccessor,
903 class DestIterator, class DestAccessor>
904 void pyramidReduceBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
905 DestIterator dul, DestIterator dlr, DestAccessor dest,
906 double centerValue = 0.4);
907 }
908 \endcode
909
910 \deprecatedAPI{pyramidReduceBurtFilter}
911 pass \ref ImageIterators and \ref DataAccessors :
912 \code
913 namespace vigra {
914 template <class SrcIterator, class SrcAccessor,
915 class DestIterator, class DestAccessor>
916 void pyramidReduceBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
917 DestIterator dul, DestIterator dlr, DestAccessor dest,
918 double centerValue = 0.4);
919 }
920 \endcode
921 use argument objects in conjunction with \ref ArgumentObjectFactories :
922 \code
923 namespace vigra {
924 template <class SrcIterator, class SrcAccessor,
925 class DestIterator, class DestAccessor>
926 void pyramidReduceBurtFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
927 triple<DestIterator, DestIterator, DestAccessor> dest,
928 double centerValue = 0.4);
929 }
930 \endcode
931 \deprecatedEnd
932
933 use a \ref vigra::ImagePyramid :
934 \code
935 namespace vigra {
936 template <class Image, class Alloc>
937 void pyramidReduceBurtFilter(ImagePyramid<Image, Alloc> & pyramid,
938 int fromLevel, int toLevel,
939 double centerValue = 0.4);
940 }
941 \endcode
942 */
doxygen_overloaded_function(template<...> void pyramidReduceBurtFilter)943 doxygen_overloaded_function(template <...> void pyramidReduceBurtFilter)
944
945 template <class SrcIterator, class SrcAccessor,
946 class DestIterator, class DestAccessor>
947 void pyramidReduceBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
948 DestIterator dul, DestIterator dlr, DestAccessor dest,
949 double centerValue = 0.4)
950 {
951 vigra_precondition(0.25 <= centerValue && centerValue <= 0.5,
952 "pyramidReduceBurtFilter(): centerValue must be between 0.25 and 0.5.");
953
954 int wold = slr.x - sul.x;
955 int wnew = dlr.x - dul.x;
956 int hold = slr.y - sul.y;
957 int hnew = dlr.y - dul.y;
958
959 vigra_precondition(wnew == (wold + 1) / 2 && hnew == (hold + 1) / 2,
960 "pyramidReduceBurtFilter(): destSize = ceil(srcSize / 2) required.");
961
962 Rational<int> samplingRatio(1,2), offset(0);
963 resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);
964
965 ArrayVector<Kernel1D<double> > kernels(1);
966 kernels[0].initExplicitly(-2, 2) = 0.25 - centerValue / 2.0, 0.25, centerValue, 0.25, 0.25 - centerValue / 2.0;
967
968 typedef typename
969 NumericTraits<typename SrcAccessor::value_type>::RealPromote
970 TmpType;
971 typedef BasicImage<TmpType> TmpImage;
972 typedef typename TmpImage::traverser TmpIterator;
973
974 BasicImage<TmpType> tmp(wnew, hold);
975
976 TmpIterator tul = tmp.upperLeft();
977
978 for(; sul.y < slr.y; ++sul.y, ++tul.y)
979 {
980 typename SrcIterator::row_iterator sr = sul.rowIterator();
981 typename TmpIterator::row_iterator tr = tul.rowIterator();
982 // FIXME: replace with reduceLineBurtFilter()
983 resamplingConvolveLine(sr, sr+wold, src, tr, tr+wnew, tmp.accessor(),
984 kernels, mapCoordinate);
985 }
986
987 tul = tmp.upperLeft();
988
989 for(; dul.x < dlr.x; ++dul.x, ++tul.x)
990 {
991 typename DestIterator::column_iterator dc = dul.columnIterator();
992 typename TmpIterator::column_iterator tc = tul.columnIterator();
993 resamplingConvolveLine(tc, tc+hold, tmp.accessor(), dc, dc+hnew, dest,
994 kernels, mapCoordinate);
995 }
996 }
997
998 template <class SrcIterator, class SrcAccessor,
999 class DestIterator, class DestAccessor>
1000 inline
pyramidReduceBurtFilter(triple<SrcIterator,SrcIterator,SrcAccessor> src,triple<DestIterator,DestIterator,DestAccessor> dest,double centerValue=0.4)1001 void pyramidReduceBurtFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1002 triple<DestIterator, DestIterator, DestAccessor> dest,
1003 double centerValue = 0.4)
1004 {
1005 pyramidReduceBurtFilter(src.first, src.second, src.third,
1006 dest.first, dest.second, dest.third, centerValue);
1007 }
1008
1009 template <class Image, class Alloc>
1010 inline
pyramidReduceBurtFilter(ImagePyramid<Image,Alloc> & pyramid,int fromLevel,int toLevel,double centerValue=0.4)1011 void pyramidReduceBurtFilter(ImagePyramid<Image, Alloc> & pyramid,
1012 int fromLevel, int toLevel,
1013 double centerValue = 0.4)
1014 {
1015 vigra_precondition(fromLevel < toLevel,
1016 "pyramidReduceBurtFilter(): fromLevel must be smaller than toLevel.");
1017 vigra_precondition(pyramid.lowestLevel() <= fromLevel && toLevel <= pyramid.highestLevel(),
1018 "pyramidReduceBurtFilter(): fromLevel and toLevel must be between the lowest and highest pyramid levels (inclusive).");
1019
1020 for(int i=fromLevel+1; i <= toLevel; ++i)
1021 pyramidReduceBurtFilter(srcImageRange(pyramid[i-1]), destImageRange(pyramid[i]), centerValue);
1022 }
1023
1024 /** \brief Two-fold up-sampling for image pyramid reconstruction.
1025
1026 This function implements the expansion by one resolution level (first signature)
1027 or across several pyramid levels (last signature) of a Gaussian pyramid as described in
1028
1029 P. Burt, E. Adelson: <i>"The Laplacian Pyramid as a Compact Image Code"</i>, IEEE Trans. Communications, 9(4):532-540, 1983
1030
1031 That is, the function first places the pixel values of the low-resolution
1032 image at the even pixel coordinates of the high-resolution image (pixels with
1033 at least one odd coordinate are zero-initialized) and then applies the
1034 interpolation filter
1035 \code
1036 [0.5 - centerValue, 0.5, 2*centerValue, 0.5, 0.5 - centerValue]
1037 \endcode
1038 to the high-resolution image. The source image shape must be
1039 <tt>src_shape = ceil(dest_shape / 2.0)</tt>.
1040 <tt>centerValue</tt> must be between 0.25 and 0.5 and determines the sharpness
1041 of the interpolation (bigger values correspond to sharper images).
1042 If <tt>fromLevel - toLevel > 1</tt> in the pyramid variant of the function,
1043 this process is repeated until <tt>toLevel</tt> is reached.
1044
1045 Typically, this functions is used in connection with a \ref vigra::ImagePyramid
1046 (last signature below) to perform several levels of expansion in one go.
1047
1048 <b> Declarations:</b>
1049
1050 <b>\#include</b> \<vigra/resampling_convolution.hxx\><br>
1051 Namespace: vigra
1052
1053 pass 2D array views:
1054 \code
1055 namespace vigra {
1056 template <class SrcIterator, class SrcAccessor,
1057 class DestIterator, class DestAccessor>
1058 void pyramidExpandBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1059 DestIterator dul, DestIterator dlr, DestAccessor dest,
1060 double centerValue = 0.4);
1061 }
1062 \endcode
1063
1064 \deprecatedAPI{pyramidExpandBurtFilter}
1065 pass \ref ImageIterators and \ref DataAccessors :
1066 \code
1067 namespace vigra {
1068 template <class SrcIterator, class SrcAccessor,
1069 class DestIterator, class DestAccessor>
1070 void pyramidExpandBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1071 DestIterator dul, DestIterator dlr, DestAccessor dest,
1072 double centerValue = 0.4);
1073 }
1074 \endcode
1075 use argument objects in conjunction with \ref ArgumentObjectFactories :
1076 \code
1077 namespace vigra {
1078 template <class SrcIterator, class SrcAccessor,
1079 class DestIterator, class DestAccessor>
1080 void pyramidExpandBurtFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1081 triple<DestIterator, DestIterator, DestAccessor> dest,
1082 double centerValue = 0.4);
1083 }
1084 \endcode
1085 \deprecatedEnd
1086
1087 use a \ref vigra::ImagePyramid :
1088 \code
1089 namespace vigra {
1090 template <class Image, class Alloc>
1091 void pyramidExpandBurtFilter(ImagePyramid<Image, Alloc> & pyramid,
1092 int fromLevel, int toLevel,
1093 double centerValue = 0.4);
1094 }
1095 \endcode
1096 */
doxygen_overloaded_function(template<...> void pyramidExpandBurtFilter)1097 doxygen_overloaded_function(template <...> void pyramidExpandBurtFilter)
1098
1099 template <class SrcIterator, class SrcAccessor,
1100 class DestIterator, class DestAccessor>
1101 void pyramidExpandBurtFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1102 DestIterator dul, DestIterator dlr, DestAccessor dest,
1103 double centerValue = 0.4)
1104 {
1105 vigra_precondition(0.25 <= centerValue && centerValue <= 0.5,
1106 "pyramidExpandBurtFilter(): centerValue must be between 0.25 and 0.5.");
1107
1108 int wold = slr.x - sul.x;
1109 int wnew = dlr.x - dul.x;
1110 int hold = slr.y - sul.y;
1111 int hnew = dlr.y - dul.y;
1112
1113 vigra_precondition(wold == (wnew + 1) / 2 && hold == (hnew + 1) / 2,
1114 "pyramidExpandBurtFilter(): oldSize = ceil(newSize / 2) required.");
1115
1116 Rational<int> samplingRatio(2), offset(0);
1117 resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);
1118
1119 ArrayVector<Kernel1D<double> > kernels(2);
1120 kernels[0].initExplicitly(-1, 1) = 0.5 - centerValue, 2.0*centerValue, 0.5 - centerValue;
1121 kernels[1].initExplicitly(-1, 0) = 0.5, 0.5;
1122
1123 typedef typename
1124 NumericTraits<typename SrcAccessor::value_type>::RealPromote
1125 TmpType;
1126 typedef BasicImage<TmpType> TmpImage;
1127 typedef typename TmpImage::traverser TmpIterator;
1128
1129 BasicImage<TmpType> tmp(wnew, hold);
1130
1131 TmpIterator tul = tmp.upperLeft();
1132
1133 for(; sul.y < slr.y; ++sul.y, ++tul.y)
1134 {
1135 typename SrcIterator::row_iterator sr = sul.rowIterator();
1136 typename TmpIterator::row_iterator tr = tul.rowIterator();
1137 // FIXME: replace with expandLineBurtFilter()
1138 resamplingConvolveLine(sr, sr+wold, src, tr, tr+wnew, tmp.accessor(),
1139 kernels, mapCoordinate);
1140 }
1141
1142 tul = tmp.upperLeft();
1143
1144 for(; dul.x < dlr.x; ++dul.x, ++tul.x)
1145 {
1146 typename DestIterator::column_iterator dc = dul.columnIterator();
1147 typename TmpIterator::column_iterator tc = tul.columnIterator();
1148 resamplingConvolveLine(tc, tc+hold, tmp.accessor(), dc, dc+hnew, dest,
1149 kernels, mapCoordinate);
1150 }
1151 }
1152
1153 template <class SrcIterator, class SrcAccessor,
1154 class DestIterator, class DestAccessor>
1155 inline
pyramidExpandBurtFilter(triple<SrcIterator,SrcIterator,SrcAccessor> src,triple<DestIterator,DestIterator,DestAccessor> dest,double centerValue=0.4)1156 void pyramidExpandBurtFilter(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1157 triple<DestIterator, DestIterator, DestAccessor> dest,
1158 double centerValue = 0.4)
1159 {
1160 pyramidExpandBurtFilter(src.first, src.second, src.third,
1161 dest.first, dest.second, dest.third, centerValue);
1162 }
1163
1164
1165 template <class Image, class Alloc>
1166 inline
pyramidExpandBurtFilter(ImagePyramid<Image,Alloc> & pyramid,int fromLevel,int toLevel,double centerValue=0.4)1167 void pyramidExpandBurtFilter(ImagePyramid<Image, Alloc> & pyramid, int fromLevel, int toLevel,
1168 double centerValue = 0.4)
1169 {
1170 vigra_precondition(fromLevel > toLevel,
1171 "pyramidExpandBurtFilter(): fromLevel must be larger than toLevel.");
1172 vigra_precondition(pyramid.lowestLevel() <= toLevel && fromLevel <= pyramid.highestLevel(),
1173 "pyramidExpandBurtFilter(): fromLevel and toLevel must be between the lowest and highest pyramid levels (inclusive).");
1174
1175 for(int i=fromLevel-1; i >= toLevel; --i)
1176 pyramidExpandBurtFilter(srcImageRange(pyramid[i+1]), destImageRange(pyramid[i]), centerValue);
1177 }
1178
1179 /** \brief Create a Laplacian pyramid.
1180
1181 This function implements the reduction across several resolution levels of
1182 a Laplacian pyramid as described in
1183
1184 P. Burt, E. Adelson: <i>"The Laplacian Pyramid as a Compact Image Code"</i>, IEEE Trans. Communications, 9(4):532-540, 1983
1185
1186 It first creates a Gaussian pyramid using \ref pyramidReduceBurtFilter(), then
1187 upsamples each level once using \ref pyramidExpandBurtFilter(), and finally
1188 stores the difference between the upsampled and original versions of
1189 each level (i.e. the Laplacian of Gaussian is approximated by a difference
1190 of Gaussian).
1191
1192 <b>\#include</b> \<vigra/resampling_convolution.hxx\><br>
1193 Namespace: vigra
1194 */
1195 template <class Image, class Alloc>
1196 inline void
pyramidReduceBurtLaplacian(ImagePyramid<Image,Alloc> & pyramid,int fromLevel,int toLevel,double centerValue=0.4)1197 pyramidReduceBurtLaplacian(ImagePyramid<Image, Alloc> & pyramid,
1198 int fromLevel, int toLevel,
1199 double centerValue = 0.4)
1200 {
1201 using namespace functor;
1202
1203 pyramidReduceBurtFilter(pyramid, fromLevel, toLevel, centerValue);
1204 for(int i=fromLevel; i < toLevel; ++i)
1205 {
1206 typename ImagePyramid<Image, Alloc>::value_type tmpImage(pyramid[i].size());
1207 pyramidExpandBurtFilter(srcImageRange(pyramid[i+1]), destImageRange(tmpImage), centerValue);
1208 combineTwoImages(srcImageRange(tmpImage), srcImage(pyramid[i]), destImage(pyramid[i]),
1209 Arg1() - Arg2());
1210 }
1211 }
1212
1213 /** \brief Reconstruct a Laplacian pyramid.
1214
1215 This function implements the reconstruction of a Gaussian pyramid
1216 across several resolution levels of a Laplacian pyramid as described in
1217
1218 P. Burt, E. Adelson: <i>"The Laplacian Pyramid as a Compact Image Code"</i>, IEEE Trans. Communications, 9(4):532-540, 1983
1219
1220 At each level starting from <tt>fromLevel</tt>, this function calls
1221 \ref pyramidExpandBurtFilter() to interpolate the image to the next highest
1222 resolution, and then adds the interpolated image to the image stored at the
1223 next level.
1224
1225 <b>\#include</b> \<vigra/resampling_convolution.hxx\><br>
1226 Namespace: vigra
1227 */
1228 template <class Image, class Alloc>
1229 inline void
pyramidExpandBurtLaplacian(ImagePyramid<Image,Alloc> & pyramid,int fromLevel,int toLevel,double centerValue=0.4)1230 pyramidExpandBurtLaplacian(ImagePyramid<Image, Alloc> & pyramid,
1231 int fromLevel, int toLevel,
1232 double centerValue = 0.4)
1233 {
1234 using namespace functor;
1235
1236 vigra_precondition(fromLevel > toLevel,
1237 "pyramidExpandBurtLaplacian(): fromLevel must be larger than toLevel.");
1238 vigra_precondition(pyramid.lowestLevel() <= toLevel && fromLevel <= pyramid.highestLevel(),
1239 "pyramidExpandBurtLaplacian(): fromLevel and toLevel must be between the lowest and highest pyramid levels (inclusive).");
1240
1241 for(int i=fromLevel-1; i >= toLevel; --i)
1242 {
1243 typename ImagePyramid<Image, Alloc>::value_type tmpImage(pyramid[i].size());
1244 pyramidExpandBurtFilter(srcImageRange(pyramid[i+1]), destImageRange(tmpImage), centerValue);
1245 combineTwoImages(srcImageRange(tmpImage), srcImage(pyramid[i]), destImage(pyramid[i]),
1246 Arg1() - Arg2());
1247 }
1248 }
1249
1250 //@}
1251
1252 } // namespace vigra
1253
1254
1255 #endif /* VIGRA_RESAMPLING_CONVOLUTION_HXX */
1256