Changeset 6985

Timestamp:

Apr 25, 2006, 2:28:15 PM (20 years ago)

Author:

Paul Price

Message:

Adding pmChipMosaic back in after bringing it up to date with the new way of doing things.

Location:

trunk/psModules/src

Files:

• astrom/Makefile.am (modified) (2 diffs)
• astrom/pmChipMosaic.c (modified) (12 diffs)
• astrom/pmChipMosaic.h (modified) (1 diff)
• psmodules.h (modified) (1 diff)

trunk/psModules/src/astrom/Makefile.am

@@ -14 +14 @@
         pmHDUUtils.c \
         pmReadout.c \
+        pmChipMosaic.c \
         pmConcepts.c \
         pmConceptsRead.c \
@@ -39 +40 @@
         pmHDUUtils.h \
         pmReadout.h \
+        pmChipMosaic.h \
         pmConcepts.h \
         pmConceptsRead.h \

trunk/psModules/src/astrom/pmChipMosaic.c

@@ -1 +1 @@
 #include <stdio.h>
 #include <assert.h>
-
 #include "pslib.h"
 #include "pmFPA.h"
-#include "pmChipMosaic.h"
-
-#define MEM_LEAKS 0
+#include "pmHDU.h"
+#include "psRegionIsBad.h"
+
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// File-static (private) functions
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+// Do two regions overlap?
+#define REGIONS_OVERLAP(region1, region2) \
+((region1->x0 > region2->x0 && region1->x0 < region2->x1) || \
+ (region1->x1 > region2->x0 && region1->x1 < region2->x1) || \
+ (region1->y0 > region2->y0 && region1->y0 < region2->y1) || \
+ (region1->y1 > region2->y0 && region1->y1 < region2->y1))
 
 // Compare a value with a maximum and minimum
@@ -17 +27 @@
 }
 
+// Update a metadata entry directly
+#define MD_UPDATE(MD, NAME, TYPE, VALUE) \
+{ \
+    psMetadataItem *item = psMetadataLookup(MD, NAME); \
+    item->data.TYPE = VALUE; \
+}
+
+// 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?
+                          )
+{
+    if (primary) {
+        *bounds = psRegionSet(INFINITY, 0, INFINITY, 0);
+    }
+    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 || psRegionIsBad(*trimsec)) {
+            psError(PS_ERR_IO, 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)) {
+            return false;
+        }
+
+        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);
+            return false;
+        }
+        if (biassecs->n == 0) {
+            // No point allocating an iterator if there's nothing there to iterate on
+            continue;
+        }
+        psListIterator *biassecsIter = psListIteratorAlloc(biassecs, PS_LIST_HEAD, false); // Iterator
+        psRegion *biassec = NULL;       // Bias section from iteration
+        while ((biassec = psListGetAndIncrement(biassecsIter))) {
+            if (psRegionIsBad(*biassec)) {
+                continue;
+            }
+            if (REGIONS_OVERLAP(biassec, bounds)) {
+                psFree(biassecsIter);
+                return false;
+            }
+        }
+        psFree(biassecsIter);
+    }
+
+    // If we've gotten this far, everything is fine.
+    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".
+                         )
+{
+    // Check that we've got the HDU in the chip or the FPA
+    if ((!chip->hdu || !chip->hdu->images) && (!chip->parent->hdu || !chip->parent->hdu->images)) {
+        return NULL;
+    }
+
+    // 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) {
+            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)) {
+        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 < 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);
+            psFree(imageBounds);
+            return NULL;
+        }
+    }
+
+    // Need to check all the other chips if the HDU is in the FPA
+    pmFPA *fpa = chip->parent;          // The parent FPA
+    if (fpa->hdu && fpa->hdu->images) {
+        psArray *chips = fpa->chips;    // Array of chips
+        for (int i = 0; i < chips->n; i++) {
+            pmChip *testChip = chips->data[i]; // The chip of interest
+            if (testChip == chip) {
+                // Already done this one
+                continue;
+            }
+            if (!chipContiguous(imageBounds, testChip, false)) {
+                psTrace(__func__, 5, "Image isn't contiguous.\n");
+                psFree(imageBounds);
+                return NULL;
+            }
+        }
+    }
+
+    return imageBounds;
+}
+
 // Mosaic multiple images, with flips, binning and offsets
