Index: trunk/stac/src/stacRead.c
===================================================================
--- trunk/stac/src/stacRead.c	(revision 3681)
+++ trunk/stac/src/stacRead.c	(revision 5743)
@@ -1,45 +1,53 @@
 #include <stdio.h>
 #include <string.h>
+#include <assert.h>
 #include "pslib.h"
 #include "stac.h"
 
-#define BUFFER 100			// Size of buffer for incrementally reading coordinates
-
-psArray *stacReadImages(psArray *filenames // The file names of the images
+#define BUFFER 100                      // Size of buffer for incrementally reading coordinates
+
+psArray *stacReadImages(psArray **headers, // The image headers, to be returned
+    psArray *filenames // The file names of the images
     )
 {
-    int nFiles = filenames->n;		// The number of input files
+    int nFiles = filenames->n;          // The number of input files
     psArray *images = psArrayAlloc(nFiles); // The input files, to be returned
-    psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read
+    assert(!headers || ! *headers || (*headers)->n == nFiles);
+    if (headers && ! *headers) {
+        *headers = psArrayAlloc(nFiles);
+    }
+
+    psRegion imageRegion = {0, 0, 0, 0}; // Region of image to read
 
     psTrace("stac.read.images", 1, "Reading input images....\n");
     for (int i = 0; i < nFiles; i++) {
-	psTrace("stac.read.images", 2, "Reading input image %s....\n",filenames->data[i]);
-	psFits *imageFile = psFitsAlloc(filenames->data[i]);
-	// We only read PHUs --- not mucking around with extensions for now
-	psImage *image = psFitsReadImage(NULL, imageFile, *imageRegion, 0);
-	if (image == NULL) {
-	    psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n",filenames->data[i]);
-	    exit(EXIT_FAILURE);
-	}
-	psTrace("stac.read.images", 4, "Image %s is %dx%d\n", filenames->data[i], image->numCols,
-		image->numRows);
-	// Convert to 32-bit floating point, in necessary
-	if (image->type.type != PS_TYPE_F32) {
-	    psTrace("stac.read.images", 3, "Converting %s to floating point in memory....",
-		    filenames->data[i]);
-	    images->data[i] = psImageCopy(NULL, image, PS_TYPE_F32);
-	    psFree(image);
-	} else {
-	    int numNaN = psImageClipNaN(image, FLT_MAX);
-	    if (numNaN) {
-		psTrace("stac.read.images", 5, "Clipped %d NaN pixels.\n", numNaN);
-	    }
-	    images->data[i] = image;
-	}
-	psFree(imageFile);
+        psTrace("stac.read.images", 2, "Reading input image %s....\n",filenames->data[i]);
+        psFits *imageFile = psFitsOpen(filenames->data[i], "r");
+        // We only read PHUs --- not mucking around with extensions for now
+        if (headers) {
+            (*headers)->data[i] = psFitsReadHeader(NULL, imageFile);
+        }
+        psImage *image = psFitsReadImage(NULL, imageFile, imageRegion, 0);
+        if (image == NULL) {
+            psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n",filenames->data[i]);
+            exit(EXIT_FAILURE);
+        }
+        psFitsClose(imageFile);
+        psTrace("stac.read.images", 4, "Image %s is %dx%d\n", filenames->data[i], image->numCols,
+                image->numRows);
+        // Convert to 32-bit floating point, in necessary
+        if (image->type.type != PS_TYPE_F32) {
+            psTrace("stac.read.images", 3, "Converting %s to floating point in memory....\n",
+                    filenames->data[i]);
+            images->data[i] = psImageCopy(NULL, image, PS_TYPE_F32);
+            psFree(image);
+        }
+        int numNaN = psImageClipNaN(image, 0.0);
+        if (numNaN) {
+            psTrace("stac.read.images", 5, "Clipped %d NaN pixels to zero.\n", numNaN);
+        }
+        images->data[i] = image;
     }
     psTrace("stac.read.images",1,"%d input images read.\n",nFiles);
-    psFree(imageRegion);
 
     return images;
@@ -50,6 +58,6 @@
     FILE *file = fopen(filename, "r");
     if (file == NULL) {
-	psLogMsg("stac.read.coords", PS_LOG_ERROR, "Cannot open coordinate file, %s\n", filename);
-	return NULL;
+        psLogMsg("stac.read.coords", PS_LOG_ERROR, "Cannot open coordinate file, %s\n", filename);
+        return NULL;
     }
 
@@ -57,15 +65,15 @@
 
     psArray *coords = psArrayAlloc(BUFFER); // The array of coordinates to be returned
-    float x, y;				// Coordinates to read
-    int num = 0;			// Number of coordinates read
+    float x, y;                         // Coordinates to read
+    int num = 0;                        // Number of coordinates read
     while (fscanf(file, "%f %f\n", &x, &y) == 2) {
-	psPlane *coord = psPlaneAlloc();// A coordinate
-	coord->x = x;
-	coord->y = y;
-	coords->data[num] = coord;
-	num++;
-	if (num % BUFFER) {
-	    coords = psArrayRealloc(coords, num + BUFFER);
-	}
+        psPlane *coord = psPlaneAlloc();// A coordinate
+        coord->x = x;
+        coord->y = y;
+        coords->data[num] = coord;
+        num++;
+        if (num % BUFFER) {
+            coords = psArrayRealloc(coords, num + BUFFER);
+        }
     }
     coords->n = num;
@@ -82,86 +90,86 @@
     FILE *mapfp = fopen(filename, "r");
     if (mapfp == NULL) {
-	psLogMsg("stac.read.map", PS_LOG_ERROR, "Cannot open map file, %s\n", filename);
-	return NULL;
+        psLogMsg("stac.read.map", PS_LOG_ERROR, "Cannot open map file, %s\n", filename);
+        return NULL;
     }
     // Read the file
     psTrace("stac.read.map", 5, "Reading map file %s....\n", filename);
-    
+
     // Format is now:
     // order
     // x coefficients
     // y coefficients
-    // 
+    //
     // where the order is 1 for linear, 2 for quadratic, 3 for cubic.
     // and the coefficients are read by the following pseudo-code:
-    // 
-    // 	for (int k = 0; k < order + 1; k++)
-    // 	    for (int j = 0; j < k; j++)
-    // 		int i = k - j;
+    //
+    //  for (int k = 0; k < order + 1; k++)
+    //      for (int j = 0; j < k; j++)
+    //          int i = k - j;
     //              read coefficient of x^i y^j
-    // 
+    //
     // This produces the ordering:
     // x^0 y^0, x^1 y^0, x^0 y^1, x^2 y^0, x^1 y^1, x^0 y^2, x^3 y^0, x^2 y^1, x^1 y^2, x^0 y^3
-    // 
+    //
     // This is, of course, for third order polynomials.
     // For lower orders, the list is truncated at the appropriate level.
 
-    int order = 0;			// Polynomial order
+    int order = 0;                      // Polynomial order
     if (fscanf(mapfp, "%d", &order) != 1) {
-	psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
-	fclose(mapfp);
-	return NULL;
-    }
-    
+        psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
+        fclose(mapfp);
+        return NULL;
+    }
+
     psTrace("stac.read.map", 5, "Polynomial order: %d\n", order);
-    
+
     psPlaneTransform *map = psPlaneTransformAlloc(order + 1, order + 1); // The transformation
     // Set coefficient masks
     for (int i = 0; i < order + 1; i++) {
-	for (int j = 0; j < order + 1; j++) {
-	    if (i + j > order + 1) {
-		map->x->mask[i][j] = 1;
-		map->y->mask[i][j] = 1;
-	    } else {
-		map->x->mask[i][j] = 0;
-		map->y->mask[i][j] = 0;
-	    }
-	}
-    }
-    
+        for (int j = 0; j < order + 1; j++) {
+            if (i + j > order + 1) {
+                map->x->mask[i][j] = 1;
+                map->y->mask[i][j] = 1;
+            } else {
+                map->x->mask[i][j] = 0;
+                map->y->mask[i][j] = 0;
+            }
+        }
+    }
+
     // Read x coefficients
     psTrace("stac.read.map", 7, "x' = \n");
     for (int k = 0; k < order + 1; k++) {
-	for (int j = 0; j < k + 1; j++) {
-	    int i = k - j;
-	    if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
-		psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
-		fclose(mapfp);
-		psFree(map);
-		return NULL;
-	    }
-	    psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
-	}
+        for (int j = 0; j < k + 1; j++) {
+            int i = k - j;
+            if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
+                psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
+                fclose(mapfp);
+                psFree(map);
+                return NULL;
+            }
+            psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
+        }
     }
     // Read y coefficients
     psTrace("stac.read.maps", 7, "y' = \n");
     for (int k = 0; k < order + 1; k++) {
-	for (int j = 0; j < k + 1; j++) {
-	    int i = k - j;
-	    if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
-		psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
-		fclose(mapfp);
-		psFree(map);
-		return NULL;
-	    }
-	    psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
-	}
-    }
-    
+        for (int j = 0; j < k + 1; j++) {
+            int i = k - j;
+            if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
+                psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
+                fclose(mapfp);
+                psFree(map);
+                return NULL;
+            }
+            psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
+        }
+    }
+
     fclose(mapfp);
 
     return map;
 }
