/** @file pmSubtraction.c
 *
 *  @author Paul Price, IfA
 *  @author GLG, MHPCC
 *
 *  @version $Revision: 1.62 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2007-09-17 21:40:22 $
 *
 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <pslib.h>

#include "pmHDU.h"                      // Required for pmFPA.h
#include "pmFPA.h"
#include "pmSubtractionStamps.h"

#include "pmSubtraction.h"

//#define TESTING

#define PIXEL_LIST_BUFFER 100           // Number of entries to add to pixel list at a time


//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Private (file-static) functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

// Return the number of polynomial terms there are for a given order
static inline long polyTerms(int order  // Order of polynomial
                             )
{
    return (order + 1) * (order + 2) / 2;
}

// Given a stamp coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j
static psImage *spatialPolyValues(int spatialOrder, // Maximum spatial polynomial order
                                  float x, float y // Normalised position of interest, [-1,1]
    )
{
    assert(spatialOrder >= 0);
    assert(x >= -1 && x <= 1);
    assert(y >= -1 && y <= 1);

    psImage *polyValues = psImageAlloc(spatialOrder + 1, spatialOrder + 1, PS_TYPE_F64); // Matrix to return
    polyValues->data.F64[0][0] = 1.0;

    double value = 1.0;
    for (int i = 1; i <= spatialOrder; i++) {
        value *= x;
        polyValues->data.F64[0][i] = value;
    }

    value = 1.0;
    for (int j = 1; j <= spatialOrder; j++) {
        value *= y;
        polyValues->data.F64[j][0] = value;
    }

    for (int j = 1; j <= spatialOrder; j++) {
        for (int i = 1; i <= spatialOrder - j; i++) {
            polyValues->data.F64[j][i] = polyValues->data.F64[j][0] * polyValues->data.F64[0][i];
        }
    }

    return polyValues;
}

// Generate the kernel to apply to the variance from the normal kernel
static psKernel *varianceKernel(psKernel *normalKernel // Normal kernel
                                )
{
    // Kernel range
    int xMin = normalKernel->xMin, xMax = normalKernel->xMax;
    int yMin = normalKernel->yMin, yMax = normalKernel->yMax;

    psKernel *kernel = psKernelAlloc(xMin, xMax, yMin, yMax); // Kernel to return

    // Take the square of the normal kernel
    double sumNormal = 0.0, sumVariance = 0.0; // Sum of the normal and variance kernels
    for (int v = yMin; v <= yMax; v++) {
        for (int u = xMin; u <= xMax; u++) {
            float value = normalKernel->kernel[v][u]; // Value of interest
            float value2 = PS_SQR(value); // Value squared
            sumNormal += value;
            sumVariance += value2;
            kernel->kernel[v][u] = value2;
        }
    }

    // Normalise so that the sum of the variance kernel is the square of the sum of the normal kernel
    // This is required to keep the relative scaling between the image and the weight map
    psBinaryOp(kernel->image, kernel->image, "*",
               psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32));

    return kernel;
}

// Generate an image of the solved kernel
static inline psKernel *solvedKernel(psKernel *kernel, // Kernel, to return
                                     const psVector *solution, // Solution to the least-squares problem
                                     const pmSubtractionKernels *kernels, // Kernel basis functions
                                     const psImage *polyValues // Spatial polynomial values
                                     )
{
    int numKernels = kernels->num;      // Number of kernel basis functions
    assert(kernels->u->n == numKernels);
    assert(kernels->v->n == numKernels);
    int spatialOrder = kernels->spatialOrder; // Maximum spatial polynomial order
    assert(solution->n == numKernels * polyTerms(spatialOrder) + polyTerms(kernels->bgOrder));

    // Ensure the subIndex for POIS kernels is what is expected
    assert((kernels->type != PM_SUBTRACTION_KERNEL_POIS &&
            kernels->type != PM_SUBTRACTION_KERNEL_SPAM &&
            kernels->type != PM_SUBTRACTION_KERNEL_FRIES &&
            kernels->type != PM_SUBTRACTION_KERNEL_GUNK &&
            kernels->type != PM_SUBTRACTION_KERNEL_RINGS) ||
           (kernels->u->data.S32[kernels->subIndex] == 0 && kernels->v->data.S32[kernels->subIndex] == 0));

    int size = kernels->size;           // Kernel half-size
    if (!kernel) {
        kernel = psKernelAlloc(-size, size, -size, size);
    }
    psImageInit(kernel->image, 0.0);

    for (int i = 0; i < numKernels; i++) {
        double value = 0.0;             // Value of this component, from the sum of multiple polynomial terms
        for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
            for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
                value += polyValues->data.F64[yOrder][xOrder] * solution->data.F64[index];
            }
        }

        switch (kernels->type) {
          case PM_SUBTRACTION_KERNEL_POIS: {
              int u = kernels->u->data.S32[i]; // Offset in x
              int v = kernels->v->data.S32[i]; // Offset in y
              kernel->kernel[v][u] += value;
              if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
                  // The (0,0) element is subtracted from most kernels to preserve photometric scaling
                  kernel->kernel[0][0] -= value;
              }
              break;
          }
          /* SPAM and FRIES use the same method */
          case PM_SUBTRACTION_KERNEL_SPAM:
          case PM_SUBTRACTION_KERNEL_FRIES: {
              int uStart = kernels->u->data.S32[i];
              int uStop = kernels->uStop->data.S32[i];
              int vStart = kernels->v->data.S32[i];
              int vStop = kernels->vStop->data.S32[i];

              // Normalising sum of kernel component to unity
              float norm = 1.0 / (float)((uStop - uStart + 1) * (vStop - vStart + 1));

              for (int v = vStart; v <= vStop; v++) {
                  for (int u = uStart; u <= uStop; u++) {
                      kernel->kernel[v][u] += norm * value;
                      if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
                          // The (0,0) element is subtracted from most kernels to preserve photometric scaling
                          kernel->kernel[0][0] -= value;
                      }
                  }
              }
              break;
          }
          case PM_SUBTRACTION_KERNEL_GUNK: {
              if (i < kernels->inner) {
                  // Using pre-calculated function
                  psKernel *preCalc = kernels->preCalc->data[i]; // Precalculated values
                  // Iterating over the kernel
                  for (int v = -size; v <= size; v++) {
                      for (int u = -size; u <= size; u++) {
                          kernel->kernel[v][u] += preCalc->kernel[v][u] * value;
                          // Photometric scaling is built into the preCalc kernel --- no subtraction!
                      }
                  }
              } else {
                  // Using delta function
                  bool subtract = (kernels->spatialOrder > 0 && i != kernels->subIndex); // Subtract (0,0)?
                  int u = kernels->u->data.S32[i]; // Offset in x
                  int v = kernels->v->data.S32[i]; // Offset in y
                  kernel->kernel[v][u] += value;
                  if (subtract) {
                      // The (0,0) element is subtracted from most kernels to preserve photometric scaling
                      kernel->kernel[0][0] -= value;
                  }
              }
              break;
          }
          case PM_SUBTRACTION_KERNEL_ISIS: {
              psKernel *preCalc = kernels->preCalc->data[i]; // Precalculated values
              // Iterating over the kernel
              for (int v = -size; v <= size; v++) {
                  for (int u = -size; u <= size; u++) {
                      kernel->kernel[v][u] += preCalc->kernel[v][u] * value;
                      // Photometric scaling is built into the preCalc kernel --- no subtraction!
                  }
              }
              break;
          }
          case PM_SUBTRACTION_KERNEL_RINGS: {
              if (i == kernels->subIndex) {
                  kernel->kernel[0][0] += value;
                  break;
              }
              psArray *preCalc = kernels->preCalc->data[i]; // Precalculated data
              psVector *uCoords = preCalc->data[0]; // u coordinates
              psVector *vCoords = preCalc->data[1]; // v coordinates
              psVector *poly = preCalc->data[2]; // Polynomial values
              int num = uCoords->n;     // Number of pixels

              for (int j = 0; j < num; j++) {
                  int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
                  kernel->kernel[v][u] += poly->data.F32[j] * value;
              }
              // Photometric scaling is built into the kernel --- no subtraction!
              break;
          }
          default:
            psAbort("Should never get here.");
        }
    }

    return kernel;
}

