Changeset 13970

Timestamp:

Jun 25, 2007, 2:27:15 PM (19 years ago)

Author:

eugene

Message:

adding fringe residual measurement in addition to fringe amplitude; clean up ppImageLoop calls to stats and fringe

Location:

trunk/ppImage/src

Files:

• Makefile.am (modified) (2 diffs)
• ppImage.h (modified) (1 diff)
• ppImageDetrendFringe.c (modified) (13 diffs)
• ppImageDetrendFringe.h (modified) (3 diffs)
• ppImageDetrendReadout.c (modified) (1 diff)
• ppImageLoop.c (modified) (4 diffs)
• ppImageOptions.c (modified) (2 diffs)
• ppImageOptions.h (modified) (1 diff)
• ppImageStats.c (added)

trunk/ppImage/src/Makefile.am

@@ -26 +26 @@
         ppImageAstrom.c \
         ppImageAddstar.c \
+        ppImageStats.c \
         ppImageDefineFile.c \
         ppImageVersion.c
@@ -48 +49 @@
         ppImageAstrom.c \
         ppImageAddstar.c \
+        ppImageStats.c \
         ppImageDefineFile.c \
         ppImageVersion.c

trunk/ppImage/src/ppImage.h

@@ -65 +65 @@
 bool ppImageDefineFile (pmConfig *config, pmFPA *input, char *filerule, char *argname, pmFPAfileType fileType, pmDetrendType detrendType);
 
+// calculate stats, including MD5
+bool ppImageStats (pmConfig *config, pmChip *chip, const pmFPAview *inputView, const ppImageOptions *options);
+
+// write stats to output file 
+bool ppImageStatsOutput (pmConfig *config, pmFPA *fpa, const ppImageOptions *options);
+
 #endif

trunk/ppImage/src/ppImageDetrendFringe.c

@@ -10 +10 @@
 #include "ppImageDetrendFringe.h"
 
-bool ppImageDetrendFringeMeasure(pmReadout *readout, pmCell *fringe, const ppImageOptions *options)
+bool ppImageDetrendFringeMeasure(pmReadout *readout, pmCell *fringe, const bool isResidual, const ppImageOptions *options)
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);
@@ -49 +49 @@
     psBinaryOp(measurements->df, measurements->df, "*", psScalarAlloc(expTime, PS_TYPE_F32));
 
+    char *scienceFringes = NULL;
+    if (isResidual) {
+        scienceFringes = psStringCopy ("FRINGE.RESIDUALS");
+    } else {
+        scienceFringes = psStringCopy ("FRINGE.MEASUREMENTS");
+    }
+
     // Science fringe measurements
-    pmFringeStats *previous = psMetadataLookupPtr(NULL, readout->parent->analysis, "FRINGE.MEASUREMENTS");
+    pmFringeStats *previous = psMetadataLookupPtr(NULL, readout->parent->analysis, scienceFringes);
     if (previous) {
         // Multiple readouts: concatenate
@@ -73 +80 @@
     }
 
