Changeset 14626 for trunk/ppStack/src/ppStackReadout.c

Timestamp:

Aug 22, 2007, 4:47:19 PM (19 years ago)

Author:

Paul Price

Message:

Changes to stacking to allow it to work on convolved images.

File:

• trunk/ppStack/src/ppStackReadout.c (modified) (10 diffs)

trunk/ppStack/src/ppStackReadout.c

@@ -15 +15 @@
 {
     // Get the recipe values
-    bool mdok;                          // Status of MD lookup
     int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Rejection iterations
     float combineRej = psMetadataLookupF32(NULL, config->arguments, "COMBINE.REJ"); // Combination threshold
-    float convolveRej = psMetadataLookupF32(NULL, config->arguments, "CONVOLVE.REJ"); // Convolution threshold
-    float extent = psMetadataLookupF32(NULL, config->arguments, "EXTENT"); // Extent of convolution
     psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask
     psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
-    psVector *seeing = psMetadataLookupPtr(&mdok, config->arguments, "SEEING"); // Seeing in each image
+    float threshold = psMetadataLookupF32(NULL, config->arguments, "THRESHOLD.MASK"); // Threshold for mask deconvolution
 
     // Get the output target
@@ -41 +38 @@
     int fileNum = 0;                    // Number of file
     float totExposure = 0.0;            // Total exposure time
-    pmReadout *templateRO = NULL;       // Template readout, for copy WCS
     psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
     psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging
@@ -53 +49 @@
 
         bool mdok;                      // Status of MD lookup
-        float seeing = psMetadataLookupF32(&mdok, inputFile->fpa->analysis, "PPSTACK.SEEING"); // Seeing FWHM
-        if (!mdok || !isfinite(seeing)) {
-            psWarning("No SEEING supplied for image %d --- set to unity.", fileNum);
-            seeing = 1.0;
-        }
         float weighting = psMetadataLookupF32(&mdok, inputFile->fpa->analysis,
                                               "PPSTACK.WEIGHTING"); // Relative weighting
@@ -69 +60 @@
             scale = 1.0;
         }
+
+        float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE"); // Exposure time
+#if 0
+        if (!isfinite(exposure)) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    "CELL.EXPOSURE is not set for input file %ld", stack->n);
+            psFree(stats);
+            psFree(rng);
+            psFree(fileIter);
+            psFree(fpaList);
+            psFree(cellList);
+            psFree(outRO);
+            psFree(stack);
+            return false;
+        }
+#endif
+        totExposure += exposure;        // Total exposure time
+        // Generate convolved version of input
+        pmReadout *convolved = pmReadoutAlloc(NULL); // Convolved version of input readout
+        if (!ppStackMatch(convolved, ro, config)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to match image %d.", fileNum);
+            psFree(stats);
+            psFree(rng);
+            psFree(fileIter);
+            psFree(fpaList);
+            psFree(cellList);
+            psFree(stack);
+            psFree(outRO);
+            psFree(convolved);
+            return false;
+        }
+
+#ifdef REJECTION_FILES
+        {
+            psString name = NULL;           // Name of image
+            psStringAppend(&name, "convolved%03d.fits", fileNum);
+            psFits *fits = psFitsOpen(name, "w");
+            psFree(name);
+            psFitsWriteImage(fits, NULL, convolved->image, 0, NULL);
+            psFitsClose(fits);
+        }
+#endif
+
+
+#if 0
+        // Background subtraction, scaling and normalisation is performed automatically by the image matching
 
         // Brain-dead background subtraction
@@ -80 +117 @@
             psFree(stack);
             psFree(outRO);
+            psFree(convolved);
             return false;
         }
@@ -99 +137 @@
             psFree(cellList);
             psFree(outRO);
+            psFree(convolved);
             psFree(stack);
             return false;
         }
-        totExposure += exposure;        // Total exposure time
+
         (void)psBinaryOp(ro->image, ro->image, "/", psScalarAlloc(exposure, PS_TYPE_F32));
         if (ro->weight) {
             (void)psBinaryOp(ro->weight, ro->weight, "/", psScalarAlloc(exposure, PS_TYPE_F32));
         }
