Changeset 29486 for branches/czw_branch/20100817/psModules/src/detrend/pmNonLinear.c

Timestamp:

Oct 19, 2010, 5:41:33 PM (16 years ago)

Author:

watersc1

Message:

Still broken with weird memory errors, but committing anyway due to JHU trip

File:

• branches/czw_branch/20100817/psModules/src/detrend/pmNonLinear.c (modified) (2 diffs)

branches/czw_branch/20100817/psModules/src/detrend/pmNonLinear.c

@@ -45 +45 @@
         } }
 
+
+# define PS_BIN_INTERPOLATE(RESULT, VECTOR, BOUNDS, BIN, VALUE) { \
+    float dX, dY, Xo, Yo, Xt;                                     \
+    if (BIN == BOUNDS->n - 1) {                                   \
+      dX = 0.5*(BOUNDS->data.F32[BIN+1] - BOUNDS->data.F32[BIN-1]);     \
+      dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN-1];             \
+      Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]);       \
+      Yo = VECTOR->data.F32[BIN];                                       \
+    } else {                                                            \
+      dX = 0.5*(BOUNDS->data.F32[BIN+2] - BOUNDS->data.F32[BIN]);       \
+      dY = VECTOR->data.F32[BIN+1] - VECTOR->data.F32[BIN];             \
+      Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]);       \
+      Yo = VECTOR->data.F32[BIN];                                       \
+    }                                                                   \
+    if (dY != 0) {                                                      \
+      Xt = (VALUE - Yo)*dX/dY + Xo;                                     \
+    } else {                                                            \
+      Xt = Xo;                                                          \
+    }                                                                   \
+    Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \
+    psTrace("pmNonLinear", 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n",  \
+            Xo, Yo, dX, dY, Xt, VALUE);                                 \
+    RESULT = Xt; }
+
 pmReadout *pmNonLinearityLookup(pmReadout *inputReadout, const psVector *inFlux, const psVector *outFlux)
 {
@@ -92 +116 @@
     return inputReadout;
 }
+
+bool pmNonLinearityApply(pmReadout *inputReadout, psArray *Ltab)
+{
+  PS_ASSERT_PTR_NON_NULL(inputReadout, false);
+  PS_ASSERT_PTR_NON_NULL(inputReadout->image, false);
+  PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, false);
+  PS_ASSERT_PTR_NON_NULL(Ltab, false);
+
+  //  psS32 tableSize = Ltab->n;
+
+  psImage *image = inputReadout->image;
+  psS32 numSamples = 39;
+  
+  psS32 numBorder  = 10;
+  // Load default data.
+  psVector *default_row_correction_fluxes = psVectorAlloc(numSamples,PS_TYPE_F32);
+  psVector *default_col_correction_fluxes = psVectorAlloc(numSamples,PS_TYPE_F32);
+  psVector *row_correction_fluxes;
+  psVector *col_correction_fluxes;
+
+  psVector *default_row_correction_factors = psVectorAlloc(numSamples,PS_TYPE_F32);
+  psVector *default_col_correction_factors = psVectorAlloc(numSamples,PS_TYPE_F32);
+  psVector *row_correction_factors;
+  psVector *col_correction_factors;
+
+  int n = 0;
+  int m = 0;
+  for (int k = 0; k < Ltab->n; k++) { // Begin load default tables
+    psMetadata *row = Ltab->data[k];    
+    if (psMetadataLookupS32(NULL,row,"POSITION") != -1) {
+      continue;
+    }
+    if (psMetadataLookupS32(NULL,row,"DIRECTION") == 0) {
+      psVectorSet(default_row_correction_fluxes,n,psMetadataLookupF32(NULL,row,"FLUX"));
+      psVectorSet(default_row_correction_factors,n,psMetadataLookupF32(NULL,row,"FACTOR"));
+      n++;
+    }
+    else {
+      psVectorSet(default_col_correction_fluxes,m,psMetadataLookupF32(NULL,row,"FLUX"));
+      psVectorSet(default_col_correction_factors,m,psMetadataLookupF32(NULL,row,"FACTOR"));
+      m++;
+    }
+  } // End load default tables
+  
+  for (int i = 0; i < image->numRows; i++) { // Loop over rows : note: problem with discontinuity here
+    n = 0;
+    if ((i < numBorder)||(image->numRows - i < numBorder)) {     // Load row correction data
+      row_correction_fluxes = psVectorAlloc(numSamples,PS_TYPE_F32);
+      row_correction_factors = psVectorAlloc(numSamples,PS_TYPE_F32);
+
+      for (int k = 0; k < Ltab->n; k++) {
+        psMetadata *row = Ltab->data[k];        
+        if ((psMetadataLookupS32(NULL,row,"DIRECTION") != 0)||
+            (psMetadataLookupS32(NULL,row,"POSITION") != i + 1)) { // image data is 0 indexed, but correction data is 1 indexed.
+          continue;
+        }
+        psVectorSet(row_correction_fluxes,n,psMetadataLookupF32(NULL,row,"FLUX"));
+        psVectorSet(row_correction_factors,n,psMetadataLookupF32(NULL,row,"FACTOR"));
+        //      psTrace("psModules.camera",6,"rows: %d\n",n);
+        n++;
+      }
+    } // End load specific row
+    else {
+      row_correction_fluxes = default_row_correction_fluxes;
+      row_correction_factors = default_row_correction_factors;
+    } // End load default row
+      
+    for (int j = 0; j < image->numCols; j++) { // Loop over columns
+      m = 0;
+
+      if ((j < numBorder)||(image->numCols - j < numBorder)) {       // Load col correction data
+        col_correction_fluxes = psVectorAlloc(numSamples,PS_TYPE_F32);
+        col_correction_factors = psVectorAlloc(numSamples,PS_TYPE_F32);
+        for (int k = 0; k < Ltab->n; k++) {
+          psMetadata *row = Ltab->data[k];      
+          if ((psMetadataLookupS32(NULL,row,"DIRECTION") != 1)||
+              (psMetadataLookupS32(NULL,row,"POSITION") != j + 1)) {
+            continue;
+          }
+          psVectorSet(col_correction_fluxes,m,psMetadataLookupF32(NULL,row,"FLUX"));
+          psVectorSet(col_correction_factors,m,psMetadataLookupF32(NULL,row,"FACTOR"));
+          //      psTrace("psModules.camera",6,"columns: %d\n",m);
+          m++;
+        }
+      } // End load specific col
+      else {
+        col_correction_fluxes = default_col_correction_fluxes;
+        col_correction_factors = default_col_correction_factors;
+      } // End load default col
+
+      // Calculate correction factor contribution for this pixel.
+      psF32 factor_row = pmNonLinearityMeasure(image->data.F32[i][j],
+                                               row_correction_fluxes,row_correction_factors);
+      psF32 factor_col = pmNonLinearityMeasure(image->data.F32[i][j],
+                                               col_correction_fluxes,col_correction_factors);
+
+      if (((i <= 9)&&(j <= 9))||((i == 9)&&(j == 5))) { // Print out if we're looking at a test case.
+        psTrace("psModules.camera",6,"Linearity: %d %d %s %f %f %f %d %d\n",i,j,
+                psMetadataLookupStr(NULL,inputReadout->parent->concepts,"CELL.NAME"),
+                image->data.F32[i][j],factor_row,factor_col,numBorder,numSamples);
+        psTrace("psModules.camera",6,"Linearity: R: %d %d %d C: %d %d %d\n",
+                i,(i < numBorder),(image->numRows - i),
+                j,(j < numBorder),(image->numCols - j));
+
+/*      psTrace("psModules.camera",6,"Linearity: V: "); */
+/*      for (int k = 0; k < numSamples; k++) { */
+/*        psTrace("psModules.camera",6,"(%f %f) ", */
+/*                col_correction_fluxes->data.F32[k],col_correction_factors->data.F32[k]); */
+/*      } */
+/*      psTrace("psModules.camera",6,"\n"); */
+      } // End Test case
+
+      // Apply correction to image data
+      image->data.F32[i][j] = image->data.F32[i][j] * ( factor_row * factor_col ) ;
+
+      // Free correction if we allocated
+      if ((j < numBorder)||(image->numCols - j < numBorder)) {
+        psFree(col_correction_fluxes);
+        psFree(col_correction_factors);
+      }
+
+    } // End loop over columns
+
+    // Free correction if we allocated
+    if ((i < numBorder)||(image->numRows - i < numBorder)) {
+      psFree(row_correction_fluxes);
+      psFree(row_correction_factors);
+    }
+
+  } // End loop over rows
+
+  psFree(default_row_correction_fluxes);
+  psFree(default_row_correction_factors);
+  psFree(default_col_correction_fluxes);
+  psFree(default_col_correction_factors);
+  
+  return(true);
+}
+
+psF32 pmNonLinearityMeasure(psF32 flux, psVector *correction_fluxes, psVector *correction_factors) {
+  //  psS32 numPixels = 0;
+  psF32 result = 0;
+  psU32 bin = 0;
+
+  bin = correction_fluxes->n - 1;
+  if (flux < correction_fluxes->data.F32[0]) {
+    return(1.0);
+  }
+  if (flux > correction_fluxes->data.F32[bin]) {
+    return(1.0);
+  }
+  
+  for (int i = 0; i < correction_fluxes->n - 1; i++) {
+    if ((flux >= correction_fluxes->data.F32[i])&&
+        (flux <  correction_fluxes->data.F32[i+1])) {
+      result = correction_factors->data.F32[i] +
+        (flux - correction_fluxes->data.F32[i]) *
+        ((correction_factors->data.F32[i+1] - correction_factors->data.F32[i]) /
+         (correction_fluxes->data.F32[i+1] - correction_fluxes->data.F32[i]));
+      continue;
+    }
+  }
+
+/*   PS_BIN_FOR_VALUE(bin,correction_fluxes,flux); */
+/*   if ((bin < 0)||(bin > correction_fluxes->n)) { */
+/*     return(1.0); */
+/*   } */
+  
+/*   PS_BIN_INTERPOLATE(result,correction_fluxes,correction_factors,bin,flux); */
+  if (!isfinite(result)) {
+    result = 1.0;
+  }
+  if (result <= 0) {
+    result = 1.0;
+  }
+  return(result);
+}
+  
+  
+