Main Page | Modules | Namespace List | Class Hierarchy | Class List | File List | Namespace Members | Class Members | File Members | Related Pages
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
