Changeset 6448 for branches/rel10_ifa/psModules/src/detrend/pmNonLinear.c

Timestamp:

Feb 17, 2006, 7:13:42 AM (20 years ago)

Author:

magnier

Message:

bulk merge of eam_rel9_p0 onto this branch

File:

• branches/rel10_ifa/psModules/src/detrend/pmNonLinear.c (modified) (5 diffs)

branches/rel10_ifa/psModules/src/detrend/pmNonLinear.c

@@ +1 @@
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the
+// one that was being worked on.
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
 /** @file  pmNonLinear.c
  *
@@ -5 +10 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-05 20:49:40 $
+ *  @version $Revision: 1.5.12.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-17 17:13:41 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -22 +27 @@
 
 #include "pmNonLinear.h"
-#include "pmSubtractBias.h"
 
-// XXX: Remove, autoconf must be
-#define PS_WARN_PTR_NON_NULL(NAME) \
-if ((NAME) == NULL) { \
-    psLogMsg(__func__, PS_LOG_WARN, "WARNING: %s is NULL.", #NAME); \
-} \
 /******************************************************************************
 pmNonLinearityLookup(): This routine will take an pmReadout image as input
 and a 1-D polynomial.  For each pixel in the input image, the polynomial will
-be evaluated at that pixels value, and the image pixel will then be set
-to
+be evaluated at that pixels value, and the image pixel will then be set to
 that value.
-*****************************************************************************/
+ *****************************************************************************/
 
-pmReadout *pmNonLinearityPolynomial(
-    pmReadout *inputReadout,
-    const psPolynomial1D *input1DPoly)
+pmReadout *pmNonLinearityPolynomial(pmReadout *inputReadout,
+                                    const psPolynomial1D *input1DPoly)
 {
     PS_ASSERT_PTR_NON_NULL(inputReadout, NULL);
@@ -45 +42 @@
     PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL);
     PS_ASSERT_PTR_NON_NULL(input1DPoly, NULL);
-    PS_WARN_PTR_NON_NULL(inputReadout->parent);
-    if (inputReadout->parent != NULL) {
-        PS_WARN_PTR_NON_NULL(inputReadout->parent->concepts);
-    }
 
-    //
-    // Determine trimmed image from metadata.
-    //
-    psImage *trimmedImg = p_psDetermineTrimmedImage(inputReadout);
+    psS32 i;
+    psS32 j;
 
