Changeset 16605 for trunk/ppStack/src/ppStackLoop.c

Timestamp:

Feb 22, 2008, 9:28:00 AM (18 years ago)

Author:

Paul Price

Message:

Merging in branch development. ppStack now works with incremental reads. Needs some cleanup, but it works.

File:

• trunk/ppStack/src/ppStackLoop.c (modified) (2 diffs)

trunk/ppStack/src/ppStackLoop.c

@@ -11 +11 @@
 #include "ppStack.h"
 
+// Here follows lists of files for activation/deactivation at various stages.  Each must be NULL-terminated.
+
+#if 0
+// All files in the system
+static char *allFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT",
+                            "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT",
+                            "PSPHOT.PSF.SAVE", "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES",
+                            "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL", "PSPHOT.BACKMDL.STDEV",
+                            "PSPHOT.BACKGND", "PSPHOT.BACKSUB", "SOURCE.PLOT.MOMENTS",
+                            "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID", "PSPHOT.INPUT.CMF",
+                            0 };
+#endif
+
+// Files required in preparation for convolution
+static char *prepareFiles[] = { "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES", 0 };
+
+// Files required for the convolution
+static char *convolveFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT", 0 };
+
+// Output files for the combination
+static char *combineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT", 0 };
+
+// Files for photometry
+static char *photFiles[] = { "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL",
+                             "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB",
+                             "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID",
+                             "PSPHOT.INPUT.CMF", 0 };
+
+//#define CONVOLVED_ALREADY               // Already have the convolution products --- testing
+
+
+// Activate/deactivate a list of files
+static void fileActivation(pmConfig *config, // Configuration
+                           char **files, // Files to turn on/off
+                           bool state   // Activation state
+    )
+{
+    assert(config);
+    assert(files);
+
+    for (int i = 0; files[i] != NULL; i++) {
+        pmFPAfileActivate(config->files, state, files[i]);
+    }
+    return;
+}
+
+// Activate/deactivate a single element for a list; return array of files
+static psArray *fileActivationSingle(pmConfig *config, // Configuration
+                                     char **files, // Files to turn on/off
+                                     bool state,   // Activation state
+                                     int num // Number of file in sequence
+                                     )
+{
+    assert(config);
+    assert(files);
+
+    psList *list = psListAlloc(NULL);   // List of files
+
+    for (int i = 0; files[i] != NULL; i++) {
+        pmFPAfile *file = pmFPAfileActivateSingle(config->files, state, files[i], num); // Activated file
+        psListAdd(list, PS_LIST_TAIL, file);
+    }
+
+    psArray *array = psListToArray(list);
+    psFree(list);
+
+    return array;
+}
+
+#if 0
+// Set the data level for a list of files
+static void fileSetDataLevel(pmConfig *config, // Configuration
+                             char **files, // Files for which to set level
+                             pmFPALevel level // Level to set
+                             )
+{
+    assert(config);
+    assert(files);
+
+    for (int i = 0; files[i] != NULL; i++) {
+        psArray *selected = pmFPAfileSelect(config->files, files[i]); // Selected files of interest
+        for (int j = 0; j < selected->n; j++) {
+            pmFPAfile *file = selected->data[j];
+            assert(file);
+            file->dataLevel = level;
+        }
+        psFree(selected);
+    }
+    return;
+}
+#endif
+
+// 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;
+    }
+    return true;
+}
+
+#ifndef CONVOLVED_ALREADY
+// Write an image to a FITS file
+static bool writeImage(const char *name, // Name of image
+                       psMetadata *header, // Header
+                       const psImage *image // Image
+                       )
+{
+    assert(name);
+    assert(image);
+
+    psFits *fits = psFitsOpen(name, "w");
+    if (!fits) {
+        psError(PS_ERR_IO, false, "Unable to open FITS file to write image.");
+        return false;
+    }
+    if (!psFitsWriteImage(fits, header, image, 0, NULL)) {
+        psError(PS_ERR_IO, false, "Unable to write FITS image.");
+        return false;
+    }
+    psFitsClose(fits);
+    return true;
+}
+#endif
+
+
 bool ppStackLoop(pmConfig *config)
 {
+    assert(config);
+
     psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
 
@@ -32 +202 @@
     }
 
-    pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSTACK.INPUT"); // Token input file
-    if (!input) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find input data!\n");
-        return false;
-    }
     pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTACK.OUTPUT"); // Output file
     if (!output) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!\n");
