Changeset 7061 for trunk/ppMerge/src/ppMergeBackground.c

Timestamp:

May 3, 2006, 5:57:35 PM (20 years ago)

Author:

Paul Price

Message:

Code coming together, doesn't yet compile

Location:

trunk/ppMerge/src

Files:

• . (modified) (1 prop)
• ppMergeBackground.c (modified) (3 diffs)

trunk/ppMerge/src/ppMergeBackground.c

@@ -5 +5 @@
 
 #include "ppMergeBackground.h"
+
+#define MAX_ITER 10                     // Maximum number of iterations for pmFlatNormalize
+#define TOLERANCE 1e-3                  // Tolerance to reach for pmFlatNormalize
+
 
 // Extract a particular statistic from the stats
@@ -29 +33 @@
 
 
-// Get the background for each chip of each input
-psImage *ppMergeBackground(const ppMergeOptions *options, // The options
-                           const pmConfig *config // The configuration
+// Get the scale and zero for each chip of each input
+bool ppMergeScaleZero(ppImage **scales, // The scales for each integration/cell
+                      ppImage **zeros, // The zeroes for each integration/cell
+                      ppMergeData *data,// The data
+                      const ppMergeOptions *options, // The options
+                      const pmConfig *config // The configuration
     )
 {
     assert(config->camera);             // Need the camera configuration
-    int numCells = 0;                   // Number of cells in the model FPA
-    {
-        // Count the cells
-        pmFPA *fpa = pmFPAConstruct(config->camera); // The FPA
-        psArray *chips = fpa->chips;        // Array of chips
-        for (int i = 0; i < chips->n; i++) {
-            pmChip *chip = chips->data[i];  // Chip of interest
-            if (!chip) {
-                continue;
-            }
-            psArray *cells = chip->cells;   // Array of cells
-            for (int j = 0; j < cells->n; j++) {
-                pmCell *cell = cells->data[j];
-                if (cell) {
-                numCells++;
-                }
-            }
-        }
-        psFree(fpa);
-    }
+
     psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names
     assert(filenames);                  // It should be here --- it's put here in ppMergeConfig
 
-    psImage *background = psImageAlloc(filenames->n, numCells, PS_TYPE_F64); // Background measurements
+    // Sanity checks
+    assert(!options->scale || scales);
+    assert(!scales || !*scales || ((*scales)->type.type == PS_TYPE_F64 &&
+                                   (*scales)->numCols == data->numCells &&
+                                   (*scales)->numRows == filenames->n));
+    assert(!options->zero || zeros);
+    assert(!zeros || !*zeros || ((*zeros)->type.type == PS_TYPE_F64 &&
+                                 (*zeros)->numCols == data->numCells &&
+                                 (*zeros)->numRows == filenames->n));
+
+    psImage *background = psImageAlloc(filenames->n, data->numCells, PS_TYPE_F64); // Background measurements
     psImageInit(background, NAN);
     psStats *bgStats = psStatsAlloc(options->background); // Statistic to measure the background
-
+    psVector *gains = NULL;             // The gains for each cell
+    psImage *exptime = NULL;            // The exposure time for each integration of each cell
+    if (options->scale) {
+        gains = psVectorAlloc(data->numCells, PS_TYPE_F64);
+    }
+    if (options->exptime) {
+        exptime = psImageAlloc(filenames->n, data->numCells, PS_TYPE_F64);
+    }
+
+    bool status = true;                 // Status of getting the scale and zero --- did everything go right?
     for (int i = 0; i < filenames->n; i++) {
         psString name = filenames->data[i]; // The name of the file
@@ -70 +77 @@
         if (!inFile) {
             psErrorPrint(PS_ERR_IO, false, "Unable to open input file %s --- ignored.\n", name);
+            status = false;
             continue;
         }
-        psMetadata *header = psFitsReadHeader(NULL, inFile); // The FITS (primary) header
-        psMetadata *format = pmConfigCameraFormatFromHeader(config, header); // The camera format
-        psFree(header);
-        if (!format) {
-            psErrorPrint(PS_ERR_IO, false, "Unable to identify camera format for input file %s --- "
-                         "ignored.\n", name);
-            continue;
-        }
-
-        pmFPA *fpa = pmFPAConstruct(config->camera); // The FPA
-        pmFPAview *view = pmFPAAddSourceFromHeader(fpa, phu, format); // View of the PHU
-        int cellNum = 0;                // Number of the cell
+
+        pmFPA *fpa = in->data[i];       // The FPA for this input
+        int cellNum = -1;               // Number of the cell
         psArray *chips = fpa->chips;    // Array of chips
-        for (int i = 0; i < chips->n; i++) {
-            pmChip *chip = chips->data[i]; // The chip of interest
+        for (int j = 0; j < chips->n; j++) {
+            pmChip *chip = chips->data[j]; // The chip of interest
             if (!chip) {
                 continue;
             }
             psArray *cells = chip->cells; // Array of cells
-            for (int j = 0; j < cells->n; j++) {
-                pmCell *cell = cells->data[j]; // The cell of interest
+            for (int k = 0; k < cells->n; k++) {
+                pmCell *cell = cells->data[k]; // The cell of interest
                 if (!cell) {
                     continue;
                 }
-
-                if (!pmCellRead(cell, inFile, config->database)) {
-                    // Nothing here
-                    psFree(cell);
-                    cellNum++;
-                    continue;
-                }
-
-                if (cell->readouts->n > 1) {
-                    psLogMsg(__func__, PS_LOG_WARN, "File %s chip %d cell %d contains more than one readout "
-                             "--- ignoring all but the first.\n", name, i, j);
-                }
-
-                pmReadout *readout = cell->readouts->data[0]; // The readout of interest
-                psImage *image = readout->image; // The pixels of interest
-                if (!image) {
-                    psFree(cell);
-                    cellNum++;
-                    continue;
-                }
-
-                psImageStats(bgStats, image, readout->mask, options->maskVal);
-                background[i][cellNum++] = getStat(bgStats, options->background);
-                psFree(cell);
+                cellNum++;
+
+                if (!pmCellReadHeader(cell, inFile)) {
+                    psLogMsg(__func__, PS_LOG_WARN, "Unable to read header for file %s chip %d cell %d --- "
+                             "ignored.\n", name, j, k);
+                    status = false;
+                }
+
+                // Normalising by the exposure time
+                if (options->exptime) {
+                    bool mdok = true;   // Status of MD lookup
+                    exptime->data.F64[i][cellNum] = psMetadataLookup(&mdok, cell->concepts, "CELL.EXPOSURE");
+                    if (!mdok || isnan(exptime->data.F64[i][cellNum])) {
+                        psLogMsg(__func__, PS_LOG_WARN, "CELL.EXPOSURE for file %s chip %d cell %d is not "
+                                 "set.\n", name, j, k);
+                        exptime->data.F64[i][cellNum] = NAN;
+                        status = false;
+                    }
+                }
+
+                // Scaling by the background
+                if (options->scale || options->zero) {
+                    if (!pmCellRead(cell, inFile, config->database)) {
+                        // Nothing here
+                        psFree(cell);
+                        continue;
+                    }
+
+                    if (cell->readouts->n > 1) {
+                        psLogMsg(__func__, PS_LOG_WARN, "File %s chip %d cell %d contains more than one "
+                                 "readout --- ignoring all but the first.\n", name, j, k);
+                        status = false;
+                    }
+
+                    pmReadout *readout = cell->readouts->data[0]; // The readout of interest
+                    psImage *image = readout->image; // The pixels of interest
+                    if (!image) {
+                        psFree(cell);
+                        continue;
+                    }
+
+                    // Get the gain
+                    if (options->scale) {
+                        bool mdok = true;   // Status of MD lookup
+                        gain->data.F32[cellNum] = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN");
+                        if (!mdok || isnan(gain->data.F32[cellNum])) {
+                            psLogMsg(__func__, PS_LOG_WARN, "CELL.GAIN for file %s chip %d cell %d is not "
+                                     "set.\n", name, j, k);
+                            gain->data.F32[cellNum] = NAN;
+                            status = false;
+                        }
+                    }
+
+                    // Get the background
+                    psImageStats(bgStats, image, readout->mask, options->maskVal);
+                    background[i][cellNum] = getStat(bgStats, options->background);
+                }
+
+                psCellFreeData(cell);
             }
-            psFree(chip);
-        }
-        psFree(fpa);
-        psFree(view);
+            psChipFreeData(chip);
+        }
+        psFPAFreeData(fpa);
         psFitsClose(inFile);
     }
     psFree(bgStats);
 
-    return background;
+    if (options->scale) {
+        // Need to normalize over the focal plane
+        bool converge = true;           // Did we converge?
+        psVector *fluxes = pmFlatNormalize(&converge, gains, background, MAX_ITER, TOLERANCE);
+        if (!converge || !fluxes) {
+            psLogMsg(__func__, PS_LOG_WARN, "Normalisation failed to converge --- continuing anyway.\n");
+            status = false;
+        }
+        psFree(gains);
+
+        *scales = psImageCopy(*scales, background, PS_TYPE_F64);
+        psImage *scalesDeref = *scales; // Dereference the pointer
+
+        for (int i = 0; i < scalesDeref->numRows; i++) {
+            psF64 bg = fluxes->data.F64[i];
+            for (int j = 0; j < scalesDeref->numCols; j++) {
+                scalesDeref->data.F64[i][j] = bg;
+            }
+        }
+    }
+
+    if (options->exptime) {
+        if (!options->scale) {
+            // Copy over the exposure times
+            psImageCopy(*scales, exptime, PS_TYPE_F64);
+        } else {
+            psBinaryOp(*scales, *scales, "*", exptime);
+        }
+    }
+
+    if (options->zero) {
+        if (!*zeros) {
+            // This is much faster than copying!
+            *zeros = psMemIncrRefCounter(background);
+        } else {
+            *zeros = psImageCopy(*zeros, background, PS_TYPE_F64);
+        }
+    }
+
+    psFree(background);
+
+    return status;
 }
+