-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
-                        )
+static psImage *imageMosaic(const psArray *source, // Images to splice in
+                            const psVector *xFlip, const psVector *yFlip, // Need to flip x and y?
+                            const psVector *xBinSource, // Binning in x of source images
+                            const psVector *yBinSource, // Binning in 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
+                           )
 {
     assert(source);
@@ -33 +235 @@
     assert(x0 && x0->type.type == PS_TYPE_S32);
     assert(y0 && y0->type.type == PS_TYPE_S32);
+    assert(xFlip->n == source->n);
+    assert(yFlip->n == source->n);
+    assert(xBinSource->n == source->n);
+    assert(yBinSource->n == source->n);
+    assert(x0->n == source->n);
+    assert(y0->n == source->n);
+
+    if (source->n == 0) {
+        return NULL;
+    }
 
     // Get the maximum extent of the mosaic image
@@ -40 +252 @@
     int yMax = - INT_MAX;
     psElemType type = 0;
+    int numImages = 0;                  // Number of images
     for (int i = 0; i < source->n; i++) {
         psImage *image = source->data[i]; // The image of interest
@@ -45 +258 @@
             continue;
         }
+        numImages++;
 
         // Only implemented for F32 and U8 images so far.
@@ -66 +280 @@
         COMPARE(x0->data.S32[i] + xParity * xBinSource->data.S32[i] * image->numCols - xParity, xMin, xMax);
         COMPARE(y0->data.S32[i] + yParity * yBinSource->data.S32[i] * image->numRows - yParity, yMin, yMax);
+    }
+    if (numImages == 0) {
+        return NULL;
     }
 
@@ -137 +354 @@
 static bool cellConcepts(pmCell *target,// Target cell
                          psArray *sources, // Source cells
-                         int xBin, int yBin // Binning
+                         int xBin, int yBin, // Binning
+                         psRegion *trimsec // The trim section
                         )
 {
@@ -143 +361 @@
     float gain       = 0.0;             // Gain
     float readnoise  = 0.0;             // Read noise
-    float saturation = 0.0;             // Saturation level
-    float bad        = 0.0;             // Bad level
+    float saturation = INFINITY;        // Saturation level
+    float bad        = -INFINITY;       // Bad level
     float exposure   = 0.0;             // Exposure time
     float darktime   = 0.0;             // Dark time
@@ -160 +378 @@
         gain       += psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN");
         readnoise  += psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");
-        saturation += psMetadataLookupF32(NULL, cell->concepts, "CELL.SATURATION");
-        bad        += psMetadataLookupF32(NULL, cell->concepts, "CELL.BAD");
         exposure   += psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE");
         darktime   += psMetadataLookupF32(NULL, cell->concepts, "CELL.DARKTIME");
-        time       += psTimeToMJD(psMetadataLookupPtr(NULL, cell->concepts, "CELL.TIME"));
+        psTime *cellTime = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TIME");
+        time       += psTimeToMJD(cellTime);
         if (i == 0) {
             timeSys = psMetadataLookupS32(NULL, cell->concepts, "CELL.TIMESYS");
@@ -172 +389 @@
             success = false;
         }
-    }
-    gain       /= (float)nCells;
-    readnoise  /= (float)nCells;
-    saturation /= (float)nCells;
-    bad        /= (float)nCells;
-    exposure   /= (float)nCells;
-    darktime   /= (float)nCells;
-    psTime *timePtr = psTimeFromMJD(time/(double)nCells);
-    timePtr = psTimeConvert(timePtr, timeSys);
-
-    // XXX *REALLY* need a generic "concept update" function that handles the type and comments transparently.
-    psMetadataAddF32(target->concepts, PS_LIST_TAIL, "CELL.GAIN", PS_META_REPLACE, "Gain (e/ADU)", gain);
-    psMetadataAddF32(target->concepts, PS_LIST_TAIL, "CELL.READNOISE", PS_META_REPLACE, "Read noise (e)", readnoise);
-    psMetadataAddF32(target->concepts, PS_LIST_TAIL, "CELL.SATURATION", PS_META_REPLACE, "Saturation level (ADU)", saturation);
-    psMetadataAddF32(target->concepts, PS_LIST_TAIL, "CELL.BAD", PS_META_REPLACE, "Bad level (ADU)", bad);
-    psMetadataAddF32(target->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE, "Exposure time (sec)", exposure);
-    psMetadataAddF32(target->concepts, PS_LIST_TAIL, "CELL.DARKTIME", PS_META_REPLACE, "Time since last CCD flush (sec)", darktime);
-    psMetadataAddPtr(target->concepts, PS_LIST_TAIL, "CELL.TIME", PS_DATA_TIME | PS_META_REPLACE, "Time of observation", timePtr);
-    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.TIMESYS", PS_META_REPLACE, "Time system", timeSys);
-    psFree(timePtr);
-
-    // Now fill in the ones I know by other means
-    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.X0", PS_META_REPLACE, "Position of (0,0) on the chip", 0);
-    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.Y0", PS_META_REPLACE, "Position of (0,0) on the chip", 0);
-    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.XPARITY", PS_META_REPLACE, "Orientation in x compared to the rest of the FPA", 1);
-    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.YPARITY", PS_META_REPLACE, "Orientation in x compared to the rest of the FPA", 1);
-    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.XBIN", PS_META_REPLACE, "Binning in x", xBin);
-    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_META_REPLACE, "Binning in x", yBin);
-    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.READDIR", PS_META_REPLACE, "Read direction (faked)", 1);
-    psRegion *trimsec = psMetadataLookupPtr(NULL, target->concepts, "CELL.TRIMSEC");
-    trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0.0;
+
+        float cellSaturation = psMetadataLookupF32(NULL, cell->concepts, "CELL.SATURATION");
+        if (cellSaturation < saturation) {
+            saturation = cellSaturation;
+        }
+        float cellBad = psMetadataLookupF32(NULL, cell->concepts, "CELL.BAD");
+        if (cellBad < bad) {
+            bad = cellBad;
+        }
+    }
+    gain      /= (float)nCells;
+    readnoise /= (float)nCells;
+    exposure  /= (float)nCells;
+    darktime  /= (float)nCells;
+    time      /= (double)nCells;
+
+    MD_UPDATE(target->concepts, "CELL.GAIN", F32, gain);
+    MD_UPDATE(target->concepts, "CELL.READNOISE", F32, readnoise);
+    MD_UPDATE(target->concepts, "CELL.SATURATION", F32, saturation);
+    MD_UPDATE(target->concepts, "CELL.BAD", F32, bad);
+    MD_UPDATE(target->concepts, "CELL.EXPOSURE", F32, exposure);
+    MD_UPDATE(target->concepts, "CELL.DARKTIME", F32, darktime);
+    MD_UPDATE(target->concepts, "CELL.TIMESYS", S32, timeSys);
+    MD_UPDATE(target->concepts, "CELL.X0", S32, 0);
+    MD_UPDATE(target->concepts, "CELL.Y0", S32, 0);
+    MD_UPDATE(target->concepts, "CELL.XPARITY", S32, 1);
+    MD_UPDATE(target->concepts, "CELL.YPARITY", S32, 1);
+    MD_UPDATE(target->concepts, "CELL.XBIN", S32, xBin);
+    MD_UPDATE(target->concepts, "CELL.YBIN", S32, yBin);
+
+    // CELL.TIME needs special care
+    {
+        psMetadataItem *timeItem = psMetadataLookup(target->concepts, "CELL.TIME");
+        psFree(timeItem->data.V);
+        timeItem->data.V = psTimeFromMJD(time);
+    }
+
+    // CELL.TRIMSEC needs special care
+    {
+        psMetadataItem *trimsecItem = psMetadataLookup(target->concepts, "CELL.TRIMSEC");
+        psFree(trimsecItem->data.V);
+        trimsecItem->data.V = psMemIncrRefCounter(trimsec);
+    }
+
 
     return success;
