Changeset 25930

Timestamp:

Oct 22, 2009, 4:54:27 PM (17 years ago)

Author:

watersc1

Message:

Updated fringe fitting code.

Corrected Bool->S32 bug for BTOOLAPP header concept.

Location:

trunk

Files:

• ppImage/src/Makefile.am (modified) (1 diff)
• ppImage/src/ppImage.h (modified) (1 diff)
• ppImage/src/ppImageDetrendFringe.c (modified) (4 diffs)
• ppImage/src/ppImageDetrendPattern.c (added)
• ppImage/src/ppImageDetrendReadout.c (modified) (1 diff)
• ppImage/src/ppImageLoop.c (modified) (1 diff)
• ppImage/src/ppImageOptions.c (modified) (1 diff)
• psModules/src/concepts/pmConcepts.c (modified) (2 diffs)
• psModules/src/concepts/pmConceptsRead.c (modified) (1 diff)
• psModules/src/concepts/pmConceptsStandard.c (modified) (2 diffs)
• psModules/src/concepts/pmConceptsStandard.h (modified) (1 diff)
• psModules/src/detrend/pmFringeStats.c (modified) (7 diffs)

trunk/ppImage/src/Makefile.am

@@ -41 +41 @@
         ppImageDetrendFringe.c \
         ppImageDetrendFree.c \
+        ppImageDetrendPattern.c \
         ppImageRebinReadout.c \
         ppImageMosaic.c \

trunk/ppImage/src/ppImage.h

@@ -159 +159 @@
 bool ppImageBurntoolMask(pmConfig *config, ppImageOptions *options, pmFPAview *view, pmReadout *mask);
 
+bool ppImageDetrendPatternApply(pmConfig *config, pmChip *chip, const pmFPAview *inputView, const ppImageOptions *options);
+
 // Record which detrend file was used for the detrending
 bool ppImageDetrendRecord(

trunk/ppImage/src/ppImageDetrendFringe.c

@@ -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
 
+        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)
-        if (cell->readouts->n > 1) {
-            psWarning ("Skipping Video Cell (%d) for ppImageDetrendFringe.c:getFringes", i);
-            continue;
-        }
-
+        // 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;
+    }
+    
     return fringes;
 }
@@ -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
@@ -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]);
         }
-        referencesCat->data[i] = pmFringeStatsConcatenate(refs, NULL, NULL);
+        referencesCat->data[i] = pmFringeStatsConcatenate(refs, NULL, NULL);
         psFree(refs);
     }

trunk/ppImage/src/ppImageDetrendReadout.c

@@ -128 +128 @@
     }
 
-    // Pattern noise correction
-    if (options->doPattern) {
-        if (!pmPatternRow(input, options->patternOrder, options->patternIter, options->patternRej,
-                          options->patternThresh, options->patternMean, options->patternStdev,
-                          options->maskValue, options->darkMask)) {
-            psFree(detview);
-            return false;
-        }
-    }
+/*     // Pattern noise correction */
+/*     if (options->doPattern) { */
+/*         if (!pmPatternRow(input, options->patternOrder, options->patternIter, options->patternRej, */
+/*                           options->patternThresh, options->patternMean, options->patternStdev, */
+/*                           options->maskValue, options->darkMask)) { */
+/*             psFree(detview); */
+/*             return false; */
+/*         } */
+/*     } */
 
     // Normalization by a single (known) constant

trunk/ppImage/src/ppImageLoop.c

@@ -163 +163 @@
         }
 