-        pmStackData *data = pmStackDataAlloc(ro, seeing, weighting); // Data to stack
+#endif
+
+        if (fileNum == 0) {
+            // Copy astrometry over
+            pmFPA *fpa = ro->parent->parent->parent; // Template FPA
+            pmHDU *hdu = fpa->hdu; // Template HDU
+            pmHDU *outHDU = outFPA->hdu; // Output HDU
+            if (!outHDU || !hdu) {
+                psWarning("Unable to find HDU at FPA level to copy astrometry.");
+            } else {
+                if (!pmAstromReadWCS(outFPA, outCell->parent, hdu->header, 1.0)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry from input FPA.");
+                    psFree(stack);
+                    psFree(outRO);
+                    return false;
+                }
+                if (!outHDU->header) {
+                    outHDU->header = psMetadataAlloc();
+                }
+                if (!pmAstromWriteWCS(outHDU->header, outFPA, outCell->parent, WCS_TOLERANCE)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to output FPA.");
+                    psFree(stack);
+                    psFree(outRO);
+                    return false;
+                }
+            }
+        }
+
+        // Don't need the original any more!
+        psFree(ro);
+
+        // Ensure there is a mask, or pmStackCombine will complain
+        if (!convolved->mask) {
+            convolved->mask = psImageAlloc(convolved->image->numCols, convolved->image->numRows,
+                                           PS_TYPE_MASK);
+            psImageInit(convolved->mask, 0);
+        }
+
+        pmStackData *data = pmStackDataAlloc(convolved, weighting); // Data to stack
+        psFree(convolved);
         psArrayAdd(stack, ARRAY_BUFFER, data);
         psFree(data);                   // Drop reference
-
-        // 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);
-        }
-
-        if (fileNum == 0) {
-            templateRO = ro;
-        }
 
         fileNum++;
@@ -159 +227 @@
 #endif
 
-    // Only perform the additional rejection if we have seeing information
-    if (seeing && !pmStackReject(stack, maskBad, extent, convolveRej)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to reject input pixels.");
-        psFree(fpaList);
-        psFree(cellList);
-        psFree(stack);
-        psFree(outRO);
-        return false;
+    for (int i = 0; i < num; i++) {
+        pmStackData *data = stack->data[i]; // Data for this image
+        pmReadout *readout = data->readout; // Readout for this image
+
+        // Extract the regions and solutions used in the image matching
+        psArray *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions
+        {
+            psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD,
+                                                               "^SUBTRACTION.REGION$"); // Iterator
+            psMetadataItem *item = NULL;// Item from iteration
+            while ((item = psMetadataGetAndIncrement(iter))) {
+                assert(item->type == PS_DATA_REGION);
+                regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V);
+            }
+            psFree(iter);
+        }
+        psArray *solutions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of solutions
+        {
+            psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD,
+                                                               "^SUBTRACTION.SOLUTION$"); // Iterator
+            psMetadataItem *item = NULL;// Item from iteration
+            while ((item = psMetadataGetAndIncrement(iter))) {
+                assert(item->type == PS_DATA_VECTOR);
+                solutions = psArrayAdd(solutions, ARRAY_BUFFER, item->data.V);
+            }
+            psFree(iter);
+        }
+        assert(regions->n == solutions->n);
+
+        pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis,
+                                                            "SUBTRACTION.KERNEL"); // Kernels
+
+        psPixels *reject = pmSubtractionDeconvolveMask(data->pixels, threshold, regions,
+                                                       solutions, kernels); // List of pixels to reject
+        psFree(data->pixels);
+        data->pixels = reject;
     }
 
@@ -195 +291 @@
     }
 
+#if 0
     // Restore image to counts using the total exposure time
     (void)psBinaryOp(outRO->image, outRO->image, "*", psScalarAlloc(totExposure, PS_TYPE_F32));
@@ -200 +297 @@
         (void)psBinaryOp(outRO->weight, outRO->weight, "*", psScalarAlloc(totExposure, PS_TYPE_F32));
     }
+#endif
     psMetadataAddF32(outCell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE,
                      "Summed exposure time (sec)", totExposure);
@@ -215 +313 @@
     psFree(cellList);
 
-    // Copy astrometry over
-    pmFPA *templateFPA = templateRO->parent->parent->parent; // Template FPA
-    pmHDU *templateHDU = templateFPA->hdu; // Template HDU
-    pmHDU *outHDU = outFPA->hdu; // Output HDU
-    if (!outHDU || !templateHDU) {
-        psWarning("Unable to find HDU at FPA level to copy astrometry.");
-    } else {
-        if (!pmAstromReadWCS(outFPA, outCell->parent, templateHDU->header, 1.0)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry from input FPA.");
-            psFree(stack);
-            psFree(outRO);
-            return false;
-        }
-        if (!outHDU->header) {
-            outHDU->header = psMetadataAlloc();
-        }
-        if (!pmAstromWriteWCS(outHDU->header, outFPA, outCell->parent, WCS_TOLERANCE)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to output FPA.");
-            psFree(stack);
-            psFree(outRO);
-            return false;
-        }
-    }
-
     psFree(stack);
     psFree(outRO);                      // Drop reference