Changeset 23992 for trunk/ppSkycell/src/ppSkycellLoop.c

Timestamp:

Apr 28, 2009, 3:31:25 PM (17 years ago)

Author:

Paul Price

Message:

It builds. Hasn't been tested yet, though.

File:

• trunk/ppSkycell/src/ppSkycellLoop.c (modified) (9 diffs)

trunk/ppSkycell/src/ppSkycellLoop.c

@@ -11 +11 @@
 #define BUFFER 16                       // Size of buffer for projections
 
-static psRegion *skycellRegion(psRegion *region, // Current region of skycell, or NULL
-                               pmAstromWCS *wcs, // World Coordinate System
+// List of input files
+static const char *inFiles[] = { "PPSKYCELL.IMAGE", "PPSKYCELL.MASK", NULL };
+
+// List of output files
+//static const char *outFiles[] = { "PPSKYCELL.JPEG1", "PPSKYCELL.JPEG2", NULL };
+
+
+static void regionMinMax(psRegion *base,// Base region; modified
+                         const psRegion *compare // Comparison region
+                         )
+{
+    base->x0 = PS_MIN(base->x0, compare->x0);
+    base->x1 = PS_MAX(base->x1, compare->x1);
+    base->y0 = PS_MIN(base->y0, compare->y0);
+    base->y1 = PS_MAX(base->y1, compare->y1);
+    return;
+}
+
+static psRegion *skycellRegion(pmAstromWCS *wcs, // World Coordinate System
                                int numCols, int numRows // Size of image
     )
 {
     psPlane *fromCoords = psPlaneAlloc(), *toCoords = psPlaneAlloc(); // Coordinates for transforms
-    if (!region) {
-        region = psRegionAlloc(INFINITY, -INFINITY, INFINITY, -INFINITY);
-    }
+    psRegion *region = psRegionAlloc(INFINITY, -INFINITY, INFINITY, -INFINITY); // Region for skycell
 
 // Limit the region using the nominated point (X,Y)
@@ -25 +40 @@
     fromCoords->x = (X); \
     fromCoords->y = (Y); \
-    psPlaneTransform(toCoords, wcs->transform, fromCoords); \
+    psPlaneTransformApply(toCoords, wcs->trans, fromCoords); \
     toCoords->x += wcs->crpix1; \
     toCoords->y += wcs->crpix2; \
@@ -33 +48 @@
     region->y1 = PS_MAX(region->y1, toCoords->y);
 
-    REGION_RANGE(0,0);                  // Lower left
-    REGION_RANGE(0,numRows);            // Upper left
-    REGION_RANGE(numCols,0);            // Lower right
-    REGION_RANGE(numCols,numRows);      // Upper right
+    REGION_RANGE(0, 0);                 // Lower left
+    REGION_RANGE(0, numRows);           // Upper left
+    REGION_RANGE(numCols, 0);           // Lower right
+    REGION_RANGE(numCols, numRows);     // Upper right
 
     return region;
+}
+
+// Activate/deactivate a single element for a list
+void fileActivationSingle(pmConfig *config, // Configuration
+                          const char **files, // Files to turn on/off
+                          bool state,   // Activation state
+                          int num // Number of file in sequence
+                          )
+{
+    assert(config);
+    for (int i = 0; files[i] != NULL; i++) {
+        pmFPAfileActivateSingle(config->files, state, files[i], num); // Activated file
+    }
+    return;
+}
+
+// Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells
+static pmFPAview *filesIterateDown(pmConfig *config // Configuration
+                                   )
+{
+    assert(config);
+
+    pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        return NULL;
+    }
+    view->chip = 0;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        return NULL;
+    }
+    view->cell = 0;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        return NULL;
+    }
+    view->readout = 0;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        return NULL;
+    }
+    return view;
+}
+
+// Iterate up the hierarchy, writing files; we can get away with this because we're working on skycells
+static bool filesIterateUp(pmConfig *config // Configuration
+                           )
+{
+    assert(config);
+
+    pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy
+    view->chip = view->cell = view->readout = 0;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        return false;
+    }
+    view->readout = -1;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        return false;
+    }
+    view->cell = -1;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        return false;
+    }
+    view->chip = -1;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        return false;
+    }
+    psFree(view);
+    return true;
 }
 
@@ -49 +130 @@
     psVector *cdelt1 = psVectorAllocEmpty(BUFFER, PS_TYPE_F64); // CDELT1 values
     psVector *cdelt2 = psVectorAllocEmpty(BUFFER, PS_TYPE_F64); // CDELT2 values
-    psArray *regions = psArrayAllocEmpty(BUFFER); // Region for projection
+    psArray *projRegions = psArrayAllocEmpty(BUFFER); // Region for projection
     int numProj = 0;                    // Number of projections
 
     psVector *target = psVectorAlloc(data->numInputs, PS_TYPE_S32); // Target for each input
