Changeset 16693

Timestamp:

Feb 27, 2008, 3:19:23 PM (18 years ago)

Author:

Paul Price

Message:

Split the combination run into two --- one for the initial combination and rejection, and the second for the final combination. Forced to do this because of the convolution of the rejected pixels --- the convolution area can go beyond the individual scans, so it's necessary to gather the list of all rejected pixels before doing the final combination run. The code seems to be working rather well now.

Location:

trunk/ppStack/src

Files:

• ppStack.h (modified) (2 diffs)
• ppStackLoop.c (modified) (2 diffs)
• ppStackReadout.c (modified) (6 diffs)

trunk/ppStack/src/ppStack.h

@@ -2 +2 @@
 #define PP_STACK_H
 
-#define PPSTACK_RECIPE "PPSTACK"
+#define PPSTACK_RECIPE "PPSTACK"        // Name of the recipe
+#define PPSTACK_REJECTED_PIXELS "PPSTACK.PIXELS" // Name of rejected pixels metadata items
 
 #include <psmodules.h>
@@ -26 +27 @@
 
 // Perform stacking on a readout
-bool ppStackReadout(const pmConfig *config,   // Configuration
+bool ppStackReadoutInitial(const pmConfig *config,   // Configuration
                     pmReadout *outRO,   // Output readout
                     const psArray *readouts, // Input readouts
                     const psArray *regions, // Array with array of regions used in each PSF matching
                     const psArray *kernels // Array with array of kernels used in each PSF matching
+    );
+
+// Perform stacking on a readout
+bool ppStackReadoutFinal(const pmConfig *config,   // Configuration
+                         pmReadout *outRO,   // Output readout
+                         const psArray *readouts, // Input readouts
+                         const psArray *rejected // Array with pixels rejected in each image
     );
 

trunk/ppStack/src/ppStackLoop.c

@@ -413 +413 @@
 #endif
 
-            if (!ppStackReadout(config, outRO, readouts, subRegions, subKernels)) {
+            if (!ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
                 psFree(readouts);
@@ -432 +432 @@
         }
 
-        psFree(readouts);
+        // Reset for the second read
+        // Extract the rejection lists
         psFree(subKernels);
         psFree(subRegions);
+        psArray *rejected = psArrayAlloc(num); // Rejected pixels
+        for (int i = 0; i < num; i++) {
+            pmReadout *ro = readouts->data[i]; // Readout of interest
+            pmReadoutFreeData(ro);
+
+            psPixels *rejects = NULL;   // Rejection list for this readout
+            psMetadataIterator *iter = psMetadataIteratorAlloc(ro->analysis, PS_LIST_HEAD,
+                                                               "^" PPSTACK_REJECTED_PIXELS "$"); // Iterator
+            psMetadataItem *item;
+            while ((item = psMetadataGetAndIncrement(iter))) {
+                psPixels *pixels = item->data.V; // Rejected pixels
+                psTrace("ppStack", 5, "Adding %ld rejected pixels to image %d", pixels->n, i);
+                rejects = psPixelsConcatenate(rejects, pixels);
+            }
+            psFree(iter);
+            psTrace("ppStack", 5, "%ld rejected pixels rejected from image %d", rejects->n, i);
+            psMetadataRemoveKey(ro->analysis, PPSTACK_REJECTED_PIXELS);
+            rejected->data[i] = rejects;
+        }
+
+
+        // Read convolutions by chunks
+        more = true;
+        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, 0) ||
+                    !pmReadoutReadChunkMask(readout, maskFits->data[i], 0, numScans, 0) ||
+                    !pmReadoutReadChunkWeight(readout, weightFits->data[i], 0, numScans, 0)) {
+                    psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i);
+                    psFree(readouts);
+                    psFree(rejected);
+                    psFree(outRO);
+                    psFree(view);
+                    return false;
+                }
+            }
+
+#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 (!ppStackReadoutFinal(config, outRO, readouts, rejected)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
+                psFree(readouts);
+                psFree(rejected);
+                psFree(outRO);
+                psFree(view);
+                return false;
+            }
+
+            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);
+
         for (int i = 0; i < num; i++) {
             psFitsClose(imageFits->data[i]);

trunk/ppStack/src/ppStackReadout.c

@@ -16 +16 @@
 #define COMBINED_FILES                  // Write combined images?
 
-static int sectionNum = 0;              // Section number; for debugging outputs
-
-
-bool ppStackReadout(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
-                    const psArray *regions, const psArray *kernels)
+
+
+bool ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
+                           const psArray *regions, const psArray *kernels)
 {
     assert(config);
@@ -29 +28 @@
     assert(readouts->n == regions->n);
     assert(regions->n == kernels->n);
+    static int sectionNum = 0;          // Section number; for debugging outputs
 
 
@@ -43 +43 @@
     int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
 
+
+    int num = readouts->n;              // Number of inputs
+    psArray *stack = psArrayAlloc(num); // Array for stacking
+
+    for (int i = 0; i < num; i++) {
+        pmReadout *ro = readouts->data[i];
+        assert(ro);
+        pmFPA *fpa = ro->parent->parent->parent; // Parent FPA
+
+        float weighting = psMetadataLookupF32(&mdok, fpa->analysis, "PPSTACK.WEIGHTING"); // Relative weight
+        if (!mdok || !isfinite(weighting)) {
+            psWarning("No WEIGHTING supplied for image %d --- set to unity.", i);
+            weighting = 1.0;
+        }
+
+        // Ensure there is a mask, or pmStackCombine will complain
+        if (!ro->mask) {
+            ro->mask = psImageAlloc(ro->image->numCols, ro->image->numRows, PS_TYPE_MASK);
+            psImageInit(ro->mask, 0);
+        }
+
+        stack->data[i] = pmStackDataAlloc(ro, weighting);
+    }
+
+    if (!pmStackCombine(outRO, stack, maskBad, maskBlank, kernelSize, iter, combineRej, useVariance, safe)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
+        psFree(stack);
+        return false;
+    }
+
+#ifdef COMBINED_FILES
+    {
+        psString name = NULL;           // Name of image
+        psStringAppend(&name, "combined_initial_%03d.fits", sectionNum);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
+        psFitsClose(fits);
+    }
+#endif
+#ifdef INSPECTION_FILES
+    for (int i = 0; i < stack->n; i++) {
+        pmStackData *data = stack->data[i]; // Data for this image
+        psImage *inspected = psPixelsToMask(NULL, data->pixels,
+                                            psRegionSet(0, outRO->image->numCols - 1,
+                                                        0, outRO->image->numRows - 1),
+                                            maskBlank);
+        psString name = NULL;           // Name of image
+        psStringAppend(&name, "inspect_%03d_%03d.fits", sectionNum, i);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, inspected, 0, NULL);
+        psFree(inspected);
+        psFitsClose(fits);
+    }
+#endif
+
+    // Reject pixels
+    for (int i = 0; i < num; i++) {
+        pmStackData *data = stack->data[i]; // Data for this image
+        pmReadout *readout = data->readout; // Readout of interest
+        int col0 = readout->col0, row0 = readout->row0; // Offset for readout
+        int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image
+
+        psRegion *valid = psRegionAlloc(col0, col0 + numCols, row0, row0 + numRows); // Valid region for rej
+        psPixels *reject = pmStackReject(data->pixels, valid, threshold, regions->data[i],
+                                         kernels->data[i]); // Pixels to reject
+        psFree(valid);
+
+        psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, PPSTACK_REJECTED_PIXELS,
+                         PS_DATA_PIXELS | PS_META_DUPLICATE_OK, "Rejected pixels from initial combination",
+                         reject);
+        psFree(reject);                 // Drop reference
+    }
+
+    psFree(stack);
+
+    sectionNum++;
+
+    return true;
+}
+
+
+
+bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
+                         const psArray *rejected)
+{
+    assert(config);
+    assert(outRO);
+    assert(readouts);
+    assert(rejected);
+    assert(readouts->n == rejected->n);
+
+    static int sectionNum = 0;          // Section number; for debugging outputs
+
+
+    // Get the recipe values
+    bool mdok;                          // Status of MD lookup
+    psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask
+    psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
+    bool useVariance = psMetadataLookupBool(&mdok, config->arguments, "VARIANCE"); // Use variance for rejection?
 
     int num = readouts->n;              // Number of inputs
@@ -110 +211 @@
         psListAdd(cellList, PS_LIST_TAIL, ro->parent);
 
-        stack->data[i] = pmStackDataAlloc(ro, weighting);
-    }
-
-    if (!pmStackCombine(outRO, stack, maskBad, maskBlank, kernelSize, iter, combineRej, useVariance, safe)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
-        psFree(fpaList);
-        psFree(cellList);
-        psFree(stack);
-        return false;
-    }
-
-#ifdef COMBINED_FILES
-    {
-        psString name = NULL;           // Name of image
-        psStringAppend(&name, "combined1_%03d.fits", sectionNum);
-        psFits *fits = psFitsOpen(name, "w");
-        psFree(name);
-        psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
-        psFitsClose(fits);
-    }
-#endif
-#ifdef INSPECTION_FILES
-    for (int i = 0; i < stack->n; i++) {
-        pmStackData *data = stack->data[i]; // Data for this image
-        psImage *inspected = psPixelsToMask(NULL, data->pixels,
-                                            psRegionSet(0, data->readout->image->numCols - 1,
-                                                        0, data->readout->image->numRows - 1),
-                                            maskBlank);
-        psString name = NULL;           // Name of image
-        psStringAppend(&name, "inspect_%03d_%03d.fits", sectionNum, i);
-        psFits *fits = psFitsOpen(name, "w");
-        psFree(name);
-        psFitsWriteImage(fits, NULL, inspected, 0, NULL);
-        psFree(inspected);
-        psFitsClose(fits);
-    }
-#endif
-
-    // Reject pixels
-    for (int i = 0; i < num; i++) {
-        pmStackData *data = stack->data[i]; // Data for this image
-        pmReadout *readout = data->readout; // Readout of interest
-        int col0 = readout->col0, row0 = readout->row0; // Offset for readout
-        int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image
-
-        psRegion *valid = psRegionAlloc(col0, col0 + numCols, row0, row0 + numRows); // Valid region for rej
-        psPixels *reject = pmStackReject(data->pixels, valid, threshold, regions->data[i],
-                                         kernels->data[i]); // Pixels to reject
-        psFree(valid);
-        psFree(data->pixels);
-        data->pixels = reject;
+        pmStackData *data = pmStackDataAlloc(ro, weighting);
+        data->pixels = psMemIncrRefCounter(rejected->data[i]);
+        stack->data[i] = data;
     }
 
 #ifdef REJECTION_FILES
