/** @file  psImageInterpolate.c
 *
 *  @brief Contains functions for interpolating an image
 *
 *  @author Robert DeSonia, MHPCC
 *  @author Ross Harman, MHPCC
 *  @author Paul Price, IfA
 *
 *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2008-06-24 02:04:55 $
 *
 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
 *
 */

#ifdef HAVE_CONFIG_H
# include "config.h"
#endif

#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <string.h>

#include "psAbort.h"
#include "psMemory.h"
#include "psError.h"
#include "psAssert.h"
#include "psString.h"
#include "psImage.h"
#include "psImageInterpolate.h"

static void imageInterpolateOptionsFree(psImageInterpolateOptions *options)
{
    // Casting away const
    psFree((psImage*)options->image);
    psFree((psImage*)options->mask);
    psFree((psImage*)options->variance);
}

psImageInterpolateOptions *psImageInterpolateOptionsAlloc(psImageInterpolateMode mode,
                                                          const psImage *image, const psImage *variance,
                                                          const psImage *mask, psMaskType maskVal,
                                                          double badImage, double badVariance,
                                                          psMaskType badMask, psMaskType poorMask,
                                                          float poorFrac)
{
    psImageInterpolateOptions *options = psAlloc(sizeof(psImageInterpolateOptions)); // Options, to return
    psMemSetDeallocator(options, (psFreeFunc)imageInterpolateOptionsFree);

    options->mode = mode;
    // Casting away const to add to options
    options->image = psMemIncrRefCounter((psImage*)image);
    options->variance = psMemIncrRefCounter((psImage*)variance);
    options->mask = psMemIncrRefCounter((psImage*)mask);
    options->maskVal = maskVal;
    options->badImage = badImage;
    options->badVariance = badVariance;
    options->badMask = badMask;
    options->poorMask = poorMask;
    options->poorFrac = poorFrac;

    return options;
}

// Interpolation engine for flat mode (nearest pixel)
static inline psImageInterpolateStatus interpolateFlat(double *imageValue, double *varianceValue,
                                                       psMaskType *maskValue, float x, float y,
                                                       const psImageInterpolateOptions *options)
{
    // Parameters have been checked by psImageInterpolate()

    const psImage *image = options->image; // Image of interest
    int xInt = round(x - 0.5 + FLT_EPSILON); // Pixel closest to point of interest in x
    int yInt = round(y - 0.5 + FLT_EPSILON); // Pixel closest to point of interest in y
    int xLast = image->numCols - 1;     // Last pixel in x
    int yLast = image->numRows - 1;     // Last pixel in y

    if (xInt < 0 || xInt > xLast || yInt < 0 || yInt > yLast) {
        // At least one pixel of the interpolation kernel is off the image
        if (imageValue) {
            *imageValue = options->badImage;
        }
        if (varianceValue) {
            *varianceValue = options->badVariance;
        }
        if (maskValue) {
            *maskValue = options->badMask;
        }

        return PS_INTERPOLATE_STATUS_OFF;
    } else {

        // Image and variance 'interpolation' according to image type
        #define FLAT_CASE(TYPE) \
          case PS_TYPE_##TYPE: { \
            if (imageValue) { \
                *imageValue = options->image->data.TYPE[yInt][xInt]; \
            } \
            if (varianceValue && options->variance) { \
                *varianceValue = options->variance->data.TYPE[yInt][xInt]; \
            } \
            break; \
        }

        switch (options->image->type.type) {
            FLAT_CASE(U8);
            FLAT_CASE(U16);
            FLAT_CASE(U32);
            FLAT_CASE(U64);
            FLAT_CASE(S8);
            FLAT_CASE(S16);
            FLAT_CASE(S32);
            FLAT_CASE(S64);
            FLAT_CASE(F32);
            FLAT_CASE(F64);
          default:
            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
                    options->image->type.type);
            return PS_INTERPOLATE_STATUS_ERROR;
        }

        if (maskValue) {
            if (options->mask) {
                *maskValue = options->mask->data.PS_TYPE_MASK_DATA[yInt][xInt];
            } else {
                *maskValue = 0;
            }
        }
    }
    return PS_INTERPOLATE_STATUS_GOOD;
}

