Changeset 15844

Timestamp:

Dec 14, 2007, 3:32:16 PM (18 years ago)

Author:

Paul Price

Message:

Using pmPSFEnvelope to generate a target PSF.

Location:

trunk/ppStack/src

Files:

• ppStack.h (modified) (1 diff)
• ppStackArguments.c (modified) (3 diffs)
• ppStackCamera.c (modified) (5 diffs)
• ppStackMatch.c (modified) (3 diffs)
• ppStackReadout.c (modified) (4 diffs)

trunk/ppStack/src/ppStack.h

@@ -38 +38 @@
                   const pmReadout *input, ///< Readout to be convolved
                   const pmReadout *sourcesRO, ///< Readout with sources
+                  const pmPSF *psf,     ///< Target PSF
                   const pmConfig *config ///< Configuration
     );

trunk/ppStack/src/ppStackArguments.c

@@ -25 +25 @@
             "\tMASK(STR):      Mask filename\n"
             "\tWEIGHT(STR)     Weight map filename\n"
+            "\tPSF(STR)        PSF filename\n"
             "\tWEIGHTING(F32): Relative weighting to be applied\n"
             "\tSCALE(F32):     Relative scaling to be applied\n",
@@ -117 +118 @@
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-threshold-mask", 0, "Threshold for mask deconvolution", NAN);
     psMetadataAddBool(arguments, PS_LIST_TAIL, "-photometry", 0, "Do photometry on stacked image?", false);
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-psf-instances", 0, "Number of instances for PSF generation", 5);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-psf-radius", 0, "Radius for PSF generation", 20.0);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-psf-model", 0, "Model name for PSF generation", "PS_MODEL_RGAUSS");
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-psf-order", 0, "Spatial order for PSF generation", 3);
 
     if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 3) {
@@ -165 +170 @@
     VALUE_ARG_RECIPE_FLOAT("-threshold-mask", "THRESHOLD.MASK", F32);
 
+    VALUE_ARG_RECIPE_INT("-psf-instances", "PSF.INSTANCES", S32, 0);
+    VALUE_ARG_RECIPE_FLOAT("-psf-radius",  "PSF.RADIUS",    F32);
+    VALUE_ARG_RECIPE_INT("-psf-order",     "PSF.ORDER",     S32, 0);
+    valueArgStr(config, arguments, "-psf-model", "PSF.MODEL", config->arguments);
+
     if (psMetadataLookupBool(NULL, arguments, "-photometry") ||
         psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {

trunk/ppStack/src/ppStackCamera.c

@@ -14 +14 @@
 {
     bool haveWeights = false;           // Do we have weight maps?
+    bool havePSFs = false;              // Do we have PSFs?
 
     psMetadata *inputs = psMetadataLookupMetadata(NULL, config->arguments, "INPUTS"); // The inputs info
@@ -39 +40 @@
         psString mask = psMetadataLookupStr(&mdok, input, "MASK"); // Name of mask
         psString weight = psMetadataLookupStr(&mdok, input, "WEIGHT"); // Name of weight map
+        psString psf = psMetadataLookupStr(&mdok, input, "PSF"); // Name of PSF
 
         float weighting = psMetadataLookupF32(&mdok, input, "WEIGHTING"); // Relative weighting
@@ -115 +117 @@
         }
 
+        // Optionally add the psf file
+        if (psf && strlen(psf) > 0) {
+            havePSFs = true;
+            psArray *psfFiles = psArrayAlloc(1); // Array of filenames for this FPA
+            psfFiles->data[0] = psMemIncrRefCounter(psf);
+            psMetadataAddArray(config->arguments, PS_LIST_TAIL, "PSF.FILENAMES", PS_META_REPLACE,
+                               "Filenames of PSF files", psfFiles);
+            psFree(psfFiles);
+
+            bool status;
+            pmFPAfile *psfFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PSPHOT.PSF.LOAD",
+                                                       "PSF.FILENAMES");
+            if (!status) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf);
+                return false;
+            }
+            if (psfFile->type != PM_FPA_FILE_PSF) {
+                psError(PS_ERR_IO, true, "PSPHOT.PSF.LOAD is not of type PSF");
+                return false;
+            }
+        }
+
         psMetadataAddF32(imageFile->fpa->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
                          "Relative weighting for image", weighting);
