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

Timestamp:

Feb 22, 2008, 9:28:00 AM (18 years ago)

Author:

Paul Price

Message:

Merging in branch development. ppStack now works with incremental reads. Needs some cleanup, but it works.

File:

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

trunk/ppStack/src/ppStackReadout.c

@@ -10 +10 @@
 #include "ppStack.h"
 
-#define ARRAY_BUFFER 16                 // Number to add to array at a time
 #define WCS_TOLERANCE 0.001             // Tolerance for WCS
 
-//#define NO_CONVOLUTION                  // Don't perform convolutions?
-//#define CONVOLUTION_FILES               // Write convolutions?
 //#define REJECTION_FILES                 // Write rejection mask?
 //#define INSPECTION_FILES                // Write inspection mask?
 
-bool ppStackReadout(pmConfig *config, const pmFPAview *view)
+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)
 {
+    assert(config);
+    assert(outRO);
+    assert(readouts);
+    assert(regions);
+    assert(kernels);
+    assert(readouts->n == regions->n);
+    assert(regions->n == kernels->n);
+
+
     // 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
@@ -27 +36 @@
     psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
     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.MODEL"); // Model for PSF
-    int psfOrder = psMetadataLookupS32(&mdok, config->arguments, "PSF.ORDER"); // Spatial order for PSF
-
-    // Get the output target
-    pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
-    pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
-    pmFPA *outFPA = outCell->parent->parent; // Output FPA
-
-    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
-    assert(num > 0);
-
-    // Get the input sources
-    psArray *stack = psArrayAllocEmpty(ARRAY_BUFFER); // The stack of inputs
-    pmReadout *sources = pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES"); // Sources for stamps
-    psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD,
-                                                           "^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 num = readouts->n;              // Number of inputs
+    psArray *stack = psArrayAlloc(num); // Array for stacking
+
+    pmCell *outCell = outRO->parent;    // Output cell
+    pmChip *outChip = outCell->parent;  // Output chip
+    pmFPA *outFPA = outChip->parent;    // Output FPA
+
     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);
-        pmFPAfile *inputFile = fileItem->data.V; // An input file
-        pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Corresponding readout
+    for (int i = 0; i < num; i++) {
+        pmReadout *ro = readouts->data[i];
+        assert(ro);
+        pmFPA *fpa = ro->parent->parent->parent; // Parent FPA
 
         bool mdok;                      // Status of MD lookup
-        float weighting = psMetadataLookupF32(&mdok, inputFile->fpa->analysis,
-                                              "PPSTACK.WEIGHTING"); // Relative weighting
+        float weighting = psMetadataLookupF32(&mdok, fpa->analysis, "PPSTACK.WEIGHTING"); // Relative weight
         if (!mdok || !isfinite(weighting)) {
-            psWarning("No WEIGHTING supplied for image %d --- set to unity.", fileNum);
+            psWarning("No WEIGHTING supplied for image %d --- set to unity.", i);
             weighting = 1.0;
         }
-        float scale = psMetadataLookupF32(&mdok, inputFile->fpa->analysis, "PPSTACK.SCALE"); // Rel. scale
-        if (!mdok || !isfinite(scale)) {
-            psWarning("No SCALE supplied for image %d --- set to unity.", fileNum);
-            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;
-        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);
-            }
-#endif // CONVOLUTION_FILES
-
-        } 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) {
+#if 0
+        if (i == 0) {
             // Copy astrometry over
-            pmFPA *fpa = ro->parent->parent->parent; // Template FPA
             pmHDU *hdu = fpa->hdu; // Template HDU
             pmHDU *outHDU = outFPA->hdu; // Output HDU
@@ -224 +78 @@
                 psWarning("Unable to find HDU at FPA level to copy astrometry.");
             } else {
-                if (!pmAstromReadWCS(outFPA, outCell->parent, hdu->header, 1.0)) {
+                if (!pmAstromReadWCS(outFPA, outChip, hdu->header, 1.0)) {
                     psErrorClear();
                     psWarning("Unable to read WCS astrometry from input FPA.");
@@ -231 +85 @@
                         outHDU->header = psMetadataAlloc();
                     }
-                    if (!pmAstromWriteWCS(outHDU->header, outFPA, outCell->parent, WCS_TOLERANCE)) {
+                    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
                         psErrorClear();
                         psWarning("Unable to write WCS astrometry to output FPA.");
@@ -238 +92 @@
             }
         }
-
-        // Don't need the original any more!
-        psFree(ro);
+#endif
 
         // 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);
-        }
-
-        psListAdd(fpaList, PS_LIST_TAIL, inputFile->fpa);
+        if (!ro->mask) {
+            ro->mask = psImageAlloc(ro->image->numCols, ro->image->numRows, PS_TYPE_MASK);
+            psImageInit(ro->mask, 0);
+        }
+
+        psListAdd(fpaList, PS_LIST_TAIL, fpa);
         psListAdd(cellList, PS_LIST_TAIL, ro->parent);
 
-        pmStackData *data = pmStackDataAlloc(convolved, weighting); // Data to stack
-        psFree(convolved);
-        psArrayAdd(stack, ARRAY_BUFFER, data);
-        psFree(data);                   // Drop reference
-
-        fileNum++;
-    }
-    psFree(fileIter);
-    psFree(stats);
-    psFree(rng);
-
-    if (fileNum == 0) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Not enough good files to combine.");
-        psFree(fpaList);
-        psFree(cellList);
-        psFree(stack);
-        psFree(outRO);
-        return false;
+        stack->data[i] = pmStackDataAlloc(ro, weighting);
     }
 
