Index: /trunk/pois/src/pois.c
===================================================================
--- /trunk/pois/src/pois.c	(revision 5741)
+++ /trunk/pois/src/pois.c	(revision 5742)
@@ -18,5 +18,5 @@
 {
     struct timeval tv;
- 
+
     gettimeofday(&tv, NULL);
     return(tv.tv_sec + 1.e-6 * tv.tv_usec);
@@ -35,73 +35,73 @@
 
     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);
+        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
-    psMetadata *header = NULL;		// Header for output image
+    psMetadata *header = NULL;          // Header for output image
 
     psRegion imageRegion = {0, 0, 0, 0};
-    psFits *reference = psFitsAlloc(config->refFile);
+    psFits *reference = psFitsOpen(config->refFile, "r");
     if (config->reverse) {
-	header = psFitsReadHeader(NULL, reference);
+        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);
-	exit(EXIT_FAILURE);
+        psErrorStackPrint(stderr, "Fatal error: unable to read %s\n", config->refFile);
+        exit(EXIT_FAILURE);
     }
     psTrace("pois", 3, "Read %s: %dx%d\n", config->refFile, refImage->numCols, refImage->numRows);
-    psFree(reference);
+    psFitsClose(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);
+        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 = psFitsOpen(config->inFile, "r");
     if (! config->reverse) {
-	header = psFitsReadHeader(NULL, input);
+        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);
-	exit(EXIT_FAILURE);
+        psErrorStackPrint(stderr, "Fatal error: unable to read %s\n", config->inFile);
+        exit(EXIT_FAILURE);
     }
     psTrace("pois", 3, "Read %s: %dx%d\n", config->inFile, inImage->numCols, inImage->numRows);
-    psFree(input);
+    psFitsClose(input);
     // 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;
+        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)) {
-	psError(PS_ERR_BAD_PARAMETER_SIZE, true,
-		"Reference and input images must have the same dimensions: %dx%d vs %dx%d\n",
-		refImage->numCols, refImage->numRows, inImage->numCols, inImage->numRows);
-	exit(EXIT_FAILURE);
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                "Reference and input images must have the same dimensions: %dx%d vs %dx%d\n",
+                refImage->numCols, refImage->numRows, inImage->numCols, inImage->numRows);
+        exit(EXIT_FAILURE);
     }
     config->xImage = refImage->numCols;
@@ -115,6 +115,6 @@
     // as the square root of the value.
     if (config->background != 0.0) {
-	(void)psBinaryOp(refImage, refImage, "+", psScalarAlloc(config->background, PS_TYPE_F32));
-	(void)psBinaryOp(inImage, inImage, "+", psScalarAlloc(config->background, PS_TYPE_F32));
+        (void)psBinaryOp(refImage, refImage, "+", psScalarAlloc(config->background, PS_TYPE_F32));
+        (void)psBinaryOp(inImage, inImage, "+", psScalarAlloc(config->background, PS_TYPE_F32));
     }
 
@@ -127,105 +127,105 @@
 #ifdef TESTING
     // Write the mask out
-    psFits *maskFile = psFitsAlloc("mask.fits");
+    psFits *maskFile = psFitsOpen("mask.fits", "w");
     if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
-	psErrorStackPrint(stderr, "Unable to write mask: mask.fits\n");
+        psErrorStackPrint(stderr, "Unable to write mask: mask.fits\n");
     }
     psTrace("pois", 1, "Mask written to mask.fits\n");
