#include <stdio.h>
#include <strings.h>
#include <assert.h>
#include "pslib.h"

#include "papStuff.h"
#include "pmFPA.h"
#include "pmFPAMorph.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;
      case PS_TYPE_F32:
        for (int y = 0; y < image->numRows; y++) {
            for (int x = 0; x < image->numCols; x++) {
                flipped->data.F32[y][size - x - 1] = image->data.F32[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;
}


#if 0
// 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;
}
#endif


// Splice an image into another image
static psRegion *spliceImage(psImage *target, // Target image onto which to splice
                             int *x0, int *y0, // Starting coordinates to splice to
                             psImage *source, // Source image from which to splice
                             const psRegion *from, // Region to splice from
                             int readdir, // Read direction
                             bool xFlip, bool yFlip // Flip x and y axes?
                        )
{
    psImage *toSplice = psImageSubset(source, *from); // Subimage, to splice into target
    if (xFlip) {
        int size = 0;           // Size of flipped image
        if (readdir == 1) {
            size = toSplice->numCols;
        } else if (readdir == 2) {
            size = target->numCols;
        }
        psImage *temp = xFlipImage(toSplice, size);
        psFree(toSplice);
        toSplice = temp;
    }
    if (yFlip) {
        int size = 0;           // Size of flipped image
        if (readdir == 1) {
            size = target->numRows;
        } else if (readdir == 2) {
            size = toSplice->numRows;
        }
        psImage *temp = yFlipImage(toSplice, size);
        psFree(toSplice);
        toSplice = temp;
    }

    psTrace(__func__, 5, "Splicing %dx%d image in at %d,%d\n", toSplice->numCols, toSplice->numRows,
            *x0, *y0);
    (void)psImageOverlaySection(target, toSplice, *x0, *y0, "=");

    if (readdir == 1) {
        *x0 += toSplice->numCols;
    } else if (readdir == 2) {
        *y0 += toSplice->numRows;
    }

    // Return the region where we pasted the image
    psRegion *outRegion = psAlloc(sizeof(psRegion));
    // Addding 1 to get FITS standard
    outRegion->x0 = *x0 + 1;
    outRegion->x1 = *x0 + toSplice->numCols;
    outRegion->y0 = *y0 + 1;
    outRegion->y1 = *y0 + toSplice->numRows;

    psFree(toSplice);
    return outRegion;
}

// Splice bias regions (overscans)
static bool spliceBias(psImage *splice, // Image to which to splice
                       int *x0, int *y0, // Starting position for splice
                       const pmCell *in, // Input cell
                       pmCell *out,     // Output cell
                       int xFlip, int yFlip, // Flip the image?
                       int readNum,     // Number of readout
                       int readdir      // Read direction
                       )
{
    psArray *readouts = in->readouts; // The readouts;
    pmReadout *readout = readouts->data[readNum]; // The readout of interest
    psList *inBiassecs = psMetadataLookupPtr(NULL, in->concepts, "CELL.BIASSEC"); // List of input biassecs
    psListIterator *inBiassecsIter = psListIteratorAlloc(inBiassecs, PS_LIST_HEAD, false);
    psRegion *inBiassec = NULL; // BIASSEC from list
    psList *outBiassecs = psListAlloc(NULL); // List of biassecs for output
    while (inBiassec = psListGetAndIncrement(inBiassecsIter)) {
        psRegion *outBiassec = spliceImage(splice, x0, y0, readout->image, inBiassec, readdir, xFlip,
                                           yFlip);
        psListAdd(outBiassecs, PS_LIST_TAIL, outBiassec);
    }
    psFree(inBiassecsIter);
    psMetadataAdd(out->concepts, PS_LIST_TAIL, "CELL.BIASSEC", PS_DATA_LIST | PS_META_REPLACE,
                  "BIASSEC from pmFPAMorph", outBiassecs);

    return true;
}

psImage *p_pmImageMosaic(const psArray *source, // Images to splice in
                         const psVector *xFlip, const psVector *yFlip, // Need to flip x and y?
                         const psVector *xBinSource, const psVector *yBinSource, // Binning in x and y of
                                                                                 // source images
                         int xBinTarget, int yBinTarget, // Binning in x and y of target images
                         const psVector *x0, const psVector *y0 // Offsets for source images on target
                         );

// Splice a list of cells together
// We need to do the following:
// 1. Copy the contents of the cells (i.e., the readouts) over from the in to the out
// 2. Splice together everything in the bucket for the image that we will write out
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->n != outCells->n) {
        psError(PS_ERR_IO, true, "Sizes of input (%d) and output (%d) lists of cells don't match.\n",
                inCells->n, outCells->n);
        return NULL;
    }

    int numCells = inCells->n;          // 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 *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?
    psVector *xBin = psVectorAlloc(numCells, PS_TYPE_S32); // Binning in x
    psVector *yBin = psVectorAlloc(numCells, PS_TYPE_S32); // Binning in y
    psVector *x0 = psVectorAlloc(numCells, PS_TYPE_S32); // Offset in x
    psVector *y0 = psVectorAlloc(numCells, PS_TYPE_S32); // Offset in y
    int xBinOut = 0;                    // Binning in x for output
    int yBinOut = 0;                    // Binning in y for output

    // We go through the cells to make sure everything's kosher, and to get the sizes.
    psListIterator *inCellsIter = psListIteratorAlloc(inCells, PS_LIST_HEAD, false); // Iterator for cells
    psListIterator *outCellsIter = psListIteratorAlloc(outCells, PS_LIST_HEAD, false); // Iterator for cells
    pmCell *inCell = NULL;              // Cell from list of input cells
    pmCell *outCell = NULL;             // Cell from list of output cells
    int cellNum = -1;                   // Cell number
    while (inCell = psListGetAndIncrement(inCellsIter) && outCell = psListGetAndIncrement(outCellsIter)) {
        cellNum++;
        psArray *inReadouts = inCell->readouts;   // The readouts comprising the cell

        // Get and check the read direction
        int cellReaddir = psMetadataLookupS32(NULL, inCell->concepts, "CELL.READDIR");
        psTrace(__func__, 10, "Checking read direction for cell %d: %d\n", cellNum, cellReaddir);
        if (cellNum == 0) {
            // First run through; set it.
            readdir = cellReaddir;
        } else if (readdir != cellReaddir) {
            psError(PS_ERR_IO, true, "Trying to splice cells with different read directions (%d vs %d). I "
                    "don't know how to do that!\n", readdir, cellReaddir);
            // XXX Clean up before returning
            return NULL;
        }

        // Get and check the number of readouts
        if (! spliced) {
            numReadouts = inReadouts->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 (inReadouts->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 and check the binning --- we can't splice to cells that have different binnings, because
        // the spliced image can't be defined.  Hence all the "to" cells must have the same binning.
        xBinIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.XBIN");
        yBinIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.YBIN");
        int xBinCheck = psMetadataLookupS32(NULL, outCell->concepts, "CELL.XBIN");
        int yBinCheck = psMetadataLookupS32(NULL, outCell->concepts, "CELL.YBIN");
        if (xBinOut == 0) {
            xBinOut = xBinCheck;
        } else if (xBinCheck != xBinOut) {
            psError(PS_ERR_IO, true, "Trying to splice to cells with different x binnings (%d vs %d)\n",
                    xBinOut, xBinCheck);
            // XXX Clean up before returning
            return NULL;
        }
        if (yBinOut == 0) {
            yBinOut = yBinCheck;
        } else if (yBinCheck != yBinOut) {
            psError(PS_ERR_IO, true, "Trying to splice to cells with different y binnings (%d vs %d)\n",
                    yBinOut, yBinCheck);
            // XXX Clean up before returning
            return NULL;
        }

        // Copy the readouts over
        psArray *outReadouts = psArrayAlloc(inReadouts->n); // Array of readouts for the output cell
        for (int i = 0; i < inReadouts->n; i++) {
            outReadouts->data[i] = psMemIncrRefCounter(inReadouts->data[i]);
        }

        // Get the offsets
        x0->data.S32[i] = psMetadataLookupS32(NULL, inCell->concepts, "CELL.X0");
        y0->data.S32[i] = psMetadataLookupS32(NULL, inCell->concepts, "CELL.Y0");

        // Get and check the parity
        int xParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.XPARITY");
        int yParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.YPARITY");
        int xParityOut = psMetadataLookupS32(NULL, outCell->concepts, "CELL.XPARITY");
        int yParityOut = psMetadataLookupS32(NULL, outCell->concepts, "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 (xParityIn != xParityOut) {
            for (int i = 0; i < outReadouts->n; i++) {
                psImage *image = outReadouts->data[i];
                outReadouts->data[i] = xFlipImage(image);
                psFree(image);
            }
            // Correct CELL.X0 for the parity shift
            psMetadataItem *x0item = psMetadataLookup(inCell->concepts, "CELL.X0");
            x0item->data.S32 -= xSize->data.S32[cellNum];
        }
        if (yParityIn != yParityOut) {
            for (int i = 0; i < outReadouts->n; i++) {
                psImage *image = outReadouts->data[i];
                outReadouts->data[i] = xFlipImage(image);
                psFree(image);
            }
            // Correct CELL.Y 0 for the parity shift
            psMetadataItem *y0item = psMetadataLookup(inCell->concepts, "CELL.Y0");
            y0item->data.S32 -= ySize->data.S32[cellNum];
        }

        //////////////////////////////////////////////////////////////////////////////////////////////////////
        // DO I WANT TO DO THIS????
        // Wouldn't it be better to splice all the images together first, and then divvy up the components?
        //////////////////////////////////////////////////////////////////////////////////////////////////////


        // Manipulate the images, if required
        for (int i = 0; i < outReadouts->n; i++) {
            psImage *image = manipulate(outReadouts

                         psImage *image = psMemIncrRefCounter(outReadouts->data[i]);
                         if (xParityIn != xParityOut) {
                             psImage *image = xFlipImage(image);
                         }


        // THIS ISN'T SO SIMPLE: we need to SET the output x0,y0 based on where the image actually gets pasted
        // Get the offsets
        x0->data.S32[i] = psMetadataLookupS32(NULL, outCell->concepts, "CELL.X0");
        y0->data.S32[i] = psMetadataLookupS32(NULL, outCell->concepts, "CELL.Y0");




// Splice a list of cells together
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->n != outCells->n) {
        psError(PS_ERR_IO, true, "Sizes of input (%d) and output (%d) lists of cells don't match.\n",
                inCells->n, outCells->n);
        return NULL;
    }

    int numCells = inCells->n;          // 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?


    // We go through the cells to make sure everything's kosher, and to get the sizes.
    psListIterator *inCellsIter = psListIteratorAlloc(inCells, PS_LIST_HEAD, false); // Iterator for cells
    psListIterator *outCellsIter = psListIteratorAlloc(outCells, PS_LIST_HEAD, false); // Iterator for cells
    pmCell *inCell = NULL;              // Cell from list of input cells
    int cellNum = 0;                    // Cell number;
    while (inCell = psListGetAndIncrement(inCellsIter)) {
        psArray *readouts = inCell->readouts;   // The readouts comprising the cell

        // Get and check the read direction
        int cellReaddir = psMetadataLookupS32(NULL, inCell->concepts, "CELL.READDIR");
        psTrace(__func__, 10, "Checking read direction for cell %d: %d\n", cellNum, cellReaddir);
        if (cellNum == 0) {
            // First run through; set it.
            readdir = cellReaddir;
        } else if (readdir != cellReaddir) {
            psError(PS_ERR_IO, true, "Trying to splice cells with different read directions (%d vs %d). I "
                    "don't know how to do that!\n", readdir, cellReaddir);
            // XXX Clean up before returning
            return NULL;
        }

        // Get and 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 and check the cell parities
        pmCell *outCell = psListGetAndIncrement(outCellsIter);  // Corresponding output cell
        int xParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.XPARITY");
        int yParityIn = psMetadataLookupS32(NULL, inCell->concepts, "CELL.YPARITY");
        int xParityOut = psMetadataLookupS32(NULL, outCell->concepts, "CELL.XPARITY");
        int yParityOut = psMetadataLookupS32(NULL, outCell->concepts, "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++) {
            pmReadout *readout = readouts->data[i]; // The readout of interest
            psImage *image = readout->image; // The image pixels

            // 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;
            }

            psRegion *trimsec = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.TRIMSEC"); // Trim section
            psList *biassecs = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.BIASSEC"); // Bias section
            psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false); // Iterator
            psRegion *biassec = NULL;   // Bias section from list iteration

            // Get the size of the spliced image
            if (readdir == 1) {
                psTrace(__func__, 7, "TRIMSEC is %.0fx%.0f\n", trimsec->x1 - trimsec->x0,
                        trimsec->y1 - trimsec->y0);
                xSize->data.S32[i] += trimsec->x1 - trimsec->x0;
                ySize->data.S32[i] = MAX(ySize->data.S32[i], trimsec->y1 - trimsec->y0);
                while (biassec = psListGetAndIncrement(biassecsIter)) {
                    psTrace(__func__, 9, "BIASSEC is %.0fx%.0f\n", biassec->x1 - biassec->x0,
                        biassec->y1 - biassec->y0);
                    xSize->data.S32[i] += biassec->x1 - biassec->x0;
                    ySize->data.S32[i] = MAX(ySize->data.S32[i], biassec->y1 - biassec->y0);
                }
            } else if (readdir == 2) {
                ySize->data.S32[i] += trimsec->y1 - trimsec->y0;
                xSize->data.S32[i] = MAX(xSize->data.S32[i], trimsec->x1 - trimsec->x0);
                while (biassec = psListGetAndIncrement(biassecsIter)) {
                    ySize->data.S32[i] += biassec->y1 - biassec->y0;
                    xSize->data.S32[i] = MAX(ySize->data.S32[i], biassec->x1 - biassec->x0);
                }
            } else {
                psError(PS_ERR_IO, true, "Invalid read direction: %d\n", readdir);
                // XXX Clean up before returning
                return NULL;
            }
            psFree(biassecsIter);
        }
    }
    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
        int x0 = 0, y0 = 0;             // Position at which the splice begins

        // Position-dependent overscans first
        if (posDep) {
            psTrace(__func__, 8, "Doing first position-dependent biases...\n");
            for (int cellNum = 0; cellNum < inCells->n / 2; cellNum++) {
                psTrace(__func__, 9, "Cell %d...\n", cellNum);
                pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest
                pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest
                spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i,
                           readdir);
            }
        }

        // Then the images
        psListIterator *inCellsIter = psListIteratorAlloc(inCells, PS_LIST_HEAD, false);// Iterator for cells
        pmCell *inCell = NULL;          // Cell from the list of cells
        psListIterator *outCellsIter = psListIteratorAlloc(outCells, PS_LIST_HEAD, false); // Iterator
        int cellNum = 0;                // Cell number
        while (inCell = psListGetAndIncrement(inCellsIter)) {
            psArray *readouts = inCell->readouts; // The readouts comprising the cell
            pmReadout *readout = readouts->data[i]; // The specific readout of interest
            pmCell *outCell = psListGetAndIncrement(outCellsIter); // Output cell
            psRegion *inTrimsec = psMetadataLookupPtr(NULL, inCell->concepts, "CELL.TRIMSEC"); // TRIMSEC
            psRegion *outTrimsec = spliceImage(splice, &x0, &y0, readout->image, inTrimsec, readdir,
                                               xFlip->data.U8[cellNum], yFlip->data.U8[cellNum]);

            // Update the output TRIMSEC
            psMetadataAdd(outCell->concepts, PS_LIST_TAIL, "CELL.TRIMSEC", PS_DATA_UNKNOWN | PS_META_REPLACE,
                             "TRIMSEC from pmFPAMorph", outTrimsec);
        }
        psFree(inCellsIter);
        psFree(outCellsIter);

        // Then the biases
        if (posDep) {
            // Need to only do half the biases
            psTrace(__func__, 8, "Doing second position-dependent biases...\n");
            for (int cellNum = inCells->n / 2; cellNum < inCells->n; cellNum++) {
                psTrace(__func__, 9, "Cell %d...\n", cellNum);
                pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest
                pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest
                spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i,
                           readdir);
            }
        } else {
            // Need to do all the biases
            for (int cellNum = 0; cellNum < inCells->n; cellNum++) {
                pmCell *inCell = psListGet(inCells, cellNum); // Input cell of interest
                pmCell *outCell = psListGet(outCells, cellNum); // Output cell of interest
                spliceBias(splice, &x0, &y0, inCell, outCell, xFlip->data.U8[i], yFlip->data.U8[i], i,
                           readdir);
            }
        }

        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

    psTrace(__func__, 1, "Copying FPA concepts over...\n");
    pap_psMetadataCopy(toFPA->concepts, fromFPA->concepts);
    pap_psMetadataCopy(toFPA->phu, fromFPA->phu);

    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;
        }

        if (! fromChip->valid) {
            toChip->valid = false;
            continue;
        }

        psTrace(__func__, 2, "Copying chip %d concepts over...\n", i);
        pap_psMetadataCopy(toChip->concepts, fromChip->concepts);

        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;
            }

            if (! fromCell->valid) {
                toCell->valid = false;
                continue;
            }