// Setup for interpolation by kernel
static inline void interpolateKernelSetup(int *xNum, int *yNum, // Size of interpolation kernel, returned
                                          int *xCentral, int *yCentral, // Central pixel of convolution
                                          float x, float y, // Coordinates of interest
                                          psImageInterpolateMode mode // Mode for interpolation
                                          )
{
    switch (mode) {
      case PS_INTERPOLATE_BILINEAR:
        *xNum = *yNum = 2;
        // Central pixel is the pixel below the point of interest
        *xCentral = floor(x - 0.5 + FLT_EPSILON);
        *yCentral = floor(y - 0.5 + FLT_EPSILON);
        break;
      case PS_INTERPOLATE_BICUBE:
      case PS_INTERPOLATE_GAUSS:
        *xNum = *yNum = 3;
        // Central pixel is the closest pixel to the point of interest
        *xCentral = x;
        *yCentral = y;
        break;
      case PS_INTERPOLATE_FLAT:
      case PS_INTERPOLATE_LANCZOS2:
      case PS_INTERPOLATE_LANCZOS3:
      case PS_INTERPOLATE_LANCZOS4:
      default:
        psAbort("Invalid interpolation mode.");
    }
}

// Generate the interpolation kernel; it should be normalised to unity
static inline void interpolateKernelGenerate(int xNum, int yNum, // Size of interpolation kernel
                                             double kernel[xNum][yNum], // Kernel, to be set
                                             int xCentral, int yCentral, // Central pixel of convolution
                                             float x, float y, // Coordinates of interest
                                             psImageInterpolateMode mode // Mode for interpolation
                                             )
{
    switch (mode) {
      case PS_INTERPOLATE_BILINEAR: {
          double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
          double yFrac = y - 0.5 - yCentral; // Fraction of pixel in y
          kernel[0][0] = (1.0 - xFrac) * (1.0 - yFrac);
          kernel[0][1] = xFrac * (1.0 - yFrac);
          kernel[1][0] = yFrac * (1.0 - xFrac);
          kernel[1][1] = xFrac * yFrac;
          break;
      }
      case PS_INTERPOLATE_BICUBE: {
          double xFrac = x - 0.5 - xCentral; // Fraction of pixel in x
          double yFrac = y - 0.5 - yCentral; // Fraction of pixel in y
          // Calculation variables
          double xxFrac = xFrac * xFrac / 6.0;
          double yyFrac = yFrac * yFrac / 6.0;
          double xyFrac = 0.25 * xFrac * yFrac;
          xFrac /= 6.0;
          yFrac /= 6.0;
          kernel[0][0] = - 1.0/9.0 - xFrac - yFrac + xxFrac + yyFrac + xyFrac;
          kernel[0][1] = 2.0/9.0 - yFrac - 2.0 * xxFrac + yyFrac;
          kernel[0][2] = - 1.0/9.0 + xFrac - yFrac + xxFrac + yyFrac - xyFrac;
          kernel[1][0] = 2.0/9.0 - xFrac + xxFrac - 2.0 * yyFrac;
          kernel[1][1] = 5.0/9.0 - 2.0 * xxFrac - 2.0 * yyFrac;
          kernel[1][2] = 2.0/9.0 + xFrac + xxFrac - 2.0 * yyFrac;
          kernel[2][0] = - 1.0/9.0 - xFrac + yFrac + xxFrac + yyFrac - xyFrac;
          kernel[2][1] = 2.0/9.0 + yFrac - 2.0 * xxFrac + yyFrac;
          kernel[2][2] = - 1.0/9.0 + xFrac + yFrac + xxFrac + yyFrac + xyFrac;
          break;
      }
      case PS_INTERPOLATE_GAUSS: {
          double xFrac = x - xCentral - 0.5; // Fraction of pixel in x
          double yFrac = y - yCentral - 0.5; // Fraction of pixel in y
          double sigma = 0.5;           // Gaussian sigma
          double norm = 0.0;            // Normalisation
          for (int j = 0, yPos = - (yNum - 1) / 2; j < yNum; j++, yPos++) {
              for (int i = 0, xPos = - (xNum - 1) / 2; i < xNum; i++, xPos++) {
                  norm += kernel[j][i] = exp(-0.5 / PS_SQR(sigma) * (PS_SQR(xPos - xFrac) +
                                                                     PS_SQR(yPos - yFrac)));
              }
          }
          norm = 1.0 / norm;
          for (int j = 0; j < yNum; j++) {
              for (int i = 0; i < xNum; i++) {
                  kernel[j][i] *= norm;
              }
          }
          break;
      }
      case PS_INTERPOLATE_FLAT:
      case PS_INTERPOLATE_LANCZOS2:
      case PS_INTERPOLATE_LANCZOS3:
      case PS_INTERPOLATE_LANCZOS4:
      default:
        psAbort("Invalid interpolation mode.");
    }
}