@@ -277 +111 @@
         psFree(cellList);
         psFree(stack);
-        psFree(outRO);
         return false;
     }
@@ -283 +116 @@
 #ifdef INSPECTION_FILES
     {
-        psFits *fits = psFitsOpen("combined.fits", "w");
+        psString name = NULL;           // Name of image
+        psStringAppend(&name, "combined_%03d.fits", sectionNum);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
         psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
         psFitsClose(fits);
@@ -291 +127 @@
         pmStackData *data = stack->data[i]; // Data for this image
         psImage *inspected = psPixelsToMask(NULL, data->pixels,
-                                            psRegionSet(0, data->readout->image->numCols,
-                                                        0, data->readout->image->numRows),
+                                            psRegionSet(0, data->readout->image->numCols - 1,
+                                                        0, data->readout->image->numRows - 1),
                                             maskBlank);
         psString name = NULL;           // Name of image
-        psStringAppend(&name, "inspect%03d.fits", i);
+        psStringAppend(&name, "inspect_%03d_%03d.fits", sectionNum, i);
         psFits *fits = psFitsOpen(name, "w");
         psFree(name);
@@ -304 +140 @@
 #endif
 
-    for (int i = 0; i < stack->n; i++) {
+    // Reject pixels
+    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 = pmStackReject(data->pixels, threshold, regions,
-                                         solutions, kernels); // List of pixels to reject
+        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;
-
-        psFree(solutions);
-        psFree(regions);
     }
 
@@ -352 +162 @@
         }
         psImage *rejected = psPixelsToMask(NULL, data->pixels,
-                                           psRegionSet(0, data->readout->image->numCols,
-                                                       0, data->readout->image->numRows),
+                                           psRegionSet(0, data->readout->image->numCols - 1,
+                                                       0, data->readout->image->numRows - 1),
                                            maskBlank);
         psString name = NULL;           // Name of image
-        psStringAppend(&name, "reject%03d.fits", i);
+        psStringAppend(&name, "reject_%03d_%03d.fits", sectionNum, i);
         psFits *fits = psFitsOpen(name, "w");
         psFree(name);
@@ -370 +180 @@
         psFree(cellList);
         psFree(stack);
-        psFree(outRO);
         return false;
-    }
-
-    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));
-        }
     }
 
@@ -397 +198 @@
     psFree(cellList);
 
-    if (photometry) {
-        pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
-        pmFPACopy(photFile->fpa, outRO->parent->parent->parent);
-
-        if (!psphotReadout(config, view)) {
-            psWarning("Unable to perform photometry on stacked image.");
-            psErrorStackPrint(stderr, "Error stack from photometry:");
-            psErrorClear();
-        }
-
-        pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
-    }
-
     psFree(stack);
-    psFree(outRO);                      // Drop reference
+
+    sectionNum++;
 
     return true;