Changeset 7454

Timestamp:

Jun 9, 2006, 12:23:43 PM (20 years ago)

Author:

Paul Price

Message:

Adding FPA mosaicking; generalised some of the chip-specific functions to support multiple calls (as for an FPA with multiple chips). Still need to update the headers. Will rename these files to pmFPAMosaic.[ch] in the near future.

Location:

trunk/psModules/src/camera

Files:

• pmChipMosaic.c (modified) (23 diffs)
• pmChipMosaic.h (modified) (1 diff)

trunk/psModules/src/camera/pmChipMosaic.c

@@ -4 +4 @@
 #include "pmFPA.h"
 #include "pmHDU.h"
-
+#include "pmChipMosaic.h"
+
+#define CELL_LIST_BUFFER 50             // Buffer size for cell lists
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -33 +35 @@
 }
 
-// Are the pixels for the chip contiguous on the HDU?
-// Work this out by examining all the CELL.TRIMSEC and CELL.BIASSEC regions for the component cells
-static bool chipContiguous(psRegion *bounds, // The bounds of the image, altered if primary==true
-                           pmChip *chip, // The chip to examine for contiguity
-                           bool primary // Is this the primary chip of interest?
-                          )
-{
-    assert(bounds);
+
+// Get the bounds for an chip's pixels on the HDU
+static bool chipBounds(psRegion *bounds, // The bounds for the chip
+                       const pmChip *chip // The chip to examine for contiguity
+                      )
+{
     assert(chip);
 
-    if (primary) {
-        *bounds = psRegionSet(INFINITY, 0, INFINITY, 0);
-    }
     psArray *cells = chip->cells;       // The array of cells
     bool mdok = true;                   // Status of MD lookup
@@ -52 +49 @@
         psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section
         if (!mdok || !trimsec || psRegionIsNaN(*trimsec)) {
-            psError(PS_ERR_IO, true, "CELL.TRIMSEC hasn't been set for cell %d.\n", i);
+            psError(PS_ERR_UNKNOWN, true, "CELL.TRIMSEC hasn't been set for cell %d.\n", i);
             return false;
         }
 
-        if (primary) {
-            if (trimsec->x0 < bounds->x0) {
-                bounds->x0 = trimsec->x0;
-            }
-            if (trimsec->x1 > bounds->x1) {
-                bounds->x1 = trimsec->x1;
-            }
-            if (trimsec->y0 < bounds->y0) {
-                bounds->y0 = trimsec->y0;
-            }
-            if (trimsec->y1 > bounds->y1) {
-                bounds->y1 = trimsec->y1;
-            }
-        } else if (REGIONS_OVERLAP(trimsec, bounds)) {
+        if (trimsec->x0 < bounds->x0) {
+            bounds->x0 = trimsec->x0;
+        }
+        if (trimsec->x1 > bounds->x1) {
+            bounds->x1 = trimsec->x1;
+        }
+        if (trimsec->y0 < bounds->y0) {
+            bounds->y0 = trimsec->y0;
+        }
+        if (trimsec->y1 > bounds->y1) {
+            bounds->y1 = trimsec->y1;
+        }
+    }
+
+    return true;
+}
+
+// Make sure the TRIMSEC doesn't overlap with the established image bounds
+static bool chipContiguousTrimsec(psRegion *bounds, // The bounds of the image, altered if primary==true
+                                  const pmChip *chip // The chip to examine for contiguity
+                                 )
+{
+    assert(bounds);
+    assert(chip);
+
+    psArray *cells = chip->cells;       // The array of cells
+    bool mdok = true;                   // Status of MD lookup
+    for (int i = 0; i < cells->n; i++) {
+        pmCell *cell = cells->data[i];  // Cell of interest
+        psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section
+        if (!mdok || !trimsec || psRegionIsNaN(*trimsec)) {
+            psError(PS_ERR_UNKNOWN, true, "CELL.TRIMSEC hasn't been set for cell %d.\n", i);
             return false;
         }
 
+        if (REGIONS_OVERLAP(trimsec, bounds)) {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+// Make sure the BIASSEC doesn't overlap with the established image bounds
+static bool chipContiguousBiassec(psRegion *bounds, // The bounds of the image, altered if primary==true
+                                  const pmChip *chip // The chip to examine for contiguity
+                                 )
+{
+    assert(bounds);
+    assert(chip);
+
+    // Check that the biases don't get in the way
+    psArray *cells = chip->cells;       // The array of cells
+    bool mdok = true;                   // Status of MD lookup
+    for (int i = 0; i < cells->n; i++) {
+        pmCell *cell = cells->data[i];  // Cell of interest
         psList *biassecs = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.BIASSEC"); // Bias sections
         if (!mdok || !biassecs) {
-            psError(PS_ERR_IO, true, "CELL.BIASSEC hasn't been set for cell %d.\n", i);
+            psError(PS_ERR_UNKNOWN, true, "CELL.BIASSEC hasn't been set for cell %d.\n", i);
             return false;
         }
@@ -100 +136 @@
 }
 
+// Are the pixels for the FPA contiguous on the HDU?
+// Work this out by examining all the CELL.TRIMSEC and CELL.BIASSEC regions for the component cells
+static bool fpaContiguous(psRegion *bounds, // The bounds of the image, returned
+                          const pmFPA *fpa // The FPA to examine for contiguity
+                         )
+{
+    assert(bounds);
+    assert(fpa);
+
+    *bounds = psRegionSet(INFINITY, 0, INFINITY, 0);
+
+    // Get the size of the pixels on the HDU
+    psArray *chips = fpa->chips;        // The array of chips
+    for (int i = 0; i < chips->n; i++) {
+        pmChip *chip = chips->data[i];  // Chip of interest
+        if (!chipBounds(bounds, chip)) {
+            return false;
+        }
+    }
+
+    // Make sure the bias regions don't get in the way of the HDU
+    for (int i = 0; i < chips->n; i++) {
+        pmChip *chip = chips->data[i];  // Chip of interest
+        if (!chipContiguousBiassec(bounds, chip)) {
+            return false;
+        }
+    }
+
+    // If we got through it all, they must all be contiguous
+    return true;
+}
+
+
+
+// Check a cell for niceness in the parity and binning
+static bool niceCellParityBinning(int *xBin, int *yBin, // Binning for cell, to be returned
+                                  const pmCell *cell // Cell to check for niceness
+                                 )
+{
+    assert(xBin);
+    assert(yBin);
+    assert(cell);
+
+    // A "nice" cell must have only a single readout
+    if (cell->readouts->n != 1) {
+        return false;
+    }
+
+    // A "nice" cell must have parity == 1
+    bool mdok = true;                   // Status of MD lookup
+    int xParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XPARITY"); // Parity in x
+    if (!mdok || xParity == 0) {
+        psError(PS_ERR_UNKNOWN, true, "CELL.XPARITY hasn't been set for cell.\n");
+        return false;
+    }
+    if (xParity != 1) {
+        return false;
+    }
+    int yParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YPARITY"); // Parity in y
+    if (!mdok || yParity == 0) {
+        psError(PS_ERR_UNKNOWN, true, "CELL.YPARITY hasn't been set for cell.\n");
+        return false;
+    }
+    if (yParity != 1) {
+        return false;
+    }
+
+    // A "nice" cell must have consistent binning
+    int xBinCell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XBIN"); // Binning in x
+    if (!mdok || xBin <= 0) {
+        psError(PS_ERR_UNKNOWN, true, "CELL.XPARITY hasn't been set for cell.\n");
+        return false;
+    }
+    int yBinCell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YBIN"); // Binning in y
+    if (!mdok || yBin <= 0) {
+        psError(PS_ERR_UNKNOWN, true, "CELL.YPARITY hasn't been set for cell.\n");
+        return false;
+    }
+    if (*xBin == 0 || *yBin == 0) {
+        *xBin = xBinCell;
+        *yBin = yBinCell;
+    } else if (xBinCell != *xBin || yBinCell != *yBin) {
+        return false;
+    }
+
+    return true;
+}
+
+
+// Check a cell for niceness in the boundaries
+static bool niceCellBounds(const pmCell *cell, // Cell to check for niceness
+                           const psRegion *imageBounds // Bounds of the image on the HDU
+                          )
+{
+    // A "nice" cell must have the (0,0) pixel at CELL.X0,CELL.Y0
+    bool mdok = true;                   // Status of MD lookup
+    int x0 = psMetadataLookupS32(&mdok, cell->concepts, "CELL.X0"); // Position of (0,0) on chip
+    if (!mdok) {
+        psError(PS_ERR_UNKNOWN, true, "CELL.X0 hasn't been set for cell.\n");
+        return false;
+    }
+    int y0 = psMetadataLookupS32(&mdok, cell->concepts, "CELL.Y0"); // Position of (0,0) on chip
+    if (!mdok) {
+        psError(PS_ERR_UNKNOWN, true, "CELL.Y0 hasn't been set for cell.\n");
+        return false;
+    }
+    pmReadout *readout = cell->readouts->data[0]; // A representative readout
+    if (!readout) {
+        return false;                   // Nothing here
+    }
+    if (x0 != readout->col0 + readout->image->col0 - (int)imageBounds->x0 ||
+            y0 != readout->row0 + readout->image->row0 - (int)imageBounds->y0) {
+        psTrace(__func__, 5, "CELL.X0,Y0 don't match: %d,%d vs %d,%d\n", x0, y0,
+                readout->col0 + readout->image->col0 - (int)imageBounds->x0,
+                readout->row0 + readout->image->row0 - (int)imageBounds->y0);
+        return false;
+    }
+
+    return true;
+}
+
 
 // Is the chip "nice"?  If so, return the region containing the chip pixels
 static psRegion *niceChip(int *xBinChip, int *yBinChip, // Binning for chip, to be returned
-                          pmChip *chip  // Chip to examine for "niceness".
+                          const pmChip *chip // Chip to examine for "niceness".
                          )
 {
@@ -116 +273 @@
 
     // Check parity and binning for component cells
-    bool mdok = true;                   // Status of MD lookup
     *xBinChip = 0;
     *yBinChip = 0;
     for (int i = 0; i < chip->cells->n; i++) {
         pmCell *cell = chip->cells->data[i]; // The cell of interest
-
-        // A "nice" chip must have only a single readout
-        if (cell->readouts->n != 1) {
+        if (!niceCellParityBinning(xBinChip, yBinChip, cell)) {
             return NULL;
         }
-
-        // A "nice" chip must have parity == 1
-        int xParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XPARITY"); // Parity in x
-        if (!mdok || xParity == 0) {
-            psError(PS_ERR_IO, true, "CELL.XPARITY hasn't been set for cell %d.\n", i);
-            return NULL;
-        }
-        if (xParity != 1) {
-            return NULL;
-        }
-        int yParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YPARITY"); // Parity in y
-        if (!mdok || yParity == 0) {
-            psError(PS_ERR_IO, true, "CELL.YPARITY hasn't been set for cell %d.\n", i);
-            return NULL;
-        }
-        if (yParity != 1) {
-            return NULL;
-        }
-
-        // A "nice" chip must have consistent binning
-        int xBin = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XBIN"); // Binning in x
-        if (!mdok || xBin <= 0) {
-            psError(PS_ERR_IO, true, "CELL.XPARITY hasn't been set for cell %d.\n", i);
-            return NULL;
-        }
-        int yBin = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YBIN"); // Binning in y
-        if (!mdok || yBin <= 0) {
-            psError(PS_ERR_IO, true, "CELL.YPARITY hasn't been set for cell %d.\n", i);
-            return NULL;
-        }
-        if (*xBinChip == 0 || *yBinChip == 0) {
-            *xBinChip = xBin;
-            *yBinChip = yBin;
-        } else if (xBin != *xBinChip || yBin != *yBinChip) {
-            return NULL;
-        }
     }
 
     // Now check that the pixels are all contiguous
-    psRegion *imageBounds = psRegionAlloc(0, 0, 0, 0); // Bound of image on HDU
-    if (!chipContiguous(imageBounds, chip, true)) {
+    psRegion *imageBounds = psRegionAlloc(INFINITY, 0, INFINITY, 0); // Bound of image on HDU
+    if (!chipBounds(imageBounds, chip) || !chipContiguousBiassec(imageBounds, chip)) {
         psTrace(__func__, 5, "Image isn't contiguous.\n");
         psFree(imageBounds);
@@ -178 +296 @@
     for (int i = 0; i < chip->cells->n; i++) {
         pmCell *cell = chip->cells->data[i]; // The cell of interest
-
-        // A "nice" chip must have the (0,0) pixel at CELL.X0,CELL.Y0
-        int x0 = psMetadataLookupS32(&mdok, cell->concepts, "CELL.X0"); // Position of (0,0) on chip
-        if (!mdok) {
-            psError(PS_ERR_IO, true, "CELL.X0 hasn't been set for cell %d.\n", i);
-            return NULL;
-        }
-        int y0 = psMetadataLookupS32(&mdok, cell->concepts, "CELL.Y0"); // Position of (0,0) on chip
-        if (!mdok) {
-            psError(PS_ERR_IO, true, "CELL.Y0 hasn't been set for cell %d.\n", i);
-            return NULL;
-        }
-        pmReadout *readout = cell->readouts->data[0]; // A representative readout
-        if (!readout) {
-            return NULL;                // Nothing here
-        }
-        if (x0 != readout->col0 + readout->image->col0 - (int)imageBounds->x0 ||
-                y0 != readout->row0 + readout->image->row0 - (int)imageBounds->y0) {
-            psTrace(__func__, 5, "CELL.X0,Y0 don't match: %d,%d vs %d,%d\n", x0, y0,
-                    readout->col0 + readout->image->col0 - (int)imageBounds->x0,
-                    readout->row0 + readout->image->row0 - (int)imageBounds->y0);
+        if (!niceCellBounds(cell, imageBounds)) {
             psFree(imageBounds);
             return NULL;
@@ -214 +312 @@
                 continue;
             }
-            if (!chipContiguous(imageBounds, testChip, false)) {
+            if (!chipContiguousTrimsec(imageBounds, testChip) ||
+                    !chipContiguousBiassec(imageBounds, testChip)) {
                 psTrace(__func__, 5, "Image isn't contiguous.\n");
+                psFree(imageBounds);
+                return NULL;
+            }
+        }
+    }
+
+    return imageBounds;
+}
+
+// Is the FPA "nice"?  If so, return the region containing the FPA pixels
+static psRegion *niceFPA(int *xBinFPA, int *yBinFPA, // Binning for FPA, to be returned
+                         const pmFPA *fpa  // FPA to examine for "niceness".
+                        )
+{
+    assert(xBinFPA);
+    assert(yBinFPA);
+    assert(fpa);
+
+    // Check that we've got the HDU in the chip or the FPA
+    if (!fpa->hdu || !fpa->hdu->images) {
+        return NULL;
+    }
+
+    // Check parity and binning for component cells
+    *xBinFPA = 0;
+    *yBinFPA = 0;
+    for (int i = 0; i < fpa->chips->n; i++) {
+        pmChip *chip = fpa->chips->data[i]; // The chip of interest
+        for (int j = 0; i < chip->cells->n; i++) {
+            pmCell *cell = chip->cells->data[j]; // The cell of interest
+            if (!niceCellParityBinning(xBinFPA, yBinFPA, cell)) {
+                return NULL;
+            }
+        }
+    }
+
+    // Now check that the pixels are all contiguous
+    psRegion *imageBounds = psRegionAlloc(0, 0, 0, 0); // Bound of image on HDU
+    if (!fpaContiguous(imageBounds, fpa)) {
+        psTrace(__func__, 5, "Image isn't contiguous.\n");
+        psFree(imageBounds);
+        return NULL;
+    }
+
+    psString region = psRegionToString(*imageBounds);
+    psTrace(__func__, 7, "Image bounds: %s\n", region);
+    psFree(region);
+
+    for (int i = 0; i < fpa->chips->n; i++) {
+        pmChip *chip = fpa->chips->data[i]; // The chip of interest
+        for (int j = 0; i < chip->cells->n; i++) {
+            pmCell *cell = chip->cells->data[j]; // The cell of interest
+            if (!niceCellBounds(cell, imageBounds)) {
                 psFree(imageBounds);
                 return NULL;
@@ -259 +411 @@
     psElemType type = 0;
     int numImages = 0;                  // Number of images
+    psTrace(__func__, 3, "Mosaicking %ld cells.\n", source->n);
     for (int i = 0; i < source->n; i++) {
         psImage *image = source->data[i]; // The image of interest
@@ -315 +468 @@
                 xFlip->data.U8[i] == 0 && yFlip->data.U8[i] == 0) {
             // Let someone else do the hard work; useful to test psImageOverlaySection if no other reason
-            psImageOverlaySection(mosaic, image, x0->data.S32[i], y0->data.S32[i], "+");
+            psImageOverlaySection(mosaic, image, x0->data.S32[i] - xMin, y0->data.S32[i] - yMin, "+");
         } else {
             // We have to do the hard work ourself
             for (int y = 0; y < image->numRows; y++) {
                 int yParity = yFlip->data.U8[i] ? -1 : 1;
-                float yTargetBase = (y0->data.S32[i] + yParity * yBinSource->data.S32[i] * y) / yBinTarget;
+                float yTargetBase = (y0->data.S32[i] + yParity * yBinSource->data.S32[i] * y - yMin) /
+                                    yBinTarget;
                 for (int x = 0; x < image->numCols; x++) {
                     int xParity = xFlip->data.U8[i] ? -1 : 1;
-                    float xTargetBase = (x0->data.S32[i] + xParity * xBinSource->data.S32[i] * x) /
+                    float xTargetBase = (x0->data.S32[i] + xParity * xBinSource->data.S32[i] * x - xMin) /
                                         xBinTarget;
 
@@ -359 +513 @@
 // Set the concepts in the new cell, based on the values in the old one
 static bool cellConcepts(pmCell *target,// Target cell
-                         psArray *sources, // Source cells
+                         psList *sources, // Source cells
                          int xBin, int yBin, // Binning
                          psRegion *trimsec // The trim section
@@ -378 +532 @@
     double time      = 0.0;             // Time of observation
     psTimeType timeSys = 0;             // Time system
+    int readdir      = 0;               // Cell read direction
 
     int nCells = 0;                     // Number of cells;
-    for (int i = 0; i < sources->n; i++) {
-        pmCell *cell = sources->data[i];// The cell of interest
+    psListIterator *sourcesIter = psListIteratorAlloc(sources, PS_LIST_HEAD, false); // Iterator for sources
+    pmCell *cell = NULL;                // Source cell from iteration
+    while ((cell = psListGetAndIncrement(sourcesIter))) {
         if (!cell) {
             continue;
@@ -393 +549 @@
         psTime *cellTime = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TIME");
         time       += psTimeToMJD(cellTime);
-        if (i == 0) {
+        if (nCells == 1) {
             timeSys = psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS");
-        } else if (timeSys != psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS")) {
-            psLogMsg(__func__, PS_LOG_ERROR, "Differing time systems in use: %d vs %d\n", timeSys,
-                     psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS"));
-            success = false;
+            readdir = psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR");
+        } else {
+            if (timeSys != psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS")) {
+                psLogMsg(__func__, PS_LOG_ERROR, "Differing time systems in use: %d vs %d\n", timeSys,
+                         psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS"));
+                success = false;
+            }
+            if (readdir != psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR")) {
+                psLogMsg(__func__, PS_LOG_ERROR, "Differing cell read directions in use: %d vs %d\n", readdir,
+                         psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR"));
+                success = false;
+            }
         }
 
@@ -410 +574 @@
         }
     }
+    psFree(sourcesIter);
+
     gain      /= (float)nCells;
     readnoise /= (float)nCells;
@@ -429 +595 @@
     MD_UPDATE(target->concepts, "CELL.XBIN", S32, xBin);
     MD_UPDATE(target->concepts, "CELL.YBIN", S32, yBin);
+    MD_UPDATE(target->concepts, "CELL.READDIR", S32, readdir);
 
     // CELL.TIME needs special care
@@ -448 +615 @@
 
 
+// Add a cell and its various properties to the arrays
+static bool addCell(psArray *images,    // Array of images
+                    psArray *masks,     // Array of masks
+                    psArray *weights,   // Array of weights
+                    psVector *x0,       // Array of X0
+                    psVector *y0,       // Array of Y0
+                    psVector *xBin,     // Array of XBIN
+                    psVector *yBin,     // Array of YBIN
+                    psVector *xFlip,    // Array indicating whether x axis should be flipped
+                    psVector *yFlip,    // Array indicating whether y axis should be flipped
+                    const pmCell *cell, // Cell to add
+                    int *xBinMin,       // The minimum x binning, returned
+                    int *yBinMin,       // The minimum y binning, returned
+                    bool chipStuff      // Worry about chip stuff as well?
+                   )
+{
+    if (!cell) {
+        return false;
+    }
+
+    // Expand the arrays and vectors to handle new data
+    long index = images->n;               // The index to use
+    if (images->n == images->nalloc) {
+        images  = psArrayRealloc(images,  index + CELL_LIST_BUFFER);
+        masks   = psArrayRealloc(masks,   index + CELL_LIST_BUFFER);
+        weights = psArrayRealloc(weights, index + CELL_LIST_BUFFER);
+        x0    = psVectorRealloc(x0,    index+ CELL_LIST_BUFFER);
+        y0    = psVectorRealloc(y0,    index+ CELL_LIST_BUFFER);
+        xBin  = psVectorRealloc(xBin,  index+ CELL_LIST_BUFFER);
+        yBin  = psVectorRealloc(yBin,  index+ CELL_LIST_BUFFER);
+        xFlip = psVectorRealloc(xFlip, index+ CELL_LIST_BUFFER);
+        yFlip = psVectorRealloc(yFlip, index+ CELL_LIST_BUFFER);
+    }
+
+    images->n = index + 1;
+    masks->n = index + 1;
+    weights->n = index + 1;
+    x0->n = index + 1;
+    y0->n = index + 1;
+    xBin->n = index + 1;
+    yBin->n = index + 1;
+    xFlip->n = index + 1;
+    yFlip->n = index + 1;
+
+    bool mdok = true;                   // Status of MD lookup
+    bool good = true;                   // Is everything good?
+
+    // Offset of the cell on the chip
+    int x0Cell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.X0");
+    if (!mdok) {
+        psError(PS_ERR_UNKNOWN, true, "CELL.X0 for cell is not set.\n");
+        good = false;
+    }
+    int y0Cell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.Y0");
+    if (!mdok) {
+        psError(PS_ERR_UNKNOWN, true, "CELL.Y0 for cell is not set.\n");
+        good = false;
+    }
+    psTrace(__func__, 5, "Cell %d: x0=%d y0=%d\n", index, x0Cell, y0Cell);
+
+    // Offset of the chip on the FPA
+    int x0Chip = 0, y0Chip = 0;
+    if (chipStuff) {
+        pmChip *chip = cell->parent;    // The parent chip
+        if (!chip) {
+            psError(PS_ERR_UNKNOWN, true, "Cell has no parent chip --- can't find CHIP.X0 and CHIP.Y0\n");
+            good = false;
+        }
+        x0Chip = psMetadataLookupS32(&mdok, chip->concepts, "CHIP.X0");
+        if (!mdok) {
+            psError(PS_ERR_UNKNOWN, true, "CHIP.X0 for chip is not set.\n");
+            good = false;
+        }
+        y0Chip = psMetadataLookupS32(&mdok, chip->concepts, "CHIP.Y0");
+        if (!mdok) {
+            psError(PS_ERR_UNKNOWN, true, "CHIP.Y0 for chip is not set.\n");
+            good = false;
+        }
+    }
+
+    // Binning
+    xBin->data.S32[index] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XBIN");
+    if (!mdok || xBin->data.S32[index] == 0) {
+        psError(PS_ERR_UNKNOWN, true, "CELL.XBIN for cell is not set.\n");
+        good = false;
+    }
+    if (xBin->data.S32[index] < *xBinMin) {
+        *xBinMin = xBin->data.S32[index];
+    }
+    yBin->data.S32[index] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YBIN");
+    if (!mdok || yBin->data.S32[index] == 0) {
+        psError(PS_ERR_UNKNOWN, true, "CELL.YBIN for cell is not set.\n");
+        good = false;
+    }
+    if (yBin->data.S32[index] < *yBinMin) {
+        *yBinMin = yBin->data.S32[index];
+    }
+
+    // Do we need to flip?
+    int xParityCell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XPARITY");
+    if (!mdok || (xParityCell != 1 && xParityCell != -1)) {
+        psError(PS_ERR_UNKNOWN, true, "CELL.XPARITY for cell is not set.\n");
+        good = false;
+    }
+    int yParityCell = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YPARITY");
+    if (!mdok || (yParityCell != 1 && yParityCell != -1)) {
+        psError(PS_ERR_UNKNOWN, true, "CELL.YPARITY for cell is not set.\n");
+        good = false;
+    }
+
+    // Parity of the chip on the FPA
+    int xParityChip = 0, yParityChip = 0;
+    if (chipStuff) {
+        pmChip *chip = cell->parent;    // The parent chip
+        xParityChip = psMetadataLookupS32(&mdok, chip->concepts, "CHIP.XPARITY");
+        if (!mdok || (xParityChip != 1 && xParityChip != -1)) {
+            psError(PS_ERR_UNKNOWN, true, "CHIP.XPARITY for chip is not set.\n");
+            good = false;
+        }
+        yParityChip = psMetadataLookupS32(&mdok, chip->concepts, "CHIP.YPARITY");
+        if (!mdok || (yParityChip != 1 && yParityChip != -1)) {
+            psError(PS_ERR_UNKNOWN, true, "CHIP.YPARITY for chip is not set.\n");
+            good = false;
+        }
+    }
+
+    // Set the flips on the basis of the parity
+    if (xParityCell * xParityChip == -1) {
+        xFlip->data.U8[index] = 1;
+    } else {
+        xFlip->data.U8[index] = 0;
+    }
+    if (yParityCell * yParityChip == -1) {
+        yFlip->data.U8[index] = 1;
+    } else {
+        yFlip->data.U8[index] = 0;
+    }
+
+    x0->data.S32[index] = x0Chip + x0Cell;
+    y0->data.S32[index] = y0Chip + y0Cell;
+
+    // Add the readout to the array of images to be mosaicked
+    psArray *readouts = cell->readouts; // The array of readouts
+    if (readouts->n > 1) {
+        psLogMsg(__func__, PS_LOG_WARN, "Cell contains more than one readout (%ld) --- only the first will "
+                 "be mosaicked.\n", readouts->n);
+    }
+    pmReadout *readout = readouts->data[0]; // The only readout we'll bother with
+
+    // The images to put into the mosaic
+    images->data[index]  = psMemIncrRefCounter(readout->image);
+    weights->data[index] = psMemIncrRefCounter(readout->weight);
+    masks->data[index]   = psMemIncrRefCounter(readout->mask);
+
+    psTrace(__func__, 9, "Added cell (%x) %ld: %d,%d; %d,%d, %d,%d.\n", cell, index,
+            x0->data.S32[index], y0->data.S32[index], xBin->data.S32[index], yBin->data.S32[index],
+            xFlip->data.U8[index], yFlip->data.U8[index]);
+
+    return true;
+}
+
+
 // Mosaic together the cells in a chip
-bool chipMosaic(psImage **mosaicImage,  // The mosaic image, to be returned
-                psImage **mosaicMask,   // The mosaic mask, to be returned
-                psImage **mosaicWeights, // The mosaic weights, to be returned
-                int *xBinChip, int *yBinChip, // The binning in x and y, to be returned
-                pmChip *chip            // Chip to mosaic
-               )
+static bool chipMosaic(psImage **mosaicImage, // The mosaic image, to be returned
+                       psImage **mosaicMask, // The mosaic mask, to be returned
+                       psImage **mosaicWeights, // The mosaic weights, to be returned
+                       int *xBinChip, int *yBinChip, // The binning in x and y, to be returned
+                       const pmChip *chip // Chip to mosaic
+                      )
 {
     assert(mosaicImage);
@@ -461 +790 @@
     assert(chip);
 
-    psArray *cells = chip->cells;       // The array of cells
-    int numCells = cells->n;            // Number of cells
-    psArray *images = psArrayAlloc(numCells); // Array of images that will be mosaicked
-    psArray *weights = psArrayAlloc(numCells); // Array of weight images to be mosaicked
-    psArray *masks = psArrayAlloc(numCells); // Array of mask images to be mosaicked
-    psVector *x0 = psVectorAlloc(numCells, PS_TYPE_S32); // Origin x coordinates
-    psVector *y0 = psVectorAlloc(numCells, PS_TYPE_S32); // Origin y coordinates
-    psVector *xBin = psVectorAlloc(numCells, PS_TYPE_S32); // Binning in x
-    psVector *yBin = psVectorAlloc(numCells, PS_TYPE_S32); // Binning in y
-    psVector *xFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Flip in x?
-    psVector *yFlip = psVectorAlloc(numCells, PS_TYPE_U8); // Flip in y?
-    images->n = weights->n = masks->n = numCells;
-    x0->n = y0->n = xBin->n = yBin->n = xFlip->n = yFlip->n = numCells;
+    psArray *images = psArrayAlloc(0); // Array of images that will be mosaicked
+    psArray *weights = psArrayAlloc(0); // Array of weight images to be mosaicked
+    psArray *masks = psArrayAlloc(0); // Array of mask images to be mosaicked
+    psVector *x0 = psVectorAlloc(0, PS_TYPE_S32); // Origin x coordinates
+    psVector *y0 = psVectorAlloc(0, PS_TYPE_S32); // Origin y coordinates
+    psVector *xBin = psVectorAlloc(0, PS_TYPE_S32); // Binning in x
+    psVector *yBin = psVectorAlloc(0, PS_TYPE_S32); // Binning in y
+    psVector *xFlip = psVectorAlloc(0, PS_TYPE_U8); // Flip in x?
+    psVector *yFlip = psVectorAlloc(0, PS_TYPE_U8); // Flip in y?
 
     // Binning for the mosaicked chip is the minimum binning allowed by the cells
@@ -480 +805 @@
 
     // Set up the required inputs
-    psTrace(__func__, 1, "Mosaicking %d cells...\n", numCells);
-    bool good = true;                   // Is everything good, well-behaved?
-    for (int i = 0; i < numCells && good; i++) {
+    bool allGood = true;                // Is everything good, well-behaved?
+    psArray *cells = chip->cells;       // The array of cells
+    for (int i = 0; i < cells->n; i++) {
         pmCell *cell = cells->data[i];  // The cell of interest
         if (!cell) {
             continue;
         }
-        bool mdok = true;               // Status of MD lookup
-        x0->data.S32[i] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.X0");
-        if (!mdok) {
-            psError(PS_ERR_IO, true, "CELL.X0 for cell %d is not set.\n", i);
-            good = false;
-        }
-        y0->data.S32[i] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.Y0");
-        if (!mdok) {
-            psError(PS_ERR_IO, true, "CELL.Y0 for cell %d is not set.\n", i);
-            good = false;
-        }
-        psTrace(__func__, 5, "Cell %d: x0=%d y0=%d\n", i, x0->data.S32[i], y0->data.S32[i]);
-        xBin->data.S32[i] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XBIN");
-        if (!mdok || xBin->data.S32[i] == 0) {
-            psError(PS_ERR_IO, true, "CELL.XBIN for cell %d is not set.\n", i);
-            good = false;
-        }
-        if (xBin->data.S32[i] < *xBinChip) {
-            *xBinChip = xBin->data.S32[i];
-        }
-        yBin->data.S32[i] = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YBIN");
-        if (!mdok || yBin->data.S32[i] == 0) {
-            psError(PS_ERR_IO, true, "CELL.YBIN for cell %d is not set.\n", i);
-            good = false;
-        }
-        if (yBin->data.S32[i] < *yBinChip) {
-            *yBinChip = yBin->data.S32[i];
-        }
-        int xParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.XPARITY");
-        if (!mdok || (xParity != 1 && xParity != -1)) {
-            psError(PS_ERR_IO, true, "CELL.XPARITY for cell %d is not set.\n", i);
-            good = false;
-        }
-        int yParity = psMetadataLookupS32(&mdok, cell->concepts, "CELL.YPARITY");
-        if (!mdok || (yParity != 1 && yParity != -1)) {
-            psError(PS_ERR_IO, true, "CELL.YPARITY for cell %d is not set.\n", i);
-            good = false;
-        }
-        if (xParity == -1) {
-            xFlip->data.U8[i] = 1;
-        } else {
-            xFlip->data.U8[i] = 0;
-        }
-        if (yParity == -1) {
-            yFlip->data.U8[i] = 1;
-        } else {
-            yFlip->data.U8[i] = 0;
-        }
-
-        psArray *readouts = cell->readouts; // The array of readouts
-        if (readouts->n > 1) {
-            psLogMsg(__func__, PS_LOG_WARN, "Cell %d contains more than one readout --- only the first will "
-                     "be mosaicked.\n", i);
-        }
-        pmReadout *readout = readouts->data[0]; // The only readout we'll bother with
-
-        // The images to put into the mosaic
-        images->data[i]  = psMemIncrRefCounter(readout->image);
-        weights->data[i] = psMemIncrRefCounter(readout->weight);
-        masks->data[i]   = psMemIncrRefCounter(readout->mask);
+        allGood |= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip,
+                           cell, xBinChip, yBinChip, false);
     }
 
     // Mosaic the images together and we're done
-    if (good) {
+    if (allGood) {
         *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0);
         *mosaicWeights = imageMosaic(weights, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0);
@@ -567 +834 @@
     psFree(y0);
 
-    return good;
+    return allGood;
+}
+
+// Mosaic together the cells in a FPA
+static bool fpaMosaic(psImage **mosaicImage, // The mosaic image, to be returned
+                      psImage **mosaicMask, // The mosaic mask, to be returned
+                      psImage **mosaicWeights, // The mosaic weights, to be returned
+                      int *xBinFPA, int *yBinFPA, // The binning in x and y, to be returned
+                      const pmFPA *fpa  // FPA to mosaic
+                     )
+{
+    assert(mosaicImage);
+    assert(mosaicMask);
+    assert(mosaicWeights);
+    assert(fpa);
+
+    psArray *images = psArrayAlloc(0); // Array of images that will be mosaicked
+    psArray *weights = psArrayAlloc(0); // Array of weight images to be mosaicked
+    psArray *masks = psArrayAlloc(0); // Array of mask images to be mosaicked
+    psVector *x0 = psVectorAlloc(0, PS_TYPE_S32); // Origin x coordinates
+    psVector *y0 = psVectorAlloc(0, PS_TYPE_S32); // Origin y coordinates
+    psVector *xBin = psVectorAlloc(0, PS_TYPE_S32); // Binning in x
+    psVector *yBin = psVectorAlloc(0, PS_TYPE_S32); // Binning in y
+    psVector *xFlip = psVectorAlloc(0, PS_TYPE_U8); // Flip in x?
+    psVector *yFlip = psVectorAlloc(0, PS_TYPE_U8); // Flip in y?
+
+    // Binning for the mosaicked chip is the minimum binning allowed by the cells
+    *xBinFPA = INT_MAX;
+    *yBinFPA = INT_MAX;
+
+    // Set up the required inputs
+    bool allGood = true;                // Is everything good, well-behaved?
+    psArray *chips = fpa->chips;        // Array of chips
+    for (int i = 0; i < chips->n; i++) {
+        pmChip *chip = chips->data[i];  // The chip of interest
+        if (!chip) {
+            continue;
+        }
+        psArray *cells = chip->cells;   // The array of cells
+        for (int j = 0; j < cells->n; j++) {
+            pmCell *cell = cells->data[j];  // The cell of interest
+            if (!cell) {
+                continue;
+            }
+            allGood |= addCell(images, masks, weights, x0, y0, xBin, yBin, xFlip, yFlip,
+                               cell, xBinFPA, yBinFPA, true);
+        }
+    }
+
+    // Mosaic the images together and we're done
+    if (allGood) {
+        *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0);
+        *mosaicWeights = imageMosaic(weights, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0);
+        *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0);
+    }
+
+    // Clean up
+    psFree(images);
+    psFree(weights);
+    psFree(masks);
+    psFree(xFlip);
+    psFree(yFlip);
+    psFree(xBin);
+    psFree(yBin);
+    psFree(x0);
+    psFree(y0);
+
+    return allGood;
 }
 
@@ -599 +933 @@
 // Once the demands of case 1 have been met, or case 2 has been performed, then we can create a cell to hold
 // the mosaic image.
-bool pmChipMosaic(pmChip *chip      // Chip whose cells will be mosaicked
+bool pmChipMosaic(pmChip *target,       // Target chip --- may contain only a single cell
+                  const pmChip *source  // Source chip whose cells will be mosaicked
                  )
 {
-    PS_ASSERT_PTR_NON_NULL(chip, false);
+    // Target exists, and has only a single cell
+    PS_ASSERT_PTR_NON_NULL(target, false);
+    PS_ASSERT_PTR_NON_NULL(target->cells, false);
+    if (target->cells->n != 1) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Target chip for mosaicking must contain a single cell.\n");
+        return false;
+    }
+    pmCell *targetCell = target->cells->data[0]; // The target cell
+    PS_ASSERT_PTR_NON_NULL(targetCell, false);
+    // Source exists
+    PS_ASSERT_PTR_NON_NULL(source, false);
+
 
     psImage *mosaicImage   = NULL;      // The mosaic image
@@ -611 +957 @@
     psRegion *chipRegion = NULL;        // Region on the HDU that corresponds to the chip
     int xBin = 0, yBin = 0;             // Binning for the chip mosaic
-    if ((chipRegion = niceChip(&xBin, &yBin, chip))) {
+    if ((chipRegion = niceChip(&xBin, &yBin, source))) {
         // Case 1 --- we need only cut out the region
         psTrace(__func__, 1, "Case 1 mosaicking: simple cut-out.\n");
-        pmHDU *hdu = chip->hdu;         // The HDU that has the pixels
+        pmHDU *hdu = source->hdu;       // The HDU that has the pixels
         if (!hdu || !hdu->images) {
-            hdu = chip->parent->hdu;
-        }
-        mosaicImage   = psMemIncrRefCounter(psImageSubset(hdu->images->data[0], *chipRegion));
+            hdu = source->parent->hdu;
+        }
+        mosaicImage = psMemIncrRefCounter(psImageSubset(hdu->images->data[0], *chipRegion));
         if (hdu->masks) {
-            mosaicMask    = psMemIncrRefCounter(psImageSubset(hdu->masks->data[0], *chipRegion));
+            mosaicMask = psMemIncrRefCounter(psImageSubset(hdu->masks->data[0], *chipRegion));
         }
         if (hdu->weights) {
@@ -628 +974 @@
         // Case 2 --- we need to mosaic by cut and paste
         psTrace(__func__, 1, "Case 2 mosaicking: cut and paste.\n");
-        if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicWeights, &xBin, &yBin, chip)) {
-            psError(PS_ERR_IO, false, "Unable to mosaic cells.\n");
+        if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicWeights, &xBin, &yBin, source)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to mosaic cells.\n");
             return false;
         }
@@ -635 +981 @@
     }
 
-    // Construct a new cell, set the concepts, and add the mosaic in
-    pmCell *newCell = pmCellAlloc(NULL, "MOSAIC"); // New cell
-    cellConcepts(newCell, chip->cells, xBin, yBin, chipRegion);
+    // Set the concepts for the target cell
+    psList *sourceCells = psArrayToList(source->cells); // List of cells
+    cellConcepts(targetCell, sourceCells, xBin, yBin, chipRegion);
+    psFree(sourceCells);
     psFree(chipRegion);
-    chip->mosaic = (struct pmCell*)newCell;
-
-    // Now make a new readout to go in the new cell
-    pmReadout *newReadout = pmReadoutAlloc(newCell); // New readout
+
+    // Now make a new readout to go in the target cell
+    pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout
     newReadout->image  = mosaicImage;
     newReadout->mask   = mosaicMask;
@@ -650 +996 @@
     return true;
 }
+
+
+// Same deal with the FPA
+bool pmFPAMosaic(pmFPA *target,         // Target FPA --- may contain only a single chip with a single cell
+                 const pmFPA *source    // FPA whose chips and cells will be mosaicked
+                )
+{
+    // Target exists, and has only a single chip with single cell
+    PS_ASSERT_PTR_NON_NULL(target, false);
+    PS_ASSERT_PTR_NON_NULL(target->chips, false);
+    if (target->chips->n != 1) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Target FPA for mosaicking must contain a single chip.\n");
+        return false;
+    }
+    pmChip *targetChip = target->chips->data[0]; // The target chip
+    PS_ASSERT_PTR_NON_NULL(targetChip, false);
+    PS_ASSERT_PTR_NON_NULL(targetChip->cells, false);
+    if (target->chips->n != 1) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Target FPA for mosaicking must contain a single cell.\n");
+        return false;
+    }
+    pmCell *targetCell = targetChip->cells->data[0]; // The target cell
+    PS_ASSERT_PTR_NON_NULL(targetCell, false);
+    // Source exists
+    PS_ASSERT_PTR_NON_NULL(source, false);
+
+    psImage *mosaicImage   = NULL;      // The mosaic image
+    psImage *mosaicMask    = NULL;      // The mosaic mask
+    psImage *mosaicWeights = NULL;      // The mosaic weights
+
+    // Find the HDU
+    psRegion *fpaRegion = NULL;         // Region on the HDU that corresponds to the FPA
+    int xBin = 0, yBin = 0;             // Binning for the FPA mosaic
+    if ((fpaRegion = niceFPA(&xBin, &yBin, source))) {
+        // Case 1 --- we need only cut out the region
+        psTrace(__func__, 1, "Case 1 mosaicking: simple cut-out.\n");
+        pmHDU *hdu = source->hdu;         // The HDU that has the pixels
+        mosaicImage = psMemIncrRefCounter(psImageSubset(hdu->images->data[0], *fpaRegion));
+        if (hdu->masks) {
+            mosaicMask = psMemIncrRefCounter(psImageSubset(hdu->masks->data[0], *fpaRegion));
+        }
+        if (hdu->weights) {
+            mosaicWeights = psMemIncrRefCounter(psImageSubset(hdu->weights->data[0], *fpaRegion));
+        }
+    } else {
+        // Case 2 --- we need to mosaic by cut and paste
+        psTrace(__func__, 1, "Case 2 mosaicking: cut and paste.\n");
+        if (!fpaMosaic(&mosaicImage, &mosaicMask, &mosaicWeights, &xBin, &yBin, source)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to mosaic chips.\n");
+            return false;
+        }
+        fpaRegion = psRegionAlloc(NAN, NAN, NAN, NAN); // We've cut and paste, so there's no valid trimsec
+    }
+
+    // Set the concepts for the target cell, and add the mosaic in
+    // First we need a list of cells
+    psList *sourceCells = psListAlloc(NULL); // List of source cells
+    psArray *chips = source->chips;        // Array of chips
+    for (long i = 0; i < chips->n; i++) {
+        pmChip *chip = chips->data[i];  // Chip of interest
+        if (!chip) {
+            continue;
+        }
+        psArray *cells = chip->cells;
+        for (long j = 0; j < cells->n; j++) {
+            pmCell *cell = cells->data[j]; // Cell of interest
+            if (!cell) {
+                continue;
+            }
+            psListAdd(sourceCells, PS_LIST_TAIL, cell);
+        }
+    }
+    cellConcepts(targetCell, sourceCells, xBin, yBin, fpaRegion);
+    psFree(sourceCells);
+    psFree(fpaRegion);
+
+    // Now make a new readout to go in the new cell
+    pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout
+    newReadout->image  = mosaicImage;
+    newReadout->mask   = mosaicMask;
+    newReadout->weight = mosaicWeights;
+    psFree(newReadout);                 // Drop reference
+
+    return true;
+}
+

trunk/psModules/src/camera/pmChipMosaic.h

@@ -2 +2 @@
 #define PM_CHIP_MOSAIC_H
 
-#include "pslib.h"
 #include "pmFPA.h"
 
 // Mosaic all cells within a chip
-bool pmChipMosaic(pmChip *chip      // Chip whose cells will be mosaicked
+bool pmChipMosaic(pmChip *target,       // Target chip --- may contain only a single cell
+                  const pmChip *source  // Source chip whose cells will be mosaicked
                  );
 
+// Mosaic all cells within an FPA
+bool pmFPAMosaic(pmFPA *target,         // Target FPA --- may contain only a single chip with a single cell
+                 const pmFPA *source    // FPA whose chips and cells will be mosaicked
+                );
+
 #endif