+        // Apply the pattern noise correction
+        if (options->doPattern) {
+          if (!ppImageDetrendPatternApply(config,chip,view,options)) {
+            ESCAPE("Unable to apply pattern corrections");
+          }
+        }
+        
         // measure various pixel-based statistics for this image
         if (!ppImagePixelStats(config, stats, options, view)) {

trunk/ppImage/src/ppImageOptions.c

@@ -249 +249 @@
 
     options->burntoolTrails = psMetadataLookupS32(&status, recipe, "BURNTOOL.TRAILS");
-    fprintf(stderr,"TRAILS: %d %d %d\n",options->burntoolTrails,psMetadataLookupS32(&status,recipe,"BURNTOOL.TRAILS"),status);
+    psTrace("psModules.detrend", 7, "burntoolTrails: %d BURNTOOL.TRAILS: %d Status: %d\n",
+            options->burntoolTrails,psMetadataLookupS32(&status,recipe,"BURNTOOL.TRAILS"),status);
     if (!status) {
       psWarning("BURNTOOL.TRAILS not found in recipe: setting to default value.\n");

trunk/psModules/src/concepts/pmConcepts.c

@@ -249 +249 @@
 CONCEPT_REGISTER_FUNCTION(S32, Enum, -1); // For enums: set default to -1
 CONCEPT_REGISTER_FUNCTION(S32, S32, 0); // For values: set default to 0
+//CONCEPT_REGISTER_FUNCTION(Bool, Bool, NULL); // For values: set default to 0
 
 static void conceptRegisterTime(const char *name, /* Name of concept */ \
@@ -333 +334 @@
         conceptRegisterF32("FPA.TELTEMP.EXTRA", "Telescope Temperatures: extra", p_pmConceptParse_TELTEMPS, NULL, NULL, false, PM_FPA_LEVEL_FPA);
         conceptRegisterF32("FPA.PON.TIME", "Power On Time", NULL, NULL, NULL, false, PM_FPA_LEVEL_FPA);
-        conceptRegisterS32("FPA.BURNTOOL.APPLIED", "Status of burntool processing", p_pmConceptParse_BTOOLAPP,NULL,NULL,false,PM_FPA_LEVEL_FPA);
+        conceptRegisterS32("FPA.BURNTOOL.APPLIED", "Status of burntool processing", p_pmConceptParse_BTOOLAPP,p_pmConceptFormat_BTOOLAPP,NULL,false,PM_FPA_LEVEL_FPA);
+        //      conceptRegisterBool("FPA.BURNTOOL.APPLIED", "Status of burntool processing", NULL, NULL, NULL, false, PM_FPA_LEVEL_FPA);
         conceptRegisterF32("FPA.EXPOSURE", "Exposure time (sec)", NULL, NULL, NULL, false, PM_FPA_LEVEL_FPA);
     }

trunk/psModules/src/concepts/pmConceptsRead.c

@@ -47 +47 @@
     case PS_DATA_F64:
         return psMetadataItemAllocF64(pattern->name, pattern->comment, psMetadataItemParseF64(concept));
+    case PS_DATA_BOOL:
+      return psMetadataItemAllocBool(pattern->name, pattern->comment, psMetadataItemParseBool(concept));
     default:
         psWarning("Concept %s (%s) is not of a standard type (%x)\n",

trunk/psModules/src/concepts/pmConceptsStandard.c

@@ -753 +753 @@
     psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Type for %s (%x) is not BOOL\n",
             concept->name, concept->type);
-    return NULL;
+    if (concept->type != PS_DATA_S32) {
+      psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Wasn't the type I'd guessed either.\n");
+      return NULL;
+    }
+    psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Looks like an S32 value? (%d)\n",
+            concept->data.S32);
+
+    if (concept->data.S32 == 0) {
+      bt_status = 1;
+    }
+    else if (concept->data.S32 == 1) {
+      bt_status = -2;
+    }
   }
 
@@ -765 +777 @@
   return psMetadataItemAllocS32(concept->name, concept->comment, bt_status);
 }  
+psMetadataItem *p_pmConceptFormat_BTOOLAPP(const psMetadataItem *concept,
+                                           pmConceptSource source,
+                                           const psMetadata *cameraFormat,
+                                           const pmFPA *fpa,
+                                           const pmChip *chip,
+                                           const pmCell *cell)
+{
+  assert(concept);
+
+  if (concept->type != PS_DATA_S32) {
+    return NULL;
+  }
+
+  if (concept->data.S32 == 0) {
+    return NULL;
+  }
+  else if (concept->data.S32 == -2) {
+    return psMetadataItemAllocBool(concept->name,concept->comment,true);
+  }
+  else if (concept->data.S32 == 1) {
+    return psMetadataItemAllocBool(concept->name,concept->comment,false);
+  }
+  else {
+    return NULL;
+  }
+
+}  
+
 
 // Get the current value of a concept

