// A list of pixel coordinates
typedef struct {
    psVector *x;			// x coordinate
    psVector *y;			// y coordinate
} psPixels;

// Convert from a list of pixels to a U8 mask image
psImage *psPixelsToMask(psImage *out,// Output mask image, or NULL
			const psPixels *pixels, // Pixels to convert
			const psRegion *region,	// Range of values; also specifies size of output image
			unsigned int maskVal // Value to give mask
			);

// Convert from a U8 mask image to a list of pixels
psPixels *psMaskToPixels(psPixels *out,	// Output pixels, or NULL
			 const psImage *mask, // Mask image to convert
			 unsigned int maskVal // Values to convert
			 );

// Concatenate a list of pixels onto another list
psPixels *psPixelsConcatenate(psPixels *out, // Output pixels, or NULL
			      const psPixels *pixels // Pixels to add to list
			      );

// Return combined image
psImage *pmCombineImages(psImage *combined, // Combined image
			 psArray **rejected, // Array of rejection masks
			 const psArray *images, // Array of input images
			 const psArray *errors, // Array of input error images
			 const psArray *masks,// Array of input masks
			 unsigned int maskVal, // Mask value
			 const psArray *pixels,	// Pixels to combine
			 int numIter,	// Number of rejection iterations
			 float sigmaClip, // Number of standard deviations at which to reject
			 const psStats *stats // Statistics to use in the combination
			 );

// Return the derivative of the transformation at the specified coordinates
psPlane *psPlaneTransformDeriv(psPlane *out, // Output value, or NULL
			       const psPlaneTransform *map, // Transformation to apply
			       const psPlane *in // Coordinates at which to evaluate the derivative
			       );

// Transform a list of pixels from one coordinate frame to another.
// Include overlaps (need to account for derivatives).
psPixels *psPixelsTransform(psPixels *out, // Output pixels, or NULL
			    const psPixels *in, // Input pixels
			    const psPlaneTransform *inToOut // Transformation map
			    );

// Transform an image --- different from current version, nothing major
psImage *psImageTransform(psImage *output, // Output image, or NULL
			  psPixels **blankPixels, // List of blank pixels
			  const psImage *input,	// Input image
			  const psImage *inputMask, // Input image mask, or NULL
			  int inputMaskVal, // Mask value
			  double exposedValue, // Value for "exposed" pixels in the transformation
                          const psPlaneTransform *outToIn, // Transformation
			  const psRegion *region, // Region of output to consider (defines image size)
			  const psPixels *pixels, // Pixels in output image to transform, or NULL for all
			  psImageInterpolateMode mode // Mode of interpolation
			  );

// Return an array of psPixels that are rejected
psArray *pmRejectPixels(const psArray *images, // Array of input images
			const psArray *pixels, // These are the pixels which were rejected in the combination
			const psPlaneTransform *inToOut, // Transformation from input to output system
			const psPlaneTransform *outToIn, // Transformation from output to input system
			float rejThreshold, // Rejection threshold
			float gradLimit // Gradient limit
			);