-    for (psS32 i=0;i<trimmedImg->numRows;i++) {
-        for (psS32 j=0;j<trimmedImg->numCols;j++) {
-            trimmedImg->data.F32[i][j] = psPolynomial1DEval(input1DPoly,
-                                         trimmedImg->data.F32[i][j]);
+    for (i=0;i<inputReadout->image->numRows;i++) {
+        for (j=0;j<inputReadout->image->numCols;j++) {
+            inputReadout->image->data.F32[i][j] = psPolynomial1DEval(input1DPoly, inputReadout->image->data.F32[i][j]);
         }
     }
@@ -71 +61 @@
 inFluxe, and the corresponding value in outFlux.  The image pixel will then
 be set to the value from outFlux.
- 
-XXX: Must assert that filename exists.  This should probably happen in
-the lookup files.
  *****************************************************************************/
-pmReadout *pmNonLinearityLookup(
-    pmReadout *inputReadout,
-    const char *filename
-)
+pmReadout *pmNonLinearityLookup(pmReadout *inputReadout,
+                                const psVector *inFlux,
+                                const psVector *outFlux)
 {
     PS_ASSERT_PTR_NON_NULL(inputReadout,NULL);
     PS_ASSERT_PTR_NON_NULL(inputReadout->image,NULL);
     PS_ASSERT_IMAGE_TYPE(inputReadout->image, PS_TYPE_F32, NULL);
-    PS_WARN_PTR_NON_NULL(inputReadout->parent);
-    if (inputReadout->parent != NULL) {
-        PS_WARN_PTR_NON_NULL(inputReadout->parent->concepts);
+    PS_ASSERT_PTR_NON_NULL(inFlux,NULL);
+    if (inFlux->n < 2) {
+        psError(PS_ERR_UNKNOWN,true, "pmNonLinearityLookup(): input vector less than 2 elements.  Returning inputReadout image.");
+        return(inputReadout);
     }
-    //
-    // Determine trimmed image from metadata.
-    //
-    psImage *trimmedImg = p_psDetermineTrimmedImage(inputReadout);
+    PS_ASSERT_PTR_NON_NULL(outFlux,NULL);
+    psS32 tableSize = inFlux->n;
+    if (inFlux->n != outFlux->n) {
+        tableSize = PS_MIN(inFlux->n, outFlux->n);
+        psLogMsg(__func__, PS_LOG_WARN,
+                 "WARNING: pmNonLinear.c: pmNonLinearityLookup(): input vectors have different sizes (%d, %d)\n", inFlux->n, outFlux->n);
+    }
+    PS_ASSERT_VECTOR_TYPE(inFlux, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_TYPE(outFlux, PS_TYPE_F32, NULL);
 
-    psLookupTable *tmpLT = psLookupTableAlloc(filename, "%f %f", 0);
-    psS32 numLines = psLookupTableRead(tmpLT);
+    psS32 i;
+    psS32 j;
+    psS32 binNum;
+    psScalar x;
     psS32 numPixels = 0;
-    if (numLines < 2) {
-        psLogMsg(__func__, PS_LOG_WARN,
-                 "WARNING: Lookup Table is too small.  Returning original pmReadout.\n");
-    } else {
-        for (psS32 i=0;i<trimmedImg->numRows;i++) {
-            for (psS32 j=0;j<trimmedImg->numCols;j++) {
-                psF64 tmpD = psLookupTableInterpolate(tmpLT, trimmedImg->data.F32[i][j], 1);
-                if (!isnan(tmpD)) {
-                    trimmedImg->data.F32[i][j] = tmpD;
-                } else {
-                    numPixels++;
-                }
+    psF32 slope;
+
+    x.type.type = PS_TYPE_F32;
+    for (i=0;i<inputReadout->image->numRows;i++) {
+        for (j=0;j<inputReadout->image->numCols;j++) {
+            x.data.F32 = inputReadout->image->data.F32[i][j];
+            binNum = p_psVectorBinDisect((psVector *)inFlux, &x);
+
+            if (binNum == -2) {
+                // We get here if x is below the table lookup range.
+                inputReadout->image->data.F32[i][j] = outFlux->data.F32[0];
+                numPixels++;
+
+            } else if (binNum == -1) {
+                // We get here if x is above the table lookup range.
+                inputReadout->image->data.F32[i][j] = outFlux->data.F32[tableSize-1];
+                numPixels++;
+
+            } else if (binNum < -2) {
+                // We get here if there was some other problem.
+                psError(PS_ERR_UNKNOWN,true, "pmNonLinearityLookup(): Could not perform p_psVectorBinDisect().  Returning inputReadout image.");
+                return(inputReadout);
+                numPixels++;
+            } else {
+                // Perform linear interpolation.
+                slope = (outFlux->data.F32[binNum+1] - outFlux->data.F32[binNum]) /
+                        (inFlux->data.F32[binNum+1]  - inFlux->data.F32[binNum]);
+                inputReadout->image->data.F32[i][j] = outFlux->data.F32[binNum] +
+                                                      ((x.data.F32 - inFlux->data.F32[binNum]) * slope);
             }
         }
-        if (numPixels > 0) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "WARNING: pmNonLinear.c: pmNonLinearityLookup(): %d pixels outside table.", numPixels);
-        }
     }
-
+    if (numPixels > 0) {
+        psLogMsg(__func__, PS_LOG_WARN,
+                 "WARNING: pmNonLinear.c: pmNonLinearityLookup(): %d pixels outside table.", numPixels);
+    }
     return(inputReadout);
 }