@@ -208 +438 @@
 
 
-
-// Mosaic a chip together into a single image
-int pmChipMosaic(pmChip *chip,// Chip to mosaic
-                 int xBinChip, int yBinChip // Binning of mosaic image in x and y
-                )
+// 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
+               )
 {
-
     psArray *cells = chip->cells;       // The array of cells
-    psArray *images = psArrayAlloc(cells->n); // Array of images that will be mosaicked
-    psArray *weights = psArrayAlloc(cells->n); // Array of weight images to be mosaicked
-    psArray *masks = psArrayAlloc(cells->n); // Array of mask images to be mosaicked
-    psVector *x0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin x coordinates
-    psVector *y0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin y coordinates
-    psVector *xBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in x
-    psVector *yBin = psVectorAlloc(cells->n, PS_TYPE_S32); // Binning in y
-    psVector *xFlip = psVectorAlloc(cells->n, PS_TYPE_U8); // Flip in x?
-    psVector *yFlip = psVectorAlloc(cells->n, PS_TYPE_U8); // Flip in y?
+    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;
+
+    // Binning for the mosaicked chip is the minimum binning allowed by the cells
+    *xBinChip = INT_MAX;
+    *yBinChip = INT_MAX;
 
     // Set up the required inputs
-    psTrace(__func__, 1, "Mosaicking %d cells...\n", cells->n);
-    for (int i = 0; i < cells->n; i++) {
+    psTrace(__func__, 1, "Mosaicking %d cells...\n", numCells);
+    bool good = true;                   // Is everything good, well-behaved?
+    for (int i = 0; i < numCells && good; i++) {
         pmCell *cell = cells->data[i];  // The cell of interest
         if (!cell) {
             continue;
         }
-        x0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
-        y0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
+        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(NULL, cell->concepts, "CELL.XBIN");
-        yBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN");
-        int xParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
-        int yParity = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
-        if (xParity == 1) {
-            xFlip->data.U8[i] = 0;
-        } else if (xParity == -1) {
+        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 {
-            psLogMsg(__func__, PS_LOG_WARN, "The x parity of cell %d is not +/- 1 (it's %d) --- "
-                     "assuming +1.\n", i, xParity);
             xFlip->data.U8[i] = 0;
         }
-        if (yParity == 1) {
-            yFlip->data.U8[i] = 0;
-        } else if (yParity == -1) {
+        if (yParity == -1) {
             yFlip->data.U8[i] = 1;
         } else {
-            psLogMsg(__func__, PS_LOG_WARN, "The y parity of cell %d is not +/- 1 (it's %d) --- "
-                     "assuming +1.\n", i, yParity);
             yFlip->data.U8[i] = 0;
         }
@@ -271 +533 @@
         masks->data[i]   = psMemIncrRefCounter(readout->mask);
     }
+
     // Mosaic the images together and we're done
-    psImage *image = p_pmImageMosaic(images, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0);
-    psImage *weight = p_pmImageMosaic(weights, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0);
-    psImage *mask = p_pmImageMosaic(masks, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0);
+    if (good) {
+        *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0);
+        *mosaicWeights = imageMosaic(weights, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0);
+        *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0);
+    }
 
     // Clean up
-    psFree(x0);
-    psFree(y0);
-    psFree(xBin);
-    psFree(yBin);
-    psFree(xFlip);
-    psFree(yFlip);
     psFree(images);
     psFree(weights);
     psFree(masks);
-    int nCells = cells->n;
-
-    // Fix up the HDU
-    if (chip->parent->hdu) {
-        psLogMsg(__func__, PS_LOG_WARN, "The original format has the entire FPA in a single extension.  "
-                 "The FPA hierarchy may be invalid following the pmChipMosaic.\n");
+    psFree(xFlip);
+    psFree(yFlip);
+    psFree(xBin);
+    psFree(yBin);
+    psFree(x0);
+    psFree(y0);
+
+    return good;
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Public functions
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+// Mosaic all the cells in a chip together.
+//
+// It is desirable to do this without using psImageOverlay (or similar) if it can be at all avoided (because
+// it's really really slow in that case).  There are therefore two cases:
+//
+// 1. The HDU is at the Chip or FPA level.  This is the fast case, and only works if the HDU is "nice", by
+// which I mean:
+//
+//    - the CELL.TRIMSECs are contiguous on the HDU image
+//    - the CELL.PARITYs are identically +1
+//    - the CELL.XBIN and CELL.YBIN are all identical
+//
+// Then we can just use psImageSubset to get the "mosaicked" chip.
+//
+//
+// 2. The HDU is at the cell level, or the above requirements are not met, in which case we mosaic the cells.
+// This is the slow case.  We need to:
+//
+//    - Throw away the bias regions
+//    - Convert all cells to common parity
+//    - Mosaic the cells into an HDU image using CELL.X0 and CELL.Y0
+//    - Update CELL.TRIMSECs
+//
+// 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
+                 )
+{
+    psImage *mosaicImage   = NULL;      // The mosaic image
+    psImage *mosaicMask    = NULL;      // The mosaic mask
+    psImage *mosaicWeights = NULL;      // The mosaic weights
+
+    // Find the HDU
+    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))) {
+        // 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
+        if (!hdu || !hdu->images) {
+            hdu = chip->parent->hdu;
+        }
+        mosaicImage   = psMemIncrRefCounter(psImageSubset(hdu->images->data[0], *chipRegion));
+        if (hdu->masks) {
+            mosaicMask    = psMemIncrRefCounter(psImageSubset(hdu->masks->data[0], *chipRegion));
+        }
+        if (hdu->weights) {
+            mosaicWeights = psMemIncrRefCounter(psImageSubset(hdu->weights->data[0], *chipRegion));
+        }
     } else {
-        if (! chip->hdu) {
-            psString chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
-            chip->hdu = p_pmHDUAlloc(chipName);
-        }
-        p_pmHDU *hdu = chip->hdu;
-        psArrayElementsFree(hdu->images);
-        psArrayElementsFree(hdu->weights);
-        psArrayElementsFree(hdu->masks);
-        hdu->images  = psArrayRealloc(hdu->images,1);
-        hdu->weights = psArrayRealloc(hdu->weights, 1);
-        hdu->masks   = psArrayRealloc(hdu->masks, 1);
-        hdu->images->data[0]  = image;
-        hdu->weights->data[0] = weight;
-        hdu->masks->data[0]   = mask;
-        psMetadataAddS32(hdu->header, PS_LIST_TAIL, "NAXIS1", PS_META_REPLACE, "Number of columns", image->numCols);
-        psMetadataAddS32(hdu->header, PS_LIST_TAIL, "NAXIS2", PS_META_REPLACE, "Number of rows", image->numRows);
-    }
-
-    // Chop off all the component cells, and construct a new one
-    pmCell *newCell = pmCellAlloc(NULL, NULL, __func__); // New cell
-    cellConcepts(newCell, cells, xBinChip, yBinChip);
-    pmChipFreeCells(chip);
-    // Have to put in the new cell manually, since we didn't want to put it in before blowing the cells away.
-    newCell->parent = chip;
-    psArrayAdd(chip->cells, 1, newCell);
-    newCell->exists = true;
-    newCell->process = true;
+        // 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");
+            return false;
+        }
+        chipRegion = psRegionAlloc(NAN, NAN, NAN, NAN); // We've cut and paste, so there's no valid trimsec
+    }
+
+    // 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);
+    psFree(chipRegion);
+    chip->mosaic = (struct pmCell*)newCell;
 
     // Now make a new readout to go in the new cell
     pmReadout *newReadout = pmReadoutAlloc(newCell); // New readout