// Interpolation engine using interpolation kernel
static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue,
                                                  psMaskType *maskValue, float x, float y,
                                                  const psImageInterpolateOptions *options)
{
    // Parameters have been checked by psImageInterpolate()

    int xNum, yNum;                     // Number of interpolation kernel pixels
    int xCentral, yCentral;             // Central pixel of the convolution
    interpolateKernelSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode);

    const psImage *image = options->image; // Image of interest
    int xLast = image->numCols - 1;     // Last pixel in x
    int yLast = image->numRows - 1;     // Last pixel in y

    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
        // At least one pixel of the interpolation kernel is off the image
        if (imageValue) {
            *imageValue = options->badImage;
        }
        if (varianceValue) {
            *varianceValue = options->badVariance;
        }
        if (maskValue) {
            *maskValue = options->badMask;
        }
        return PS_INTERPOLATE_STATUS_OFF;
    }

    double kernel[yNum][xNum];          // Interpolation kernel for straight interpolation
    interpolateKernelGenerate(xNum, yNum, kernel, xCentral, yCentral, x, y, options->mode);

    // Image interpolation, according to image type
    #define KERNEL_IMAGE_CASE(TYPE) \
        case PS_TYPE_##TYPE: { \
            for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
                for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
                    *imageValue += kernel[j][i] * image->data.TYPE[yPix][xPix]; \
                } \
            } \
            break; \
        }

    // Calculate the value for the image
    if (imageValue) {
        *imageValue = 0.0;
        switch (image->type.type) {
            KERNEL_IMAGE_CASE(U8);
            KERNEL_IMAGE_CASE(U16);
            KERNEL_IMAGE_CASE(U32);
            KERNEL_IMAGE_CASE(U64);
            KERNEL_IMAGE_CASE(S8);
            KERNEL_IMAGE_CASE(S16);
            KERNEL_IMAGE_CASE(S32);
            KERNEL_IMAGE_CASE(S64);
            KERNEL_IMAGE_CASE(F32);
            KERNEL_IMAGE_CASE(F64);
          default:
            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
                    image->type.type);
            return PS_INTERPOLATE_STATUS_ERROR;
        }
    }

    // Check the mask value
    const psImage *mask = options->mask; // Image mask
    psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_GOOD; // Status of interpolation
    if (maskValue) {
        *maskValue = 0;
        if (mask) {
            psMaskType maskVal = options->maskVal; // Mask value
            double badContrib = 0.0;        // Amount of kernel on bad pixels
            double totContrib = 0.0;        // Total amount of kernel
            bool badPix = false;            // Is there a bad pixel?
            for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) {
                for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) {
                    if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) {
                        badContrib += fabs(kernel[j][i]);
                        badPix = true;
                    }
                    *maskValue |= mask->data.PS_TYPE_MASK_DATA[yPix][xPix];
                    totContrib += fabs(kernel[j][i]);
                }
            }

            if (badPix) {
                if (badContrib / totContrib >= options->poorFrac) {
                    *maskValue |= options->badMask;
                    status = PS_INTERPOLATE_STATUS_BAD;
                } else {
                    *maskValue |= options->poorMask;
                    status = PS_INTERPOLATE_STATUS_POOR;
                }
            }
        }
    }

    // Finally, the variance
    const psImage *variance = options->variance; // Image variance
    if (varianceValue && variance) {
        *varianceValue = 0.0;

        // Variance interpolation, according to image type
        #define KERNEL_VARIANCE_CASE(TYPE) \
            case PS_TYPE_##TYPE: { \
                double sumKernel2 = 0.0; /* Sum of kernel squares */ \
                for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
                    for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
                        double kernel2 = PS_SQR(kernel[j][i]); /* Kernel squared */ \
                        sumKernel2 += kernel2; \
                        *varianceValue += kernel2 * variance->data.TYPE[yPix][xPix]; \
                    } \
                } \
                *varianceValue /= sumKernel2; /* Normalise so that sum of kernel squares is unity */ \
                break; \
            }

        switch (variance->type.type) {
            KERNEL_VARIANCE_CASE(U8);
            KERNEL_VARIANCE_CASE(U16);
            KERNEL_VARIANCE_CASE(U32);
            KERNEL_VARIANCE_CASE(U64);
            KERNEL_VARIANCE_CASE(S8);
            KERNEL_VARIANCE_CASE(S16);
            KERNEL_VARIANCE_CASE(S32);
            KERNEL_VARIANCE_CASE(S64);
            KERNEL_VARIANCE_CASE(F32);
            KERNEL_VARIANCE_CASE(F64);
          default:
            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
                    image->type.type);
            return PS_INTERPOLATE_STATUS_ERROR;
        }
    }

    return status;
}