-        return false;
-    }
-
-    pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy
-    pmHDU *lastHDU = NULL;              // Last HDU that was updated
-
-    // Iterate over the FPA hierarchy
-    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-        return false;
-    }
-
-    pmChip *chip;                       // Chip of interest
-    while ((chip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) {
-        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-            return false;
-        }
-
-        pmCell *cell;                // Cell of interest
-        while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
-            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!");
+        return false;
+    }
+    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
+    int numScans = psMetadataLookupS32(NULL, config->arguments, "ROWS"); // Number of scans to read at once
+    psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
+    int overlap = 2 * psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Overlap by kernel size between consecutive scans
+
+    // Preparation iteration: Load the sources, and get a target PSF model
+    pmReadout *sources = NULL;          // Readout with sources to use for PSF matching
+    pmPSF *targetPSF = NULL;            // Target PSF
+    {
+        pmFPAfileActivate(config->files, false, NULL);
+        fileActivation(config, prepareFiles, true);
+        pmFPAview *view = filesIterateDown(config);
+        if (!view) {
+            return false;
+        }
+
+        // We want to hang on to the 'sources' even when its host FPA is blown away
+        sources = psMemIncrRefCounter(pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES"));
+        if (!sources) {
+            psError(PS_ERR_UNKNOWN, true, "Unable to find sources.");
+            psFree(view);
+            return false;
+        }
+
+        // Generate list of PSFs
+        psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD,
+                                                               "^PPSTACK.INPUT$"); // Iterator
+        psMetadataItem *fileItem; // Item from iteration
+        psList *psfList = psListAlloc(NULL); // List of PSFs for PSF envelope
+        int numCols = 0, numRows = 0;   // Size of image
+        while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
+            assert(fileItem->type == PS_DATA_UNKNOWN);
+            pmFPAfile *inputFile = fileItem->data.V; // An input file
+            pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF
+            pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF
+            if (!psf) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to find PSF.");
+                psFree(view);
+                psFree(sources);
+                psFree(fileIter);
+                psFree(psfList);
                 return false;
             }
-
-            // Put version information into the header
+            psListAdd(psfList, PS_LIST_TAIL, psf);
+
+            pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest
             pmHDU *hdu = pmHDUFromCell(cell);
-            if (hdu && hdu != lastHDU) {
-                if (!hdu->header) {
-                    hdu->header = psMetadataAlloc();
-                }
-                ppStackVersionMetadata(hdu->header);
-                lastHDU = hdu;
-            }
-
-            pmReadout *readout;         // Readout of interest
-            while ((readout = pmFPAviewNextReadout(view, input->fpa, 1))) {
-                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+            assert(hdu && hdu->header);
+            int naxis1 = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); // Number of columns
+            int naxis2 = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows
+            if (naxis1 <= 0 || naxis2 <= 0) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to determine size of image from PSF.");
+                psFree(view);
+                psFree(sources);
+                psFree(fileIter);
+                psFree(psfList);
+                return false;
+            }
+            if (numCols == 0 && numRows == 0) {
+                numCols = naxis1;
+                numRows = naxis2;
+            }
+        }
+        psFree(fileIter);
+        psFree(view);
+
+        targetPSF = ppStackPSF(config, numCols, numRows, psfList);
+        psFree(psfList);
+        if (!targetPSF) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
+            psFree(sources);
+            return false;
+        }
+
+        filesIterateUp(config);
+    }
+
+    const char *suffix = "conv.fits";   // Suffix for convolved images; ultimately this will be from recipe
+    const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
+    assert(outName);
+    psArray *imageNames = psArrayAlloc(num);
+    psArray *maskNames = psArrayAlloc(num);
+    psArray *weightNames = psArrayAlloc(num);
+    for (int i = 0; i < num; i++) {
+        psString imageName = NULL, maskName = NULL, weightName = NULL; // Names for convolved images
+        psStringAppend(&imageName, "%s.im-%d.%s", outName, i, suffix);
+        psStringAppend(&maskName, "%s.mk-%d.%s", outName, i, suffix);
+        psStringAppend(&weightName, "%s.wt-%d.%s", outName, i, suffix);
+        imageNames->data[i] = imageName;
+        maskNames->data[i] = maskName;
+        weightNames->data[i] = weightName;
+    }
+
+    // Generate convolutions and write them to disk
+    psArray *cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again
+    psArray *subKernels = psArrayAlloc(num); // Subtraction kernels --- required in the stacking
+    psArray *subRegions = psArrayAlloc(num); // Subtraction regions --- required in the stacking
+    for (int i = 0; i < num; i++) {
+        pmFPAfileActivate(config->files, false, NULL);
+        psArray *files = fileActivationSingle(config, convolveFiles, true, i);
+        pmFPAview *view = filesIterateDown(config);
+        if (!view) {
+            psFree(sources);
+            psFree(targetPSF);
+            return false;
+        }
+
+        pmReadout *readout = pmFPAviewThisReadout(view, ((pmFPAfile*)files->data[0])->fpa); // Input readout
+        psFree(view);
+
+#ifndef CONVOLVED_ALREADY
+        // Background subtraction, scaling and normalisation is performed automatically by the image matching
+        psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
+        if (!ppStackMatch(readout, &regions, &kernels, sources, targetPSF, config)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to match image %d --- ignoring.", i);
+            psFree(sources);
+            psFree(targetPSF);
+            return false;
+        }
+
+        subRegions->data[i] = regions;
+        subKernels->data[i] = kernels;
+
+        // Write the temporary convolved files
+        pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
+        assert(hdu);
+        writeImage(imageNames->data[i],  hdu->header, readout->image);
+        writeImage(maskNames->data[i],   hdu->header, readout->mask);
+        writeImage(weightNames->data[i], hdu->header, readout->weight);
+#endif
+
+        cells->data[i] = psMemIncrRefCounter(readout->parent);
+        filesIterateUp(config);
+    }
+    psFree(sources);
+    psFree(targetPSF);
+
+    // Stack the convolved files
+    {
+        pmFPAfileActivate(config->files, false, NULL);
+        fileActivation(config, combineFiles, true);
+        if (psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY")) {
+            fileActivation(config, photFiles, true);
+        }
+        pmFPAview *view = filesIterateDown(config);
+        if (!view) {
+            psFree(cells);
+            psFree(subKernels);
+            psFree(subRegions);
+            return false;
+        }
+        pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
+        pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
+
+        psArray *readouts = psArrayAlloc(num); // Readouts for convolved images
+        for (int i = 0; i < num; i++) {
+            readouts->data[i] = pmReadoutAlloc(cells->data[i]); // Readout into which to read
+        }
+        psFree(cells);
+
+        // FITS files
+        psArray *imageFits  = psArrayAlloc(num);
+        psArray *maskFits   = psArrayAlloc(num);
+        psArray *weightFits = psArrayAlloc(num);
+        for (int i = 0; i < num; i++) {
+            imageFits->data[i] = psFitsOpen(imageNames->data[i], "r");
+            maskFits->data[i] = psFitsOpen(maskNames->data[i], "r");
+            weightFits->data[i] = psFitsOpen(weightNames->data[i], "r");
+        }
+
+        // Read convolutions by chunks
+        bool more = true;               // More to read?
+        for (int numChunk = 0; more; numChunk++) {
+            for (int i = 0; i < num; i++) {
+                pmReadout *readout = readouts->data[i];
+                assert(readout);
+
+                if (!pmReadoutReadChunk(readout, imageFits->data[i], 0, numScans, overlap) ||
+                    !pmReadoutReadChunkMask(readout, maskFits->data[i], 0, numScans, overlap) ||
+                    !pmReadoutReadChunkWeight(readout, weightFits->data[i], 0, numScans, overlap)) {
+                    psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i);
+                    psFree(readouts);
+                    psFree(subKernels);
+                    psFree(subRegions);
+                    psFree(outRO);
+                    psFree(view);
                     return false;
                 }
-
-                // Perform the analysis
-                if (!ppStackReadout(config, view)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
-                    return false;
-                }
-
-                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-                    return false;
-                }
-            }
-
-            // Perform statistics on the cell
-            if (stats) {
-                ppStatsFPA(stats, output->fpa, view, maskBlank, config);
-            }
-
-            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+            }
+
+#ifdef TESTING
+            {
+                pmReadout *ro = readouts->data[0];
+                psTrace("ppStack", 1, "Stack: [%d:%d,%d:%d]\n", ro->col0, ro->col0 + ro->image->numCols,
+                        ro->row0, ro->row0 + ro->image->numRows);
+            }
+#endif
+
+            if (!ppStackReadout(config, outRO, readouts, subRegions, subKernels)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
+                psFree(readouts);
+                psFree(subKernels);
+                psFree(subRegions);
+                psFree(outRO);
+                psFree(view);
                 return false;
             }
