Changeset 5786 for trunk/archive/scripts/src/phase2/pmChipMosaic.c

Timestamp:

Dec 13, 2005, 5:47:35 PM (21 years ago)

Author:

Paul Price

Message:

Upgrading for psLib-0.9.0

File:

• trunk/archive/scripts/src/phase2/pmChipMosaic.c (modified) (3 diffs)

trunk/archive/scripts/src/phase2/pmChipMosaic.c

@@ -1 +1 @@
 #include <stdio.h>
+#include <assert.h>
 
 #include "pslib.h"
@@ -14 +15 @@
     }
 
-psImage *pmChipMosaic(pmChip *chip,     // Chip to mosaic
-                      int xBinChip, int yBinChip // Binning of mosaic image in x and y
+// 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
     )
 {
-    psArray *cells = chip->cells;       // The array of cells
-    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 *xParity = psVectorAlloc(cells->n, PS_TYPE_S32); // Parity in x
-    psVector *yParity = psVectorAlloc(cells->n, PS_TYPE_S32); // Parity in y
-    psArray *trimsec = psArrayAlloc(cells->n); // Trim section
+    assert(source);
+    assert(xFlip && xFlip->type.type == PS_TYPE_U8);
+    assert(yFlip && yFlip->type.type == PS_TYPE_U8);
+    assert(xBinSource && xBinSource->type.type == PS_TYPE_S32);
+    assert(yBinSource && yBinSource->type.type == PS_TYPE_S32);
+    assert(x0 && x0->type.type == PS_TYPE_S32);
+    assert(y0 && y0->type.type == PS_TYPE_S32);
 
     // Get the maximum extent of the mosaic image
@@ -32 +37 @@
     int yMin = INT_MAX;
     int yMax = 0;
-    for (int i = 0; i < cells->n; i++) {
-        pmCell *cell = cells->data[i];  // The cell of interest
-        x0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
-        y0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
-        xBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN");
-        yBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN");
-        xParity->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
-        yParity->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
-        trimsec->data[i] = psMemIncrRefCounter(psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"));
-        psTrace(__func__, 7, "Cell %d trimsec: [%.0f:%.0f,%.0f:%.0f]\n", i, ((psRegion*)trimsec->data[i])->x0,
-                ((psRegion*)trimsec->data[i])->x1, ((psRegion*)trimsec->data[i])->y0,
-                ((psRegion*)trimsec->data[i])->y1);
-        // Size of cell in x and y
-        int nx = (int)(((psRegion*)trimsec->data[i])->x1 - (int)((psRegion*)trimsec->data[i])->x0);
-        int ny = (int)(((psRegion*)trimsec->data[i])->y1 - (int)((psRegion*)trimsec->data[i])->y0);
-        psTrace(__func__, 5, "Extent of cell %d: %d -> %d , %d -> %d\n", i, x0->data.S32[i],
-                x0->data.S32[i] + xParity->data.S32[i] * xBin->data.S32[i] * nx, y0->data.S32[i],
-                y0->data.S32[i] + yParity->data.S32[i] * yBin->data.S32[i] * ny);
+    for (int i = 0; i < source->n; i++) {
+        psImage *image = source->data[i]; // The image of interest
 
-        COMPARE(x0->data.S32[i], xMin, xMax);
-        COMPARE(y0->data.S32[i], yMin, yMax);
-        // Subtract the parity to get the inclusive limit (not exclusive)
-        COMPARE(x0->data.S32[i] + xParity->data.S32[i] * xBin->data.S32[i] * nx - xParity->data.S32[i], xMin, xMax);
-        COMPARE(y0->data.S32[i] + yParity->data.S32[i] * yBin->data.S32[i] * ny - yParity->data.S32[i], yMin, yMax);
+        assert(image->type.type == PS_TYPE_F32); // Only implemented for F32 images so far.
+
+        // Size of cell in x and y
+        int xParity = xFlip->data.U8[i] ? -1 : 1;
+        int yParity = yFlip->data.U8[i] ? -1 : 1;
+        psTrace(__func__, 5, "Extent of cell %d: %d -> %d , %d -> %d\n", i, x0->data.S32[i],
+                x0->data.S32[i] + xParity * xBinSource->data.S32[i] * image->numCols, y0->data.S32[i],
+                y0->data.S32[i] + yParity * yBinSource->data.S32[i] * image->numRows);
+
+        COMPARE(x0->data.S32[i], xMin, xMax);
+        COMPARE(y0->data.S32[i], yMin, yMax);
+        // Subtract the parity to get the inclusive limit (not exclusive)
+        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);
     }
 
     // Set up the image
     // Since both upper and lower values are inclusive, we need to add one to the size
-    float xSize = (float)(xMax - xMin + 1) / (float)xBinChip;
+    float xSize = (float)(xMax - xMin + 1) / (float)xBinTarget;
     if (xSize - (int)xSize > 0) {
-        xSize += 1;
+        xSize += 1;
     }
-    float ySize = (float)(yMax - yMin + 1) / (float)yBinChip;
+    float ySize = (float)(yMax - yMin + 1) / (float)yBinTarget;
     if (ySize - (int)ySize > 0) {
-        ySize += 1;
+        ySize += 1;
     }
 
-    psTrace(__func__, 3, "Mosaicked chip will be %dx%d\n", (int)xSize, (int)ySize);
+    psTrace(__func__, 3, "Spliced image will be %dx%d\n", (int)xSize, (int)ySize);
     psImage *mosaic = psImageAlloc((int)xSize, (int)ySize, PS_TYPE_F32); // The mosaic image
     psImageInit(mosaic, 0.0);
 
-    // Set the image
+    // Next pass through the images to do the mosaicking
+    for (int i = 0; i < source->n; i++) {
+        psImage *image = source->data[i]; // The image of interest
+        if (xBinSource->data.S32[i] == xBinTarget && yBinSource->data.S32[i] == yBinTarget &&
+            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], "+");
+        } 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;
+                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) /
+                        xBinTarget;
+
+                    // In case the original image is binned but the mosaic is not, we need to fill in the
+                    // values in the mosaic.
+                    for (int j = 0; j < yBinSource->data.S32[i]; j++) {
+                        int yTarget = (int)(yTargetBase + yParity * (float)j / (float)yBinTarget);
+                        for (int i = 0; i < xBinSource->data.S32[i]; i++) {
+                            int xTarget = (int)(xTargetBase + xParity * (float)i / (float)xBinTarget);
+
+                            mosaic->data.F32[yTarget][xTarget] += image->data.F32[y][x];
+                        }
+                    } // Iterating over mosaic image for binned input image
+                }
+            } // Iterating over input image
+        }
+    }
+
+    return mosaic;
+}
+
+psImage *pmChipMosaic(pmChip *chip,     // Chip to mosaic
+                      int xBinChip, int yBinChip // Binning of mosaic image in x and y
+    )
+{
+    psArray *cells = chip->cells;       // The array of cells
+    psArray *images = psArrayAlloc(cells->n); // Array of images that will 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?
+
+    // Set up the required inputs
     for (int i = 0; i < cells->n; i++) {
-        pmCell *cell = cells->data[i];  // The cell of interest
-        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);
-        }
-        psImage *image = ((pmReadout*)readouts->data[0])->image; // The image to put into the mosaic
-        psImage *trimmed = psImageSubset(image, *(psRegion*)(trimsec->data[i])); // Trimmed image (no overscan)
+        pmCell *cell = cells->data[i];  // The cell of interest
+        x0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
+        y0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
+        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) {
+            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) {
+            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;
+        }
 