// Generate Lanczos interpolation kernels
static void lanczos(double values[],    // Interpolation kernel to generate
                    int num,            // Number of values in the kernel
                    float frac          // Sub-pixel position
    )
{
    // XXX: Instead of generating a convolution kernel that does no shifting, need to set a boolean that says
    // we can do an exact shift.
    if (fabs(frac) < DBL_EPSILON) {
        // No real shift
        for (int i = 0; i < (num - 1) / 2; i++) {
            values[i] = 0.0;
        }
        values[(num - 1) / 2] = 1.0;
        for (int i = (num - 1) / 2 + 1; i < num; i++) {
            values[i] = 0.0;
        }
    } else {
        double norm4 = 0.0;             // Normalisation
        double norm1 = 2.0 / PS_SQR(M_PI); // Normalisation for laczos
        double norm2 = M_PI * 4.0 / (float)num; // Normalisation for sinc function 1
        double norm3 =M_PI_2 * 4.0 / (float)num; // Normalisation for sinc function 2
        double pos = - (num - 1)/2 - frac;  // Position of interest
        for (int i = 0; i < num; i++, pos += 1.0) {
            norm4 += values[i] = norm1 * sin(pos * norm2) * sin(pos * norm3) / PS_SQR(pos);
        }
        norm4 = 1.0 / norm4;
        for (int i = 0; i < num; i++) {
            values[i] *= norm4;
        }
    }

    return;
}

// Setup for interpolation by separable kernels
static inline void interpolateSeparateSetup(int *xNum, int *yNum, // Size of interpolation kernel, returned
                                            int *xCentral, int *yCentral, // Central pixel of convolution
                                            float x, float y, // Coordinates of interest
                                            psImageInterpolateMode mode // Mode for interpolation
                                            )
{
    // Central pixel is the pixel below the point of interest
    *xCentral = floor(x - 0.5);
    *yCentral = floor(y - 0.5);
    switch (mode) {
      case PS_INTERPOLATE_LANCZOS2:
        *xNum = *yNum = 4;
        break;
      case PS_INTERPOLATE_LANCZOS3:
        *xNum = *yNum = 6;
        break;
      case PS_INTERPOLATE_LANCZOS4:
        *xNum = *yNum = 8;
        break;
      case PS_INTERPOLATE_FLAT:
      case PS_INTERPOLATE_BILINEAR:
      case PS_INTERPOLATE_BICUBE:
      case PS_INTERPOLATE_GAUSS:
      default:
        psAbort("Invalid interpolation mode.");
    }
}

// Generate the interpolation kernels for separable case; they should be normalised to unity
static inline void interpolateSeparateGenerate(int xNum, int yNum, // Size of interpolation kernel
                                               double xKernel[xNum], double yKernel[yNum], // Kernels
                                               int xCentral, int yCentral, // Central pixel of convolution
                                               float x, float y, // Coordinates of interest
                                               psImageInterpolateMode mode // Mode for interpolation
                                               )
{
    // XXX Could put in an "exact shift" (i.e., xFrac = 0.0) version
    switch (mode) {
      case PS_INTERPOLATE_LANCZOS2:
      case PS_INTERPOLATE_LANCZOS3:
      case PS_INTERPOLATE_LANCZOS4: {
          double xFrac = x - xCentral - 0.5; // Fraction of pixel in x
          lanczos(xKernel, xNum, xFrac);
          double yFrac = y - yCentral - 0.5; // Fraction of pixel in y
          lanczos(yKernel, yNum, yFrac);
          break;
      }
      case PS_INTERPOLATE_FLAT:
      case PS_INTERPOLATE_BILINEAR:
      case PS_INTERPOLATE_BICUBE:
      case PS_INTERPOLATE_GAUSS:
      default:
        psAbort("Invalid interpolation mode.");
    }
}

