Changeset 27839 for branches/simmosaic_branches/ppImage/src/ppImageDetrendFringe.c

Timestamp:

May 3, 2010, 8:45:22 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simmosaic_branches

Files:

• . (modified) (1 prop)
• ppImage/src/ppImageDetrendFringe.c (modified) (14 diffs)

branches/simmosaic_branches/ppImage/src/ppImageDetrendFringe.c

@@ -14 +14 @@
     psArray *references = psMemIncrRefCounter(psMetadataLookupPtr(NULL, fringe->analysis, "FRINGE.MEASUREMENTS"));
     if (!references) {
-        psError(PS_ERR_IO, false, "Unable to find fringe references.\n");
-        return false;
+        psError(PS_ERR_IO, false, "Unable to find fringe references.\n");
+        return false;
     }
 
@@ -40 +40 @@
     char *scienceFringes = NULL;
     if (isResidual) {
-        scienceFringes = psStringCopy ("FRINGE.RESIDUALS");
+        scienceFringes = psStringCopy ("FRINGE.RESIDUALS");
     } else {
-        scienceFringes = psStringCopy ("FRINGE.MEASUREMENTS");
+        scienceFringes = psStringCopy ("FRINGE.MEASUREMENTS");
     }
 
@@ -84 +84 @@
     psArray *cells = chip->cells;       // Component cells
     psArray *fringes = psArrayAlloc(cells->n); // Fringes, to return
+    int video_cell_zero = 0;
+
     for (int i = 0; i < cells->n; i++) {
         fringes->data[i] = NULL;
@@ -89 +91 @@
         pmCell *cell = cells->data[i];  // Cell of interest
 
-        // XXX for now, skip the video cells (cell->readouts->n > 1)
-        if (cell->readouts->n > 1) {
-            psWarning ("Skipping Video Cell (%d) for ppImageDetrendFringe.c:getFringes", i);
-            continue;
-        }
+        psTrace("psModules.detrend",7,"Readouts: Cell %d chip: %ld\n",i,cell->readouts->n);
+        // XXX for now, skip the video cells (cell->readouts->n > 1)
+        // CZW: This mess creates a fake set of fringe stats by stealing the previous one.
+        // We let the fitting code know that this is all lies by scaling the weights by a crazy amount.
+
+        if ( (cell->readouts->n > 1) ) {
+          psTrace("psModules.detrend",7,"Should be skipping scichip: %d\n",i);
+          psWarning ("Skipping Video Cell (%d) for ppImageDetrendFringe.c:getFringes", i);
+
+          if (i == 0) {
+            video_cell_zero = 1;
+          }
+          else {
+            pmFringeStats *prevFringe = fringes->data[i-1];
+            pmFringeStats *fringe     = pmFringeStatsAlloc(prevFringe->regions);
+            for (int j = 0; j < fringe->regions->nRequested; j++) {
+              fringe->f->data.F32[j] = prevFringe->f->data.F32[j];
+              fringe->df->data.F32[j] = prevFringe->df->data.F32[j] / 1e6;
+            }
+            fringes->data[i] = fringe;
+          }
+
+          continue;
+        }
 
         fringes->data[i] = psMemIncrRefCounter(psMetadataLookupPtr(NULL, cell->analysis, source));
+    }
+
+    if (video_cell_zero == 1) {
+      pmFringeStats *prevFringe = fringes->data[1];
+      pmFringeStats *fringe     = pmFringeStatsAlloc(prevFringe->regions);
+      for (int j = 0; j < fringe->regions->nRequested; j++) {
+        fringe->f->data.F32[j] = NAN;
+        fringe->df->data.F32[j] = 1.0;
+      }
+      fringes->data[0] = fringe;
     }
 
@@ -112 +143 @@
     psArray *science = NULL;
     if (isResidual) {
-        science = getFringes(scienceChip, "FRINGE.RESIDUALS"); // Fringe residuals on science chip
+        science = getFringes(scienceChip,  "FRINGE.RESIDUALS"); // Fringe residuals on science chip
     } else {
-        science = getFringes(scienceChip, "FRINGE.MEASUREMENTS"); // Fringe measurements on science chip
+        science = getFringes(scienceChip, "FRINGE.MEASUREMENTS"); // Fringe measurements on science chip
     }
 
@@ -124 +155 @@
     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
