Changeset 23992

Timestamp:

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

Author:

Paul Price

Message:

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

Location:

trunk/ppSkycell/src

Files:

• Makefile.am (modified) (1 diff)
• ppSkycell.c (modified) (3 diffs)
• ppSkycell.h (modified) (1 diff)
• ppSkycellArguments.c (modified) (2 diffs)
• ppSkycellCamera.c (modified) (2 diffs)
• ppSkycellLoop.c (modified) (9 diffs)
• ppSkycellVersion.c (modified) (2 diffs)

trunk/ppSkycell/src/Makefile.am

@@ -31 +31 @@
         ppSkycell.c             \
         ppSkycellArguments.c    \
+        ppSkycellCamera.c       \
+        ppSkycellData.c         \
         ppSkycellLoop.c         \
-        ppSkycellReadout.c      \
         ppSkycellVersion.c            
 

trunk/ppSkycell/src/ppSkycell.c

@@ -5 +5 @@
 #include <stdio.h>
 #include <pslib.h>
+#include <psmodules.h>
+
+#include "ppSkycell.h"
 
 int main(int argc, char *argv[])
@@ -11 +14 @@
     if (!data) {
         psErrorStackPrint(stderr, "Unable to initialise.");
-        exit(PS_EXIT_CONFIG_ERROR);
+        return PS_EXIT_CONFIG_ERROR;
     }
 
@@ -17 +20 @@
         psErrorStackPrint(stderr, "Unable to parse arguments.");
         psFree(data);
-        exit(PS_EXIT_CONFIG_ERROR);
+        return PS_EXIT_CONFIG_ERROR;
     }
 
+    if (!ppSkycellCamera(data)) {
+        psErrorStackPrint(stderr, "Unable to parse camera configuration.");
+        psFree(data);
+        return PS_EXIT_CONFIG_ERROR;
+    }
 
+    if (!ppSkycellLoop(data)) {
+        psErrorStackPrint(stderr, "Unable to process data.");
+        psFree(data);
+        return PS_EXIT_DATA_ERROR;
+    }
 
+    psFree(data);
+
+    return PS_EXIT_SUCCESS;
+}
+

trunk/ppSkycell/src/ppSkycell.h

@@ -2 +2 @@
 #define PP_SKYCELL_H
 
+#include <pslib.h>
+#include <psmodules.h>
+
+#define PPSKYCELL_RECIPE "PPSKYCELL"    // Recipe name
+
+// Data for processing
 typedef struct {
-    psString inName;                    // Input filename
-    psString outName;                   // Output name
+    psString imagesName;                // Filename with images
+    psString masksName;                 // Filename with masks
+    psString outRoot;                   // Output root name
     int numInputs;                      // Number of inputs
+    psImageMaskType maskVal;            // Value to mask
     int bin1, bin2;                     // Binning factors
     pmConfig *config;                   // Configuration
 } ppSkycellData;
 
+/// Initialise data for processing
+ppSkycellData *ppSkycellDataInit(int *argc, char *argv[] // Command-line arguments
+    );
+
+/// Parse command-line arguments
+bool ppSkycellArguments(ppSkycellData *data, // Data for processing
+                        int argc, char *argv[] // Command-line arguments
+    );
+
+/// Parse camera configurations
+bool ppSkycellCamera(ppSkycellData *data // Data for processing
+    );
+
+/// Loop over input data, processing
+bool ppSkycellLoop(ppSkycellData *data // Data for processing
+    );
+
+
+#endif

trunk/ppSkycell/src/ppSkycellArguments.c