// The image combination looks something like this:
psImage *imageCombination(const psArray *inputs, // Input detector images
			  psArray *inputMask,	// Input mask images
			  const psArray *inputsErr, // Input error images
			  const psArray *skyToDetector, // Maps from sky coordinates to detector coordinates
			  const psRegion *combineRegion, // Sky coordinate pixels to combine
			  int numIter,	// Number of iterations
			  float rejThreshold, // Threshold for rejection
			  float gradLimit // Limit for gradient
	      )
{
    psArray *transformed = psArrayAlloc(nImages); // Array of transformed images
    psArray *transformedErr = psArrayAlloc(nImages); // Array of transformed error images
    psArray *transformedMask = psArrayAlloc(nImages); // Array of masks for transformed images

    for (int i = 0; i < nImages; i++) {
	psPixels *blanks = NULL;	// List of blank pixels
	transformed->data[i] = psImageTransform(NULL, &blanks, inputs->data[i],
						inputMask->data[i], inputMaskVal, NAN, skyToDetector,
						combineRegion, NULL, PS_INTERPOLATE_BILINEAR);
	transformedErr->data[i] = psImageTransform(NULL, NULL, inputsErr->data[i], inputMask->data[i],
						   inputMaskVal, NAN, skyToDetector, combineRegion, NULL,
						   PS_INTERPOLATE_BILINEAR_VARIANCE);
	psImage *skyImage = transformed->data[i]; // Dereference the transformed image
	psRegion *blankRegion = psRegionAlloc(0, 0, skyImage->numCols, skyImage->numRows); // Size of
											   // transformed
											   // image
	transformedMask->data[i] = psPixelsToMask(NULL, blanks, blankRegion, PS_MASK_BLANK);
	psFree(blankRegion);
	psFree(blanks);
    }

    psArray *rejected = NULL;		// Array of rejected pixel lists
    psStats *combineStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistic to use in doing the combination
    psImage *combined = pmCombineImages(NULL, &rejected, transformed, transformedErr, transformedMask, 0,
					NULL, numIter, sigmaClip, combineStats); // Combined image
    psArray *bad = pmRejectPixels(inputs, rejected, NULL, skyToDetector, rejThreshold, gradLimit); // Bad pix
    psPixels *combinePixels = NULL;	// Pixels to combine
    for (int i = 0; i < nImages; i++) {
	psPixels *badSource = psPixelsTransform(NULL, bad->data[i], skyToDetector); // Bad pixels on the input
	psImage *badMask = psPixelsToMask(NULL, badSource, PS_MASK_COSMICRAY); // Mask image for the input
	(void)psBinaryOp(inputMask->data[i], inputMask->data[i], "|", badMask);	// Put CRs into original mask
	psFree(badSource);
	psFree(badMask);

	combinePixels = psPixelsConcatenate(redo, bad->data[i]);

	// Update transformed image
	psPixels *blanks = NULL;	// List of blank pixels
	transformed->data[i] = psImageTransform(transformed->data[i], &blanks, inputs->data[i],
						inputMask->data[i], inputMaskVal | PS_MASK_COSMICRAY, NAN,
						skyToDetector, combineRegion, bad->data[i],
						PS_INTERPOLATE_BILINEAR);
	transformedErr->data[i] = psImageTransform(transformedErr->data[i], NULL, inputsErr->data[i],
						   inputMask->data[i], inputMaskVal | PS_MASK_COSMICRAY,
						   NAN, skyToDetector, combineRegion, bad->data[i],
						   PS_INTERPOLATE_BILINEAR_VARIANCE);
	psImage *skyImage = transformed->data[i]; // Dereference the transformed image
	psRegion *blankRegion = psRegionAlloc(0, 0, skyImage->numCols, skyImage->numRows); // Size of
											   // transformed
											   // image
	transformedMask->data[i] = psPixelsToMask(transformedMask->data[i], blanks, blankRegion,
						  PS_MASK_BLANK);
	psFree(blankRegion);
	psFree(blanks);
    }
    psFree(bad);

    // Combine with no rejection
    combined = pmCombineImages(combined, NULL, transformed, transformedErr, transformedMask,
			       PS_MASK_BLANK, combinePixels, 0, 0.0, combineStats);
    psFree(combineStats);
    psFree(combinePixels);
    psFree(transformed);
    psFree(transformedErr);
    psFree(transformedMask);

    return combined;
}


//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Image subtraction APIs                                                                                   //
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

typedef enum {
    PM_STAMP_USED,			// Use this stamp
    PM_STAMP_RESET,			// This stamp has been rejected --- reset
    PM_STAMP_RECALC,			// Having been reset, this stamp needs to be recalculated
    PM_STAMP_NONE			// No stamp in this region
} pmStampStatus;

typedef struct {
    int x, y;				// Position
    psImage *matrix;			// Associated matrix
    psVector *vector;			// Associated vector
    psStampStatus status;		// Status of stamp
} pmStamp;

// There's no reason why we can't do *both* ISIS and POIS --- it only takes a couple of changes.
// Let's leave in the flexibility.

typedef enum {
    PM_SUBTRACTION_KERNEL_POIS,		// POIS kernel --- delta functions
    PM_SUBTRACTION_KERNEL_ISIS		// ISIS kernel --- gaussians modified by polynomials
} pmSubtractionKernelType;

