Changeset 16367 for branches/pap_branch_080207/ppStack/src/ppStackLoop.c

Timestamp:

Feb 7, 2008, 6:37:16 PM (18 years ago)

Author:

Paul Price

Message:

Making progress with refactoring of ppStack for optimisation. This version of the code does NOT yet work.

File:

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

branches/pap_branch_080207/ppStack/src/ppStackLoop.c

@@ -10 +10 @@
 
 #include "ppStack.h"
+
+// Here follows lists of files for activation/deactivation at various stages.  Each must be NULL-terminated.
+
+// All files in the system
+static char *allFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT",
+                            "PPSTACK.OUTCONV", "PPSTACK.OUTCONV.MASK", "PPSTACK.OUTCONV.WEIGHT",
+                            "PPSTACK.INCONV", "PPSTACK.INCONV.MASK", "PPSTACK.INCONV.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 };
+
+// Files required in preparation for convolution
+static char *prepareFiles[] = { "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES", 0 };
+
+// Files for convolution; these should be turned on one by one
+static char *convolutionFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT",
+                                    "PPSTACK.OUTCONV", "PPSTACK.OUTCONV.MASK", "PPSTACK.OUTCONV.WEIGHT",
+                                    0 };
+
+// Input files for the combination
+static char *inCombineFiles = { "PPSTACK.INCONV", "PPSTACK.INCONV.MASK", "PPSTACK.INCONV.WEIGHT", 0 };
+
+// Output files for the combination and photometry
+static char *outCombineFiles = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT",
+                                 "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 };
+
+
+// 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;
+}
+
+// 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;
+}
+
+
+
+// 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;
+}
+
 
 bool ppStackLoop(pmConfig *config)
@@ -32 +177 @@
     }
 
-    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");
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!");
         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)) {
+    // Preparation iteration: Load the sources, and get a target PSF model
+    psArray *sources = NULL;            // 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;
+        }
+
+        pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES"); // Readout
+        // We want to hang on to the 'sources' even when its host FPA is blown away
+        sources = psMemIncrRefCounter(psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.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
+        int fileNum = 0;        // Number of file
+        psList *psfList = psListAlloc(NULL); // List of PSFs for PSF envelope
+        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, ro->parent->parent->analysis, "PSPHOT.PSF"); // PSF
+            if (!psf) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to find PSF.");
+                psFree(view);
+                psFree(sources);
                 return false;
             }
