Changeset 7061 for trunk/ppMerge

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)
• .cvsignore (added)
• Makefile.am (modified) (1 diff)
• ppMerge.c (modified) (2 diffs)
• ppMerge.h (modified) (1 diff)
• ppMergeBackground.c (modified) (3 diffs)
• ppMergeBackground.h (modified) (1 diff)
• ppMergeCheckInputs.c (modified) (2 diffs)
• ppMergeCheckInputs.h (modified) (1 diff)
• ppMergeCombine.c (added)
• ppMergeCombine.h (added)
• ppMergeData.c (added)
• ppMergeData.h (added)
• ppMergeOptions.c (modified) (2 diffs)

trunk/ppMerge/src/Makefile.am

@@ -4 +4 @@
 ppImage_LDADD = $(PSLIB_LIBS) $(PSMODULE_LIBS)
 ppImage_SOURCES = \
-        ppConfig.c \
-        ppFile.c \
-        ppMem.c \
-        ppImage.c \
-        ppImageConfig.c \
-        ppImageData.c \
-        ppImageDetrendBias.c \
-        ppImageDetrendCell.c \
-        ppImageDetrendMask.c \
-        ppImageDetrendNonLinear.c \
-        ppImageLoop.c \
-        ppImageOptions.c \
-        ppImageParseCamera.c \
-        ppImageParseDetrend.c \
-        ppImageWeights.c
-#       ppImageLoadPixels.c
-#       ppDetrendFlat.c
-#       ppImageOutput.c
-#       ppImageDetrendPedestal.c
-#       ppImageOutput.c
-#       ppImagePhot.c
+        ppMerge.c               \
+        ppMergeBackground.c     \
+        ppMergeCheckInputs.c    \
+        ppMergeCombine.c        \
+        ppMergeConfig.c         \
+        ppMergeData.c           \
+        ppMergeOptions.c
 
 
 noinst_HEADERS = \
-        ppConfig.h \
-        ppFile.h \
-        ppMem.h \
-        ppImage.h \
-        ppImageData.h \
-        ppImageDetrend.h \
-        ppImageOptions.h
+        ppMerge.h               \
+        ppMergeBackground.h     \
+        ppMergeCheckInputs.h    \
+        ppMergeCombine.h        \
+        ppMergeConfig.h         \
+        ppMergeData.h           \
+        ppMergeOptions.h
 
 clean-local:

trunk/ppMerge/src/ppMerge.c

@@ -1 +1 @@
 #include <stdio.h>
-#include "pslib.h"
-#include "pmConcepts.h"
+#include <pslib.h>
+#include <psmodules.h>
+
 #include "ppMerge.h"
 #include "ppMem.h"
@@ -13 +14 @@
 
     // Parse the configuration and arguments
-    ppImageConfig(config, argc, argv);
-
-    // Open the input image, output image, output mask
-    // Get camera configuration from header if not already defined
+    // Open the input image(s)
+    // Determine camera, format from header if not already defined
     // Construct camera in preparation for reading
-    ppImageParseCamera(data, config);
+    pmConfig *config = ppMergeConfig(&argc, argv);
 
     // Set various tasks (define optional operations)
-    ppImageOptionsParse(data, options, config);
+    ppMergeOptions *options = ppMergeOptionsParse(config);
 
-    // open detrend images, load headers, optionally load pixels
-    ppImageParseDetrend(data, options, config);
+    // Check the inputs
+    ppMergeData *data = ppMergeCheckInputs(options, config);
+    if (!data) {
+        psError("Not enough valid input files.\n");
+        exit(EXIT_FAILURE);
+    }
 
-    // Image Arithmetic Loop
-    ppImageLoop(data, options, config);
+    // Measure the background in each image
+    psImage *scale = NULL;              // The scalings
+    psImage *zero = NULL;               // The zeros
+    ppMergeScaleZero(&scale, &zero, data, options, config);
+
+    // Do the combination and write
+    ppMergeCombine(scale, zero, data, options, config);
 
     // Cleaning up
     psFree(data);
+    psFree(scale);
+    psFree(zero);
     psFree(options);
     psFree(config);
+
     psTimerStop();
     psTraceReset();
+    psErrorClear();
     pmConceptsDone();
+    pmConfigDone();
+
     ppMemCheck();
 

trunk/ppMerge/src/ppMerge.h

@@ -1 @@
-# include <stdio.h>
-# include <strings.h>
-# include <glob.h>
+#ifndef PP_MERGE_H
+#define PP_MERGE_H
 
-# include "pslib.h"
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
 
-# include "psAdditionals.h"
+#define RECIPE "PPMERGE"                // Name of the recipe to use
 
