Index: trunk/pois/src/pois.c
===================================================================
--- trunk/pois/src/pois.c	(revision 3813)
+++ 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);
 }
+