+    psArray *imageRegions = psArrayAlloc(data->numInputs); // Region for image
 
     for (int i = 0; i < data->numInputs; i++) {
@@ -68 +150 @@
         int numRows = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows
 
-        double xLow = wcs->crpix1 / fabs(wcs->crpix1) + wcs->cdelt2 / fabs(wcs->crpix2)
-        double xHigh = numCols + wcs->crpix1 / fabs(wcs->crpix1) + numRows * wcs->cdelt2 / fabs(wcs->crpix2);
-
         pmAstromWCS *wcs = pmAstromWCSfromHeader(hdu->header); // World Coordinate System
         if (!wcs) {
@@ -76 +155 @@
             return false;
         }
+
+        psRegion *region = imageRegions->data[i] = skycellRegion(wcs, numCols, numRows); // Region of image
 
         bool found = false;             // Found a projection?
@@ -83 +164 @@
                 wcs->cdelt1 == cdelt1->data.F64[j] &&
                 wcs->cdelt2 == cdelt1->data.F64[j]) {
-                skycellRegion(regions->data[j], wcs, numCols, numRows);
+                regionMinMax(projRegions->data[j], region);
                 target->data.S32[i] = j;
+                skycellRegion(wcs, numCols, numRows);
                 found = true;
             }
@@ -94 +176 @@
             psVectorAppend(cdelt1, wcs->cdelt1);
             psVectorAppend(cdelt2, wcs->cdelt2);
-            psArrayAdd(regions, regions->n, skycellRegion(NULL, wcs, numCols, numRows))
+            psArrayAdd(projRegions, projRegions->n, region);
             target->data.S32[i] = numProj;
 
@@ -101 +183 @@
     }
 
+    pmFPAfileActivate(data->config->files, false, NULL);
+
+    for (int i = 0; i < numProj; i++) {
+        psRegion *projRegion = projRegions->data[i]; // Region for skycell projection
+        int xSize = projRegion->x1 - projRegion->x0; // Size of unbinned image
+        int ySize = projRegion->y1 - projRegion->y0; // Size of unbinned image
+        // Size of binned image 1
+        int numCols1 = xSize / (float)data->bin1 + 0.5, numRows1 = ySize / (float)data->bin1 + 0.5;
+        // Size of binned image 2
+        int numCols2 = numCols1 / (float)data->bin2 + 0.5, numRows2 = numRows2 / (float)data->bin1 + 0.5;
+
+        psImage *image1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); // Binned image
+        psImage *image2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); // Binned image
+        psImageInit(image1, 0);
+        psImageInit(image2, 0);
+
+        for (int j = 0; j < data->numInputs; j++) {
+            if (target->data.S32[j] != i) {
+                continue;
+            }
+
+            fileActivationSingle(data->config, inFiles, true, j);
+            pmFPAview *view = filesIterateDown(data->config); // View to readout
+            if (!view) {
+                psError(psErrorCodeLast(), false, "Unable to iterate down.");
+                // XXX Cleanup
+                return false;
+            }
+
+            pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", j);
+            pmReadout *inRO = pmFPAviewThisReadout(view, file->fpa); // Readout with input
+            psFree(view);
+
+            pmReadout *bin1RO = pmReadoutAlloc(NULL), *bin2RO = pmReadoutAlloc(NULL); // Binned readouts
+            if (!pmReadoutRebin(bin1RO, inRO, data->maskVal, data->bin1, data->bin1)) {
+                psError(psErrorCodeLast(), false, "Unable to rebin image");
+                // XXX Cleanup
+                return false;
+            }
+            if (!pmReadoutRebin(bin2RO, bin1RO, data->maskVal, data->bin2, data->bin2)) {
+                psError(psErrorCodeLast(), false, "Unable to rebin image");
+                // XXX Cleanup
+                return false;
+            }
+
+            psRegion *imageRegion = imageRegions->data[j]; // Region for image
+            // Offsets for image on skycell
+            int xOffset1 = (imageRegion->x0 - projRegion->x0) / (float)data->bin1 + 0.5;
+            int yOffset1 = (imageRegion->y0 - projRegion->y0) / (float)data->bin1 + 0.5;
+            int xOffset2 = xOffset1 / (float)data->bin2 + 0.5, yOffset2 = yOffset1 / (float)data->bin2 + 0.5;
+
+            // XXX Completely neglecting rotations
+            // The skycells are divided up neatly with them all having the same orientation
+            psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "=");
+            psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "=");
+
+            psFree(bin1RO);
+            psFree(bin2RO);
+            filesIterateUp(data->config);
+        }
+
+        {
+            psString filename = NULL;   // Filename for image
+            psStringAppend(&filename, "skycell_%d_1.fits", i);
+            psFits *fits = psFitsOpen(filename, "w");
+            psFree(filename);
+            psFitsWriteImage(fits, NULL, image1, 0, NULL);
+            psFitsClose(fits);
+        }
+
+        {
+            psString filename = NULL;   // Filename for image
+            psStringAppend(&filename, "skycell_%d_2.fits", i);
+            psFits *fits = psFitsOpen(filename, "w");
+            psFree(filename);
+            psFitsWriteImage(fits, NULL, image2, 0, NULL);
+            psFitsClose(fits);
+        }
+
+        psFree(image1);
+        psFree(image2);
+    }
+
+    // XXX Cleanup
+    return true;
+}