Changeset 29833 for trunk

Timestamp:

Nov 25, 2010, 9:45:07 PM (15 years ago)

Author:

watersc1

Message:

Merge of czw_branch/20100817 into trunk. This includes the non-linearity correction application code.

Location:

trunk

Files:

• . (modified) (1 prop)
• dbconfig (modified) (1 prop)
• ippScripts (modified) (1 prop)
• ippTasks (modified) (1 prop)
• ippTools (modified) (1 prop)
• ippconfig/gpc1/ppImage.config (modified) (1 diff)
• ippconfig/recipes/filerules-mef.mdc (modified) (1 diff)
• ippconfig/recipes/filerules-simple.mdc (modified) (1 diff)
• ippconfig/recipes/filerules-split.mdc (modified) (1 diff)
• ppImage/src/ppImage.h (modified) (1 diff)
• ppImage/src/ppImageArguments.c (modified) (2 diffs)
• ppImage/src/ppImageDefineFile.c (modified) (4 diffs)
• ppImage/src/ppImageDetrendFree.c (modified) (2 diffs)
• ppImage/src/ppImageDetrendNonLinear.c (modified) (2 diffs)
• ppImage/src/ppImageDetrendReadout.c (modified) (3 diffs)
• ppImage/src/ppImageLoop.c (modified) (1 diff)
• ppImage/src/ppImageOptions.c (modified) (3 diffs)
• ppImage/src/ppImageParseCamera.c (modified) (1 diff)
• psLib/src/types/psMetadata.c (modified) (1 diff)
• psModules/src/camera/pmFPAfile.c (modified) (2 diffs)
• psModules/src/camera/pmFPAfile.h (modified) (1 diff)
• psModules/src/camera/pmFPAfileFitsIO.c (modified) (1 diff)
• psModules/src/camera/pmFPAfileIO.c (modified) (6 diffs)
• psModules/src/detrend/pmBias.c (modified) (3 diffs)
• psModules/src/detrend/pmBias.h (modified) (1 diff)
• psModules/src/detrend/pmDetrendDB.c (modified) (1 diff)
• psModules/src/detrend/pmDetrendDB.h (modified) (1 diff)
• psModules/src/detrend/pmNonLinear.c (modified) (3 diffs)
• psModules/src/detrend/pmNonLinear.h (modified) (1 diff)

trunk/ippconfig/gpc1/ppImage.config

@@ -262 +262 @@
   NOISEMAP METADATA
   END   
+  LINEARITY METADATA
+  END
 END
 

trunk/ippconfig/recipes/filerules-mef.mdc

@@ -63 +63 @@
 PPIMAGE.FRINGE          INPUT    @DETDB        CHIP       FRINGE
 PPIMAGE.SHUTTER         INPUT    @DETDB        CHIP       IMAGE     
+PPIMAGE.LINEARITY       INPUT    @DETDB        CHIP       LINEARITY
 
 ## Files used by ppMerge

trunk/ippconfig/recipes/filerules-simple.mdc

@@ -28 +28 @@
 PPIMAGE.FRINGE            INPUT    @DETDB        CHIP       FRINGE
 PPIMAGE.SHUTTER           INPUT    @DETDB        CHIP       IMAGE     
+PPIMAGE.LINEARITY         INPUT    @DETDB        CHIP       LINEARITY
                           
 ## Files used by ppMerge  

trunk/ippconfig/recipes/filerules-split.mdc

@@ -46 +46 @@
 PPIMAGE.FRINGE            INPUT    @DETDB        CHIP       FRINGE
 PPIMAGE.SHUTTER           INPUT    @DETDB        CHIP       IMAGE
+PPIMAGE.LINEARITY         INPUT    @DETDB        CHIP       LINEARITY
 
 ## Files used by ppMerge

trunk/ppImage/src/ppImage.h