trunk/psModules/src/concepts/pmConceptsStandard.h

@@ -146 +146 @@
     const pmCell *cell ///< Cell for concept, or NULL
     );
+psMetadataItem *p_pmConceptFormat_BTOOLAPP(
+    const psMetadataItem *concept, ///< Concept to format
+    pmConceptSource source, ///< Source for concept
+    const psMetadata *cameraFormat, ///< Camera format definition
+    const pmFPA *fpa, ///< FPA for concept, or NULL
+    const pmChip *chip, ///< Chip for concept, or NULL
+    const pmCell *cell ///< Cell for concept, or NULL
+    );
 
 

trunk/psModules/src/detrend/pmFringeStats.c

@@ -11 +11 @@
 #include "pmFPA.h"
 #include "pmFringeStats.h"
+
+#include "psPolynomialMD.h"
+#include "psMinimizePolyFit.h"
+#include "psVector.h"
+
 
 // Future optimisations for speed:
@@ -331 +336 @@
         dfPt[i] = 1.0 / medianSd->sampleStdev;
 
-        psTrace("psModules.detrend", 7, "[%d:%d,%d:%d]: %f %f\n", (int)region.x0, (int)region.x1,
-                (int)region.y0, (int)region.y1, fPt[i], dfPt[i]);
+        psTrace("psModules.detrend", 7, "[%d:%d,%d:%d]: %f %f : %s\n", (int)region.x0, (int)region.x1,
+                (int)region.y0, (int)region.y1, fPt[i], dfPt[i], readout->parent->hdu->extname);
     }
     psFree(sky);
@@ -840 +845 @@
             }
         }
-        B->data.F64[i] = vector;
+        B->data.F64[i] = vector;
     }
 
@@ -966 +971 @@
             }
         }
+        
     }
 
@@ -988 +994 @@
                 psTrace("psModules.detrend", 9, "Masking region %d because not finite in fringe %d.\n", j, i);
             }
-        }
-    }
-
+            else if (fabs(fringe->f->data.F32[j]) > 0.1) {
+              mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 1;
+              psTrace("psModules.detrend", 9, "Masking region %d because too large fringe %d.\n", j, i);
+            }
+            // Mask bad points in the science data as well.
+            if ((i == 0) && (!isfinite(science->f->data.F32[j]))) {
+                mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 1;
+                psTrace("psModules.detrend", 9, "Masking region %d because not finite in science fringe %d.\n", j, i);
+            }         
+            psTrace("psModules.detrend", 7, "F %f %f %f %d\n",
+                    fringe->f->data.F32[j], science->f->data.F32[j],
+                    1 / science->df->data.F32[j],(int) mask->data.PS_TYPE_VECTOR_MASK_DATA[j]);
+        }
+    }
+    
+/*     // Begin switch from old outlier removal and fitting code. */
+
+    psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1);
+
+    pmFringeStats *fringe = fringes->data[0];
+    psVector *errors = psVectorAlloc(science->df->n,PS_TYPE_F32);
+    for (int j = 0; j < errors->n; j++) {
+      errors->data.F32[j] = 1 / science->df->data.F32[j];
+    }
+    psVectorFitPolynomial1D(poly,mask,0xff,science->f,errors,fringe->f);
+
+    for (int i = 0; i <= poly->nX; i++) {
+      scale->coeff->data.F32[i] = poly->coeff[i];
+      psTrace("psModules.detrend",7,"COEFFS: %d %g %g %g\n",i,scale->coeff->data.F32[i],poly->coeff[i],poly->coeffErr[i]);
+    }
+
+    psFree(poly);
+    //    psFree(fringe);
+    psFree(errors);
+
+    psFree(median);
+    psFree(diff);
+    return scale;
+    // End switch from old code.
+    
+
+    
 # if (0)
     // Write fringe data to file for a test
@@ -1019 +1064 @@
         }
         scale->coeff->data.F32[0] -= median->sampleMedian;
-        scale->coeff->data.F32[i] = 0.0;
+        if (i != 0) {
+          scale->coeff->data.F32[i] = 0.0;
+        }
     }
     psFree(median);
@@ -1051 +1098 @@
 
     return scale;
+    //# endif
 }