-        if (xBin->data.S32[i] == xBinChip && yBin->data.S32[i] == yBinChip && xParity->data.S32[i] == 1 &&
-            yParity->data.S32[i] == 1) {
-            // Let someone else do the hard work; useful to test psImageOverlaySection if no other reason
-            psImageOverlaySection(mosaic, trimmed, x0->data.S32[i], y0->data.S32[i], "+");
-        } else {
-            // We have to do the hard work ourself
-            for (int y = 0; y < trimmed->numRows; y++) {
-                float yTargetBase = (y0->data.S32[i] + yParity->data.S32[i] * yBin->data.S32[i] * y) /
-                    yBinChip;
-                for (int x = 0; x < trimmed->numCols; x++) {
-                    float xTargetBase = (x0->data.S32[i] + xParity->data.S32[i] * xBin->data.S32[i] * x) /
-                        xBinChip;
+        // Trim the image to get rid of the overscan
+        psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC");
+        psTrace(__func__, 7, "Cell %d trimsec: [%.0f:%.0f,%.0f:%.0f]\n", i, trimsec->x0, trimsec->x1,
+                trimsec->y0, trimsec->y1);
+        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);
+        }
+        psImage *image = ((pmReadout*)readouts->data[0])->image; // The image to put into the mosaic
+        images->data[i] = psImageSubset(image, *trimsec); // Trimmed image
+    }
 
-                    // In case the original image is binned but the mosaic is not, we need to fill in the
-                    // values in the mosaic.
-                    for (int j = 0; j < yBin->data.S32[i]; j++) {
-                        int yTarget = (int)(yTargetBase + yParity->data.S32[i] * (float)j / (float)yBinChip);
-                        for (int i = 0; i < xBin->data.S32[i]; i++) {
-                            int xTarget = (int)(xTargetBase +
-                                                xParity->data.S32[i] * (float)i / (float)xBinChip);
+    // Mosaic the images together and we're done
+    psImage *mosaic = p_pmImageMosaic(images, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0);
 
-                            mosaic->data.F32[yTarget][xTarget] += trimmed->data.F32[y][x];
-                        }
-                    } // Iterating over mosaic image for binned input image
-                }
-            } // Iterating over input image
-        }
-        psFree(trimmed);
-    } // Iterating over cells
-
+    // Clean up
     psFree(x0);
     psFree(y0);
     psFree(xBin);
     psFree(yBin);
-    psFree(xParity);
-    psFree(yParity);
-    psFree(trimsec);
+    psFree(xFlip);
+    psFree(yFlip);
+    psFree(images);
 
     return mosaic;