Index: /trunk/pois/Makefile.am
===================================================================
--- /trunk/pois/Makefile.am	(revision 5716)
+++ /trunk/pois/Makefile.am	(revision 5717)
@@ -1,3 +1,3 @@
-SUBDIRS = src swig
+SUBDIRS = src #swig
 
 CLEANFILES = *.pyc *~ core core.*
Index: /trunk/pois/configure.ac
===================================================================
--- /trunk/pois/configure.ac	(revision 5716)
+++ /trunk/pois/configure.ac	(revision 5717)
@@ -15,12 +15,12 @@
 AC_PROG_LIBTOOL
 
-PKG_CHECK_MODULES(PSLIB, pslib >= 0.5.0) 
+PKG_CHECK_MODULES(PSLIB, pslib >= 0.8.0) 
 
-dnl is this the best was to setup recursive CFLAGS?
-pois_CFLAGS="-Wall -g -O3 -std=c99"
+dnl Set CFLAGS for build
+pois_CFLAGS="-Wall -Werror -g -O3 -std=c99"
 AC_SUBST([pois_CFLAGS])
 
 AC_C_BIGENDIAN([PS_BIGENDIAN=1],[PS_BIGENDIAN=0])
-		
+
 AC_PROG_SWIG(1.3.24)
 SWIG_MULTI_MODULE_SUPPORT
@@ -29,5 +29,5 @@
 dnl lifted from autoconf 2.59's configure.ac
 AC_PATH_PROG([PERL], perl, no)
-AC_SUBST([PERL])dnl
+AC_SUBST([PERL])
 if test "$PERL" = no; then
   AC_MSG_ERROR([perl is not found])
Index: /trunk/pois/src/Makefile.am
===================================================================
--- /trunk/pois/src/Makefile.am	(revision 5716)
+++ /trunk/pois/src/Makefile.am	(revision 5717)
@@ -1,4 +1,7 @@
-AM_CFLAGS = $(pois_CFLAGS) $(PSLIB_CFLAGS)
+AM_CFLAGS = $(PSLIB_CFLAGS) $(pois_CFLAGS)
 AM_LDFLAGS = $(PSLIB_LIBS)
+
+pois_CPPFLAGS = $(PSLIB_CFLAGS) $(pois_CFLAGS)
+pois_LDFLAGS = $(PSLIB_LIBS)
 
 include_HEADERS = pois.h poisErrorCodes.h
@@ -19,18 +22,20 @@
 libpois_la_SOURCES = \
 	pois.h \
-	poisParseConfig.c \
-	poisKernelBasisFunctions.c \
+	poisAddPenalty.c \
+	poisCalculateDeviations.c \
 	poisCalculateEquation.c \
+	poisCheckKernel.c \
 	poisConvolveImage.c \
 	poisExtractKernel.c \
+	poisFindStamps.c \
+	poisKernelBasisFunctions.c \
 	poisMakeMask.c \
-	poisFindStamps.c \
-	poisCalculateDeviations.c \
-	poisAddPenalty.c \
+	poisReadStamps.c \
+	poisRejectStamps.c \
+	poisParseConfig.c \
+	poisSolveEquation.c \
 	poisStamp.c \
-	poisRejectStamps.c \
-	poisSolveEquation.c \
-	poisMaskOps.c \
 	errorCodes.c
+###	poisMaskOps.c
 
 #
@@ -49,2 +54,6 @@
 	-$(RM) testout.fits stamp*.fits kernel*.fits
 	./pois -v test1.fits test2.fits testout.fits
+
+# Tags for emacs
+tags:
+	etags `find . -name \*.[ch] -print`
Index: /trunk/pois/src/pois.c
===================================================================
--- /trunk/pois/src/pois.c	(revision 5716)
+++ /trunk/pois/src/pois.c	(revision 5717)
@@ -7,4 +7,5 @@
 
 #include <stdio.h>
+#include <string.h>
 #include <sys/time.h>
 #include "pslib.h"
@@ -25,34 +26,44 @@
 int main(int argc, char **argv)
 {
-    poisErrorRegister();
-
     double startTime = getTime();
 
     // Set trace levels
     psTraceSetLevel(".", 0);
-    psTraceSetLevel("pois", 10);
-    psTraceSetLevel("pois.config", 10);
-    psTraceSetLevel("pois.time", 10);
-    psTraceSetLevel("pois.solution", 0);
-    psTraceSetLevel("pois.kernelBasisFunctions", 10);
-    psTraceSetLevel("pois.calculateEquation", 10);
-    psTraceSetLevel("pois.convolveImage", 10);
-    psTraceSetLevel("pois.findStamps", 10);
-    psTraceSetLevel("pois.parseConfig", 10);
-    psTraceSetLevel("pois.makeMask", 10);
-    psTraceSetLevel("pois.extractKernel", 10);
-    psTraceSetLevel("pois.calculateDeviations", 10);
-
-    // Set logging level
-    psLogSetLevel(9);
 
     // Parse the command line
     poisConfig *config = poisParseConfig(argc, argv); // Configuration values
 
+    if (config->verbose) {
+	psTraceSetLevel("pois", 10);
+	psTraceSetLevel("pois.config", 10);
+	psTraceSetLevel("pois.time", 10);
+	psTraceSetLevel("pois.solution", 0);
+	psTraceSetLevel("pois.kernelBasisFunctions", 10);
+	psTraceSetLevel("pois.calculateEquation", 10);
+	psTraceSetLevel("pois.convolveImage", 10);
+	psTraceSetLevel("pois.findStamps", 10);
+	psTraceSetLevel("pois.parseConfig", 10);
+	psTraceSetLevel("pois.makeMask", 10);
+	psTraceSetLevel("pois.extractKernel", 10);
+	psTraceSetLevel("pois.calculateDeviations", 10);
+	psTraceSetLevel("pois.checkKernel", 10);
+	if (config->stampFile) {
+	    psTraceSetLevel("pois.checkStamp", 10);
+	} else {
+	    psTraceSetLevel("pois.checkStamp", 0);
+	}
+	// Set logging level
+	psLogSetLevel(9);
+    }
+
     // Read inputs
-    psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0);
+    psMetadata *header = NULL;		// Header for output image
+
+    psRegion imageRegion = {0, 0, 0, 0};
     psFits *reference = psFitsAlloc(config->refFile);