@@ -129 +153 @@
     if (psMetadataLookup(config->arguments, "WEIGHT.FILENAMES")) {
         psMetadataRemoveKey(config->arguments, "WEIGHT.FILENAMES");
+    }
+    if (psMetadataLookup(config->arguments, "PSF.FILENAMES")) {
+        psMetadataRemoveKey(config->arguments, "PSF.FILENAMES");
     }
 
@@ -214 +241 @@
     }
 
+    // Output PSF
+    if (havePSFs) {
+        pmFPAfile *outPSF = pmFPAfileDefineOutput(config, output->fpa, "PSPHOT.PSF.SAVE");
+        if (!outPSF) {
+            psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.PSF.SAVE"));
+            return false;
+        }
+        if (outPSF->type != PM_FPA_FILE_PSF) {
+            psError(PS_ERR_IO, true, "PSPHOT.PSF.SAVE is not of type PSF");
+            return false;
+        }
+        outPSF->save = true;
+    }
+
     // Sources for use as stamps
     bool status = false;                // Found the file?

trunk/ppStack/src/ppStackMatch.c

@@ -12 +12 @@
 
 bool ppStackMatch(pmReadout *output, const pmReadout *input, const pmReadout *sourcesRO,
-                  const pmConfig *config)
+                  const pmPSF *psf, const pmConfig *config)
 {
     // Look up appropriate values from the ppSub recipe
@@ -46 +46 @@
     // These values are specified specifically for stacking
     const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS"); // Filename for stamps
-    float target = psMetadataLookupF32(NULL, config->arguments, "TARGET"); // Target PSF width
 
     psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
@@ -66 +65 @@
         }
     }
-    pmReadout *fake = pmReadoutFakeFromSources(input->image->numCols, input->image->numRows,
-                                               sources, target, powf(10.0, -0.4 * maxMag));
+
+    pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
+
+    if (!pmReadoutFakeFromSources(fake, input->image->numCols, input->image->numRows, sources, NULL, NULL,
+                                  psf, powf(10.0, -0.4 * maxMag), 0, false)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
+        psFree(fake);
+        psFree(optWidths);
+        return false;
+    }
 
 #ifdef TESTING

trunk/ppStack/src/ppStackReadout.c

@@ -28 +28 @@
     float threshold = psMetadataLookupF32(NULL, config->arguments, "THRESHOLD.MASK"); // Threshold for mask deconvolution
     bool photometry = psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY"); // Perform photometry?
+    int psfInstances = psMetadataLookupS32(&mdok, config->arguments, "PSF.INSTANCES"); // Number of instances for PSF
+    float psfRadius = psMetadataLookupF32(NULL, config->arguments, "PSF.RADIUS"); // Radius for PSF
+    const char *psfModel = psMetadataLookupStr(NULL, config->arguments, "PSF.MODE"); // Model for PSF
+    int psfOrder = psMetadataLookupS32(&mdok, config->arguments, "PSF.ORDER"); // Spatial order for PSF
 
     // Get the output target
@@ -43 +47 @@
                                                            "^PPSTACK.INPUT$"); // Iterator over input files
     psMetadataItem *fileItem;           // Item from iteration