-    psMetadataAdd(readout->parent->analysis, PS_LIST_TAIL, "FRINGE.MEASUREMENTS",
+    psMetadataAdd(readout->parent->analysis, PS_LIST_TAIL, scienceFringes,
                   PS_DATA_UNKNOWN | PS_META_REPLACE, "Fringe measurements", measurements);
     psFree(measurements);
-
+    psFree(scienceFringes);
     psFree(references);
 
@@ -84 +91 @@
 
 // Pull the fringes out of the cell analysis FRINGE.MEASUREMENTS for a chip
-static psArray *getFringes(const pmChip *chip // Chip of interest
-    )
+// XXX need some error checks
+static psArray *getFringes(const pmChip *chip, const char *source)
 {
     psArray *cells = chip->cells;       // Component cells
@@ -91 +98 @@
     for (int i = 0; i < cells->n; i++) {
         pmCell *cell = cells->data[i];  // Cell of interest
-        fringes->data[i] = psMemIncrRefCounter(psMetadataLookupPtr(NULL, cell->analysis,
-                                                                   "FRINGE.MEASUREMENTS"));
+        fringes->data[i] = psMemIncrRefCounter(psMetadataLookupPtr(NULL, cell->analysis, source));
     }
 
@@ -101 +107 @@
 // Solve the fringe system: we have science fringe measurements for each cell, and an array of reference
 // fringe measurements for each cell.  Need to concatenate these together first, and then solve.
-bool ppImageDetrendFringeSolve(pmChip *scienceChip, const pmChip *refChip, const ppImageOptions *options)
+bool ppImageDetrendFringeSolve(pmChip *scienceChip, const pmChip *refChip, const bool isResidual, const ppImageOptions *options)
 {
     PS_ASSERT_PTR_NON_NULL(scienceChip, NULL);
@@ -107 +113 @@
     PS_ASSERT_PTR_NON_NULL(options, NULL);
 
-    psArray *science = getFringes(scienceChip); // Fringe measurements on science chip
+    psArray *science = NULL;
+    if (isResidual) {
+        science = getFringes(scienceChip, "FRINGE.RESIDUALS"); // Fringe residuals on science chip
+    } else {
+        science = getFringes(scienceChip, "FRINGE.MEASUREMENTS"); // Fringe measurements on science chip
+    }
+
     pmFringeStats *scienceCat = pmFringeStatsConcatenate(science, NULL, NULL); // Science fringes
     psFree(science);
@@ -113 +125 @@
     // Need to transform the array of cells each with an array of fringes --> array of fringes for the chip as
     // a whole
-    psArray *references = getFringes(refChip); // Fringe measurements on reference chip
+    psArray *references = getFringes(refChip, "FRINGE.MEASUREMENTS"); // Fringe measurements on reference chip
     int numRefs = ((psArray*)references->data[0])->n; // Number of reference fringes
     psArray *referencesCat = psArrayAlloc(numRefs); // Reference fringes
@@ -132 +144 @@
                                                    options->fringeIter, options->fringeKeep);
 
-    psMetadataAdd(scienceChip->analysis, PS_LIST_TAIL, "FRINGE.SOLUTION", PS_DATA_UNKNOWN,
-                  "Fringe solution", solution);
+    if (isResidual) {
+        psMetadataAdd(scienceChip->analysis, PS_LIST_TAIL, "FRINGE.RESIDUAL.SOLUTION", PS_DATA_UNKNOWN, "Fringe solution", solution);
+    } else {
+        psMetadataAdd(scienceChip->analysis, PS_LIST_TAIL, "FRINGE.SOLUTION", PS_DATA_UNKNOWN, "Fringe solution", solution);
+    }
+
+# if (0)
+    // write the fringe amplitude or residual amplitude to the header
+    // XXX this is measured per cell, but we only have headers per chip
+    pmHDU *hdu = pmHDUFromCell(science);// HDU  of interest
+    for (int i = 0; i < solution->nFringeFrames; i++) {
+        // write metadata header value
+        psString keyword = NULL;
+        if (isResidual) {
+            psStringAppend (&keyword, "FRES_%02dV", i);
+        } else {
+            psStringAppend (&keyword, "FRNG_%02dV", i);
+        }
+        psMetadataAddF32(hdu->header, PS_LIST_TAIL, keyword, PS_META_REPLACE, "Fringe Amplitude", solution->coeff->data.F32[i + 1]);
+        psFree (keyword);
+
+        keyword = NULL;
+        if (isResidual) {
+            psStringAppend (&keyword, "FRES_%02dE", i);
+        } else {
+            psStringAppend (&keyword, "FRNG_%02dE", i);
+        }
+        psMetadataAddF32(hdu->header, PS_LIST_TAIL, keyword, PS_META_REPLACE, "Fringe Amplitude error", solution->coeffErr->data.F32[i + 1]);
+        psFree (keyword);
+    }
+# endif
+
     psFree(solution);
     psFree(scienceCat);
@@ -142 +184 @@
 
 
-bool ppImageDetrendFringeGenerate(pmCell *science, pmCell *fringes)
+bool ppImageDetrendFringeGenerate(pmCell *science, pmCell *fringes, const ppImageOptions *options)
 {
     PS_ASSERT_PTR_NON_NULL(science, false);
@@ -191 +233 @@
         psFree(md5string);
 
+# if (0)
         // write metadata header value
+        // XXX this is measured per cell, but we only have headers per chip
         psString keyword = NULL;
-        psStringAppend (&keyword, "FRNG_%02d", i);
+        psStringAppend (&keyword, "FRNG_%02dV", i);
         psMetadataAddF32(hdu->header, PS_LIST_TAIL, keyword, PS_META_REPLACE, "Fringe Amplitude", solution->coeff->data.F32[i + 1]);
         psFree (keyword);
 
         keyword = NULL;
-        psStringAppend (&keyword, "FRNG_%02dD", i);
+        psStringAppend (&keyword, "FRNG_%02dE", i);
         psMetadataAddF32(hdu->header, PS_LIST_TAIL, keyword, PS_META_REPLACE, "Fringe Amplitude error", solution->coeffErr->data.F32[i + 1]);
         psFree (keyword);
+# endif
+
     }
     if (expTime != 1.0) {
@@ -219 +265 @@
             return false;
         }
+
+        // measure residual fringe amplitude. results go to FRINGE.RESIDUALS
+        ppImageDetrendFringeMeasure (readout, fringes, true, options);
     }
     psFree(sumFringe);
@@ -233 +282 @@
 }
 