@@ -155 +155 @@
 bool ppImageDetrendBias(pmReadout *inputReadout, pmReadout *bias, pmReadout *dark, ppImageOptions *options);
 
-bool ppImageDetrendNonLinear(pmReadout *input, ppImageOptions *options);
+//bool ppImageDetrendNonLinear(pmReadout *input, ppImageOptions *options);
+bool ppImageDetrendNonLinear(pmReadout *input, pmFPAview *linearity, pmConfig *config);
 bool ppImageDetrendNonLinearLookup(pmReadout *input, psMetadataItem *dataItem);
 bool ppImageDetrendNonLinearPolynomial(pmReadout *input, psMetadataItem *dataItem);

trunk/ppImage/src/ppImageArguments.c

@@ -24 +24 @@
     fprintf(stderr, "\t-mask/-masklist: Mask image.\n");
     fprintf(stderr, "\t-fringe/-fringelist: Fringe image and data.\n");
+    fprintf(stderr, "\t-linearity/-linearlist: linearity correction file.\n");
     fprintf(stderr, "\n");
     exit (2);
@@ -120 +121 @@
     pmConfigFileSetsMD (config->arguments, &argc, argv, "MASK", "-mask", "-masklist");
     pmConfigFileSetsMD (config->arguments, &argc, argv, "FRINGE", "-fringe", "-fringelist");
+    pmConfigFileSetsMD (config->arguments, &argc, argv, "LINEARITY", "-linearity", "-linearlist");
 
     // chip selection is used to limit chips to be processed

trunk/ppImage/src/ppImageDefineFile.c

@@ -15 +15 @@
         file = pmFPAfileDefineFromArgs(&status, config, filerule, argname);
         if (!status) {
-            psError(PS_ERR_UNKNOWN, false, "failed to load file definition");
+            psError(PS_ERR_UNKNOWN, false, "failed to load file definition ARG LIST");
             return false;
         }
@@ -23 +23 @@
         file = pmFPAfileDefineFromRun(&status, NULL, config, filerule);
         if (!status) {
-            psError(PS_ERR_UNKNOWN, false, "failed to load file definition");
+            psError(PS_ERR_UNKNOWN, false, "failed to load file definition RUN");
             return false;
         }
@@ -31 +31 @@
         file = pmFPAfileDefineFromConf(&status, config, filerule);
         if (!status) {
-            psError(PS_ERR_UNKNOWN, false, "failed to load file definition");
+            psError(PS_ERR_UNKNOWN, false, "failed to load file definition CONFIG");
             return false;
         }
@@ -39 +39 @@
         file = pmFPAfileDefineFromDetDB(&status, config, filerule, input, detrendType);
         if (!status) {
-            psError(PS_ERR_UNKNOWN, false, "failed to load file definition");
+            psError(PS_ERR_UNKNOWN, false, "failed to load file definition DETREND");
             return false;
         }

trunk/ppImage/src/ppImageDetrendFree.c

@@ -13 +13 @@
     "PPIMAGE.FLAT",
     "PPIMAGE.SHUTTER",
+    "PPIMAGE.LINEARITY",
     NULL
 };
@@ -23 +24 @@
 
         pmFPAfile *file = psMetadataLookupPtr(&status, config->files, detrendTypes[i]); // File of interest
+        psTrace("pmFPAfileFree",1,"Working on %s\n",detrendTypes[i]);
         if (!file) continue; // not all detrends are used in any given run
 

trunk/ppImage/src/ppImageDetrendNonLinear.c

@@ -44 +44 @@
 }
 