@@ -29 +29 @@
 {
     fprintf(stderr, "\nPan-STARRS skycell JPEGifier\n\n");
-    fprintf(stderr, "Usage: %s INPUT.list OUTPUT_ROOT\n\n",
+    fprintf(stderr, "Usage: %s -images INPUT.list [-masks MASK.list] OUTPUT_ROOT\n\n",
             program);
     fprintf(stderr, "\n");
@@ -45 +45 @@
 bool ppSkycellArguments(ppSkycellData *data, int argc, char *argv[])
 {
-    assert(config);
+    assert(data);
+    assert(data->config);
 
     psMetadata *arguments = psMetadataAlloc(); // Command-line arguments
-    if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 3) {
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-images", 0, "Filename with input images", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-masks", 0, "Filename with input masks", NULL);
+    if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 2) {
         usage(argv[0], arguments, data);
     }
 
-    data->inName = psStringCopy(argv[1]);
-    data->outRoot = psStringCopy(argv[2]);
+    data->imagesName = psMetadataLookupStr(NULL, arguments, "-images");
+    data->masksName = psMetadataLookupStr(NULL, arguments, "-masks");
+    data->outRoot = psStringCopy(argv[1]);
 
     psTrace("ppSkycell", 1, "Done reading command-line arguments\n");
     psFree(arguments);
+
+    PS_ASSERT_STRING_NON_EMPTY(data->imagesName, false);
+    PS_ASSERT_STRING_NON_EMPTY(data->outRoot, false);
 
     return true;

trunk/ppSkycell/src/ppSkycellCamera.c

@@ -9 +9 @@
 #include "ppSkycell.h"
 
+/// Read list of images
+static psArray *fileList(const char *filename // Filename
+    )
+{
+    psString input = psSlurpFilename(filename);
+    if (!input) {
+        psError(PS_ERR_IO, false, "Unable to read %s", filename);
+        return false;
+    }
+    psArray *inputs = psStringSplitArray(input, "\n", false); // Input filenames
+    psFree(input);
+    if (!inputs || inputs->n == 0) {
+        psError(PS_ERR_IO, false, "Unable to read filenames from %s", filename);
+        psFree(inputs);
+        return NULL;
+    }
+    return inputs;
+}
+
 /// Add a single filename to the arguments as an array, so that it can be used with pmFPAfileBindFromArgs, etc
-static void fileList(const char *file, // The symbolic name for the file
-                     const char *name, // The name of the file
-                     const char *comment, // Description of the file
-                     pmConfig *config // Configuration
+static void fileArguments(const char *file, // The symbolic name for the file
+                          const char *name, // The name of the file
+                          const char *comment, // Description of the file
+                          pmConfig *config // Configuration
     )
 {
     psArray *files = psArrayAlloc(1); // Array with file names
     files->data[0] = psStringCopy(name);
+    if (psMetadataLookup(config->arguments, file)) {
+        psMetadataRemoveKey(config->arguments, file);
+    }
     psMetadataAddArray(config->arguments, PS_LIST_TAIL, file, 0, comment, files);
     psFree(files);
@@ -27 +49 @@
     )
 {
-    psString input = psSlurpFilename(data->inName);
-    if (!input) {
-        psError(PS_ERR_IO, false, "Unable to read %s", data->inName);
+    psArray *images = fileList(data->imagesName); // Image names
+    if (!images) {
+        psError(psErrorCodeLast(), false, "No images provided.");
         return false;
     }
+    data->numInputs = images->n;
 
-    psArray *inputs = psStringSplitArray(input, "\n", false); // Input filenames
-    psFree(input);
-    if (!inputs || inputs->n == 0) {
-        psError(PS_ERR_IO, false, "Unable to read filenames from %s", data->inName);
-        psFree(inputs);
-        return false;
-    }
-
-    data->numInputs = inputs->n;
-    psMetadataAddArray(config->arguments, PS_LIST_TAIL, "INPUTS", 0, "Input files", inputs);
-    psFree(inputs);
-    psMetadataAddString(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Output root", data->outRoot);
-
-    pmFPAfile *inFile = NULL;           // Representative input file
-    for (int i = 0; i < data->numInputs; i++) {
-        bool found = false;             // Found file definition?
-        pmFPAfile *file = pmFPAfileDefineSingleFromArgs(&found, data->config, "PPSKYCELL.INPUT",
-                                                        "INPUTS", i); // Input file
-        if (!file || !found) {
-            psError(psErrorCodeLast(), false, "Unable to define input %d\n", i);
+    psArray *masks = NULL;              // Mask names
+    if (data->masksName) {
+        masks = fileList(data->masksName);
+        if (!masks) {
+            psError(psErrorCodeLast(), false, "No masks provided.");
+            psFree(images);
             return false;
         }
-        if (!inFile) {
-            inFile = file;
+        if (masks->n != data->numInputs) {
+            psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Number of images (%ld) and masks (%ld) do not match",
+                    images->n, masks->n);
+            psFree(images);
+            psFree(masks);
+            return false;
         }
     }
 
-    pmFPAfile *outFile = pmFPAfileDefineFromFile(data->config, inFile, data->bin1, data->bin1,
-                                                 "PPSKYCELL.JPEG1"); // Output file
-    if (!outFile) {
+    psMetadataAddStr(data->config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Output root", data->outRoot);
+
+    for (int i = 0; i < data->numInputs; i++) {
+        bool status = false;             // Status of file definition
+        fileArguments("IMAGE", images->data[i], "Name of the image", data->config);
+        pmFPAfile *image = pmFPAfileDefineFromArgs(&status, data->config, "PPSKYCELL.IMAGE", "IMAGE"); // File
+        if (!status || !image) {
+            psError(PS_ERR_IO, false, "Failed to build file from PPSKYCELL.IMAGE");
+            // XXX Cleanup
+            return false;
+        }
+
+        if (data->masksName) {
+            fileArguments("MASK", masks->data[i], "Name of the mask", data->config);
+            if (!pmFPAfileBindFromArgs(&status, image, data->config, "PPSKYCELL.MASK", "MASK") || !status) {
+                psError(PS_ERR_IO, false, "Failed to build file from PPSKYCELL.MASK");
+                // XXX Cleanup
+                return false;
+            }
+        }
+    }
+
+    if (!pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.JPEG1")) {
         psError(psErrorCodeLast(), false, "Unable to define output.");
         return false;
     }
-    if (!pmFPAfileDefineFromFile(data->config, outFile, data->bin2, data->bin2, "PPSKYCELL.JPEG2")) {
+    if (!pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.JPEG2")) {
         psError(psErrorCodeLast(), false, "Unable to define output.");
         return false;

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

trunk/ppSkycell/src/ppSkycellVersion.c

@@ -18 +18 @@
 #include <pslib.h>
 #include <psmodules.h>
-#include <ppStats.h>
 
 #include "ppSkycell.h"
@@ -78 +77 @@
     psLibVersionHeader(header);
     psModulesVersionHeader(header);
-    ppStatsVersionHeader(header);
 
     psString version = ppSkycellVersion(); // Software version