+    int fileNum = 0;                    // Number of file
+    psList *psfList = psListAlloc(NULL); // List of PSFs for PSF envelope
+    pmPSF *outPSF = NULL;               // Ouptut PSF
+    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
+        pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Corresponding readout
+        pmCell *cell = ro->parent;      // The cell
+        pmChip *chip = cell->parent;    // The chip: holds the PSF
+
+        bool mdok;                      // Status of MD lookup
+        pmPSF *psf = psMetadataLookupPtr(&mdok, chip->analysis, "PSPHOT.PSF");
+        if (mdok && psf) {
+            psListAdd(psfList, PS_LIST_TAIL, psf);
+            if (numCols == 0 && numRows == 0) {
+                numCols = ro->image->numCols;
+                numRows = ro->image->numRows;
+            }
+        }
+    }
+
+    bool convolve = false;              // Convolve the input images?
+    if (psfList->n > 0) {
+        psArray *psfArray = psListToArray(psfList); // Array of PSFs
+        outPSF = pmPSFEnvelope(numCols, numRows, psfArray, psfInstances, psfRadius, psfModel,
+                               psfOrder, psfOrder);
+        psFree(psfArray);
+        if (!outPSF) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
+            // XXX Free stuff
+            return false;
+        }
+        convolve = true;
+        PS_ASSERT_PTR_NON_NULL(sources, false);
+    }
+
+    // Iterate through again to get the convolved images (or not)
+
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics
     psRandom *rng = psRandomAlloc(0, PS_RANDOM_TAUS); // Random number generator
-    int fileNum = 0;                    // Number of file
     float totExposure = 0.0;            // Total exposure time
     psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
     psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging
+
+    psMetadataIteratorSet(fileIter, PS_LIST_HEAD);
+    fileNum = 0;
     while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
         assert(fileItem->type == PS_DATA_UNKNOWN);
@@ -85 +130 @@
 
         // Generate convolved version of input