-bool ppImageDetrendNonLinear(pmReadout *input, ppImageOptions *options) {
+
+bool ppImageDetrendNonLinear(pmReadout *input, pmFPAview *detview, pmConfig  *config) {
+    bool status;
+
+    pmFPAfile *linearity_file = psMetadataLookupPtr(&status,config->files,"PPIMAGE.LINEARITY");
+    psFits *linearity_fits = linearity_file->fits;
+
+    char *extname = psMetadataLookupStr(&status,input->parent->concepts,"CELL.NAME");
+    if (!extname) {
+        psError(PS_ERR_IO, false, "missing CELL.NAME in concepts");
+        return(false);
+    }
+
+    if (!psFitsMoveExtName(linearity_fits,extname)) {
+        psError(PS_ERR_IO, false, "Unable to move to non-linearity table %s", extname);
+        return(false);
+    }
+  
+    psArray *table = psFitsReadTable(linearity_fits);
+    if (!table) {
+        psError(PS_ERR_IO, false, "Unable to read non-linearity table.\n");
+        return(false);
+    }
+
+    // It might be better to pack lookup table here...
+    // Why? I only use that lookup table once for the single cell it matches. 
+  
+    if (!pmNonLinearityApply(input,table)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to apply non-linearity corrections.\n");
+        psFree (table);
+        return(false);
+    }       
+    psFree (table);
+
+    return true;
+}
+
+bool ppImageDetrendNonLinear_Original(pmReadout *input, ppImageOptions *options) {
 
     psMetadataItem *concept;
@@ -59 +96 @@
 
       case PS_DATA_METADATA:
-        // XXX EAM: this is somewhat confusing : let's wrap in a function when i understand it
-
         // Go looking for the value in the hierarchy
         concept = psMetadataLookup(cell->concepts, options->nonLinearSource);

trunk/ppImage/src/ppImageDetrendReadout.c

@@ -51 +51 @@
     }
 
-
-# if 0
+    // Subtract the overscan
+    if (options->doOverscan) {
+      if (!pmOverscanSubtract (input, options->overscan)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to subtract overscan.");
+        psFree(detview);
+        return false;
+      }
+    }
+
     // Non-linearity correction
     if (options->doNonLin) {
-        ppImageDetrendNonLinear(detrend->input, input, options);
-    }
-# endif
+      //      linearity = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.LINEARITY");
+      if (!ppImageDetrendNonLinear(input,detview,config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to correct NonLinearity");
+        psFree(detview);
+        return(false);
+      }
+      /*       ppImageDetrendNonLinear(detrend->input, input, options); */
+    }
 
     // set up the dark and bias
@@ -77 +89 @@
     }
 