// Subtract the (0,0) element to preserve photometric scaling
static void convolveSub(psKernel *convolved, // Convolved image
                        const psKernel *image, // Image being convolved
                        int footprint  // Size of region of interest
                        )
{
    // Can't use psBinaryOp because the images are of different size
    for (int y = -footprint; y <= footprint; y++) {
        for (int x = -footprint; x <= footprint; x++) {
            convolved->kernel[y][x] -= image->kernel[y][x];
        }
    }
    return;
}

// Generate the convolution given some offset
static psKernel *convolveOffset(const psKernel *image, // Image to convolve (a kernel for convenience)
                                int u, int v, // Offset to apply
                                int footprint // Size of region of interest
                                )
{
    psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image
    int numBytes = (2 * footprint + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
    for (int y = -footprint; y <= footprint; y++) {
        // Convolution with a delta function is just the value specified by the offset
        memcpy(&convolved->kernel[y][-footprint], &image->kernel[y - v][-footprint - u], numBytes);
    }
    return convolved;
}

// Generate the convolved reference image
static psKernel *convolveRef(const pmSubtractionKernels *kernels, // Kernel basis functions
                             int index, // Kernel basis function index
                             const psKernel *image, // Image to convolve (a kernel for convenience)
                             int footprint // Size of region of interest
                             )
{
    switch (kernels->type) {
      case PM_SUBTRACTION_KERNEL_POIS: {
          int u = kernels->u->data.S32[index]; // Offset in x
          int v = kernels->v->data.S32[index]; // Offset in y
          psKernel *convolved = convolveOffset(image, u, v, footprint); // Convolved image
          if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
              convolveSub(convolved, image, footprint);
          }
          return convolved;
      }
        // Method for SPAM and FRIES is the same
      case PM_SUBTRACTION_KERNEL_SPAM:
      case PM_SUBTRACTION_KERNEL_FRIES: {
          psKernel *convolved = psKernelAlloc(-footprint, footprint,
                                              -footprint, footprint); // Convolved image
          int uStart = kernels->u->data.S32[index];
          int uStop = kernels->uStop->data.S32[index];
          int vStart = kernels->v->data.S32[index];
          int vStop = kernels->vStop->data.S32[index];
          float norm = 1.0 / (uStop - uStart + 1) * (vStop - vStart + 1); // Normalisation
          for (int y = -footprint; y <= footprint; y++) {
              for (int x = -footprint; x <= footprint; x++) {
                  double sum = 0.0;
                  for (int v = vStart; v <= vStop; v++) {
                      for (int u = uStart; u <= uStop; u++) {
                          sum += image->kernel[y - v][x - u];
                      }
                  }
                  convolved->kernel[y][x] = norm * sum;
              }
          }
          if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
              convolveSub(convolved, image, footprint);
          }
          return convolved;
      }
      case PM_SUBTRACTION_KERNEL_GUNK: {
          if (index < kernels->inner) {
              // Photometric scaling is already built in to the precalculated kernel
              return p_pmSubtractionConvolveStampPrecalc(image, kernels->preCalc->data[index]);
          }
          // Using delta function
          int u = kernels->u->data.S32[index]; // Offset in x
          int v = kernels->v->data.S32[index]; // Offset in y
          psKernel *convolved = convolveOffset(image, u, v, footprint); // Convolved image
          if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
              convolveSub(convolved, image, footprint);
          }
          return convolved;
      }
      case PM_SUBTRACTION_KERNEL_ISIS: {
          // Photometric scaling is already built in to the precalculated kernel
          return p_pmSubtractionConvolveStampPrecalc(image, kernels->preCalc->data[index]);
      }
      case PM_SUBTRACTION_KERNEL_RINGS: {
          psKernel *convolved = psKernelAlloc(-footprint, footprint,
                                              -footprint, footprint); // Convolved image
          if (index == kernels->subIndex) {
              // Simply copying over the image
              return convolveOffset(image, 0, 0, footprint);
          }

          psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data
          psVector *uCoords = preCalc->data[0]; // u coordinates
          psVector *vCoords = preCalc->data[1]; // v coordinates
          psVector *poly = preCalc->data[2]; // Polynomial values
          int num = uCoords->n;         // Number of pixels
          for (int y = -footprint; y <= footprint; y++) {
              for (int x = -footprint; x <= footprint; x++) {
                  double sum = 0.0;             // Accumulated sum from convolution
                  for (int j = 0; j < num; j++) {
                      int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
                      sum += image->kernel[y - v][x - u] * poly->data.F32[j];
                  }
                  convolved->kernel[y][x] = sum;
                  // Photometric scaling is built into the kernel --- no subtraction!
              }
          }
          return convolved;
      }
      default:
        psAbort("Should never get here.");
    }
    return NULL;
}