// Interpolation engine for separable interpolation kernels (either for good reasons or for practical reasons)
static psImageInterpolateStatus interpolateSeparate(double *imageValue, double *varianceValue,
                                                    psMaskType *maskValue, float x, float y,
                                                    const psImageInterpolateOptions *options)
{
    // Parameters have been checked by psImageInterpolate()

    int xNum, yNum;                     // Number of interpolation kernel pixels
    int xCentral, yCentral; // Central pixel of the convolution
    interpolateSeparateSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode);

    const psImage *image = options->image; // Image of interest
    int xLast = image->numCols - 1;     // Last pixel in x
    int yLast = image->numRows - 1;     // Last pixel in y

    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
        // At least one pixel of the interpolation kernel is off the image
        if (imageValue) {
            *imageValue = options->badImage;
        }
        if (varianceValue) {
            *varianceValue = options->badVariance;
        }
        if (maskValue) {
            *maskValue = options->badMask;
        }
        return PS_INTERPOLATE_STATUS_OFF;
    }

    double xKernel[xNum], yKernel[yNum]; // Interpolation kernels
    interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, xCentral, yCentral, x, y, options->mode);

    // Image interpolation, according to image type
    #define SEPARATE_IMAGE_CASE(TYPE) \
      case PS_TYPE_##TYPE: { \
        for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
            double xInterpValue = 0.0; /* Interpolation in x */ \
            for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
                xInterpValue += xKernel[i] * image->data.TYPE[yPix][xPix]; \
            } \
            *imageValue += yKernel[j] * xInterpValue; /* Interpolating in y */ \
        } \
        break; \
      }

    // Calculate the value for the image
    if (imageValue) {
        *imageValue = 0.0;
        switch (image->type.type) {
            SEPARATE_IMAGE_CASE(U8);
            SEPARATE_IMAGE_CASE(U16);
            SEPARATE_IMAGE_CASE(U32);
            SEPARATE_IMAGE_CASE(U64);
            SEPARATE_IMAGE_CASE(S8);
            SEPARATE_IMAGE_CASE(S16);
            SEPARATE_IMAGE_CASE(S32);
            SEPARATE_IMAGE_CASE(S64);
            SEPARATE_IMAGE_CASE(F32);
            SEPARATE_IMAGE_CASE(F64);
          default:
            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
                    image->type.type);
            return PS_INTERPOLATE_STATUS_ERROR;
        }
    }

    // Check the mask value
    const psImage *mask = options->mask; // Image mask
    psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_GOOD; // Status of interpolation
    if (maskValue) {
        *maskValue = 0;
        if (mask) {
            psMaskType maskVal = options->maskVal; // Mask value
            double badContrib = 0.0;        // Amount of kernel on bad pixels
            double totContrib = 0.0;        // Total amount of kernel
            bool badPix = false;            // Number of bad pixels
            for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) {
                // Interpolation in x
                double xBadContrib = 0.0;   // Amount of x kernel on bad pixels
                double xTotContrib = 0.0;   // Total amount of x kernel
                for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) {
                    if (mask->data.PS_TYPE_MASK_DATA[yPix][xPix] & maskVal) {
                        xBadContrib += fabs(xKernel[i]);
                        badPix = true;
                    }
                    *maskValue |= mask->data.PS_TYPE_MASK_DATA[yPix][xPix];
                    xTotContrib += fabs(xKernel[i]);
                }
                // Interpolating in y
                badContrib += fabs(yKernel[j]) * xBadContrib;
                totContrib += fabs(yKernel[j]) * xTotContrib;
            }

            if (badPix) {
                if (badContrib / totContrib >= options->poorFrac) {
                    *maskValue |= options->badMask;
                    status = PS_INTERPOLATE_STATUS_BAD;
                } else {
                    *maskValue |= options->poorMask;
                    status = PS_INTERPOLATE_STATUS_POOR;
                }
            }
        }
    }

    // Finally, the variance
    const psImage *variance = options->variance; // Image variance
    if (varianceValue && variance) {
        *varianceValue = 0.0;

        // Variance interpolation, according to image type
        #define SEPARATE_VARIANCE_CASE(TYPE) \
          case PS_TYPE_##TYPE: { \
            double ySumKernel2 = 0.0; /* Sum of kernel squared in y */ \
            for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
                double xSumKernel2 = 0.0; /* Sum of kernel squared in x */ \
                double xInterpValue = 0.0; /* Interpolation in x */ \
                for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
                    double kernel2 = PS_SQR(xKernel[i]); /* Kernel squared */ \
                    xSumKernel2 += kernel2; \
                    xInterpValue += kernel2 * variance->data.TYPE[yPix][xPix]; \
                } \
                double kernel2 = PS_SQR(yKernel[j]); /* Kernel squared */ \
                ySumKernel2 += xSumKernel2 * kernel2; \
                *varianceValue += xInterpValue * kernel2; /* Interpolating in y */ \
            } \
            *varianceValue /= ySumKernel2; \
            break; \
          }

        switch (variance->type.type) {
            SEPARATE_VARIANCE_CASE(U8);
            SEPARATE_VARIANCE_CASE(U16);
            SEPARATE_VARIANCE_CASE(U32);
            SEPARATE_VARIANCE_CASE(U64);
            SEPARATE_VARIANCE_CASE(S8);
            SEPARATE_VARIANCE_CASE(S16);
            SEPARATE_VARIANCE_CASE(S32);
            SEPARATE_VARIANCE_CASE(S64);
            SEPARATE_VARIANCE_CASE(F32);
            SEPARATE_VARIANCE_CASE(F64);
          default:
            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unrecognised type for image: %x",
                    variance->type.type);
            return PS_INTERPOLATE_STATUS_ERROR;
        }
    }

    return status;
}



psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue, psMaskType *maskValue,
                                            float x, float y, const psImageInterpolateOptions *options)
{
    PS_ASSERT_PTR_NON_NULL(options, PS_INTERPOLATE_STATUS_ERROR);

    const psImage *image = options->image; // Image to interpolate
    const psImage *mask = options->mask; // Mask to interpolate
    const psImage *variance = options->variance; // Variance to interpolate

    PS_ASSERT_IMAGE_NON_NULL(image, PS_INTERPOLATE_STATUS_ERROR);
    if (varianceValue && variance) {
        PS_ASSERT_IMAGE_NON_NULL(variance, PS_INTERPOLATE_STATUS_ERROR);
        PS_ASSERT_IMAGE_TYPE(variance, image->type.type, PS_INTERPOLATE_STATUS_ERROR);
        psAssert(image->numCols == variance->numCols && image->numRows == variance->numRows,
                 "Image and variance sizes");
    }
    if (maskValue && mask) {
        PS_ASSERT_IMAGE_NON_NULL(mask, PS_INTERPOLATE_STATUS_ERROR);
        PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, PS_INTERPOLATE_STATUS_ERROR);
        psAssert(image->numCols == mask->numCols && image->numRows == mask->numRows, "Image and mask sizes");
    }

    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(options->poorFrac, 0.0, PS_INTERPOLATE_STATUS_ERROR);

    switch (options->mode) {
      case PS_INTERPOLATE_FLAT:
        return interpolateFlat(imageValue, varianceValue, maskValue, x, y, options);
      case PS_INTERPOLATE_BILINEAR:
      case PS_INTERPOLATE_BICUBE:
      case PS_INTERPOLATE_GAUSS:
        return interpolateKernel(imageValue, varianceValue, maskValue, x, y, options);
      case PS_INTERPOLATE_LANCZOS2:
      case PS_INTERPOLATE_LANCZOS3:
      case PS_INTERPOLATE_LANCZOS4:
        return interpolateSeparate(imageValue, varianceValue, maskValue, x, y, options);
      default:
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                _("Specified interpolation method (%d) is not supported."),
                options->mode);
        return PS_INTERPOLATE_STATUS_ERROR;
    }

    psAbort("Should never reach here.");
    return PS_INTERPOLATE_STATUS_ERROR;
}


static float varianceFactorFlat(float x, float y, const psImageInterpolateOptions *options)
{
    // There's no smearing
    return 1.0;
}

