Changeset 14626

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.

Location:

trunk/ppStack/src

Files:

• Makefile.am (modified) (1 diff)
• ppStack.h (modified) (1 diff)
• ppStackArguments.c (modified) (5 diffs)
• ppStackCamera.c (modified) (2 diffs)
• ppStackMatch.c (added)
• ppStackReadout.c (modified) (10 diffs)

trunk/ppStack/src/Makefile.am

@@ -10 +10 @@
         ppStackLoop.c           \
         ppStackReadout.c        \
-        ppStackVersion.c
+        ppStackVersion.c        \
+        ppStackMatch.c
 
 noinst_HEADERS =                \

trunk/ppStack/src/ppStack.h

@@ -34 +34 @@
     );
 
+/// Convolve image to match specified seeing
+bool ppStackMatch(pmReadout *output,    ///< Convolved readout
+                  const pmReadout *input, ///< Readout to be convolved
+                  const pmConfig *config ///< Configuration
+    );
+
 #endif

trunk/ppStack/src/ppStackArguments.c

@@ -19 +19 @@
 {
     fprintf(stderr, "\nPan-STARRS Image combination\n\n");
-    fprintf(stderr, "Usage: %s INPUTS.mdc OUTPUT_ROOT\n"
+    fprintf(stderr, "Usage: %s INPUTS.mdc OUTPUT_ROOT -target FWHM -stamps FILENAME\n"
             "where INPUTS.mdc contains various METADATAs, each with:\n"
             "\tIMAGE(STR):     Image filename\n"
             "\tMASK(STR):      Mask filename\n"
             "\tWEIGHT(STR)     Weight map filename\n"
-            "\tSEEING(F32):    Seeing FWHM (pixels)\n"
             "\tWEIGHTING(F32): Relative weighting to be applied\n"
             "\tSCALE(F32):     Relative scaling to be applied\n",
@@ -103 +102 @@
 
     psMetadata *arguments = psMetadataAlloc(); // Command-line arguments
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-stamps", 0, "Stamps file with x,y,flux per line", NULL);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-target", 0, "Target PSF FWHM", NAN);
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL);
     psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Number of rejection iterations", 0);
@@ -110 +111 @@
     psMetadataAddU8(arguments,  PS_LIST_TAIL, "-mask-bad", 0, "Mask value for bad pixels", 0);
     psMetadataAddU8(arguments,  PS_LIST_TAIL, "-mask-blank", 0, "Mask value for blank region", 0);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-threshold-mask", 0, "Threshold for mask deconvolution", NAN);
 
     // XXX I want to get this from the recipe as well.  (same for everything else...)
@@ -117 +119 @@
         usage(argv[0], arguments, config);
     }
+
+    const char *stampsName = psMetadataLookupStr(NULL, arguments, "-stamps"); // Name of stamps file
+    if (!stampsName || strlen(stampsName) == 0) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Stamps file not specified with -stamps.");
+        goto ERROR;
+    }
+    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "STAMPS", 0, "Stamps file", stampsName);
+
+    float target = psMetadataLookupF32(NULL, arguments, "-target"); // Target PSF width
+    if (!isfinite(target)) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Target PSF width not specified with -target.");
+        goto ERROR;
+    }
+    psMetadataAddF32(config->arguments, PS_LIST_TAIL, "TARGET", 0, "Target PSF width", target);
+
 
     unsigned int numBad = 0;                     // Number of bad lines
@@ -138 +155 @@
     }
 
-    VALUE_ARG_RECIPE_INT("-iter",           "ITER",         S32, 0);
-    VALUE_ARG_RECIPE_FLOAT("-combine-rej",  "COMBINE.REJ",  F32);
-    VALUE_ARG_RECIPE_FLOAT("-convolve-rej", "CONVOLVE.REJ", F32);
-    VALUE_ARG_RECIPE_FLOAT("-extent",       "EXTENT",       F32);
-    VALUE_ARG_RECIPE_MASK("-mask-bad",      "MASK.BAD");
-    VALUE_ARG_RECIPE_MASK("-mask-blank",    "MASK.BLANK");
+    VALUE_ARG_RECIPE_INT("-iter",             "ITER",           S32, 0);
+    VALUE_ARG_RECIPE_FLOAT("-combine-rej",    "COMBINE.REJ",    F32);
+    VALUE_ARG_RECIPE_FLOAT("-convolve-rej",   "CONVOLVE.REJ",   F32);
+    VALUE_ARG_RECIPE_FLOAT("-extent",         "EXTENT",         F32);
+    VALUE_ARG_RECIPE_MASK("-mask-bad",        "MASK.BAD");
+    VALUE_ARG_RECIPE_MASK("-mask-blank",      "MASK.BLANK");
+    VALUE_ARG_RECIPE_FLOAT("-threshold-mask", "THRESHOLD.MASK", F32);
 
     psTrace("ppStack", 1, "Done reading command-line arguments\n");

trunk/ppStack/src/ppStackCamera.c

@@ -40 +40 @@
         psString weight = psMetadataLookupStr(&mdok, input, "WEIGHT"); // Name of weight map
 
-        float seeing = psMetadataLookupF32(&mdok, input, "SEEING"); // Seeing FWHM
-        if (!mdok || isnan(seeing) || seeing == 0.0) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks SEEING of type F32", item->name);
-            psFree(iter);
-            return false;
-        }
-
         float weighting = psMetadataLookupF32(&mdok, input, "WEIGHTING"); // Relative weighting
         if (!mdok || isnan(weighting) || weighting == 0.0) {
@@ -122 +115 @@
         }
 
-        psMetadataAddF32(imageFile->fpa->analysis, PS_LIST_TAIL, "PPSTACK.SEEING", 0,
-                         "Seeing for image", seeing);
         psMetadataAddF32(imageFile->fpa->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
                          "Relative weighting for image", weighting);

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