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

Timestamp:

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

Author:

Paul Price

Message:

Using pmPSFEnvelope to generate a target PSF.

File:

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

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);