static float varianceFactorKernel(float x, float y, const psImageInterpolateOptions *options)
{
    int xNum, yNum;                     // Number of interpolation kernel pixels
    int xCentral, yCentral;             // Central pixel of the convolution
    interpolateKernelSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode);

    const psImage *image = options->image; // Image of interest
    int xLast = image->numCols - 1;     // Last pixel in x
    int yLast = image->numRows - 1;     // Last pixel in y

    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
        // At least one pixel of the interpolation kernel is off the image
        return NAN;
    }

    double kernel[yNum][xNum];          // Interpolation kernel for straight interpolation
    interpolateKernelGenerate(xNum, yNum, kernel, xCentral, yCentral, x, y, options->mode);

    double sumKernel2 = 0.0;            // Sum of kernel squares
    for (int j = 0; j < yNum; j++) {
        for (int i = 0; i < xNum; i++) {
            sumKernel2 += PS_SQR(kernel[j][i]);
        }
    }

    return sumKernel2;
}

static float varianceFactorSeparate(float x, float y, const psImageInterpolateOptions *options)
{
    int xNum, yNum;                     // Number of interpolation kernel pixels
    int xCentral, yCentral;             // Central pixel of the convolution
    interpolateSeparateSetup(&xNum, &yNum, &xCentral, &yCentral, x, y, options->mode);

    const psImage *image = options->image; // Image of interest
    int xLast = image->numCols - 1;     // Last pixel in x
    int yLast = image->numRows - 1;     // Last pixel in y

    if (xCentral - (xNum - 1) / 2 < 0 || xCentral + xNum / 2 > xLast ||
        yCentral - (yNum - 1) / 2 < 0 || yCentral + yNum / 2 > yLast) {
        // At least one pixel of the interpolation kernel is off the image
        return NAN;
    }

    double xKernel[xNum], yKernel[yNum]; // Interpolation kernels for separable interpolation
    interpolateSeparateGenerate(xNum, yNum, xKernel, yKernel, xCentral, yCentral, x, y, options->mode);

    double ySumKernel2 = 0.0;           // Sum of kernel squared in y
    for (int j = 0; j < yNum; j++) {
        double xSumKernel2 = 0.0;       // Sum of kernel squared in x
        for (int i = 0; i < xNum; i++) {
            xSumKernel2 += PS_SQR(xKernel[i]);
        }
        ySumKernel2 += xSumKernel2 * PS_SQR(yKernel[j]);
    }

    return ySumKernel2;
}

float psImageInterpolateVarianceFactor(float x, float y, const psImageInterpolateOptions *options)
{
    PS_ASSERT_PTR_NON_NULL(options, PS_INTERPOLATE_STATUS_ERROR);

    switch (options->mode) {
      case PS_INTERPOLATE_FLAT:
        return varianceFactorFlat(x, y, options);
      case PS_INTERPOLATE_BILINEAR:
      case PS_INTERPOLATE_BICUBE:
      case PS_INTERPOLATE_GAUSS:
        return varianceFactorKernel(x, y, options);
      case PS_INTERPOLATE_LANCZOS2:
      case PS_INTERPOLATE_LANCZOS3:
      case PS_INTERPOLATE_LANCZOS4:
        return varianceFactorSeparate(x, y, options);
      default:
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                _("Specified interpolation method (%d) is not supported."),
                options->mode);
        return NAN;
    }

    psAbort("Should never reach here.");
    return NAN;
}


psImageInterpolateMode psImageInterpolateModeFromString(const char *name)
{
    PS_ASSERT_STRING_NON_EMPTY(name, PS_INTERPOLATE_NONE);

    if (!strcasecmp(name, "FLAT"))     return PS_INTERPOLATE_FLAT;
    if (!strcasecmp(name, "BILINEAR")) return PS_INTERPOLATE_BILINEAR;
    if (!strcasecmp(name, "BICUBE"))   return PS_INTERPOLATE_BICUBE;
    if (!strcasecmp(name, "GAUSS"))    return PS_INTERPOLATE_GAUSS;
    if (!strcasecmp(name, "LANCZOS2")) return PS_INTERPOLATE_LANCZOS2;
    if (!strcasecmp(name, "LANCZOS3")) return PS_INTERPOLATE_LANCZOS3;
    if (!strcasecmp(name, "LANCZOS4")) return PS_INTERPOLATE_LANCZOS4;

    psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Unknown interpolation type %s"), name);
    return PS_INTERPOLATE_NONE;
}