-    // Want the readouts to contain a subimage, but that subimage is the whole image.
-    // This preserves the relationship there was before, where freeing the parent frees the child.
-    psRegion entire = {0.0, 0.0, 0.0, 0.0};
-    newReadout->image = psMemIncrRefCounter(psImageSubset(image, entire));
-    newReadout->weight = psMemIncrRefCounter(psImageSubset(weight, entire));
-    newReadout->mask = psMemIncrRefCounter(psImageSubset(mask, entire));
-    // Drop references
-    psFree(newReadout);
-    psFree(newCell);
-
-    // Well, we've stuffed around with the camera configuration, so it's no longer valid...
-    #if 0
-
-    psFree(chip->parent->camera);
-    chip->parent->camera = NULL;
-    #endif
-
-    return nCells;
-}
-
-
-int pmFPAMosaicCells(pmFPA *fpa,        // FPA
-                     int xBinChip, int yBinChip // Binning of mosaic image in x and y
-                    )
-{
-    assert(fpa);
-
-    int numChips = 0;
-    psArray *chips = fpa->chips;        // Component chips
-    for (int i = 0; i < chips->n; i++) {
-        pmChip *chip = chips->data[i];  // The chip of interest
-        if (! chip || ! chip->exists || ! chip->process) {
-            continue;
-        }
-
-        if (pmChipMosaic(chip, xBinChip, yBinChip) > 0) {
-            numChips++;
-        }
-    }
-
-    return numChips;
-
-}
+    newReadout->image  = mosaicImage;
+    newReadout->mask   = mosaicMask;
+    newReadout->weight = mosaicWeights;
+    psFree(newReadout);                 // Drop reference
+
+    return true;
+}

trunk/psModules/src/astrom/pmChipMosaic.h

@@ -5 +5 @@
 #include "pmFPA.h"
 
-// Mosaic multiple images, with flips, binning and offsets
-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
-                        );
-
-// Mosaic a chip together into a single cell with single readout
-int pmChipMosaic(pmChip *chip,// Chip to mosaic
-                 int xBinChip, int yBinChip // Binning of mosaic image in x and y
-                );
-
-// Mosaic all the cells in all (valid) chips together (neglecting the overscans); return number of chips
-int pmFPAMosaicCells(pmFPA *fpa,        // FPA
-                     int xBinChip, int yBinChip // Binning of mosaic image in x and y
-                    );
+// Mosaic all cells within a chip
+bool pmChipMosaic(pmChip *chip      // Chip whose cells will be mosaicked
+                 );
 
 #endif

trunk/psModules/src/psmodules.h

@@ -37 +37 @@
 #include <pmFPAWrite.h>
 #include <pmFPA_JPEG.h>
-
 #include <pmReadout.h>
-// #include <pmChipMosaic.h>
+#include <pmChipMosaic.h>
 // #include <pmFPAAstrometry.h>