// Convolve an image using FFT
static void convolveFFT(psImage *target,// Place the result in here
                        const psImage *image, // Image to convolve
                        const psImage *mask, // Mask image
                        const psKernel *kernel, // Kernel by which to convolve
                        psRegion region,// Region of interest
                        float background, // Background to add
                        int size        // Size of (square) kernel
                        )
{
    psRegion border = psRegionSet(region.x0 - size, region.x1 + size,
                                  region.y0 - size, region.y1 + size); // Add a border

    // Casting away const so psImageSubset can add the child
    psImage *subImage = psImageSubset((psImage*)image, border); // Subimage to convolve
    psImage *subMask = mask ? psImageSubset((psImage*)mask, border) : NULL; // Subimage mask
    psImage *convolved = psImageConvolveFFT(subImage, subMask, PM_SUBTRACTION_MASK_REF,
                                            kernel, 0.0); // Convolution
    psFree(subImage);
    psFree(subMask);

    // Now, we have to chop off the borders, and stick it in where it belongs
    psImage *subConv = psImageSubset(convolved, psRegionSet(size, -size, size, -size)); // Cut off the edges
    psImage *subTarget = psImageSubset(target, region); // Target for this subregion
    if (background != 0.0) {
        psBinaryOp(subTarget, subConv, "+", psScalarAlloc(background, PS_TYPE_F32));
    } else {
        int numBytes = subTarget->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
        for (int y = 0; y < subTarget->numRows; y++) {
            memcpy(subTarget->data.F32[y], subConv->data.F32[y], numBytes);
        }
    }
    psFree(subConv);
    psFree(convolved);
    psFree(subTarget);
    return;
}

// Convolve an image directly
static void convolveDirect(psImage *target, // Put the result here
                           const psImage *image, // Image to convolve
                           const psKernel *kernel, // Kernel by which to convolve
                           psRegion region,// Region of interest
                           float background, // Background to add
                           int size        // Size of (square) kernel
                           )
{
    for (int y = region.y0; y < region.y1; y++) {
        for (int x = region.x0; x < region.x1; x++) {
            target->data.F32[y][x] = background;
            for (int v = -size; v <= size; v++) {
                for (int u = -size; u <= size; u++) {
                    target->data.F32[y][x] += kernel->kernel[v][u] *
                        image->data.F32[y - v][x - u];
                }
            }
        }
    }
    return;
}