-    psFree(maskFile);
-#endif
-
-    psArray *stamps = NULL;		// Array of stamps
-    psVector *solution = NULL;		// Solution vector
-
-    FILE *stampsFP = NULL;		// File pointer for stamps
+    psFitsClose(maskFile);
+#endif
+
+    psArray *stamps = NULL;             // Array of 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);
-	}
+        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
-    bool badStamps = true;		// 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);
-
-	// Find stamps
-	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++) {
-	    poisStamp *stamp = stamps->data[s];	// Stamp of interest
-	    if (stamp->status != POIS_STAMP_BAD) {
-		numStamps++;
-	    }
-	}
-	psTrace("pois.time", 1, "%d stamps found at %f sec\n", numStamps, getTime() - startTime);
-	if (numStamps == 0) {
-	    fprintf(stderr, "No stamps found.  Check input parameters.\n");
-	    exit(EXIT_FAILURE);
-	}
-
-	// Calculate equation
-	(void)poisCalculateEquation(stamps, refImage, inImage, kernelBasisFunctions, config);
-	psTrace("pois.time", 1, "Equation calculated at %f sec\n", getTime() - startTime);
-
-	// Solve the equation
-	solution = poisSolveEquation(solution, stamps, config);
-	psTrace("pois.time", 1, "Matrix equation solved at %f sec\n", getTime() - startTime);
-
-#ifdef TESTING
-	// Print solution
-	psTrace("pois.solution", 5, "Solution:\n");
-	for (int i = 0; i < solution->n - 1; i++) {
-	    poisKernelBasis *kernel = kernelBasisFunctions->data[i]; // The i-th kernel basis function
-	    psTrace("pois.solution", 7, "%d: (%d,%d) x^%d y^%d --> %f\n", i, kernel->u, kernel->v,
-		    kernel->xOrder, kernel->yOrder, solution->data.F64[i]);
-	}
-#endif
-
-#ifdef TESTING
-	// Save the kernel postage stamp
-	char kernelName[MAXCHAR];		// File name for kernel
-	snprintf(kernelName, MAXCHAR, "kernel%d.fits", iterNum);
-	psImage *kernelImage = poisExtractKernel(solution, kernelBasisFunctions, 0.0, 0.0, config);
-	psFits *kernelFile = psFitsAlloc(kernelName);
-	if (!psFitsWriteImage(kernelFile, NULL, kernelImage, 0)) {
-	    psErrorStackPrint(stderr, "Unable to write kernel: %s\n", kernelName);
-	}
-	psTrace("pois", 1, "Kernel written to %s\n", kernelName);
-	psFree(kernelFile);
-	psFree(kernelImage);
-#endif
-
-	// Calculate deviations for the stamps
-	psVector *deviations = poisCalculateDeviations(NULL, stamps, refImage, inImage, mask,
-						       kernelBasisFunctions, solution, config);
-	badStamps = poisRejectStamps(stamps, mask, deviations, config);
-
-	psFree(deviations);
+        psTrace("pois", 1, "Iteration %d...\n", iterNum);
+
+        // Find stamps
+        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++) {
+            poisStamp *stamp = stamps->data[s]; // Stamp of interest
+            if (stamp->status != POIS_STAMP_BAD) {
+                numStamps++;
+            }
+        }
+        psTrace("pois.time", 1, "%d stamps found at %f sec\n", numStamps, getTime() - startTime);
+        if (numStamps == 0) {
+            fprintf(stderr, "No stamps found.  Check input parameters.\n");
+            exit(EXIT_FAILURE);
+        }
+
+        // Calculate equation
+        (void)poisCalculateEquation(stamps, refImage, inImage, kernelBasisFunctions, config);
+        psTrace("pois.time", 1, "Equation calculated at %f sec\n", getTime() - startTime);
+
+        // Solve the equation
+        solution = poisSolveEquation(solution, stamps, config);
+        psTrace("pois.time", 1, "Matrix equation solved at %f sec\n", getTime() - startTime);
+
+#ifdef TESTING
+        // Print solution
+        psTrace("pois.solution", 5, "Solution:\n");
+        for (int i = 0; i < solution->n - 1; i++) {
+            poisKernelBasis *kernel = kernelBasisFunctions->data[i]; // The i-th kernel basis function
+            psTrace("pois.solution", 7, "%d: (%d,%d) x^%d y^%d --> %f\n", i, kernel->u, kernel->v,
+                    kernel->xOrder, kernel->yOrder, solution->data.F64[i]);
+        }
+#endif
+
+#ifdef TESTING
+        // Save the kernel postage stamp
+        char kernelName[MAXCHAR];               // File name for kernel
+        snprintf(kernelName, MAXCHAR, "kernel%d.fits", iterNum);
+        psImage *kernelImage = poisExtractKernel(solution, kernelBasisFunctions, 0.0, 0.0, config);
+        psFits *kernelFile = psFitsOpen(kernelName, "w");
+        if (!psFitsWriteImage(kernelFile, NULL, kernelImage, 0)) {
+            psErrorStackPrint(stderr, "Unable to write kernel: %s\n", kernelName);
+        }
+        psTrace("pois", 1, "Kernel written to %s\n", kernelName);
+        psFitsClose(kernelFile);
+        psFree(kernelImage);
+#endif
+
+        // Calculate deviations for the stamps
+        psVector *deviations = poisCalculateDeviations(NULL, stamps, refImage, inImage, mask,
+                                                       kernelBasisFunctions, solution, config);
+        badStamps = poisRejectStamps(stamps, mask, deviations, config);
+
+        psFree(deviations);
     } // Iterating for good solution
 
     // If there was rejection on the last round, re-solve the equation using only the good stamps
     if (badStamps) {
-	solution = poisSolveEquation(solution, stamps, config);
+        solution = poisSolveEquation(solution, stamps, config);
     }
 
     if (stampsFP) {
-	fclose(stampsFP);
+        fclose(stampsFP);
     }
 
     // Undo the mask for bad stamp area
     for (int j = 0; j < mask->numRows; j++) {
-	for (int i = 0; i < mask->numCols; i++) {
-	    if (mask->data.U8[j][i] & POIS_MASK_STAMP) {
-		mask->data.U8[j][i] &= ~POIS_MASK_STAMP;
-	    }
-	}
+        for (int i = 0; i < mask->numCols; i++) {
+            if (mask->data.U8[j][i] & POIS_MASK_STAMP) {
+                mask->data.U8[j][i] &= ~POIS_MASK_STAMP;
+            }
+        }
     }
 