-    //psMetadata *refHeader = psFitsReadHeader(NULL, reference);
-    psImage *refImage = psFitsReadImage(NULL, reference, *imageRegion, 0);
+    if (config->reverse) {
+	header = psFitsReadHeader(NULL, reference);
+    }
+    psImage *refImage = psFitsReadImage(NULL, reference, imageRegion, 0);
     if (refImage == NULL) {
 	psErrorStackPrint(stderr, "Fatal error: unable to read %s\n", config->refFile);
@@ -61,8 +72,17 @@
     psTrace("pois", 3, "Read %s: %dx%d\n", config->refFile, refImage->numCols, refImage->numRows);
     psFree(reference);
+    // Convert to 32-bit floating point, in necessary
+    if (refImage->type.type != PS_TYPE_F32) {
+	psTrace("pois", 3, "Converting %s to floating point in memory....\n", config->refFile);
+	psImage *temp = psImageCopy(NULL, refImage, PS_TYPE_F32);
+	psFree(refImage);
+	refImage = temp;
+    }
 
     psFits *input = psFitsAlloc(config->inFile);
-    //psMetadata *inHeader = psFitsReadHeader(NULL, input);
-    psImage *inImage = psFitsReadImage(NULL, input, *imageRegion, 0);
+    if (! config->reverse) {
+	header = psFitsReadHeader(NULL, input);
+    }
+    psImage *inImage = psFitsReadImage(NULL, input, imageRegion, 0);
     if (inImage == NULL) {
 	psErrorStackPrint(stderr, "Fatal error: unable to read %s\n", config->inFile);
@@ -71,5 +91,11 @@
     psTrace("pois", 3, "Read %s: %dx%d\n", config->inFile, inImage->numCols, inImage->numRows);
     psFree(input);
-    psFree(imageRegion);
+    // Convert to 32-bit floating point, in necessary
+    if (inImage->type.type != PS_TYPE_F32) {
+	psTrace("pois", 3, "Converting %s to floating point in memory....\n", config->inFile);
+	psImage *temp = psImageCopy(NULL, inImage, PS_TYPE_F32);
+	psFree(inImage);
+	inImage = temp;
+    }
 
     if ((refImage->numCols != inImage->numCols) || (refImage->numRows != inImage->numRows)) {
@@ -102,5 +128,5 @@
     // Write the mask out
     psFits *maskFile = psFitsAlloc("mask.fits");
-    if (!psFitsWriteImage(maskFile, NULL, mask, 0, NULL)) {
+    if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
 	psErrorStackPrint(stderr, "Unable to write mask: mask.fits\n");
     }
@@ -110,16 +136,27 @@
 
     psArray *stamps = NULL;		// Array of stamps
-    psVector *stampMask = NULL;		// Mask for stamps
     psVector *solution = NULL;		// Solution vector
 
+    FILE *stampsFP = NULL;		// File pointer for stamps
+    if (config->stampFile) {
+	psTrace("pois", 5, "Opening stamp file, %s\n", config->stampFile);
+	stampsFP = fopen(config->stampFile, "r");
+	if (! stampsFP) {
+	    psError(PS_ERR_IO, true, "Unable to open stamps file, %s!\n", config->stampFile);
+	    exit(EXIT_FAILURE);
+	}
+    }
+
     // Iterate for a good solution
-    psList *badStamps = NULL;		// Do we have bad stamps, such that we need to continue to iterate?
+    bool badStamps = true;		// Do we have bad stamps, such that we need to continue to iterate?
     for (int iterNum = 0; iterNum < config->numIter && badStamps; iterNum++) {
 	psTrace("pois", 1, "Iteration %d...\n", iterNum);
 
-	psFree(badStamps);		// left over from last iteration
-	
 	// Find stamps
-	stamps = poisFindStamps(stamps, refImage, mask, config);
+	if (config->stampFile) {
+	    stamps = poisReadStamps(stamps, &stampsFP, refImage, mask, config);
+	} else {
+	    stamps = poisFindStamps(stamps, refImage, mask, config);
+	}
 	int numStamps = 0;
 	for (int s = 0; s < stamps->n; s++) {
@@ -159,5 +196,5 @@
 	psImage *kernelImage = poisExtractKernel(solution, kernelBasisFunctions, 0.0, 0.0, config);
 	psFits *kernelFile = psFitsAlloc(kernelName);
-	if (!psFitsWriteImage(kernelFile, NULL, kernelImage, 0, NULL)) {
+	if (!psFitsWriteImage(kernelFile, NULL, kernelImage, 0)) {
 	    psErrorStackPrint(stderr, "Unable to write kernel: %s\n", kernelName);
 	}
@@ -177,6 +214,9 @@
     // If there was rejection on the last round, re-solve the equation using only the good stamps
     if (badStamps) {
-	psFree(badStamps);
 	solution = poisSolveEquation(solution, stamps, config);
+    }
+
+    if (stampsFP) {
+	fclose(stampsFP);
     }
 
@@ -185,8 +225,11 @@
 	for (int i = 0; i < mask->numCols; i++) {
 	    if (mask->data.U8[j][i] & POIS_MASK_STAMP) {
-		mask->data.U8[j][i] = POIS_MASK_OK;
-	    }
-	}
-    }
+		mask->data.U8[j][i] &= ~POIS_MASK_STAMP;
+	    }
+	}
+    }
+
+    // Check the result
+    (void)poisCheckKernel(solution, kernelBasisFunctions, config);
 
     // Convolve the input image
@@ -199,5 +242,5 @@
     snprintf(convName, MAXCHAR, "%s.conv", config->outFile);
     psFits *convFile = psFitsAlloc(convName);
-    if (!psFitsWriteImage(convFile, NULL, convImage, 0, NULL)) {
+    if (!psFitsWriteImage(convFile, NULL, convImage, 0)) {
 	psErrorStackPrint(stderr, "Unable to write convolved image: %s\n", convName);
     }
@@ -208,5 +251,9 @@
     // Do the subtraction
     psImage *subImage = psImageAlloc(config->xImage, config->yImage, PS_TYPE_F32);
-    (void)psBinaryOp(subImage, inImage, "-", convImage);
+    if (config->reverse) {
+	(void)psBinaryOp(subImage, convImage, "-", inImage);
+    } else {
+	(void)psBinaryOp(subImage, inImage, "-", convImage);
+    }
     psTrace("pois.time", 1, "Subtraction completed at %f sec\n", getTime() - startTime);
 
@@ -231,16 +278,16 @@
 	    }
 	}
-	
-	for (int y = 0; y < config->yKernel; y++) {
-	    for (int x = 0; x < subImage->numCols; x++) {
-		subImage->data.F32[y][x] = 0.0;
-	    }
-	}
-	for (int y = subImage->numRows - config->yKernel; y < subImage->numRows; y++) {
-	    for (int x = 0; x < subImage->numCols; x++) {
-		subImage->data.F32[y][x] = 0.0;
-	    }
-	}
-    }
+    }
+    for (int y = 0; y < config->yKernel; y++) {
+	for (int x = 0; x < subImage->numCols; x++) {
+	    subImage->data.F32[y][x] = 0.0;
+	}
+    }
+    for (int y = subImage->numRows - config->yKernel; y < subImage->numRows; y++) {
+	for (int x = 0; x < subImage->numCols; x++) {
+	    subImage->data.F32[y][x] = 0.0;
+	}
+    }
+
 
     // Mask out bad pixels
@@ -253,6 +300,23 @@
     }
 
+    if (config->mask) {
+	char maskName[80];		// Name of mask image
+	sprintf(maskName, "%s.mask", config->outFile);
+	psFits *maskFile = psFitsAlloc(maskName);
+        if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
+            psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName);
+        }
+        psTrace("pois", 1, "Mask image written to %s\n", maskName);
+        psFree(maskFile);
+    }
+
+
+    int numNaN = psImageClipNaN(subImage, 0.0);
+    if (numNaN > 0) {
+	printf("Masked %d NAN pixels to zero.\n", numNaN);
+    }
+
     psFits *subFile = psFitsAlloc(config->outFile);
-    if (!psFitsWriteImage(subFile, NULL, subImage, 0, NULL)) {
+    if (!psFitsWriteImage(subFile, header, subImage, 0)) {
 	psErrorStackPrint(stderr, "Unable to write subtracted image: %s\n", config->outFile);
     }
@@ -261,11 +325,12 @@
     psTrace("pois.time", 1, "I/O completed at %f sec\n", getTime() - startTime);
 
-
     // Clean up
-    (void)psFree(kernelBasisFunctions);
-    (void)psFree(refImage);
-    (void)psFree(inImage);
-    (void)psFree(convImage);
-    (void)psFree(subImage);
-    (void)psFree(config);
+    psFree(header);
+    psFree(kernelBasisFunctions);
+    psFree(refImage);
+    psFree(inImage);
+    psFree(convImage);
+    psFree(subImage);
+    psFree(config);
 }
+
Index: /trunk/pois/src/pois.h
===================================================================
--- /trunk/pois/src/pois.h	(revision 5716)
+++ /trunk/pois/src/pois.h	(revision 5717)
@@ -3,5 +3,4 @@
 
 #include "pslib.h"
-#include "poisErrorCodes.h"
 
 #define abs(x) ((x) >= 0 ? (x) : (-(x)))
@@ -21,5 +20,5 @@
     int nsx, nsy;			// Number of stamps in x and y
     float refSat, inSat;		// Saturation values for reference and input images
-    float bad;				// Value below which pixels are bad
+    float refBad, inBad;		// Bad values for reference and input images
     float background;			// Value to add to the background in images to get above zero
     int spatialOrder;			// Order of spatial variations in the kernel
@@ -27,4 +26,7 @@
     float sigmaRej;			// Limit (in std dev) for rejection
     float penalty;			// Penalty value
+    char *stampFile;			// Filename for file with stamps
+    bool reverse;			// Reverse the subtraction?
+    bool mask;				// Output the mask?
 } poisConfig;
 
@@ -64,4 +66,11 @@
 void poisStampFree(poisStamp *stamp);
 
+// Return true if the stamp is OK
+bool poisCheckStamp(const psImage *image, // Image to check for threshold
+		    const psImage *mask, // Mask to check for bad pixels
+		    int x, int y,	// Pixel coordinates
+		    const poisConfig *config // Configuration
+    );
+
 /************************************************************************************************************/
 
@@ -78,9 +87,7 @@
 
 // Print usage information
-//void help(void);
+void help(void);
 
 // Parse the command line arguments