// Mark a pixel as blank in the image, mask and weight
static inline void markBlank(psImage *image, // Image to mark as blank
                             psImage *mask, // Mask to mark as blank (or NULL)
                             psImage *weight, // Weight map to mark as blank (or NULL)
                             int x, int y, // Coordinates to mark blank
                             psMaskType blank // Blank mask value
    )
{
    image->data.F32[y][x] = NAN;
    if (mask) {
        mask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
    }
    if (weight) {
        weight->data.F32[y][x] = NAN;
    }
    return;
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Semi-public functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

// Generate the convolution given a precalculated kernel
psKernel *p_pmSubtractionConvolveStampPrecalc(const psKernel *image, const psKernel *kernel)
{
    PS_ASSERT_KERNEL_NON_NULL(image, NULL);
    PS_ASSERT_KERNEL_NON_NULL(kernel, NULL);

    psImage *conv = psImageConvolveFFT(image->image, NULL, 0, kernel, 0.0); // Convolved image
    int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image
    psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version
    psFree(conv);
    return convolved;
}


//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Public functions
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

psImage *pmSubtractionMask(const psImage *refMask, const psImage *inMask, psMaskType maskVal,
                           int size, int footprint)
{
    PS_ASSERT_IMAGE_NON_NULL(refMask, NULL);
    PS_ASSERT_IMAGE_TYPE(refMask, PS_TYPE_MASK, NULL);
    if (inMask) {
        PS_ASSERT_IMAGE_NON_NULL(inMask, NULL);
        PS_ASSERT_IMAGE_TYPE(inMask, PS_TYPE_MASK, NULL);
        PS_ASSERT_IMAGES_SIZE_EQUAL(inMask, refMask, NULL);
    }
    PS_ASSERT_INT_NONNEGATIVE(size, NULL);
    PS_ASSERT_INT_NONNEGATIVE(footprint, NULL);

    // Size of the images
    int numCols = refMask->numCols;
    int numRows = refMask->numRows;

    // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask
    psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask
    psImageInit(mask, 0);

    // Dereference for convenience
    psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA;
    psMaskType **inData = NULL;
    if (inMask) {
        inData = inMask->data.PS_TYPE_MASK_DATA;
    }
    psMaskType **refData = refMask->data.PS_TYPE_MASK_DATA;

    // Block out a border around the edge of the image

    // Bottom stripe
    for (int y = 0; y < PS_MIN(size + footprint, numRows); y++) {
        for (int x = 0; x < numCols; x++) {
            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
        }
    }
    // Either side
    for (int y = PS_MIN(size + footprint, numRows); y < numRows - size - footprint; y++) {
        for (int x = 0; x < PS_MIN(size + footprint, numCols); x++) {
            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
        }
        for (int x = PS_MAX(numCols - size - footprint, 0); x < numCols; x++) {
            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
        }
    }
    // Top stripe
    for (int y = PS_MAX(numRows - size - footprint, 0); y < numRows; y++) {
        for (int x = 0; x < numCols; x++) {
            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
        }
    }

    // Mask the bad pixels, and around them
    for (int y = 0; y < numRows; y++) {
        for (int x = 0; x < numCols; x++) {
            if (inData && inData[y][x] & maskVal) {
                maskData[y][x] |= PM_SUBTRACTION_MASK_INPUT;
                // Block out the entire stamp footprint around this pixel
                for (int v = PS_MAX(y - footprint, 0); v <= PS_MIN(y + footprint, numRows - 1); v++) {
                    for (int u = PS_MAX(x - footprint, 0); u <= PS_MIN(x + footprint, numCols - 1); u++) {
                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
                    }
                }
            }
            if (refData[y][x] & maskVal) {
                maskData[y][x] |= PM_SUBTRACTION_MASK_REF;

                // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with
                // this bad pixel (within 'size').  Then we want to block out with the FOOTPRINT mask
                // everything within a footprint's distance of those (within 'footprint').

                // Bottom stripe
                for (int v = PS_MAX(y - footprint - size, 0); v < y - size; v++) {
                    for (int u = PS_MAX(x - footprint - size, 0);
                         u <= PS_MIN(x + footprint + size, numCols - 1); u++) {
                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
                    }
                }
                // Middle
                for (int v = PS_MAX(y - size, 0); v <= PS_MIN(y + size, numRows - 1); v++) {
                    // Left side
                    for (int u = PS_MAX(x - footprint - size, 0); u < x - size; u++) {
                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
                    }
                    // Centre
                    for (int u = PS_MAX(x - size, 0); u <= PS_MIN(x + size, numCols - 1); u++) {
                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_CONVOLVE;
                    }
                    // Right side
                    for (int u = x + size + 1; u <= PS_MIN(x + footprint + size, numCols - 1); u++) {
                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
                    }
                }
                // Top stripe
                for (int v = y + size + 1; v <= PS_MIN(y + footprint + size, numRows - 1); v++) {
                    for (int u = PS_MAX(x - footprint - size, 0);
                         u <= PS_MIN(x + footprint + size, numCols - 1); u++) {
                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
                    }
                }
            }
        }
    }

    return mask;
}

bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, const pmSubtractionKernels *kernels, int footprint)
{
    PS_ASSERT_PTR_NON_NULL(stamp, false);
    PS_ASSERT_KERNEL_NON_NULL(stamp->reference, false);
    PS_ASSERT_PTR_NON_NULL(kernels, false);
    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false);
    PS_ASSERT_INT_NONNEGATIVE(footprint, false);

    if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Stamp not marked for calculation.");
        return false;
    }

    if (stamp->convolutions) {
        // Already done
        return true;
    }

    stamp->convolutions = psArrayAlloc(kernels->num);
    psKernel *reference = stamp->reference; // Reference postage stamp

    for (int i = 0; i < kernels->num; i++) {
        stamp->convolutions->data[i] = convolveRef(kernels, i, reference, footprint);
    }

    return true;
}


bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
                                    int footprint)
{
    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    PS_ASSERT_PTR_NON_NULL(kernels, false);
    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false);
    PS_ASSERT_INT_NONNEGATIVE(footprint, false);

    int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation
    int numKernels = kernels->num;      // Number of kernel basis functions
    int numSpatial = polyTerms(spatialOrder); // Number of spatial variations
    int numBackground = polyTerms(kernels->bgOrder); // Number of background terms

    // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the
    // number of coefficients for the spatial polynomial, and a constant background offset.
    int numParams = numKernels * numSpatial + numBackground;
    int bgIndex = numParams - numBackground; // Index in matrix for the background

    // We iterate over each stamp, allocate the matrix and vectors if
    // necessary, and then calculate those matrix/vectors.
    for (int i = 0; i < stamps->num; i++) {
        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
        if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
            continue;
        }

        psImage *stampMatrix = stamp->matrix; // Least-squares matrix for this stamp
        psVector *stampVector = stamp->vector; // Least-squares vector for this stamp

        if (!stampMatrix) {
            stamp->matrix = stampMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
        }
        assert(stampMatrix->type.type == PS_TYPE_F64);
        assert(stampMatrix->numCols == numParams && stampMatrix->numRows == numParams);
        psImageInit(stampMatrix, 0.0);

        if (!stampVector) {
            stamp->vector = stampVector = psVectorAlloc(numParams, PS_TYPE_F64);
        }
        assert(stampVector->type.type == PS_TYPE_F64);
        assert(stampVector->n == numParams);
        psVectorInit(stampVector, 0.0);

        psKernel *input = stamp->input; // Input postage stamp
        psKernel *weight = stamp->weight; // Weight map postage stamp

        // Generate convolutions of the reference
        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i);
            return NULL;
        }

        psArray *convolutions = stamp->convolutions; // Convolutions of the reference for each kernel
        psImage *polyValues = spatialPolyValues(spatialOrder, stamp->xNorm, stamp->yNorm); // Polynomial terms

        // Generate least-squares vector and matrix
        bool bad = false;           // Are there bad values?
        for (int j = 0; j < numKernels; j++) {
            psKernel *jConv = convolutions->data[j]; // Convolution for j-th element

            // Generate matrix
            for (int k = j; k < numKernels; k++) {
                psKernel *kConv = convolutions->data[k]; // Convolution for k-th element
                double sumCC = 0.0; // Sum of the convolution products
                for (int y = - footprint; y <= footprint; y++) {
                    for (int x = - footprint; x <= footprint; x++) {
                        sumCC += jConv->kernel[y][x] * kConv->kernel[y][x] / weight->kernel[y][x];
                    }
                }
                if (!isfinite(sumCC)) {
                    bad = true;
                }
                // Generate the pseudo-convolutions from the spatial polynomial terms
                for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) {
                    for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder;
                         jxOrder++, jIndex += numKernels) {
                        for (int kyOrder = 0, kIndex = k; kyOrder <= spatialOrder; kyOrder++) {
                            for (int kxOrder = 0; kxOrder <= spatialOrder - kyOrder;
                                 kxOrder++, kIndex += numKernels) {
                                double convPoly = sumCC * polyValues->data.F64[jyOrder][jxOrder] *
                                    polyValues->data.F64[kyOrder][kxOrder]; // Temporary value
                                stampMatrix->data.F64[jIndex][kIndex] = convPoly;
                                stampMatrix->data.F64[kIndex][jIndex] = convPoly;
                            }
                        }
                    }
                }
            }

            // Vector and background term for matrix
            double sumC = 0.0;      // Sum of the convolution
            double sumIC = 0.0;     // Sum of the convolution/input product
            for (int y = - footprint; y <= footprint; y++) {
                for (int x = - footprint; x <= footprint; x++) {
                    double convProduct = jConv->kernel[y][x] / weight->kernel[y][x]; // Convolution / noise^2
                    sumC += convProduct;
                    sumIC += input->kernel[y][x] * convProduct;
                }
            }
            if (!isfinite(sumC) || !isfinite(sumIC)) {
                bad = true;
            }
            for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) {
                for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder;
                     jxOrder++, jIndex += numKernels) {
                    stampVector->data.F64[jIndex] = sumIC * polyValues->data.F64[jyOrder][jxOrder];

                    double convPoly = sumC * polyValues->data.F64[jyOrder][jxOrder]; // Value
                    stampMatrix->data.F64[jIndex][bgIndex] = convPoly;
                    stampMatrix->data.F64[bgIndex][jIndex] = convPoly;
                }
            }

        }
        psFree(polyValues);

        // Background only terms
        double sum1 = 0.0;              // Sum of the weighting
        double sumI = 0.0;              // Sum of the input
        for (int y = - footprint; y <= footprint; y++) {
            for (int x = - footprint; x <= footprint; x++) {
                double invNoise2 = 1.0 / weight->kernel[y][x]; // Inverse noise, squared
                sum1 += invNoise2;
                sumI += input->kernel[y][x] * invNoise2;
            }
        }
        if (!isfinite(sum1) || !isfinite(sumI)) {
            bad = true;
        }
        stampMatrix->data.F64[bgIndex][bgIndex] = sum1;
        stampVector->data.F64[bgIndex] = sumI;

        if (bad) {
            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
            psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d) because of bad equation\n",
                    i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
        } else {
            stamp->status = PM_SUBTRACTION_STAMP_USED;
        }

        if (psTraceGetLevel("psModules.imcombine.equation") >= 10) {
            psString matrixName = NULL;
            psStringAppend(&matrixName, "matrix%d.fits", i);
            psFits *matrixFile = psFitsOpen(matrixName, "w");
            psFree(matrixName);
            psFitsWriteImage(matrixFile, NULL, stampMatrix, 0, NULL);
            psFitsClose(matrixFile);
        }

    }

    return true;
}