@@ -239,11 +239,11 @@
 #ifdef TESTING
     // Save the convolved image
-    char convName[MAXCHAR];		// Filename for convolved image
+    char convName[MAXCHAR];             // Filename for convolved image
     snprintf(convName, MAXCHAR, "%s.conv", config->outFile);
-    psFits *convFile = psFitsAlloc(convName);
+    psFits *convFile = psFitsOpen(convName, "w");
     if (!psFitsWriteImage(convFile, NULL, convImage, 0)) {
-	psErrorStackPrint(stderr, "Unable to write convolved image: %s\n", convName);
-    }
-    psFree(convFile);
+        psErrorStackPrint(stderr, "Unable to write convolved image: %s\n", convName);
+    }
+    psFitsClose(convFile);
     psTrace("pois", 1, "Convolved image written to %s\n", convName);
 #endif
@@ -252,7 +252,7 @@
     psImage *subImage = psImageAlloc(config->xImage, config->yImage, PS_TYPE_F32);
     if (config->reverse) {
-	(void)psBinaryOp(subImage, convImage, "-", inImage);
+        (void)psBinaryOp(subImage, convImage, "-", inImage);
     } else {
-	(void)psBinaryOp(subImage, inImage, "-", convImage);
+        (void)psBinaryOp(subImage, inImage, "-", convImage);
     }
     psTrace("pois.time", 1, "Subtraction completed at %f sec\n", getTime() - startTime);
@@ -260,8 +260,8 @@
 #ifdef TESTING
     {
-	psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
-	(void)psImageStats(stats, subImage, mask, 0);
-	psTrace("pois", 5, "Subtracted image statistics: %f %f\n", stats->sampleMean, stats->sampleStdev);
-	psFree(stats);
+        psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+        (void)psImageStats(stats, subImage, mask, 0);
+        psTrace("pois", 5, "Subtracted image statistics: %f %f\n", stats->sampleMean, stats->sampleStdev);
+        psFree(stats);
     }
 #endif
