Index: trunk/psLib/src/imageops/Makefile.am =================================================================== --- trunk/psLib/src/imageops/Makefile.am (revision 12738) +++ trunk/psLib/src/imageops/Makefile.am (revision 12741) @@ -8,4 +8,5 @@ psImageConvolve.c \ psImageGeomManip.c \ + psImageInterpolate.c \ psImagePixelExtract.c \ psImagePixelManip.c \ @@ -22,4 +23,5 @@ psImageConvolve.h \ psImageGeomManip.h \ + psImageInterpolate.h \ psImagePixelExtract.h \ psImagePixelManip.h \ Index: trunk/psLib/src/imageops/psImageConvolve.c =================================================================== --- trunk/psLib/src/imageops/psImageConvolve.c (revision 12738) +++ trunk/psLib/src/imageops/psImageConvolve.c (revision 12741) @@ -7,6 +7,6 @@ /// @author Eugene Magnier, IfA /// -/// @version $Revision: 1.50 $ $Name: not supported by cvs2svn $ -/// @date $Date: 2007-03-02 22:19:21 $ +/// @version $Revision: 1.51 $ $Name: not supported by cvs2svn $ +/// @date $Date: 2007-04-04 22:42:02 $ /// /// Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -102,5 +102,5 @@ const psVector *xShifts, const psVector *yShifts, - bool tRelative, + float totalTime, bool xyRelative) { @@ -113,4 +113,9 @@ PS_ASSERT_VECTOR_TYPE(xShifts, PS_TYPE_S32, NULL); PS_ASSERT_VECTOR_TYPE(yShifts, PS_TYPE_S32, NULL); + + if (isnan(totalTime)) { + // It's more expensive to check for NAN than 0.0 + totalTime = 0.0; + } // If there are no shifts, the kernel is just a 1 at 0,0 @@ -153,5 +158,5 @@ } - if (tRelative) { + if (totalTime <= 0) { tSum += tShifts->data.F32[i]; } @@ -160,8 +165,8 @@ psTrace("psLib.imageops", 5, "Kernel range: %d:%d,%d:%d\n", xMin, xMax, yMin, yMax); - if (!tRelative) { + if (totalTime > 0) { // Then the total time is simply the final value // NB: We assume the counter starts at zero! - tSum = tShifts->data.F32[tShifts->n - 1]; + tSum = totalTime; } @@ -181,5 +186,5 @@ } float t = tShifts->data.F32[i]; - if (!tRelative) { + if (totalTime > 0) { t -= tLast; tLast = tShifts->data.F32[i]; Index: trunk/psLib/src/imageops/psImageConvolve.h =================================================================== --- trunk/psLib/src/imageops/psImageConvolve.h (revision 12738) +++ trunk/psLib/src/imageops/psImageConvolve.h (revision 12741) @@ -5,6 +5,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-03-02 22:19:21 $ + * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-04-04 22:42:02 $ * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii */ @@ -98,9 +98,9 @@ /// psKernel *psKernelGenerate( - const psVector *tShifts, ///< list of time shifts (F32) - const psVector *xShifts, ///< list of x-axis shifts (S32) - const psVector *yShifts, ///< list of y-axis shifts (S32) - bool tRelative, ///< Are times relative (durations) or absolute? - bool xyRelative ///< Are x,y positions relative (shifts) or absolute? + const psVector *tShifts, ///< list of time shifts (F32) + const psVector *xShifts, ///< list of x-axis shifts (S32) + const psVector *yShifts, ///< list of y-axis shifts (S32) + float totalTime, ///< Total time (relative times if negative) + bool xyRelative ///< Are x,y positions relative (shifts) or absolute? ); Index: trunk/psLib/src/imageops/psImageGeomManip.c =================================================================== --- trunk/psLib/src/imageops/psImageGeomManip.c (revision 12738) +++ trunk/psLib/src/imageops/psImageGeomManip.c (revision 12741) @@ -10,6 +10,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-03-14 00:39:50 $ + * @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-04-04 22:42:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -33,5 +33,5 @@ #include "psMemory.h" #include "psAssert.h" - +#include "psImageInterpolate.h" #include "psCoord.h" @@ -200,12 +200,4 @@ } - if (mode > PS_INTERPOLATE_LANCZOS4_VARIANCE ) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Specified interpolation mode, %d, is unsupported."), - mode); - psFree(out); - return NULL; - } - // create an output image of the same size // and type @@ -213,4 +205,7 @@ outCols = in->numCols * scale; invScale = 1.0f / (float)scale; + + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(mode, in, NULL, NULL, 0, + NAN, NAN, 0, 0, 0); #define PSIMAGE_RESAMPLE_CASE(TYPE) \ @@ -221,5 +216,12 @@ float inRow = (float)row * invScale; \ for (psS32 col=0;col PS_MAX_##TYPE || \ @@ -515,67 +523,44 @@ outRow = out->data.TYPE[y]; \ for (psS32 x = 0; x < outCols; x++) { \ - outRow[x] = p_psImagePixelInterpolate##MODE##_##TYPE(input,inX,inY,NULL,0,exposed); \ + double value; \ + if (!psImageInterpolate(&value, NULL, NULL, inX, inY, interp)) { \ + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); \ + psFree(out); \ + psFree(interp); \ + return NULL; \ + } \ + outRow[x] = value; \ inX += cosT; \ inY -= sinT; \ } \ } \ - } - - #define PSIMAGE_ROTATE_ARBITRARY_CASE(MODE) \ - case PS_INTERPOLATE_##MODE: \ - switch (type) { \ - case PS_TYPE_U8: \ - PSIMAGE_ROTATE_ARBITRARY_LOOP(U8,MODE); \ break; \ - case PS_TYPE_U16: \ - PSIMAGE_ROTATE_ARBITRARY_LOOP(U16,MODE); \ - break; \ - case PS_TYPE_U32: /* Not a requirement */ \ - PSIMAGE_ROTATE_ARBITRARY_LOOP(U32,MODE); \ - break; \ - case PS_TYPE_U64: /* Not a requirement */ \ - PSIMAGE_ROTATE_ARBITRARY_LOOP(U64,MODE); \ - break; \ - case PS_TYPE_S8: \ - PSIMAGE_ROTATE_ARBITRARY_LOOP(S8,MODE); \ - break; \ - case PS_TYPE_S16: \ - PSIMAGE_ROTATE_ARBITRARY_LOOP(S16,MODE); \ - break; \ - case PS_TYPE_S32: /* Not a requirement */ \ - PSIMAGE_ROTATE_ARBITRARY_LOOP(S32,MODE); \ - break; \ - case PS_TYPE_S64: /* Not a requirement */ \ - PSIMAGE_ROTATE_ARBITRARY_LOOP(S64,MODE); \ - break; \ - case PS_TYPE_F32: \ - PSIMAGE_ROTATE_ARBITRARY_LOOP(F32,MODE); \ - break; \ - case PS_TYPE_F64: \ - PSIMAGE_ROTATE_ARBITRARY_LOOP(F64,MODE); \ - break; \ - default: { \ - char* typeStr; \ - PS_TYPE_NAME(typeStr,type); \ - psError(PS_ERR_BAD_PARAMETER_TYPE, true, \ - _("Specified psImage type, %s, is not supported."), \ - typeStr); \ - psFree(out); \ - out = NULL; \ - } \ - } \ - break; - - switch (mode) { - PSIMAGE_ROTATE_ARBITRARY_CASE(FLAT); - PSIMAGE_ROTATE_ARBITRARY_CASE(BILINEAR); - PSIMAGE_ROTATE_ARBITRARY_CASE(BILINEAR_VARIANCE); - default: - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Specified interpolation mode, %d, is unsupported."), - mode); - psFree(out); - out = NULL; - } + } + + switch (type) { + PSIMAGE_ROTATE_ARBITRARY_LOOP(U8); + PSIMAGE_ROTATE_ARBITRARY_LOOP(U16); + PSIMAGE_ROTATE_ARBITRARY_LOOP(U32); + PSIMAGE_ROTATE_ARBITRARY_LOOP(U64); + PSIMAGE_ROTATE_ARBITRARY_LOOP(S8); + PSIMAGE_ROTATE_ARBITRARY_LOOP(S16); + PSIMAGE_ROTATE_ARBITRARY_LOOP(S32); + PSIMAGE_ROTATE_ARBITRARY_LOOP(S64); + PSIMAGE_ROTATE_ARBITRARY_LOOP(F32); + PSIMAGE_ROTATE_ARBITRARY_LOOP(F64); + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + _("Specified psImage type, %s, is not supported."), + typeStr); + psFree(out); + psFree(interp); + out = NULL; + } + } + + psFree(interp); + } @@ -695,5 +680,8 @@ out = psImageRecycle(out, outCols, outRows, type); - #define PSIMAGE_SHIFT_CASE(MODE,TYPE) \ + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(mode, input, NULL, NULL, 0, + exposed, NAN, 0, 0, 0.0); + + #define PSIMAGE_SHIFT_CASE(TYPE) \ case PS_TYPE_##TYPE: \ if (exposed < PS_MIN_##TYPE || \ @@ -719,50 +707,39 @@ for (psS32 col=0;coltype.type) { \ - PSIMAGE_SHIFT_CASE(MODE,U8); \ - PSIMAGE_SHIFT_CASE(MODE,U16); \ - PSIMAGE_SHIFT_CASE(MODE,U32); /* Not a requirement */ \ - PSIMAGE_SHIFT_CASE(MODE,U64); /* Not a requirement */ \ - PSIMAGE_SHIFT_CASE(MODE,S8); \ - PSIMAGE_SHIFT_CASE(MODE,S16); \ - PSIMAGE_SHIFT_CASE(MODE,S32); /* Not a requirement */ \ - PSIMAGE_SHIFT_CASE(MODE,S64); /* Not a requirement */ \ - PSIMAGE_SHIFT_CASE(MODE,F32); \ - PSIMAGE_SHIFT_CASE(MODE,F64); \ - \ - default: { \ - char* typeStr; \ - PS_TYPE_NAME(typeStr,type); \ - psError(PS_ERR_BAD_PARAMETER_TYPE, true, \ - _("Specified psImage type, %s, is not supported."), \ - typeStr); \ - psFree(out); \ - out = NULL; \ - } \ - } \ - break; - - // EAM: added BICUBE - switch (mode) { - PSIMAGE_SHIFT_ARBITRARY_CASE(FLAT); - PSIMAGE_SHIFT_ARBITRARY_CASE(BILINEAR); - PSIMAGE_SHIFT_ARBITRARY_CASE(BILINEAR_VARIANCE); - PSIMAGE_SHIFT_ARBITRARY_CASE(BICUBE); - default: - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Specified interpolation mode, %d, is unsupported."), - mode); - psFree(out); - out = NULL; - } - + switch (input->type.type) { + PSIMAGE_SHIFT_CASE(U8); + PSIMAGE_SHIFT_CASE(U16); + PSIMAGE_SHIFT_CASE(U32); + PSIMAGE_SHIFT_CASE(U64); + PSIMAGE_SHIFT_CASE(S8); + PSIMAGE_SHIFT_CASE(S16); + PSIMAGE_SHIFT_CASE(S32); + PSIMAGE_SHIFT_CASE(S64); + PSIMAGE_SHIFT_CASE(F32); + PSIMAGE_SHIFT_CASE(F64); + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, _("Specified psImage type, %s, is not supported."), + typeStr); + psFree(out); + psFree(interp); + return NULL; + } + } + + psFree(interp); return out; } @@ -866,10 +843,14 @@ // loop through the output image using the domain above and transform - // each output pixel to input coordinates and use psImagePixelInterpolate + // each output pixel to input coordinates and use psImageInterpolate // to determine the pixel value. psPlane outPosition; psPlane* inPosition = NULL; - #define PSIMAGE_TRANSFORM_DOTRANSFORM(TYPE,MODE) \ + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(mode, input, NULL, inputMask, + inputMaskVal, NAN, NAN, 0, 0, 0.0); + + + #define PSIMAGE_TRANSFORM_DOTRANSFORM(TYPE) \ /* apply the transform to get the position in the input image */ \ inPosition = psPlaneTransformApply(inPosition, outToIn, &outPosition); \ @@ -882,7 +863,11 @@ } \ /* interpolate the cooresponding input pixel to get the output pixel value. */ \ - ps##TYPE value = p_psImagePixelInterpolate##MODE##_##TYPE(input, \ - inPosition->x, inPosition->y, \ - inputMask, inputMaskVal, NAN); \ + double value; \ + if (!psImageInterpolate(&value, NULL, NULL, inPosition->x, inPosition->y, interp)) { \ + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); \ + psFree(output); \ + psFree(interp); \ + return NULL; \ + } \ /* psFree(inPosition); */\ if (isnan(value)) { \ @@ -893,5 +878,6 @@ } \ - #define PSIMAGE_TRANSFORM_LOOP(TYPE, MODE) { \ + #define PSIMAGE_TRANSFORM_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ for (int row = 0; row < numRows; row++) { \ outPosition.y = row+row0; \ @@ -899,55 +885,23 @@ for (int col = 0; col < numCols; col++) { \ outPosition.x = col+col0; \ - PSIMAGE_TRANSFORM_DOTRANSFORM(TYPE,MODE) \ + PSIMAGE_TRANSFORM_DOTRANSFORM(TYPE) \ outputData[col] = value; \ } \ } \ - } - - #define PSIMAGE_TRANSFORM_FROMLIST(TYPE, MODE) { \ - int n = pixels->n; \ - for (int i= 0; i < n; i++) { \ - int x = pixels->data[i].x; \ - int y = pixels->data[i].y; \ - if (x >= col0 && x < col1 && y >= row0 && y < row1) { \ - outPosition.x = x; \ - outPosition.y = y; \ - PSIMAGE_TRANSFORM_DOTRANSFORM(TYPE,MODE) \ - output->data.TYPE[y][x] = value; \ - } \ - } \ - } - - #define PSIMAGE_TRANSFORM_CASE(MODE) \ -case PS_INTERPOLATE_##MODE: \ - switch (type) { \ - case PS_TYPE_F32: \ - PSIMAGE_TRANSFORM_LOOP(F32,MODE); \ break; \ - case PS_TYPE_F64: \ - PSIMAGE_TRANSFORM_LOOP(F64,MODE); \ - break; \ - default: { \ - char* typeStr; \ - PS_TYPE_NAME(typeStr,type); \ - psError(PS_ERR_BAD_PARAMETER_TYPE, true, \ - _("Specified psImage type, %s, is not supported."), \ - typeStr); \ - psFree(output); \ - return NULL; \ - } \ - } \ - break; - - switch (mode) { - PSIMAGE_TRANSFORM_CASE(FLAT); - PSIMAGE_TRANSFORM_CASE(BILINEAR); - PSIMAGE_TRANSFORM_CASE(BILINEAR_VARIANCE); - default: - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Specified interpolation mode, %d, is unsupported."), - mode); - psFree(output); - return NULL; + } + + switch (type) { + PSIMAGE_TRANSFORM_CASE(F32); + PSIMAGE_TRANSFORM_CASE(F64); + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, _("Specified psImage type, %s, is not supported."), + typeStr); + psFree(output); + psFree(inPosition); + return NULL; + } } Index: trunk/psLib/src/imageops/psImageGeomManip.h =================================================================== --- trunk/psLib/src/imageops/psImageGeomManip.h (revision 12738) +++ trunk/psLib/src/imageops/psImageGeomManip.h (revision 12741) @@ -6,6 +6,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-03-14 00:39:50 $ + * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-04-04 22:42:02 $ * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii */ @@ -17,4 +17,5 @@ #include "psImage.h" +#include "psImageInterpolate.h" #include "psCoord.h" #include "psStats.h" @@ -158,5 +159,5 @@ * coordinates in the input image of a pixel in the output image — note that * this is the reverse of what might be naively expected, but it is what is - * required in order to use psImagePixelInterpolate. If the pixels array is + * required in order to use psImageInterpolate. If the pixels array is * non-NULL, it shall consist of psPixelCoords, and only those pixels in the * output image shall be transformed; otherwise, the entire image is Index: trunk/psLib/src/imageops/psImageInterpolate.c =================================================================== --- trunk/psLib/src/imageops/psImageInterpolate.c (revision 12741) +++ trunk/psLib/src/imageops/psImageInterpolate.c (revision 12741) @@ -0,0 +1,600 @@ +/** @file psImageInterpolate.c + * + * @brief Contains functions for interpolating an image + * + * @author Robert DeSonia, MHPCC + * @author Ross Harman, MHPCC + * @author Paul Price, IfA + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-04-04 22:42:02 $ + * + * Copyright 2004-2007 Institute for Astronomy, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +# include "config.h" +#endif + +#include +#include +#include +#include + +#include "psAbort.h" +#include "psMemory.h" +#include "psError.h" +#include "psAssert.h" +#include "psString.h" +#include "psImage.h" +#include "psImageInterpolate.h" + +static void imageInterpolateOptionsFree(psImageInterpolateOptions *options) +{ + // Casting away const + psFree((psImage*)options->image); + psFree((psImage*)options->mask); + psFree((psImage*)options->variance); +} + +psImageInterpolateOptions *psImageInterpolateOptionsAlloc(psImageInterpolateMode mode, + const psImage *image, const psImage *variance, + const psImage *mask, psMaskType maskVal, + double badImage, double badVariance, + psMaskType badMask, psMaskType poorMask, + float poorFrac) +{ + psImageInterpolateOptions *options = psAlloc(sizeof(psImageInterpolateOptions)); // Options, to return + psMemSetDeallocator(options, (psFreeFunc)imageInterpolateOptionsFree); + + options->mode = mode; + // Casting away const to add to options + options->image = psMemIncrRefCounter((psImage*)image); + options->variance = psMemIncrRefCounter((psImage*)variance); + options->mask = psMemIncrRefCounter((psImage*)mask); + options->maskVal = maskVal; + options->badImage = badImage; + options->badVariance = badVariance; + options->badMask = badMask; + options->poorMask = poorMask; + options->poorFrac = poorFrac; + + return options; +} + +// Interpolation engine for flat mode (nearest pixel) +static inline bool interpolateFlat(double *imageValue, double *varianceValue, psMaskType *maskValue, + float x, float y, const psImageInterpolateOptions *options) +{ + // Parameters have been checked by psImageInterpolate() + + const psImage *image = options->image; // Image of interest + int xInt = round(x - 0.5 + FLT_EPSILON); // Pixel closest to point of interest in x + int yInt = round(y - 0.5 + FLT_EPSILON); // Pixel closest to point of interest in y + int xLast = image->numCols - 1; // Last pixel in x + int yLast = image->numRows - 1; // Last pixel in y + + if (xInt < 0 || xInt > xLast || yInt < 0 || yInt > yLast) { + // At least one pixel of the interpolation kernel is off the image + if (imageValue) { + *imageValue = options->badImage; + } + if (varianceValue) { + *varianceValue = options->badVariance; + } + if (maskValue) { + *maskValue = options->badMask; + } + } else { + + // Image and variance 'interpolation' according to image type + #define FLAT_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + if (imageValue) { \ + *imageValue = options->image->data.TYPE[yInt][xInt]; \ + } \ + if (varianceValue) { \ + *varianceValue = options->variance->data.TYPE[yInt][xInt]; \ + } \ + break; \ + } + + switch (options->image->type.type) { + FLAT_CASE(U8); + FLAT_CASE(U16); + FLAT_CASE(U32); + FLAT_CASE(U64); + FLAT_CASE(S8); + FLAT_CASE(S16); + FLAT_CASE(S32); + FLAT_CASE(S64); + FLAT_CASE(F32); + FLAT_CASE(F64); + default: + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", + options->image->type.type); + return false; + } + + if (maskValue) { + *maskValue = options->mask->data.PS_TYPE_MASK_DATA[yInt][xInt]; + } + } + return true; +} + +// Interpolation engine using interpolation kernel +static bool interpolateKernel(double *imageValue, double *varianceValue, psMaskType *maskValue, + float x, float y, const psImageInterpolateOptions *options) +{ + // Parameters have been checked by psImageInterpolate() + + int xNum, yNum; // Number of interpolation kernel pixels + switch (options->mode) { + case PS_INTERPOLATE_BILINEAR: + xNum = yNum = 2; + break; + case PS_INTERPOLATE_BICUBE: + case PS_INTERPOLATE_GAUSS: + xNum = yNum = 3; + break; + case PS_INTERPOLATE_FLAT: + case PS_INTERPOLATE_LANCZOS2: + case PS_INTERPOLATE_LANCZOS3: + case PS_INTERPOLATE_LANCZOS4: + default: + psAbort("Invalid interpolation mode."); + } + + const psImage *image = options->image; // Image of interest + int xFloor = floor(x - 0.5 + FLT_EPSILON); // Pixel below point of interest in x + int yFloor = floor(y - 0.5 + FLT_EPSILON); // Pixel below point of interest in y + int xLast = image->numCols - 1; // Last pixel in x + int yLast = image->numRows - 1; // Last pixel in y + + if (xFloor - (xNum - 1) / 2 < 0 || xFloor + xNum / 2 > xLast || + yFloor - (yNum - 1) / 2 < 0 || yFloor + yNum / 2 > yLast) { + // At least one pixel of the interpolation kernel is off the image + if (imageValue) { + *imageValue = options->badImage; + } + if (varianceValue) { + *varianceValue = options->badVariance; + } + if (maskValue) { + *maskValue = options->badMask; + } + return true; + } + double kernel[yNum][xNum]; // Interpolation kernel for straight interpolation + switch (options->mode) { + case PS_INTERPOLATE_BILINEAR: { + double xFrac = x - 0.5 - xFloor; // Fraction of pixel in x + double yFrac = y - 0.5 - yFloor; // Fraction of pixel in y + kernel[0][0] = (1.0 - xFrac) * (1.0 - yFrac); + kernel[0][1] = xFrac * (1.0 - yFrac); + kernel[1][0] = yFrac * (1.0 - xFrac); + kernel[1][1] = xFrac * yFrac; + break; + } + case PS_INTERPOLATE_BICUBE: { + double xFrac = x - 0.5 - xFloor; // Fraction of pixel in x + double yFrac = y - 0.5 - yFloor; // Fraction of pixel in y + // Calculation variables + double xxFrac = xFrac * xFrac / 6.0; + double yyFrac = yFrac * yFrac / 6.0; + double xyFrac = 0.25 * xFrac * yFrac; + xFrac /= 6.0; + yFrac /= 6.0; + kernel[0][0] = - 1.0/9.0 - xFrac - yFrac + xxFrac + yyFrac + xyFrac; + kernel[0][1] = 2.0/9.0 - yFrac - 2.0 * xxFrac + yyFrac; + kernel[0][2] = - 1.0/9.0 + xFrac - yFrac + xxFrac + yyFrac - xyFrac; + kernel[1][0] = 2.0/9.0 - xFrac + xxFrac - 2.0 * yyFrac; + kernel[1][1] = 5.0/9.0 - 2.0 * xxFrac - 2.0 * yyFrac; + kernel[1][2] = 2.0/9.0 + xFrac + xxFrac - 2.0 * yyFrac; + kernel[2][0] = - 1.0/9.0 - xFrac + yFrac + xxFrac + yyFrac - xyFrac; + kernel[2][1] = 2.0/9.0 + yFrac - 2.0 * xxFrac + yyFrac; + kernel[2][2] = - 1.0/9.0 + xFrac + yFrac + xxFrac + yyFrac + xyFrac; + break; + } + case PS_INTERPOLATE_GAUSS: { + double xFrac = x - 0.5 - xFloor; // Fraction of pixel in x + double yFrac = y - 0.5 - yFloor; // Fraction of pixel in y + double xGaussFrac = 2.0 * erf((double)xNum / 4.0) - 1.0; // Fraction of Gaussian in x + double yGaussFrac = 2.0 * erf((double)yNum / 4.0) - 1.0; // Fraction of Gaussian in x + double norm = 1.0 / (double)xNum / (double)yNum * + (1.0 - xGaussFrac) / xGaussFrac * (1.0 - yGaussFrac) / yGaussFrac; // Normalisation + for (int j = 0, yPos = - (yNum - 1) / 2; j < yNum; j++, yPos++) { + for (int i = 0, xPos = - (xNum - 1) / 2; i < xNum; i++, xPos++) { + kernel[j][i] = norm * exp(-0.5 * (PS_SQR(xPos - xFrac) + PS_SQR(yPos - yFrac))); + } + } + break; + } + case PS_INTERPOLATE_FLAT: + case PS_INTERPOLATE_LANCZOS2: + case PS_INTERPOLATE_LANCZOS3: + case PS_INTERPOLATE_LANCZOS4: + default: + psAbort("Invalid interpolation mode."); + } + + // Image interpolation, according to image type + #define KERNEL_IMAGE_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) { \ + for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) { \ + value += values[j][i] = kernel[j][i] * image->data.TYPE[yPix][xPix]; \ + } \ + } \ + break; \ + } + + // Calculate the value for the image + double value = 0.0; // Value to return + double values[yNum][xNum]; // Values in image + switch (image->type.type) { + KERNEL_IMAGE_CASE(U8); + KERNEL_IMAGE_CASE(U16); + KERNEL_IMAGE_CASE(U32); + KERNEL_IMAGE_CASE(U64); + KERNEL_IMAGE_CASE(S8); + KERNEL_IMAGE_CASE(S16); + KERNEL_IMAGE_CASE(S32); + KERNEL_IMAGE_CASE(S64); + KERNEL_IMAGE_CASE(F32); + KERNEL_IMAGE_CASE(F64); + default: + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", + image->type.type); + return false; + } + if (imageValue) { + *imageValue = value; + } + + // Check the mask value + if (maskValue) { + const psImage *mask = options->mask; // Image mask + psMaskType maskVal = options->maskVal; // Mask value + double badValue = 0.0; // Amount of flux in bad pixels + int badPix = 0; // Number of bad pixels + *maskValue = 0; + for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) { + for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) { + if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { + badValue += values[j][i]; + badPix++; + } + *maskValue |= mask->data.PS_TYPE_MASK_DATA[yPix][xPix]; + } + } + + if (badPix > 0) { + if (fabs(badValue) >= fabs(value) * options->poorFrac) { + *maskValue |= options->badMask; + } else { + *maskValue |= options->poorMask; + } + } + } + + // Finally, the variance + if (varianceValue) { + const psImage *variance = options->variance; // Image variance + *varianceValue = 0.0; + + // Variance interpolation, according to image type + #define KERNEL_VARIANCE_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) { \ + for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) { \ + *varianceValue += PS_SQR(kernel[j][i]) * variance->data.TYPE[yPix][xPix]; \ + } \ + } \ + break; \ + } + + switch (variance->type.type) { + KERNEL_VARIANCE_CASE(U8); + KERNEL_VARIANCE_CASE(U16); + KERNEL_VARIANCE_CASE(U32); + KERNEL_VARIANCE_CASE(U64); + KERNEL_VARIANCE_CASE(S8); + KERNEL_VARIANCE_CASE(S16); + KERNEL_VARIANCE_CASE(S32); + KERNEL_VARIANCE_CASE(S64); + KERNEL_VARIANCE_CASE(F32); + KERNEL_VARIANCE_CASE(F64); + default: + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", + image->type.type); + return false; + } + } + + return true; +} + + +// Generate Lanczos interpolation kernels +static void lanczos(double values[], // Interpolation kernel to generate + int num, // Number of values in the kernel + float frac // Sub-pixel position + ) +{ + // XXX: Instead of generating a convolution kernel that does no shifting, need to set a boolean that says + // we can do an exact shift. + if (fabs(frac) < DBL_EPSILON) { + // No real shift + for (int i = 0; i < (num - 1) / 2; i++) { + values[i] = 0.0; + } + values[(num - 1) / 2] = 1.0; + for (int i = (num - 1) / 2 + 1; i < num; i++) { + values[i] = 0.0; + } + } else { + double norm1 = 2.0 / PS_SQR(M_PI); // Normalisation for laczos + double norm2 = 2.0 / (num / 2); // Normalisation for sinc functions + double pos = frac - (num - 1)/2; // Position of interest + for (int i = 0; i < num; i++, pos += 1.0) { + values[i] = norm1 * sin(M_PI * pos * norm2) * sin(M_PI_2 * pos * norm2) / PS_SQR(pos); + } + } + + return; +} + +// Interpolation engine for separable interpolation kernels (either for good reasons or for practical reasons) +static bool interpolateSeparate(double *imageValue, double *varianceValue, psMaskType *maskValue, + float x, float y, const psImageInterpolateOptions *options) +{ + // Parameters have been checked by psImageInterpolate() + + int xNum, yNum; // Number of interpolation kernel pixels + switch (options->mode) { + case PS_INTERPOLATE_LANCZOS2: + xNum = yNum = 4; + break; + case PS_INTERPOLATE_LANCZOS3: + xNum = yNum = 6; + break; + case PS_INTERPOLATE_LANCZOS4: + xNum = yNum = 8; + break; + case PS_INTERPOLATE_FLAT: + case PS_INTERPOLATE_BILINEAR: + case PS_INTERPOLATE_BICUBE: + case PS_INTERPOLATE_GAUSS: + default: + psAbort("Invalid interpolation mode."); + } + + const psImage *image = options->image; // Image of interest + int xFloor = floor(x - 0.5 + FLT_EPSILON); // Pixel below point of interest in x + int yFloor = floor(y - 0.5 + FLT_EPSILON); // Pixel below point of interest in y + int xLast = image->numCols - 1; // Last pixel in x + int yLast = image->numRows - 1; // Last pixel in y + + if (xFloor - (xNum - 1) / 2 < 0 || xFloor + xNum / 2 > xLast || + yFloor - (yNum - 1) / 2 < 0 || yFloor + yNum / 2 > yLast) { + // At least one pixel of the interpolation kernel is off the image + if (imageValue) { + *imageValue = options->badImage; + } + if (varianceValue) { + *varianceValue = options->badVariance; + } + if (maskValue) { + *maskValue = options->badMask; + } + return true; + } + +// bool xExact, yExact; // Is the shift exactly on? + double xKernel[xNum], yKernel[yNum];// Interpolation kernels in x and y + switch (options->mode) { + case PS_INTERPOLATE_LANCZOS2: + case PS_INTERPOLATE_LANCZOS3: + case PS_INTERPOLATE_LANCZOS4: { + double xFrac = x - 0.5 - xFloor; // Fraction of pixel in x +#if 0 + if (fabs(xFrac) < DBL_EPSILON) { + xExact = true; + } else { +#endif + lanczos(xKernel, xNum, xFrac); +#if 0 + xExact = false; + } +#endif + double yFrac = y - 0.5 - yFloor; // Fraction of pixel in y +#if 0 + if (fabs(yFrac) < DBL_EPSILON) { + yExact = true; + } else { +#endif + lanczos(yKernel, yNum, yFrac); +#if 0 + yExact = false; + } +#endif + break; + } + case PS_INTERPOLATE_FLAT: + case PS_INTERPOLATE_BILINEAR: + case PS_INTERPOLATE_BICUBE: + case PS_INTERPOLATE_GAUSS: + default: + psAbort("Invalid interpolation mode."); + } + + // Image interpolation, according to image type + #define SEPARATE_IMAGE_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) { \ + double xInterpValue = 0.0; /* Interpolation in x */ \ + for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) { \ + xInterpValue += values[j][i] = xKernel[i] * image->data.TYPE[yPix][xPix]; \ + } \ + value += yKernel[j] * xInterpValue; /* Interpolating in y */ \ + } \ + break; \ + } + + // Calculate the value for the image + double value = 0.0; // Value to return + double values[yNum][xNum]; // Values of interest + switch (image->type.type) { + SEPARATE_IMAGE_CASE(U8); + SEPARATE_IMAGE_CASE(U16); + SEPARATE_IMAGE_CASE(U32); + SEPARATE_IMAGE_CASE(U64); + SEPARATE_IMAGE_CASE(S8); + SEPARATE_IMAGE_CASE(S16); + SEPARATE_IMAGE_CASE(S32); + SEPARATE_IMAGE_CASE(S64); + SEPARATE_IMAGE_CASE(F32); + SEPARATE_IMAGE_CASE(F64); + default: + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", + image->type.type); + return false; + } + if (imageValue) { + *imageValue = value; + } + + // Check the mask value + if (maskValue) { + const psImage *mask = options->mask; // Image mask + psMaskType maskVal = options->maskVal; // Mask value + double badValue = 0.0; // Amount of flux in bad pixels + int badPix = 0; // Number of bad pixels + *maskValue = 0; + for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) { + // Interpolation in x + double xInterpValue = 0.0; + for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) { + if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { + xInterpValue += values[j][i]; + badPix++; + } + *maskValue |= mask->data.PS_TYPE_MASK_DATA[yPix][xPix]; + } + // Interpolating in y + badValue += yKernel[j] * xInterpValue; + } + + if (badPix > 0) { + if (fabs(badValue) >= fabs(value) * options->poorFrac) { + *maskValue |= options->badMask; + } else { + *maskValue |= options->poorMask; + } + } + } + + // Finally, the variance + if (varianceValue) { + const psImage *variance = options->variance; // Image variance + *varianceValue = 0.0; + + // Variance interpolation, according to image type + #define SEPARATE_VARIANCE_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + for (int j = 0, yPix = yFloor - (yNum - 1) / 2; j < yNum; j++, yPix++) { \ + double xInterpValue = 0.0; /* Interpolation in x */ \ + for (int i = 0, xPix = xFloor - (xNum - 1) / 2; i < xNum; i++, xPix++) { \ + xInterpValue += PS_SQR(xKernel[i]) * variance->data.TYPE[yPix][xPix]; \ + } \ + *varianceValue += xInterpValue * PS_SQR(yKernel[j]); /* Interpolating in y */ \ + } \ + break; \ + } + + switch (variance->type.type) { + SEPARATE_VARIANCE_CASE(U8); + SEPARATE_VARIANCE_CASE(U16); + SEPARATE_VARIANCE_CASE(U32); + SEPARATE_VARIANCE_CASE(U64); + SEPARATE_VARIANCE_CASE(S8); + SEPARATE_VARIANCE_CASE(S16); + SEPARATE_VARIANCE_CASE(S32); + SEPARATE_VARIANCE_CASE(S64); + SEPARATE_VARIANCE_CASE(F32); + SEPARATE_VARIANCE_CASE(F64); + default: + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", + variance->type.type); + return false; + } + } + + return true; +} + + + +bool psImageInterpolate(double *imageValue, double *varianceValue, psMaskType *maskValue, + float x, float y, const psImageInterpolateOptions *options) +{ + PS_ASSERT_PTR_NON_NULL(options, false); + + const psImage *image = options->image; // Image to interpolate + const psImage *mask = options->mask; // Mask to interpolate + const psImage *variance = options->variance; // Variance to interpolate + + PS_ASSERT_IMAGE_NON_NULL(image, false); + if (varianceValue) { + PS_ASSERT_IMAGE_NON_NULL(variance, false); + PS_ASSERT_IMAGE_TYPE(variance, image->type.type, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(variance, image, false); + } + if (maskValue) { + PS_ASSERT_IMAGE_NON_NULL(mask, false); + PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(mask, image, false); + } + + PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(options->poorFrac, 0.0, false); + + switch (options->mode) { + case PS_INTERPOLATE_FLAT: + return interpolateFlat(imageValue, varianceValue, maskValue, x, y, options); + case PS_INTERPOLATE_BILINEAR: + case PS_INTERPOLATE_BICUBE: + case PS_INTERPOLATE_GAUSS: + return interpolateKernel(imageValue, varianceValue, maskValue, x, y, options); + case PS_INTERPOLATE_LANCZOS2: + case PS_INTERPOLATE_LANCZOS3: + case PS_INTERPOLATE_LANCZOS4: + return interpolateSeparate(imageValue, varianceValue, maskValue, x, y, options); + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + _("Specified interpolation method (%d) is not supported."), + options->mode); + return false; + } + + psAbort("Should never reach here."); + return false; +} + +psImageInterpolateMode psImageInterpolateModeFromString(const char *name) +{ + PS_ASSERT_STRING_NON_EMPTY(name, PS_INTERPOLATE_NONE); + + if (!strcasecmp(name, "FLAT")) return PS_INTERPOLATE_FLAT; + if (!strcasecmp(name, "BILINEAR")) return PS_INTERPOLATE_BILINEAR; + if (!strcasecmp(name, "BICUBE")) return PS_INTERPOLATE_BICUBE; + if (!strcasecmp(name, "GAUSS")) return PS_INTERPOLATE_GAUSS; + if (!strcasecmp(name, "LANCZOS2")) return PS_INTERPOLATE_LANCZOS2; + if (!strcasecmp(name, "LANCZOS3")) return PS_INTERPOLATE_LANCZOS3; + if (!strcasecmp(name, "LANCZOS4")) return PS_INTERPOLATE_LANCZOS4; + + psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Unknown interpolation type %s"), name); + return PS_INTERPOLATE_NONE; +} Index: trunk/psLib/src/imageops/psImageInterpolate.h =================================================================== --- trunk/psLib/src/imageops/psImageInterpolate.h (revision 12741) +++ trunk/psLib/src/imageops/psImageInterpolate.h (revision 12741) @@ -0,0 +1,79 @@ +/* @file psImageInterpolate.h + * @brief Functions for interpolating an image + * + * @author Robert DeSonia, MHPCC + * @author Ross Harman, MHPCC + * @author Joshua Hoblitt, University of Hawaii + * @author Paul Price, Institute for Astronomy + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-04-04 22:42:02 $ + * Copyright 2004-2007 Institute for Astronomy, University of Hawaii + */ + +#ifndef PS_IMAGE_PIXEL_INTERPOLATE_H +#define PS_IMAGE_PIXEL_INTERPOLATE_H + + +/// Enumeration of options in interpolation +typedef enum { + PS_INTERPOLATE_NONE, ///< no interpolate defined (error state) + PS_INTERPOLATE_FLAT, ///< Flat interpolation (nearest pixel) + PS_INTERPOLATE_BILINEAR, ///< Bilinear interpolation + PS_INTERPOLATE_BICUBE, ///< Bicubic interpolation with 3x3 region + PS_INTERPOLATE_GAUSS, ///< Gaussian inteprolation with 3x3 region + PS_INTERPOLATE_LANCZOS2, ///< Sinc interpolation with 4x4 pixel kernel + PS_INTERPOLATE_LANCZOS3, ///< Sinc interpolation with 6x6 pixel kernel + PS_INTERPOLATE_LANCZOS4, ///< Sinc interpolation with 8x8 pixel kernel +} psImageInterpolateMode; + + +/// Options for general interpolation. +/// +/// We stuff in here all the constant values when doing interpolation, so that not all of it has to be pushed +/// onto the stack in the middle of a tight loop. For this reason, even the image, mask and variance map are +/// included. +typedef struct { + psImageInterpolateMode mode; // Interpolation mode + const psImage *image; // Input image for interpolation + const psImage *variance; // Variance image for interpolation + const psImage *mask; // Mask image for interpolation + psMaskType maskVal; // Value to mask + double badImage; // Image value if x,y location is not good + double badVariance; // Variance value if x,y location is not good + psMaskType badMask; // Mask value to give bad pixels + psMaskType poorMask; // Mask value to give poor pixels + float poorFrac; // Fraction of flux in bad pixels before output is marked bad +} psImageInterpolateOptions; + + +/// Allocator +psImageInterpolateOptions *psImageInterpolateOptionsAlloc(psImageInterpolateMode mode, // Interpolation mode + const psImage *image, // Input image + const psImage *variance, // Variance image + const psImage *mask, // Mask image + psMaskType maskVal, // Value to mask + double badImage, // Value for image if bad + double badVariance, // Value for variance if bad + psMaskType badMask, // Mask value for bad pixels + psMaskType poorMask, // Mask value for poor pixels + float poorFrac // Fraction of flux for question + ); + + + +/// Interpolate image pixel value given floating point coordinates. +bool psImageInterpolate(double *imageValue, ///< Return value for image + double *varianceValue, ///< Return value for variance + psMaskType *maskValue, ///< Return value for mask + float x, float y, ///< Location to which to interpolate + const psImageInterpolateOptions *options ///< Options for interpolation + ); + + +// Return the appropriate interpolation mode given a char string name for that mode +psImageInterpolateMode psImageInterpolateModeFromString(const char *name // Mode name + ); + + +#endif Index: trunk/psLib/src/imageops/psImagePixelExtract.c =================================================================== --- trunk/psLib/src/imageops/psImagePixelExtract.c (revision 12738) +++ trunk/psLib/src/imageops/psImagePixelExtract.c (revision 12741) @@ -8,6 +8,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.31 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-03-14 00:39:50 $ + * @version $Revision: 1.32 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-04-04 22:42:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -23,8 +23,8 @@ #include "psMemory.h" #include "psVector.h" +#include "psError.h" +#include "psImage.h" +#include "psImageInterpolate.h" #include "psImagePixelExtract.h" -#include "psError.h" - - #define VECTOR_STORE_ROW_CASE(TYPE) \ @@ -679,4 +679,7 @@ float dY = (endRow - startRow) / (float)(nSamples-1); + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(mode, input, NULL, mask, maskVal, + 0, 0, 0, 0, 0); + #define LINEAR_CUT_CASE(TYPE) \ case PS_TYPE_##TYPE: { \ @@ -692,5 +695,12 @@ cutRowsData[i] = y; \ } \ - outData[i] = psImagePixelInterpolate(input,x,y,mask,maskVal,0,mode); \ + double value; \ + if (!psImageInterpolate(&value, NULL, NULL, x, y, interp)) { \ + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); \ + psFree(interp); \ + psFree(out); \ + return NULL; \ + } \ + outData[i] = value; \ } \ } \ @@ -709,15 +719,17 @@ LINEAR_CUT_CASE(F32); LINEAR_CUT_CASE(F64); - - default: { - char* typeStr; - PS_TYPE_NAME(typeStr,input->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psImage type, %s, is not supported."), - typeStr); - psFree(out); - out = NULL; - } - } + default: { + char* typeStr; + PS_TYPE_NAME(typeStr,input->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + _("Specified psImage type, %s, is not supported."), + typeStr); + psFree(interp); + psFree(out); + out = NULL; + } + } + + psFree(interp); return out; Index: trunk/psLib/src/mathtypes/psImage.c =================================================================== --- trunk/psLib/src/mathtypes/psImage.c (revision 12738) +++ trunk/psLib/src/mathtypes/psImage.c (revision 12741) @@ -9,6 +9,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.127 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-03-22 00:11:08 $ + * @version $Revision: 1.128 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-04-04 22:42:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -572,82 +572,4 @@ } -double psImagePixelInterpolate(const psImage* input, - float x, - float y, - const psImage* mask, - psMaskType maskVal, - double unexposedValue, - psImageInterpolateMode mode) -{ - PS_ASSERT_IMAGE_NON_NULL(input, unexposedValue); - - #define PSIMAGE_PIXEL_INTERPOLATE_CASE(TYPE) \ -case PS_TYPE_##TYPE: \ - switch (mode) { \ - case PS_INTERPOLATE_FLAT: \ - return p_psImagePixelInterpolateFLAT_##TYPE( \ - input, \ - x, \ - y, \ - mask, \ - maskVal, \ - unexposedValue); \ - break; \ - case PS_INTERPOLATE_BILINEAR: \ - return p_psImagePixelInterpolateBILINEAR_##TYPE( \ - input, \ - x, \ - y, \ - mask, \ - maskVal, \ - unexposedValue); \ - break; \ - case PS_INTERPOLATE_BILINEAR_VARIANCE: \ - return p_psImagePixelInterpolateBILINEAR_VARIANCE_##TYPE( \ - input, \ - x, \ - y, \ - mask, \ - maskVal, \ - unexposedValue); \ - break; \ - case PS_INTERPOLATE_BICUBE: \ - return p_psImagePixelInterpolateBICUBE_##TYPE( \ - input, \ - x, \ - y, \ - mask, \ - maskVal, \ - unexposedValue); \ - break; \ - default: \ - psError(PS_ERR_BAD_PARAMETER_VALUE,true, \ - _("Specified interpolation method (%d) is not supported."), \ - mode); \ - } \ - break; - - switch (input->type.type) { - PSIMAGE_PIXEL_INTERPOLATE_CASE(U8); - PSIMAGE_PIXEL_INTERPOLATE_CASE(U16); - PSIMAGE_PIXEL_INTERPOLATE_CASE(U32); - PSIMAGE_PIXEL_INTERPOLATE_CASE(U64); - PSIMAGE_PIXEL_INTERPOLATE_CASE(S8); - PSIMAGE_PIXEL_INTERPOLATE_CASE(S16); - PSIMAGE_PIXEL_INTERPOLATE_CASE(S32); - PSIMAGE_PIXEL_INTERPOLATE_CASE(S64); - PSIMAGE_PIXEL_INTERPOLATE_CASE(F32); - PSIMAGE_PIXEL_INTERPOLATE_CASE(F64); - default: { - char* typeStr; - PS_TYPE_NAME(typeStr,input->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE,true, - _("Specified psImage type, %s, is not supported."), - typeStr); - } - } - - return unexposedValue; -} psF64 p_psImageGetElementF64(psImage* image, @@ -699,249 +621,2 @@ } } - -#define PSIMAGE_PIXEL_INTERPOLATE_FLAT(TYPE,RETURNTYPE) \ -inline RETURNTYPE p_psImagePixelInterpolateFLAT_##TYPE( \ - const psImage* input, \ - float x, \ - float y, \ - const psImage* mask, \ - psMaskType maskVal, \ - RETURNTYPE unexposedValue) \ -{ \ - psS32 intX = (psS32) round((psF64)(x) - 0.5 + FLT_EPSILON); \ - psS32 intY = (psS32) round((psF64)(y) - 0.5 + FLT_EPSILON); \ - psS32 lastX = input->numCols - 1; \ - psS32 lastY = input->numRows - 1; \ - \ - if ((intX < 0) || \ - (intX > lastX) || \ - (intY < 0) || \ - (intY > lastY) || \ - ( (mask!=NULL) && \ - ((mask->data.PS_TYPE_MASK_DATA[intY][intX] & maskVal) != 0) ) ) { \ - return unexposedValue; \ - } \ - \ - return input->data.TYPE[intY][intX]; \ -} - -PSIMAGE_PIXEL_INTERPOLATE_FLAT(U8,psF64) -PSIMAGE_PIXEL_INTERPOLATE_FLAT(U16,psF64) -PSIMAGE_PIXEL_INTERPOLATE_FLAT(U32,psF64) -PSIMAGE_PIXEL_INTERPOLATE_FLAT(U64,psF64) -PSIMAGE_PIXEL_INTERPOLATE_FLAT(S8,psF64) -PSIMAGE_PIXEL_INTERPOLATE_FLAT(S16,psF64) -PSIMAGE_PIXEL_INTERPOLATE_FLAT(S32,psF64) -PSIMAGE_PIXEL_INTERPOLATE_FLAT(S64,psF64) -PSIMAGE_PIXEL_INTERPOLATE_FLAT(F32,psF64) -PSIMAGE_PIXEL_INTERPOLATE_FLAT(F64,psF64) - -#define PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(TYPE, RETURNTYPE, SUFFIX, FRACFUNC) \ -inline RETURNTYPE p_psImagePixelInterpolateBILINEAR_##SUFFIX( \ - const psImage* input, \ - float x, \ - float y, \ - const psImage* mask, \ - psMaskType maskVal, \ - RETURNTYPE unexposedValue) \ -{ \ - int floorX = floor((x) - 0.5); \ - int floorY = floor((y) - 0.5); \ - float fracX = x - 0.5 - floorX; \ - float fracY = y - 0.5 - floorY; \ - int lastX = input->numCols - 1; \ - int lastY = input->numRows - 1; \ - ps##TYPE V00 = 0; \ - ps##TYPE V01 = 0; \ - ps##TYPE V10 = 0; \ - ps##TYPE V11 = 0; \ - bool valid00; \ - bool valid01; \ - bool valid10; \ - bool valid11; \ - \ - if (floorY >= 0 && floorY <= lastY) { \ - if (floorX >= 0 && floorX <= lastX) { \ - V00 = input->data.TYPE[floorY][floorX]; \ - valid00 = (mask == NULL) || \ - ((mask->data.PS_TYPE_MASK_DATA[floorY][floorX] & maskVal) == 0); \ - } else { \ - valid00 = false; \ - } \ - if (floorX >= -1 && floorX < lastX) { \ - V10 = input->data.TYPE[floorY][floorX+1]; \ - valid10 = (mask == NULL) || \ - ((mask->data.PS_TYPE_MASK_DATA[floorY][floorX+1] & maskVal) == 0); \ - } else { \ - valid10 = false; \ - } \ - } else { \ - valid00 = false; \ - valid10 = false; \ - } \ - if (floorY >= -1 && floorY < lastY) { \ - if (floorX >= 0 && floorX <= lastX) { \ - V01 = input->data.TYPE[floorY+1][floorX]; \ - valid01 = (mask == NULL) || \ - ((mask->data.PS_TYPE_MASK_DATA[floorY+1][floorX] & maskVal) == 0); \ - } else { \ - valid01 = false; \ - } \ - if (floorX >= -1 && floorX < lastX) { \ - V11 = input->data.TYPE[floorY+1][floorX+1]; \ - valid11 = (mask == NULL) || \ - ((mask->data.PS_TYPE_MASK_DATA[floorY+1][floorX+1] & maskVal) == 0); \ - } else { \ - valid11 = false; \ - } \ - } else { \ - valid01 = false; \ - valid11 = false; \ - } \ - \ - /* cover likely case of all pixels being valid more efficiently */ \ - if (valid00 && valid10 && valid01 && valid11) { \ - /* formula from the ADD */ \ - return V00*FRACFUNC((1.0-fracX)*(1.0-fracY)) + V10*FRACFUNC(fracX*(1.0-fracY)) + \ - V01*FRACFUNC(fracY*(1.0-fracX)) + V11*FRACFUNC(fracX*fracY); \ - } \ - \ - /* OK, at least one pixel is not valid - need to do it piecemeal */ \ - \ - RETURNTYPE V0 = 0.0; \ - bool valid0 = true; \ - if (valid00 && valid10) { \ - V0 = V00*FRACFUNC(1-fracX)+V10*FRACFUNC(fracX); \ - } else if (valid00) { \ - V0 = V00; \ - } else if (valid10) { \ - V0 = V10; \ - } else { \ - valid0 = false; \ - } \ - \ - RETURNTYPE V1 = 0.0; \ - bool valid1 = true; \ - if (valid01 && valid11) { \ - V1 = V01*FRACFUNC(1-fracX)+V11*FRACFUNC(fracX); \ - } else if (valid01) { \ - V1 = V01; \ - } else if (valid11) { \ - V1 = V11; \ - } else { \ - valid1 = false; \ - } \ - \ - if (valid0 && valid1) { \ - return V0*FRACFUNC(1-fracY) + V1*FRACFUNC(fracY); \ - } else if (valid0) { \ - return V0; \ - } else if (valid1) { \ - return V1; \ - } \ - \ - return unexposedValue; \ -} - -// XXX this would be much faster to use a 3x3 kernel to determine the shift -// the same function can be used for all equivalent (kernel-based) shifts... -#define PSIMAGE_PIXEL_INTERPOLATE_BICUBE(TYPE, RETURNTYPE, SUFFIX, FRACFUNC) \ -inline RETURNTYPE p_psImagePixelInterpolateBICUBE_##SUFFIX( \ - const psImage* input, \ - float x, \ - float y, \ - const psImage* mask, \ - psMaskType maskVal, \ - RETURNTYPE unexposedValue) \ -{ \ - int floorX = floor(x); \ - int floorY = floor(y); \ - psF64 fracX = x - floorX - 0.5; \ - psF64 fracY = y - floorY - 0.5; \ - psS32 lastX = input->numCols - 1; \ - psS32 lastY = input->numRows - 1; \ - if (floorX < 1) return unexposedValue; \ - if (floorY < 1) return unexposedValue; \ - if (floorX >= lastX) return unexposedValue; \ - if (floorY >= lastY) return unexposedValue; \ - \ - /* XXX use care for masked and boundary pixels */ \ - psF64 Vmm = input->data.TYPE[floorY-1][floorX-1]; \ - psF64 Vom = input->data.TYPE[floorY-1][floorX+0]; \ - psF64 Vpm = input->data.TYPE[floorY-1][floorX+1]; \ - psF64 Vmo = input->data.TYPE[floorY+0][floorX-1]; \ - psF64 Voo = input->data.TYPE[floorY+0][floorX+0]; \ - psF64 Vpo = input->data.TYPE[floorY+0][floorX+1]; \ - psF64 Vmp = input->data.TYPE[floorY+1][floorX-1]; \ - psF64 Vop = input->data.TYPE[floorY+1][floorX+0]; \ - psF64 Vpp = input->data.TYPE[floorY+1][floorX+1]; \ - \ - psF64 Vxm = Vmm + Vmo + Vmp; \ - psF64 Vxp = Vpm + Vpo + Vpp; \ - psF64 Vym = Vmm + Vom + Vpm; \ - psF64 Vyp = Vmp + Vop + Vpp; \ - psF64 Vo = Vym + Vyp + Vmo + Voo + Vpo; \ - \ - psF64 Z_00 = Vo*(5.0/9.0) - (Vxp + Vxm)/3.0 - (Vyp + Vym)/3.0; \ - \ - psF64 Z_10 = (Vxp - Vxm)/6.0; \ - psF64 Z_01 = (Vyp - Vym)/6.0; \ - psF64 Z_20 = (Vxp + Vxm)/2.0 - Vo/3.0; \ - psF64 Z_02 = (Vyp + Vym)/2.0 - Vo/3.0; \ - psF64 Z_11 = (Vpp + Vmm - Vpm - Vmp)/4.0; \ - \ - psF64 value = Z_00 + Z_10*fracX + Z_01*fracY + Z_20*fracX*fracX + Z_11*fracX*fracY + Z_02*fracY*fracY; \ - return value; \ -} - -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U8,psF64,U8,) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U16,psF64,U16,) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U32,psF64,U32,) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U64,psF64,U64,) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S8,psF64,S8,) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S16,psF64,S16,) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S32,psF64,S32,) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S64,psF64,S64,) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(F32,psF64,F32,) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(F64,psF64,F64,) - -PSIMAGE_PIXEL_INTERPOLATE_BICUBE(U8,psF64,U8,) -PSIMAGE_PIXEL_INTERPOLATE_BICUBE(U16,psF64,U16,) -PSIMAGE_PIXEL_INTERPOLATE_BICUBE(U32,psF64,U32,) -PSIMAGE_PIXEL_INTERPOLATE_BICUBE(U64,psF64,U64,) -PSIMAGE_PIXEL_INTERPOLATE_BICUBE(S8,psF64,S8,) -PSIMAGE_PIXEL_INTERPOLATE_BICUBE(S16,psF64,S16,) -PSIMAGE_PIXEL_INTERPOLATE_BICUBE(S32,psF64,S32,) -PSIMAGE_PIXEL_INTERPOLATE_BICUBE(S64,psF64,S64,) -PSIMAGE_PIXEL_INTERPOLATE_BICUBE(F32,psF64,F32,) -PSIMAGE_PIXEL_INTERPOLATE_BICUBE(F64,psF64,F64,) - -// Variance Version -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U8,psF64,VARIANCE_U8,PS_SQR) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U16,psF64,VARIANCE_U16,PS_SQR) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U32,psF64,VARIANCE_U32,PS_SQR) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U64,psF64,VARIANCE_U64,PS_SQR) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S8,psF64,VARIANCE_S8,PS_SQR) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S16,psF64,VARIANCE_S16,PS_SQR) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S32,psF64,VARIANCE_S32,PS_SQR) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S64,psF64,VARIANCE_S64,PS_SQR) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(F32,psF64,VARIANCE_F32,PS_SQR) -PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(F64,psF64,VARIANCE_F64,PS_SQR) - -psImageInterpolateMode psImageInterpolateModeFromString (char *name) { - - if (!strcasecmp(name, "FLAT")) return PS_INTERPOLATE_FLAT; - if (!strcasecmp(name, "BILINEAR")) return PS_INTERPOLATE_BILINEAR; - if (!strcasecmp(name, "BICUBE")) return PS_INTERPOLATE_BICUBE; - if (!strcasecmp(name, "GAUSS")) return PS_INTERPOLATE_GAUSS; - if (!strcasecmp(name, "LANCZOS2")) return PS_INTERPOLATE_LANCZOS2; - if (!strcasecmp(name, "LANCZOS3")) return PS_INTERPOLATE_LANCZOS3; - if (!strcasecmp(name, "LANCZOS4")) return PS_INTERPOLATE_LANCZOS4; - if (!strcasecmp(name, "BILINEAR_VARIANCE")) return PS_INTERPOLATE_BILINEAR_VARIANCE; - if (!strcasecmp(name, "LANCZOS2_VARIANCE")) return PS_INTERPOLATE_LANCZOS2_VARIANCE; - if (!strcasecmp(name, "LANCZOS3_VARIANCE")) return PS_INTERPOLATE_LANCZOS3_VARIANCE; - if (!strcasecmp(name, "LANCZOS4_VARIANCE")) return PS_INTERPOLATE_LANCZOS4_VARIANCE; - - psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Unknown interpolate type %s"), name); - return PS_INTERPOLATE_NONE; -} Index: trunk/psLib/src/mathtypes/psImage.h =================================================================== --- trunk/psLib/src/mathtypes/psImage.h (revision 12738) +++ trunk/psLib/src/mathtypes/psImage.h (revision 12741) @@ -9,6 +9,6 @@ * @author Joshua Hoblitt, University of Hawaii * - * @version $Revision: 1.91 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-03-22 00:11:08 $ + * @version $Revision: 1.92 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-04-04 22:42:02 $ * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii */ @@ -24,23 +24,4 @@ #include "psArray.h" #include "psConstants.h" - -/** enumeration of options in interpolation - * - */ -typedef enum { - PS_INTERPOLATE_NONE, ///< no interpolate defined (error state) - PS_INTERPOLATE_FLAT, ///< 'flat' interpolation (nearest pixel) - PS_INTERPOLATE_BILINEAR, ///< bi-linear interpolation - PS_INTERPOLATE_BICUBE, ///< bi-cubic interpolation with 3x3 region (EAM) - PS_INTERPOLATE_GAUSS, ///< bi-cubic interpolation with 3x3 region (EAM) - PS_INTERPOLATE_LANCZOS2, ///< Sinc interpolation with 4x4 pixel kernel - PS_INTERPOLATE_LANCZOS3, ///< Sinc interpolation with 6x6 pixel kernel - PS_INTERPOLATE_LANCZOS4, ///< Sinc interpolation with 8x8 pixel kernel - PS_INTERPOLATE_BILINEAR_VARIANCE, ///< Variance version of PS_INTERPOLATE_BILINEAR - PS_INTERPOLATE_LANCZOS2_VARIANCE, ///< Variance version of PS_INTERPOLATE_LANCZOS2 - PS_INTERPOLATE_LANCZOS3_VARIANCE, ///< Variance version of PS_INTERPOLATE_LANCZOS3 - PS_INTERPOLATE_LANCZOS4_VARIANCE ///< Variance version of PS_INTERPOLATE_LANCZOS4 - // PS_INTERPOLATE_NUM_MODES ///< enum end-marker; does not coorespond to a interpolation mode -} psImageInterpolateMode; @@ -226,54 +207,5 @@ -/** Interpolate image pixel value given floating point coordinates. - * - * @return double Pixel value interpolated from image or unexposedValue if - * given x,y doesn't coorespond to a valid image location - */ -double psImagePixelInterpolate( - const psImage* input, ///< input image for interpolation - float x, ///< column location to derive value of - float y, ///< row location ot derive value of - const psImage* mask, ///< if not NULL, the mask of the input image - psMaskType maskVal, ///< the mask value - double unexposedValue, ///< return value if x,y location is not in image. - psImageInterpolateMode mode ///< interpolation mode -); - -// return the mode equivalent to a char string name -psImageInterpolateMode psImageInterpolateModeFromString (char *name); - - -#define PIXEL_INTERPOLATE_FCN_PROTOTYPE(SUFFIX, RETURNTYPE) \ -inline RETURNTYPE p_psImagePixelInterpolate##SUFFIX( \ - const psImage* input, /**< input image for interpolation */ \ - float x, /**< column location to derive value of */ \ - float y, /**< row location ot derive value of */ \ - const psImage* mask, /**< if not NULL, the mask of the input image */ \ - psMaskType maskVal, /**< the mask value */ \ - RETURNTYPE unexposedValue /**< return value if x,y location is not in image. */ \ - ); - -#define PIXEL_INTERPOLATE_FCNS(MODE) \ -PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_U8,psF64) \ -PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_U16,psF64) \ -PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_U32,psF64) \ -PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_U64,psF64) \ -PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_S8,psF64) \ -PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_S16,psF64) \ -PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_S32,psF64) \ -PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_S64,psF64) \ -PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_F32,psF64) \ -PIXEL_INTERPOLATE_FCN_PROTOTYPE(MODE##_F64,psF64) - -#ifndef SWIG -PIXEL_INTERPOLATE_FCNS(FLAT) -PIXEL_INTERPOLATE_FCNS(BILINEAR) -PIXEL_INTERPOLATE_FCNS(BILINEAR_VARIANCE) -PIXEL_INTERPOLATE_FCNS(BICUBE) -#endif // ! SWIG - -#undef PIXEL_INTERPOLATE_FCN_PROTOTYPE -#undef PIXEL_INTERPOLATE_FCNS + /***************************************************************************** Index: trunk/psLib/src/pslib_strict.h =================================================================== --- trunk/psLib/src/pslib_strict.h (revision 12738) +++ trunk/psLib/src/pslib_strict.h (revision 12741) @@ -9,6 +9,6 @@ * @author Eric Van Alst, MHPCC * -* @version $Revision: 1.28 $ $Name: not supported by cvs2svn $ -* @date $Date: 2007-03-27 02:43:22 $ +* @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ +* @date $Date: 2007-04-04 22:42:02 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -44,4 +44,5 @@ #include "psRegion.h" +#include "psImageInterpolate.h" #include "psImageConvolve.h" #include "psImageGeomManip.h" Index: trunk/psLib/test/mathtypes/tap_psImageInterpolate.c =================================================================== --- trunk/psLib/test/mathtypes/tap_psImageInterpolate.c (revision 12738) +++ trunk/psLib/test/mathtypes/tap_psImageInterpolate.c (revision 12741) @@ -8,5 +8,5 @@ int main (void) { - plan_tests(47); + plan_tests(88); // diag("psImageInterpolate() tests"); @@ -19,25 +19,34 @@ // generate simple image (x ramp) - psImage *image = psImageAlloc(32, 32, PS_TYPE_F32); - ok(image != NULL, "psImage successfully allocated"); - skip_start(image == NULL, 5, "Skipping tests because psImageAlloc() failed"); - - image->data.F32[10][10] = 1; - - // center of pixels is 0.5, 0.5 - float value; - - value = psImagePixelInterpolate (image, 10.5, 10.5, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float (value, 1.0, "pixel center value - %f", value); + psImage *image = psImageAlloc(32, 32, PS_TYPE_F64); + ok(image != NULL, "psImage successfully allocated"); + skip_start(image == NULL, 5, "Skipping tests because psImageAlloc() failed"); + + psImageInit(image, 0.0); + image->data.F64[10][10] = 1; + + // center of pixels is 0.5, 0.5 + double value; + + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, + image, NULL, NULL, 0, 0.0, 0.0, + 0, 0, 0.0); + ok(interp, "Interpolation options set"); + + ok(psImageInterpolate(&value, NULL, NULL, 10.5, 10.5, interp), "Interpolation"); + is_double (value, 1.0, "pixel center value - %f", value); // diag ("why do I need to have tolerances of 4epsilon or so??"); - value = psImagePixelInterpolate (image, 10.9, 10.5, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float_tol (value, 0.6, 4.0*FLT_EPSILON, "pixel value - %.20f", value); - - value = psImagePixelInterpolate (image, 10.5, 10.9, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float_tol (value, 0.6, 4.0*FLT_EPSILON, "pixel value - %.20f", value); - - value = psImagePixelInterpolate (image, 10.1, 10.5, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float_tol (value, 0.6, 4.0*FLT_EPSILON, "pixel value - %.20f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 10.9, 10.5, interp), "Interpolation"); + is_double_tol (value, 0.6, 4.0*FLT_EPSILON, "pixel value - %.20f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 10.5, 10.9, interp), "Interpolation"); + is_double_tol (value, 0.6, 4.0*FLT_EPSILON, "pixel value - %.20f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 10.1, 10.5, interp), "Interpolation"); + is_double_tol (value, 0.6, 4.0*FLT_EPSILON, "pixel value - %.20f", value); + + psFree(interp); skip_end(); @@ -54,5 +63,63 @@ // generate simple image (x ramp) - psImage *image = psImageAlloc(32, 32, PS_TYPE_F32); + psImage *image = psImageAlloc(32, 32, PS_TYPE_F64); + ok(image != NULL, "psImage successfully allocated"); + skip_start(image == NULL, 5, "Skipping tests because psImageAlloc() failed"); + + for (int j = 0; j < image->numRows; j++) { + for (int i = 0; i < image->numCols; i++) { + image->data.F64[j][i] = i + 0.5; + } + } + + // center of pixels is 0.5, 0.5 + double value; + + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, + image, NULL, NULL, 0, 0.0, 0.0, + 0, 0, 0.0); + ok(interp, "Interpolation options set"); + + ok(psImageInterpolate(&value, NULL, NULL, 2.5, 2.5, interp), "Interpolation"); + is_double_tol (value, 2.5, 5.0e-8, "pixel center value - %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.2, 2.5, interp), "Interpolation"); + is_double_tol (value, 2.2, 5.0e-8, "coord: 2.2, 2.5, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.8, 2.5, interp), "Interpolation"); + is_double_tol (value, 2.8, 5.0e-8, "coord: 2.8, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.8, 2.2, interp), "Interpolation"); + is_double_tol (value, 2.8, 5.0e-8, "coord: 2.8, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.8, 2.8, interp), "Interpolation"); + is_double_tol (value, 2.8, 5.0e-8, "coord: 2.8, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 0.8, 2.8, interp), "Interpolation"); + is_double_tol (value, 0.8, 5.0e-8, "coord: 0.8, value: %f", value); + + // no extrapolation + ok(psImageInterpolate(&value, NULL, NULL, 0.3, 2.8, interp), "Interpolation"); + is_double (value, 0.0, "coord: 0.3, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, -0.2, 2.8, interp), "Interpolation"); + is_double (value, 0.0, "coord: -0.2, value: %f", value); + + psFree(interp); + + skip_end(); + + psFree(image); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + // very simple tests: no mask, bilinear mode, yramp image only + { + psMemId id = psMemGetId(); + +// diag ("interpolate a y-ramp: "); + + // generate simple image (y ramp) + psImage *image = psImageAlloc(32, 32, PS_TYPE_F64); ok(image != NULL, "psImage successfully allocated"); skip_start(image == NULL, 5, "Skipping tests because psImageAlloc() failed"); @@ -61,50 +128,47 @@ { for (int i = 0; i < image->numCols; i++) { - image->data.F32[j][i] = i + 0.5; - } - } - - // center of pixels is 0.5, 0.5 - float value; - - value = psImagePixelInterpolate (image, 2.5, 2.5, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float (value, 2.5, "pixel center value - %f", value); - - value = psImagePixelInterpolate (image, 2.2, 2.5, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float (value, 2.2, "coord: 2.2, 2.5, value: %f", value); - - value = psImagePixelInterpolate (image, 2.8, 2.5, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float (value, 2.8, "coord: 2.8, value: %f", value); - - value = psImagePixelInterpolate (image, 2.8, 2.2, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float (value, 2.8, "coord: 2.8, value: %f", value); - - value = psImagePixelInterpolate (image, 2.8, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float (value, 2.8, "coord: 2.8, value: %f", value); - - value = psImagePixelInterpolate (image, 0.8, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float (value, 0.8, "coord: 0.8, value: %f", value); - - // no extrapolation - value = psImagePixelInterpolate (image, 0.3, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float (value, 0.5, "coord: 0.3, value: %f", value); - - value = psImagePixelInterpolate (image, -0.2, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float (value, 0.5, "coord: -0.2, value: %f", value); - - skip_end(); - - psFree(image); - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - // very simple tests: no mask, bilinear mode, yramp image only - { - psMemId id = psMemGetId(); - -// diag ("interpolate a y-ramp: "); - - // generate simple image (y ramp) - psImage *image = psImageAlloc(32, 32, PS_TYPE_F32); + image->data.F64[j][i] = j + 0.5; + } + } + + // center of pixels is 0.5, 0.5 + double value; + + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, + image, NULL, NULL, 0, 0.0, 0.0, + 0, 0, 0.0); + ok(interp, "Interpolation options set"); + + ok(psImageInterpolate(&value, NULL, NULL, 2.5, 2.5, interp), "Interpolation"); + is_double_tol (value, 2.5, 5.0e-8, "pixel center value - %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.2, 2.2, interp), "Interpolation"); + is_double_tol (value, 2.2, 5.0e-8, "coord: 2.2, 2.2, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.5, 2.8, interp), "Interpolation"); + is_double_tol (value, 2.8, 5.0e-8, "coord: 2.8, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.2, 2.8, interp), "Interpolation"); + is_double_tol (value, 2.8, 5.0e-8, "coord: 2.8, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.8, 2.8, interp), "Interpolation"); + is_double_tol (value, 2.8, 5.0e-8, "coord: 2.8, value: %f", value); + + psFree(interp); + + skip_end(); + + psFree(image); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + // very simple tests: no mask, bicube mode, xramp image only + { + psMemId id = psMemGetId(); + +// diag ("interpolate an x-ramp (bicube)"); + + // generate simple image (x ramp) + psImage *image = psImageAlloc(32, 32, PS_TYPE_F64); ok(image != NULL, "psImage successfully allocated"); skip_start(image == NULL, 5, "Skipping tests because psImageAlloc() failed"); @@ -113,40 +177,59 @@ { for (int i = 0; i < image->numCols; i++) { - image->data.F32[j][i] = j + 0.5; - } - } - - // center of pixels is 0.5, 0.5 - float value; - - value = psImagePixelInterpolate (image, 2.5, 2.5, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float (value, 2.5, "pixel center value - %f", value); - - value = psImagePixelInterpolate (image, 2.2, 2.2, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float (value, 2.2, "coord: 2.2, 2.5, value: %f", value); - - value = psImagePixelInterpolate (image, 2.5, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float (value, 2.8, "coord: 2.8, value: %f", value); - - value = psImagePixelInterpolate (image, 2.2, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float (value, 2.8, "coord: 2.8, value: %f", value); - - value = psImagePixelInterpolate (image, 2.8, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); - is_float (value, 2.8, "coord: 2.8, value: %f", value); - - skip_end(); - - psFree(image); - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - // very simple tests: no mask, bicube mode, xramp image only - { - psMemId id = psMemGetId(); - -// diag ("interpolate an x-ramp (bicube)"); - - // generate simple image (x ramp) - psImage *image = psImageAlloc(32, 32, PS_TYPE_F32); + image->data.F64[j][i] = i + 0.5; + } + } + + // center of pixels is 0.5, 0.5 + double value; + + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BICUBE, + image, NULL, NULL, 0, 0.0, 0.0, + 0, 0, 0.0); + ok(interp, "Interpolation options set"); + + ok(psImageInterpolate(&value, NULL, NULL, 2.5, 2.5, interp), "Interpolation"); + is_double_tol (value, 2.5, 5.0e-8, "coord; 2.5, 2.5, value - %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.2, 2.5, interp), "Interpolation"); + is_double_tol (value, 2.2, 5.0e-8, "coord: 2.2, 2.5, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.8, 2.5, interp), "Interpolation"); + is_double_tol (value, 2.8, 5.0e-8, "coord: 2.8, 2.5, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.8, 2.2, interp), "Interpolation"); + is_double_tol (value, 2.8, 5.0e-8, "coord: 2.8, 2.2, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.8, 2.8, interp), "Interpolation"); + is_double_tol (value, 2.8, 5.0e-8, "coord: 2.8, 2.8, value: %f", value); + +// diag ("coords outside of nominal range (1 < x < Nx - 2) return 'uncover'"); + + // no extrapolation: these return the 'uncover' value + ok(psImageInterpolate(&value, NULL, NULL, 0.8, 2.8, interp), "Interpolation"); + is_double (value, 0.0, "coord: 0.8, 2.8, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 0.3, 2.8, interp), "Interpolation"); + is_double (value, 0.0, "coord: 0.3, 2.8, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, -0.2, 2.8, interp), "Interpolation"); + is_double (value, 0.0, "coord: -0.2, 2.8, value: %f", value); + + psFree(interp); + + skip_end(); + + psFree(image); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + // very simple tests: no mask, bilinear mode, yramp image only + { + psMemId id = psMemGetId(); + +// diag ("interpolate a y-ramp (bicube)"); + + // generate simple image (y ramp) + psImage *image = psImageAlloc(32, 32, PS_TYPE_F64); ok(image != NULL, "psImage successfully allocated"); skip_start(image == NULL, 5, "Skipping tests because psImageAlloc() failed"); @@ -155,52 +238,47 @@ { for (int i = 0; i < image->numCols; i++) { - image->data.F32[j][i] = i + 0.5; - } - } - - // center of pixels is 0.5, 0.5 - float value; - - value = psImagePixelInterpolate (image, 2.5, 2.5, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 2.5, "coord; 2.5, 2.5, value - %f", value); - - value = psImagePixelInterpolate (image, 2.2, 2.5, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 2.2, "coord: 2.2, 2.5, value: %f", value); - - value = psImagePixelInterpolate (image, 2.8, 2.5, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 2.8, "coord: 2.8, 2.5, value: %f", value); - - value = psImagePixelInterpolate (image, 2.8, 2.2, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 2.8, "coord: 2.8, 2.2, value: %f", value); - - value = psImagePixelInterpolate (image, 2.8, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 2.8, "coord: 2.8, 2.8, value: %f", value); - -// diag ("coords outside of nominal range (1 < x < Nx - 2) return 'uncover'"); - - // no extrapolation: these return the 'uncover' value - value = psImagePixelInterpolate (image, 0.8, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 0.0, "coord: 0.8, 2.8, value: %f", value); - - value = psImagePixelInterpolate (image, 0.3, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 0.0, "coord: 0.3, 2.8, value: %f", value); - - value = psImagePixelInterpolate (image, -0.2, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 0.0, "coord: -0.2, 2.8, value: %f", value); - - skip_end(); - - psFree(image); - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - // very simple tests: no mask, bilinear mode, yramp image only - { - psMemId id = psMemGetId(); - -// diag ("interpolate a y-ramp (bicube)"); - - // generate simple image (y ramp) - psImage *image = psImageAlloc(32, 32, PS_TYPE_F32); + image->data.F64[j][i] = j + 0.5; + } + } + + // center of pixels is 0.5, 0.5 + double value; + + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BICUBE, + image, NULL, NULL, 0, 0.0, 0.0, + 0, 0, 0.0); + ok(interp, "Interpolation options set"); + + ok(psImageInterpolate(&value, NULL, NULL, 2.5, 2.5, interp), "Interpolation"); + is_double_tol (value, 2.5, 5.0e-8, "pixel center value - %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.2, 2.2, interp), "Interpolation"); + is_double_tol (value, 2.2, 5.0e-8, "coord: 2.2, 2.5, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.5, 2.8, interp), "Interpolation"); + is_double_tol (value, 2.8, 5.0e-8, "coord: 2.8, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.2, 2.8, interp), "Interpolation"); + is_double_tol (value, 2.8, 5.0e-8, "coord: 2.8, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.8, 2.8, interp), "Interpolation"); + is_double_tol (value, 2.8, 5.0e-8, "coord: 2.8, value: %f", value); + + psFree(interp); + + skip_end(); + + psFree(image); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + // very simple tests: no mask, bilinear mode, x,y 2nd order shape + { + psMemId id = psMemGetId(); + +// diag ("interpolate a quadratic shape (bicube)"); + + // generate simple image (x ramp) + psImage *image = psImageAlloc(32, 32, PS_TYPE_F64); ok(image != NULL, "psImage successfully allocated"); skip_start(image == NULL, 5, "Skipping tests because psImageAlloc() failed"); @@ -209,67 +287,32 @@ { for (int i = 0; i < image->numCols; i++) { - image->data.F32[j][i] = j + 0.5; - } - } - - // center of pixels is 0.5, 0.5 - float value; - - value = psImagePixelInterpolate (image, 2.5, 2.5, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 2.5, "pixel center value - %f", value); - - value = psImagePixelInterpolate (image, 2.2, 2.2, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 2.2, "coord: 2.2, 2.5, value: %f", value); - - value = psImagePixelInterpolate (image, 2.5, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 2.8, "coord: 2.8, value: %f", value); - - value = psImagePixelInterpolate (image, 2.2, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 2.8, "coord: 2.8, value: %f", value); - - value = psImagePixelInterpolate (image, 2.8, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 2.8, "coord: 2.8, value: %f", value); - - skip_end(); - - psFree(image); - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - // very simple tests: no mask, bilinear mode, x,y 2nd order shape - { - psMemId id = psMemGetId(); - -// diag ("interpolate a quadratic shape (bicube)"); - - // generate simple image (x ramp) - psImage *image = psImageAlloc(32, 32, PS_TYPE_F32); - ok(image != NULL, "psImage successfully allocated"); - skip_start(image == NULL, 5, "Skipping tests because psImageAlloc() failed"); - - for (int j = 0; j < image->numRows; j++) - { - for (int i = 0; i < image->numCols; i++) { - image->data.F32[j][i] = 0.25*PS_SQR(i + 0.5) + j + 0.5; - } - } - - // center of pixels is 0.5, 0.5 - float value; - - value = psImagePixelInterpolate (image, 2.5, 2.5, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 4.0625, "pixel center value - %f", value); - - value = psImagePixelInterpolate (image, 2.2, 2.2, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 3.41, "coord: 2.2, 2.5, value: %f", value); - - value = psImagePixelInterpolate (image, 2.5, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 4.3625002, "coord: 2.5, 2.8, value: %f", value); - - value = psImagePixelInterpolate (image, 2.2, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 4.010000229, "coord: 2.2, 2.8, value: %f", value); - - value = psImagePixelInterpolate (image, 2.8, 2.8, NULL, 0, 0.0, PS_INTERPOLATE_BICUBE); - is_float (value, 4.75999975, "coord: 2.8, 2.8, value: %f", value); + image->data.F64[j][i] = 0.25*PS_SQR(i + 0.5) + j + 0.5; + } + } + + // center of pixels is 0.5, 0.5 + double value; + + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BICUBE, + image, NULL, NULL, 0, 0.0, 0.0, + 0, 0, 0.0); + ok(interp, "Interpolation options set"); + + ok(psImageInterpolate(&value, NULL, NULL, 2.5, 2.5, interp), "Interpolation"); + is_double_tol (value, 4.0625, 2.0e-7, "pixel center value - %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.2, 2.2, interp), "Interpolation"); + is_double_tol (value, 3.41, 2.0e-7, "coord: 2.2, 2.5, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.5, 2.8, interp), "Interpolation"); + is_double_tol (value, 4.3625, 2.0e-7, "coord: 2.5, 2.8, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.2, 2.8, interp), "Interpolation"); + is_double_tol (value, 4.01, 2.0e-7, "coord: 2.2, 2.8, value: %f", value); + + ok(psImageInterpolate(&value, NULL, NULL, 2.8, 2.8, interp), "Interpolation"); + is_double_tol (value, 4.76, 2.0e-7, "coord: 2.8, 2.8, value: %f", value); + + psFree(interp); skip_end();