psVector *pmSubtractionSolveEquation(psVector *solution, const pmSubtractionStampList *stamps)
{
    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL);

    // Check inputs, while summing the stamp matrices and vectors
    long numParams = -1;                // Number of parameters
    for (int i = 0; i < stamps->num; i++) {
        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
        PS_ASSERT_PTR_NON_NULL(stamp, NULL);
        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
            continue;
        }
        if (numParams == -1) {
            numParams = stamp->vector->n;
        } else {
            PS_ASSERT_VECTOR_SIZE(stamp->vector, numParams, NULL);
        }
        PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, NULL);
        PS_ASSERT_IMAGE_TYPE(stamp->matrix, PS_TYPE_F64, NULL);
        PS_ASSERT_IMAGE_SIZE(stamp->matrix, numParams, numParams, NULL);
    }
    if (numParams == -1) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No suitable stamps found.");
        return NULL;
    }

    if (solution) {
        PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL);
        PS_ASSERT_VECTOR_SIZE(solution, numParams, NULL);
    }

    // Accumulate the least-squares matricies and vectors
    psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); // Combined matrix for all stamps
    psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); // Combined vector for all stamps
    psVectorInit(sumVector, 0.0);
    psImageInit(sumMatrix, 0.0);
    for (int i = 0; i < stamps->num; i++) {
        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
        if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
            (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
            (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
        }
    }

    psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
    psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix);
    psFree(sumMatrix);
    if (!luMatrix) {
        psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
        psFree(sumVector);
        psFree(luMatrix);
        psFree(permutation);
        return NULL;
    }

    solution = psMatrixLUSolve(solution, luMatrix, sumVector, permutation);
    psFree(sumVector);
    psFree(luMatrix);
    psFree(permutation);
    if (!solution) {
        psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
        return NULL;
    }

    if (psTraceGetLevel("psModules.imcombine") >= 7) {
        for (int i = 0; i < solution->n; i++) {
            printf("       Solution %d: %f\n", i, solution->data.F64[i]);
        }
    }

    return solution;
}


int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, psImage *subMask, const psVector *solution,
                              int footprint, float sigmaRej, const pmSubtractionKernels *kernels)
{
    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, -1);
    PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1);
    PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1);
    PS_ASSERT_VECTOR_NON_NULL(solution, -1);
    PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, -1);
    PS_ASSERT_PTR_NON_NULL(kernels, -1);
    PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) +
                          polyTerms(kernels->bgOrder), -1);
    PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, -1);
    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, -1);

    // Measure deviations
    psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps
    double totalSquareDev = 0.0;        // Total square deviation from zero
    int numStamps = 0;                  // Number of used stamps
    int numKernels = kernels->num;      // Number of kernels
    int spatialOrder = kernels->spatialOrder; // Order of kernel spatial variations
    double devNorm = 1.0 / PS_SQR(2 * footprint + 1); // Normalisation for deviations
    {
        float background = solution->data.F64[solution->n-1]; // The difference in background

        for (int i = 0; i < stamps->num; i++) {
            pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
            if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
                continue;
            }

            psImage *polyValues = spatialPolyValues(kernels->spatialOrder, stamp->xNorm,
                                                    stamp->yNorm); // Polynomial terms

            psKernel *input = stamp->input; // Input postage stamp
            psKernel *weight = stamp->weight; // Weight postage stamp
            psArray *convolutions = stamp->convolutions; // Convolutions of reference image for each kernel
#ifdef TESTING
            psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint);
#endif
            float deviation = 0.0;      // Deviation for this stamp
            for (int y = - footprint; y <= footprint; y++) {
                for (int x = - footprint; x <= footprint; x++) {
                    float conv = background; // The value of the convolution
                    for (int j = 0; j < numKernels; j++) {
                        psKernel *convolution = convolutions->data[j]; // Convolution of reference
                        double polynomial = 0.0; // Value of the polynomial
                        for (int yOrder = 0, index = j; yOrder <= spatialOrder; yOrder++) {
                            for (int xOrder = 0; xOrder <= spatialOrder - yOrder;
                                 xOrder++, index += numKernels) {
                                polynomial += solution->data.F64[index] *
                                    polyValues->data.F64[yOrder][xOrder];
                            }
                        }
                        conv += convolution->kernel[y][x] * polynomial;
                    }
                    float diff = input->kernel[y][x] - conv;
                    deviation += PS_SQR(diff) / weight->kernel[y][x];
#ifdef TESTING
                    residual->kernel[y][x] = diff;
#endif
                }
            }
            psFree(polyValues);

            deviations->data.F32[i] = devNorm * deviation;
            psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
                    i, (int)stamp->x, (int)stamp->y, deviations->data.F32[i]);
            totalSquareDev += PS_SQR(deviations->data.F32[i]);
            numStamps++;

