Index: /trunk/archive/modules/include/phase4stuff.h
===================================================================
--- /trunk/archive/modules/include/phase4stuff.h	(revision 3433)
+++ /trunk/archive/modules/include/phase4stuff.h	(revision 3433)
@@ -0,0 +1,369 @@
+// 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);