-# include "pmAstrometry.h"
-# include "pmReadout.h"
-# include "pmConfig.h"
-# include "pmFPAConstruct.h"
-# include "pmFPARead.h"
-# include "pmFPAConceptsGet.h"
-# include "pmFPAWrite.h"
-
-# include "pmReadoutCombine.h"
-
-# define RECIPE "MERGE" // Name of the recipe to use
-
-// XXX : same as in ppImage
-typedef enum {
-    PP_LOAD_NONE,
-    PP_LOAD_FPA,
-    PP_LOAD_CHIP,
-    PP_LOAD_CELL,
-} ppImageLoadDepth;
-
-// XXX : same as in ppImage
-typedef struct {
-    psMetadata *site;
-    psMetadata *camera;
-    psMetadata *recipe;
-    psMetadata *arguments;
-    psDB       *database;
-} ppConfig;
-
-typedef struct {
-    char *filename;
-    pmFPA *fpa;
-    psFits *fits;
-    psMetadata *header;
-    double zero;
-    double scale;
-} ppFPA; 
-
-typedef struct {
-    psArray *input;
-    ppFPA *output;
-    ppFPA *process;
-    ppFPA *mask;
-} ppData;
-
-typedef struct {
-    ppImageLoadDepth imageLoadDepth;
-    bool doMask;                        // apply a pixel mask before stacking
-    pmCombineParams *combineParams;
-    bool applyZeroScale;
-    double gain;
-    double readnoise;
-} ppOptions;
-
-bool ppMergeConfig (ppConfig *config, int argc, char **argv);
-bool ppMergeLoop (ppData *data, ppOptions *options, ppConfig *config);
-bool ppMergeOptions (ppData *data, ppOptions *options, ppConfig *config);
-bool ppMergeOutput (ppData *data, ppOptions *options, ppConfig *config);
-bool ppMergeParseCamera (ppData *data, ppConfig *config);
-bool ppMergeParseDetrend (ppData *data, ppOptions *options, ppConfig *config);
-
-bool ppMergeCell (pmCell *output, pmCell *mask, psArray *cellList, ppOptions *options, ppConfig *config);
-
-// XXX : these functions are identical to the ppImage equivalents
-pmReadout* ppDetrendSelectFirst (pmCell *cell, char *name, bool doThis);
-bool ppFPAOpen (ppFPA *fpa, psMetadata *camera, char *name, bool doThis);
-bool ppMergeLoadPixels (ppFPA *input, ppFPA *process, psDB *db, int nChip, int nCell);

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;
 }
+

trunk/ppMerge/src/ppMergeBackground.h

@@ -2 +2 @@
 #define PP_MERGE_BACKGROUND_H
 
-// 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
     );
 

trunk/ppMerge/src/ppMergeCheckInputs.c

@@ -5 +5 @@
 
 #include "ppMergeCheckInputs.h"
+#include "ppMergeData.h"
 
 // Check input files to make sure everything's consistent
-bool ppCheckInputs(const pmConfig *config // Configuration
+ppMergeData *ppMergeCheckInputs(ppMergeOptions *options, // Options
+                                const pmConfig *config // Configuration
     )
 {
+    ppMergeData *data = ppMergeDataAlloc(); // The data, to return
+
+    // Output file
+    psString outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // The output file name
+    assert(outName);                    // It should be there!
+    data->outFile = psFitsOpen(outName, "w"); // Output FITS file
+    if (!data->outFile) {
+        // There's no point in continuing if we can't open the output
+        psErrorStackPrint(stderr, "Can't open output image: %s\n", outName);
+        exit(EXIT_FAILURE);
+    }
+
     psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names
     assert(filenames);
@@ -27 +41 @@
         }
         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);
+        psFitsClose(inFile);
+
+        // The formats must be identical.  The chief reason for this is so that we know what output format to
+        // use.  I guess one could specify a different output format on the command line, but how do we
+        // generate a PHU for that?  Perhaps we could revisit this restriction in the future (construct an
+        // FPAview from the specified camera format configuration, and use pmFPAAddSourceFromView), but for
+        // now it's less hassle just to limit the output format to be the input format.
+        if (!data->format) {
+            data->format = pmConfigCameraFormatFromHeader(config, header);
+            psFree(header);
+            if (!options->format) {
+                psLogMsg(__func__, PS_LOG_WARN, "Unable to identify camera format for input file %s --- "
+                             "ignored.\n", name);
+                // Kick it out
+                psFree(header);
+                data->in->data[i] = NULL;
+                continue;
+            }
+        } else if (!pmConfigValidateCameraFormat(options->format, header)) {
+            psLogMsg(__func__, PS_LOG_WARN, "Input file %s doesn't match camera format --- ignored.\n", name);
             // Kick it out
-            psFree(filenames->data[i]);
-            filenames->data[i] = NULL;
+            psFree(header);
+            data->in->data[i] = NULL;
             continue;
         }
-        psFitsClose(inFile);
+
+        data->in->data[i] = pmFPAConstruct(config->camera);
+        pmFPAview *view = pmFPAAddSourceFromHeader(data->in->data[i], header, data->format);
+        psFree(view);
+
+        // Use the first valid input as the basis for the output --- including the header
+        if (!data->out) {
+            data->out = pmFPAConstruct(config->camera);
+            pmFPAview *view = pmFPAAddSourceFromHeader(data->out, header, data->format);
+            psFree(view);
+        }
+
         numGood++;
     }
 
+    // Count the cells
+    int numCells = 0;           // Number of cells in the output FPA
+    psArray *chips = data->out->chips; // Array of chips in output
+    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++;
+                }
+        }
+    }
+    data->numCells = numCells;
+
     return (numGood > 1);
 }
+
+

trunk/ppMerge/src/ppMergeCheckInputs.h

@@ -3 +3 @@
 
 // Check input files to make sure everything's consistent
-bool ppCheckInputs(const pmConfig *config // Configuration
+bool ppMergeCheckInputs(ppMergeOptions *options, // Options
+                        const pmConfig *config // Configuration
     );
 

trunk/ppMerge/src/ppMergeOptions.c

@@ -7 +7 @@
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
+#if 0
 // Free function
 static void mergeOptionsFree(ppMergeOptions *options // Options to free
     )
 {
-    psFree(options->combine);
 }
+#endif
 
 // Allocator
@@ -18 +19 @@
 {
     ppMergeOptions *options = psAlloc(sizeof(ppMergeOptions)); // The options, to return
-    psMemSetDeallocator(options, (psFreeFunc)mergeOptionsFree);
+    // psMemSetDeallocator(options, (psFreeFunc)mergeOptionsFree);
 
     options->rows = 0;