/* $Id: cpl_image_resample_body.h,v 1.14 2007/11/13 14:27:40 yjung Exp $ * * This file is part of the ESO Common Pipeline Library * Copyright (C) 2001-2004 European Southern Observatory * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ /* Type dependent macros */ #if CPL_CLASS == CPL_CLASS_DOUBLE #define CPL_TYPE double #define CPL_TYPE_T CPL_TYPE_DOUBLE #define CPL_TYPE_INTERPOLATE cpl_image_get_interpolated_double #elif CPL_CLASS == CPL_CLASS_FLOAT #define CPL_TYPE float #define CPL_TYPE_T CPL_TYPE_FLOAT #define CPL_TYPE_INTERPOLATE cpl_image_get_interpolated_float #elif CPL_CLASS == CPL_CLASS_INT #define CPL_TYPE int #define CPL_TYPE_T CPL_TYPE_INT #define CPL_TYPE_INTERPOLATE cpl_image_get_interpolated_int #else #undef CPL_TYPE #undef CPL_TYPE_T #define CPL_TYPE_INTERPOLATE #endif #if defined CPL_OPERATION && CPL_OPERATION == CPL_IMAGE_RESAMPLE_SUBSAMPLE case CPL_TYPE_T: { CPL_TYPE * pi; CPL_TYPE * po; out_im = cpl_image_new(new_nx, new_ny, CPL_TYPE_T) ; pi = (CPL_TYPE*)image->pixels ; po = (CPL_TYPE*)out_im->pixels ; for (j=0 ; jnx ; po[i+j*new_nx] = pi[pos] ; } } } break ; #elif defined CPL_OPERATION && CPL_OPERATION == CPL_IMAGE_WARP case CPL_TYPE_T: /* Double loop on the output image */ for (j=0 ; j < out->ny; j++) { for (i=0 ; i < out->nx; i++) { /* Compute the original source for this pixel */ pos = i + j * out->nx ; x = i+1 - pdeltax[pos] ; y = j+1 - pdeltay[pos] ; /* Compute the new value */ value = CPL_TYPE_INTERPOLATE(in, in->nx, in->ny, xtabsperpix, ytabsperpix, x, y, yweight, pxprof, xradius, pyprof, yradius, sqyradius, sqyxratio, &confidence); if (confidence > 0) cpl_image_set(out, i+1, j+1, value); else { pbad[i + j * out->nx] = CPL_BINARY_1; hasbad = 1; } } } break; #elif defined CPL_OPERATION && CPL_OPERATION == CPL_IMAGE_WARP_POLYNOMIAL case CPL_TYPE_T: /* Double loop on the output image */ for (j=0, pval[1] = 1; j < out->ny; j++, pval[1]++) { for (i=0, pval[0] = 1; i < out->nx; i++, pval[0]++) { /* The CPL-calls here cannot fail */ /* Compute the original source for this pixel */ x = cpl_polynomial_eval(poly_u, val); y = cpl_polynomial_eval(poly_v, val); value = CPL_TYPE_INTERPOLATE(in, in->nx, in->ny, xtabsperpix, ytabsperpix, x, y, yweight, pxprof, xradius, pyprof, yradius, sqyradius, sqyxratio, &confidence); if (confidence > 0) cpl_image_set(out, i+1, j+1, value); else { pbad[i + j * out->nx] = CPL_BINARY_1; hasbad = 1; } } } break; #else #ifdef ADDTYPE #undef ADDTYPE #endif #define ADDTYPE(a) CONCAT2X(a,CPL_TYPE) /*----------------------------------------------------------------------------*/ /** @brief Interpolate a pixel @see cpl_image_get_interpolated() For efficiency reasons the caller must ensure that all arguments are valid. Specifically, this function must be called with an image of the correct pixel type. */ /*----------------------------------------------------------------------------*/ inline static double ADDTYPE(cpl_image_get_interpolated)( const cpl_image * source, int nx, int ny, double xtabsperpix, double ytabsperpix, double xpos, double ypos, double * yweight, const double * pxprof, double xradius, const double * pyprof, double yradius, double sqyradius, double sqyxratio, double * pconfid) { double sum = 0; double wsum = 0; double wconf = 0; double wabs = 0; double xdel; double sqrest; double xweight; double ydel; double value, weight; const CPL_TYPE * pixel = source->pixels; const int xfirst = ceil(xpos - xradius); const int xlast = floor(xpos + xradius); const int yfirst = ceil(ypos - yradius); const int ylast = floor(ypos + yradius); cpl_flops skipped = 0; int x, y; #if 0 assert( xradius > 0 ); assert( yradius > 0 ); assert( source != NULL ); assert( yweight != NULL ); assert( pxprof != NULL ); assert( pyprof != NULL ); assert( pconfid != NULL ); #endif if (xfirst > xlast || yfirst > ylast) { /* For sufficiently small radii no pixels can be reached */ *pconfid = 0; return 0; } /* Generate weights for all points in the Y-direction */ for (y = yfirst; y <= ylast; y++) yweight[y-yfirst] = pyprof[(int)(fabs(y - ypos) * ytabsperpix+0.5)]; /* Iterate through all the points in the X-direction */ for (x = xfirst; x <= xlast; x++) { xdel = x - xpos; sqrest = sqyradius - xdel*xdel*sqyxratio; xweight = pxprof[(int)(fabs(xdel) * xtabsperpix+0.5)]; /* Iterate through all the points in the Y-direction */ for (y = yfirst; y <= ylast; y++) { ydel = y - ypos; if (ydel*ydel > sqrest) { skipped += 8; continue; } /* (x,y) is within the ellipse with semi-major and -minor (xradius, yradius) at (xpos, ypos) */ /* The pixel weight is the product of the two 1D-profiles */ weight = xweight * yweight[y-yfirst]; /* Sum of all absolute weights in the inclusion area */ wconf += fabs(weight); if (x < 1 || y < 1 || x > nx || y > ny) { skipped += 5; continue; } if (cpl_image_is_rejected(source, x, y)) { skipped += 5; continue; } value = pixel[(x-1) + (y-1) * nx]; /* Sum of weigths actually used */ wsum += weight; /* Sum of absolute weigths actually used */ wabs += fabs(weight); /* Weighted pixel sum */ sum += value * weight; } } /* The confidence is the ratio of the absolute weights used over the sum of all absolute weights inside the inclusion area */ *pconfid = wconf > 0 ? wabs / wconf : 0; cpl_tools_add_flops( 3 * ( ylast - yfirst + 1) + 6 * ( xlast - xfirst + 1) + 11* ( ylast - yfirst + 1) * ( xlast - xfirst + 1) - skipped ); return wsum > 0 ? sum / wsum : 0; } #endif #undef CPL_TYPE #undef CPL_TYPE_T #undef CPL_TYPE_INTERPOLATE