-
-            // Put version information into the header
-            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)) {
-                    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;
+            psListAdd(psfList, PS_LIST_TAIL, psf);
+        }
+        psFree(fileIter);
+        psFree(view);
+
+        psArray *psfArray = psListToArray(psfList); // Array of PSFs
+        psFree(psfList);
+        targetPSF = pmPSFEnvelope(numCols, numRows, psfArray, psfInstances, psfRadius, psfModel,
+                                  psfOrder, psfOrder);
+        psFree(psfArray);
+        if (!targetPSF) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
+            psFree(sources);
+            return false;
+        }
+
+        filesIterateUp(config);
+    }
+
+
+    // Generate convolutions and write them to disk
+    psArray *cells = psArrayAlloc(numFiles);    // Cells for convolved images --- a handle for reading again
+    for (int i = 0; i < numFiles; i++) {
+        pmFPAfileActivate(config->files, false, NULL);
+        psArray *files = fileActivationSingle(config, prepareFiles, true, i);
+        pmFPAview *view = filesIterateDown(config);
+        if (!view) {
+            return false;
+        }
+
+        pmReadout *readout = pmFPAviewThisReadout(view, files->data[0]); // Input readout
+        psFree(view);
+
+        // Background subtraction, scaling and normalisation is performed automatically by the image matching
+        if (!ppStackMatch(ro, sources, outPSF, config)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to match image %d --- ignoring.", i);
+            psFree(convolved);
+            return false;
+        }
+
+#ifdef CONVOLUTION_FILES
+        if (readout->image) {
+            psString name = NULL;           // Name of image
+            psStringAppend(&name, "convolved%03d_image.fits", i);
+            psFits *fits = psFitsOpen(name, "w");
+            psFree(name);
+            psFitsWriteImage(fits, NULL, readout->image, 0, NULL);
+            psFitsClose(fits);
+        }
+        if (readout->mask) {
+            psString name = NULL;           // Name of image
+            psStringAppend(&name, "convolved%03d_mask.fits", i);
+            psFits *fits = psFitsOpen(name, "w");
+            psFree(name);
+            psFitsWriteImage(fits, NULL, readout->mask, 0, NULL);
+            psFitsClose(fits);
+        }
+        if (readout->weight) {
+            psString name = NULL;           // Name of image
+            psStringAppend(&name, "convolved%03d_weight.fits", i);
+            psFits *fits = psFitsOpen(name, "w");
+            psFree(name);
+            psFitsWriteImage(fits, NULL, readout->weight, 0, NULL);
+            psFitsClose(fits);
+        }
+#endif // CONVOLUTION_FILES
+
+        cells->data[i] = psMemIncrRefCounter(readout->parent);
+        filesIterateUp(config);
+    }
+
+    // Set files to read at the readout level --- then when we iterate down, the pmFPAfileIO stuff should take
+    // care of the basic things, and we can just worry about grabbing the readouts by chunks
+    fileSetDataLevel(config, inCombineFiles, PM_FPA_LEVEL_CHUNK);
+    pmFPAfileActivate(config->files, false, NULL);
+    fileActivation(inCombineFiles, true);
+    fileActivation(outCombineFiles, true);
+
+    {
+        pmFPAview *view = filesIterateDown(config);
+        if (!view) {
+            psFree(cells);
+            return false;
+        }
+        pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSTACK.OUTPUT"); // Output readout
+        psFree(view);
+
+        psArray *readouts = psArrayAlloc(numFiles); // Readouts for convolved images
+        for (int i = 0; i < numFiles; i++) {
+            readouts->data[i] = pmReadoutAlloc(cells->data[i]); // Readout into which to read
+        }
+
+        // Read convolutions by chunks
+        int numRead = 0;                    // Number of inputs read
+        int numChunk = 0;                   // Chunk number
+        bool more = true;               // More to read?
+        for (int numChunk = 0; more; numChunk++) {
+            for (int i = 0; i < numFiles; i++) {
+                pmReadout *readout = readouts->data[i];
+                assert(readout);
+
+                // Files, containing the FITS handle
+                pmFPAfile *imageFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV", i);
+                pmFPAfile *maskFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.MASK", i);
+                pmFPAfile *weightFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.WEIGHT", i);
+
+                // FITS handles
+                psFits *fitsImage = imageFile->fits;
+                psFits *fitsMask = maskFile->fits;
+                psFits *fitsWeight = weightFile->fits;
+
+                if (!pmReadoutReadChunk(readout, fitsImage, 0, numScans, overlap) ||
+                    !pmReadoutReadChunkMask(readout, fitsMask, 0, numScans, overlap) ||
+                    !pmReadoutReadChunkWeight(readout, fitsWeight, 0, numScans, overlap)) {
+                    psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i);
+                    psFree(readouts);
+                    psFree(cells);
                 }
             }
 
-            // Perform statistics on the cell
-            if (stats) {
-                ppStatsFPA(stats, output->fpa, view, maskBlank, config);
-            }
-
-            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+            if (!ppStackReadout(config, outRO, readouts)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
                 return false;
             }
-        }
-
-        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-            return false;
-        }
-    }
-
-    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-        return false;
+
+            for (int i = 0; i < numFiles && more; i++) {
+                pmReadout *readout = readouts->data[i];
+                assert(readout);
+
+                // Files, containing the FITS handle
+                pmFPAfile *imageFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV", i);
+                pmFPAfile *maskFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.MASK", i);
+                pmFPAfile *weightFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.WEIGHT", i);
+
+                // FITS handles
+                psFits *fitsImage = imageFile->fits;
+                psFits *fitsMask = maskFile->fits;
+                psFits *fitsWeight = weightFile->fits;
+
+                more &= pmReadoutMore(readout, fitsImage, 0, numScans, overlap);
+                more &= pmReadoutMoreMask(readout, fitsMask, 0, numScans, overlap);
+                more &= pmReadoutMoreWeight(readout, fitsWeight, 0, numScans, overlap);
+            }
+        }
+
+        // Get rid of the input files
+        fileActivation(outCombineFiles, false);
+        filesIterateUp(config);
+
+        ppStackPhotometry(config, outRO);
+
+        // Put version information into the header
+        pmHDU *hdu = pmHDUFromCell(outRO->parent);
+        if (!hdu) {
+            psError(PS_ERR_UNKNOWN, true, "Unable to find HDU for output.");
+            return false;
+        }
+        if (!hdu->header) {
+            hdu->header = psMetadataAlloc();
+        }
+        ppStackVersionMetadata(hdu->header);
+
+        // Statistics on output
+        if (stats) {
+            ppStatsFPA(stats, output->fpa, view, maskBlank, config);
+        }
+
+        // Write out the output files
+        filesActivation(outCombineFiles, true);
+        filesIterateUp(config);
     }