#ifdef TESTING
            int xParityIn = psMetadataLookupS32(NULL, fromCell->concepts, "CELL.XPARITY");
            int yParityIn = psMetadataLookupS32(NULL, fromCell->concepts, "CELL.YPARITY");
            int xParityOut = psMetadataLookupS32(NULL, toCell->concepts, "CELL.XPARITY");
            int yParityOut = psMetadataLookupS32(NULL, toCell->concepts, "CELL.YPARITY");
            psTrace(__func__, 3, "Chip %d, cell %d: in (%d,%d) out (%d,%d)\n", i, j, xParityIn, yParityIn,
                    xParityOut, yParityOut);
#endif

#ifdef TESTING
            psMetadataPrint(toCell->concepts, 7);
#endif

            // Copy the concepts over

            // Need to be a little tricky here --- some concepts are already set from the camera configuration
            // file (e.g., CELL.NAME), and we want to preserve these, not replace it with the old value.
            psTrace(__func__, 3, "Copying cell %d concepts over...\n", j);
            psMetadata *newConcepts = pap_psMetadataCopy(NULL, fromCell->concepts);
            pap_psMetadataCopy(newConcepts, toCell->concepts); // This preserves any already existing concepts
            psFree(toCell->concepts);
            toCell->concepts = newConcepts;

            // Now, we need to check if the following concepts for the target cell are defined statically:
            // CELL.XPARITY, CELL.YPARITY, CELL.XBIN, CELL.YBIN, CELL.X0, CELL.Y0
            // 1. If they are specified by the header, we can do what we want (because we can set the header).
            // 2. If they are not specified by the header, we use those values.
            const char *checkConcepts[4] = { "CELL.XPARITY",
                                             "CELL.YPARITY",
                                             "CELL.XBIN",
                                             "CELL.YBIN",
                                             "CELL.X0",
                                             "CELL.Y0" }; // Concepts that we need to check
            bool mdok = false;          // Result of MD lookup
            psMetadata *translation = psMetadataLookupMD(&mdok, toFPA->camera, "TRANSLATION"); // Header
                                        // translation table
            if (!mdok || ! translation) {
                psError(PS_ERR_IO, true, "Unable to find TRANSLATION in camera configuration!\n");
                return false;
            }
            for (int c = 0; c < 4; c++) {
                psString headerName = psMetadataLookupString(&mdok, translation, checkConcepts[c]); // Name of
                                        // header holding the concept, or NULL
                if (!mdok || strlen(headerName) <= 0) {
                    p_pmCellIngestConcept(toFPA, toChip, toCell, db, checkConcepts[c]);
                }
            }