typedef struct {
    pmSubtractionKernelType type;	// Type of kernel --- allowing the use of multiple kernels
    int u, v;				// Offset (for POIS) or polynomial order (for ISIS)
    int xOrder, yOrder;			// Spatial polynomial order (for all)
    float sigma;			// Width of Gaussian (for ISIS)
    psImage *preCalc;			// Pre-calculated kernel (to accelerate ISIS; don't use for POIS)
} pmSubtractionKernel;

// Create the kernel basis functions
psArray *pmSubtractionKernelsGeneratePOIS(int size, // Half-size of kernel (full size is 2*size+1)
					  int spatialOrder // Maximum spatial order
					  );
psArray *pmSubtractionKernelsGenerateISIS(const psVector *sigmas, // Gaussian widths
					  const psVector *orders, // Polynomial orders for the gaussians
					  int size, // Half-size of kernel (full size is 2*size+1)
					  int spatialOrder	// Maximum spatial order
					  );

psImage *psImageGrowMask(psImage *out,	// Output mask, or NULL
			 const psImage *in, // Input mask
			 unsigned int maskVal, // Value to identify
			 unsigned int growRadius, // Radius by which to grow
			 unsigned int growValue	// Value to OR with grown pixels
			 );

// Return an array of stamps
psArary *pmSubtractionFindStamps(psArray *stamps, // Output stamps, or NULL
				 const psImage *image, // Image for which to find stamps
				 const psImage *mask, // Mask
				 unsigned int maskVal, // Value for mask
				 float threshold, // Threshold for stamps in the image
				 int xNum, int yNum, // Number of stamps in x and y
				 int border // Border around image to ignore (should be size of kernel)
				 );

// Calculate the matrix equations for each kernel basis function
int pmSubtractionCalculateEquation(psArray *stamps, // The stamps for which to calculate the equation
				   const psImage *reference, // Reference image
				   const psImage *input, // Input image
				   const psArray *kernels, // The kernel basis functions
				   int footprint // Half-size of region over which to calculate equation
				   );

// Solve the matrix equation (sum the stamps, solve)
psVector *pmSubtractionSolveEquation(psVector *solution,	// Solution vector, or NULL
				     const psArray *stamps // Array of stamps
				     );

// Convolve an image by the kernel
psImage *pmSubtractionConvolveImage(psImage *out,// Output image, or NULL
				    const psImage *in, // Input image, to be convolved
				    const psImage *mask, // Mask image
				    const psVector *solution, // The solution vector (the coefficients)
				    const psArray *kernels, // The kernel parameters
				    int footprint // Size of region over which spatial variation is small (to optimise)
				    );

// Reject stamps, based on the solution; return TRUE if rejection performed
bool pmSubtractionRejectStamps(psArray *stamps, // Array of stamps to check for rejection
			       psImage *mask, // Mask image
			       unsigned int badStampMaskVal, // Value to use in mask for bad stamp
			       int footprint, // Region to mask if stamp is bad
			       float sigmaRej, // Number of RMS deviations above zero at which to reject
			       const psImage *refImage, // Reference image
			       const psImage *inImage, // Input image
			       const psVector *solution, // Solution vector
			       const psArray *kernelParams // Array of kernel parameters
			       );

// Create a kernel image from the solution vector
psImage *pmSubtractionKernelImage(psImage *out,	// Output image, or NULL
				  const psVector *solution, // Vector containing the kernel elements
				  const psArray *kernelParams, // The kernel parameters
				  float x, float y // The relative (-1 to 1) x position for which to get the
						   // kernel
				  );


