Main Page | Modules | Namespace List | Class Hierarchy | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages | Examples | Search for 
hugin_base / vigra_ext

Interpolators.h

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 -*-
00027 #ifndef VIGRA_EXT_INTERPOLATORS_H
00028 #define VIGRA_EXT_INTERPOLATORS_H
00029 
00030 #include <iostream>
00031 #include <iomanip>
00032 
00033 #include <math.h>
00034 #include <hugin_math/hugin_math.h>
00035 #include <algorithm>
00036 
00037 #include <vigra/accessor.hxx>
00038 #include <vigra/diff2d.hxx>
00039 
00040 using std::endl;
00041 
00042 namespace vigra_ext {
00043 
00044 // Some locally needed math functions
00045 
00046 static double   sinc            ( double x );
00047 static double   cubic01         ( double x );
00048 static double   cubic12         ( double x );
00049 
00050 static double sinc( double x )
00051 {
00052     x *= M_PI;
00053     if(x != 0.0)
00054         return(sin(x) / x);
00055     return(1.0);
00056 }
00057 
00058 
00059 // Cubic polynomial with parameter A
00060 // A = -1: sharpen; A = - 0.5 homogeneous
00061 // make sure x >= 0
00062 static const double     A(-0.75);
00063 
00064 // 0 <= x < 1
00065 static double cubic01( double x )
00066 {
00067     return      (( A + 2.0 )*x - ( A + 3.0 ))*x*x +1.0;
00068 }
00069 // 1 <= x < 2
00070 
00071 static double cubic12( double x )
00072 {
00073     return      (( A * x - 5.0 * A ) * x + 8.0 * A ) * x - 4.0 * A;
00074 
00075 }
00076 
00078 enum Interpolator {
00079     INTERP_CUBIC = 0,
00080     INTERP_SPLINE_16,
00081     INTERP_SPLINE_36,
00082     INTERP_SINC_256,
00083     INTERP_SPLINE_64,
00084     INTERP_BILINEAR,
00085     INTERP_NEAREST_NEIGHBOUR,
00086     INTERP_SINC_1024
00087 };
00088 
00089 
00096 struct interp_nearest
00097 {
00098     // size of neighbourhood
00099     static const int size = 2;
00100 
00101     void calc_coeff(double x, double * w) const
00102         {
00103             w[1] = (x >= 0.5) ? 1 : 0;
00104             w[0] = (x < 0.5) ? 1 : 0;
00105         }
00106 
00107     void emitGLSL(std::ostringstream& oss) const {
00108         oss << "    return (i == 0.0) ? float(f < 0.5) : float(f >= 0.5);" << endl;
00109     }
00110 };
00111 
00112 
00114 struct interp_bilin
00115 {
00116     // size of neighbourhood
00117     static const int size = 2;
00118 
00119     void calc_coeff(double x, double * w) const
00120         {
00121             w[1] = x;
00122             w[0] = 1.0-x;
00123         }
00124 
00125     void emitGLSL(std::ostringstream& oss) const {
00126         oss << "    return abs(i + f - 1.0);" << endl;
00127     }
00128 };
00129 
00131 struct interp_cubic
00132 {
00133     // size of neighbourhood
00134     static const int size = 4;
00135 
00137     void calc_coeff(double x, double * w) const
00138         {
00139             w[3] = cubic12( 2.0 - x );
00140             w[2] = cubic01( 1.0 - x );
00141             w[1] = cubic01( x );
00142             w[0] = cubic12( x + 1.0 );
00143         }
00144 
00145     void emitGLSL(std::ostringstream& oss) const {
00146         oss << "    float A = " << A << ";" << endl
00147             << "    float c = abs(i - 1.0);" << endl
00148             << "    float m = (i > 1.0) ? -1.0 : 1.0;" << endl
00149             << "    float p = c + m * f;" << endl
00150             << "    if (i == 1.0 || i == 2.0) {" << endl
00151             << "        return (( A + 2.0 )*p - ( A + 3.0 ))*p*p + 1.0;" << endl
00152             << "    } else {" << endl
00153             << "        return (( A * p - 5.0 * A ) * p + 8.0 * A ) * p - 4.0 * A;" << endl
00154             << "    }" << endl;
00155     }
00156 };
00157 
00159 struct interp_spline16
00160 {
00161     // size of neighbourhood
00162     static const int size = 4;
00163 
00165     void calc_coeff(double x, double * w) const
00166         {
00167             w[3] = ( ( 1.0/3.0  * x - 1.0/5.0 ) * x -   2.0/15.0 ) * x;
00168             w[2] = ( ( 6.0/5.0 - x     ) * x +   4.0/5.0 ) * x;
00169             w[1] = ( ( x - 9.0/5.0 ) * x -   1.0/5.0     ) * x + 1.0;
00170             w[0] = ( ( -1.0/3.0 * x + 4.0/5.0     ) * x -   7.0/15.0 ) * x;
00171         }
00172 
00173     void emitGLSL(std::ostringstream& oss) const {
00174         oss << "    return (i > 1.0) ? (i == 3.0) ? (( ( 1.0/3.0  * f - 1.0/5.0 ) * f -   2.0/15.0 ) * f)" << endl
00175             << "                                  : (( ( 6.0/5.0 - f     ) * f +   4.0/5.0 ) * f)" << endl
00176             << "                     : (i == 1.0) ? (( ( f - 9.0/5.0 ) * f -   1.0/5.0     ) * f + 1.0)" << endl
00177             << "                                  : (( ( -1.0/3.0 * f + 4.0/5.0     ) * f -   7.0/15.0 ) * f);" << endl;
00178     }
00179 };
00180 
00182 struct interp_spline36
00183 {
00184     // size of neighbourhood
00185     static const int size = 6;
00186 
00193     void calc_coeff(double x, double* w) const
00194         {
00195             w[5] = ( ( -  1.0/11.0  * x +  12.0/ 209.0 ) * x +   7.0/ 209.0  ) * x;
00196             w[4] = ( (    6.0/11.0  * x -  72.0/ 209.0 ) * x -  42.0/ 209.0  ) * x;
00197             w[3] = ( ( - 13.0/11.0  * x + 288.0/ 209.0 ) * x + 168.0/ 209.0  ) * x;
00198             w[2] = ( (   13.0/11.0  * x - 453.0/ 209.0 ) * x -   3.0/ 209.0  ) * x + 1.0;
00199             w[1] = ( ( -  6.0/11.0  * x + 270.0/ 209.0 ) * x - 156.0/ 209.0  ) * x;
00200             w[0] = ( (    1.0/11.0  * x -  45.0/ 209.0 ) * x +  26.0/ 209.0  ) * x;
00201         }
00202 
00203     void emitGLSL(std::ostringstream& oss) const {
00204         oss << "    return (i > 3.0) ? (i == 5.0) ? (( ( -  1.0/11.0  * f +  12.0/ 209.0 ) * f +   7.0/ 209.0  ) * f)" << endl
00205             << "                                  : (( (    6.0/11.0  * f -  72.0/ 209.0 ) * f -  42.0/ 209.0  ) * f)" << endl
00206             << "                     : (i > 1.0) ? (i == 3.0) ? (( ( - 13.0/11.0  * f + 288.0/ 209.0 ) * f + 168.0/ 209.0  ) * f)" << endl
00207             << "                                              : (( (   13.0/11.0  * f - 453.0/ 209.0 ) * f -   3.0/ 209.0  ) * f + 1.0)" << endl
00208             << "                                 : (i == 1.0) ? (( ( -  6.0/11.0  * f + 270.0/ 209.0 ) * f - 156.0/ 209.0  ) * f)" << endl
00209             << "                                              : (( (    1.0/11.0  * f -  45.0/ 209.0 ) * f +  26.0/ 209.0  ) * f);" << endl;
00210     }
00211 };
00212 
00213 
00215 struct interp_spline64
00216 {
00217     // size of neighbourhood
00218     static const int size = 8;
00219 
00221     void calc_coeff(double x, double * w) const
00222         {
00223             w[7] = ((  1.0/41.0 * x -   45.0/2911.0) * x -   26.0/2911.0) * x;
00224             w[6] = ((- 6.0/41.0 * x +  270.0/2911.0) * x +  156.0/2911.0) * x;
00225             w[5] = (( 24.0/41.0 * x - 1080.0/2911.0) * x -  624.0/2911.0) * x;
00226             w[4] = ((-49.0/41.0 * x + 4050.0/2911.0) * x + 2340.0/2911.0) * x;
00227             w[3] = (( 49.0/41.0 * x - 6387.0/2911.0) * x -    3.0/2911.0) * x + 1.0;
00228             w[2] = ((-24.0/41.0 * x + 4032.0/2911.0) * x - 2328.0/2911.0) * x;
00229             w[1] = ((  6.0/41.0 * x - 1008.0/2911.0) * x +  582.0/2911.0) * x;
00230             w[0] = ((- 1.0/41.0 * x +  168.0/2911.0) * x -   97.0/2911.0) * x;
00231         }
00232 
00233     void emitGLSL(std::ostringstream& oss) const {
00234         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
00235             << "                                              : (((- 6.0/41.0 * f +  270.0/2911.0) * f +  156.0/2911.0) * f)" << endl
00236             << "                                 : (i == 5.0) ? ((( 24.0/41.0 * f - 1080.0/2911.0) * f -  624.0/2911.0) * f)" << endl
00237             << "                                              : (((-49.0/41.0 * f + 4050.0/2911.0) * f + 2340.0/2911.0) * f)" << endl
00238             << "                     : (i > 1.0) ? (i == 3.0) ? ((( 49.0/41.0 * f - 6387.0/2911.0) * f -    3.0/2911.0) * f + 1.0)" << endl
00239             << "                                              : (((-24.0/41.0 * f + 4032.0/2911.0) * f - 2328.0/2911.0) * f)" << endl
00240             << "                                 : (i == 1.0) ? (((  6.0/41.0 * f - 1008.0/2911.0) * f +  582.0/2911.0) * f)" << endl
00241             << "                                              : (((- 1.0/41.0 * f +  168.0/2911.0) * f -   97.0/2911.0) * f);" << endl;
00242     }
00243 };
00244 
00246 template <int size_>
00247 struct interp_sinc
00248 {
00249     // size of neighbourhood
00250     static const int size = size_;
00251 
00253     void calc_coeff(double x, double * w) const
00254         {
00255             int idx;
00256             double xadd;
00257             for( idx = 0, xadd = size / 2 - 1.0 + x;
00258                  idx < size / 2;
00259                  xadd-=1.0)
00260             {
00261                 w[idx++] = sinc( xadd ) * sinc( xadd / ( size / 2 ));
00262             }
00263             for( xadd = 1.0 - x;
00264                  idx < size;
00265                  xadd+=1.0)
00266             {
00267                 w[idx++] = sinc( xadd ) * sinc( xadd / ( size / 2 ));
00268             }
00269         }
00270 
00271     void emitGLSL(std::ostringstream& oss) const {
00272         oss << "    float c = (i < " << (size/2.0) << ") ? 1.0 : -1.0;" << endl
00273             << "    float x = c * (" << (size/2 - 1.0) << " - i + f);" << endl
00274             << "    vec2 xpi = vec2(x, x / " << (size/2.0) << ") * " << M_PI << ";" << endl
00275             << "    vec2 xsin = sin(xpi);" << endl
00276             << "    vec2 result = vec2(1.0, 1.0);" << endl
00277             << "    if (xpi.x != 0.0) result.x = xsin.x / xpi.x;" << endl
00278             << "    if (xpi.y != 0.0) result.y = xsin.y / xpi.y;" << endl
00279             << "    return result.x * result.y;" << endl;
00280     }
00281 };
00282 
00283 
00288 template <typename SrcImageIterator, typename SrcAccessor,
00289           typename INTERPOLATOR>
00290 class ImageInterpolator
00291 {
00292 public:
00293     typedef typename SrcAccessor::value_type PixelType;
00294     // dummy mask type to be compatible to algorithms expecting a ImageMaskInterpolator object
00295     typedef typename vigra::UInt8 MaskType;
00296 private:
00297     typedef typename vigra::NumericTraits<PixelType>::RealPromote RealPixelType;
00298 
00299     SrcImageIterator m_sIter;
00300     SrcAccessor m_sAcc;
00301     int m_w;
00302     int m_h;
00303     bool m_warparound;
00304 
00305     INTERPOLATOR m_inter;
00306 
00307 public:
00309     ImageInterpolator(vigra::triple<SrcImageIterator, SrcImageIterator,SrcAccessor> const & src,
00310                           INTERPOLATOR & inter,
00311                           bool warparound)
00312     : m_sIter(src.first),
00313       m_sAcc(src.third),
00314       m_w(src.second.x - src.first.x),
00315       m_h(src.second.y - src.first.y),
00316       m_warparound(warparound),
00317       m_inter(inter)
00318     {
00319     }
00320 
00324     ImageInterpolator(SrcImageIterator src_upperleft,
00325                           SrcImageIterator src_lowerright,
00326                           SrcAccessor sa,
00327                           INTERPOLATOR & inter,
00328                           bool warparound)
00329     : m_sIter(src_upperleft),
00330       m_sAcc(sa),
00331       m_w(src_lowerright.x - src_upperleft.x),
00332       m_h(src_lowerright.y - src_upperleft.y),
00333       m_warparound(warparound),
00334       m_inter(inter)
00335     {
00336     }
00337 
00339     bool operator()(double x, double y,
00340                     PixelType & result, MaskType & mask) const
00341     {
00342         mask = 255;
00343         return operator()(x,y, result);
00344     }
00345 
00347     bool operator()(double x, double y,
00348                     PixelType & result) const
00349     {
00350 
00351         // skip all further interpolation if we cannot interpolate anything
00352         if (x < -INTERPOLATOR::size/2 || x > m_w + INTERPOLATOR::size/2) return false;
00353         if (y < -INTERPOLATOR::size/2 || y > m_h + INTERPOLATOR::size/2) return false;
00354 
00355         double t = floor(x);
00356         double dx = x - t;
00357         int srcx = int(t);
00358         t = floor(y);
00359         double dy = y - t;
00360         int srcy = int(t);
00361 
00362         if ( srcx > INTERPOLATOR::size/2 && srcx < m_w -INTERPOLATOR::size/2 &&
00363              srcy > INTERPOLATOR::size/2 && srcy < m_h - INTERPOLATOR::size/2)
00364         {
00365             return interpolateNoMaskInside(srcx, srcy, dx, dy, result);
00366         }
00367 
00368         double wx[INTERPOLATOR::size];
00369         double wy[INTERPOLATOR::size];
00370 
00371         // calculate x interpolation coefficients
00372         m_inter.calc_coeff(dx, wx);
00373         m_inter.calc_coeff(dy, wy);
00374 
00375         RealPixelType p(vigra::NumericTraits<RealPixelType>::zero());
00376         double weightsum = 0.0;
00377         for (int ky = 0; ky < INTERPOLATOR::size; ky++) {
00378             int bounded_ky = srcy + 1 + ky - INTERPOLATOR::size/2;
00379 
00380             // Boundary condition: do not replicate top and bottom
00381             if (bounded_ky < 0 || bounded_ky >= m_h) {
00382                 continue;
00383             }
00384 
00385             for (int kx = 0; kx < INTERPOLATOR::size; kx++) {
00386                 int bounded_kx = srcx + 1 + kx - INTERPOLATOR::size/2;
00387 
00388                 if (m_warparound) {
00389                     // Boundary condition: wrap around the image.
00390                     if (bounded_kx < 0) 
00391                         bounded_kx += m_w;
00392                     if (bounded_kx >= m_w) 
00393                         bounded_kx -= m_w;
00394                 } else {
00395                     // Boundary condition: replicate first and last column.
00396                     //                if (srcx + kx < 0) bounded_kx -= (srcx + kx);
00397                     //                if (srcx + kx >= src_w) bounded_kx -= (srcx + kx - (src_w - 1));
00398                     // Boundary condition: do not replicate left and right
00399                     if (bounded_kx < 0) 
00400                         continue;
00401                     if (bounded_kx >= m_w)
00402                         continue;
00403                 }
00404 
00405                 // check mask
00406                 double f = wx[kx]*wy[ky];
00407                 p += f * m_sAcc(m_sIter, vigra::Diff2D(bounded_kx, bounded_ky));
00408                 weightsum += f;
00409             }
00410         }
00411 
00412         // force a certain weight
00413         if (weightsum <= 0.2) return false;
00414         // Adjust filter for any ignored transparent pixels.
00415         if (weightsum != 1.0) p /= weightsum;
00416 
00417         result = vigra::detail::RequiresExplicitCast<PixelType>::cast(p);
00418         return true;
00419     }
00420 
00421 
00423     bool interpolateNoMaskInside(int srcx, int srcy, double dx, double dy,
00424                                     PixelType & result) const
00425     {
00426         double w[INTERPOLATOR::size];
00427         RealPixelType resX[INTERPOLATOR::size];
00428 
00429         // calculate x interpolation coefficients
00430         m_inter.calc_coeff(dx, w);
00431 
00432         RealPixelType p;
00433 
00434         // first pass of separable filter, x pass
00435         vigra::Diff2D offset(srcx - INTERPOLATOR::size/2 + 1,
00436                              srcy - INTERPOLATOR::size/2 + 1);
00437         SrcImageIterator ys(m_sIter + offset);
00438         for (int ky = 0; ky < INTERPOLATOR::size; ky++, ++(ys.y)) {
00439             p = vigra::NumericTraits<RealPixelType>::zero();
00440             typename SrcImageIterator::row_iterator xs(ys.rowIterator());
00441             //SrcImageIterator xs(ys);
00442             for (int kx = 0; kx < INTERPOLATOR::size; kx++, ++xs) {
00443                 p += w[kx] * m_sAcc(xs);
00444             }
00445             resX[ky] = p;
00446         }
00447 
00448         // y pass.
00449         m_inter.calc_coeff(dy, w);
00450         p = vigra::NumericTraits<RealPixelType>::zero();
00451         for (int ky = 0; ky < INTERPOLATOR::size; ky++) {
00452             p += w[ky] * resX[ky];
00453         }
00454 
00455         result = vigra::detail::RequiresExplicitCast<PixelType>::cast(p);
00456         return true;
00457     }
00458 
00459     void emitGLSL(std::ostringstream& oss) const {
00460         m_inter.emitGLSL(oss);
00461     }
00462 
00463 };
00464 
00465 
00471 template <typename SrcImageIterator, typename SrcAccessor,
00472           typename MaskIterator, typename MaskAccessor,
00473           typename INTERPOLATOR>
00474 class ImageMaskInterpolator
00475 {
00476 public:
00477     typedef typename SrcAccessor::value_type PixelType;
00478     typedef typename MaskAccessor::value_type MaskType;
00479 private:
00480     typedef typename vigra::NumericTraits<PixelType>::RealPromote RealPixelType;
00481 
00482     SrcImageIterator m_sIter;
00483     SrcAccessor m_sAcc;
00484     MaskIterator m_mIter;
00485     MaskAccessor m_mAcc;
00486     int m_w;
00487     int m_h;
00488     bool m_warparound;
00489 
00490     INTERPOLATOR m_inter;
00491 
00492 public:
00493 
00495     ImageMaskInterpolator(vigra::triple<SrcImageIterator, SrcImageIterator,SrcAccessor> const & src,
00496                           std::pair<MaskIterator, MaskAccessor> mask,
00497                           INTERPOLATOR & inter,
00498                           bool warparound)
00499     : m_sIter(src.first),
00500       m_sAcc(src.third),
00501       m_mIter(mask.first),
00502       m_mAcc(mask.second),
00503       m_w(src.second.x - src.first.x),
00504       m_h(src.second.y - src.first.y),
00505       m_warparound(warparound),
00506       m_inter(inter)
00507     {
00508     }
00509 
00513     ImageMaskInterpolator(SrcImageIterator src_upperleft,
00514                           SrcImageIterator src_lowerright,
00515                           SrcAccessor sa,
00516                           MaskIterator mask_upperleft,
00517                           MaskAccessor ma,
00518                           INTERPOLATOR & inter,
00519                           bool warparound)
00520     : m_sIter(src_upperleft),
00521       m_sAcc(sa),
00522       m_mIter(mask_upperleft),
00523       m_mAcc(ma),
00524       m_w(src_lowerright.x - src_upperleft.x),
00525       m_h(src_lowerright.y - src_upperleft.y),
00526       m_warparound(warparound),
00527       m_inter(inter)
00528     {
00529     }
00530 #if 0
00531 
00550 //    bool operator()(float x, float y,
00551     // this is slower than the full version, thanks to the normalized interpolation (with masks).
00552     bool interpolateSeperable(float x, float y,
00553                      PixelType & result) const
00554     {
00555 
00556         // skip all further interpolation if we cannot interpolate anything
00557         if (x < -INTERPOLATOR::size/2 || x > m_w + INTERPOLATOR::size/2) return false;
00558         if (y < -INTERPOLATOR::size/2 || y > m_h + INTERPOLATOR::size/2) return false;
00559 
00560         double t = floor(x);
00561         double dx = x - t;
00562         int srcx = int(t);
00563         t = floor(y);
00564         double dy = y - t;
00565         int srcy = int(t);
00566 
00567 
00568         double w[INTERPOLATOR::size];
00569 
00570         double weightsX[INTERPOLATOR::size];
00571         PixelType resX[INTERPOLATOR::size];
00572 
00573         // calculate x interpolation coefficients
00574         m_inter.calc_coeff(dx, w);
00575 
00576         // first pass of separable filter
00577 
00578         for (int ky = 0; ky < INTERPOLATOR::size; ky++) {
00579             int bounded_ky = srcy + 1 + ky - INTERPOLATOR::size/2;
00580 
00581             // Boundary condition: replicate top and bottom rows.
00582             //                 if (srcy + ky < 0) bounded_ky -= (srcy + ky);
00583             //                 if (srcy + ky >= src_h) bounded_ky -= (srcy + ky - (src_h - 1));
00584 
00585             // Boundary condition: do not replicate top and bottom
00586             if (bounded_ky < 0 || bounded_ky >= m_h) {
00587                 weightsX[ky] = 0;
00588                 resX[ky] = 0;
00589                 continue;
00590             }
00591 
00592             RealPixelType p(vigra::NumericTraits<RealPixelType>::zero());
00593             double weightsum = 0.0;
00594 
00595             for (int kx = 0; kx < INTERPOLATOR::size; kx++) {
00596                 int bounded_kx = srcx + 1 + kx - INTERPOLATOR::size/2;
00597 
00598                 if (m_warparound) {
00599                     // Boundary condition: wrap around the image.
00600                     if (bounded_kx < 0) 
00601                         bounded_kx += m_w;
00602                     if (bounded_kx >= m_w) 
00603                         bounded_kx -= m_w;
00604                 } else {
00605                     // Boundary condition: replicate first and last column.
00606                     //                if (srcx + kx < 0) bounded_kx -= (srcx + kx);
00607                     //                if (srcx + kx >= src_w) bounded_kx -= (srcx + kx - (src_w - 1));
00608                     // Boundary condition: do not replicate left and right
00609                     if (bounded_kx < 0) 
00610                         continue;
00611                     if (bounded_kx >= m_w)
00612                         continue;
00613                 }
00614                 if (m_mIter(bounded_kx, bounded_ky)) {
00615                     // check mask
00616                     p += w[kx] * m_sIter(bounded_kx, bounded_ky);
00617                     weightsum += w[kx];
00618                 }
00619             }
00620             weightsX[ky] = weightsum;
00621             resX[ky] = p;
00622         }
00623 
00624         // y pass.
00625         m_inter.calc_coeff(dy, w);
00626         RealPixelType p(vigra::NumericTraits<RealPixelType>::zero());
00627         double weightsum = 0.0;
00628         for (int ky = 0; ky < INTERPOLATOR::size; ky++) {
00629             weightsum += weightsX[ky] * w[ky];
00630             p += w[ky] * resX[ky];
00631         }
00632 
00633         if (weightsum == 0.0) return false;
00634         // Adjust filter for any ignored transparent pixels.
00635         if (weightsum != 1.0) p /= weightsum;
00636 
00637         result = vigra::detail::RequiresExplicitCast<PixelType>::cast(p);
00638         return true;
00639     }
00640 #endif
00641 
00661     bool operator()(double x, double y,
00662                     PixelType & result, MaskType & mask) const
00663     {
00664 
00665         // skip all further interpolation if we cannot interpolate anything
00666         if (x < -INTERPOLATOR::size/2 || x > m_w + INTERPOLATOR::size/2) return false;
00667         if (y < -INTERPOLATOR::size/2 || y > m_h + INTERPOLATOR::size/2) return false;
00668 
00669         double t = floor(x);
00670         double dx = x - t;
00671         int srcx = int(t);
00672         t = floor(y);
00673         double dy = y - t;
00674         int srcy = int(t);
00675 
00676         if ( srcx > INTERPOLATOR::size/2 && srcx < m_w -INTERPOLATOR::size/2 &&
00677              srcy > INTERPOLATOR::size/2 && srcy < m_h - INTERPOLATOR::size/2)
00678         {
00679             return interpolateInside(srcx, srcy, dx, dy, result, mask);
00680         }
00681 
00682         double wx[INTERPOLATOR::size];
00683         double wy[INTERPOLATOR::size];
00684 
00685         // calculate x interpolation coefficients
00686         m_inter.calc_coeff(dx, wx);
00687         m_inter.calc_coeff(dy, wy);
00688 
00689         // first pass of separable filter
00690 
00691         RealPixelType p(vigra::NumericTraits<RealPixelType>::zero());
00692         double m = 0; 
00693         double weightsum = 0.0;
00694         for (int ky = 0; ky < INTERPOLATOR::size; ky++) {
00695             int bounded_ky = srcy + 1 + ky - INTERPOLATOR::size/2;
00696 
00697             // Boundary condition: do not replicate top and bottom
00698             if (bounded_ky < 0 || bounded_ky >= m_h) {
00699                 continue;
00700             }
00701 
00702             for (int kx = 0; kx < INTERPOLATOR::size; kx++) {
00703                 int bounded_kx = srcx + 1 + kx - INTERPOLATOR::size/2;
00704 
00705                 if (m_warparound) {
00706                     // Boundary condition: wrap around the image.
00707                     if (bounded_kx < 0) 
00708                         bounded_kx += m_w;
00709                     if (bounded_kx >= m_w) 
00710                         bounded_kx -= m_w;
00711                 } else {
00712                     // Boundary condition: replicate first and last column.
00713                     //                if (srcx + kx < 0) bounded_kx -= (srcx + kx);
00714                     //                if (srcx + kx >= src_w) bounded_kx -= (srcx + kx - (src_w - 1));
00715                     // Boundary condition: do not replicate left and right
00716                     if (bounded_kx < 0) 
00717                         continue;
00718                     if (bounded_kx >= m_w)
00719                         continue;
00720                 }
00721 
00722                 MaskType cmask = m_mIter(bounded_kx, bounded_ky);
00723                 if (cmask) {
00724                     // check mask
00725                     double f = wx[kx]*wy[ky];
00726                     // TODO: check if this is good, influences the HDR stitching masks
00727                     m += f * cmask;
00728                     p += f * m_sAcc(m_sIter, vigra::Diff2D(bounded_kx, bounded_ky));
00729                     weightsum += f;
00730                 }
00731             }
00732         }
00733 
00734         // force a certain weight
00735         if (weightsum <= 0.2) return false;
00736         // Adjust filter for any ignored transparent pixels.
00737         if (weightsum != 1.0) {
00738             p /= weightsum;
00739             m /= weightsum;
00740         }
00741 
00742         mask = vigra::detail::RequiresExplicitCast<MaskType>::cast(m);
00743         result = vigra::detail::RequiresExplicitCast<PixelType>::cast(p);
00744         return true;
00745     }
00746 
00747 
00749     bool interpolateInside(int srcx, int srcy, double dx, double dy,
00750                                     PixelType & result, MaskType & mask) const
00751     {
00752 
00753         double wx[INTERPOLATOR::size];
00754         double wy[INTERPOLATOR::size];
00755 
00756         // calculate x interpolation coefficients
00757         m_inter.calc_coeff(dx, wx);
00758         m_inter.calc_coeff(dy, wy);
00759 
00760         RealPixelType p(vigra::NumericTraits<RealPixelType>::zero());
00761         double weightsum = 0.0;
00762         double m = 0.0;
00763         vigra::Diff2D offset(srcx - INTERPOLATOR::size/2 + 1,
00764                              srcy - INTERPOLATOR::size/2 + 1);
00765         SrcImageIterator ys(m_sIter + offset);
00766         MaskIterator yms(m_mIter + offset);
00767         for (int ky = 0; ky < INTERPOLATOR::size; ky++, ++(ys.y), ++(yms.y)) {
00768 //            int bounded_ky = srcy + 1 + ky - INTERPOLATOR::size/2;
00769             typename SrcImageIterator::row_iterator xs(ys.rowIterator());
00770             typename MaskIterator::row_iterator xms(yms.rowIterator());
00771             for (int kx = 0; kx < INTERPOLATOR::size; kx++, ++xs, ++xms) {
00772 //                int bounded_kx = srcx + 1 + kx - INTERPOLATOR::size/2;
00773 
00774                 MaskType cmask = *xms;
00775                 if (cmask) {
00776                     // check mask
00777                     double f = wx[kx]*wy[ky];
00778                     // TODO: check if this is good, influences the HDR stitching masks
00779                     m += f * cmask;
00780                     p += f * m_sAcc(xs);
00781                     weightsum += f;
00782                 }
00783             }
00784         }
00785 
00786         // force a certain weight
00787         if (weightsum <= 0.2) return false;
00788         // Adjust filter for any ignored transparent pixels.
00789         if (weightsum != 1.0) {
00790             p /= weightsum;
00791             m /= weightsum;
00792         }
00793 
00794         result = vigra::detail::RequiresExplicitCast<PixelType>::cast(p);
00795         mask = vigra::detail::RequiresExplicitCast<MaskType>::cast(m);
00796         return true;
00797     }
00798 
00799 };
00800 
00801 /********************************************************/
00802 /*                                                      */
00803 /*                  InterpolatingAccessor               */
00804 /*                                                      */
00805 /********************************************************/
00806 
00856 template <class ACCESSOR, class VALUETYPE, class INTERPOLATOR>
00857 class InterpolatingAccessor
00858 {
00859 public:
00862     typedef VALUETYPE value_type;
00863 
00866     InterpolatingAccessor(ACCESSOR a, INTERPOLATOR inter)
00867         : a_(a), inter_x(inter), inter_y(inter)
00868     {}
00869 
00886     template <class ITERATOR>
00887     value_type operator()(ITERATOR const & i, float x, float y) const
00888     {
00889         int ix = int(x);
00890         int iy = int(y);
00891         float dx = x - ix;
00892         float dy = y - iy;
00893         double wx[INTERPOLATOR::size];
00894         double wy[INTERPOLATOR::size];
00895 
00896         // promote value_type for multiplication
00897         typename vigra::NumericTraits<value_type>::RealPromote
00898             ret (vigra::NumericTraits<value_type>::zero());
00899 
00900         // calculate interpolation coefficients
00901         inter_x.calc_coeff(dx, wx);
00902         inter_y.calc_coeff(dy, wy);
00903 
00904         ITERATOR ys(i + vigra::Diff2D(ix - inter_x.size/2 + 1,
00905                                       iy - inter_y.size/2 + 1));
00906         for(int y = 0; y < inter_y.size; ++y, ++ys.y) {
00907             ITERATOR xs(ys);
00908             for(int x = 0; x < inter_x.size; x++, ++xs.x) {
00909                 ret += wx[x] * wy[y] * a_(xs);
00910             }
00911         }
00912         return vigra::detail::RequiresExplicitCast<value_type>::cast(ret);
00913     }
00914 
00936     template <class ITERATOR, class ALPHAITERATOR, class ALPHAACCESSOR>
00937     bool operator()(ITERATOR const & i, std::pair<ALPHAITERATOR, ALPHAACCESSOR> const & alpha,
00938                     float x, float y, value_type & result) const
00939     {
00940         int ix = int(x);
00941         int iy = int(y);
00942         float dx = x - ix;
00943         float dy = y - iy;
00944         double wx[INTERPOLATOR::size];
00945         double wy[INTERPOLATOR::size];
00946 
00947         // promote value_type for multiplication
00948         typename vigra::NumericTraits<value_type>::RealPromote
00949             ret (vigra::NumericTraits<value_type>::zero());
00950 
00951         // calculate interpolation coefficients
00952         inter_x.calc_coeff(dx, wx);
00953         inter_y.calc_coeff(dy, wy);
00954 
00955         ITERATOR ys(i + vigra::Diff2D(ix - inter_x.size/2 + 1,
00956                                       iy - inter_y.size/2 + 1));
00957         ALPHAITERATOR ays(alpha.first + vigra::Diff2D(ix - inter_x.size/2 + 1,
00958                                       iy - inter_y.size/2 + 1));
00959         for(int y = 0; y < inter_y.size; ++y, ++ys.y, ++ays.y) {
00960             ITERATOR xs(ys);
00961             ALPHAITERATOR axs(ays);
00962             for(int x = 0; x < inter_x.size; x++, ++xs.x, ++axs.x) {
00963                 if (alpha.second(axs) <= 0 ) {
00964                     return false;
00965                 }
00966                 ret += wx[x] * wy[y] * a_(xs);
00967             }
00968         }
00969         result = vigra::detail::RequiresExplicitCast<value_type>::cast(ret);
00970         return true;
00971     }
00972 
00973 
00974 private:
00975     ACCESSOR a_;
00976     INTERPOLATOR inter_x, inter_y;
00977 };
00978 
00979 } // namespace
00980 
00981 #endif // VIGRA_EXT_INTERPOLATORS_H
Generated on Fri May 24 01:25:39 2013 for Hugintrunk by  doxygen 1.3.9.1