@@ -269,23 +269,23 @@
     // Mask out the edges
     for (int y = config->yKernel; y < subImage->numRows - config->yKernel; y++) {
-	for (int x = config->xKernel; x < subImage->numCols - config->xKernel; x++) {
-
-	    for (int x = 0; x < config->xKernel; x++) {
-		subImage->data.F32[y][x] = 0.0;
-	    }
-	    for (int x = subImage->numCols - config->xKernel; x < subImage->numCols; x++) {
-		subImage->data.F32[y][x] = 0.0;
-	    }
-	}
+        for (int x = config->xKernel; x < subImage->numCols - config->xKernel; x++) {
+
+            for (int x = 0; x < config->xKernel; x++) {
+                subImage->data.F32[y][x] = 0.0;
+            }
+            for (int x = subImage->numCols - config->xKernel; 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 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 x = 0; x < subImage->numCols; x++) {
+            subImage->data.F32[y][x] = 0.0;
+        }
     }
 
@@ -293,20 +293,20 @@
     // Mask out bad pixels
     for (int y = config->yKernel; y < subImage->numRows - config->yKernel; y++) {
-	for (int x = config->xKernel; x < subImage->numCols - config->xKernel; x++) {
-	    if (mask->data.U8[y][x]) {
-		subImage->data.F32[y][x] = 0.0;
-	    }
-	}
+        for (int x = config->xKernel; x < subImage->numCols - config->xKernel; x++) {
+            if (mask->data.U8[y][x]) {
+                subImage->data.F32[y][x] = 0.0;
+            }
+        }
     }
 
     if (config->mask) {
-	char maskName[80];		// Name of mask image
-	sprintf(maskName, "%s.mask", config->outFile);
-	psFits *maskFile = psFitsAlloc(maskName);
+        char maskName[80];              // Name of mask image
+        sprintf(maskName, "%s.mask", config->outFile);
+        psFits *maskFile = psFitsOpen(maskName, "w");
         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);
+        psFitsClose(maskFile);
     }
 
@@ -314,12 +314,12 @@
     int numNaN = psImageClipNaN(subImage, 0.0);
     if (numNaN > 0) {
-	printf("Masked %d NAN pixels to zero.\n", numNaN);
-    }
-
-    psFits *subFile = psFitsAlloc(config->outFile);
+        printf("Masked %d NAN pixels to zero.\n", numNaN);
+    }
+
+    psFits *subFile = psFitsOpen(config->outFile, "w");
     if (!psFitsWriteImage(subFile, header, subImage, 0)) {
-	psErrorStackPrint(stderr, "Unable to write subtracted image: %s\n", config->outFile);
-    }
-    psFree(subFile);
+        psErrorStackPrint(stderr, "Unable to write subtracted image: %s\n", config->outFile);
+    }
+    psFitsClose(subFile);
     psTrace("pois", 1, "Subtracted image written to %s\n", config->outFile);
     psTrace("pois.time", 1, "I/O completed at %f sec\n", getTime() - startTime);
Index: /trunk/pois/src/poisCalculateDeviations.c
===================================================================
--- /trunk/pois/src/poisCalculateDeviations.c	(revision 5741)
+++ /trunk/pois/src/poisCalculateDeviations.c	(revision 5742)
@@ -6,12 +6,12 @@
 #define MAXCHAR 80
 