// The whole image subtraction process looks something like this:
psImage *subtractImages(psImage **subtractionMask, // Mask for subtraction (output)
			const psImage *reference, // Reference image
			const psImage *refMask, // Mask for reference image
			const psImage *input,	// Input image
			const psImage *inMask,// Mask for input image
			unsigned int maskVal, // Value to be masked
			pmSubtractionKernelType kernelType, // Type of kernel
			int kernelHalfSize, // Half the kernel size (full size is 2*kernelHalfSize + 1)
			const psVector *sigmas, // Sigmas for ISIS Gaussians
			const psVector *polyOrders, // Polynomial orders for ISIS Gaussians
			int spatialOrder, // Maximum spatial order
			float stampThreshold, // Threshold for stamps
			int nStampsX, int nStampsY, // Number of stamps in x and y
			int stampSize,	// Half size of stamp footprint
			int numIter,	// Number of iterations on the stamps
			float sigmaRej	// Rejection threshold
			)
{
    // Mask around bad pixels in the reference image.  There are two cases to worry about:
    // 1. Bad pixels within the kernel, which will affect the subtracted image
    // 2. Bad pixels within the stamp, which affects the calculation of the kernel
    psImage *subMask = psImageGrowMask(NULL, refMask, maskVal, kernelHalfSize, PS_MASK_NEAR_BAD);
    (void)psImageGrowMask(subMask, refMask, maskVal, stampSize, PS_MASK_BAD_STAMP);
    // Add in the mask for the input image.  Don't need to grow this, since it isn't convolved.
    (void)psBinaryOp(subMask, subMask, "|", inMask);

    // Generate kernel basis functions
    psArray *kernels = NULL;		// Array of kernel basis functions
    switch (kernelType) {
      case PM_SUBTRACTION_KERNEL_POIS:
	// Create the kernel basis functions
	kernels = pmSubtractionKernelsGeneratePOIS(kernelHalfSize, spatialOrder);
	break;
      case PM_SUBTRACTION_KERNEL_ISIS:
	kernels = pmSubtractionKernelsGenerateISIS(sigmas, polyOrders, kernelHalfSize, spatialOrder);
	break;
      default:
	barf();
    }

    psArray *stamps = NULL;		// Array of stamps
    psVector *kernelCoeffs = NULL;	// Coefficients for the kernels
    bool rejected = true;		// Did we reject a stamp in the last iteration?

    // Iterate for a solution
    for (int iter = 0; iter < numIter && rejected; iter++) {

	// Find stamps
	stamps = pmSubtractionFindStamps(stamps, reference, subMask, maskVal | PS_MASK_BAD_STAMP,
					 stampThreshold, nStampsX, nStampsY, stampSize, kernelHalfSize);

	// Generate and solve matrix equations
	(void)pmSubtractionCalculateEquation(stamps, reference, input, kernels, stampSize);
	kernelCoeffs = pmSubtractionSolveEquation(kernelCoeffs, stamps);

	// Reject bad stamps
	rejected = pmSubtractionRejectStamps(stamps, subMask, PS_MASK_BAD_STAMP, stampSize, sigmaRej,
					     reference, input, kernelCoeffs, kernels);
    }

    // Convolve the reference image
    psImage *referenceConvolved = pmSubtractionConvolveImage(NULL, reference, subMask, kernelCoeffs, kernels);
    // Subtract
    psImage *subtracted = (psImage*)psBinaryOp(NULL, input, "-", referenceConvolved);

    // What does the kernel look like?
    psImage *kernelImage = pmSubtractionKernelImage(NULL, kernelCoeffs, kernels, 0.0, 0.0);
    // Check/save kernel image, print statistics....

    psFree(referenceConvolved);
    psFree(stamps);
    psFree(kernels);
    psFree(kernelCoeffs);

    *subtractionMask = subMask;
    return subtracted;
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// This stuff is still under development.  We don't need it immediately, but need to think about it for the //
// long run.                                                                                                //
//////////////////////////////////////////////////////////////////////////////////////////////////////////////


// A tangent plane on the sky, which may be a weird shape (triangle, hexagon, silly healpix thing).
typedef struct {
    double ra, dec;			// RA and Dec of tangent point
    double plateScale;			// Plate scale (arcsec/pixel), identical in x and y
    psArray *rows;			// The data, stored as an array of vectors
    psVector *offsets;			// Integer offset for row
} psSkyPlane;

psImage *psSkyPlaneToImage(psImage *out,	// Output image, or NULL
			   psImage **mask, // Output mask image, or NULL
			   const psSkyPlane *sky // Input sky
			   );

psSkyPlane *psImageToSkyPlane(psSkyPlane *out, // Output sky, or NULL
			      const psImage *image, // Input image
			      const psImage *mask	// Mask image, or NULL
			      );

// Return sky planes within a region on the sky.
// Could be broken into two logical steps: identify sky planes, then read/copy/retrieve
psArray *psSkyPlaneRetrieve(double ra, double dec, // Coordinates around which to search
			    int raSize, int decSize // Size around which to search
			    );

// Rough
astrometry *psImageSetWCSFromTP(psImage *image, double ra, double dec, double plateScale);