-    // Bias, dark and overscan subtraction are all merged.
-    if (options->doBias || options->doOverscan) {
-        if (!pmBiasSubtract(input, options->overscan, bias, oldDark, view)) {
+    // Bias and temperature-independent-dark subtraction are merged.
+    if (options->doBias) {
+        if (!pmBiasSubtract(input, bias, oldDark, view)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to subtract bias.");
             psFree(detview);
@@ -85 +97 @@
         }
     }
-
+    
     // Weight on the basis of pixel value needs to be done after the overscan has been subtracted
     if (options->doVarianceBuild) {

trunk/ppImage/src/ppImageLoop.c

@@ -152 +152 @@
             ESCAPE("Unable to free detrend images");
         }
-
+    
         // Apply the fringe correction
         if (options->doFringe) {

trunk/ppImage/src/ppImageOptions.c

@@ -8 +8 @@
 {
     psFree(options->overscan);
-    psFree(options->nonLinearData);
-    psFree(options->nonLinearSource);
+    // psFree(options->nonLinearData);
+    // psFree(options->nonLinearSource);
 }
 
@@ -130 +130 @@
         psMetadataItem *dataItem = psMetadataLookup(recipe, "NONLIN.DATA");
         if (! dataItem) {
-            psLogMsg(__func__, PS_LOG_ERROR, "Non-linearity correction desired, but unable to "
-                     "find NONLIN.DATA in recipe %s.", RECIPE_NAME);
+            psLogMsg("ppImage", PS_LOG_ERROR, "Non-linearity correction desired, but unable to find NONLIN.DATA in recipe %s.", RECIPE_NAME);
             exit(EXIT_FAILURE);
         }
@@ -147 +146 @@
             // This is a menu; we need the key
           case PS_DATA_METADATA:
-            {
-                bool status;
-                options->nonLinearSource = psMetadataLookupStr(&status, recipe, "NONLIN.SOURCE");
-                if (! status || ! options->nonLinearSource) {
-                    psLogMsg(__func__, PS_LOG_ERROR, "Non-linearity correction desired, but unable to "
-                            "find NONLIN.SOURCE in recipe %s.", RECIPE_NAME);
-                    exit(EXIT_FAILURE);
-                }
-            }
+            options->nonLinearSource = psMetadataLookupStr(&status, recipe, "NONLIN.SOURCE");
+            if (! status || ! options->nonLinearSource) {
+                psLogMsg("ppImage", PS_LOG_ERROR, "Non-linearity correction desired, but unable to find NONLIN.SOURCE in recipe %s.", RECIPE_NAME);
+                exit(EXIT_FAILURE);
+            }
             break;
           default:
-            psLogMsg(__func__, PS_LOG_ERROR, "Non-linearity correction desired, but "
-                    "NONLIN.DATA is of invalid type in recipe %s.", RECIPE_NAME);
+            psLogMsg("ppImage", PS_LOG_ERROR, "Non-linearity correction desired, but NONLIN.DATA is of invalid type in recipe %s.", RECIPE_NAME);
             exit(EXIT_FAILURE);
         }

trunk/ppImage/src/ppImageParseCamera.c

@@ -73 +73 @@
             return NULL;
         }
+    }
+
+    if (options->doNonLin) {
+      if (!ppImageDefineFile(config, input->fpa, "PPIMAGE.LINEARITY", "LINEARITY",
+                             PM_FPA_FILE_LINEARITY, PM_DETREND_TYPE_LINEARITY)) {
+        psError(PS_ERR_IO, false, "Can't find a non-linearity correction source");
+        psFree(options);
+        return NULL;
+      }
     }
 

trunk/psLib/src/types/psMetadata.c

@@ -109 +109 @@
     type = metadataItem->type;
 
-
+    //    fprintf(stderr,"%s %s %d\n",metadataItem->name,metadataItem->comment,type);
     psMemDecrRefCounter(metadataItem->name);
     psMemDecrRefCounter(metadataItem->comment);

trunk/psModules/src/camera/pmFPAfile.c

@@ -36 +36 @@
         return;
     }
-
+    psTrace ("pmFPAfileFree", 5, "freeing %s %ld\n", file->name,(psU64) file->fits);
     psAssert(!fpaFileFreeStrict || file->fits == NULL, "File %s wasn't closed.", file->name);
 
@@ -523 +523 @@
     if (!strcasecmp(type, "HEADER"))     {
         return PM_FPA_FILE_HEADER;
+    }
+    if (!strcasecmp(type, "LINEARITY"))  {
+      return PM_FPA_FILE_LINEARITY;
     }
     // deprecate this?

trunk/psModules/src/camera/pmFPAfile.h

@@ -50 +50 @@
     PM_FPA_FILE_SRCTEXT,
     PM_FPA_FILE_PATTERN,
+    PM_FPA_FILE_LINEARITY,
 } pmFPAfileType;
 

trunk/psModules/src/camera/pmFPAfileFitsIO.c

@@ -145 +145 @@
       case PM_FPA_FILE_FRINGE:
       case PM_FPA_FILE_DARK:
+      case PM_FPA_FILE_LINEARITY:
       case PM_FPA_FILE_CMP:
       case PM_FPA_FILE_CMF:

trunk/psModules/src/camera/pmFPAfileIO.c

@@ -233 +233 @@
       case PM_FPA_FILE_JPEG:
       case PM_FPA_FILE_KAPA:
+      case PM_FPA_FILE_LINEARITY:
         break;
       default:
@@ -560 +561 @@
       case PM_FPA_FILE_ASTROM_MODEL:
       case PM_FPA_FILE_ASTROM_REFSTARS:
+      case PM_FPA_FILE_LINEARITY:
         psTrace ("psModules.camera", 5, "closing %s (%s) (%d:%d:%d)\n", file->filename, file->name, view->chip, view->cell, view->readout);
         status = psFitsClose (file->fits);
@@ -575 +577 @@
       case PM_FPA_FILE_JPEG:
       case PM_FPA_FILE_KAPA:
+
         break;
 
@@ -619 +622 @@
       case PM_FPA_FILE_FRINGE:
       case PM_FPA_FILE_DARK:
+        //
         status = pmFPAviewFreeData(view, file);
         break;
@@ -636 +640 @@
       case PM_FPA_FILE_JPEG:
       case PM_FPA_FILE_KAPA:
+      case PM_FPA_FILE_LINEARITY:
         psTrace ("psModules.camera", 5, "nothing to free for %s (%s)\n", file->filename, file->name);
         return true;
@@ -791 +796 @@
       case PM_FPA_FILE_ASTROM_MODEL:
       case PM_FPA_FILE_ASTROM_REFSTARS:
+      case PM_FPA_FILE_LINEARITY:
         psTrace ("psModules.camera", 5, "opening %s (%s) (%d:%d:%d)\n",
                  file->filename, file->name, view->chip, view->cell, view->readout);

trunk/psModules/src/detrend/pmBias.c

@@ -163 +163 @@
 }
 
-bool pmBiasSubtract(pmReadout *in, pmOverscanOptions *overscanOpts,
+bool pmBiasSubtract(pmReadout *in, 
                     pmReadout *bias, pmReadout *dark, const pmFPAview *view)
 {
@@ -184 +184 @@
 
     pmHDU *hdu = pmHDUFromReadout(in);  // HDU of interest
-
-    if (!pmOverscanSubtract (in, overscanOpts)) {
-        return false;
-    }
 
     // Bias frame subtraction
@@ -247 +243 @@
     psString timeString = psTimeToISO(time); // String with time
     psFree(time);
-    psStringPrepend(&timeString, "Overscan/bias/dark processing completed at ");
+    psStringPrepend(&timeString, "Bias/dark processing completed at ");
     psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
                      timeString, "");

trunk/psModules/src/detrend/pmBias.h

@@ -21 +21 @@
 /// bias (if non-NULL) and dark (if non-NULL) scaled by the CELL.DARKTIME concept.
 bool pmBiasSubtract(pmReadout *in,      ///< Input readout, to be overscan/bias/dark corrected
-                    pmOverscanOptions *overscanOpts, ///< Options for overscan subtraction, or NULL
                     pmReadout *bias, ///< Bias to subtract, or NULL
                     pmReadout *dark, ///< Dark to scale and subtract, or NULL

trunk/psModules/src/detrend/pmDetrendDB.c

@@ -106 +106 @@
         DETREND_STRING_CASE(ASTROM);
         DETREND_STRING_CASE(NOISEMAP);
+        DETREND_STRING_CASE(LINEARITY);
     default:
         return NULL;

trunk/psModules/src/detrend/pmDetrendDB.h

@@ -36 +36 @@
     PM_DETREND_TYPE_ASTROM,
     PM_DETREND_TYPE_NOISEMAP,
+    PM_DETREND_TYPE_LINEARITY,
 } pmDetrendType;
 

trunk/psModules/src/detrend/pmNonLinear.c

@@ -10 +10 @@
 #include "pmNonLinear.h"
 
+psF32 pmNonLinearityMeasureNoisy(psF32 flux, psVector *correction_fluxes, psVector *correction_factors);
 pmReadout *pmNonLinearityPolynomial(pmReadout *inputReadout, const psPolynomial1D *input1DPoly)
 {
@@ -27 +28 @@
 
 // set the bin closest to the corresponding value.  
-#define PS_BIN_FOR_VALUE(RESULT, VECTOR, VALUE) { \
-        psVectorBinaryDisectResult result; \
-        psScalar tmpScalar; \
-        tmpScalar.type.type = PS_TYPE_F32; \
-        tmpScalar.data.F32 = (VALUE); \
-        RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar); \
-        switch (result) { \
-          case PS_BINARY_DISECT_PASS: \
-            break; \
-          case PS_BINARY_DISECT_OUTSIDE_RANGE: \
-            numPixels ++; \
-            break; \
-          case PS_BINARY_DISECT_INVALID_INPUT: \
-          case PS_BINARY_DISECT_INVALID_TYPE: \
-            psAbort ("programming error"); \
-            break; \
+#define PS_BIN_FOR_VALUE(RESULT, VECTOR, VALUE) {                       \
+        psVectorBinaryDisectResult result;                              \
+        psScalar tmpScalar;                                             \
+        tmpScalar.type.type = PS_TYPE_F32;                              \
+        tmpScalar.data.F32 = (VALUE);                                   \
+        RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar);    \
+        switch (result) {                                               \
+          case PS_BINARY_DISECT_PASS:                                   \
+            break;                                                      \
+          case PS_BINARY_DISECT_OUTSIDE_RANGE:                          \
+            numPixels ++;                                               \
+            break;                                                      \
+          case PS_BINARY_DISECT_INVALID_INPUT:                          \
+          case PS_BINARY_DISECT_INVALID_TYPE:                           \
+            psAbort ("programming error");                              \
+            break;                                                      \
         } }