-psVector *poisCalculateDeviations(psVector *deviations,	// Output array of deviations, or NULL
-				  psArray *stamps, // Array of stamps
-				  psImage *refImage, // Reference image
-				  psImage *inImage, // Input image
-				  psImage *mask, // Mask image
-				  psArray *kernelParams, // Array of kernel parameters
-				  psVector *solution, // Solution vector
-				  poisConfig *config // Configuration
+psVector *poisCalculateDeviations(psVector *deviations, // Output array of deviations, or NULL
+                                  psArray *stamps, // Array of stamps
+                                  psImage *refImage, // Reference image
+                                  psImage *inImage, // Input image
+                                  psImage *mask, // Mask image
+                                  psArray *kernelParams, // Array of kernel parameters
+                                  psVector *solution, // Solution vector
+                                  poisConfig *config // Configuration
     )
 {
@@ -32,63 +32,63 @@
 
     if (!deviations) {
-	deviations = psVectorAlloc(stamps->n, PS_TYPE_F32);
+        deviations = psVectorAlloc(stamps->n, PS_TYPE_F32);
     }
 
-    int footprint = config->footprint;	// Size of stamp footprint
+    int footprint = config->footprint;  // Size of stamp footprint
     int xSize = footprint + config->xKernel; // (Half) size of subimage in x
     int ySize = footprint + config->yKernel; // (Half) size of subimage in y
-    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);	// Statistics
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics
 
     psImage *subStamp = psImageAlloc(2 * xSize, 2 * ySize, PS_TYPE_F32); // Subtraction of stamp
     for (int s = 0; s < stamps->n; s++) {
-	poisStamp *stamp = stamps->data[s]; // The coordinates of the stamp of interest
-	int x = stamp->x;		// Stamp x coord
-	int y = stamp->y;		// Stamp y coord
-	if (stamp->status == POIS_STAMP_USED) {
-	    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);
-	    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);
-	    psImage *tempMask = psImageCopy(NULL, maskStampTrim, PS_TYPE_U8);
-	    
-	    (void)psImageStats(stats, tempImage, tempMask, POIS_MASK_BAD | POIS_MASK_NEAR_BAD);
-	    
-	    psFree(tempImage);
-	    psFree(tempMask);
-	    
-	    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]);
-	    
+        poisStamp *stamp = stamps->data[s]; // The coordinates of the stamp of interest
+        int x = stamp->x;               // Stamp x coord
+        int y = stamp->y;               // Stamp y coord
+        if (stamp->status == POIS_STAMP_USED) {
+            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);
+            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);
+            psImage *tempMask = psImageCopy(NULL, maskStampTrim, PS_TYPE_U8);
+
+            (void)psImageStats(stats, tempImage, tempMask, POIS_MASK_BAD | POIS_MASK_NEAR_BAD);
+
+            psFree(tempImage);
+            psFree(tempMask);
+
+            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]);
+
 #ifdef TESTING
-	    char stampName[MAXCHAR];		// File name for stamp
-	    snprintf(stampName, MAXCHAR, "stamp%d.fits", s);
-	    psFits *stampFile = psFitsAlloc(stampName);
-	    if (!psFitsWriteImage(stampFile, NULL, subStampTrim, 0)) {
-		psErrorStackPrint(stderr, "Unable to write stamp: %s\n", stampName);
-	    }
-	    psFree(stampFile);
+            char stampName[MAXCHAR];            // File name for stamp
+            snprintf(stampName, MAXCHAR, "stamp%d.fits", s);
+            psFits *stampFile = psFitsOpen(stampName, "w");
+            if (!psFitsWriteImage(stampFile, NULL, subStampTrim, 0)) {
+                psErrorStackPrint(stderr, "Unable to write stamp: %s\n", stampName);
+            }
+            psFitsClose(stampFile);
 #endif
 
 #if 0
-	    psFree(convRefStamp);
-	    psFree(maskStampTrim);
-	    psFree(subStampTrim);
-	    psFree(maskStamp);
-	    psFree(refStamp);
-	    psFree(inStamp);
+            psFree(convRefStamp);
+            psFree(maskStampTrim);
+            psFree(subStampTrim);
+            psFree(maskStamp);
+            psFree(refStamp);
+            psFree(inStamp);
 #endif
-	}
+        }
     }