-    for (int i = 0; i < stack->n; i++) {
-        pmStackData *data = stack->data[i]; // Data for this image
-        if (!data) {
-            continue;
-        }
-        psImage *rejected = psPixelsToMask(NULL, data->pixels,
-                                           psRegionSet(0, data->readout->image->numCols - 1,
-                                                       0, data->readout->image->numRows - 1),
-                                           maskBlank);
-        psString name = NULL;           // Name of image
-        psStringAppend(&name, "reject_%03d_%03d.fits", sectionNum, i);
-        psFits *fits = psFitsOpen(name, "w");
-        psFree(name);
-        psFitsWriteImage(fits, NULL, rejected, 0, NULL);
-        psFree(rejected);
-        psFitsClose(fits);
-    }
-#endif
-
-    if (!pmStackCombine(outRO, stack, maskBad, maskBlank, kernelSize, 0, combineRej, useVariance, false)) {
+    if (sectionNum == 0) {
+        for (int i = 0; i < stack->n; i++) {
+            pmStackData *data = stack->data[i]; // Data for this image
+            if (!data) {
+                continue;
+            }
+            psImage *reject = psPixelsToMask(NULL, data->pixels,
+                                             psRegionSet(0, outRO->image->numCols - 1,
+                                                         0, outRO->image->numRows - 1),
+                                             maskBlank);
+            psString name = NULL;           // Name of image
+            psStringAppend(&name, "reject_%03d.fits", i);
+            psFits *fits = psFitsOpen(name, "w");
+            psFree(name);
+            psFitsWriteImage(fits, NULL, reject, 0, NULL);
+            psFree(reject);
+            psFitsClose(fits);
+        }
+    }
+#endif
+
+    if (!pmStackCombine(outRO, stack, maskBad, maskBlank, 0, 0, NAN, useVariance, false)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
         psFree(fpaList);
@@ -194 +249 @@
     {
         psString name = NULL;           // Name of image
-        psStringAppend(&name, "combined2_%03d.fits", sectionNum);
+        psStringAppend(&name, "combined_final_%03d.fits", sectionNum);
         psFits *fits = psFitsOpen(name, "w");
         psFree(name);
@@ -219 +274 @@
     psFree(stack);
 
-    sectionNum++;
-
     return true;
 }