-        pmReadout *convolved = pmReadoutAlloc(NULL); // Convolved version of input readout
-#ifndef NO_CONVOLUTION
-        if (!ppStackMatch(convolved, ro, sources, config)) {
-            psWarning("Unable to match image %d --- ignoring.", fileNum);
-            psErrorClear();
-            psFree(convolved);
-            // XXX Free the bad image so it's not taking up good memory!
-            continue;
-        }
+        pmReadout *convolved;
+        if (convolve) {
+            convolved = pmReadoutAlloc(NULL); // Convolved version of input readout
+            // Background subtraction, scaling and normalisation is performed automatically by the image
+            // matching
+            if (!ppStackMatch(convolved, ro, sources, outPSF, config)) {
+                psWarning("Unable to match image %d --- ignoring.", fileNum);
+                psErrorClear();
+                psFree(convolved);
+                // XXX Free the bad image so it's not taking up good memory!
+                continue;
+            }
 
 #ifdef CONVOLUTION_FILES
-        if (convolved->image) {
-            psString name = NULL;           // Name of image
-            psStringAppend(&name, "convolved%03d_image.fits", fileNum);
-            psFits *fits = psFitsOpen(name, "w");
-            psFree(name);
-            psFitsWriteImage(fits, NULL, convolved->image, 0, NULL);
-            psFitsClose(fits);
-        }
-        if (convolved->mask) {
-            psString name = NULL;           // Name of image
-            psStringAppend(&name, "convolved%03d_mask.fits", fileNum);
-            psFits *fits = psFitsOpen(name, "w");
-            psFree(name);
-            psFitsWriteImage(fits, NULL, convolved->mask, 0, NULL);
-            psFitsClose(fits);
-        }
-        if (convolved->weight) {
-            psString name = NULL;           // Name of image
-            psStringAppend(&name, "convolved%03d_weight.fits", fileNum);
-            psFits *fits = psFitsOpen(name, "w");
-            psFree(name);
-            psFitsWriteImage(fits, NULL, convolved->weight, 0, NULL);
-            psFitsClose(fits);
-        }
+            if (convolved->image) {
+                psString name = NULL;           // Name of image
+                psStringAppend(&name, "convolved%03d_image.fits", fileNum);
+                psFits *fits = psFitsOpen(name, "w");
+                psFree(name);
+                psFitsWriteImage(fits, NULL, convolved->image, 0, NULL);
+                psFitsClose(fits);
+            }
+            if (convolved->mask) {
+                psString name = NULL;           // Name of image
+                psStringAppend(&name, "convolved%03d_mask.fits", fileNum);
+                psFits *fits = psFitsOpen(name, "w");
+                psFree(name);
+                psFitsWriteImage(fits, NULL, convolved->mask, 0, NULL);
+                psFitsClose(fits);
+            }
+            if (convolved->weight) {
+                psString name = NULL;           // Name of image
+                psStringAppend(&name, "convolved%03d_weight.fits", fileNum);
+                psFits *fits = psFitsOpen(name, "w");
+                psFree(name);
+                psFitsWriteImage(fits, NULL, convolved->weight, 0, NULL);
+                psFitsClose(fits);
+            }
 #endif // CONVOLUTION_FILES
 
-#else
-        convolved = psMemIncrRefCounter(ro);
-        sources = NULL;
-#endif // NO_CONVOLVUTION
-
-
-#if 0
-        // Background subtraction, scaling and normalisation is performed automatically by the image matching
-
-        // Brain-dead background subtraction
-        if (!psImageBackground(stats, ro->image, ro->mask, maskBad, rng)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to get image background on image %ld", stack->n);
-            psFree(stats);
-            psFree(rng);
-            psFree(fileIter);
-            psFree(fpaList);
-            psFree(cellList);
-            psFree(stack);
-            psFree(outRO);
-            psFree(convolved);
-            return false;
-        }
-        psTrace("ppStack", 3, "Background for image %d is %f\n", fileNum, stats->robustMedian);
-        (void)psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(stats->robustMedian, PS_TYPE_F32));
-
-        // Apply scaling
-        (void)psBinaryOp(ro->image, ro->image, "*", psScalarAlloc(1.0 / scale, PS_TYPE_F32));
-
-        // Normalise each input by the exposure time
-        float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE"); // Exposure time
-        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(convolved);
-            psFree(stack);
-            return false;
-        }
-
-        (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));
-        }
-#endif
+        } else {
+            // No convolution --- just use the unconvolved images as the "convolved" images, with some tweaks
+            convolved = psMemIncrRefCounter(ro);
+            sources = NULL;
+
+            // Brain-dead background subtraction
+            if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskBad, rng)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to get image background on image %ld", stack->n);
+                psFree(stats);
+                psFree(rng);
+                psFree(fileIter);
+                psFree(fpaList);
+                psFree(cellList);
+                psFree(stack);
+                psFree(outRO);
+                psFree(convolved);
+                return false;
+            }
+            psTrace("ppStack", 3, "Background for image %d is %f\n", fileNum, stats->robustMedian);
+            (void)psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(stats->robustMedian, PS_TYPE_F32));
+
+            // Apply scaling
+            (void)psBinaryOp(ro->image, ro->image, "*", psScalarAlloc(1.0 / scale, PS_TYPE_F32));
+
+            // Normalise each input by the exposure time
+            float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE");// Exposure time
+            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(convolved);
+                psFree(stack);
+                return false;
+            }
+
+            (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));
+            }
+        }
 
         if (fileNum == 0) {
@@ -330 +374 @@
     }
 
-#if 0
-    // Restore image to counts using the total exposure time
-    (void)psBinaryOp(outRO->image, outRO->image, "*", psScalarAlloc(totExposure, PS_TYPE_F32));
-    if (outRO->weight) {
-        (void)psBinaryOp(outRO->weight, outRO->weight, "*", psScalarAlloc(totExposure, PS_TYPE_F32));
-    }
-#endif
+    if (!convolve) {
+        // Restore image to counts using the total exposure time
+        (void)psBinaryOp(outRO->image, outRO->image, "*", psScalarAlloc(totExposure, PS_TYPE_F32));
+        if (outRO->weight) {
+            (void)psBinaryOp(outRO->weight, outRO->weight, "*", psScalarAlloc(totExposure, PS_TYPE_F32));
+        }
+    }
+
     psMetadataAddF32(outCell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE,
                      "Summed exposure time (sec)", totExposure);