-poisConfig *restrict poisConfigAlloc(void);
-
 poisConfig *poisParseConfig(int argc,	// Number of command-line arguments
 			    char **argv	// Command-line arguments
@@ -111,5 +118,5 @@
 
 // Find stamps
-psArray *poisFindStamps(psArray *stampsIO, // Existing list of stamps, or NULL
+psArray *poisFindStamps(psArray *stamps, // Existing list of stamps, or NULL
 			const psImage *image, // Image for which to find stamps
 			const psImage *mask, // Mask image
@@ -158,5 +165,5 @@
 
 // Calculate deviations from the best fit for the stamps
-psVector *poisCalculateDeviations(psVector *deviationsIO, // Output array of deviations, or NULL
+psVector *poisCalculateDeviations(psVector *deviations,	// Output array of deviations, or NULL
 				  psArray *stamps, // Array of stamps
 				  psImage *refImage, // Reference image
@@ -173,8 +180,8 @@
 
 // Reject stamps, based on the deviations
-psList *restrict poisRejectStamps(psArray *stamps,	// Array of stamps
-				  psImage *mask,	// Mask image
-				  const psVector *deviations, // Vector of deviations for the stamps
-				  const poisConfig *config // Configuration
+bool poisRejectStamps(psArray *stamps,	// Array of stamps
+		      psImage *mask,	// Mask image
+		      const psVector *deviations, // Vector of deviations for the stamps
+		      const poisConfig *config // Configuration
     );
 
@@ -184,5 +191,5 @@
 
 // Solve the matrix equation
-psVector *poisSolveEquation(psVector *solutionIO,	// Solution vector, or NULL
+psVector *poisSolveEquation(psVector *solution,	// Solution vector, or NULL
 			    const psArray *stamps, // Array of stamps
 			    const poisConfig *config
@@ -203,14 +210,25 @@
 /************************************************************************************************************/
 
-psImage *poisImageSetVal(psImage *restrict IO,	     // Input image
-			 const psC64 val	     // set to this value
-    );
-
-psImage *poisImageSetValInMask(psImage *IO,			// Image to update, or NULL
-			       const psImage *restrict in,	// Input image
-			       const psImage *mask,		// mask for image
-			       const psC64 val,			// set to this value
-			       const unsigned long bits		// set pixels where mask & bits != 0
-    );
+// poisReadStamps.c
+
+// Read stamps from a file
+psArray *poisReadStamps(psArray *stamps,// Stamps
+			FILE **stampFP, // File pointer for stamps file
+			const psImage *image, // Image for which to find stamps
+			const psImage *mask, // Mask image
+			const poisConfig *config // Configuration
+    );
+
+/************************************************************************************************************/
+
+// poisCheckKernel.c
+
+// Check kernel for power at large radii
+bool poisCheckKernel(const psVector *solution,// Kernel solution
+		     const psArray *kernels, // Kernel basis functions
+		     const poisConfig *config // Configuration
+    );
+
+/************************************************************************************************************/
 
 #endif
Index: /trunk/pois/src/poisCalculateDeviations.c
===================================================================
--- /trunk/pois/src/poisCalculateDeviations.c	(revision 5716)
+++ /trunk/pois/src/poisCalculateDeviations.c	(revision 5717)
@@ -46,18 +46,17 @@
 	int y = stamp->y;		// Stamp y coord
 	if (stamp->status == POIS_STAMP_USED) {
-	    psImage *refStamp = psImageSubset(refImage, x - xSize, y - ySize, x + xSize, y + ySize);
-	    psImage *inStamp = psImageSubset(inImage, x - xSize, y - ySize, x + xSize, y + ySize);
-	    psImage *maskStamp = psImageSubset(mask, x - xSize, y - ySize, x + xSize, y + ySize);
+	    psRegion stampRegion = {x - xSize, x + xSize, y - ySize, y + ySize};
+	    psImage *refStamp = psImageSubset(refImage, stampRegion);
+	    psImage *inStamp = psImageSubset(inImage, stampRegion);
+	    psImage *maskStamp = psImageSubset(mask, stampRegion);
 	    psImage *convRefStamp = poisConvolveImage(refStamp, maskStamp, solution, kernelParams, config);
 	    // Calculate chi^2
 	    (void)psBinaryOp(subStamp, inStamp, "-", convRefStamp);
+	    (void)psBinaryOp(subStamp, subStamp, "*", subStamp);
 	    (void)psBinaryOp(subStamp, subStamp, "/", inStamp);
-	    (void)psBinaryOp(subStamp, subStamp, "*", subStamp);
-	    psImage *subStampTrim = psImageSubset(subStamp, config->xKernel, config->yKernel,
-						  config->xKernel + 2 * footprint,
-						  config->yKernel + 2 * footprint);
-	    psImage *maskStampTrim = psImageSubset(maskStamp, config->xKernel, config->yKernel,
-						   config->xKernel + 2 * footprint,
-						   config->yKernel + 2 * footprint);
+	    psRegion stampTrim = { config->xKernel, config->xKernel + 2 * footprint, config->yKernel,
+				   config->yKernel + 2 * footprint };
+	    psImage *subStampTrim = psImageSubset(subStamp, stampTrim);
+	    psImage *maskStampTrim = psImageSubset(maskStamp, stampTrim);
 	    // Copy image to workaround bug 305
 	    psImage *tempImage = psImageCopy(NULL, subStampTrim, PS_TYPE_F32);
@@ -69,5 +68,5 @@
 	    psFree(tempMask);
 	    
-	    deviations->data.F32[s] = stats->sampleMean * (float)footprint * (float)footprint * 4.0;
+	    deviations->data.F32[s] = sqrtf(stats->sampleMean / 2.0);
 	    psTrace("pois.calculateDeviations", 5, "Deviation of stamp %d (%d,%d) is %f\n", s, x, y,
 		    deviations->data.F32[s]);
@@ -77,5 +76,5 @@
 	    snprintf(stampName, MAXCHAR, "stamp%d.fits", s);
 	    psFits *stampFile = psFitsAlloc(stampName);
-	    if (!psFitsWriteImage(stampFile, NULL, subStampTrim, 0, NULL)) {
+	    if (!psFitsWriteImage(stampFile, NULL, subStampTrim, 0)) {
 		psErrorStackPrint(stderr, "Unable to write stamp: %s\n", stampName);
 	    }
@@ -83,4 +82,5 @@
 #endif
 
+#if 0
 	    psFree(convRefStamp);
 	    psFree(maskStampTrim);
@@ -89,4 +89,5 @@
 	    psFree(refStamp);
 	    psFree(inStamp);
+#endif
 	}
     }
Index: /trunk/pois/src/poisCalculateEquation.c
===================================================================
--- /trunk/pois/src/poisCalculateEquation.c	(revision 5716)
+++ /trunk/pois/src/poisCalculateEquation.c	(revision 5717)
@@ -7,8 +7,8 @@
 /* Calculate the matrix and vector for the matrix equation */
 int poisCalculateEquation(psArray *stamps, // The stamps
-			  const psImage *refImage, // The reference image
-			  const psImage *inImage, // The input image
-			  const psArray *kernelParams, // Kernel parameters
-			  const poisConfig *config // Configuration
+                          const psImage *refImage, // The reference image
+                          const psImage *inImage, // The input image
+                          const psArray *kernelParams, // Kernel parameters
+                          const poisConfig *config // Configuration
     )
 {
@@ -23,150 +23,148 @@
     // in the background levels.
     int numSolveParams = kernelParams->n + 1;
-    int bgIndex = kernelParams->n;	// Index in matrix for the background (it's the last one, after the
-					// kernel parameters)
+    int bgIndex = kernelParams->n;      // Index in matrix for the background (it's the last one, after the
+                                        // kernel parameters)
     float **reference = refImage->data.F32; // The reference pixels
-    float **input = inImage->data.F32;	// The input pixels
-    int nPix = 0;			// Number of pixels used in calculation (for interest)
+    float **input = inImage->data.F32;  // The input pixels
+    int nPix = 0;                       // Number of pixels used in calculation (for interest)
     float nx = 0.5 * (float)refImage->numCols; // Half-size of input image in x
     float ny = 0.5 * (float)refImage->numRows; // Half-size of input image in y
     int numKernelParams = kernelParams->n; // Number of kernel parameters
 
-    psDPolynomial2D *poly = psDPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
-						 PS_POLYNOMIAL_ORD); // Polynomial for evaluation
+    psPolynomial2D *poly = psPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
+                                               PS_POLYNOMIAL_ORD); // Polynomial for evaluation
     for (int j = 0; j < config->spatialOrder + 1; j++) {
-	for (int i = 0; i < config->spatialOrder + 1; i++) {
-	    poly->coeff[j][i] = 1.0;
-	    poly->mask[j][i] = 1;	// Mask all coefficients; unmask to evaluate
-	}
+        for (int i = 0; i < config->spatialOrder + 1; i++) {
+            poly->coeff[j][i] = 1.0;
+            poly->mask[j][i] = 1;       // Mask all coefficients; unmask to evaluate
+        }
     }
 
     psImage *polyValues = psImageAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
-				       PS_TYPE_F64); // Evaluations for each of the polynomial orders
+                                       PS_TYPE_F64); // Evaluations for each of the polynomial orders
 
     /* Iterate over the stamps */
     for (int s = 0; s < stamps->n; s++) {
-	poisStamp *stamp = stamps->data[s]; // The stamp
+        poisStamp *stamp = stamps->data[s]; // The stamp
 
-	if (stamp->status == POIS_STAMP_RECALC) {
-	    psTrace("pois.calculateEquation", 5, "Calculating for stamp %d: %d,%d\n", s, stamp->x, stamp->y);
+        if (stamp->status == POIS_STAMP_RECALC) {
+            psTrace("pois.calculateEquation", 5, "Calculating for stamp %d: %d,%d\n", s, stamp->x, stamp->y);
 
-	    // Initialise the matrix and vector if required
-	    if (! stamp->matrix) {
-		psTrace("pois.calculateEquation", 9, "Allocating matrix: %dx%d\n", numSolveParams,
-			numSolveParams);
-		psImage *matrix = psImageAlloc(numSolveParams, numSolveParams, PS_TYPE_F64);
-		for (int j = 0; j < numSolveParams; j++) {
-		    for (int i = 0; i < numSolveParams; i++) {
-			matrix->data.F64[j][i] = 0.0;
-		    }
-		}
-		stamp->matrix = matrix;
-	    }
-	    if (! stamp->vector) {
-		psTrace("pois.calculateEquation", 9, "Allocating vector: %d\n", numSolveParams);
-		psVector *vector = psVectorAlloc(numSolveParams, PS_TYPE_F64);
-		for (int i = 0; i < numSolveParams; i++) {
-		    vector->data.F64[i] = 0.0;
-		}
-		stamp->vector = vector;
-	    }
-	    
-	    // Dereference, for convenience
-	    psImage *matrix = stamp->matrix;
-	    psVector *vector = stamp->vector;
-	    
-	    // Will assert on the type, since can't do much about it if it's wrong.
-	    assert(matrix->type.type == PS_TYPE_F64);
-	    assert(vector->type.type == PS_TYPE_F64);
-	    
-	    // Evaluate the polynomial for each order
-	    for (int j = 0; j < config->spatialOrder + 1; j++) {
-		for (int i = 0; i < config->spatialOrder + 1 - j; i++) {
-		    poly->mask[j][i] = 0;
-		    polyValues->data.F64[j][i] = psDPolynomial2DEval(poly, ((double)stamp->x - nx)/nx,
-								     ((double)stamp->y - ny)/ny);
-		    poly->mask[j][i] = 1;
-		}
-	    }
-	    
-	    // Iterate over the pixels
-	    for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) {
-		for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) {
-		    float invNoise2 = 1.0/reference[y][x]; // The inverse of the noise, squared.
-		    nPix++;
-		    
-		    /* Iterate over the first convolution */
-		    for (int k1 = 0; k1 < numKernelParams; k1++) {
-			poisKernelBasis *kernel1 = kernelParams->data[k1]; // Kernel basis #1
-			
-			int u1 = kernel1->u; // Offset in x
-			int v1 = kernel1->v; // Offset in y
-			int i1 = kernel1->xOrder; // Polynomial order in x
-			int j1 = kernel1->yOrder; // Polynomial order in y
-			// First convolution
-			float conv1 = polyValues->data.F64[j1][i1] * reference[y - v1][x - u1];
+            // Initialise the matrix and vector if required
+            if (! stamp->matrix) {
+                psTrace("pois.calculateEquation", 9, "Allocating matrix: %dx%d\n", numSolveParams,
+                        numSolveParams);
+                stamp->matrix = psImageAlloc(numSolveParams, numSolveParams, PS_TYPE_F64);
+            }
+            if (! stamp->vector) {
+                psTrace("pois.calculateEquation", 3, "Allocating vector: %d\n", numSolveParams);
+                stamp->vector = psVectorAlloc(numSolveParams, PS_TYPE_F64);
+            }
 
-			// Assuming that the first kernel component is 0 order in x and y, and 0 offset
-			if (k1 != 0) {
-			    conv1 -= reference[y][x];
-			}
-			
-			/* Iterate over the second convolution */
-			for (int k2 = k1; k2 < numKernelParams; k2++) {
-			    poisKernelBasis *kernel2 = kernelParams->data[k2]; // Kernel basis #2
-			    
-			    int u2 = kernel2->u; // Offset in x
-			    int v2 = kernel2->v; // Offset in y
-			    int i2 = kernel2->xOrder; // Polynomial order in x
-			    int j2 = kernel2->yOrder; // Polynomial order in y
-			    // Second convolution
-			    float conv2 = polyValues->data.F64[j2][i2] * reference[y - v2][x - u2];
-			    
-			    // Assuming that the first kernel component is 0 order in x and y, and 0 offset
-			    if (k2 != 0) {
-				conv2 -= reference[y][x];
-			    }
+            // Dereference, for convenience
+            psImage *matrix = stamp->matrix;
+            psVector *vector = stamp->vector;
 
-			    // Add into the matrix element
-			    matrix->data.F64[k1][k2] += conv1 * // First convolution
-				conv2 * // Second convolution
-				invNoise2; // Divide by sigma^2, ~ poisson.
-			} // Second convolution
-			
-			// Add into the vector element
-			vector->data.F64[k1] += input[y][x] * // Input image, no convolution
-			    conv1 *	// Convolved reference image
-			    invNoise2; // Divide by sigma^2, ~ poisson.
-			
-			/* Background term */
-			matrix->data.F64[k1][bgIndex] += conv1 * // Convolved reference image
-			    invNoise2; // Divide by sigma^2, ~ poisson.
-			
-		    } // First convolution
-		    
-		    /* Background only terms */
-		    matrix->data.F64[bgIndex][bgIndex] += invNoise2;
-		    vector->data.F64[bgIndex] += input[y][x] * invNoise2;
-		}
-	    } // Iterating over pixels in stamp
+            // Will assert on the type, since can't do much about it if it's wrong.
+            assert(matrix->type.type == PS_TYPE_F64);
+            assert(vector->type.type == PS_TYPE_F64);
 
-	    // Fill in other side of symmetric matrix
-	    for (int k1 = 0; k1 < kernelParams->n; k1++) {
-		for (int k2 = 0; k2 < k1; k2++) {
-		    matrix->data.F64[k1][k2] = matrix->data.F64[k2][k1];
-		}
-		matrix->data.F64[bgIndex][k1] = matrix->data.F64[k1][bgIndex];
-	    }
+            // Initialise the matrix and vector
+            for (int j = 0; j < numSolveParams; j++) {
+                for (int i = 0; i < numSolveParams; i++) {
+                    matrix->data.F64[j][i] = 0.0;
+                }
+                vector->data.F64[j] = 0.0;
+            }
 
-	    for (int k = 0; k < kernelParams->n; k++) {
-		poisKernelBasis *kernel = kernelParams->data[k]; // Kernel basis #2
-		int u = kernel->u; // Offset in x
-		int v = kernel->v; // Offset in y
-		matrix->data.F64[k][k] += config->penalty * (float)(u*u + v*v);
-	    }
+            // Evaluate the polynomial for each order
+            for (int j = 0; j < config->spatialOrder + 1; j++) {
+                for (int i = 0; i < config->spatialOrder + 1 - j; i++) {
+                    poly->mask[j][i] = 0;
+                    polyValues->data.F64[j][i] = psPolynomial2DEval(poly, ((double)stamp->x - nx)/nx,
+                                                                    ((double)stamp->y - ny)/ny);
+                    poly->mask[j][i] = 1;
+                }
+            }
 
-	    stamp->status = POIS_STAMP_USED; // We've calculated it now, so don't need to do so again.
+            // Iterate over the pixels
+            for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) {
+                for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) {
+                    float invNoise2 = 1.0/reference[y][x]; // The inverse of the noise, squared.
+                    nPix++;
 
-	} // Legitimate stamps
+                    /* Iterate over the first convolution */
+                    for (int k1 = 0; k1 < numKernelParams; k1++) {
+                        poisKernelBasis *kernel1 = kernelParams->data[k1]; // Kernel basis #1
+
+                        int u1 = kernel1->u; // Offset in x
+                        int v1 = kernel1->v; // Offset in y
+                        int i1 = kernel1->xOrder; // Polynomial order in x
+                        int j1 = kernel1->yOrder; // Polynomial order in y
+                        // First convolution
+                        float conv1 = polyValues->data.F64[j1][i1] * reference[y - v1][x - u1];
+
+                        // Assuming that the first kernel component is 0 order in x and y, and 0 offset
+                        if (k1 != 0) {
+                            conv1 -= reference[y][x];
+                        }
+
+                        /* Iterate over the second convolution */
+                        for (int k2 = k1; k2 < numKernelParams; k2++) {
+                            poisKernelBasis *kernel2 = kernelParams->data[k2]; // Kernel basis #2
+
+                            int u2 = kernel2->u; // Offset in x
+                            int v2 = kernel2->v; // Offset in y
+                            int i2 = kernel2->xOrder; // Polynomial order in x
+                            int j2 = kernel2->yOrder; // Polynomial order in y
+                            // Second convolution
+                            float conv2 = polyValues->data.F64[j2][i2] * reference[y - v2][x - u2];
+
+                            // Assuming that the first kernel component is 0 order in x and y, and 0 offset
+                            if (k2 != 0) {
+                                conv2 -= reference[y][x];
+                            }
+
+                            // Add into the matrix element
+                            matrix->data.F64[k1][k2] += conv1 * // First convolution
+                                conv2 * // Second convolution
+                                invNoise2; // Divide by sigma^2, ~ poisson.
+                        } // Second convolution
+
+                        // Add into the vector element
+                        vector->data.F64[k1] += input[y][x] * // Input image, no convolution
+                            conv1 *     // Convolved reference image
+                            invNoise2; // Divide by sigma^2, ~ poisson.
+
+                        /* Background term */
+                        matrix->data.F64[k1][bgIndex] += conv1 * // Convolved reference image
+                            invNoise2; // Divide by sigma^2, ~ poisson.
+
+                    } // First convolution
+
+                    /* Background only terms */
+                    matrix->data.F64[bgIndex][bgIndex] += invNoise2;
+                    vector->data.F64[bgIndex] += input[y][x] * invNoise2;
+                }
+            } // Iterating over pixels in stamp
+
+            // Fill in other side of symmetric matrix
+            for (int k1 = 0; k1 < kernelParams->n; k1++) {
+                for (int k2 = 0; k2 < k1; k2++) {
+                    matrix->data.F64[k1][k2] = matrix->data.F64[k2][k1];
+                }
+                matrix->data.F64[bgIndex][k1] = matrix->data.F64[k1][bgIndex];
+            }
+
+            for (int k = 0; k < kernelParams->n; k++) {
+                poisKernelBasis *kernel = kernelParams->data[k]; // Kernel basis #2
+                int u = kernel->u; // Offset in x
+                int v = kernel->v; // Offset in y
+                matrix->data.F64[k][k] += config->penalty * (float)(u*u + v*v);
+            }
+
+            stamp->status = POIS_STAMP_USED; // We've calculated it now, so don't need to do so again.
+
+        } // Legitimate stamps
 
     } // Iterating over stamps
Index: /trunk/pois/src/poisConvolveImage.c
===================================================================
--- /trunk/pois/src/poisConvolveImage.c	(revision 5716)
+++ /trunk/pois/src/poisConvolveImage.c	(revision 5717)
@@ -7,91 +7,90 @@
 /* Convolve an image by the kernel */
 psImage *poisConvolveImage(const psImage *input, // Input image, to be convolved
-			   const psImage *mask,	// Mask image
-			   const psVector *solution, // The solution vector, containing the coefficients
-			   const psArray *kernelParams, // The kernel parameters
-			   const poisConfig *config // The configuration
+                           const psImage *mask, // Mask image
+                           const psVector *solution, // The solution vector, containing the coefficients
+                           const psArray *kernelParams, // The kernel parameters
+                           const poisConfig *config // The configuration
     )
 {
-    assert(solution->n == kernelParams->n + 1);	// Extra parameter for background
+    assert(solution->n == kernelParams->n + 1); // Extra parameter for background
     assert(solution->type.type == PS_TYPE_F64);
     assert(input->numCols == mask->numCols && input->numRows == mask->numRows);
     assert(mask->type.type == PS_TYPE_U8);
     assert(input->type.type == PS_TYPE_F32);
-    
-    int nx = input->numCols;		// The size of the image in x
-    int ny = input->numRows;		// The size of the image in y
-    int xKernel = config->xKernel;	// The half-size of the kernel in x
-    int yKernel = config->yKernel;	// The half-size of the kernel in y
-    float background = solution->data.F64[solution->n - 1]; // The difference in background
-    int numParams = kernelParams->n;	// Number of kernel parameters
+
+    int nx = input->numCols;            // The size of the image in x
+    int ny = input->numRows;            // The size of the image in y
+    float nxHalf = 0.5 * (float)nx;     // Half the size of the image
+    float nyHalf = 0.5 * (float)ny;     // Half the size of the image
+    int xKernel = config->xKernel;      // The half-size of the kernel in x
+    int yKernel = config->yKernel;      // The half-size of the kernel in y
+    float background = solution->data.F64[solution->n-1]; // The difference in background
+    int numParams = kernelParams->n;    // Number of kernel parameters
 
     psImage *convolved = psImageAlloc(nx, ny, PS_TYPE_F32); // The convolved image, to be returned
 
-    psDPolynomial2D *poly = psDPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
-						 PS_POLYNOMIAL_ORD); // Polynomial
+    psPolynomial2D *poly = psPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
+                                               PS_POLYNOMIAL_ORD ); // Polynomial
     for (int i = 0; i < config->spatialOrder + 1; i++) {
-	for (int j = 0; j < config->spatialOrder + 1; j++) {
-	    poly->coeff[i][j] = 1.0;
-	    poly->mask[i][j] = 1;	// Mask all coefficients; unmask to evaluate
-	}
+        for (int j = 0; j < config->spatialOrder + 1; j++) {
+            poly->coeff[i][j] = 1.0;
+            poly->mask[i][j] = 1;       // Mask all coefficients; unmask to evaluate
+        }
     }
 
     // Convolve only the middle part
     for (int y = yKernel; y < ny - yKernel; y++) {
-	for (int x = xKernel; x < nx - xKernel; x++) {
-	    if (! mask->data.U8[y][x] & (POIS_MASK_BAD | POIS_MASK_NEAR_BAD)) {
-		double conv = background; // Convolved value
+        for (int x = xKernel; x < nx - xKernel; x++) {
+            if (! mask->data.U8[y][x] & (POIS_MASK_BAD | POIS_MASK_NEAR_BAD)) {
+                double conv = background; // Convolved value
 
-#if 0
-		float imageX = (x - nxHalf) / nxHalf; // Normalised position in x
-		float imageY = (y - nyHalf) / nyHalf; // Normalised position in y
-#endif
-		
-		// Iterate over the kernel basis functions
-		for (int k = 0; k < numParams; k++) {
-		    poisKernelBasis *kernel = kernelParams->data[k]; // The kernel basis function
-		    int u = kernel->u;
-		    int v = kernel->v;
-		    
-#if 0
-		    int xOrder = kernel->xOrder;
-		    int yOrder = kernel->yOrder;
-		    // Evaluate the polynomial
-		    poly->mask[xOrder][yOrder] = 0;
-		    double polyVal = psDPolynomial2DEval(poly, imageX, imageY);
-		    poly->mask[xOrder][yOrder] = 1;
-#else
-		    double polyVal = 1.0;
-#endif
+                float imageX = (x - nxHalf) / nxHalf; // Normalised position in x
+                float imageY = (y - nyHalf) / nyHalf; // Normalised position in y
 
-		    if (k == 0) {
-			conv += solution->data.F64[k] * input->data.F32[y - v][x - u] * polyVal;
-		    } else {
-			conv += solution->data.F64[k] *
-			    (input->data.F32[y - v][x - u] * polyVal - input->data.F32[y][x]);
-		    }
-		}
-		convolved->data.F32[y][x] = (float)conv;
-	    }
-	}
-	    
-	/* Pad the rest with zeros */
-	for (int x = 0; x < xKernel; x++) {
-	    convolved->data.F32[y][x] = 0.0;
-	}
-	for (int x = nx - xKernel; x < nx; x++) {
-	    convolved->data.F32[y][x] = 0.0;
-	}
+                // Iterate over the kernel basis functions
+                for (int k = 0; k < numParams; k++) {
+                    poisKernelBasis *kernel = kernelParams->data[k]; // The kernel basis function
+                    int u = kernel->u;
+                    int v = kernel->v;
+                    int xOrder = kernel->xOrder;
+                    int yOrder = kernel->yOrder;
+
+                    double polyVal = 1.0;
+                    if (xOrder != 0 || yOrder != 0) {
+                        // Evaluate the polynomial
+                        poly->mask[xOrder][yOrder] = 0;
+                        polyVal = psPolynomial2DEval(poly, imageX, imageY);
+                        poly->mask[xOrder][yOrder] = 1;
+                    }
+
+                    if (k == 0) {
+                        conv += solution->data.F64[k] * input->data.F32[y - v][x - u] * polyVal;
+                    } else {
+                        conv += solution->data.F64[k] *
+                            (input->data.F32[y - v][x - u] * polyVal - input->data.F32[y][x]);
+                    }
+                }
+                convolved->data.F32[y][x] = (float)conv;
+            }
+        }
+
+        /* Pad the rest with zeros */
+        for (int x = 0; x < xKernel; x++) {
+            convolved->data.F32[y][x] = 0.0;
+        }
+        for (int x = nx - xKernel; x < nx; x++) {
+            convolved->data.F32[y][x] = 0.0;
+        }
     }
 
     for (int y = 0; y < yKernel; y++) {
-	for (int x = 0; x < nx; x++) {
-	    convolved->data.F32[y][x] = 0.0;
-	}
+        for (int x = 0; x < nx; x++) {
+            convolved->data.F32[y][x] = 0.0;
+        }
     }
     for (int y = ny - yKernel; y < ny; y++) {
-	for (int x = 0; x < nx; x++) {
-	    convolved->data.F32[y][x] = 0.0;
-	}
+        for (int x = 0; x < nx; x++) {
+            convolved->data.F32[y][x] = 0.0;
+        }
     }
 