#ifdef TESTING
            psMetadataPrint(toCell->concepts, 7);
#endif

            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);

            // Copy the cell contents over
            toCell->readouts = psMemIncrRefCounter(fromCell->readouts);

            if (toCell->hdu && strlen(toCell->hdu->extname) > 0) {
                // Splice the component cells
                p_pmHDU *hdu = toCell->hdu;
                if (! hdu->header) {
                    hdu->header = psMetadataAlloc();
                }
                hdu->images = spliceCells(targetCells, sourceCells, positionDependent);
                // Purge the lists
                while (psListRemove(targetCells, PS_LIST_TAIL));
                while (psListRemove(sourceCells, PS_LIST_TAIL));
            }

        }

        if (toChip->hdu && strlen(toChip->hdu->extname) > 0) {
            // Splice the component cells
            p_pmHDU *hdu = toChip->hdu;
            if (! hdu->header) {
                hdu->header = psMetadataAlloc();
            }
            hdu->images = spliceCells(targetCells, sourceCells, positionDependent);
            // Purge the lists
            while (psListRemove(targetCells, PS_LIST_TAIL));
            while (psListRemove(sourceCells, PS_LIST_TAIL));
        }

    }

    if (toFPA->hdu && strlen(toFPA->hdu->extname) > 0) {
        // Splice the component cells
        p_pmHDU *hdu = toFPA->hdu;
        if (! hdu->header) {
            hdu->header = psMetadataAlloc();
        }
        hdu->images = spliceCells(targetCells, sourceCells, positionDependent);
        // Purge the lists
        while (psListRemove(targetCells, PS_LIST_TAIL));
        while (psListRemove(sourceCells, PS_LIST_TAIL));
    }

    return true;
}