-        }
-
-        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-            return false;
-        }
-    }
-
-    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-        return false;
-    }
-
-    psFree(view);
+
+            for (int i = 0; i < num && more; i++) {
+                pmReadout *readout = readouts->data[i];
+                assert(readout);
+                more &= pmReadoutMore(readout, imageFits->data[i], 0, numScans);
+                more &= pmReadoutMoreMask(readout, maskFits->data[i], 0, numScans);
+                more &= pmReadoutMoreWeight(readout, weightFits->data[i], 0, numScans);
+            }
+        }
+
+        psFree(readouts);
+        psFree(subKernels);
+        psFree(subRegions);
+        for (int i = 0; i < num; i++) {
+            psFitsClose(imageFits->data[i]);
+            psFitsClose(maskFits->data[i]);
+            psFitsClose(weightFits->data[i]);
+        }
+
+        if (psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY") &&
+            !ppStackPhotometry(config, outRO, view)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output.");
+            psFree(outRO);
+            psFree(view);
+            return false;
+        }
+
+        // Statistics on output
+        if (stats) {
+            ppStatsFPA(stats, outCell->parent->parent, view, maskBlank, config);
+        }
+
+        // Put version information into the header
+        pmHDU *hdu = pmHDUFromCell(outRO->parent);
+        if (!hdu) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output.");
+            return false;
+        }
+        if (!hdu->header) {
+            hdu->header = psMetadataAlloc();
+        }
+        ppStackVersionMetadata(hdu->header);
+
+        // Write out the output files
+        fileActivation(config, combineFiles, true);
+        filesIterateUp(config);
+
+        psFree(outRO);
+        psFree(view);
+    }
+
 
     // Write out summary statistics