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

Timestamp:

Jun 1, 2006, 10:16:15 AM (20 years ago)

Author:

Paul Price

Message:

Working version for flat-fielding. Need to add extenal scale/zero input, mask saturated pixels, test on dark, bias, fringe.

File:

• trunk/ppMerge/src/ppMergeBackground.c (modified) (9 diffs)

trunk/ppMerge/src/ppMergeBackground.c

@@ -49 +49 @@
     // Sanity checks
     assert(!options->scale || scales);
-    assert(!scales || !*scales || ((*scales)->type.type == PS_TYPE_F64 &&
+    assert(!scales || !*scales || ((*scales)->type.type == PS_TYPE_F32 &&
                                    (*scales)->numCols == data->numCells &&
                                    (*scales)->numRows == filenames->n));
     assert(!options->zero || zeros);
-    assert(!zeros || !*zeros || ((*zeros)->type.type == PS_TYPE_F64 &&
+    assert(!zeros || !*zeros || ((*zeros)->type.type == PS_TYPE_F32 &&
                                  (*zeros)->numCols == data->numCells &&
                                  (*zeros)->numRows == filenames->n));
 
-    psImage *background = psImageAlloc(filenames->n, data->numCells, PS_TYPE_F64); // Background measurements
+    psImage *background = psImageAlloc(data->numCells, filenames->n, PS_TYPE_F32); // Background measurements
     psImageInit(background, NAN);
     psStats *bgStats = psStatsAlloc(options->background); // Statistic to measure the background
@@ -63 +63 @@
     psImage *exptime = NULL;            // The exposure time for each integration of each cell
     if (options->scale) {
-        gains = psVectorAlloc(data->numCells, PS_TYPE_F64);
+        gains = psVectorAlloc(data->numCells, PS_TYPE_F32);
+        gains->n = data->numCells;
     }
     if (options->exptime) {
-        exptime = psImageAlloc(filenames->n, data->numCells, PS_TYPE_F64);
+        exptime = psImageAlloc(filenames->n, data->numCells, PS_TYPE_F32);
     }
 
@@ -75 +76 @@
             continue;
         }
+        psTrace(__func__, 9, "Opening %s to get background...\n", name);
         psFits *inFile = psFitsOpen(filenames->data[i], "r"); // The FITS file to read
         if (!inFile) {
@@ -107 +109 @@
                 if (options->exptime) {
                     bool mdok = true;   // Status of MD lookup
-                    exptime->data.F64[i][cellNum] = psMetadataLookupF32(&mdok, cell->concepts,
+                    exptime->data.F32[i][cellNum] = psMetadataLookupF32(&mdok, cell->concepts,
                                                                         "CELL.EXPOSURE");
-                    if (!mdok || isnan(exptime->data.F64[i][cellNum])) {
+                    if (!mdok || isnan(exptime->data.F32[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;
+                        exptime->data.F32[i][cellNum] = NAN;
                         status = false;
                     }
@@ -151 +153 @@
 
                     // Get the background
-                    psImageStats(bgStats, image, readout->mask, options->combine->maskVal);
-                    background->data.F64[i][cellNum] = getStat(bgStats, options->background);
+                    int sampleSize = (image->numCols * image->numRows) / options->sample; // Size of sample
+                    psVector *sample = psVectorAlloc(sampleSize, PS_TYPE_F32); // Sample of the image
+                    sample->n = sampleSize;
+                    psVector *sampleMask = NULL; // Mask for sample
+                    if (readout->mask) {
+                        sampleMask = psVectorAlloc(sampleSize, PS_TYPE_U8);
+                        sampleMask->n = sampleSize;
+                    }
+                    psImage *mask = readout->mask; // The mask image
+                    for (long i = 0; i < sampleSize; i++) {
+                        int j = i * options->sample; // Index into image
+                        int x = j % image->numCols; // x index
+                        int y = j / image->numCols; // y index
+                        sample->data.F32[i] = image->data.F32[y][x];
+                        if (readout->mask) {
+                            sampleMask->data.U8[i] = mask->data.U8[y][x];
+                        }
+                    }
+                    psVectorStats(bgStats, sample, sampleMask, NULL, options->combine->maskVal);
+                    psFree(sample);
+                    psFree(sampleMask);
+                    background->data.F32[i][cellNum] = getStat(bgStats, options->background);
+                    psTrace(__func__, 3, "Background for %s, cell %d is %f\n", name, cellNum,
+                            background->data.F32[i][cellNum]);
                 }
 
@@ -167 +191 @@
         // Need to normalize over the focal plane
         bool converge = true;           // Did we converge?
+        if (psTraceGetLevel(__func__) > 9) {
+            for (int i = 0; i < gains->n; i++) {
+                psTrace(__func__, 10, "Gain for cell %d is %f\n", i, gains->data.F32[i]);
+            }
+        }
         psVector *fluxes = pmFlatNormalize(&converge, gains, background, MAX_ITER, TOLERANCE);
         if (!converge || !fluxes) {
@@ -174 +203 @@
         psFree(gains);
 
-        *scales = psImageCopy(*scales, background, PS_TYPE_F64);
+        *scales = psImageCopy(*scales, background, PS_TYPE_F32);
         psImage *scalesDeref = *scales; // Dereference the pointer
 
         for (int i = 0; i < scalesDeref->numRows; i++) {
-            psF64 bg = fluxes->data.F64[i];
+            psF32 bg = fluxes->data.F32[i];
             for (int j = 0; j < scalesDeref->numCols; j++) {
-                scalesDeref->data.F64[i][j] = bg;
+                scalesDeref->data.F32[i][j] = bg;
             }
         }
@@ -188 +217 @@
         if (!options->scale) {
             // Copy over the exposure times
-            psImageCopy(*scales, exptime, PS_TYPE_F64);
+            psImageCopy(*scales, exptime, PS_TYPE_F32);
         } else {
             psBinaryOp(*scales, *scales, "*", exptime);
@@ -199 +228 @@
             *zeros = psMemIncrRefCounter(background);
         } else {
-            *zeros = psImageCopy(*zeros, background, PS_TYPE_F64);
+            *zeros = psImageCopy(*zeros, background, PS_TYPE_F32);
+        }
+    }
+
+    // Diagnostic stuff
+    if (scales && *scales && psTraceGetLevel(__func__) > 9) {
+        psImage *scalesDeref = *scales; // Dereference the pointer
+        for (int i = 0; i < scalesDeref->numRows; i++) {
+            for (int j = 0; j < scalesDeref->numCols; j++) {
+                psTrace(__func__, 9, "Scale for exposure %d, cell %d is: %f\n", i, j,
+                        scalesDeref->data.F32[i][j]);
+            }
+        }
+    }
+    if (zeros && *zeros && psTraceGetLevel(__func__) > 9) {
+        psImage *zerosDeref = *zeros; // Dereference the pointer
+        for (int i = 0; i < zerosDeref->numRows; i++) {
+            for (int j = 0; j < zerosDeref->numCols; j++) {
+                psTrace(__func__, 9, "Zero for exposure %d, cell %d is: %f\n", i, j,
+                        zerosDeref->data.F32[i][j]);
+            }
         }
     }