Hugin trunk 0.1
Loading...
Searching...
No Matches
Interpolators.h
Go to the documentation of this file.
1// -*- c-basic-offset: 4 -*-
27#ifndef VIGRA_EXT_INTERPOLATORS_H
28#define VIGRA_EXT_INTERPOLATORS_H
29
30#include <iostream>
31#include <iomanip>
32
33#include <math.h>
35#include <algorithm>
36
37#include <vigra/accessor.hxx>
38#include <vigra/diff2d.hxx>
39
40using std::endl;
41
42namespace vigra_ext {
43
44// Some locally needed math functions
45
46static double sinc ( double x );
47static double cubic01 ( double x );
48static double cubic12 ( double x );
49
50static double sinc( double x )
51{
52 x *= M_PI;
53 if(x != 0.0)
54 return(sin(x) / x);
55 return(1.0);
56}
57
58
59// Cubic polynomial with parameter A
60// A = -1: sharpen; A = - 0.5 homogeneous
61// make sure x >= 0
62static const double A(-0.75);
63
64// 0 <= x < 1
65static double cubic01( double x )
66{
67 return (( A + 2.0 )*x - ( A + 3.0 ))*x*x +1.0;
68}
69// 1 <= x < 2
70
71static double cubic12( double x )
72{
73 return (( A * x - 5.0 * A ) * x + 8.0 * A ) * x - 4.0 * A;
74
75}
76
93
94
102{
103 // size of neighbourhood
104 static const int size = 2;
105
106 void calc_coeff(double x, double * w) const
107 {
108 w[1] = (x >= 0.5) ? 1 : 0;
109 w[0] = (x < 0.5) ? 1 : 0;
110 }
111
112 void emitGLSL(std::ostringstream& oss) const {
113 oss << " return (i == 0.0) ? float(f < 0.5) : float(f >= 0.5);" << endl;
114 }
115};
116
117
120{
121 // size of neighbourhood
122 static const int size = 2;
123
124 void calc_coeff(double x, double * w) const
125 {
126 w[1] = x;
127 w[0] = 1.0-x;
128 }
129
130 void emitGLSL(std::ostringstream& oss) const {
131 oss << " return abs(i + f - 1.0);" << endl;
132 }
133};
134
137{
138 // size of neighbourhood
139 static const int size = 4;
140
142 void calc_coeff(double x, double * w) const
143 {
144 w[3] = cubic12( 2.0 - x );
145 w[2] = cubic01( 1.0 - x );
146 w[1] = cubic01( x );
147 w[0] = cubic12( x + 1.0 );
148 }
149
150 void emitGLSL(std::ostringstream& oss) const {
151 oss << " float A = " << A << ";" << endl
152 << " float c = abs(i - 1.0);" << endl
153 << " float m = (i > 1.0) ? -1.0 : 1.0;" << endl
154 << " float p = c + m * f;" << endl
155 << " if (i == 1.0 || i == 2.0) {" << endl
156 << " return (( A + 2.0 )*p - ( A + 3.0 ))*p*p + 1.0;" << endl
157 << " } else {" << endl
158 << " return (( A * p - 5.0 * A ) * p + 8.0 * A ) * p - 4.0 * A;" << endl
159 << " }" << endl;
160 }
161};
162
165{
166 // size of neighbourhood
167 static const int size = 4;
168
170 void calc_coeff(double x, double * w) const
171 {
172 w[3] = ( ( 1.0/3.0 * x - 1.0/5.0 ) * x - 2.0/15.0 ) * x;
173 w[2] = ( ( 6.0/5.0 - x ) * x + 4.0/5.0 ) * x;
174 w[1] = ( ( x - 9.0/5.0 ) * x - 1.0/5.0 ) * x + 1.0;
175 w[0] = ( ( -1.0/3.0 * x + 4.0/5.0 ) * x - 7.0/15.0 ) * x;
176 }
177
178 void emitGLSL(std::ostringstream& oss) const {
179 oss << " return (i > 1.0) ? (i == 3.0) ? (( ( 1.0/3.0 * f - 1.0/5.0 ) * f - 2.0/15.0 ) * f)" << endl
180 << " : (( ( 6.0/5.0 - f ) * f + 4.0/5.0 ) * f)" << endl
181 << " : (i == 1.0) ? (( ( f - 9.0/5.0 ) * f - 1.0/5.0 ) * f + 1.0)" << endl
182 << " : (( ( -1.0/3.0 * f + 4.0/5.0 ) * f - 7.0/15.0 ) * f);" << endl;
183 }
184};
185
188{
189 // size of neighbourhood
190 static const int size = 6;
191
198 void calc_coeff(double x, double* w) const
199 {
200 w[5] = ( ( - 1.0/11.0 * x + 12.0/ 209.0 ) * x + 7.0/ 209.0 ) * x;
201 w[4] = ( ( 6.0/11.0 * x - 72.0/ 209.0 ) * x - 42.0/ 209.0 ) * x;
202 w[3] = ( ( - 13.0/11.0 * x + 288.0/ 209.0 ) * x + 168.0/ 209.0 ) * x;
203 w[2] = ( ( 13.0/11.0 * x - 453.0/ 209.0 ) * x - 3.0/ 209.0 ) * x + 1.0;
204 w[1] = ( ( - 6.0/11.0 * x + 270.0/ 209.0 ) * x - 156.0/ 209.0 ) * x;
205 w[0] = ( ( 1.0/11.0 * x - 45.0/ 209.0 ) * x + 26.0/ 209.0 ) * x;
206 }
207
208 void emitGLSL(std::ostringstream& oss) const {
209 oss << " return (i > 3.0) ? (i == 5.0) ? (( ( - 1.0/11.0 * f + 12.0/ 209.0 ) * f + 7.0/ 209.0 ) * f)" << endl
210 << " : (( ( 6.0/11.0 * f - 72.0/ 209.0 ) * f - 42.0/ 209.0 ) * f)" << endl
211 << " : (i > 1.0) ? (i == 3.0) ? (( ( - 13.0/11.0 * f + 288.0/ 209.0 ) * f + 168.0/ 209.0 ) * f)" << endl
212 << " : (( ( 13.0/11.0 * f - 453.0/ 209.0 ) * f - 3.0/ 209.0 ) * f + 1.0)" << endl
213 << " : (i == 1.0) ? (( ( - 6.0/11.0 * f + 270.0/ 209.0 ) * f - 156.0/ 209.0 ) * f)" << endl
214 << " : (( ( 1.0/11.0 * f - 45.0/ 209.0 ) * f + 26.0/ 209.0 ) * f);" << endl;
215 }
216};
217
218
221{
222 // size of neighbourhood
223 static const int size = 8;
224
226 void calc_coeff(double x, double * w) const
227 {
228 w[7] = (( 1.0/41.0 * x - 45.0/2911.0) * x - 26.0/2911.0) * x;
229 w[6] = ((- 6.0/41.0 * x + 270.0/2911.0) * x + 156.0/2911.0) * x;
230 w[5] = (( 24.0/41.0 * x - 1080.0/2911.0) * x - 624.0/2911.0) * x;
231 w[4] = ((-49.0/41.0 * x + 4050.0/2911.0) * x + 2340.0/2911.0) * x;
232 w[3] = (( 49.0/41.0 * x - 6387.0/2911.0) * x - 3.0/2911.0) * x + 1.0;
233 w[2] = ((-24.0/41.0 * x + 4032.0/2911.0) * x - 2328.0/2911.0) * x;
234 w[1] = (( 6.0/41.0 * x - 1008.0/2911.0) * x + 582.0/2911.0) * x;
235 w[0] = ((- 1.0/41.0 * x + 168.0/2911.0) * x - 97.0/2911.0) * x;
236 }
237
238 void emitGLSL(std::ostringstream& oss) const {
239 oss << " return (i > 3.0) ? (i > 5.0) ? (i == 7.0) ? ((( 1.0/41.0 * f - 45.0/2911.0) * f - 26.0/2911.0) * f)" << endl
240 << " : (((- 6.0/41.0 * f + 270.0/2911.0) * f + 156.0/2911.0) * f)" << endl
241 << " : (i == 5.0) ? ((( 24.0/41.0 * f - 1080.0/2911.0) * f - 624.0/2911.0) * f)" << endl
242 << " : (((-49.0/41.0 * f + 4050.0/2911.0) * f + 2340.0/2911.0) * f)" << endl
243 << " : (i > 1.0) ? (i == 3.0) ? ((( 49.0/41.0 * f - 6387.0/2911.0) * f - 3.0/2911.0) * f + 1.0)" << endl
244 << " : (((-24.0/41.0 * f + 4032.0/2911.0) * f - 2328.0/2911.0) * f)" << endl
245 << " : (i == 1.0) ? ((( 6.0/41.0 * f - 1008.0/2911.0) * f + 582.0/2911.0) * f)" << endl
246 << " : (((- 1.0/41.0 * f + 168.0/2911.0) * f - 97.0/2911.0) * f);" << endl;
247 }
248};
249
251template <int size_>
253{
254 // size of neighbourhood
255 static const int size = size_;
256
258 void calc_coeff(double x, double * w) const
259 {
260 int idx;
261 double xadd;
262 for( idx = 0, xadd = size / 2 - 1.0 + x;
263 idx < size / 2;
264 xadd-=1.0)
265 {
266 w[idx++] = sinc( xadd ) * sinc( xadd / ( size / 2 ));
267 }
268 for( xadd = 1.0 - x;
269 idx < size;
270 xadd+=1.0)
271 {
272 w[idx++] = sinc( xadd ) * sinc( xadd / ( size / 2 ));
273 }
274 }
275
276 void emitGLSL(std::ostringstream& oss) const {
277 oss << " float c = (i < " << (size/2.0) << ") ? 1.0 : -1.0;" << endl
278 << " float x = c * (" << (size/2 - 1.0) << " - i + f);" << endl
279 << " vec2 xpi = vec2(x, x / " << (size/2.0) << ") * " << M_PI << ";" << endl
280 << " vec2 xsin = sin(xpi);" << endl
281 << " vec2 result = vec2(1.0, 1.0);" << endl
282 << " if (xpi.x != 0.0) result.x = xsin.x / xpi.x;" << endl
283 << " if (xpi.y != 0.0) result.y = xsin.y / xpi.y;" << endl
284 << " return result.x * result.y;" << endl;
285 }
286};
287
288
293template <typename SrcImageIterator, typename SrcAccessor,
294 typename INTERPOLATOR>
296{
297public:
298 typedef typename SrcAccessor::value_type PixelType;
299 // dummy mask type to be compatible to algorithms expecting a ImageMaskInterpolator object
300 typedef typename vigra::UInt8 MaskType;
301private:
302 typedef typename vigra::NumericTraits<PixelType>::RealPromote RealPixelType;
303
306 int m_w;
307 int m_h;
309
311
312public:
314 ImageInterpolator(vigra::triple<SrcImageIterator, SrcImageIterator,SrcAccessor> const & src,
316 bool warparound)
317 : m_sIter(src.first),
319 m_w(src.second.x - src.first.x),
320 m_h(src.second.y - src.first.y),
323 {
324 }
325
342
344 bool operator()(double x, double y,
345 PixelType & result, MaskType & mask) const
346 {
347 mask = 255;
348 return operator()(x,y, result);
349 }
350
352 bool operator()(double x, double y,
353 PixelType & result) const
354 {
355
356 // skip all further interpolation if we cannot interpolate anything
357 if (x < -INTERPOLATOR::size/2 || x > m_w + INTERPOLATOR::size/2) return false;
358 if (y < -INTERPOLATOR::size/2 || y > m_h + INTERPOLATOR::size/2) return false;
359
360 double t = floor(x);
361 double dx = x - t;
362 int srcx = int(t);
363 t = floor(y);
364 double dy = y - t;
365 int srcy = int(t);
366
367 if ( srcx > INTERPOLATOR::size/2 && srcx < m_w -INTERPOLATOR::size/2 &&
368 srcy > INTERPOLATOR::size/2 && srcy < m_h - INTERPOLATOR::size/2)
369 {
370 return interpolateNoMaskInside(srcx, srcy, dx, dy, result);
371 }
372
373 double wx[INTERPOLATOR::size];
374 double wy[INTERPOLATOR::size];
375
376 // calculate x interpolation coefficients
377 m_inter.calc_coeff(dx, wx);
378 m_inter.calc_coeff(dy, wy);
379
380 RealPixelType p(vigra::NumericTraits<RealPixelType>::zero());
381 double weightsum = 0.0;
382 for (int ky = 0; ky < INTERPOLATOR::size; ky++) {
383 int bounded_ky = srcy + 1 + ky - INTERPOLATOR::size/2;
384
385 // Boundary condition: do not replicate top and bottom
387 continue;
388 }
389
390 for (int kx = 0; kx < INTERPOLATOR::size; kx++) {
391 int bounded_kx = srcx + 1 + kx - INTERPOLATOR::size/2;
392
393 if (m_warparound) {
394 // Boundary condition: wrap around the image.
395 if (bounded_kx < 0)
396 bounded_kx += m_w;
397 if (bounded_kx >= m_w)
398 bounded_kx -= m_w;
399 } else {
400 // Boundary condition: replicate first and last column.
401 // if (srcx + kx < 0) bounded_kx -= (srcx + kx);
402 // if (srcx + kx >= src_w) bounded_kx -= (srcx + kx - (src_w - 1));
403 // Boundary condition: do not replicate left and right
404 if (bounded_kx < 0)
405 continue;
406 if (bounded_kx >= m_w)
407 continue;
408 }
409
410 // check mask
411 double f = wx[kx]*wy[ky];
412 p += f * m_sAcc(m_sIter, vigra::Diff2D(bounded_kx, bounded_ky));
413 weightsum += f;
414 }
415 }
416
417 // force a certain weight
418 if (weightsum <= 0.2) return false;
419 // Adjust filter for any ignored transparent pixels.
420 if (weightsum != 1.0) p /= weightsum;
421
422 result = vigra::detail::RequiresExplicitCast<PixelType>::cast(p);
423 return true;
424 }
425
426
428 bool interpolateNoMaskInside(int srcx, int srcy, double dx, double dy,
429 PixelType & result) const
430 {
431 double w[INTERPOLATOR::size];
432 RealPixelType resX[INTERPOLATOR::size];
433
434 // calculate x interpolation coefficients
435 m_inter.calc_coeff(dx, w);
436
438
439 // first pass of separable filter, x pass
440 vigra::Diff2D offset(srcx - INTERPOLATOR::size/2 + 1,
441 srcy - INTERPOLATOR::size/2 + 1);
442 SrcImageIterator ys(m_sIter + offset);
443 for (int ky = 0; ky < INTERPOLATOR::size; ky++, ++(ys.y)) {
444 p = vigra::NumericTraits<RealPixelType>::zero();
445 typename SrcImageIterator::row_iterator xs(ys.rowIterator());
446 //SrcImageIterator xs(ys);
447 for (int kx = 0; kx < INTERPOLATOR::size; kx++, ++xs) {
448 p += w[kx] * m_sAcc(xs);
449 }
450 resX[ky] = p;
451 }
452
453 // y pass.
454 m_inter.calc_coeff(dy, w);
455 p = vigra::NumericTraits<RealPixelType>::zero();
456 for (int ky = 0; ky < INTERPOLATOR::size; ky++) {
457 p += w[ky] * resX[ky];
458 }
459
460 result = vigra::detail::RequiresExplicitCast<PixelType>::cast(p);
461 return true;
462 }
463
464 void emitGLSL(std::ostringstream& oss) const {
465 m_inter.emitGLSL(oss);
466 }
467
468};
469
470
476template <typename SrcImageIterator, typename SrcAccessor,
477 typename MaskIterator, typename MaskAccessor,
478 typename INTERPOLATOR>
480{
481public:
482 typedef typename SrcAccessor::value_type PixelType;
483 typedef typename MaskAccessor::value_type MaskType;
484private:
485 typedef typename vigra::NumericTraits<PixelType>::RealPromote RealPixelType;
486
490 MaskAccessor m_mAcc;
491 int m_w;
492 int m_h;
494
496
497public:
498
500 ImageMaskInterpolator(vigra::triple<SrcImageIterator, SrcImageIterator,SrcAccessor> const & src,
501 std::pair<MaskIterator, MaskAccessor> mask,
503 bool warparound)
504 : m_sIter(src.first),
506 m_mIter(mask.first),
507 m_mAcc(mask.second),
508 m_w(src.second.x - src.first.x),
509 m_h(src.second.y - src.first.y),
512 {
513 }
514
535#if 0
555// bool operator()(float x, float y,
556 // this is slower than the full version, thanks to the normalized interpolation (with masks).
557 bool interpolateSeperable(float x, float y,
558 PixelType & result) const
559 {
560
561 // skip all further interpolation if we cannot interpolate anything
562 if (x < -INTERPOLATOR::size/2 || x > m_w + INTERPOLATOR::size/2) return false;
563 if (y < -INTERPOLATOR::size/2 || y > m_h + INTERPOLATOR::size/2) return false;
564
565 double t = floor(x);
566 double dx = x - t;
567 int srcx = int(t);
568 t = floor(y);
569 double dy = y - t;
570 int srcy = int(t);
571
572
573 double w[INTERPOLATOR::size];
574
575 double weightsX[INTERPOLATOR::size];
576 PixelType resX[INTERPOLATOR::size];
577
578 // calculate x interpolation coefficients
579 m_inter.calc_coeff(dx, w);
580
581 // first pass of separable filter
582
583 for (int ky = 0; ky < INTERPOLATOR::size; ky++) {
584 int bounded_ky = srcy + 1 + ky - INTERPOLATOR::size/2;
585
586 // Boundary condition: replicate top and bottom rows.
587 // if (srcy + ky < 0) bounded_ky -= (srcy + ky);
588 // if (srcy + ky >= src_h) bounded_ky -= (srcy + ky - (src_h - 1));
589
590 // Boundary condition: do not replicate top and bottom
592 weightsX[ky] = 0;
593 resX[ky] = 0;
594 continue;
595 }
596
597 RealPixelType p(vigra::NumericTraits<RealPixelType>::zero());
598 double weightsum = 0.0;
599
600 for (int kx = 0; kx < INTERPOLATOR::size; kx++) {
601 int bounded_kx = srcx + 1 + kx - INTERPOLATOR::size/2;
602
603 if (m_warparound) {
604 // Boundary condition: wrap around the image.
605 if (bounded_kx < 0)
606 bounded_kx += m_w;
607 if (bounded_kx >= m_w)
608 bounded_kx -= m_w;
609 } else {
610 // Boundary condition: replicate first and last column.
611 // if (srcx + kx < 0) bounded_kx -= (srcx + kx);
612 // if (srcx + kx >= src_w) bounded_kx -= (srcx + kx - (src_w - 1));
613 // Boundary condition: do not replicate left and right
614 if (bounded_kx < 0)
615 continue;
616 if (bounded_kx >= m_w)
617 continue;
618 }
620 // check mask
621 p += w[kx] * m_sIter(bounded_kx, bounded_ky);
622 weightsum += w[kx];
623 }
624 }
626 resX[ky] = p;
627 }
628
629 // y pass.
630 m_inter.calc_coeff(dy, w);
631 RealPixelType p(vigra::NumericTraits<RealPixelType>::zero());
632 double weightsum = 0.0;
633 for (int ky = 0; ky < INTERPOLATOR::size; ky++) {
634 weightsum += weightsX[ky] * w[ky];
635 p += w[ky] * resX[ky];
636 }
637
638 if (weightsum == 0.0) return false;
639 // Adjust filter for any ignored transparent pixels.
640 if (weightsum != 1.0) p /= weightsum;
641
642 result = vigra::detail::RequiresExplicitCast<PixelType>::cast(p);
643 return true;
644 }
645#endif
646
666 bool operator()(double x, double y,
667 PixelType & result, MaskType & mask) const
668 {
669
670 // skip all further interpolation if we cannot interpolate anything
671 if (x < -INTERPOLATOR::size/2 || x > m_w + INTERPOLATOR::size/2) return false;
672 if (y < -INTERPOLATOR::size/2 || y > m_h + INTERPOLATOR::size/2) return false;
673
674 double t = floor(x);
675 double dx = x - t;
676 int srcx = int(t);
677 t = floor(y);
678 double dy = y - t;
679 int srcy = int(t);
680
681 if ( srcx > INTERPOLATOR::size/2 && srcx < m_w -INTERPOLATOR::size/2 &&
682 srcy > INTERPOLATOR::size/2 && srcy < m_h - INTERPOLATOR::size/2)
683 {
684 return interpolateInside(srcx, srcy, dx, dy, result, mask);
685 }
686
687 double wx[INTERPOLATOR::size];
688 double wy[INTERPOLATOR::size];
689
690 // calculate x interpolation coefficients
691 m_inter.calc_coeff(dx, wx);
692 m_inter.calc_coeff(dy, wy);
693
694 // first pass of separable filter
695
696 RealPixelType p(vigra::NumericTraits<RealPixelType>::zero());
697 double m = 0;
698 double weightsum = 0.0;
699 for (int ky = 0; ky < INTERPOLATOR::size; ky++) {
700 int bounded_ky = srcy + 1 + ky - INTERPOLATOR::size/2;
701
702 // Boundary condition: do not replicate top and bottom
704 continue;
705 }
706
707 for (int kx = 0; kx < INTERPOLATOR::size; kx++) {
708 int bounded_kx = srcx + 1 + kx - INTERPOLATOR::size/2;
709
710 if (m_warparound) {
711 // Boundary condition: wrap around the image.
712 if (bounded_kx < 0)
713 bounded_kx += m_w;
714 if (bounded_kx >= m_w)
715 bounded_kx -= m_w;
716 } else {
717 // Boundary condition: replicate first and last column.
718 // if (srcx + kx < 0) bounded_kx -= (srcx + kx);
719 // if (srcx + kx >= src_w) bounded_kx -= (srcx + kx - (src_w - 1));
720 // Boundary condition: do not replicate left and right
721 if (bounded_kx < 0)
722 continue;
723 if (bounded_kx >= m_w)
724 continue;
725 }
726
728 if (cmask) {
729 // check mask
730 double f = wx[kx]*wy[ky];
731 // TODO: check if this is good, influences the HDR stitching masks
732 m += f * cmask;
733 p += f * m_sAcc(m_sIter, vigra::Diff2D(bounded_kx, bounded_ky));
734 weightsum += f;
735 }
736 }
737 }
738
739 // force a certain weight
740 if (weightsum <= 0.2) return false;
741 // Adjust filter for any ignored transparent pixels.
742 if (weightsum != 1.0) {
743 p /= weightsum;
744 m /= weightsum;
745 }
746
747 mask = vigra::detail::RequiresExplicitCast<MaskType>::cast(m);
748 result = vigra::detail::RequiresExplicitCast<PixelType>::cast(p);
749 return true;
750 }
751
752
754 bool interpolateInside(int srcx, int srcy, double dx, double dy,
755 PixelType & result, MaskType & mask) const
756 {
757
758 double wx[INTERPOLATOR::size];
759 double wy[INTERPOLATOR::size];
760
761 // calculate x interpolation coefficients
762 m_inter.calc_coeff(dx, wx);
763 m_inter.calc_coeff(dy, wy);
764
765 RealPixelType p(vigra::NumericTraits<RealPixelType>::zero());
766 double weightsum = 0.0;
767 double m = 0.0;
768 vigra::Diff2D offset(srcx - INTERPOLATOR::size/2 + 1,
769 srcy - INTERPOLATOR::size/2 + 1);
770 SrcImageIterator ys(m_sIter + offset);
771 MaskIterator yms(m_mIter + offset);
772 for (int ky = 0; ky < INTERPOLATOR::size; ky++, ++(ys.y), ++(yms.y)) {
773// int bounded_ky = srcy + 1 + ky - INTERPOLATOR::size/2;
774 typename SrcImageIterator::row_iterator xs(ys.rowIterator());
775 typename MaskIterator::row_iterator xms(yms.rowIterator());
776 for (int kx = 0; kx < INTERPOLATOR::size; kx++, ++xs, ++xms) {
777// int bounded_kx = srcx + 1 + kx - INTERPOLATOR::size/2;
778
779 MaskType cmask = *xms;
780 if (cmask) {
781 // check mask
782 double f = wx[kx]*wy[ky];
783 // TODO: check if this is good, influences the HDR stitching masks
784 m += f * cmask;
785 p += f * m_sAcc(xs);
786 weightsum += f;
787 }
788 }
789 }
790
791 // force a certain weight
792 if (weightsum <= 0.2) return false;
793 // Adjust filter for any ignored transparent pixels.
794 if (weightsum != 1.0) {
795 p /= weightsum;
796 m /= weightsum;
797 }
798
799 result = vigra::detail::RequiresExplicitCast<PixelType>::cast(p);
800 mask = vigra::detail::RequiresExplicitCast<MaskType>::cast(m);
801 return true;
802 }
803
804};
805
806/********************************************************/
807/* */
808/* InterpolatingAccessor */
809/* */
810/********************************************************/
811
861template <class ACCESSOR, class VALUETYPE, class INTERPOLATOR>
863{
864public:
868
874
891 template <class ITERATOR>
892 value_type operator()(ITERATOR const & i, float x, float y) const
893 {
894 int ix = int(x);
895 int iy = int(y);
896 float dx = x - ix;
897 float dy = y - iy;
898 double wx[INTERPOLATOR::size];
899 double wy[INTERPOLATOR::size];
900
901 // promote value_type for multiplication
902 typename vigra::NumericTraits<value_type>::RealPromote
903 ret (vigra::NumericTraits<value_type>::zero());
904
905 // calculate interpolation coefficients
906 inter_x.calc_coeff(dx, wx);
907 inter_y.calc_coeff(dy, wy);
908
909 ITERATOR ys(i + vigra::Diff2D(ix - inter_x.size/2 + 1,
910 iy - inter_y.size/2 + 1));
911 for(int y = 0; y < inter_y.size; ++y, ++ys.y) {
912 ITERATOR xs(ys);
913 for(int x = 0; x < inter_x.size; x++, ++xs.x) {
914 ret += wx[x] * wy[y] * a_(xs);
915 }
916 }
917 return vigra::detail::RequiresExplicitCast<value_type>::cast(ret);
918 }
919
941 template <class ITERATOR, class ALPHAITERATOR, class ALPHAACCESSOR>
942 bool operator()(ITERATOR const & i, std::pair<ALPHAITERATOR, ALPHAACCESSOR> const & alpha,
943 float x, float y, value_type & result) const
944 {
945 int ix = int(x);
946 int iy = int(y);
947 float dx = x - ix;
948 float dy = y - iy;
949 double wx[INTERPOLATOR::size];
950 double wy[INTERPOLATOR::size];
951
952 // promote value_type for multiplication
953 typename vigra::NumericTraits<value_type>::RealPromote
954 ret (vigra::NumericTraits<value_type>::zero());
955
956 // calculate interpolation coefficients
957 inter_x.calc_coeff(dx, wx);
958 inter_y.calc_coeff(dy, wy);
959
960 ITERATOR ys(i + vigra::Diff2D(ix - inter_x.size/2 + 1,
961 iy - inter_y.size/2 + 1));
962 ALPHAITERATOR ays(alpha.first + vigra::Diff2D(ix - inter_x.size/2 + 1,
963 iy - inter_y.size/2 + 1));
964 for(int y = 0; y < inter_y.size; ++y, ++ys.y, ++ays.y) {
965 ITERATOR xs(ys);
967 for(int x = 0; x < inter_x.size; x++, ++xs.x, ++axs.x) {
968 if (alpha.second(axs) <= 0 ) {
969 return false;
970 }
971 ret += wx[x] * wy[y] * a_(xs);
972 }
973 }
974 result = vigra::detail::RequiresExplicitCast<value_type>::cast(ret);
975 return true;
976 }
977
978
979private:
982};
983
984} // namespace
985
986#endif // VIGRA_EXT_INTERPOLATORS_H
#define M_PI
"wrapper" for efficient interpolation access to an image
bool operator()(double x, double y, PixelType &result, MaskType &mask) const
Interpolate without mask, but return dummy alpha value nevertheless.
ImageInterpolator(SrcImageIterator src_upperleft, SrcImageIterator src_lowerright, SrcAccessor sa, INTERPOLATOR &inter, bool warparound)
Construct interpolator for specific image.
bool operator()(double x, double y, PixelType &result) const
Interpolate without mask.
ImageInterpolator(vigra::triple< SrcImageIterator, SrcImageIterator, SrcAccessor > const &src, INTERPOLATOR &inter, bool warparound)
Construct interpolator for an given image.
vigra::NumericTraits< PixelType >::RealPromote RealPixelType
SrcAccessor::value_type PixelType
bool interpolateNoMaskInside(int srcx, int srcy, double dx, double dy, PixelType &result) const
Interpolate without boundary check and mask.
void emitGLSL(std::ostringstream &oss) const
"wrapper" for efficient interpolation access to an image
bool operator()(double x, double y, PixelType &result, MaskType &mask) const
Interpolate the data item at a non-integer position x, y.
MaskAccessor::value_type MaskType
bool interpolateInside(int srcx, int srcy, double dx, double dy, PixelType &result, MaskType &mask) const
Interpolate without boundary check.
ImageMaskInterpolator(SrcImageIterator src_upperleft, SrcImageIterator src_lowerright, SrcAccessor sa, MaskIterator mask_upperleft, MaskAccessor ma, INTERPOLATOR &inter, bool warparound)
Construct interpolator for specific image.
ImageMaskInterpolator(vigra::triple< SrcImageIterator, SrcImageIterator, SrcAccessor > const &src, std::pair< MaskIterator, MaskAccessor > mask, INTERPOLATOR &inter, bool warparound)
Construct interpolator for an given image.
vigra::NumericTraits< PixelType >::RealPromote RealPixelType
SrcAccessor::value_type PixelType
interpolation at non-integer positions.
VALUETYPE value_type
the iterators' pixel type
InterpolatingAccessor(ACCESSOR a, INTERPOLATOR inter)
init from given accessor
value_type operator()(ITERATOR const &i, float x, float y) const
Interpolate the data item at a non-integer position x, y.
bool operator()(ITERATOR const &i, std::pair< ALPHAITERATOR, ALPHAACCESSOR > const &alpha, float x, float y, value_type &result) const
Interpolate the data item at a non-integer position x, y.
misc math function & classes used by other parts of the program
static double sinc(double x)
static const double A(-0.75)
static double cubic01(double x)
Interpolator
enum with all interpolation methods
@ INTERP_NEAREST_NEIGHBOUR
static double cubic12(double x)
simple bilinear interpolation
static const int size
void calc_coeff(double x, double *w) const
void emitGLSL(std::ostringstream &oss) const
cubic interpolation
static const int size
void emitGLSL(std::ostringstream &oss) const
void calc_coeff(double x, double *w) const
initialize weights for given x
several classes to calculate interpolator weights,
void emitGLSL(std::ostringstream &oss) const
void calc_coeff(double x, double *w) const
sinc interpolation, with variable width
static const int size
void calc_coeff(double x, double *w) const
initialize weights for given offset x
void emitGLSL(std::ostringstream &oss) const
spline16 interpolation
void emitGLSL(std::ostringstream &oss) const
void calc_coeff(double x, double *w) const
initialize weights for given x
spline36 interpolation
void emitGLSL(std::ostringstream &oss) const
void calc_coeff(double x, double *w) const
calculate weights for given offset x.
spline64 interpolation
void emitGLSL(std::ostringstream &oss) const
void calc_coeff(double x, double *w) const
initialize weights for given offset x
std::vector< deghosting::BImagePtr > threshold(const std::vector< deghosting::FImagePtr > &inputImages, const double threshold, const uint16_t flags)
Threshold function used for creating alpha masks for images.
Definition threshold.h:41