-    
+
 
 
@@ -169,20 +177,20 @@
     )
 {
-    int nFiles = filenames->n;		// The number of input files
+    int nFiles = filenames->n;          // The number of input files
     psArray *maps = psArrayAlloc(nFiles); // The maps, to be returned
-    char mapfile[MAXCHAR];		// Filename of map
-    
+    char mapfile[MAXCHAR];              // Filename of map
+
     psTrace("stac.read.maps",1,"Reading maps....\n");
     for (int i = 0; i < nFiles; i++) {
-	if (strlen(filenames->data[i]) > MAXCHAR - 4) {
-	    psLogMsg("stac.read.maps",PS_LOG_ERROR,"Filename %s is too long.\n",filenames->data[i]);
-	    exit(EXIT_FAILURE);
-	}
-	// Read the file
-	sprintf(mapfile,"%s.map",filenames->data[i]);
-	maps->data[i] = stacReadMap(mapfile);
-	if (maps->data[i] == NULL) {
-	    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map: %s\n", mapfile);
-	}
+        if (strlen(filenames->data[i]) > MAXCHAR - 4) {
+            psLogMsg("stac.read.maps",PS_LOG_ERROR,"Filename %s is too long.\n",filenames->data[i]);
+            exit(EXIT_FAILURE);
+        }
+        // Read the file
+        sprintf(mapfile,"%s.map",filenames->data[i]);
+        maps->data[i] = stacReadMap(mapfile);
+        if (maps->data[i] == NULL) {
+            psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map: %s\n", mapfile);
+        }
 
     }