-    for (int i = 0; i < numRefs; i++) { // Iterate over fringes
-        psArray *refs = psArrayAlloc(references->n); // Array of fringes for each cell
-        for (int j = 0; j < references->n; j++) { // Iterate over cells
-            psArray *ref = references->data[j]; // Array of references for this cell
+    psArray *referencesCat = psArrayAlloc(numRefs);   // Reference fringes
+    for (int i = 0; i < numRefs; i++) {               // Iterate over fringes
+        psArray *refs = psArrayAlloc(references->n);  // Array of fringes for each cell
+        for (int j = 0; j < references->n; j++) {     // Iterate over cells
+            psArray *ref = references->data[j];       // Array of references for this cell
+
             refs->data[j] = psMemIncrRefCounter(ref->data[i]);
         }
@@ -142 +174 @@
 
     if (isResidual) {
-        psMetadataAdd(scienceChip->analysis, PS_LIST_TAIL, "FRINGE.RESIDUAL.SOLUTION", PS_DATA_UNKNOWN, "Fringe solution", solution);
+        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);
+        psMetadataAdd(scienceChip->analysis, PS_LIST_TAIL, "FRINGE.SOLUTION", PS_DATA_UNKNOWN, "Fringe solution", solution);
     }
 
@@ -152 +184 @@
     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);
+        // 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
@@ -192 +224 @@
     }
     assert(fringes->readouts->n == solution->nFringeFrames);
+
+    psImageMaskType maskVal = options->flatMask;
 
     bool mdok;                          // Status of MD lookup
@@ -231 +265 @@
 
 # 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_%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_%02dE", i);
-        psMetadataAddF32(hdu->header, PS_LIST_TAIL, keyword, PS_META_REPLACE, "Fringe Amplitude error", solution->coeffErr->data.F32[i + 1]);
-        psFree (keyword);
+        // write metadata header value
+        // XXX this is measured per cell, but we only have headers per chip
+        psString keyword = NULL;
+        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_%02dE", i);
+        psMetadataAddF32(hdu->header, PS_LIST_TAIL, keyword, PS_META_REPLACE, "Fringe Amplitude error", solution->coeffErr->data.F32[i + 1]);
+        psFree (keyword);
 # endif
 
@@ -257 +291 @@
         }
 
+        // subtract fringe and update mask if fringe value is NAN
+        for (int iy = 0; iy < readout->image->numRows; iy++) {
+            for (int ix = 0; ix < readout->image->numCols; ix++) {
+                readout->image->data.F32[iy][ix] -= sumFringe->data.F32[iy][ix];
+                if (!isfinite(sumFringe->data.F32[iy][ix]) && readout->mask) {
+                    readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= maskVal;
+                }
+            }
+        }
+
         // XXX: Make generic, so subregions may be subtracted as well
-        if (!psBinaryOp(readout->image, readout->image, "-", sumFringe)) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to subtract fringe.\n");
-            return false;
-        }
-
-        // measure residual fringe amplitude. results go to FRINGE.RESIDUALS
-        ppImageDetrendFringeMeasure (readout, fringes, true, options);
+        // if (!psBinaryOp(readout->image, readout->image, "-", sumFringe)) {
+        //     psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to subtract fringe.\n");
+        //     return false;
+        // }
+
+        // measure residual fringe amplitude. results go to FRINGE.RESIDUALS
+        ppImageDetrendFringeMeasure (readout, fringes, true, options);
     }
     psFree(sumFringe);
@@ -276 +320 @@
     psFree(timeString);
 
-    return sumFringe;
+    return true;
 }
 
@@ -294 +338 @@
     pmChip *fringe = pmFPAfileThisChip(config->files, view, "PPIMAGE.FRINGE");
     if (!fringe) {
-        psError(PS_ERR_UNKNOWN, false, "missing fringe reference data.\n");
-        psFree (view);
-        return false;
+        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;
+        psError(PS_ERR_UNKNOWN, false, "failed to solve the fringe system.\n");
+        psFree (view);
+        return false;
     }
 
@@ -309 +353 @@
     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;
-        }
-
-        // XXX for now, skip the video cells (cell->readouts->n > 1)
-        if (cell->readouts->n > 1) {
-          psWarning ("Skipping Video Cell for ppImageDetrendFringeApply");
-          continue;
-        }
-
-        if (!ppImageDetrendFringeGenerate(cell, fringeCell, options)) {
-            psError(PS_ERR_UNKNOWN, false, "failed to apply fringe image.\n");
-            psFree (view);
-            return false;
-        }
+        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;
+        }
+
+        // XXX for now, skip the video cells (cell->readouts->n > 1)
+        if (cell->readouts->n > 1) {
+          psWarning ("Skipping Video Cell for ppImageDetrendFringeApply");
+          continue;
+        }
+
+        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;
+        psError(PS_ERR_UNKNOWN, false, "failed to solve the residual fringe system.\n");
+        psFree (view);
+        return false;
     }