Index: /trunk/pois/src/poisExtractKernel.c
===================================================================
--- /trunk/pois/src/poisExtractKernel.c	(revision 5716)
+++ /trunk/pois/src/poisExtractKernel.c	(revision 5717)
@@ -30,6 +30,6 @@
     }
 
-    psDPolynomial2D *poly = psDPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
-						 PS_POLYNOMIAL_ORD); // Polynomial for evaluation
+    psPolynomial2D *poly = psPolynomial2DAlloc(config->spatialOrder + 1, config->spatialOrder + 1,
+					       PS_POLYNOMIAL_ORD); // Polynomial for evaluation
     for (int j = 0; j < config->spatialOrder + 1; j++) {
 	for (int i = 0; i < config->spatialOrder + 1; i++) {
@@ -51,5 +51,5 @@
 	// Evaluate polynomial
 	poly->mask[xOrder][yOrder] = 0;
-	double evaluation = solution->data.F64[k] * psDPolynomial2DEval(poly, x, y);
+	double evaluation = solution->data.F64[k] * psPolynomial2DEval(poly, x, y);
 	poly->mask[xOrder][yOrder] = 1;
 
Index: /trunk/pois/src/poisFindStamps.c
===================================================================
--- /trunk/pois/src/poisFindStamps.c	(revision 5716)
+++ /trunk/pois/src/poisFindStamps.c	(revision 5717)
@@ -7,11 +7,11 @@
 
 psArray *poisFindStamps(psArray *stamps, // Existing list of stamps, or NULL
-			const psImage *image, // Image for which to find stamps
-			const psImage *mask, // Mask image
-			const poisConfig *config // Configuration values
+                        const psImage *image, // Image for which to find stamps
+                        const psImage *mask, // Mask image
+                        const poisConfig *config // Configuration values
     )
 {
-    int nx = config->xImage;		// Size of image in x
-    int ny = config->yImage;		// Size of image in y
+    int nx = config->xImage;            // Size of image in x
+    int ny = config->yImage;            // Size of image in y
 
     assert(image->numCols == nx && image->numRows == ny);
@@ -21,20 +21,18 @@
     assert(!stamps || stamps->n == config->nsx * config->nsy);
 
-    float threshold = config->threshold;// Stamp threshold
-    int footprint = config->footprint;	// Footprint for stamps
+    int footprint = config->footprint;  // Footprint for stamps
     int borderx = footprint + config->xKernel; // Border around image in x
     int bordery = footprint + config->yKernel; // Border around image in y
-    int nsx = config->nsx;		// Number of stamps in x
-    int nsy = config->nsy;		// Number of stamps in y
+    int nsx = config->nsx;              // Number of stamps in x
+    int nsy = config->nsy;              // Number of stamps in y
 
-    psF32 **pixels = image->data.F32;	// The image pixels
-    psU8 **maskpix = mask->data.U8; // The mask pixels
+    psF32 **pixels = image->data.F32;   // The image pixels
 
     if (stamps == NULL) {
-	stamps = psArrayAlloc(nsx * nsy);
-					// Initialise; should be done by psArrayAlloc
-	for (int i = 0; i < nsx * nsy; i++) {
-	    stamps->data[i] = NULL;
-	}
+        stamps = psArrayAlloc(nsx * nsy);
+        // Initialise
+        for (int i = 0; i < nsx * nsy; i++) {
+            stamps->data[i] = poisStampAlloc();
+        }
     }
 
@@ -42,57 +40,42 @@
     int num = 0;
     for (int i = 0; i < nsx; i++) {
-	for (int j = 0; j < nsy; j++) {
-	    poisStamp *stamp = stamps->data[num]; // The stamp
-	    if (stamp == NULL) {
-		stamp = stamps->data[num] = poisStampAlloc();
-	    }
+        for (int j = 0; j < nsy; j++) {
+            poisStamp *stamp = stamps->data[num]; // The stamp
 
-	    // Only find a new stamp if we need to
-	    if (stamp->status == POIS_STAMP_RESET) {
-		// Find maximum non-masked value in the image section, but don't include a footprint around
-		// the edge
-		float max = - INFINITY;		// Negative infinity
-		int bestx = 0, besty = 0;	// Position of maximum
-		
-		for (int y = bordery + j * (ny - 2.0 * bordery) / nsy; 
-		     y < bordery + (j + 1) * (ny - 2.0 * bordery) / nsy; y++) {
-		    for (int x = borderx + i * (nx - 2.0 * borderx) / nsx;
-			 x < borderx + (i + 1) * (nx - 2.0 * borderx) / nsx; x++) {
-			if (pixels[y][x] > max && pixels[y][x] >= threshold) {
-			    // Check the footprint
-			    int ok = 1;
-			    for (int v = y - footprint; v <= y + footprint && ok; v++) {
-				for (int u = x - footprint; u <= x + footprint && ok; u++) {
-				    if (maskpix[v][u] &
-					(POIS_MASK_BAD | POIS_MASK_NEAR_BAD | POIS_MASK_STAMP)) {
-					ok = 0;
-				    }
-				}
-			    }
-			    if (ok) {
-				max = pixels[y][x];
-				bestx = x;
-				besty = y;
-			    }
-			} // Done checking the candidate pixel
-		    }
-		} // Done iterating over this image section
-		
-		// Store the stamp
-		if (max != - INFINITY) {
-		    stamp->x = bestx;
-		    stamp->y = besty;
-		    stamp->status = POIS_STAMP_RECALC;
-		} else {
-		    // No stamp in this section
-		    stamp->status = POIS_STAMP_BAD;
-		}
-	    } // Finding a new stamp
+            // Only find a new stamp if we need to
+            if ((stamp->x == 0 && stamp->y == 0) || stamp->status == POIS_STAMP_RESET) {
+                // Find maximum non-masked value in the image section, but don't include a footprint around
+                // the edge
+                float max = - INFINITY;         // Negative infinity
+                int bestx = 0, besty = 0;       // Position of maximum
 
-	    psTrace("pois.findStamps", 7, "Stamp %d: %d,%d\n", num, stamp->x, stamp->y);
-	    num++;
-	}
+                for (int y = bordery + j * (ny - 2.0 * bordery) / nsy;
+                     y < bordery + (j + 1) * (ny - 2.0 * bordery) / nsy; y++) {
+                    for (int x = borderx + i * (nx - 2.0 * borderx) / nsx;
+                         x < borderx + (i + 1) * (nx - 2.0 * borderx) / nsx; x++) {
+                        if (pixels[y][x] > max && poisCheckStamp(image, mask, x, y, config)) {
+                            max = pixels[y][x];
+                            bestx = x;
+                            besty = y;
+                        }
+                    }
+                } // Done iterating over this image section
+
+                // Store the stamp
+                if (max != - INFINITY) {
+                    stamp->x = bestx;
+                    stamp->y = besty;
+                    stamp->status = POIS_STAMP_RECALC;
+                } else {
+                    // No stamp in this section
+                    stamp->status = POIS_STAMP_BAD;
+                }
+            } // Finding a new stamp
+
+            psTrace("pois.findStamps", 7, "Stamp %d: %d,%d\n", num, stamp->x, stamp->y);
+            num++;
+        }
     } // Done iterating over sections
 
     return stamps;
-}				    
+}
Index: /trunk/pois/src/poisKernelBasisFunctions.c
===================================================================
--- /trunk/pois/src/poisKernelBasisFunctions.c	(revision 5716)
+++ /trunk/pois/src/poisKernelBasisFunctions.c	(revision 5717)
@@ -36,5 +36,5 @@
 	    array->data[num] = kernel;
 
-	    psTrace("pois.kernelBasisFunctions", 8, "--> %d: (%d,%d) x^%d y^%d\n", i, kernel->u,
+	    psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", i, kernel->u,
 		    kernel->v, kernel->xOrder, kernel->yOrder);
 	    num++;
@@ -64,5 +64,5 @@
 			array->data[num] = kernel;
 
-			psTrace("pois.kernelBasisFunctions", 8, "--> %d: (%d,%d) x^%d y^%d\n", num,
+			psTrace("pois.kernelBasisFunctions", 3, "--> %d: (%d,%d) x^%d y^%d\n", num,
 				kernel->u, kernel->v, kernel->xOrder, kernel->yOrder);
 			num++;
Index: /trunk/pois/src/poisMakeMask.c
===================================================================
--- /trunk/pois/src/poisMakeMask.c	(revision 5716)
+++ /trunk/pois/src/poisMakeMask.c	(revision 5717)
@@ -41,5 +41,5 @@
 	    // Mask all around something that's saturated in the reference, since it gets convolved
 	    // But only mask single pixel in the input image, since it doesn't get convolved
-	    if (refImage->data.F32[y][x] >= config->refSat || refImage->data.F32[y][x] <= config->bad ||
+	    if (refImage->data.F32[y][x] >= config->refSat || refImage->data.F32[y][x] <= config->refBad ||
 		isnan(refImage->data.F32[y][x])) {
 		mask->data.U8[y][x] |= POIS_MASK_BAD;
@@ -53,5 +53,5 @@
 		    }
 		}
-	    } else if (inImage->data.F32[y][x] >= config->inSat || inImage->data.F32[y][x] <= config->bad ||
+	    } else if (inImage->data.F32[y][x] >= config->inSat || inImage->data.F32[y][x] <= config->inBad ||
 		isnan(inImage->data.F32[y][x])) {
 		mask->data.U8[y][x] |= POIS_MASK_BAD;
Index: /trunk/pois/src/poisParseConfig.c
===================================================================
--- /trunk/pois/src/poisParseConfig.c	(revision 5716)
+++ /trunk/pois/src/poisParseConfig.c	(revision 5717)
@@ -4,8 +4,8 @@
 #include "pois.h"
 
-static void help(void)
+void help(void)
 {
     fprintf (stderr, "POIS: Pan-STARRS (or Pricey's) Optimal Image Subtraction\n"
-	     "Usage: pois [-h] [-v] [-k X Y] [-t T] [-s S1 S2] [-b BAD] [-B BG] [-o N] [-n NSTAMPS] [-i ITER] [-r REJ] REF IN OUT\n"
+	     "Usage: pois [-h] [-v] [-k X Y] [-t T] [-s S1 S2] [-b B1 B2] [-B BG] [-o N] [-n NSTAMPS] [-S STAMPS] [-i ITER] [-r REJ] [-m] REF IN OUT\n"
 	     "where\n"
 	     "\t-h         Help (this info)\n"
@@ -14,10 +14,13 @@
 	     "\t-t T       Threshold for stamps on reference image (0)\n"
 	     "\t-s S1 S2   Saturation level for reference (S1) and input (S2) (50000)\n"
-	     "\t-b BAD     Bad level for both images (0)\n"
+	     "\t-b B1 B2   Bad level for both images (0)\n"
 	     "\t-B BG      Level to add to background to get above zero (0)\n"
 	     "\t-o ORDER   Order of spatial variations in the kernel (0)\n"
 	     "\t-n NSTAMPS Number of stamps in each dimension (5)\n"
+	     "\t-S STAMPS  File from which to read stamps\n"
+	     "\t-R         Reverse the sense of the subtraction (REF - IN)\n"
 	     "\t-i ITER    Number of rejection iterations (10)\n"
 	     "\t-r REJ     Rejection level in standard deviations (2.5)\n"
+	     "\t-m         Output mask frame (OUT.mask)\n"
 	     "\tREF        Reference image\n"
 	     "\tIN         Input image (to be matched to REF)\n"
@@ -26,7 +29,17 @@
 }
 
-poisConfig *restrict poisConfigAlloc(void)
+poisConfig *poisParseConfig(int argc, // Number of command-line arguments
+			    char **argv // Command-line arguments
+    )
 {
-    poisConfig *config = (poisConfig *) psAlloc(sizeof(poisConfig));
+    poisConfig *config;			// Configuration values
+
+    /* Variables for getopt */
+    int opt;   /* Option, from getopt */
+    extern char *optarg;   /* Argument accompanying switch */
+    extern int optind;   /* getopt variables */
+
+    /* Initialise configuration */
+    config = (poisConfig *) psAlloc(sizeof(poisConfig));
     config->verbose = 0;
     config->refFile = NULL;
@@ -39,5 +52,6 @@
     config->refSat = 60000.0;
     config->inSat = 60000.0;
-    config->bad = 0.0;
+    config->refBad = 0.0;
+    config->inBad = 0.0;
     config->background = 0.0;
     config->spatialOrder = 0;
@@ -47,24 +61,10 @@
     config->sigmaRej = 2.5;
     config->penalty = 0.0;
-
-    return config;
-}
-
-
-poisConfig *poisParseConfig(int argc, // Number of command-line arguments
-			    char **argv // Command-line arguments
-    )
-{
-    poisConfig *config  = poisConfigAlloc(); // Configuration values
-
-    /* Variables for getopt */
-    int opt;   /* Option, from getopt */
-    extern char *optarg;   /* Argument accompanying switch */
-    extern int optind;   /* getopt variables */
-
-    /* Initialise configuration */
+    config->stampFile = NULL;
+    config->reverse = false;
+    config->mask = false;
 
     /* Parse command-line arguments using getopt */
-    while ((opt = getopt(argc, argv, "hvk:f:q:t:s:o:b:B:n:i:r:p:")) != -1) {
+    while ((opt = getopt(argc, argv, "hvRmk:f:q:t:s:o:b:B:n:i:r:p:S:")) != -1) {
         switch (opt) {
 	  case 'h':
@@ -115,8 +115,13 @@
 	    break;
 	  case 'b':
-	    if (sscanf(optarg, "%f", &config->bad) != 1) {
-		help();
+            if ((sscanf(argv[optind-1], "%f", &config->refBad) != 1) ||
+                (sscanf(argv[optind++], "%f", &config->inBad) != 1)) {
+                /*
+                  Note: incrementing optind, so I can read more than
+                  one parameter.
+                */
+                help();
                 exit(EXIT_FAILURE);
-	    }
+            }
 	    break;
 	  case 'B':
@@ -151,4 +156,14 @@
             }
 	    break;
+	  case 'S':
+	    config->stampFile = argv[optind-1];
+	    config->threshold = -INFINITY;
+	    break;
+	  case 'R':
+	    config->reverse = true;
+	    break;
+	  case 'm':
+            config->mask = true;
+            break;
 	  default:
 	    help();
Index: /trunk/pois/src/poisReadStamps.c
===================================================================
--- /trunk/pois/src/poisReadStamps.c	(revision 5716)
+++ /trunk/pois/src/poisReadStamps.c	(revision 5717)
@@ -3,59 +3,56 @@
 #include "pois.h"
 
-#define MAXLINE 128			// Maximum length of a line in the stamps file
+#define MAXLINE 128                     // Maximum length of a line in the stamps file
 
 
 psArray *poisReadStamps(psArray *stamps,// Stamps
-			FILE **stampFP, // File pointer for stamps file
-			const psImage *image, // Image for which to find stamps
-			const psImage *mask, // Mask image
-			const poisConfig *config // Configuration
+                        FILE **stampFP, // File pointer for stamps file
+                        const psImage *image, // Image for which to find stamps
+                        const psImage *mask, // Mask image
+                        const poisConfig *config // Configuration
     )
 {
-    FILE *fp = *stampFP;		// File pointer
-
-    int num = 0;			// Number of stamps
+    FILE *fp = *stampFP;                // File pointer
 
     if (stamps == NULL) {
-	stamps = psArrayAlloc(config->nsx);
-	// Initialise
-	for (int i = 0; i < config->nsx; i++) {
-	    stamps->data[i] = poisStampAlloc();
-	}
+        stamps = psArrayAlloc(config->nsx);
+        // Initialise
+        for (int i = 0; i < config->nsx; i++) {
+            stamps->data[i] = poisStampAlloc();
+        }
     }
 
     for (int i = 0; i < stamps->n; i++) {
-	poisStamp *stamp = stamps->data[i];
+        poisStamp *stamp = stamps->data[i];
 
-	if (stamp->status == POIS_STAMP_RESET) {
-	    if (fp) {
-		float x, y;		// Position of stamp
-		int xPixel, yPixel;	// Integer position
-		int result = 0;		// Result of fscanf
-		bool badPixels = false;	// Result of check on stamp footprint
-		char line[MAXLINE];	// The line from the stamps file
-		char *lineBack;		// The line back from fgets
-		do {
-		    lineBack = fgets(line, MAXLINE, fp);
-		    result = sscanf(line, "%f %f", &x, &y);
+        if (stamp->status == POIS_STAMP_RESET) {
+            if (fp) {
+                float x, y;             // Position of stamp
+                int xPixel, yPixel;     // Integer position
+                int result = 0;         // Result of fscanf
+                char line[MAXLINE];     // The line from the stamps file
+                char *lineBack;         // The line back from fgets
+                do {
+                    lineBack = fgets(line, MAXLINE, fp);
+                    result = sscanf(line, "%f %f", &x, &y);
 
-		    // Add 0.5 for rounding, and subtract 1.0 for unit-offset to zero-offset
-		    xPixel = (int)(x - 0.5);
-		    yPixel = (int)(y - 0.5);
-		} while (result == 2 && lineBack && ! poisCheckStamp(image, mask, xPixel, yPixel, config));
-		if (result == 2 && lineBack) {
-		    stamp->x = xPixel;
-		    stamp->y = yPixel;
-		    stamp->status = POIS_STAMP_RECALC;
-		    psTrace("pois.readStamps", 7, "Stamp %d: %d,%d\n", i, stamp->x, stamp->y);
-		} else {
-		    // We've run out of stamps in the file --- probably got EOF
-		    stamp->status = POIS_STAMP_BAD;
-		}
-	    } else {
-		// Stamp file isn't open!
-		stamp->status = POIS_STAMP_BAD;
-	    }
-	}
+                    // Add 0.5 for rounding, and subtract 1.0 for unit-offset to zero-offset
+                    xPixel = (int)(x - 0.5);
+                    yPixel = (int)(y - 0.5);
+                } while (result == 2 && lineBack && ! poisCheckStamp(image, mask, xPixel, yPixel, config));
+                if (result == 2 && lineBack) {
+                    stamp->x = xPixel;
+                    stamp->y = yPixel;
+                    stamp->status = POIS_STAMP_RECALC;
+                    psTrace("pois.readStamps", 7, "Stamp %d: %d,%d\n", i, stamp->x, stamp->y);
+                } else {
+                    // We've run out of stamps in the file --- probably got EOF
+                    stamp->status = POIS_STAMP_BAD;
+                }
+            } else {
+                // Stamp file isn't open!
+                stamp->status = POIS_STAMP_BAD;
+            }
+        }
     }
 
Index: /trunk/pois/src/poisRejectStamps.c
===================================================================
--- /trunk/pois/src/poisRejectStamps.c	(revision 5716)
+++ /trunk/pois/src/poisRejectStamps.c	(revision 5717)
@@ -6,8 +6,8 @@
 #define SQUARE(x) ((x)*(x))
 
-psList *restrict poisRejectStamps(psArray *stamps, // Array of stamps
-				  psImage *mask,	// Mask image
-				  const psVector *deviations, // Vector of deviations for the stamps
-				  const poisConfig *config // Configuration
+bool poisRejectStamps(psArray *stamps,  // Array of stamps
+                      psImage *mask,    // Mask image
+                      const psVector *deviations, // Vector of deviations for the stamps
+                      const poisConfig *config // Configuration
     )
 {
@@ -21,56 +21,46 @@
     assert(deviations->type.type == PS_TYPE_F32);
 
-    psVector *devMask = psVectorAlloc(stamps->n, PS_TYPE_U8); // Mask for statistics
+    bool masked = false;                // Have we masked any stamps?
 
     // Don't want the standard deviation, but rather the deviation from zero
-    double meanDev = 0.0;
+    double sum = 0.0;
     int numDev = 0;
     for (int i = 0; i < deviations->n; i++) {
-	poisStamp *stamp = stamps->data[i];	// The stamp
-	// Only interested in the stamps that we're actually using
-	if (stamp != NULL && stamp->status == POIS_STAMP_USED) {
-	    meanDev += SQUARE(deviations->data.F32[i]);
-	    numDev++;
-	}
+        poisStamp *stamp = stamps->data[i];     // The stamp
+        // Only interested in the stamps that we're actually using
+        if (stamp->status == POIS_STAMP_USED) {
+            sum += SQUARE(deviations->data.F32[i]);
+            numDev++;
+        }
     }
-    float rmsDev = sqrt(meanDev/numDev);
-    float limit = rmsDev * config->sigmaRej;
-    psTrace("pois.rejectStamps", 2, "RMS deviation: %f\n", rmsDev);
+    float meanDev = sqrtf(sum / (float)numDev);
+    float limit = meanDev * config->sigmaRej;
+    psTrace("pois.rejectStamps", 2, "Mean RMS deviation: %f\n", meanDev);
     psTrace("pois.rejectStamps", 2, "Rejection limit: %f\n", limit);
 
     // Reject stamps
-    psList *rejected = NULL;		// rejected stamps
     for (int s = 0; s < deviations->n; s++) {
-	poisStamp *stamp = stamps->data[s];
-	if (stamp == NULL) {
-	    continue;
-	}
-	
-	if (stamp->status == POIS_STAMP_USED && fabsf(deviations->data.F32[s]) > limit) {
-	    psTrace("pois.rejectStamps", 1, "Rejecting stamp %d (%d,%d): %f\n", s, stamp->x, stamp->y,
-		    deviations->data.F32[s]);
-	    
-	    // Mask out the stamp in the image so you don't find it again
-	    for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) {
-		for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) {
-		    mask->data.U8[y][x] |= POIS_MASK_STAMP;
-		}
-	    }
-	    
-	    // Mark stamp for replacement
-	    stamps->data[s] = NULL;
-	    
-	    if (rejected == NULL) {
-		rejected = psListAlloc(stamp);
-	    } else {
-		bool success = psListAdd(rejected, PS_LIST_TAIL, stamp);
-		assert(success);
-	    }
-	    psFree(stamp);
-	} // Bad stamps
+        poisStamp *stamp = stamps->data[s];
+        if (stamp->status == POIS_STAMP_USED &&
+            (fabsf(deviations->data.F32[s]) > limit || ! isfinite(deviations->data.F32[s]))) {
+            masked = true;
+            psTrace("pois.rejectStamps", 1, "Rejecting stamp %d (%d,%d): %f\n", s, stamp->x, stamp->y,
+                    deviations->data.F32[s]);
+
+            // Mask out the stamp in the image so you don't find it again
+            for (int y = stamp->y - config->footprint; y < stamp->y + config->footprint; y++) {
+                for (int x = stamp->x - config->footprint; x < stamp->x + config->footprint; x++) {
+                    mask->data.U8[y][x] |= POIS_MASK_STAMP;
+                }
+            }
+
+            // Set stamp for replacement
+            stamp->x = 0;
+            stamp->y = 0;
+            stamp->status = POIS_STAMP_RESET;
+
+        } // Bad stamps
     } // Iterating over stamps
 
-    psFree(devMask);
-
-    return rejected;
+    return masked;
 }
Index: /trunk/pois/src/poisSolveEquation.c
===================================================================
--- /trunk/pois/src/poisSolveEquation.c	(revision 5716)
+++ /trunk/pois/src/poisSolveEquation.c	(revision 5717)
@@ -39,15 +39,4 @@
 
 	if (stamp->status == POIS_STAMP_USED) {
-
-#ifdef TESTING
-	    printf("Matrix %d:\n", s);
-	    for (int i = 0; i < matrix->numCols; i++) {
-		for (int j = 0; j < matrix->numRows; j++) {
-		    printf("%f ", matrix->data.F64[i][j]);
-		}
-		printf("\n");
-	    }
-#endif
-
 	    (void)psBinaryOp(matrix, matrix, "+", stampMatrix);
 	    (void)psBinaryOp(vector, vector, "+", stampVector);
@@ -61,4 +50,14 @@
     psTrace("pois.solveEquation", 1, "Background difference is %f\n", solution->data.F64[solution->n - 1]);
 
+#if 0
+    printf("Matrix:\n");
+    for (int i = 0; i < matrix->numCols; i++) {
+	for (int j = 0; j < matrix->numRows; j++) {
+	    printf("%f ", matrix->data.F64[i][j]);
+	}
+	printf("\n");
+    }
+#endif
+
     // Clean up
     psFree(luMatrix);
Index: /trunk/pois/src/poisStamp.c
===================================================================
--- /trunk/pois/src/poisStamp.c	(revision 5716)
+++ /trunk/pois/src/poisStamp.c	(revision 5717)
@@ -11,5 +11,5 @@
     stamp->matrix = NULL;
     stamp->vector = NULL;
-    psMemSetDeallocator(stamp, (psFreeFcn)poisStampFree);
+    psMemSetDeallocator(stamp, (psFreeFunc)poisStampFree);
 
     return stamp;
@@ -18,5 +18,44 @@
 void poisStampFree(poisStamp *stamp)
 {
-    psFree(stamp->matrix);
-    psFree(stamp->vector);
+    if (stamp->matrix) {
+	psFree(stamp->matrix);
+    }
+    if (stamp->vector) {
+	psFree(stamp->vector);
+    }
+    psFree(stamp);
 }
+
+
+// Return true if the stamp is OK
+bool poisCheckStamp(const psImage *image, // Image to check for threshold
+		    const psImage *mask, // Mask to check for bad pixels
+		    int x, int y,	// Pixel coordinates
+		    const poisConfig *config // Configuration
+    )
+{
+    if (x < config->footprint + config->xKernel ||
+	y < config->footprint + config->yKernel ||
+	x > config->xImage - config->footprint - config->xKernel ||
+	y > config->yImage - config->footprint - config->yKernel ||
+	image->data.F32[y][x] < config->threshold) {
+	psTrace("pois.checkStamp", 9, "Rejecting stamp at %d,%d (border/threshold)\n", x, y);
+	return false;
+    }
+
+    // Check the footprint
+    bool ok = true;			// Is the footprint OK?
+    int footprint = config->footprint;	// Footprint size
+    for (int v = y - footprint; v <= y + footprint && ok; v++) {
+	for (int u = x - footprint; u <= x + footprint && ok; u++) {
+	    if (mask->data.U8[v][u] &
+		(POIS_MASK_BAD | POIS_MASK_NEAR_BAD | POIS_MASK_STAMP)) {
+		psTrace("pois.checkStamp", 9, "Rejecting stamp at %d,%d (bad footprint)\n", x, y);
+		ok = false;
+	    }
+	}
+    }
+    
+    return ok;
+}
+