+
+
+# 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 +117 @@
     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);
+
+    psTimerStart ("nonlinear");
+
+    psS32 numSamples = 39;
+    psS32 numBorder  = 10;
+
+    //  psS32 tableSize = Ltab->n;
+
+    psImage *image = inputReadout->image;
+
+    // Load default data.
+    psVector *default_row_correction_fluxes = psVectorAlloc(numSamples,PS_TYPE_F32);
+    psVector *default_col_correction_fluxes = psVectorAlloc(numSamples,PS_TYPE_F32);
+
+    psVector *default_row_correction_factors = psVectorAlloc(numSamples,PS_TYPE_F32);
+    psVector *default_col_correction_factors = psVectorAlloc(numSamples,PS_TYPE_F32);
+
+    // pre-allocate the correction vectors
+    psVector *row_correction_fluxes  = NULL;
+    psVector *col_correction_fluxes  = NULL;
+    psVector *row_correction_factors = NULL;
+    psVector *col_correction_factors = NULL;
+
+    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
+  
+    psLogMsg ("psModules", PS_LOG_MINUTIA, "load default data from table: %f sec\n", psTimerMark ("nonlinear"));
+
+    // pre-allocate arrays with the correction vectors for the borders
+    psArray *x_lo_flux = psArrayAlloc(numBorder);
+    psArray *x_hi_flux = psArrayAlloc(numBorder);
+    psArray *y_lo_flux = psArrayAlloc(numBorder);
+    psArray *y_hi_flux = psArrayAlloc(numBorder);
+
+    psArray *x_lo_fact = psArrayAlloc(numBorder);
+    psArray *x_hi_fact = psArrayAlloc(numBorder);
+    psArray *y_lo_fact = psArrayAlloc(numBorder);
+    psArray *y_hi_fact = psArrayAlloc(numBorder);
+    for (int i = 0; i < numBorder; i++) {
+        // pre-allocate the correction vectors
+        x_lo_flux->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
+        x_hi_flux->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
+        y_lo_flux->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
+        y_hi_flux->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
+
+        x_lo_fact->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
+        x_hi_fact->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
+        y_lo_fact->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
+        y_hi_fact->data[i] = psVectorAllocEmpty(numSamples,PS_TYPE_F32);
+    }   
+
+    // parse out the full table:
+    for (int k = 0; k < Ltab->n; k++) {
+        psMetadata *row = Ltab->data[k];        
+        int dir = psMetadataLookupS32(NULL,row,"DIRECTION");
+        int pos = psMetadataLookupS32(NULL,row,"POSITION");
+
+        psVector *fluxVector = NULL;
+        psVector *factVector = NULL;
+
+        int seq = -1;
+        if ((dir == 0) && (pos < image->numCols/2)) {
+          seq = pos ;
+            if (seq < 0) continue;
+            fluxVector = x_lo_flux->data[seq];
+            factVector = x_lo_fact->data[seq];
+        } 
+        if ((dir == 0) && (pos > image->numCols/2)) {
+            seq = pos + numBorder - image->numCols;
+            if (seq >= image->numCols) continue;
+            fluxVector = x_hi_flux->data[seq];
+            factVector = x_hi_fact->data[seq];
+        }
+        if ((dir == 1) && (pos < image->numRows/2)) {
+            seq = pos;
+            if (seq < 0) continue;
+            fluxVector = y_lo_flux->data[seq];
+            factVector = y_lo_fact->data[seq];
+        } 
+        if ((dir == 1) && (pos > image->numRows/2)) {
+            seq = pos + numBorder - image->numRows;
+            if (seq >= image->numRows) continue;
+            fluxVector = y_hi_flux->data[seq];
+            factVector = y_hi_fact->data[seq];
+        }
+
+        float flux = psMetadataLookupF32(NULL,row,"FLUX");
+        float factor = psMetadataLookupF32(NULL,row,"FACTOR");
+
+        psVectorAppend(fluxVector, flux);
+        psVectorAppend(factVector, factor);
+    }
+    psLogMsg ("psModules", PS_LOG_MINUTIA, "load border data from table: %f sec\n", psTimerMark ("nonlinear"));
+
+    for (int i = 0; i < image->numRows; i++) { // Loop over rows : note: problem with discontinuity here
+        row_correction_fluxes = NULL;
+        if (i < numBorder) {
+            row_correction_fluxes  = y_lo_flux->data[i];
+            row_correction_factors = y_lo_fact->data[i];
+        }
+        if (i > image->numRows - numBorder) {
+            row_correction_fluxes  = y_hi_flux->data[i + numBorder - image->numRows];
+            row_correction_factors = y_hi_fact->data[i + numBorder - image->numRows];
+        }
+        if (row_correction_fluxes == NULL) {
+          row_correction_factors = default_row_correction_factors;
+            row_correction_fluxes = default_row_correction_fluxes;
+        }
+
+        for (int j = 0; j < image->numCols; j++) { // Loop over columns
+            col_correction_fluxes = NULL;
+            if (j < numBorder) {
+                col_correction_fluxes  = x_lo_flux->data[j];
+                col_correction_factors = x_lo_fact->data[j];
+            }
+            if (j > image->numCols - numBorder) {
+                col_correction_fluxes  = x_hi_flux->data[j + numBorder - image->numCols];
+                col_correction_factors = x_hi_fact->data[j + numBorder - image->numCols];
+            }
+            if (col_correction_fluxes == NULL) {
+              col_correction_factors = default_col_correction_factors;
+              col_correction_fluxes = default_col_correction_fluxes;
+            }
+
+            // 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 (0)
+            if (((i == 200)&&(j == 200))||((i == 9)&&(j == 5))) { // Print out if we're looking at a test case.
+              psF32 factor_row = pmNonLinearityMeasureNoisy(image->data.F32[i][j], row_correction_fluxes,row_correction_factors);
+              psF32 factor_col = pmNonLinearityMeasureNoisy(image->data.F32[i][j], col_correction_fluxes,col_correction_factors);
+              
+                psTrace("psModules.nonlin",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.nonlin",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.nonlin",6,"Linearity: V: ");
+                for (int k = 0; k < numSamples; k++) {
+                    psTrace("psModules.nonlin",6,"(%f %f) (%f %f) DDDD> (%f %f) (%f %f)",
+                            col_correction_fluxes->data.F32[k],col_correction_factors->data.F32[k],
+                            row_correction_fluxes->data.F32[k],row_correction_factors->data.F32[k],
+                            default_col_correction_fluxes->data.F32[k],default_col_correction_factors->data.F32[k],
+                            default_row_correction_fluxes->data.F32[k],default_row_correction_factors->data.F32[k]);
+                }
+                psTrace("psModules.nonlin",6,"\n");
+            } // End Test case
+#endif
+            // Apply correction to image data
+#if (0)
+            if (((i == 200)&&(j == 200))||((i == 9)&&(j == 5))) { // Print out if we're looking at a test case.
+              psTrace("psModules.nonlin",4,"Applied Linearity Correction: %s %d %d : %f %f -> %f\n",
+                      psMetadataLookupStr(NULL,inputReadout->parent->concepts,"CELL.NAME"),
+                      i,j,image->data.F32[i][j],(factor_row + factor_col) / 2.0,
+                      image->data.F32[i][j] + (factor_row + factor_col) / 2.0);
+            }
+# endif
+            image->data.F32[i][j] = image->data.F32[i][j] - ( factor_row + factor_col ) / 2.0;
+
+        } // End loop over columns
+    } // End loop over rows
+    psLogMsg ("psModules", PS_LOG_MINUTIA, "apply correction: %f sec\n", psTimerMark ("nonlinear"));
+
+    psFree(x_lo_flux);
+    psFree(x_hi_flux);
+    psFree(y_lo_flux);
+    psFree(y_hi_flux);
+
+    psFree(x_lo_fact);
+    psFree(x_hi_fact);
+    psFree(y_lo_fact);
+    psFree(y_hi_fact);
+
+    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(0.0);
+    }
+/*     if (flux > correction_fluxes->data.F32[bin]) { */
+/*      return(0.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 = 0.0;
+    }
+/*     if (result <= 0) { */
+/*      result = 1.0; */
+/*     } */
+    return(result);
+}
+
+psF32 pmNonLinearityMeasureNoisy(psF32 flux, psVector *correction_fluxes, psVector *correction_factors) {
+    //  psS32 numPixels = 0;
+    psF32 result = 0;
+    psU32 bin = 0;
+
+    bin = correction_fluxes->n - 1;
+    psTrace("psModules.nonlin",6,"NLMN: %f %d %f %f\n",flux,bin,correction_fluxes->data.F32[0],correction_fluxes->data.F32[bin]);
+    if (flux < correction_fluxes->data.F32[0]) {
+        return(0.0);
+    }
+/*     if (flux > correction_fluxes->data.F32[bin]) { */
+/*      return(0.0); */
+/*     } */
+
+    for (int i = 0; i < correction_fluxes->n - 1; i++) {
+      psTrace("psModules.nonlin",6,"NLMN: %f %d %f %f %f %f\n",flux,i,correction_fluxes->data.F32[i],correction_fluxes->data.F32[i],
+              correction_factors->data.F32[i],correction_factors->data.F32[i+1]);
+        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 = 0.0;
+    }
+/*     if (result <= 0) { */
+/*      result = 1.0; */
+/*     } */
+    return(result);
+}
+
+  
+  

trunk/psModules/src/detrend/pmNonLinear.h

@@ -31 +31 @@
                                 const psVector *outFlux ///< Table column with output fluxes
                                );
-
+bool pmNonLinearityApply(pmReadout *inputReadout, psArray *Ltab);
+psF32 pmNonLinearityMeasure(psF32 flux, psVector *correction_fluxes, psVector *correction_factors);
 /// @}
 #endif