Index: trunk/psLib/src/imageops/Makefile.am =================================================================== --- trunk/psLib/src/imageops/Makefile.am (revision 12699) +++ 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 12699) +++ 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 12699) +++ 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 12699) +++ 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 12699) +++ 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 12699) +++ 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;