+bool ppImageDetrendFringeApply (pmConfig *config, pmChip *chip, const pmFPAview *inputView, const ppImageOptions *options) {
+
+    pmCell *cell = NULL;
+
+    assert (options->doFringe); // do not call if not needed
+    assert (inputView->chip != -1);
+    assert (inputView->cell == -1);
+    assert (inputView->readout == -1);
+
+    pmFPAview *view = pmFPAviewAlloc(0); // View for local processing
+    *view = *inputView;
+
+    // select the reference chip
+    pmChip *fringe = pmFPAfileThisChip(config->files, view, "PPIMAGE.FRINGE");
+    if (!fringe) {
+        psError(PS_ERR_UNKNOWN, false, "missing fringe reference data.\n");
+        psFree (view);
+        return false;
+    }
+
+    // Solve the fringe system
+    if (!ppImageDetrendFringeSolve(chip, fringe, false, options)) {
+        psError(PS_ERR_UNKNOWN, false, "failed to solve the fringe system.\n");
+        psFree (view);
+        return false;
+    }
+
+    // Go back over the cells to apply the fringe correction
+    view->cell = view->readout = -1;
+    while ((cell = pmFPAviewNextCell(view, chip->parent, 1)) != NULL) {
+        if (!cell->process || !cell->file_exists) {
+            continue;
+        }
+
+        // Apply the fringe correction
+        psTrace("ppImage", 3, "Applying fringe correction...\n");
+        pmCell *fringeCell = pmFPAfileThisCell(config->files, view, "PPIMAGE.FRINGE");
+        if (!fringeCell) {
+            psError(PS_ERR_UNKNOWN, false, "missing fringe reference data.\n");
+            psFree (view);
+            return false;
+        }
+
+        if (!ppImageDetrendFringeGenerate(cell, fringeCell, options)) {
+            psError(PS_ERR_UNKNOWN, false, "failed to apply fringe image.\n");
+            psFree (view);
+            return false;
+        }
+    }
+
+    // Solve the residual fringe system
+    if (!ppImageDetrendFringeSolve(chip, fringe, true, options)) {
+        psError(PS_ERR_UNKNOWN, false, "failed to solve the residual fringe system.\n");
+        psFree (view);
+        return false;
+    }
+
+    psFree (view);
+    return true;
+}

trunk/ppImage/src/ppImageDetrendFringe.h

@@ -9 +9 @@
 bool ppImageDetrendFringeMeasure(pmReadout *readout, // Readout to measure
                                  pmCell *fringe, // Fringe cell (each readout is a different component)
+                                 const bool isResidual,
                                  const ppImageOptions *options // Options
     );
@@ -15 +16 @@
 bool ppImageDetrendFringeSolve(pmChip *scienceChip, // Chip with science
                                const pmChip *refChip, // Chip with reference fringes
+                               const bool isResidual,
                                const ppImageOptions *options // Options
     );
@@ -20 +22 @@
 // Generate fringe frame
 bool ppImageDetrendFringeGenerate(pmCell *science, // Science cell
-                                  pmCell *fringes // Fringe cell, one readout per fringe component
+                                  pmCell *fringes, // Fringe cell, one readout per fringe component
+                                  const ppImageOptions *options // Options
     );
 
 
+bool ppImageDetrendFringeApply (pmConfig *config, // config
+                                pmChip *chip, // science chip
+                                const pmFPAview *inputView, // current view
+                                const ppImageOptions *options // options
+);
+
 #endif

trunk/ppImage/src/ppImageDetrendReadout.c

