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

Timestamp:

Feb 8, 2008, 11:43:52 AM (18 years ago)

Author:

Paul Price

Message:

Code compiles and is ready for preliminary testing. Some features (e.g., rejection based on initial combination) are disabled for now, marked with XXX in the code.

File:

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

branches/pap_branch_080207/ppStack/src/ppStackLoop.c

@@ -13 +13 @@
 // 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",
@@ -23 +24 @@
                             "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 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 };
+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 };
+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 };
 
 
@@ -103 +100 @@
 
 
-
 // Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells
 static pmFPAview *filesIterateDown(pmConfig *config // Configuration
@@ -158 +154 @@
 bool ppStackLoop(pmConfig *config)
 {
+    assert(config);
+
     psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
 
@@ -182 +180 @@
         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
+
+
+
+    // XXX Will need to adjust this when we're doing the convolution
+    int overlap = 0;                    // Overlap between consecutive scans
+
+
+
 
     // Preparation iteration: Load the sources, and get a target PSF model
-    psArray *sources = NULL;            // Sources to use for PSF matching
+    pmReadout *sources = NULL;          // Readout with sources to use for PSF matching
     pmPSF *targetPSF = NULL;            // Target PSF
     {
@@ -194 +202 @@
         }
 
-        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"));
+        sources = psMemIncrRefCounter(pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES"));
         if (!sources) {
             psError(PS_ERR_UNKNOWN, true, "Unable to find sources.");
@@ -207 +214 @@
                                                                "^PPSTACK.INPUT$"); // Iterator
         psMetadataItem *fileItem; // Item from iteration
-        int fileNum = 0;        // Number of file
         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, ro->parent->parent->analysis, "PSPHOT.PSF"); // 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;
             }
             psListAdd(psfList, PS_LIST_TAIL, psf);
+
+            pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest
+            pmHDU *hdu = pmHDUFromCell(cell);
+            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);
 
-        psArray *psfArray = psListToArray(psfList); // Array of PSFs
+        targetPSF = ppStackPSF(config, numCols, numRows, psfList);
         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.");
@@ -241 +265 @@
 
     // 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++) {
+    psArray *cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again
+    for (int i = 0; i < num; i++) {
         pmFPAfileActivate(config->files, false, NULL);
         psArray *files = fileActivationSingle(config, prepareFiles, true, i);
         pmFPAview *view = filesIterateDown(config);
         if (!view) {
+            psFree(sources);
+            psFree(targetPSF);
             return false;
         }
@@ -254 +280 @@
 
         // Background subtraction, scaling and normalisation is performed automatically by the image matching
-        if (!ppStackMatch(ro, sources, outPSF, config)) {
+        if (!ppStackMatch(readout, sources, targetPSF, config)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to match image %d --- ignoring.", i);
-            psFree(convolved);
+            psFree(sources);
+            psFree(targetPSF);
             return false;
         }
@@ -290 +317 @@
         filesIterateUp(config);
     }
+    psFree(sources);
+    psFree(targetPSF);
 
     // Set files to read at the readout level --- then when we iterate down, the pmFPAfileIO stuff should take
@@ -295 +324 @@
     fileSetDataLevel(config, inCombineFiles, PM_FPA_LEVEL_CHUNK);
     pmFPAfileActivate(config->files, false, NULL);
-    fileActivation(inCombineFiles, true);
-    fileActivation(outCombineFiles, true);
+    fileActivation(config, inCombineFiles, true);
+    fileActivation(config, outCombineFiles, true);
 
     {
@@ -304 +333 @@
             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++) {
+        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);
 
         // 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++) {
+            for (int i = 0; i < num; i++) {
                 pmReadout *readout = readouts->data[i];
                 assert(readout);
@@ -336 +364 @@
                     psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i);
                     psFree(readouts);
-                    psFree(cells);
+                    psFree(outRO);
+                    psFree(view);
+                    return false;
                 }
             }
 
+#if 0
             if (!ppStackReadout(config, outRO, readouts)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
+                psFree(readouts);
+                psFree(outRO);
+                psFree(view);
                 return false;
             }
-
-            for (int i = 0; i < numFiles && more; i++) {
+#else
+            printf("Stack...\n");
+#endif
+
+            for (int i = 0; i < num && more; i++) {
                 pmReadout *readout = readouts->data[i];
                 assert(readout);
@@ -364 +401 @@
             }
         }
+        psFree(readouts);
 
         // Get rid of the input files
-        fileActivation(outCombineFiles, false);
+        fileActivation(config, outCombineFiles, false);
         filesIterateUp(config);
 
-        ppStackPhotometry(config, outRO);
+        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, true, "Unable to find HDU for output.");
+            psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output.");
             return false;
         }
@@ -382 +431 @@
         ppStackVersionMetadata(hdu->header);
 
-        // Statistics on output
-        if (stats) {
-            ppStatsFPA(stats, output->fpa, view, maskBlank, config);
-        }
-
         // Write out the output files
-        filesActivation(outCombineFiles, true);
+        fileActivation(config, outCombineFiles, true);
         filesIterateUp(config);
-    }
-
-    psFree(view);
+
+        psFree(outRO);
+        psFree(view);
+    }
+
 
     // Write out summary statistics