IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 21, 2008, 4:10:37 PM (18 years ago)
Author:
Paul Price
Message:

Reworking psImageInterpolate so that kernels are pre-calculated. Renamed a few things in the process: psImageInterpolateOptions --> psImageInterpolation (shorter), PS_INTERPOLATE_BICUBE --> PS_INTERPOLATE_BIQUADRATIC (accurate). I've moved all the interpolation kernels into 1D so that I can pre-calculate them with a reasonable amount of memory. This might be a mistake for the BIQUADRATIC kernel (doesn't have an obvious 1D kernel function), so I might revert just that one soon. Haven't checked if the speed is faster yet, but I checked something in that broke the build, so I need to check this in. Tests pass.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/imageops/psImageInterpolate.c

    r20207 r20306  
    77 *  @author Paul Price, IfA
    88 *
    9  *  @version $Revision: 1.22 $ $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 $
    1111 *
    1212 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    3131#include "psImageInterpolate.h"
    3232
    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'
     38static 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
     49static 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
     59static 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
     77static 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
     91static 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
     123static 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
     138static 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
     171static void imageInterpolationFree(psImageInterpolation *interp)
    34172{
    35173    // 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
     182psImageInterpolation *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
    66235
    67236// Interpolation engine for flat mode (nearest pixel)
    68237static inline psImageInterpolateStatus interpolateFlat(double *imageValue, double *varianceValue,
    69238                                                       psMaskType *maskValue, float x, float y,
    70                                                        const psImageInterpolateOptions *options)
     239                                                       const psImageInterpolation *interp)
    71240{
    72241    // Parameters have been checked by psImageInterpolate()
    73242
    74     const psImage *image = options->image; // Image of interest
     243    const psImage *image = interp->image; // Image of interest
    75244    int xInt = round(x - 0.5 + FLT_EPSILON); // Pixel closest to point of interest in x
    76245    int yInt = round(y - 0.5 + FLT_EPSILON); // Pixel closest to point of interest in y
    … …  
    81250        // At least one pixel of the interpolation kernel is off the image
    82251        if (imageValue) {
    83             *imageValue = options->badImage;
     252            *imageValue = interp->badImage;
    84253        }
    85254        if (varianceValue) {
    86             *varianceValue = options->badVariance;
     255            *varianceValue = interp->badVariance;
    87256        }
    88257        if (maskValue) {
    89             *maskValue = options->badMask;
     258            *maskValue = interp->badMask;
    90259        }
    91260
    … …  
    97266          case PS_TYPE_##TYPE: { \
    98267            if (imageValue) { \
    99                 *imageValue = options->image->data.TYPE[yInt][xInt]; \
     268                *imageValue = interp->image->data.TYPE[yInt][xInt]; \
    100269            } \
    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]; \
    103272            } \
    104273            break; \
    105274        }
    106275
    107         switch (options->image->type.type) {
     276        switch (interp->image->type.type) {
    108277            FLAT_CASE(U8);
    109278            FLAT_CASE(U16);
    … …  
    118287          default:
    119288            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
    120                     options->image->type.type);
     289                    interp->image->type.type);
    121290            return PS_INTERPOLATE_STATUS_ERROR;
    122291        }
    123292
    124293        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];
    127296            } else {
    128297                *maskValue = 0;
    … …  
    133302}
    134303
    135 // Setup for interpolation by kernel
    136 static inline void interpolateKernelSetup(int *xNum, int *yNum, // Size of interpolation kernel, returned
    137                                           int *xCentral, int *yCentral, // Central pixel of convolution
    138                                           float x, float y, // Coordinates of interest
    139                                           psImageInterpolateMode mode // Mode for interpolation
    140                                           )
    141 {
    142     switch (mode) {
    143       case PS_INTERPOLATE_BILINEAR:
    144         *xNum = *yNum = 2;
    145         // Central pixel is the pixel below the point of interest
    146         *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 interest
    153         *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 kernel
    166 // No need to normalise to unity
    167 static inline void interpolateKernelGenerate(int xNum, int yNum, // Size of interpolation kernel
    168                                              float kernel[xNum][yNum], // Kernel, to be set
    169                                              bool *xExact, bool *yExact, // Exact shift?
    170                                              int xCentral, int yCentral, // Central pixel of convolution
    171                                              float x, float y, // Coordinates of interest
    172                                              psImageInterpolateMode mode, // Mode for interpolation
    173                                              bool allowExact // Allow exact shifts?
    174                                              )
    175 {
    176     float xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
    177     float yFrac = y - 0.5 - yCentral; // Fraction of pixel in y
    178     *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 x
    191           float yFrac = y - 0.5 - yCentral; // Fraction of pixel in y
    192           // Calculation variables
    193           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 x
    211           float yFrac = y - yCentral - 0.5; // Fraction of pixel in y
    212           float sigma = 0.5;           // Gaussian sigma
    213           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 }
    230304
    231305// Can't interpolate: kernel is off the image
    232306#define INTERPOLATE_OFF() { \
    233         *imageValue = options->badImage; \
     307        *imageValue = interp->badImage; \
    234308        if (varianceValue) { \
    235             *varianceValue = options->badVariance; \
     309            *varianceValue = interp->badVariance; \
    236310        } \
    237311        if (maskValue) { \
    238             *maskValue = options->badMask; \
     312            *maskValue = interp->badMask; \
    239313        } \
    240314        return PS_INTERPOLATE_STATUS_OFF; \
    241315    }
    242316
    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
    312325static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue,
    313326                                                  psMaskType *maskValue, float x, float y,
    314                                                   const psImageInterpolateOptions *options)
     327                                                  const psImageInterpolation *interp)
    315328{
    316329    // Parameters have been checked by psImageInterpolate()
    317330
    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    }
    332455
    333456    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
    337458    double sumImage = 0.0;              // Sum of image multiplied by kernel
    338459    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) {
    345465    // 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
    352476        if (!yExact) {
    353             // Bottom rows
    354477            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();
    357481                }
     482                INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();
    358483            }
    359484        }
     485        // Two sides
    360486        if (!xExact) {
    361             // Two sides
    362487            for (int j = jMin; j < jMax; j++) {
     488                INTERPOLATE_SEPARATE_SETUP_OFFIMAGE_COL();
    363489                for (int i = 0; i < iMin; i++) {
    364                     INTERPOLATE_KERNEL_ADD_OFFIMAGE();
     490                    INTERPOLATE_SEPARATE_ADD_OFFIMAGE_COL();
    365491                }
    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();
    368494                }
     495                INTERPOLATE_SEPARATE_ADD_OFFIMAGE_ROW();
    369496            }
    370497        }
     498        // Top rows
    371499        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();
    376507                }
    377508            }
    … …  
    384515          if (haveMask) { \
    385516              /* 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; \
    392537                      } 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; \
    397541                      } \
    398542                  } \
     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; \
    399549              } \
    400               \
    401550          } else { \
    402551              /* 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; \
    411570                  } \
     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; \
    412576              } \
    413577          } \
    414578      } else if (haveMask) { \
    415579          /* 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; \
    422597                  } else { \
    423                       sumImage += kernelValue * image->data.TYPE[yPix][xPix]; \
    424                       sumKernel += kernelValue; \
    425                       sumKernel2 += kernelValue2; \
     598                      xSumImage += kernelValue * *imageData; \
     599                      xSumKernel += kernelValue; \
    426600                  } \
    427601              } \
     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; \
    428607          } \
    429       } else { \
     608      } else {\
    430609          /* 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; \
    436621              } \
     622              float kernelValue = *yKernelData; /* Value of kernel in y */ \
     623              sumImage += kernelValue * xSumImage; \
     624              sumKernel += kernelValue * xSumKernel; \
    437625          } \
    438626      } \
    … …  
    466654        // Completely good pixel
    467655        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) ) {
    469657        // Some pixels masked: poor pixel
    470658        if (haveMask && maskValue) {
    471             *maskValue |= options->poorMask;
     659            *maskValue |= interp->poorMask;
    472660        }
    473661        status = PS_INTERPOLATE_STATUS_POOR;
    … …  
    475663        // Many pixels (or a few important pixels) masked: bad pixel
    476664        if (haveMask && maskValue) {
    477             *maskValue |= options->badMask;
     665            *maskValue |= interp->badMask;
    478666        }
    479667        status = PS_INTERPOLATE_STATUS_BAD;
    480668    }
    481669
     670    psFree(xKernelNew);
     671    psFree(yKernelNew);
     672    psFree(xKernel2New);
     673    psFree(yKernel2New);
     674
    482675    return status;
    483676}
    484677
    485678
    486 // Generate Lanczos interpolation kernels
    487 static void lanczos(float values[],    // Interpolation kernel to generate
    488                     bool *exact,        // Exact shift?
    489                     int num,            // Number of values in the kernel
    490                     float frac,         // Sub-pixel position
    491                     bool allowExact     // Allow exact shift?
    492     )
    493 {
    494     if (allowExact && fabs(frac) < DBL_EPSILON) {
    495         *exact = true;
    496         // No real shift
    497         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 laczos
    507         float norm2 = M_PI * 4.0 / (float)num; // Normalisation for sinc function 1
    508         float norm3 = M_PI_2 * 4.0 / (float)num; // Normalisation for sinc function 2
    509         float pos = - (num - 1)/2 - frac;  // Position of interest
    510         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 kernels
    523 static inline void interpolateSeparateSetup(int *xNum, int *yNum, // Size of interpolation kernel, returned
    524                                             int *xCentral, int *yCentral, // Central pixel of convolution
    525                                             float x, float y, // Coordinates of interest
    526                                             psImageInterpolateMode mode // Mode for interpolation
    527                                             )
    528 {
    529     // Central pixel is the pixel below the point of interest
    530     *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 unity
    552 static inline void interpolateSeparateGenerate(int xNum, int yNum, // Size of interpolation kernel
    553                                                float xKernel[xNum], float yKernel[yNum], // Kernels
    554                                                bool *xExact, bool *yExact, // Exact shifts?
    555                                                int xCentral, int yCentral, // Central pixel of convolution
    556                                                float x, float y, // Coordinates of interest
    557                                                psImageInterpolateMode mode, // Mode for interpolation
    558                                                bool allowExact // Allow an exact shift?
    559                                                )
    560 {
    561     // XXX Could put in an "exact shift" (i.e., xFrac = 0.0) version
    562     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 x
    567           lanczos(xKernel, xExact, xNum, xFrac, allowExact);
    568           float yFrac = y - yCentral - 0.5; // Fraction of pixel in y
    569           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 pixels
    589     int xCentral, yCentral; // Central pixel of the convolution
    590     interpolateSeparateSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode);
    591 
    592     const psImage *image = options->image; // Image of interest
    593     const psImage *mask = options->mask; // Image mask
    594     const psImage *variance = options->variance; // Image variance
    595 
    596     INTERPOLATE_BOUNDS();
    597 
    598     float xKernel[xNum], yKernel[yNum]; // Interpolation kernels
    599     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 kernel
    604     double sumKernel2 = 0.0;            // Sum of the kernel-squared
    605     double sumBad = 0.0;                // Sum of bad kernel-squared contributions
    606     psMaskType maskVal = options->maskVal; // Value to mask
    607     double sumImage = 0.0;              // Sum of image multiplied by kernel
    608     double sumVariance = 0.0;           // Sum of variance multiplied by kernel-squared
    609 
    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 image
    616 #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 image
    628     if (offImage) {
    629         // Bottom rows
    630         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 sides
    640         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 rows
    653         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 interpolation
    787     *imageValue = sumImage / sumKernel;
    788     if (wantVariance) {
    789         *varianceValue = sumVariance / sumKernel2;
    790     }
    791     if (sumBad == 0) {
    792         // Completely good pixel
    793         status = PS_INTERPOLATE_STATUS_GOOD;
    794     } else if (sumBad < PS_SQR(options->poorFrac) * (sumBad + sumKernel2) ) {
    795         // Some pixels masked: poor pixel
    796         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 pixel
    802         if (haveMask && maskValue) {
    803             *maskValue |= options->badMask;
    804         }
    805         status = PS_INTERPOLATE_STATUS_BAD;
    806     }
    807 
    808     return status;
    809 }
    810 
    811 
    812679
    813680psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue, psMaskType *maskValue,
    814                                             float x, float y, const psImageInterpolateOptions *options)
     681                                            float x, float y, const psImageInterpolation *interp)
    815682{
    816683    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 interpolate
    820     const psImage *mask = options->mask; // Mask to interpolate
    821     const psImage *variance = options->variance; // Variance to interpolate
     684    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
    822689
    823690    PS_ASSERT_IMAGE_NON_NULL(image, PS_INTERPOLATE_STATUS_ERROR);
    … …  
    834701    }
    835702
    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) {
    839706      case PS_INTERPOLATE_FLAT:
    840         return interpolateFlat(imageValue, varianceValue, maskValue, x, y, options);
     707        return interpolateFlat(imageValue, varianceValue, maskValue, x, y, interp);
    841708      case PS_INTERPOLATE_BILINEAR:
    842       case PS_INTERPOLATE_BICUBE:
     709      case PS_INTERPOLATE_BIQUADRATIC:
    843710      case PS_INTERPOLATE_GAUSS:
    844         return interpolateKernel(imageValue, varianceValue, maskValue, x, y, options);
    845711      case PS_INTERPOLATE_LANCZOS2:
    846712      case PS_INTERPOLATE_LANCZOS3:
    847713      case PS_INTERPOLATE_LANCZOS4:
    848         return interpolateSeparate(imageValue, varianceValue, maskValue, x, y, options);
     714        return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp);
    849715      default:
    850716        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    851717                _("Specified interpolation method (%d) is not supported."),
    852                 options->mode);
     718                interp->mode);
    853719        return PS_INTERPOLATE_STATUS_ERROR;
    854720    }
    … …  
    859725
    860726
    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) {
     727static 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
     750float psImageInterpolateVarianceFactor(float x, float y, psImageInterpolateMode mode)
     751{
     752    switch (mode) {
    951753      case PS_INTERPOLATE_FLAT:
    952         return varianceFactorFlat(x, y, options);
     754        // No smearing by design
     755        return 1.0;
    953756      case PS_INTERPOLATE_BILINEAR:
    954       case PS_INTERPOLATE_BICUBE:
     757      case PS_INTERPOLATE_BIQUADRATIC:
    955758      case PS_INTERPOLATE_GAUSS:
    956         return varianceFactorKernel(x, y, options);
    957759      case PS_INTERPOLATE_LANCZOS2:
    958760      case PS_INTERPOLATE_LANCZOS3:
    959761      case PS_INTERPOLATE_LANCZOS4:
    960         return varianceFactorSeparate(x, y, options);
     762        return varianceFactorKernel(x, y, mode);
    961763      default:
    962764        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    963765                _("Specified interpolation method (%d) is not supported."),
    964                 options->mode);
     766                mode);
    965767        return NAN;
    966768    }
    … …  
    977779    if (!strcasecmp(name, "FLAT"))     return PS_INTERPOLATE_FLAT;
    978780    if (!strcasecmp(name, "BILINEAR")) return PS_INTERPOLATE_BILINEAR;
    979     if (!strcasecmp(name, "BICUBE"))   return PS_INTERPOLATE_BICUBE;
     781    if (!strcasecmp(name, "BIQUADRATIC")) return PS_INTERPOLATE_BIQUADRATIC;
    980782    if (!strcasecmp(name, "GAUSS"))    return PS_INTERPOLATE_GAUSS;
    981783    if (!strcasecmp(name, "LANCZOS2")) return PS_INTERPOLATE_LANCZOS2;
Note: See TracChangeset for help on using the changeset viewer.