@@ -88 +88 @@
     if (options->doFringe) {
         pmCell *fringe = pmFPAfileThisCell(config->files, detview, "PPIMAGE.FRINGE");
-        if (!ppImageDetrendFringeMeasure(input, fringe, options)) {
+        if (!ppImageDetrendFringeMeasure(input, fringe, false, options)) {
             return false;
         }

trunk/ppImage/src/ppImageLoop.c

@@ -3 +3 @@
 #endif
 
-#include <ppStats.h>
 #include "ppImage.h"
 #include "ppImageDetrendFringe.h"
@@ -10 +9 @@
 bool ppImageLoop (pmConfig *config, ppImageOptions *options) {
 
-    bool mdok;                      // Status of MD lookup
     bool status;
     pmChip *chip;
     pmCell *cell;
     pmReadout *readout;
-    psMetadata *stats = NULL;
-
-    const char *statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS"); // Filename for statistics
-    FILE *statsFile = NULL;             // File stream for statistics
-    if (mdok && statsName && strlen(statsName) > 0) {
-        psString resolved = pmConfigConvertFilename(statsName, config, true); // Resolved filename
-        statsFile = fopen(resolved, "w");
-        if (!statsFile) {
-            psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
-            psFree(resolved);
-            return false;
-        }
-        stats = psMetadataAlloc();
-        psFree(resolved);
-    }
 
     pmFPAfile *input = psMetadataLookupPtr(&status, config->files, "PPIMAGE.INPUT");
@@ -96 +79 @@
         }
 
-        // Solve the fringe system
+        // Apply the fringe correction
         if (options->doFringe) {
-            pmChip *fringe = pmFPAfileThisChip(config->files, view, "PPIMAGE.FRINGE");
-            if (!ppImageDetrendFringeSolve(chip, fringe, options)) {
+            if (!ppImageDetrendFringeApply (config, chip, view, options)) {
                 psFree (view);
                 return false;
             }
-        }
-
-        // Go back over the cells to apply the fringe and do statistics
-        view->cell = view->readout = -1;
-        while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
-            if (!cell->process || !cell->file_exists) {
-                continue;
-            }
-
-            // Apply the fringe correction
-            if (options->doFringe) {
-                psTrace("ppImage", 3, "Applying fringe correction...\n");
-                pmCell *fringeCell = pmFPAfileThisCell(config->files, view, "PPIMAGE.FRINGE");
-                if (!ppImageDetrendFringeGenerate(cell, fringeCell)) {
-                    psFree (view);
-                    return false;
-                }
-            }
-
-            // Perform statistics on the detrended cell
-            if (stats) {
-                bool mdok;              // Status of MD lookup
-                pmFPAfile *output = psMetadataLookupPtr(&mdok, config->files, "PPIMAGE.OUTPUT");
-                if (!mdok || !output) {
-                    psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find file PPIMAGE.OUTPUT.");
-                    psFree (view);
-                    return false;
-                }
-
-                if (!ppStats(stats, output->fpa, view,
-                             options->satMask | options->badMask | options->maskValue,
-                             config)) {
-                    psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate stats for image.");
-                    psFree(stats);
-                    psFree(view);
-                    return false;
-                }
-
-                if (options->doFringe && !ppStatsFringe(stats, chip, "FRINGE.SOLUTION")) {
-                    psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to extract fringe solution for image.");
-                    psFree(stats);
-                    psFree(view);
-                    return false;
-                }
-            }
-
-            // Add MD5 information for cell
-            pmHDU *hdu = pmHDUFromCell(cell); // HDU that owns the cell
-            while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) {
-                const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
-                const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME");
-
-                psString headerName = NULL; // Header name for MD5
-                psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
-
-                psVector *md5 = psImageMD5(readout->image); // md5 hash
-                psString md5string = psMD5toString(md5); // String
-                psFree(md5);
-                psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE,
-                                 "Image MD5", md5string);
-                psFree(md5string);
-                psFree(headerName);
-            }
-        }
+        }
+        
+        // measure various statistics for this image
+        if (!ppImageStats (config, chip, view, options)) {
+            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to measures stats for image");
+            psFree (view);
+            return false;
+        }
 
         if (!ppImageMosaicChip(config, options, view, "PPIMAGE.CHIP", "PPIMAGE.OUTPUT")) {
@@ -214 +140 @@
 
     // Write out summary statistics
-    if (stats) {
-        char *statsMDC = psMetadataConfigFormat(stats);
-        if (!statsMDC || strlen(statsMDC) == 0) {
-            psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n");
-        } else {
-            fprintf(statsFile, "%s", statsMDC);
-        }
-        psFree(statsMDC);
-        fclose(statsFile);
-        psFree(stats);
+    if (!ppImageStatsOutput (config, input->fpa, options)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to write statistics file.\n");
+        psFree (view);
+        return false;
     }
 

trunk/ppImage/src/ppImageOptions.c

@@ -39 +39 @@
     options->doAstromChip = false;      // Astrometry
     options->doAstromMosaic = false;    // Astrometry
+
+    options->doStats    = false;        // Measure and save image statistics
 
     options->BaseFITS   = false;        // create binned image (scale 1)
@@ -187 +189 @@
     options->doShutter = psMetadataLookupBool(NULL, recipe, "SHUTTER");
 
+    options->doStats = false;
+    char *statsName = psMetadataLookupStr(&status, config->arguments, "STATS"); // Filename for statistics
+    if (statsName) {
+        options->doStats = true;
+    }
+
     // binned image options
     options->xBin1 = psMetadataLookupS32(&status, recipe, "BIN1.XBIN");

trunk/ppImage/src/ppImageOptions.h

@@ -29 +29 @@
     void *nonLinearSource;
 
+    bool doStats;
+
     bool BaseFITS;
     bool ChipFITS;