/* $Id: cpl_image_resample.c,v 1.29 2007/11/14 09:07:47 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 */ /* * $Author: yjung $ * $Date: 2007/11/14 09:07:47 $ * $Revision: 1.29 $ * $Name: $ */ #ifdef HAVE_CONFIG_H #include #endif /*----------------------------------------------------------------------------- Includes -----------------------------------------------------------------------------*/ #include #include #include #include #include #include #include "cpl_tools.h" #include "cpl_memory.h" #include "cpl_image_resample.h" #include "cpl_mask.h" #include "cpl_image_defs.h" /*----------------------------------------------------------------------------- Defines -----------------------------------------------------------------------------*/ #define CPL_IMAGE_RESAMPLE_SUBSAMPLE 1 #define CPL_IMAGE_WARP 2 #define CPL_IMAGE_WARP_POLYNOMIAL 3 #define CONCAT(a,b) a ## _ ## b #define CONCAT2X(a,b) CONCAT(a,b) /*----------------------------------------------------------------------------*/ /** * @defgroup cpl_image_resample Images functionalities that need resampling * * This module provides functions to warp images. * * @par Synopsis: * @code * #include "cpl_image_resample.h" * @endcode */ /*----------------------------------------------------------------------------*/ /**@{*/ /*----------------------------------------------------------------------------- Private function prototypes -----------------------------------------------------------------------------*/ inline static double cpl_image_get_interpolated_double(const cpl_image *, int, int, double, double, double, double, double *, const double *, double, const double *, double, double, double, double *); inline static double cpl_image_get_interpolated_float(const cpl_image *, int, int, double, double, double, double, double *, const double *, double, const double *, double, double, double, double *); inline static double cpl_image_get_interpolated_int(const cpl_image *, int, int, double, double, double, double, double *, const double *, double, const double *, double, double, double, double *); /*----------------------------------------------------------------------------- Function codes -----------------------------------------------------------------------------*/ #define CPL_OPERATION CPL_IMAGE_WARP /*----------------------------------------------------------------------------*/ /** @brief Warp an image @param out Pre-allocated image to hold the result @param in Image to warp. @param deltax The x shifts for each pixel @param deltay The y shifts for each pixel @param xprofile Interpolation weight as a function of the distance in X @param xradius Positive inclusion radius in the X-dimension @param yprofile Interpolation weight as a function of the distance in Y @param yradius Positive inclusion radius in the Y-dimension @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error @see cpl_image_get_interpolated() 'out' and 'in' can have different sizes. The 'out' image pixels positions correspond to the lower left pixels positions in the 'in' image. In other words, the ID transformation would write the lower left part of 'in' in 'out'. deltax and deltay must be of type CPL_TYPE_DOUBLE and of the same size as out. For each pixel, the deltax and deltay are such that: @verbatim X = x + deltax Y = y + deltay @endverbatim where (X,Y) is the position in the 'out' image, (x,y) is the position in the 'in' image FIXME: Flux-calibration is missing. Beware that extreme transformations may lead to blank images. The input image type can be CPL_TYPE_INT, CPL_TYPE_FLOAT and CPL_TYPE_DOUBLE. Examples of profiles and radius are: @verbatim xprofile = cpl_vector_new(CPL_KERNEL_DEF_SAMPLES) ; cpl_vector_fill_kernel_profile(profile, CPL_KERNEL_DEFAULT, CPL_KERNEL_DEF_WIDTH) ; xradius = CPL_KERNEL_DEF_WIDTH ; @endverbatim Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL - CPL_ERROR_ILLEGAL_INPUT if the input images sizes are incompatible or if the delta images are not of type CPL_TYPE_DOUBLE - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_warp( cpl_image * out, const cpl_image * in, const cpl_image * deltax, const cpl_image * deltay, const cpl_vector * xprofile, double xradius, const cpl_vector * yprofile, double yradius) { double x, y ; cpl_mask * bad ; cpl_binary * pbad; double value, confidence; int hasbad ; const double sqxradius = xradius * xradius; const double sqyradius = yradius * yradius; double sqyxratio; double xtabsperpix; double ytabsperpix; const double * pxprof; const double * pyprof; double * yweight; int ixprolen; int iyprolen; int pos ; const double * pdeltax ; const double * pdeltay ; int i, j; /* Check entries */ cpl_ensure_code(out, CPL_ERROR_NULL_INPUT); cpl_ensure_code(in, CPL_ERROR_NULL_INPUT); cpl_ensure_code(deltax, CPL_ERROR_NULL_INPUT); cpl_ensure_code(deltay, CPL_ERROR_NULL_INPUT); cpl_ensure_code(cpl_image_get_size_x(deltax) == cpl_image_get_size_x(out), CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(cpl_image_get_size_y(deltax) == cpl_image_get_size_y(out), CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(cpl_image_get_size_x(deltay) == cpl_image_get_size_x(out), CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(cpl_image_get_size_y(deltay) == cpl_image_get_size_y(out), CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(cpl_image_get_type(deltax) == CPL_TYPE_DOUBLE, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(cpl_image_get_type(deltay) == CPL_TYPE_DOUBLE, CPL_ERROR_ILLEGAL_INPUT); ixprolen = cpl_vector_get_size(xprofile); cpl_ensure_code(ixprolen > 0, CPL_ERROR_ILLEGAL_INPUT); iyprolen = cpl_vector_get_size(yprofile); cpl_ensure_code(iyprolen > 0, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(xradius > 0, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(yradius > 0, CPL_ERROR_ILLEGAL_INPUT); /* Initialise */ bad = NULL ; hasbad = 0 ; pxprof = cpl_vector_get_data_const(xprofile); pyprof = cpl_vector_get_data_const(yprofile); sqyxratio = sqyradius/sqxradius; xtabsperpix = (ixprolen-1)/xradius; ytabsperpix = (iyprolen-1)/yradius; yweight = cpl_malloc((1 + (int)(2* yradius)) * sizeof(double)); /* Access the offsets */ pdeltax = cpl_image_get_data_double_const(deltax) ; pdeltay = cpl_image_get_data_double_const(deltay) ; /* Create the bad pixels mask */ bad = cpl_mask_new(out->nx, out->ny); pbad = cpl_mask_get_data(bad); switch (in->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_resample_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_resample_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_resample_body.h" #undef CPL_CLASS default: cpl_error_set(cpl_func, CPL_ERROR_TYPE_MISMATCH); } cpl_free(yweight); /* Handle bad pixels */ if (hasbad) cpl_image_reject_from_mask(out, bad); cpl_mask_delete(bad); return CPL_ERROR_NONE; } #undef CPL_OPERATION #define CPL_OPERATION CPL_IMAGE_WARP_POLYNOMIAL /*----------------------------------------------------------------------------*/ /** @brief Warp an image according to a 2d polynomial transformation. @param out Pre-allocated image to hold the result @param in Image to warp. @param poly_u 2D Polynomial transform in U. @param poly_v 2D Polynomial transform in V. @param xprofile Interpolation weight as a function of the distance in X @param xradius Positive inclusion radius in the X-dimension @param yprofile Interpolation weight as a function of the distance in Y @param yradius Positive inclusion radius in the Y-dimension @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error @see cpl_image_get_interpolated() 'out' and 'in' can have different sizes. The 'out' image pixels positions correspond to the lower left pixels positions in the 'in' image. In other words, the ID transformation would write the lower left part of 'in' in 'out'. Warp an image according to a polynomial transform. Provide two 2d cpl_polynomial poly_u and poly_v such as: @verbatim x = cpl_polynomial_eval(pu, (u, v)) y = cpl_polynomial_eval(pv, (u, v)) @endverbatim Attention! The polynomials define a reverse transform. (u,v) are coordinates in the warped image and (x,y) are coordinates in the original image. The transform you provide is used to compute from the warped image, which pixels contributed in the original image. FIXME: Flux-calibration is missing. Beware that extreme transformations may lead to blank images. The input image type can be CPL_TYPE_INT, CPL_TYPE_FLOAT and CPL_TYPE_DOUBLE. The two images may have different dimensions and types. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL - CPL_ERROR_ILLEGAL_INPUT if the polynomial dimensions are not 2 - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_warp_polynomial( cpl_image * out, const cpl_image * in, const cpl_polynomial * poly_u, const cpl_polynomial * poly_v, const cpl_vector * xprofile, double xradius, const cpl_vector * yprofile, double yradius) { double x, y; cpl_vector * val = NULL; double * pval; cpl_mask * bad = NULL; cpl_binary * pbad; double value, confidence; int i, j; int hasbad = 0; const double sqxradius = xradius * xradius; const double sqyradius = yradius * yradius; double sqyxratio; double xtabsperpix; double ytabsperpix; const double * pxprof; const double * pyprof; double * yweight; int ixprolen; int iyprolen; /* Check entries */ cpl_ensure_code(out, CPL_ERROR_NULL_INPUT); cpl_ensure_code(in, CPL_ERROR_NULL_INPUT); cpl_ensure_code(poly_u, CPL_ERROR_NULL_INPUT); cpl_ensure_code(poly_v, CPL_ERROR_NULL_INPUT); cpl_ensure_code(cpl_polynomial_get_dimension(poly_u) == 2, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(cpl_polynomial_get_dimension(poly_v) == 2, CPL_ERROR_ILLEGAL_INPUT); ixprolen = cpl_vector_get_size(xprofile); cpl_ensure(ixprolen > 0, cpl_error_get_code(), -3); iyprolen = cpl_vector_get_size(yprofile); cpl_ensure(iyprolen > 0, cpl_error_get_code(), -4); cpl_ensure(xradius > 0, CPL_ERROR_ILLEGAL_INPUT, -5); cpl_ensure(yradius > 0, CPL_ERROR_ILLEGAL_INPUT, -6); pxprof = cpl_vector_get_data_const(xprofile); pyprof = cpl_vector_get_data_const(yprofile); sqyxratio = sqyradius/sqxradius; xtabsperpix = (ixprolen-1)/xradius; ytabsperpix = (iyprolen-1)/yradius; yweight = cpl_malloc((1 + (int)(2* yradius)) * sizeof(double)); val = cpl_vector_new(2); pval = cpl_vector_get_data(val); bad = cpl_mask_new(out->nx, out->ny); pbad = cpl_mask_get_data(bad); switch (in->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_resample_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_resample_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_resample_body.h" #undef CPL_CLASS default: cpl_error_set(cpl_func, CPL_ERROR_TYPE_MISMATCH); } cpl_vector_delete(val); cpl_free(yweight); if (hasbad) cpl_image_reject_from_mask(out, bad); cpl_mask_delete(bad); return CPL_ERROR_NONE; } #undef CPL_OPERATION #define CPL_OPERATION CPL_IMAGE_RESAMPLE_SUBSAMPLE /*----------------------------------------------------------------------------*/ /** @brief Sub-sample an image @param image The image to subsample @param xstep Take every xstep pixel in x @param ystep Take every ystep pixel in y @return The newly allocated sub-sampled image or NULL in error case @see cpl_image_extract xstep and ystep shall be positive. Currently xstep must equal ystep. step represents the sampling step in x and y: a step of 2 will create an image with a qaurter of the pixels of the input image. Currently, both the X- and Y-size of the image must be divisible by the step. image type can be CPL_TYPE_INT, CPL_TYPE_FLOAT and CPL_TYPE_DOUBLE. The returned image must be deallocated using cpl_image_delete(). Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL - CPL_ERROR_ILLEGAL_INPUT if xstep, ystep are not as requested - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_extract_subsample( const cpl_image * image, int xstep, int ystep) { cpl_image * out_im; int new_nx, new_ny; int pos; int i, j; /* Check entries */ cpl_ensure(image != NULL, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(xstep == ystep, CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(xstep>0, CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(image->nx >= xstep, CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(image->ny >= xstep, CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(!(image->nx % xstep), CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(!(image->ny % xstep), CPL_ERROR_ILLEGAL_INPUT, NULL); /* Initialize */ new_nx = image->nx/xstep; new_ny = image->ny/xstep; /* Switch on the image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_resample_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_resample_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_resample_body.h" #undef CPL_CLASS default: cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL); } /* Sub sample the bad pixel map */ if (image->bpm != NULL) out_im->bpm = cpl_mask_extract_subsample(image->bpm, xstep, ystep) ; return out_im; } #undef CPL_OPERATION /*----------------------------------------------------------------------------*/ /** @brief Interpolate a pixel @param source Interpolation source @param xpos Pixel x floating-point position (FITS convention) @param ypos Pixel y floating-point position (FITS convention) @param xprofile Interpolation weight as a function of the distance in X @param xradius Positive inclusion radius in the X-dimension @param yprofile Interpolation weight as a function of the distance in Y @param yradius Positive inclusion radius in the Y-dimension @param pconfid Confidence level of the interpolated value (range 0 to 1) @return The interpolated pixel value, or undefined on error @see cpl_image_get() If the X- and Y-radii are identical the area of inclusion is a circle, otherwise it is an ellipse, with the larger of the two radii as the semimajor axis and the other as the semiminor axis. The radii are only required to be positive. However, for small radii, especially radii less than 1/sqrt(2), (xpos, ypos) may be located such that no source pixels are included in the interpolation, causing the interpolated pixel value to be undefined. The X- and Y-profiles can be generated with cpl_vector_fill_kernel_profile(profile, radius). For profiles generated with cpl_vector_fill_kernel_profile() it is important to use the same radius both there and in cpl_image_get_interpolated(). A good profile length is CPL_KERNEL_DEF_SAMPLES, using radius CPL_KERNEL_DEF_WIDTH. On error *pconfid is negative (unless pconfid is NULL). Otherwise, if *pconfid is zero, the interpolated pixel-value is undefined. Otherwise, if *pconfid is less than 1, the area of inclusion is close to the image border or contains rejected pixels. The input image type can be CPL_TYPE_INT, CPL_TYPE_FLOAT and CPL_TYPE_DOUBLE. Here is an example of a simple image unwarping (with error-checking omitted for brevity): const double xyradius = CPL_KERNEL_DEF_WIDTH; cpl_vector * xyprofile = cpl_vector_new(CPL_KERNEL_DEF_SAMPLES); cpl_image * unwarped = cpl_image_new(nx, ny, CPL_TYPE_DOUBLE); cpl_vector_fill_kernel_profile(xyprofile, CPL_KERNEL_DEFAULT, xyradius); for (iv = 1; iv <= ny; iv++) { for (iu = 1; iu <= nx; iu++) { double confidence; const double x = my_unwarped_x(); const double y = my_unwarped_y(); const double value = cpl_image_get_interpolated(warped, x, y, xyprofile, xyradius, xyprofile, xyradius, &confidence); if (confidence > 0) cpl_image_set(unwarped, iu, iv, value); else cpl_image_reject(unwarped, iu, iv); } } cpl_vector_delete(xyprofile); Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL - CPL_ERROR_ILLEGAL_INPUT if xradius, xprofile, yprofile and yradius are not as requested - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ double cpl_image_get_interpolated(const cpl_image * source, double xpos, double ypos, const cpl_vector * xprofile, double xradius, const cpl_vector * yprofile, double yradius, double * pconfid) { const double * pxprof = cpl_vector_get_data_const(xprofile); const double * pyprof = cpl_vector_get_data_const(yprofile); const double sqxradius = xradius * xradius; const double sqyradius = yradius * yradius; double sqyxratio; double xtabsperpix, ytabsperpix; const int nx = cpl_image_get_size_x(source); const int ny = cpl_image_get_size_y(source); const int ixprolen = cpl_vector_get_size(xprofile); const int iyprolen = cpl_vector_get_size(yprofile); double * yweight; double value; cpl_ensure(pconfid !=NULL, CPL_ERROR_NULL_INPUT, -1); *pconfid = -1; cpl_ensure(nx > 0, cpl_error_get_code(), -2); cpl_ensure(ixprolen > 0, cpl_error_get_code(), -3); cpl_ensure(iyprolen > 0, cpl_error_get_code(), -4); cpl_ensure(xradius > 0, CPL_ERROR_ILLEGAL_INPUT, -5); cpl_ensure(yradius > 0, CPL_ERROR_ILLEGAL_INPUT, -6); xtabsperpix = (ixprolen-1)/xradius; ytabsperpix = (iyprolen-1)/yradius; sqyxratio = sqyradius/sqxradius; yweight = cpl_malloc((1 + (int)(2* yradius)) * sizeof(double)); switch (source->type) { case CPL_TYPE_DOUBLE: value = cpl_image_get_interpolated_double(source, nx, ny, xtabsperpix, ytabsperpix, xpos, ypos, yweight, pxprof, xradius, pyprof, yradius, sqyradius, sqyxratio, pconfid); break; case CPL_TYPE_FLOAT: value = cpl_image_get_interpolated_float(source, nx, ny, xtabsperpix, ytabsperpix, xpos, ypos, yweight, pxprof, xradius, pyprof, yradius, sqyradius, sqyxratio, pconfid); break; case CPL_TYPE_INT: value = cpl_image_get_interpolated_int(source, nx, ny, xtabsperpix, ytabsperpix, xpos, ypos, yweight, pxprof, xradius, pyprof, yradius, sqyradius, sqyxratio, pconfid); break; default: value = 0; cpl_error_set(cpl_func, CPL_ERROR_TYPE_MISMATCH); } cpl_free(yweight); return value; } /**@}*/ #define CPL_CLASS CPL_CLASS_DOUBLE #undef CPL_OPERATION #include "cpl_image_resample_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #undef CPL_OPERATION #include "cpl_image_resample_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #undef CPL_OPERATION #include "cpl_image_resample_body.h" #undef CPL_CLASS