Changeset 20306
- Timestamp:
- Oct 21, 2008, 4:10:37 PM (18 years ago)
- Location:
- trunk/psLib
- Files:
-
- 4 edited
-
src/imageops/psImageInterpolate.c (modified) (12 diffs)
-
src/imageops/psImageInterpolate.h (modified) (4 diffs)
-
src/imageops/psImagePixelExtract.c (modified) (2 diffs)
-
test/imageops/tap_psImageInterpolate2.c (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/imageops/psImageInterpolate.c
r20207 r20306 7 7 * @author Paul Price, IfA 8 8 * 9 * @version $Revision: 1.2 2$ $Name: not supported by cvs2svn $10 * @date $Date: 2008-10- 16 22:55:45$9 * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2008-10-22 02:10:37 $ 11 11 * 12 12 * Copyright 2004-2007 Institute for Astronomy, University of Hawaii … … 31 31 #include "psImageInterpolate.h" 32 32 33 static void imageInterpolateOptionsFree(psImageInterpolateOptions *options) 33 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 34 // Interpolation kernels 35 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 36 37 // Kernel sizes as a function of interpolation mode: probably faster than a 'switch' 38 static int kernelSizes[] = { 0, // NONE 39 0, // FLAT 40 2, // BILINEAR 41 3, // BIQUADRATIC 42 3, // GAUSS 43 4, // LANCZOS2 44 6, // LANCZOS3 45 8 // LANCZOS4 46 }; 47 48 /// Generate a linear interpolation kernel 49 static inline void interpolationKernelBilinear(psF32 *kernel, // Kernel vector to populate 50 float frac // Fraction of pixel 51 ) 52 { 53 kernel[0] = 1.0 - frac; 54 kernel[1] = frac; 55 } 56 57 /// Generate a biquadratic interpolation kernel 58 /// This reduces Gene's original made-up 2D kernel to 1D 59 static inline void interpolationKernelBiquadratic(psF32 *kernel, // Kernel vector to populate 60 float frac // Fraction of pixel 61 ) 62 { 63 float x = frac / 6.0; 64 float xx = frac * x; 65 kernel[0] = -x + xx; 66 kernel[1] = 1.0/3.0 - 2.0 * xx; 67 kernel[2] = x + xx; 68 } 69 70 #if 0 71 /// Generate a bicubic interpolation kernel 72 /// "cubic Hermite spline" has a=-1/2 for: 73 /// W(x) = 1 when x == 0 74 /// = (a+2)|x|^3 - (a+3)|x|^2 + 1 when 0 < |x| < 1 75 /// = a|x|^3 - 5a|x|^2 + 8a|x| - 4a when 1 < |x| < 2 76 /// = 0 otherwise 77 static inline void interpolationKernelBicubic(psF32 *kernel, // Kernel vector to populate 78 float frac // Fraction of pixel 79 ) 80 { 81 float frac2 = PS_SQR(frac); 82 float frac3 = frac2 * frac; 83 kernel[0] = 0.5 * frac + frac2 - 0.5 * frac3; // x = -(1 + frac) 84 kernel[1] = 1.5 * frac3 - 2.5 * frac2 + 1.0; // x = -frac 85 kernel[2] = 0.5 * frac + 2.0 * frac2 - 1.5 * frac3; // x = 1 - frac 86 kernel[3] = -0.5 * frac2 - frac3; // x = 2 - frac 87 } 88 #endif 89 90 // Generate Lanczos interpolation kernel 91 static inline void interpolationKernelLanczos(psF32 *kernel, // Kernel vector to populate 92 int size, // Size of kernel 93 float frac // Fraction of pixel 94 ) 95 { 96 if (fabs(frac) < FLT_EPSILON) { 97 // No real shift 98 for (int i = 0; i < (size - 1) / 2; i++) { 99 kernel[i] = 0.0; 100 } 101 kernel[(size - 1) / 2] = 1.0; 102 for (int i = (size - 1) / 2 + 1; i < size; i++) { 103 kernel[i] = 0.0; 104 } 105 return; 106 } 107 108 float norm1 = 2.0 / PS_SQR(M_PI); // Normalisation for laczos 109 float norm2 = M_PI * 4.0 / (float)size; // Normalisation for sinc function 1 110 float norm3 = M_PI_2 * 4.0 / (float)size; // Normalisation for sinc function 2 111 float pos = - (size - 1)/2 - frac; // Position of interest 112 for (int i = 0; i < size; i++, pos += 1.0) { 113 if (fabs(pos) < FLT_EPSILON) { 114 kernel[i] = 1.0; 115 } else { 116 kernel[i] = norm1 * sinf(pos * norm2) * sinf(pos * norm3) / PS_SQR(pos); 117 } 118 } 119 return; 120 } 121 122 // Generate Gauss interpolation kernel 123 static inline void interpolationKernelGauss(psF32 *kernel, // Kernel vector to populate 124 int size, // Size of kernel 125 float sigma, // Width of kernel 126 float frac // Fraction of pixel 127 ) 128 { 129 for (int i = 0, pos = - (size - 1) / 2; i < size; i++, pos++) { 130 kernel[i] = expf(-0.5 / PS_SQR(sigma) * PS_SQR(pos - frac)); 131 } 132 // XXX Normalise? 133 return; 134 } 135 136 137 // Generate interpolation kernel 138 static inline void interpolationKernel(psF32 *kernel, // Kernel vector to populate 139 float frac, // Fraction of pixel 140 psImageInterpolateMode mode // Mode of interpolation 141 ) 142 { 143 psAssert(frac >= 0.0 && frac < 1.0, "Pixel fraction is %f", frac); 144 switch (mode) { 145 case PS_INTERPOLATE_FLAT: 146 // Nothing to pre-compute 147 break; 148 case PS_INTERPOLATE_BILINEAR: 149 interpolationKernelBilinear(kernel, frac); 150 break; 151 case PS_INTERPOLATE_BIQUADRATIC: 152 interpolationKernelBiquadratic(kernel, frac); 153 break; 154 case PS_INTERPOLATE_GAUSS: 155 interpolationKernelGauss(kernel, kernelSizes[mode], 0.5, frac); 156 break; 157 case PS_INTERPOLATE_LANCZOS2: 158 case PS_INTERPOLATE_LANCZOS3: 159 case PS_INTERPOLATE_LANCZOS4: 160 interpolationKernelLanczos(kernel, kernelSizes[mode], frac); 161 break; 162 default: 163 psAbort("Unsupported interpolation mode: %x", mode); 164 } 165 166 return; 167 } 168 169 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 170 171 static void imageInterpolationFree(psImageInterpolation *interp) 34 172 { 35 173 // Casting away const 36 psFree((psImage*)options->image); 37 psFree((psImage*)options->mask); 38 psFree((psImage*)options->variance); 39 } 40 41 psImageInterpolateOptions *psImageInterpolateOptionsAlloc(psImageInterpolateMode mode, 42 const psImage *image, const psImage *variance, 43 const psImage *mask, psMaskType maskVal, 44 double badImage, double badVariance, 45 psMaskType badMask, psMaskType poorMask, 46 float poorFrac) 47 { 48 psImageInterpolateOptions *options = psAlloc(sizeof(psImageInterpolateOptions)); // Options, to return 49 psMemSetDeallocator(options, (psFreeFunc)imageInterpolateOptionsFree); 50 51 options->mode = mode; 52 // Casting away const to add to options 53 options->image = psMemIncrRefCounter((psImage*)image); 54 options->variance = psMemIncrRefCounter((psImage*)variance); 55 options->mask = psMemIncrRefCounter((psImage*)mask); 56 options->maskVal = maskVal; 57 options->badImage = badImage; 58 options->badVariance = badVariance; 59 options->badMask = badMask; 60 options->poorMask = poorMask; 61 options->poorFrac = poorFrac; 62 options->shifting = true; 63 64 return options; 65 } 174 psFree((psImage*)interp->image); 175 psFree((psImage*)interp->mask); 176 psFree((psImage*)interp->variance); 177 psFree((psImage*)interp->kernel); 178 psFree((psImage*)interp->kernel2); 179 psFree((psVector*)interp->sumKernel2); 180 } 181 182 psImageInterpolation *psImageInterpolationAlloc(psImageInterpolateMode mode, 183 const psImage *image, const psImage *variance, 184 const psImage *mask, psMaskType maskVal, 185 double badImage, double badVariance, 186 psMaskType badMask, psMaskType poorMask, 187 float poorFrac, int numKernels) 188 { 189 psImageInterpolation *interp = psAlloc(sizeof(psImageInterpolation)); // Options, to return 190 psMemSetDeallocator(interp, (psFreeFunc)imageInterpolationFree); 191 192 interp->mode = mode; 193 // Casting away const to add to interp 194 interp->image = psMemIncrRefCounter((psImage*)image); 195 interp->variance = psMemIncrRefCounter((psImage*)variance); 196 interp->mask = psMemIncrRefCounter((psImage*)mask); 197 interp->maskVal = maskVal; 198 interp->badImage = badImage; 199 interp->badVariance = badVariance; 200 interp->badMask = badMask; 201 interp->poorMask = poorMask; 202 interp->poorFrac = poorFrac; 203 interp->shifting = true; 204 205 // Pre-calculate interpolation kernels 206 psImage *kernel = NULL; // Kernel 207 psImage *kernel2 = NULL; // Kernel^2 208 psVector *sumKernel2 = NULL; // Sum of kernel^2 209 if (numKernels > 0) { 210 int size = kernelSizes[mode]; // Size of kernel 211 212 kernel = psImageAlloc(size, numKernels, PS_TYPE_F32); // Kernel 213 kernel2 = psImageAlloc(size, numKernels, PS_TYPE_F32); // Kernel^2 214 sumKernel2 = psVectorAlloc(numKernels, PS_TYPE_F32); // Sum of kernel^2 215 216 for (int i = 0; i < numKernels; i++) { 217 float frac = i / (float)numKernels; // Fraction of shift 218 interpolationKernel(kernel->data.F32[i], frac, mode); 219 220 float sum = 0.0; // Sum of kernel^2 221 for (int j = 0; j < size; j++) { 222 sum += kernel2->data.F32[i][j] = PS_SQR(kernel->data.F32[i][j]); 223 } 224 sumKernel2->data.F32[i] = sum; 225 } 226 } 227 interp->kernel = kernel; 228 interp->kernel2 = kernel2; 229 interp->sumKernel2 = sumKernel2; 230 interp->numKernels = numKernels; 231 232 return interp; 233 } 234 66 235 67 236 // Interpolation engine for flat mode (nearest pixel) 68 237 static inline psImageInterpolateStatus interpolateFlat(double *imageValue, double *varianceValue, 69 238 psMaskType *maskValue, float x, float y, 70 const psImageInterpolat eOptions *options)239 const psImageInterpolation *interp) 71 240 { 72 241 // Parameters have been checked by psImageInterpolate() 73 242 74 const psImage *image = options->image; // Image of interest243 const psImage *image = interp->image; // Image of interest 75 244 int xInt = round(x - 0.5 + FLT_EPSILON); // Pixel closest to point of interest in x 76 245 int yInt = round(y - 0.5 + FLT_EPSILON); // Pixel closest to point of interest in y … … 81 250 // At least one pixel of the interpolation kernel is off the image 82 251 if (imageValue) { 83 *imageValue = options->badImage;252 *imageValue = interp->badImage; 84 253 } 85 254 if (varianceValue) { 86 *varianceValue = options->badVariance;255 *varianceValue = interp->badVariance; 87 256 } 88 257 if (maskValue) { 89 *maskValue = options->badMask;258 *maskValue = interp->badMask; 90 259 } 91 260 … … 97 266 case PS_TYPE_##TYPE: { \ 98 267 if (imageValue) { \ 99 *imageValue = options->image->data.TYPE[yInt][xInt]; \268 *imageValue = interp->image->data.TYPE[yInt][xInt]; \ 100 269 } \ 101 if (varianceValue && options->variance) { \102 *varianceValue = options->variance->data.TYPE[yInt][xInt]; \270 if (varianceValue && interp->variance) { \ 271 *varianceValue = interp->variance->data.TYPE[yInt][xInt]; \ 103 272 } \ 104 273 break; \ 105 274 } 106 275 107 switch ( options->image->type.type) {276 switch (interp->image->type.type) { 108 277 FLAT_CASE(U8); 109 278 FLAT_CASE(U16); … … 118 287 default: 119 288 psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x", 120 options->image->type.type);289 interp->image->type.type); 121 290 return PS_INTERPOLATE_STATUS_ERROR; 122 291 } 123 292 124 293 if (maskValue) { 125 if ( options->mask) {126 *maskValue = options->mask->data.PS_TYPE_MASK_DATA[yInt][xInt];294 if (interp->mask) { 295 *maskValue = interp->mask->data.PS_TYPE_MASK_DATA[yInt][xInt]; 127 296 } else { 128 297 *maskValue = 0; … … 133 302 } 134 303 135 // Setup for interpolation by kernel136 static inline void interpolateKernelSetup(int *xNum, int *yNum, // Size of interpolation kernel, returned137 int *xCentral, int *yCentral, // Central pixel of convolution138 float x, float y, // Coordinates of interest139 psImageInterpolateMode mode // Mode for interpolation140 )141 {142 switch (mode) {143 case PS_INTERPOLATE_BILINEAR:144 *xNum = *yNum = 2;145 // Central pixel is the pixel below the point of interest146 *xCentral = floor(x - 0.5 + FLT_EPSILON);147 *yCentral = floor(y - 0.5 + FLT_EPSILON);148 break;149 case PS_INTERPOLATE_BICUBE:150 case PS_INTERPOLATE_GAUSS:151 *xNum = *yNum = 3;152 // Central pixel is the closest pixel to the point of interest153 *xCentral = x;154 *yCentral = y;155 break;156 case PS_INTERPOLATE_FLAT:157 case PS_INTERPOLATE_LANCZOS2:158 case PS_INTERPOLATE_LANCZOS3:159 case PS_INTERPOLATE_LANCZOS4:160 default:161 psAbort("Invalid interpolation mode.");162 }163 }164 165 // Generate the interpolation kernel166 // No need to normalise to unity167 static inline void interpolateKernelGenerate(int xNum, int yNum, // Size of interpolation kernel168 float kernel[xNum][yNum], // Kernel, to be set169 bool *xExact, bool *yExact, // Exact shift?170 int xCentral, int yCentral, // Central pixel of convolution171 float x, float y, // Coordinates of interest172 psImageInterpolateMode mode, // Mode for interpolation173 bool allowExact // Allow exact shifts?174 )175 {176 float xFrac = x - 0.5 - xCentral; // Fraction of pixel in x177 float yFrac = y - 0.5 - yCentral; // Fraction of pixel in y178 *xExact = (allowExact && fabs(xFrac) < DBL_EPSILON);179 *yExact = (allowExact && fabs(yFrac) < DBL_EPSILON);180 181 switch (mode) {182 case PS_INTERPOLATE_BILINEAR: {183 kernel[0][0] = (1.0 - xFrac) * (1.0 - yFrac);184 kernel[0][1] = xFrac * (1.0 - yFrac);185 kernel[1][0] = yFrac * (1.0 - xFrac);186 kernel[1][1] = xFrac * yFrac;187 break;188 }189 case PS_INTERPOLATE_BICUBE: {190 float xFrac = x - 0.5 - xCentral; // Fraction of pixel in x191 float yFrac = y - 0.5 - yCentral; // Fraction of pixel in y192 // Calculation variables193 double xxFrac = xFrac * xFrac / 6.0;194 double yyFrac = yFrac * yFrac / 6.0;195 double xyFrac = 0.25 * xFrac * yFrac;196 xFrac /= 6.0;197 yFrac /= 6.0;198 kernel[0][0] = - 1.0/9.0 - xFrac - yFrac + xxFrac + yyFrac + xyFrac;199 kernel[0][1] = 2.0/9.0 - yFrac - 2.0 * xxFrac + yyFrac;200 kernel[0][2] = - 1.0/9.0 + xFrac - yFrac + xxFrac + yyFrac - xyFrac;201 kernel[1][0] = 2.0/9.0 - xFrac + xxFrac - 2.0 * yyFrac;202 kernel[1][1] = 5.0/9.0 - 2.0 * xxFrac - 2.0 * yyFrac;203 kernel[1][2] = 2.0/9.0 + xFrac + xxFrac - 2.0 * yyFrac;204 kernel[2][0] = - 1.0/9.0 - xFrac + yFrac + xxFrac + yyFrac - xyFrac;205 kernel[2][1] = 2.0/9.0 + yFrac - 2.0 * xxFrac + yyFrac;206 kernel[2][2] = - 1.0/9.0 + xFrac + yFrac + xxFrac + yyFrac + xyFrac;207 break;208 }209 case PS_INTERPOLATE_GAUSS: {210 float xFrac = x - xCentral - 0.5; // Fraction of pixel in x211 float yFrac = y - yCentral - 0.5; // Fraction of pixel in y212 float sigma = 0.5; // Gaussian sigma213 for (int j = 0, yPos = - (yNum - 1) / 2; j < yNum; j++, yPos++) {214 for (int i = 0, xPos = - (xNum - 1) / 2; i < xNum; i++, xPos++) {215 kernel[j][i] = expf(-0.5 / PS_SQR(sigma) * (PS_SQR(xPos - xFrac) +216 PS_SQR(yPos - yFrac)));217 }218 }219 break;220 }221 case PS_INTERPOLATE_FLAT:222 case PS_INTERPOLATE_LANCZOS2:223 case PS_INTERPOLATE_LANCZOS3:224 case PS_INTERPOLATE_LANCZOS4:225 default:226 psAbort("Invalid interpolation mode.");227 }228 return;229 }230 304 231 305 // Can't interpolate: kernel is off the image 232 306 #define INTERPOLATE_OFF() { \ 233 *imageValue = options->badImage; \307 *imageValue = interp->badImage; \ 234 308 if (varianceValue) { \ 235 *varianceValue = options->badVariance; \309 *varianceValue = interp->badVariance; \ 236 310 } \ 237 311 if (maskValue) { \ 238 *maskValue = options->badMask; \312 *maskValue = interp->badMask; \ 239 313 } \ 240 314 return PS_INTERPOLATE_STATUS_OFF; \ 241 315 } 242 316 243 // Set up and check interpolation bounds 244 // This macro defines many useful values 245 #define INTERPOLATE_BOUNDS() \ 246 int xLast = image->numCols - 1, yLast = image->numRows - 1; /* Last pixels of image */ \ 247 /* Start and stop of kernel on image */ \ 248 int xStart = xCentral - (xNum - 1) / 2, xStop = xCentral + xNum / 2; \ 249 int yStart = yCentral - (yNum - 1) / 2, yStop = yCentral + yNum / 2; \ 250 if (xStart > xLast || xStop < 0 || yStart > yLast || yStop < 0) { \ 251 /* Interpolation kernel is entirely off the image */ \ 252 INTERPOLATE_OFF(); \ 253 } \ 254 int xMin, xMax, yMin, yMax; /* Minimum and maximum valid pixels on image */ \ 255 int iMin, iMax, jMin, jMax; /* Minimum and maximum valid pixels on kernel */ \ 256 bool offImage = false; /* At least one pixel of the kernel is off the image */ \ 257 if (xStart < 0) { \ 258 xMin = 0; \ 259 iMin = -xStart; \ 260 offImage = true; \ 261 } else { \ 262 xMin = xStart; \ 263 iMin = 0; \ 264 } \ 265 if (xStop > xLast) { \ 266 xMax = xLast; \ 267 iMax = xLast - xStart + 1; \ 268 offImage = true; \ 269 } else { \ 270 xMax = xStop; \ 271 iMax = xNum; \ 272 } \ 273 if (yStart < 0) { \ 274 yMin = 0; \ 275 jMin = -yStart; \ 276 offImage = true; \ 277 } else { \ 278 yMin = yStart; \ 279 jMin = 0; \ 280 } \ 281 if (yStop > yLast) { \ 282 yMax = yLast; \ 283 jMax = yLast - yStart + 1; \ 284 offImage = true; \ 285 } else { \ 286 yMax = yStop; \ 287 jMax = yNum; \ 288 } 289 290 // Refine limits if the interpolation is exact 291 #define INTERPOLATE_EXACT() { \ 292 if (xExact) { \ 293 if (iMin > 0 || iMax < 1) { \ 294 INTERPOLATE_OFF(); /* Returns */ \ 295 } \ 296 iMin = (xNum - 1) / 2; \ 297 iMax = iMin + 1; \ 298 } \ 299 if (yExact) { \ 300 if (jMin > 0 || jMax < 1) { \ 301 INTERPOLATE_OFF(); /* Returns */ \ 302 } \ 303 jMin = (yNum - 1) / 2; \ 304 jMax = jMin + 1; \ 305 } \ 306 if (xExact && yExact) { \ 307 return interpolateFlat(imageValue, varianceValue, maskValue, x, y, options); \ 308 } \ 309 } 310 311 // Interpolation engine using interpolation kernel 317 // Set up the kernel parameters; defines some useful values 318 #define INTERPOLATE_KERNEL_SETUP(X, Y) \ 319 /* Central pixel is the pixel below the point of interest */ \ 320 int xCentral = floor((X) - 0.5), yCentral = floor((Y) - 0.5); /* Central pixel */ \ 321 float xFrac = (X) - xCentral - 0.5, yFrac = (Y) - yCentral - 0.5; /* Fraction of pixel */ 322 323 324 // Interpolation engine for (separable) interpolation kernels 312 325 static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue, 313 326 psMaskType *maskValue, float x, float y, 314 const psImageInterpolat eOptions *options)327 const psImageInterpolation *interp) 315 328 { 316 329 // Parameters have been checked by psImageInterpolate() 317 330 318 int xNum, yNum; // Number of interpolation kernel pixels 319 int xCentral, yCentral; // Central pixel of the convolution 320 interpolateKernelSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode); 321 322 const psImage *image = options->image; // Image of interest 323 const psImage *mask = options->mask; // Image mask 324 const psImage *variance = options->variance; // Image variance 325 326 INTERPOLATE_BOUNDS(); 327 328 float kernel[yNum][xNum]; // Interpolation kernel for straight interpolation 329 bool xExact, yExact; // Exact shifts? 330 interpolateKernelGenerate(xNum, yNum, kernel, &xExact, &yExact, xCentral, yCentral, 331 x, y, options->mode, options->shifting); 331 psImageInterpolateMode mode = interp->mode; // Mode of interpolation 332 const psImage *image = interp->image; // Image of interest 333 const psImage *mask = interp->mask; // Image mask 334 const psImage *variance = interp->variance; // Image variance 335 psMaskType maskVal = interp->maskVal; // Value to mask 336 bool wantVariance = variance && varianceValue; // Does the user want the variance value? 337 bool haveMask = mask && maskVal; // Does the user want the variance value? 338 339 // Kernel basics 340 int size = kernelSizes[mode]; // Size of kernel 341 INTERPOLATE_KERNEL_SETUP(x, y); 342 343 // Extent of the kernel on the image 344 bool xExact = fabsf(xFrac) < FLT_EPSILON, yExact = fabsf(yFrac) < FLT_EPSILON; // Are shifts exact? 345 if (xExact && yExact) { 346 // Both shifts are exact 347 return interpolateFlat(imageValue, varianceValue, maskValue, x, y, interp); 348 } 349 350 int xLast = image->numCols - 1, yLast = image->numRows - 1; // Last pixels of image 351 int xStart = xCentral - (size - 1) / 2, xStop = xCentral + size / 2; // Start and stop of kernel on image 352 int yStart = yCentral - (size - 1) / 2, yStop = yCentral + size / 2; // Start and stop of kernel on image 353 if (xStart > xLast || xStop < 0 || yStart > yLast || yStop < 0) { 354 // Interpolation kernel is entirely off the image 355 INTERPOLATE_OFF(); // Returns 356 } 357 int xMin, xMax, yMin, yMax; // Minimum and maximum valid pixels on image 358 int iMin, iMax, jMin, jMax; // Minimum and maximum valid pixels on kernel 359 bool offImage = false; // At least one pixel of the kernel is off the image 360 361 // XXX Haven't got single dimension exact shifts implemented properly yet 362 // XXX When it is ready, note that the check below limits 363 xExact = yExact = false; 364 if (xExact) { 365 if (iMin > 0 || iMax < 1) { 366 INTERPOLATE_OFF(); // Returns 367 } 368 iMin = (size - 1) / 2; 369 iMax = iMin + 1; 370 } else { 371 if (xStart < 0) { 372 xMin = 0; 373 iMin = -xStart; 374 offImage = true; 375 } else { 376 xMin = xStart; 377 iMin = 0; 378 } 379 if (xStop > xLast) { 380 xMax = xLast; 381 iMax = xLast - xStart + 1; 382 offImage = true; 383 } else { 384 xMax = xStop; 385 iMax = size; 386 } 387 } 388 if (yExact) { 389 if (jMin > 0 || jMax < 1) { 390 INTERPOLATE_OFF(); // Returns 391 } 392 jMin = (size - 1) / 2; 393 jMax = jMin + 1; 394 } else { 395 if (yStart < 0) { 396 yMin = 0; 397 jMin = -yStart; 398 offImage = true; 399 } else { 400 yMin = yStart; 401 jMin = 0; 402 } 403 if (yStop > yLast) { 404 yMax = yLast; 405 jMax = yLast - yStart + 1; 406 offImage = true; 407 } else { 408 yMax = yStop; 409 jMax = size; 410 } 411 } 412 413 // Get the appropriate kernels 414 psVector *xKernelNew = NULL, *yKernelNew = NULL; // New kernels if not pre-calculated 415 psVector *xKernel2New = NULL, *yKernel2New = NULL;// New kernel^2 if not pre-calculated 416 psF32 *xKernel, *yKernel; // Kernel values 417 psF32 *xKernel2, *yKernel2; // Kernel^2 values 418 float xSumKernel2, ySumKernel2; // Sum of kernel^2 419 int numKernels = interp->numKernels; // Number of pre-calculated kernels 420 if (numKernels > 0) { 421 const psImage *kernel = interp->kernel; // Pre-calculated kernels 422 const psImage *kernel2 = interp->kernel2; // Pre-calculated kernel^2 423 const psVector *sumKernel2 = interp->sumKernel2; // Pre-calculated kernel^2 424 int xIndex = xFrac * numKernels + 0.5; // Index of closest pre-calculated kernel 425 xKernel = kernel->data.F32[xIndex]; 426 xKernel2 = kernel2->data.F32[xIndex]; 427 xSumKernel2 = sumKernel2->data.F32[xIndex]; 428 int yIndex = yFrac * numKernels + 0.5; // Index of closest pre-calculated kernel 429 yKernel = kernel->data.F32[yIndex]; 430 yKernel2 = kernel2->data.F32[yIndex]; 431 ySumKernel2 = sumKernel2->data.F32[yIndex]; 432 } else { 433 xKernelNew = psVectorAlloc(size, PS_TYPE_F32); 434 xKernel = xKernelNew->data.F32; 435 interpolationKernel(xKernel, xFrac, mode); 436 yKernelNew = psVectorAlloc(size, PS_TYPE_F32); 437 yKernel = yKernelNew->data.F32; 438 interpolationKernel(yKernel, yFrac, mode); 439 440 if (haveMask || wantVariance || offImage) { 441 xKernel2New = psVectorAlloc(size, PS_TYPE_F32); 442 xKernel2 = xKernel2New->data.F32; 443 interpolationKernel(xKernel2, xFrac, mode); 444 yKernel2New = psVectorAlloc(size, PS_TYPE_F32); 445 yKernel2 = yKernel2New->data.F32; 446 interpolationKernel(yKernel2, yFrac, mode); 447 xSumKernel2 = 0.0; 448 ySumKernel2 = 0.0; 449 for (int i = 0; i < size; i++) { 450 xSumKernel2 += xKernel2[i] = PS_SQR(xKernel2[i]); 451 ySumKernel2 += yKernel2[i] = PS_SQR(yKernel2[i]); 452 } 453 } 454 } 332 455 333 456 float sumKernel = 0.0; // Sum of the kernel 334 double sumKernel2 = 0.0; // Sum of the kernel-squared 335 float sumBad = 0.0; // Sum of bad kernel-squared contributions 336 psMaskType maskVal = options->maskVal; // Value to mask 457 double sumBad = 0.0; // Sum of bad kernel-squared contributions 337 458 double sumImage = 0.0; // Sum of image multiplied by kernel 338 459 double sumVariance = 0.0; // Sum of variance multiplied by kernel-squared 339 340 bool wantVariance = variance && varianceValue; // Does the user want the variance value? 341 bool haveMask = mask && maskVal; // Does the user want the variance value? 342 343 INTERPOLATE_EXACT(); 344 460 float sumKernel2 = xSumKernel2 * ySumKernel2; // Sum of kernel^2 461 462 // Measure kernel contribution from outside image 463 464 if (offImage) { 345 465 // Add contributions in an area outside the image 346 #define INTERPOLATE_KERNEL_ADD_OFFIMAGE() { \ 347 sumBad += PS_SQR(kernel[j][i]); \ 348 } 349 350 // Measure kernel contribution from outside image 351 if (offImage) { 466 #define INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL() \ 467 float xSumBad = 0.0; 468 #define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL() { \ 469 xSumBad += xKernel2[i]; \ 470 } 471 #define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW() { \ 472 sumBad += yKernel2[j] * xSumBad; \ 473 } 474 475 // Bottom rows 352 476 if (!yExact) { 353 // Bottom rows354 477 for (int j = 0; j < jMin; j++) { 355 for (int i = 0; i < xNum; i++) { 356 INTERPOLATE_KERNEL_ADD_OFFIMAGE(); 478 INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL(); 479 for (int i = 0; i < size; i++) { 480 INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); 357 481 } 482 INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW(); 358 483 } 359 484 } 485 // Two sides 360 486 if (!xExact) { 361 // Two sides362 487 for (int j = jMin; j < jMax; j++) { 488 INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL(); 363 489 for (int i = 0; i < iMin; i++) { 364 INTERPOLATE_ KERNEL_ADD_OFFIMAGE();490 INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); 365 491 } 366 for (int i = iMax + 1; i < xNum; i++) {367 INTERPOLATE_ KERNEL_ADD_OFFIMAGE();492 for (int i = iMax; i < size; i++) { 493 INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); 368 494 } 495 INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW(); 369 496 } 370 497 } 498 // Top rows 371 499 if (!yExact) { 372 // Top rows 373 for (int j = jMax + 1; j < yNum; j++) { 374 for (int i = 0; i < xNum; i++) { 375 INTERPOLATE_KERNEL_ADD_OFFIMAGE(); 500 for (int j = jMax; j < size; j++) { 501 if (!xExact) { 502 INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL(); 503 for (int i = 0; i < size; i++) { 504 INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL(); 505 } 506 INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW(); 376 507 } 377 508 } … … 384 515 if (haveMask) { \ 385 516 /* Variance and mask */ \ 386 for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ 387 for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ 388 float kernelValue = kernel[j][i]; /* Value of kernel */ \ 389 float kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ 390 if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \ 391 sumBad += kernelValue2; \ 517 const psF32 *yKernelData = yKernel; \ 518 const psF32 *yKernel2Data = yKernel2; \ 519 for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++, yKernelData++, yKernel2Data++) { \ 520 double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ 521 double xSumVariance = 0.0; /* Sum of image multiplied by kernel-squared in x */ \ 522 float xSumKernel = 0.0; /* Sum of kernel in x */ \ 523 float xSumBad = 0.0; /* Sum of bad kernel-squared in x */ \ 524 /* Dereferenced versions of inputs */ \ 525 const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \ 526 const ps##TYPE *varianceData = &variance->data.TYPE[yPix][xMin]; \ 527 const psMaskType *maskData = &mask->data.PS_TYPE_MASK_DATA[yPix][xMin]; \ 528 const psF32 *xKernelData = xKernel; \ 529 const psF32 *xKernel2Data = xKernel2; \ 530 for (int i = iMin, xPix = xMin; i < iMax; \ 531 i++, xPix++, imageData++, varianceData++, maskData++, \ 532 xKernelData++, xKernel2Data++) { \ 533 float kernelValue = *xKernelData; /* Value of kernel in x */ \ 534 float kernelValue2 = *xKernel2Data; /* Square of kernel in x */ \ 535 if (*maskData & maskVal) { \ 536 xSumBad += kernelValue2; \ 392 537 } else { \ 393 sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ 394 sumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \ 395 sumKernel += kernelValue; \ 396 sumKernel2 += kernelValue2; \ 538 xSumImage += kernelValue * *imageData; \ 539 xSumVariance += kernelValue2 * *varianceData; \ 540 xSumKernel += kernelValue; \ 397 541 } \ 398 542 } \ 543 float kernelValue = *yKernel; /* Value of kernel in y */ \ 544 float kernelValue2 = *yKernel2Data; /* Value of kernel-squared in y */ \ 545 sumImage += kernelValue * xSumImage; \ 546 sumVariance += kernelValue2 * xSumVariance; \ 547 sumBad += kernelValue2 * xSumBad; \ 548 sumKernel += kernelValue * xSumKernel; \ 399 549 } \ 400 \401 550 } else { \ 402 551 /* Variance, no mask */ \ 403 for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ 404 for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ 405 float kernelValue = kernel[j][i]; /* Value of kernel */ \ 406 float kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ 407 sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ 408 sumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \ 409 sumKernel += kernelValue; \ 410 sumKernel2 += kernelValue2; \ 552 const psF32 *yKernelData = yKernel; \ 553 const psF32 *yKernel2Data = yKernel2; \ 554 for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++, yKernelData++, yKernel2Data++) { \ 555 double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ 556 double xSumVariance = 0.0; /* Sum of image multiplied by kernel-squared in x */ \ 557 float xSumKernel = 0.0; /* Sum of kernel in x */ \ 558 /* Dereferenced versions of inputs */ \ 559 const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \ 560 const ps##TYPE *varianceData = &variance->data.TYPE[yPix][xMin]; \ 561 const psF32 *xKernelData = xKernel; \ 562 const psF32 *xKernel2Data = xKernel2; \ 563 for (int i = iMin, xPix = xMin; i < iMax; \ 564 i++, xPix++, imageData++, varianceData++, xKernelData++, xKernel2Data++) { \ 565 float kernelValue = *xKernelData; /* Value of kernel */ \ 566 float kernelValue2 = *xKernel2Data; /* Square of kernel */ \ 567 xSumImage += kernelValue * *imageData; \ 568 xSumVariance += kernelValue2 * *varianceData; \ 569 xSumKernel += kernelValue; \ 411 570 } \ 571 float kernelValue = *yKernelData; /* Value of kernel in y */ \ 572 float kernelValue2 = *yKernel2Data; /* Value of kernel-squared in y */ \ 573 sumImage += kernelValue * xSumImage; \ 574 sumVariance += kernelValue2 * xSumVariance; \ 575 sumKernel += kernelValue * xSumKernel; \ 412 576 } \ 413 577 } \ 414 578 } else if (haveMask) { \ 415 579 /* Mask, no variance */ \ 416 for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ 417 for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ 418 float kernelValue = kernel[j][i]; /* Value of kernel */ \ 419 float kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \ 420 if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \ 421 sumBad += kernelValue2; \ 580 const psF32 *yKernelData = yKernel; \ 581 const psF32 *yKernel2Data = yKernel2; \ 582 for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++, yKernelData++, yKernel2Data++) { \ 583 double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ 584 float xSumKernel = 0.0; /* Sum of kernel in x */ \ 585 float xSumBad = 0.0; /* Sum of bad kernel-squared in x */ \ 586 /* Dereferenced versions of inputs */ \ 587 const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \ 588 const psMaskType *maskData = &mask->data.PS_TYPE_MASK_DATA[yPix][xMin]; \ 589 const psF32 *xKernelData = xKernel; \ 590 const psF32 *xKernel2Data = xKernel2; \ 591 for (int i = iMin, xPix = xMin; i < iMax; \ 592 i++, xPix++, imageData++, maskData++, xKernelData++, xKernel2Data++) { \ 593 float kernelValue = *xKernelData; /* Value of kernel */ \ 594 float kernelValue2 = *xKernel2Data; /* Value of kernel-squared */ \ 595 if (*maskData & maskVal) { \ 596 xSumBad += kernelValue2; \ 422 597 } else { \ 423 sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ 424 sumKernel += kernelValue; \ 425 sumKernel2 += kernelValue2; \ 598 xSumImage += kernelValue * *imageData; \ 599 xSumKernel += kernelValue; \ 426 600 } \ 427 601 } \ 602 float kernelValue = *yKernelData; /* Value of kernel in y */ \ 603 float kernelValue2 = *yKernel2Data; /* Value of kernel-squared in y */ \ 604 sumImage += kernelValue * xSumImage; \ 605 sumBad += kernelValue2 * xSumBad; \ 606 sumKernel += kernelValue * xSumKernel; \ 428 607 } \ 429 } else { 608 } else {\ 430 609 /* Neither variance nor mask */ \ 431 for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \ 432 for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \ 433 float kernelValue = kernel[j][i]; /* Value of kernel */ \ 434 sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \ 435 sumKernel += kernelValue; \ 610 const psF32 *yKernelData = yKernel; \ 611 for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++, yKernelData++) { \ 612 double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \ 613 float xSumKernel = 0.0; /* Sum of kernel in x */ \ 614 /* Dereferenced versions of inputs */ \ 615 const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \ 616 const psF32 *xKernelData = xKernel; \ 617 for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++, imageData++, xKernelData++) { \ 618 float kernelValue = *xKernelData; /* Value of kernel */ \ 619 xSumImage += kernelValue * *imageData; \ 620 xSumKernel += kernelValue; \ 436 621 } \ 622 float kernelValue = *yKernelData; /* Value of kernel in y */ \ 623 sumImage += kernelValue * xSumImage; \ 624 sumKernel += kernelValue * xSumKernel; \ 437 625 } \ 438 626 } \ … … 466 654 // Completely good pixel 467 655 status = PS_INTERPOLATE_STATUS_GOOD; 468 } else if (sumBad < PS_SQR( options->poorFrac) * (sumBad + sumKernel2)) {656 } else if (sumBad < PS_SQR(interp->poorFrac) * (sumBad + sumKernel2) ) { 469 657 // Some pixels masked: poor pixel 470 658 if (haveMask && maskValue) { 471 *maskValue |= options->poorMask;659 *maskValue |= interp->poorMask; 472 660 } 473 661 status = PS_INTERPOLATE_STATUS_POOR; … … 475 663 // Many pixels (or a few important pixels) masked: bad pixel 476 664 if (haveMask && maskValue) { 477 *maskValue |= options->badMask;665 *maskValue |= interp->badMask; 478 666 } 479 667 status = PS_INTERPOLATE_STATUS_BAD; 480 668 } 481 669 670 psFree(xKernelNew); 671 psFree(yKernelNew); 672 psFree(xKernel2New); 673 psFree(yKernel2New); 674 482 675 return status; 483 676 } 484 677 485 678 486 // Generate Lanczos interpolation kernels487 static void lanczos(float values[], // Interpolation kernel to generate488 bool *exact, // Exact shift?489 int num, // Number of values in the kernel490 float frac, // Sub-pixel position491 bool allowExact // Allow exact shift?492 )493 {494 if (allowExact && fabs(frac) < DBL_EPSILON) {495 *exact = true;496 // No real shift497 for (int i = 0; i < (num - 1) / 2; i++) {498 values[i] = 0.0;499 }500 values[(num - 1) / 2] = 1.0;501 for (int i = (num - 1) / 2 + 1; i < num; i++) {502 values[i] = 0.0;503 }504 } else {505 *exact = false;506 float norm1 = 2.0 / PS_SQR(M_PI); // Normalisation for laczos507 float norm2 = M_PI * 4.0 / (float)num; // Normalisation for sinc function 1508 float norm3 = M_PI_2 * 4.0 / (float)num; // Normalisation for sinc function 2509 float pos = - (num - 1)/2 - frac; // Position of interest510 for (int i = 0; i < num; i++, pos += 1.0) {511 if (fabs(pos) < DBL_EPSILON) {512 values[i] = 1.0;513 } else {514 values[i] = norm1 * sinf(pos * norm2) * sinf(pos * norm3) / PS_SQR(pos);515 }516 }517 }518 519 return;520 }521 522 // Setup for interpolation by separable kernels523 static inline void interpolateSeparateSetup(int *xNum, int *yNum, // Size of interpolation kernel, returned524 int *xCentral, int *yCentral, // Central pixel of convolution525 float x, float y, // Coordinates of interest526 psImageInterpolateMode mode // Mode for interpolation527 )528 {529 // Central pixel is the pixel below the point of interest530 *xCentral = floor(x - 0.5);531 *yCentral = floor(y - 0.5);532 switch (mode) {533 case PS_INTERPOLATE_LANCZOS2:534 *xNum = *yNum = 4;535 break;536 case PS_INTERPOLATE_LANCZOS3:537 *xNum = *yNum = 6;538 break;539 case PS_INTERPOLATE_LANCZOS4:540 *xNum = *yNum = 8;541 break;542 case PS_INTERPOLATE_FLAT:543 case PS_INTERPOLATE_BILINEAR:544 case PS_INTERPOLATE_BICUBE:545 case PS_INTERPOLATE_GAUSS:546 default:547 psAbort("Invalid interpolation mode.");548 }549 }550 551 // Generate the interpolation kernels for separable case; they should be normalised to unity552 static inline void interpolateSeparateGenerate(int xNum, int yNum, // Size of interpolation kernel553 float xKernel[xNum], float yKernel[yNum], // Kernels554 bool *xExact, bool *yExact, // Exact shifts?555 int xCentral, int yCentral, // Central pixel of convolution556 float x, float y, // Coordinates of interest557 psImageInterpolateMode mode, // Mode for interpolation558 bool allowExact // Allow an exact shift?559 )560 {561 // XXX Could put in an "exact shift" (i.e., xFrac = 0.0) version562 switch (mode) {563 case PS_INTERPOLATE_LANCZOS2:564 case PS_INTERPOLATE_LANCZOS3:565 case PS_INTERPOLATE_LANCZOS4: {566 float xFrac = x - xCentral - 0.5; // Fraction of pixel in x567 lanczos(xKernel, xExact, xNum, xFrac, allowExact);568 float yFrac = y - yCentral - 0.5; // Fraction of pixel in y569 lanczos(yKernel, yExact, yNum, yFrac, allowExact);570 break;571 }572 case PS_INTERPOLATE_FLAT:573 case PS_INTERPOLATE_BILINEAR:574 case PS_INTERPOLATE_BICUBE:575 case PS_INTERPOLATE_GAUSS:576 default:577 psAbort("Invalid interpolation mode.");578 }579 }580 581 // Interpolation engine for separable interpolation kernels (either for good reasons or for practical reasons)582 static psImageInterpolateStatus interpolateSeparate(double *imageValue, double *varianceValue,583 psMaskType *maskValue, float x, float y,584 const psImageInterpolateOptions *options)585 {586 // Parameters have been checked by psImageInterpolate()587 588 int xNum, yNum; // Number of interpolation kernel pixels589 int xCentral, yCentral; // Central pixel of the convolution590 interpolateSeparateSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode);591 592 const psImage *image = options->image; // Image of interest593 const psImage *mask = options->mask; // Image mask594 const psImage *variance = options->variance; // Image variance595 596 INTERPOLATE_BOUNDS();597 598 float xKernel[xNum], yKernel[yNum]; // Interpolation kernels599 bool xExact, yExact; // Exact shift?600 interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, &xExact, &yExact,601 xCentral, yCentral, x, y, options->mode, options->shifting);602 603 float sumKernel = 0.0; // Sum of the kernel604 double sumKernel2 = 0.0; // Sum of the kernel-squared605 double sumBad = 0.0; // Sum of bad kernel-squared contributions606 psMaskType maskVal = options->maskVal; // Value to mask607 double sumImage = 0.0; // Sum of image multiplied by kernel608 double sumVariance = 0.0; // Sum of variance multiplied by kernel-squared609 610 bool wantVariance = variance && varianceValue; // Does the user want the variance value?611 bool haveMask = mask && maskVal; // Does the user want the variance value?612 613 INTERPOLATE_EXACT();614 615 // Add contributions in an area outside the image616 #define INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL() \617 float xSumBad = 0.0;618 619 #define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL() { \620 xSumBad += PS_SQR(xKernel[i]); \621 }622 623 #define INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW() { \624 sumBad += PS_SQR(yKernel[j]) * xSumBad; \625 }626 627 // Measure kernel contribution from outside image628 if (offImage) {629 // Bottom rows630 if (!yExact) {631 for (int j = 0; j < jMin; j++) {632 INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();633 for (int i = 0; i < xNum; i++) {634 INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();635 }636 INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();637 }638 }639 // Two sides640 if (!xExact) {641 for (int j = jMin; j < jMax; j++) {642 INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();643 for (int i = 0; i < iMin; i++) {644 INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();645 }646 for (int i = iMax + 1; i < xNum; i++) {647 INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();648 }649 INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();650 }651 }652 // Top rows653 if (!yExact) {654 for (int j = jMax + 1; j < yNum; j++) {655 if (!xExact) {656 INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();657 for (int i = 0; i < xNum; i++) {658 INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();659 }660 INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();661 }662 }663 }664 }665 666 #define INTERPOLATE_SEPARATE_CASE(TYPE) \667 case PS_TYPE_##TYPE: { \668 if (wantVariance) { \669 if (haveMask) { \670 /* Variance and mask */ \671 for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \672 double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \673 double xSumVariance = 0.0; /* Sum of image multiplied by kernel-squared in x */ \674 float xSumKernel = 0.0; /* Sum of kernel in x */ \675 double xSumKernel2 = 0.0; /* Sum of kernel-squared in x */ \676 float xSumBad = 0.0; /* Sum of bad kernel-squared in x */ \677 /* Dereferenced versions of inputs */ \678 const ps##TYPE *imageData = &image->data.TYPE[yPix][xMin]; \679 const ps##TYPE *varianceData = &variance->data.TYPE[yPix][xMin]; \680 const psMaskType *maskData = &mask->data.PS_TYPE_MASK_DATA[yPix][xMin]; \681 for (int i = iMin, xPix = xMin; i < iMax; \682 i++, xPix++, imageData++, varianceData++, maskData++) { \683 float kernelValue = xKernel[i]; /* Value of kernel in x */ \684 double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel in x */ \685 if (*maskData & maskVal) { \686 xSumBad += kernelValue2; \687 } else { \688 xSumImage += kernelValue * *imageData; \689 xSumVariance += kernelValue2 * *varianceData; \690 xSumKernel += kernelValue; \691 xSumKernel2 += kernelValue2; \692 } \693 } \694 float kernelValue = yKernel[j]; /* Value of kernel in y */ \695 double kernelValue2 = PS_SQR(kernelValue); /* Value of kernel-squared in y */ \696 sumImage += kernelValue * xSumImage; \697 sumVariance += kernelValue2 * xSumVariance; \698 sumBad += kernelValue2 * xSumBad; \699 sumKernel += kernelValue * xSumKernel; \700 sumKernel2 += kernelValue2 * xSumKernel2; \701 } \702 } else { \703 /* Variance, no mask */ \704 for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \705 double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \706 double xSumVariance = 0.0; /* Sum of image multiplied by kernel-squared in x */ \707 float xSumKernel = 0.0; /* Sum of kernel in x */ \708 double xSumKernel2 = 0.0; /* Sum of kernel-squared in x */ \709 for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \710 float kernelValue = xKernel[i]; /* Value of kernel */ \711 double kernelValue2 = PS_SQR(kernelValue); /* Square of kernel */ \712 xSumImage += kernelValue * image->data.TYPE[yPix][xPix]; \713 xSumVariance += kernelValue2 * variance->data.TYPE[yPix][xPix]; \714 xSumKernel += kernelValue; \715 xSumKernel2 += kernelValue2; \716 } \717 float kernelValue = yKernel[j]; /* Value of kernel in y */ \718 double kernelValue2 = PS_SQR(kernelValue); /* Value of kernel-squared in y */ \719 sumImage += kernelValue * xSumImage; \720 sumVariance += kernelValue2 * xSumVariance; \721 sumKernel += kernelValue * xSumKernel; \722 sumKernel2 += kernelValue2 * xSumKernel2; \723 } \724 } \725 } else if (haveMask) { \726 /* Mask, no variance */ \727 for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \728 double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \729 float xSumKernel = 0.0; /* Sum of kernel in x */ \730 double xSumKernel2 = 0.0; /* Sum of kernel-squared in x */ \731 float xSumBad = 0.0; /* Sum of bad kernel-squared in x */ \732 for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \733 float kernelValue = xKernel[i]; /* Value of kernel */ \734 double kernelValue2 = PS_SQR(kernelValue); /* Value of kernel-squared */ \735 if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) { \736 xSumBad += kernelValue2; \737 } else { \738 xSumImage += kernelValue * image->data.TYPE[yPix][xPix]; \739 xSumKernel += kernelValue; \740 xSumKernel2 += kernelValue2; \741 } \742 } \743 float kernelValue = yKernel[j]; /* Value of kernel in y */ \744 double kernelValue2 = PS_SQR(kernelValue); /* Value of kernel-squared in y */ \745 sumImage += kernelValue * xSumImage; \746 sumBad += kernelValue2 * xSumBad; \747 sumKernel += kernelValue * xSumKernel; \748 sumKernel2 += kernelValue2 * xSumKernel2; \749 } \750 } else {\751 /* Neither variance nor mask */ \752 for (int j = jMin, yPix = yMin; j < jMax; j++, yPix++) { \753 double xSumImage = 0.0; /* Sum of image multiplied by kernel in x */ \754 float xSumKernel = 0.0; /* Sum of kernel in x */ \755 for (int i = iMin, xPix = xMin; i < iMax; i++, xPix++) { \756 float kernelValue = xKernel[i]; /* Value of kernel */ \757 xSumImage += kernelValue * image->data.TYPE[yPix][xPix]; \758 xSumKernel += kernelValue; \759 } \760 float kernelValue = yKernel[j]; /* Value of kernel in y */ \761 sumImage += kernelValue * xSumImage; \762 sumKernel += kernelValue * xSumKernel; \763 } \764 } \765 } \766 break;767 768 769 switch (image->type.type) {770 INTERPOLATE_SEPARATE_CASE(U8);771 INTERPOLATE_SEPARATE_CASE(U16);772 INTERPOLATE_SEPARATE_CASE(U32);773 INTERPOLATE_SEPARATE_CASE(U64);774 INTERPOLATE_SEPARATE_CASE(S8);775 INTERPOLATE_SEPARATE_CASE(S16);776 INTERPOLATE_SEPARATE_CASE(S32);777 INTERPOLATE_SEPARATE_CASE(S64);778 INTERPOLATE_SEPARATE_CASE(F32);779 INTERPOLATE_SEPARATE_CASE(F64);780 default:781 psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",782 image->type.type);783 return PS_INTERPOLATE_STATUS_ERROR;784 }785 786 psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_ERROR; // Status of interpolation787 *imageValue = sumImage / sumKernel;788 if (wantVariance) {789 *varianceValue = sumVariance / sumKernel2;790 }791 if (sumBad == 0) {792 // Completely good pixel793 status = PS_INTERPOLATE_STATUS_GOOD;794 } else if (sumBad < PS_SQR(options->poorFrac) * (sumBad + sumKernel2) ) {795 // Some pixels masked: poor pixel796 if (haveMask && maskValue) {797 *maskValue |= options->poorMask;798 }799 status = PS_INTERPOLATE_STATUS_POOR;800 } else {801 // Many pixels (or a few important pixels) masked: bad pixel802 if (haveMask && maskValue) {803 *maskValue |= options->badMask;804 }805 status = PS_INTERPOLATE_STATUS_BAD;806 }807 808 return status;809 }810 811 812 679 813 680 psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue, psMaskType *maskValue, 814 float x, float y, const psImageInterpolat eOptions *options)681 float x, float y, const psImageInterpolation *interp) 815 682 { 816 683 PS_ASSERT_PTR_NON_NULL(imageValue, PS_INTERPOLATE_STATUS_ERROR); 817 PS_ASSERT_PTR_NON_NULL( options, PS_INTERPOLATE_STATUS_ERROR);818 819 const psImage *image = options->image; // Image to interpolate820 const psImage *mask = options->mask; // Mask to interpolate821 const psImage *variance = options->variance; // Variance to interpolate684 PS_ASSERT_PTR_NON_NULL(interp, PS_INTERPOLATE_STATUS_ERROR); 685 686 const psImage *image = interp->image; // Image to interpolate 687 const psImage *mask = interp->mask; // Mask to interpolate 688 const psImage *variance = interp->variance; // Variance to interpolate 822 689 823 690 PS_ASSERT_IMAGE_NON_NULL(image, PS_INTERPOLATE_STATUS_ERROR); … … 834 701 } 835 702 836 PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL( options->poorFrac, 0.0, PS_INTERPOLATE_STATUS_ERROR);837 838 switch ( options->mode) {703 PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(interp->poorFrac, 0.0, PS_INTERPOLATE_STATUS_ERROR); 704 705 switch (interp->mode) { 839 706 case PS_INTERPOLATE_FLAT: 840 return interpolateFlat(imageValue, varianceValue, maskValue, x, y, options);707 return interpolateFlat(imageValue, varianceValue, maskValue, x, y, interp); 841 708 case PS_INTERPOLATE_BILINEAR: 842 case PS_INTERPOLATE_BI CUBE:709 case PS_INTERPOLATE_BIQUADRATIC: 843 710 case PS_INTERPOLATE_GAUSS: 844 return interpolateKernel(imageValue, varianceValue, maskValue, x, y, options);845 711 case PS_INTERPOLATE_LANCZOS2: 846 712 case PS_INTERPOLATE_LANCZOS3: 847 713 case PS_INTERPOLATE_LANCZOS4: 848 return interpolate Separate(imageValue, varianceValue, maskValue, x, y, options);714 return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp); 849 715 default: 850 716 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 851 717 _("Specified interpolation method (%d) is not supported."), 852 options->mode);718 interp->mode); 853 719 return PS_INTERPOLATE_STATUS_ERROR; 854 720 } … … 859 725 860 726 861 static float varianceFactorFlat(float x, float y, const psImageInterpolateOptions *options) 862 { 863 // There's no smearing 864 return 1.0; 865 } 866 867 static float varianceFactorKernel(float x, float y, const psImageInterpolateOptions *options) 868 { 869 int xNum, yNum; // Number of interpolation kernel pixels 870 int xCentral, yCentral; // Central pixel of the convolution 871 interpolateKernelSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode); 872 873 const psImage *image = options->image; // Image of interest 874 int xLast = image->numCols - 1; // Last pixel in x 875 int yLast = image->numRows - 1; // Last pixel in y 876 877 if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast || 878 yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) { 879 // At least one pixel of the interpolation kernel is off the image 880 return NAN; 881 } 882 883 float kernel[yNum][xNum]; // Interpolation kernel for straight interpolation 884 bool xExact, yExact; // Exact shift? 885 // Pretend we're not shifting pixels (i.e., we don't care about any exact shifting) because on the off 886 // chance that a shift is integral, we don't want to get a trivial kernel back, but something that sort of 887 // reflects what's going on. 888 interpolateKernelGenerate(xNum, yNum, kernel, &xExact, &yExact, xCentral, yCentral, 889 x, y, options->mode, false); 890 891 float sumKernel = 0.0; // Sum of kernel 892 double sumKernel2 = 0.0; // Sum of kernel squares 893 for (int j = 0; j < yNum; j++) { 894 for (int i = 0; i < xNum; i++) { 895 float value = kernel[j][i]; // Value of kernel 896 sumKernel += value; 897 sumKernel2 += PS_SQR(value); 898 } 899 } 900 901 return sumKernel2 / PS_SQR(sumKernel); 902 } 903 904 static float varianceFactorSeparate(float x, float y, const psImageInterpolateOptions *options) 905 { 906 int xNum, yNum; // Number of interpolation kernel pixels 907 int xCentral, yCentral; // Central pixel of the convolution 908 interpolateSeparateSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode); 909 910 const psImage *image = options->image; // Image of interest 911 int xLast = image->numCols - 1; // Last pixel in x 912 int yLast = image->numRows - 1; // Last pixel in y 913 914 if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast || 915 yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) { 916 // At least one pixel of the interpolation kernel is off the image 917 return NAN; 918 } 919 920 float xKernel[xNum], yKernel[yNum]; // Interpolation kernels for separable interpolation 921 bool xExact, yExact; // Exact shift? 922 // Pretend we're not shifting pixels (i.e., we don't care about any exact shifting) because on the off 923 // chance that a shift is integral, we don't want to get a trivial kernel back, but something that sort of 924 // reflects what's going on. 925 interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, &xExact, &yExact, 926 xCentral, yCentral, x, y, options->mode, false); 927 928 float ySumKernel = 0.0; // Sum of kernel in y 929 double ySumKernel2 = 0.0; // Sum of kernel squared in y 930 for (int j = 0; j < yNum; j++) { 931 float xSumKernel = 0.0; // Sum of kernel in x 932 double xSumKernel2 = 0.0; // Sum of kernel squared in x 933 for (int i = 0; i < xNum; i++) { 934 float value = xKernel[i]; // Value of kernel 935 xSumKernel += value; 936 xSumKernel2 += PS_SQR(value); 937 } 938 float value = yKernel[j]; // Value of kernel 939 ySumKernel += xSumKernel * value; 940 ySumKernel2 += xSumKernel2 * PS_SQR(value); 941 } 942 943 return ySumKernel2 / PS_SQR(ySumKernel); 944 } 945 946 float psImageInterpolateVarianceFactor(float x, float y, const psImageInterpolateOptions *options) 947 { 948 PS_ASSERT_PTR_NON_NULL(options, PS_INTERPOLATE_STATUS_ERROR); 949 950 switch (options->mode) { 727 static float varianceFactorKernel(float x, float y, psImageInterpolateMode mode) 728 { 729 int size = kernelSizes[mode]; // Size of kernel 730 731 // Kernel basics 732 INTERPOLATE_KERNEL_SETUP(x, y); 733 734 psF32 xKernel[size], yKernel[size]; // Interpolation kernels 735 interpolationKernel(xKernel, xFrac, mode); 736 interpolationKernel(yKernel, yFrac, mode); 737 738 float xSumKernel = 0.0, ySumKernel = 0.0; // Sum of kernel 739 float xSumKernel2 = 0.0, ySumKernel2 = 0.0; // Sum of kernel^2 740 for (int i = 0; i < size; i++) { 741 xSumKernel += xKernel[i]; 742 xSumKernel2 += PS_SQR(xKernel[i]); 743 ySumKernel += yKernel[i]; 744 ySumKernel2 += PS_SQR(yKernel[i]); 745 } 746 747 return (xSumKernel2 * ySumKernel2) / PS_SQR(xSumKernel * ySumKernel); 748 } 749 750 float psImageInterpolateVarianceFactor(float x, float y, psImageInterpolateMode mode) 751 { 752 switch (mode) { 951 753 case PS_INTERPOLATE_FLAT: 952 return varianceFactorFlat(x, y, options); 754 // No smearing by design 755 return 1.0; 953 756 case PS_INTERPOLATE_BILINEAR: 954 case PS_INTERPOLATE_BI CUBE:757 case PS_INTERPOLATE_BIQUADRATIC: 955 758 case PS_INTERPOLATE_GAUSS: 956 return varianceFactorKernel(x, y, options);957 759 case PS_INTERPOLATE_LANCZOS2: 958 760 case PS_INTERPOLATE_LANCZOS3: 959 761 case PS_INTERPOLATE_LANCZOS4: 960 return varianceFactor Separate(x, y, options);762 return varianceFactorKernel(x, y, mode); 961 763 default: 962 764 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 963 765 _("Specified interpolation method (%d) is not supported."), 964 options->mode);766 mode); 965 767 return NAN; 966 768 } … … 977 779 if (!strcasecmp(name, "FLAT")) return PS_INTERPOLATE_FLAT; 978 780 if (!strcasecmp(name, "BILINEAR")) return PS_INTERPOLATE_BILINEAR; 979 if (!strcasecmp(name, "BI CUBE")) return PS_INTERPOLATE_BICUBE;781 if (!strcasecmp(name, "BIQUADRATIC")) return PS_INTERPOLATE_BIQUADRATIC; 980 782 if (!strcasecmp(name, "GAUSS")) return PS_INTERPOLATE_GAUSS; 981 783 if (!strcasecmp(name, "LANCZOS2")) return PS_INTERPOLATE_LANCZOS2; -
trunk/psLib/src/imageops/psImageInterpolate.h
r19140 r20306 7 7 * @author Paul Price, Institute for Astronomy 8 8 * 9 * @version $Revision: 1. 5$ $Name: not supported by cvs2svn $10 * @date $Date: 2008- 08-21 01:12:44$9 * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2008-10-22 02:10:37 $ 11 11 * Copyright 2004-2007 Institute for Astronomy, University of Hawaii 12 12 */ … … 18 18 /// Enumeration of options in interpolation 19 19 typedef enum { 20 PS_INTERPOLATE_NONE , ///< no interpolate defined (error state)21 PS_INTERPOLATE_FLAT, ///< Flat interpolation (nearest pixel)22 PS_INTERPOLATE_BILINEAR, ///< Bilinear interpolation23 PS_INTERPOLATE_BI CUBE, ///< Bicubic interpolation with 3x3 region24 PS_INTERPOLATE_GAUSS, ///< Gaussian inteprolation with 3x3 region25 PS_INTERPOLATE_LANCZOS2, ///< Sinc interpolation with 4x4 pixel kernel26 PS_INTERPOLATE_LANCZOS3, ///< Sinc interpolation with 6x6 pixel kernel27 PS_INTERPOLATE_LANCZOS4, ///< Sinc interpolation with 8x8 pixel kernel20 PS_INTERPOLATE_NONE = 0, ///< No interpolate defined (error state) 21 PS_INTERPOLATE_FLAT, ///< Flat interpolation (nearest pixel) 22 PS_INTERPOLATE_BILINEAR, ///< Bilinear interpolation 23 PS_INTERPOLATE_BIQUADRATIC, ///< Biquadratic interpolation with 3x3 region 24 PS_INTERPOLATE_GAUSS, ///< Gaussian inteprolation with 3x3 region 25 PS_INTERPOLATE_LANCZOS2, ///< Sinc interpolation with 4x4 pixel kernel 26 PS_INTERPOLATE_LANCZOS3, ///< Sinc interpolation with 6x6 pixel kernel 27 PS_INTERPOLATE_LANCZOS4, ///< Sinc interpolation with 8x8 pixel kernel 28 28 } psImageInterpolateMode; 29 29 … … 54 54 float poorFrac; ///< Fraction of flux in bad pixels before output is marked bad 55 55 bool shifting; ///< Shifting images? Don't interpolate if the shift is exact. 56 } psImageInterpolateOptions; 56 int numKernels; ///< Number of pre-calculated kernels 57 const psImage *kernel, *kernel2; ///< 1D interpolation kernel and kernel^2 (row) for each spacing 58 const psVector *sumKernel2; ///< Sum of kernel^2 for each spacing 59 } psImageInterpolation; 57 60 58 61 59 62 /// Allocator 60 psImageInterpolateOptions *psImageInterpolateOptionsAlloc(psImageInterpolateMode mode, // Interpolation mode 61 const psImage *image, // Input image 62 const psImage *variance, // Variance image 63 const psImage *mask, // Mask image 64 psMaskType maskVal, // Value to mask 65 double badImage, // Value for image if bad 66 double badVariance, // Value for variance if bad 67 psMaskType badMask, // Mask value for bad pixels 68 psMaskType poorMask, // Mask value for poor pixels 69 float poorFrac // Fraction of flux for question 63 psImageInterpolation *psImageInterpolationAlloc( 64 psImageInterpolateMode mode, // Interpolation mode 65 const psImage *image, // Input image 66 const psImage *variance, // Variance image 67 const psImage *mask, // Mask image 68 psMaskType maskVal, // Value to mask 69 double badImage, // Value for image if bad 70 double badVariance, // Value for variance if bad 71 psMaskType badMask, // Mask value for bad pixels 72 psMaskType poorMask, // Mask value for poor pixels 73 float poorFrac, // Fraction of flux for question 74 int numKernels // Number of interpolation kernels to pre-calculate 70 75 ) PS_ATTR_MALLOC; 71 76 72 77 73 74 78 /// Interpolate image pixel value given floating point coordinates. 75 psImageInterpolateStatus psImageInterpolate(double *imageValue, ///< Return value for image 76 double *varianceValue, ///< Return value for variance 77 psMaskType *maskValue, ///< Return value for mask 78 float x, float y, ///< Location to which to interpolate 79 const psImageInterpolateOptions *options ///< Options 79 psImageInterpolateStatus psImageInterpolate( 80 double *imageValue, ///< Return value for image 81 double *varianceValue, ///< Return value for variance 82 psMaskType *maskValue, ///< Return value for mask 83 float x, float y, ///< Location to which to interpolate 84 const psImageInterpolation *options ///< Options 80 85 ); 81 82 86 83 87 // Return the appropriate interpolation mode given a char string name for that mode … … 91 95 /// factor. 92 96 float psImageInterpolateVarianceFactor(float x, float y, ///< Position of interest 93 const psImageInterpolateOptions *options ///< Interpolation options97 psImageInterpolateMode mode ///< Interpolation mode 94 98 ); 95 99 -
trunk/psLib/src/imageops/psImagePixelExtract.c
r12741 r20306 8 8 * @author Robert DeSonia, MHPCC 9 9 * 10 * @version $Revision: 1.3 2$ $Name: not supported by cvs2svn $11 * @date $Date: 200 7-04-04 22:42:02$10 * @version $Revision: 1.33 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2008-10-22 02:10:37 $ 12 12 * 13 13 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 679 679 float dY = (endRow - startRow) / (float)(nSamples-1); 680 680 681 psImageInterpolat eOptions *interp = psImageInterpolateOptionsAlloc(mode, input, NULL, mask, maskVal,682 681 psImageInterpolation *interp = psImageInterpolationAlloc(mode, input, NULL, mask, maskVal, 682 0, 0, 0, 0, 0, 0); 683 683 684 684 #define LINEAR_CUT_CASE(TYPE) \ -
trunk/psLib/test/imageops/tap_psImageInterpolate2.c
r19129 r20306 79 79 psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 12345); 80 80 81 psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(mode, image, 82 variance, mask, 0, NAN, NAN, 83 0, 0, 0.0); 81 psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, variance, mask, 0, NAN, NAN, 82 0, 0, 0.0, 0); 84 83 ok(interp, "Interpolation options set"); 85 84 skip_start(!interp, 2 * num, "Interpolation options failed"); … … 93 92 psMaskType maskVal; 94 93 95 psImageInterpolateStatus status = psImageInterpolate(&imageVal, &varianceVal, &maskVal, x, y, interp); 96 ok(status != PS_INTERPOLATE_STATUS_ERROR, "Interpolation at %.2f,%.2f (%x)", x, y, status); 97 skip_start(status == PS_INTERPOLATE_STATUS_ERROR, 1, "Interpolation failed"); 94 psImageInterpolateStatus interpOK = psImageInterpolate(&imageVal, &varianceVal, &maskVal, 95 x, y, interp); 96 ok(interpOK != PS_INTERPOLATE_STATUS_ERROR, "Interpolation at %.2f,%.2f (%x)", x, y, interpOK); 97 skip_start(interpOK == PS_INTERPOLATE_STATUS_ERROR, 1, "Interpolation failed"); 98 98 99 int xCentral, yCentral; // Central pixel of interpolation 100 switch (mode) { 101 case PS_INTERPOLATE_BICUBE: 102 case PS_INTERPOLATE_GAUSS: 103 xCentral = x; 104 yCentral = y; 105 break; 106 default: 107 xCentral = floor(x - 0.5); 108 yCentral = floor(y - 0.5); 109 break; 110 } 99 int xCentral = x - 0.5, yCentral = y - 0.5; // Central pixel of interpolation 111 100 112 101 if (xCentral - (xKernel - 1) / 2 < 0 || xCentral + xKernel / 2 > xSize - 1 || 113 102 yCentral - (yKernel - 1) / 2 < 0 || yCentral + yKernel / 2 > ySize - 1) { 114 ok(status == PS_INTERPOLATE_STATUS_BAD || status == PS_INTERPOLATE_STATUS_POOR, 115 "Interpolation at border"); 103 ok(interpOK == PS_INTERPOLATE_STATUS_OFF || interpOK == PS_INTERPOLATE_STATUS_BAD || 104 interpOK == PS_INTERPOLATE_STATUS_POOR, 105 "Interpolation at %.2f,%.2f (border): %x", x, y, interpOK); 116 106 } else { 117 is_double_tol(imageVal, imageFunc(x, y), tol, "Interpolation = %f vs %f ",118 imageVal, imageFunc(x, y) );107 is_double_tol(imageVal, imageFunc(x, y), tol, "Interpolation = %f vs %f (%d)", 108 imageVal, imageFunc(x, y), interpOK); 119 109 } 120 110 … … 146 136 check(16, 16, PS_TYPE_F64, 32, PS_INTERPOLATE_BILINEAR, 2, 2, 1.0e-6); // 66 tests 147 137 148 check(16, 16, PS_TYPE_F32, 32, PS_INTERPOLATE_BI CUBE, 3, 3, 1.0e-6); // 66 tests149 check(16, 16, PS_TYPE_F64, 32, PS_INTERPOLATE_BI CUBE, 3, 3, 1.0e-6); // 66 tests138 check(16, 16, PS_TYPE_F32, 32, PS_INTERPOLATE_BIQUADRATIC, 3, 3, 1.0e-6); // 66 tests 139 check(16, 16, PS_TYPE_F64, 32, PS_INTERPOLATE_BIQUADRATIC, 3, 3, 1.0e-6); // 66 tests 150 140 151 141 check(16, 16, PS_TYPE_F32, 32, PS_INTERPOLATE_GAUSS, 3, 3, 1.0e-3); // 66 tests
Note:
See TracChangeset
for help on using the changeset viewer.