#ifdef TESTING
            {
                psString filename = NULL;
                psStringAppend(&filename, "resid_%03d.fits", i);
                psFits *fits = psFitsOpen(filename, "w");
                psFree(filename);
                psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
                psFitsClose(fits);
            }
            psFree(residual);
#endif

        }
    }

    float limit = sigmaRej * sqrt(totalSquareDev / numStamps); // Limit on maximum deviation
    psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit);

    int numRejected = 0;
    for (int i = 0; i < stamps->num; i++) {
        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
        if (stamp->status == PM_SUBTRACTION_STAMP_USED && deviations->data.F32[i] > limit) {
            // Mask out the stamp in the image so you it's not found again
            psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i,
                    (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
            numRejected++;
            for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
                for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) {
                    subMask->data.PS_TYPE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ;
                }
            }

            // Set stamp for replacement
            stamp->x = 0;
            stamp->y = 0;
            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
            // Recalculate convolutions
            psFree(stamp->convolutions);
            stamp->convolutions = NULL;
        }
    }

    psFree(deviations);

    return numRejected;
}

psImage *pmSubtractionKernelImage(const psVector *solution, const pmSubtractionKernels *kernels,
                                  float x, float y
                                 )
{
    PS_ASSERT_VECTOR_NON_NULL(solution, NULL);
    PS_ASSERT_PTR_NON_NULL(kernels, NULL);
    PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) +
                          polyTerms(kernels->bgOrder), NULL);
    PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL);
    PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL);
    PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, false);
    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL);

    // Precalulate polynomial values
    psImage *polyValues = spatialPolyValues(kernels->spatialOrder, x, y);

    // The appropriate kernel
    psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues);

    psFree(polyValues);

    psImage *image = psMemIncrRefCounter(kernel->image); // Image of the kernel
    psFree(kernel);

    return image;
}

psArray *pmSubtractionKernelSolutions(const psVector *solution, const pmSubtractionKernels *kernels,
                                      float x, float y
                                      )
{
    PS_ASSERT_VECTOR_NON_NULL(solution, NULL);
    PS_ASSERT_PTR_NON_NULL(kernels, NULL);
    PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) +
                          polyTerms(kernels->bgOrder), NULL);
    PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL);
    PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL);
    PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, false);
    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL);

    psArray *images = psArrayAlloc(solution->n - 1); // Images of each kernel to return
    psVector *fakeSolution = psVectorAlloc(solution->n, PS_TYPE_F64); // Fake solution vector
    psVectorInit(fakeSolution, 0.0);

    for (int i = 0; i < solution->n - 1; i++) {
        fakeSolution->data.F64[i] = solution->data.F64[i];
        images->data[i] = pmSubtractionKernelImage(fakeSolution, kernels, x, y);
        fakeSolution->data.F64[i] = 0.0;
    }

    psFree(fakeSolution);

    return images;
}

bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask,
                           const psImage *inImage, const psImage *inWeight, const psImage *subMask,
                           psMaskType blank, const psRegion *region, const psVector *solution,
                           const pmSubtractionKernels *kernels, bool useFFT)
{
    PS_ASSERT_IMAGE_NON_NULL(inImage, false);
    PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, false);
    if (subMask) {
        PS_ASSERT_IMAGE_NON_NULL(subMask, false);
        PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false);
        PS_ASSERT_IMAGES_SIZE_EQUAL(subMask, inImage, false);
    }
    if (inWeight) {
        PS_ASSERT_IMAGE_NON_NULL(inWeight, false);
        PS_ASSERT_IMAGE_TYPE(inWeight, PS_TYPE_F32, false);
        PS_ASSERT_IMAGES_SIZE_EQUAL(inWeight, inImage, false);
    }
    PS_ASSERT_PTR_NON_NULL(outImage, false);
    if (*outImage) {
        PS_ASSERT_IMAGE_NON_NULL(*outImage, false);
        PS_ASSERT_IMAGES_SIZE_EQUAL(*outImage, inImage, false);
        PS_ASSERT_IMAGE_TYPE(*outImage, PS_TYPE_F32, false);
    }
    if (outMask && *outMask) {
        PS_ASSERT_IMAGE_NON_NULL(*outMask, false);
        PS_ASSERT_IMAGES_SIZE_EQUAL(*outMask, inImage, false);
        PS_ASSERT_IMAGE_TYPE(*outMask, PS_TYPE_MASK, false);
        PS_ASSERT_IMAGE_NON_NULL(subMask, false);
    }
    if (outWeight && *outWeight) {
        PS_ASSERT_IMAGE_NON_NULL(*outWeight, false);
        PS_ASSERT_IMAGES_SIZE_EQUAL(*outWeight, inImage, false);
        PS_ASSERT_IMAGE_TYPE(*outWeight, PS_TYPE_F32, false);
        PS_ASSERT_IMAGE_NON_NULL(inWeight, false);
    }
    PS_ASSERT_VECTOR_NON_NULL(solution, false);
    PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, false);
    PS_ASSERT_PTR_NON_NULL(kernels, false);
    PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) +
                          polyTerms(kernels->bgOrder), -1);
    PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, false);
    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false);
    if (region) {
        if (psRegionIsNaN(*region)) {
            psString string = psRegionToString(*region);
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string);
            psFree(string);
            return false;
        }
        if (region->x0 < 0 || region->x1 > inImage->numCols ||
            region->y0 < 0 || region->y1 > inImage->numRows) {
            psString string = psRegionToString(*region);
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)",
                    string, inImage->numCols, inImage->numRows);
            psFree(string);
            return false;
        }
    }

    int numBackground = polyTerms(kernels->bgOrder); // Number of background terms
    float background = solution->data.F64[solution->n - numBackground]; // The difference in background
    int numCols = inImage->numCols, numRows = inImage->numRows; // Image dimensions

    psImage *convImage = *outImage;    // Convolved image
    if (!convImage) {
        *outImage = psImageAlloc(numCols, numRows, PS_TYPE_F32);
        convImage = *outImage;
    }
    psImage *convMask = NULL;           // Convolved mask image
    if (outMask && subMask) {
        if (!*outMask) {
            *outMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
        }
        convMask = *outMask;
        psImageInit(convMask, 0);
    }
    psImage *convWeight = NULL;         // Convolved weight (variance) image
    if (outWeight && inWeight) {
        if (!*outWeight) {
            *outWeight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
        }
        convWeight = *outWeight;
        psImageInit(convWeight, 0.0);
    }

    int size = kernels->size;           // Half-size of kernel
    int fullSize = 2 * size + 1;        // Full size of kernel
    psImage *polyValues = NULL;         // Pre-calculated polynomial values

    psKernel *kernelImage = NULL;       // Kernel for the image
    psKernel *kernelWeight = NULL;      // Kernel for the weight map

    // Get region for convolution: [xMin:xMax,yMin:yMax]
    int xMin = size, xMax = numCols - size;
    int yMin = size, yMax = numRows - size;
    if (region) {
        xMin = PS_MAX(region->x0, xMin);
        xMax = PS_MIN(region->x1, xMax);
        yMin = PS_MAX(region->y0, yMin);
        yMax = PS_MIN(region->y1, yMax);
    }

    for (int j = yMin; j < yMax; j += fullSize) {
        int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest
        for (int i = xMin; i < xMax; i += fullSize) {
            int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest

            // Only generate polynomial values every kernel footprint, since we have already assumed
            // (with the stamps) that it does not vary rapidly on this scale.
            psFree(polyValues);
            polyValues = spatialPolyValues(kernels->spatialOrder,
                                           2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols,
                                           2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows);

            kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues);
            if (inWeight) {
                kernelWeight = varianceKernel(kernelImage);
            }


            if (useFFT) {
                // Use Fast Fourier Transform to do the convolution
                // This provides a big speed-up for large kernels
                convolveFFT(convImage, inImage, subMask, kernelImage, psRegionSet(i, xSubMax, j, ySubMax),
                            background, size);
                if (inWeight) {
                    convolveFFT(convWeight, inWeight, subMask, kernelWeight,
                                psRegionSet(i, xSubMax, j, ySubMax), 0.0, size);
                }
            } else {
                convolveDirect(convImage, inImage, kernelImage, psRegionSet(i, xSubMax, j, ySubMax),
                               background, size);
                if (inWeight) {
                    convolveDirect(convWeight, inWeight, kernelWeight, psRegionSet(i, xSubMax, j, ySubMax),
                                   0.0, size);
                }
            }

            // Propagate the mask
            for (int y = j; y < ySubMax; y++) {
                for (int x = i; x < xSubMax; x++) {
                    if (subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] &
                                    (PM_SUBTRACTION_MASK_INPUT | PM_SUBTRACTION_MASK_CONVOLVE))) {
                        convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
                        convImage->data.F32[y][x] = NAN;
                        if (inWeight) {
                            convWeight->data.F32[y][x] = NAN;
                        }
                    }
                }
            }

        }
    }
    psFree(kernelImage);
    psFree(kernelWeight);
    psFree(polyValues);

    return true;
}

bool pmSubtractionBorder(psImage *image, psImage *weight, psImage *mask,
                         int size, psMaskType blank)
{
    PS_ASSERT_IMAGE_NON_NULL(image, false);
    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
    if (mask) {
        PS_ASSERT_IMAGE_NON_NULL(mask, false);
        PS_ASSERT_IMAGES_SIZE_EQUAL(mask, image, false);
        PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, false);
    }
    if (weight) {
        PS_ASSERT_IMAGE_NON_NULL(weight, false);
        PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image, false);
        PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
    }

    int numCols = image->numCols, numRows = image->numRows; // Image dimensions

    for (int y = size; y < numRows - size; y++) {
        for (int x = 0; x < size; x++) {
            markBlank(image, mask, weight, x, y, blank);
        }
        for (int x = numCols - size; x < numCols; x++) {
            markBlank(image, mask, weight, x, y, blank);
        }
    }
    for (int y = 0; y < size; y++) {
        for (int x = 0; x < numCols; x++) {
            markBlank(image, mask, weight, x, y, blank);
        }
    }
    for (int y = numRows - size; y < numRows; y++) {
        for (int x = 0; x < numCols; x++) {
            markBlank(image, mask, weight, x, y, blank);
        }
    }

    return true;
}


