Index: /trunk/archive/scripts/src/phase2/pmFPAMorph.c
===================================================================
--- /trunk/archive/scripts/src/phase2/pmFPAMorph.c	(revision 5106)
+++ /trunk/archive/scripts/src/phase2/pmFPAMorph.c	(revision 5106)
@@ -0,0 +1,468 @@
+#include <stdio.h>
+#include <string.h>
+#include <assert.h>
+#include "pslib.h"
+#include "papmodule.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+
+
+// Flip the chip in x
+static psImage *xFlipImage(psImage *image, // Image to be flipped
+			   int size	// Size of image in flip dimension
+    )
+{
+    psImage *flipped = psImageAlloc(size, image->numRows, image->type.type); // Flipped image
+    switch (image->type.type) {
+      case PS_TYPE_U16:
+	for (int y = 0; y < image->numRows; y++) {
+	    for (int x = 0; x < image->numCols; x++) {
+		flipped->data.U16[y][size - x - 1] = image->data.U16[y][x];
+	    }
+	}
+	break;
+      default:
+	psError(PS_ERR_IO, true, "Image type %x not yet supported for flip.\n", image->type.type);
+    }
+
+    return flipped;
+}
+
+// Flip the chip in y
+static psImage *yFlipImage(psImage *image, // Image to be flipped
+			   int size	// Size of image in flip dimension
+    )
+{
+    psImage *flipped = psImageAlloc(image->numCols, size, image->type.type); // Flipped image
+    switch (image->type.type) {
+      case PS_TYPE_U16:
+	for (int y = 0; y < image->numRows; y++) {
+	    for (int x = 0; x < image->numCols; x++) {
+		flipped->data.U16[size - y - 1][x] = image->data.U16[y][x];
+	    }
+	}
+	break;
+      default:
+	psError(PS_ERR_IO, true, "Image type %x not yet supported for flip.\n", image->type.type);
+    }
+
+    return flipped;
+}
+
+// Return a list of the region strings and the method for interpreting them
+static psList *getRegions(psString *method, // Method for evaluating the regions
+			  psString methodValues // The METHOD:VALUES string
+    )
+{
+    char *values = strchr(methodValues, ':'); // The VALUES
+    *method = psStringNCopy(methodValues, strlen(methodValues) - strlen(values)); // The METHOD
+    psList *regionStrings = NULL;	// List of the region strings
+    if (strncmp(*method, "HEADER", 6) == 0 || strncmp(*method, "HD", 2) == 0) {
+	regionStrings = papSplit(values, ", ");
+    } else if (strncmp(*method, "VALUE", 5) == 0 || strncmp(*method, "VAL", 3) == 0) {
+	regionStrings = papSplit(values, "; ");
+    }
+
+    return regionStrings;
+}
+
+
+static psArray *spliceCells(psList *outCells, // List of target cells (required for parity info)
+			    psList *inCells, // List of cells to splice together
+			    bool posDep	// Position dependent placement of overscans?
+    )
+{
+    assert(outCells);
+    assert(inCells);
+
+    psTrace(__func__, 5, "Splicing together all cells in the bucket.\n");
+
+    if (inCells->size != outCells->size) {
+	psError(PS_ERR_IO, true, "Sizes of input (%d) and output (%d) lists of cells don't match.\n",
+		inCells->size, outCells->size);
+	return NULL;
+    }
+
+    int numCells = inCells->size;	// Number of cells
+    int readdir = 0;			// Read direction for CCD; currently unknown
+    int numReadouts = 0;		// Number of readouts in a cell
+    psArray *spliced = NULL;		// Array of spliced readouts
+    psElemType type = 0;		// Image type (U16, S32, F32, etc)
+
+    psVector *xSize = NULL;		// Size of spliced image in x
+    psVector *ySize = NULL;		// Size of spliced image in y
+
+    psVector *xFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Do we need to flip in x?
+    psVector *yFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Do we need to flip in y?
+
+    psListIterator *inCellsIter = psListIteratorAlloc(inCells, PS_LIST_HEAD, false); // Iterator for cells
+    papCell *inCell = NULL;		// Cell from list of input cells
+    psListIterator *outCellsIter = psListIteratorAlloc(outCells, PS_LIST_HEAD, false); // Iterator for cells
+    int cellNum = 0;			// Cell number; 
+    while (inCell = psListGetAndIncrement(inCellsIter)) {
+	psArray *readouts = inCell->readouts;	// The readouts comprising the cell
+
+	// Check the read direction
+	int cellReaddir = psMetadataLookupS32(NULL, inCell->concepts, "CELL.READDIR");
+	if (readdir == 0) {
+	    readdir = cellReaddir;
+	} else if (readdir != cellReaddir) {
+	    psError(PS_ERR_IO, true, "Trying to splice cells with different read directions (%d vs %d)\n",
+		    readdir, cellReaddir);
+	    // XXX Clean up before returning
+	    return NULL;
+	}
+	// Check the number of readouts
+	if (! spliced) {
+	    numReadouts = readouts->n;
+	    spliced = psArrayAlloc(numReadouts);
+	    xSize = psVectorAlloc(numReadouts, PS_TYPE_S32);
+	    ySize = psVectorAlloc(numReadouts, PS_TYPE_S32);
+	    for (int i = 0; i < numReadouts; i++) {
+		xSize->data.S32[i] = 0;
+		ySize->data.S32[i] = 0;
+	    }
+	} else if (readouts->n != numReadouts) {
+	    psError(PS_ERR_IO, true, "Trying to splice cells with different number of reads (%d vs %d)\n",
+		    numReadouts, readouts->n);
+	    // XXX Clean up before returning
+	    return NULL;
+	}
+
+	// Get the cell parities
+	papCell *outCell = psListGetAndIncrement(outCellsIter);	// Corresponding output cell
+	int xParityIn = psMetadataLookupS32(NULL, inCell, "CELL.XPARITY");
+	int yParityIn = psMetadataLookupS32(NULL, inCell, "CELL.YPARITY");
+	int xParityOut = psMetadataLookupS32(NULL, outCell, "CELL.XPARITY");
+	int yParityOut = psMetadataLookupS32(NULL, outCell, "CELL.YPARITY");
+	if (xParityIn == 0 || xParityOut == 0 || yParityIn == 0 || yParityOut == 0) {
+	    psError(PS_ERR_IO, false, "Unable to determine parity of cell!\n");
+	    // XXX Clean up before returning
+	    return NULL;
+	}
+
+	if (yParityIn == yParityOut) {
+	    yFlip->data.U8[cellNum] = 0;
+	} else {
+	    psTrace(__func__, 7, "Need to flip y.\n");
+	    yFlip->data.U8[cellNum] = 1;
+	}
+	if (xParityIn == xParityOut) {
+	    xFlip->data.U8[cellNum] = 0;
+	} else {
+	    psTrace(__func__, 7, "Need to flip x.\n");
+	    xFlip->data.U8[cellNum] = 1;
+	}
+	cellNum++;
+
+	// Calculate the sizes of the spliced images
+	for (int i = 0; i < numReadouts; i++) {
+	    papReadout *readout = readouts->data[i]; // The readout of interest
+	    psImage *image = readout->image; // The image pixels
+	    psList *overscans = readout->overscans; // The overscan regions
+
+	    // Check image type
+	    if (type == 0) {
+		type = image->type.type;
+	    } else if (image->type.type != type) {
+		psError(PS_ERR_IO, true, "Trying to splice readouts with different image types (%d vs %d)\n",
+			type, image->type.type);
+		// XXX Clean up before returning
+		return NULL;
+	    }
+
+	    // Get the size of the spliced image
+	    if (readdir == 1) {
+		psTrace(__func__, 9, "Image size: %dx%d\n", image->numCols, image->numRows);
+		xSize->data.S32[i] += image->numCols;
+		ySize->data.S32[i] = MAX(ySize->data.S32[i], image->numRows);
+	    } else if (readdir == 2) {
+		psTrace(__func__, 9, "Image size: %dx%d\n", image->numCols, image->numRows);
+		ySize->data.S32[i] += image->numRows;
+		xSize->data.S32[i] = MAX(xSize->data.S32[i], image->numCols);
+	    } else {
+		psError(PS_ERR_IO, true, "Invalid read direction: %d\n", readdir);
+		// XXX Clean up before returning
+		return NULL;
+	    }
+
+	    psListIterator *overscansIter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false);
+	    psImage *overscan = NULL;	// Overscan image
+	    while (overscan = psListGetAndIncrement(overscansIter)) {
+
+		// Check image type
+		if (type == 0) {
+		    type = image->type.type;
+		} else if (image->type.type != type) {
+		    psError(PS_ERR_IO, true, "Trying to splice readouts with different image types "
+			    "(%d vs %d)\n", type, image->type.type);
+		    // XXX Clean up before returning
+		    return NULL;
+		}
+
+		// Get the size of the spliced image
+		if (readdir == 1) {
+		    psTrace(__func__, 9, "Overscan size: %dx%d\n", image->numCols, image->numRows);
+		    xSize->data.S32[i] += overscan->numCols;
+		    ySize->data.S32[i] = MAX(ySize->data.S32[i], overscan->numRows);
+		} else if (readdir == 2) {
+		    psTrace(__func__, 9, "Overscan size: %dx%d\n", image->numCols, image->numRows);
+		    ySize->data.S32[i] += overscan->numRows;
+		    xSize->data.S32[i] = MAX(xSize->data.S32[i], overscan->numCols);
+		}
+	    }
+
+	}
+    }
+    psFree(inCellsIter);
+    psFree(outCellsIter);
+
+    // Make sure all the readouts have the same size
+    int numRows = ySize->data.S32[0];	// Number of rows for spliced image
+    int numCols = xSize->data.S32[0];	// Number of columns for spliced image
+    for (int i = 1; i < numReadouts; i++) {
+	if (xSize->data.S32[i] != numCols || ySize->data.S32[i] != numRows) {
+	    psError(PS_ERR_IO, true, "Spliced readouts would have different sizes: %dx%d vs %dx%d\n",
+		    numCols, numRows, xSize->data.S32[i], ySize->data.S32[i]);
+	}
+    }
+
+    psTrace(__func__, 5, "Spliced image has dimensions %dx%d\n", numCols, numRows);
+
+    // Now we have done the requisite checks, and know the sizes; we just go through and splice together
+    for (int i = 0; i < numReadouts; i++) {
+	psImage *splice = psImageAlloc(numCols, numRows, type); // The spliced image
+
+	psListIterator *inCellsIter = psListIteratorAlloc(inCells, PS_LIST_HEAD, false);// Iterator for cells
+	papCell *inCell = NULL;		// Cell from the list of cells
+	psListIterator *outCellsIter = psListIteratorAlloc(outCells, PS_LIST_HEAD, false); // Iterator
+	int cellNum = 0;		// Cell number
+	int x0 = 0, y0 = 0;		// Position at which the splice begins
+	while (inCell = psListGetAndIncrement(inCellsIter)) {
+	    psArray *readouts = inCell->readouts; // The readouts comprising the cell
+	    papReadout *readout = readouts->data[i]; // The specific readout of interest
+	    psImage *image = readout->image; // The image pixels
+	    psList *overscans = readout->overscans; // The overscan regions
+	    papCell *outCell = psListGetAndIncrement(outCellsIter); // Output cell
+	    
+	    psImage *toSplice = psMemIncrRefCounter(image); // Image to splice in
+	    if (xFlip->data.U8[cellNum]) {
+		int size = 0;		// Size of flipped image
+		if (readdir == 1) {
+		    size = toSplice->numCols;
+		} else if (readdir == 2) {
+		    size = numCols;
+		}
+		psImage *temp = xFlipImage(toSplice, size);
+		psFree(toSplice);
+		toSplice = temp;
+	    }
+	    if (yFlip->data.U8[cellNum]) {
+		int size = 0;		// Size of flipped image
+		if (readdir == 1) {
+		    size = numRows;
+		} else if (readdir == 2) {
+		    size = toSplice->numRows;
+		}
+		psImage *temp = yFlipImage(toSplice, size);
+		psFree(toSplice);
+		toSplice = temp;
+	    }
+
+	    (void)psImageOverlaySection(splice, toSplice, x0, y0, "=");
+
+	    if (readdir == 1) {
+		x0 += toSplice->numCols;
+	    } else if (readdir == 2) {
+		y0 += toSplice->numRows;
+	    }
+
+	    // Update the TRIMSEC for the output cell
+	    psString trimsecString = psMetadataLookupString(NULL, outCell, "CELL.TRIMSEC");
+	    psString trimsecMethod = NULL; // Method of determining trimsec
+	    psList *trimsecs = getRegions(&trimsecMethod, trimsecString); // List of trimsecs
+	    if (strncmp(trimsecMethod, "HEADER", 6) != 0 && strncmp(trimsecMethod, "HD", 2) != 0) {
+		psError(PS_ERR_IO, true, "Can only splice cells if the CELL.TRIMSEC (= %s) is determined "
+			"from the header.\n", trimsecString);
+		// XXX Clean up before returning
+		return NULL;
+	    }
+	    if (trimsecs->size != 1) {
+		psError(PS_ERR_IO, true, "Multiple headers specified for CELL.TRIMSEC: %s\n", trimsecString);
+		// XXX Clean up before returning
+		return NULL;
+	    }
+	    psString trimsecName = trimsecs->head->data; // Grab header for TRIMSEC off the list
+	    psFree(trimsecs);
+
+	    // TRIMSEC region; add 1 to get FITS standard
+	    psRegion trimsecRegion = {x0 + 1, x0 + toSplice->numCols, y0 + 1, y0 + toSplice->numRows};
+	    psString trimsecValue = psRegionToString(trimsecRegion);
+	    psMetadataAddStr(header, PS_LIST_TAIL, trimsecName, 0, "TRIMSEC from pmFPAMorph", trimsecValue);
+	    // XXX Does psMetadataAddStr overwrite a previous value?
+
+	    // Repeat the above for each overscan (flip, overlay, update BIASSEC)
+	    // Prepare the BIASSEC
+	    psString biassecString = psMetadataLookupString(NULL, outCell, "CELL.BIASSEC");
+	    psString biassecMethod = NULL; // Method of determining biassec
+	    psList *biassecs = getRegions(&biassecMethod, biassecString); // List of biassecs
+	    if (strncmp(biassecMethod, "HEADER", 6) != 0 && strncmp(biassecMethod, "HD", 2) != 0) {
+		psError(PS_ERR_IO, true, "Can only splice cells if the CELL.BIASSEC (= %s) is determined "
+			"from the header.\n", biassecString);
+		// XXX Clean up before returning
+		return NULL;
+	    }
+	    if (biassecs->size != overscans->size);
+
+	    psListIterator *overscansIter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false);
+	    psImage *overscan = NULL;	// Overscan from list
+	    psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false);
+	    psString biassecName = NULL;// FITS header keyword for biassec
+	    while ((overscan = psListGetAndIncrement(overscansIter)) &&
+		   (biassecName = psListGetAndIncrement(biassecsIter))) {
+		// Splice the overscan on
+		psImage *toSplice = psMemIncrRefCounter(overscan); // Image to splice in
+		if (xFlip->data.U8[cellNum]) {
+		    int size = 0;	// Size of flipped image
+		    if (readdir == 1) {
+			size = toSplice->numCols;
+		    } else if (readdir == 2) {
+			size = numCols;
+		    }
+		    psImage *temp = xFlipImage(toSplice, size);
+		    psFree(toSplice);
+		    toSplice = temp;
+		}
+		if (yFlip->data.U8[cellNum]) {
+		    int size = 0;	// Size of flipped image
+		    if (readdir == 1) {
+			size = numRows;
+		    } else if (readdir == 2) {
+			size = toSplice->numRows;
+		    }
+		    psImage *temp = yFlipImage(toSplice, size);
+		    psFree(toSplice);
+		    toSplice = temp;
+		}
+		
+		(void)psImageOverlaySection(splice, toSplice, x0, y0, "=");
+
+		if (readdir == 1) {
+		    x0 += toSplice->numCols;
+		} else if (readdir == 2) {
+		    y0 += toSplice->numRows;
+		}
+
+		// Update the header
+		psRegion biassecRegion = {x0 + 1, x0 + toSplice->numCols, y0 + 1, y0 + toSplice->numRows};
+		psString biassecValue = psRegionToString(biassecRegion);
+		psMetadataAddStr(header, PS_LIST_TAIL, biassecName, 0, "BIASSEC from pmFPAMorph",
+				 biassecValue);
+	    }
+	    psFree(overscansIter);
+	    psFree(biassecsIter);
+	    psFree(biassecs);
+	} // Iterating over input cells
+
+	spliced->data[i] = splice;
+    } // Iterating over readouts
+
+    return spliced;
+}
+
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+bool pmFPAMorph(pmFPA *toFPA,		// FPA structure to which to morph
+		pmFPA *fromFPA,		// FPA structure from which to morph
+		bool positionDependent, // Is the position of the overscan dependent on the position?
+		int chipNum,		// Chip number, in case the scopes are different
+		int cellNum		// Cell number, in case the scopes are different
+    )
+{
+    psList *targetCells = psListAlloc(NULL); // List of target cells
+    psList *sourceCells = psListAlloc(NULL); // List of source cells
+
+    psArray *toChips = toFPA->chips;	// Array of chips
+    psArray *fromChips = fromFPA->chips;// Array of chips    
+
+    for (int i = 0; i < toChips->n; i++) {
+	pmChip *toChip = toChips->data[i];
+	pmChip *fromChip = NULL;
+	// Select the correct chip
+	if (toChips->n == 1 && fromChips->n > chipNum) {
+	    fromChip = fromChips->data[chipNum];
+	} else if (fromChips->n == 1 && toChips->n > chipNum) {
+	    fromChip = fromChips->data[0];
+	} else if (toChips->n == fromChips->n) {
+	    fromChip = fromChips->data[i];
+	} else {
+	    psError(PS_ERR_IO, true, "Unable to discern intended chip.\n");
+	    return false;
+	}
+	
+	psArray *toCells = toChip->cells; // Array of cells
+	psArray *fromCells = fromChip->cells; // Array of cells
+
+	for (int j = 0; j < toCells->n; j++) {
+	    pmCell *toCell = toCells->data[j];
+	    pmCell *fromCell = NULL;
+	    // Select the correct cell
+	    if (toCells->n == 1 && fromCells->n > cellNum) {
+		fromCell = fromCells->data[cellNum];
+	    } else if (fromCells->n == 1 && toCells->n > cellNum) {
+		fromCell = fromCells->data[0];
+	    } else if (toCells->n == fromCells->n) {
+		fromCell = fromCells->data[j];
+	    } else {
+		psError(PS_ERR_IO, true, "Unable to discern intended cell.\n");
+		return false;
+	    }
+
+	    psTrace(__func__, 5, "Putting chip %d, cell %d in the bucket.\n", i, j);
+	    psListAdd(targetCells, PS_LIST_TAIL, toCell);
+	    psListAdd(sourceCells, PS_LIST_TAIL, fromCell);
+
+	    int xParityIn = psMetadataLookupS32(NULL, fromCell, "CELL.XPARITY");
+	    int yParityIn = psMetadataLookupS32(NULL, fromCell, "CELL.YPARITY");
+	    int xParityOut = psMetadataLookupS32(NULL, toCell, "CELL.XPARITY");
+	    int yParityOut = psMetadataLookupS32(NULL, toCell, "CELL.YPARITY");
+	    psTrace(__func__, 3, "Chip %d, cell %d: in (%d,%d) out (%d,%d)\n", i, j, xParityIn, yParityIn,
+		    xParityOut, yParityOut);
+
+	    // Copy the cell contents over
+	    toCell->readouts = psMemIncrRefCounter(fromCell->readouts);
+	    
+	    if (toCell->private && strlen(toCell->private->extname) > 0) {
+		// Splice the component cells
+		p_pmHDU *hdu = toCell->private;
+		if (! hdu->header) {
+		    hdu->header = psMetadataAlloc();
+		}
+		hdu->pixels = spliceCells(hdu->header, targetCells, sourceCells, positionDependent);
+		// Purge the lists
+		while (psListRemove(targetCells, PS_LIST_TAIL));
+		while (psListRemove(sourceCells, PS_LIST_TAIL));
+	    }
+
+	}
+
+	if (toChip->private && strlen(toChip->private->extname) > 0) {
+	    // Splice the component cells
+	    p_pmHDU *hdu = toChip->private;
+	    if (! hdu->header) {
+		hdu->header = psMetadataAlloc();
+	    }
+	    hdu->pixels = spliceCells(hdu->header, targetCells, sourceCells);
+	    // Purge the lists
+	    while (psListRemove(targetCells, PS_LIST_TAIL));
+	    while (psListRemove(sourceCells, PS_LIST_TAIL));
+	}
+
+    }
+
+    return true;
+}
+
Index: /trunk/archive/scripts/src/phase2/pmFPAMorph.h
===================================================================
--- /trunk/archive/scripts/src/phase2/pmFPAMorph.h	(revision 5106)
+++ /trunk/archive/scripts/src/phase2/pmFPAMorph.h	(revision 5106)
@@ -0,0 +1,12 @@
+#ifndef PM_FPA_MORPH_H
+#define PM_FPA_MORPH_H
+
+// Morph one focal plane representation into another
+bool pmFPAMorph(papFPA *toFPA,		// FPA structure to which to morph
+		papFPA *fromFPA,	// FPA structure from which to morph
+		int chipNum,		// Chip number, in case the scopes are different
+		int cellNum		// Cell number, in case the scopes are different
+    );
+
+
+#endif