Index: trunk/ippconfig/gpc1/ppImage.config =================================================================== --- trunk/ippconfig/gpc1/ppImage.config (revision 29831) +++ trunk/ippconfig/gpc1/ppImage.config (revision 29833) @@ -262,4 +262,6 @@ NOISEMAP METADATA END + LINEARITY METADATA + END END Index: trunk/ippconfig/recipes/filerules-mef.mdc =================================================================== --- trunk/ippconfig/recipes/filerules-mef.mdc (revision 29831) +++ trunk/ippconfig/recipes/filerules-mef.mdc (revision 29833) @@ -63,4 +63,5 @@ PPIMAGE.FRINGE INPUT @DETDB CHIP FRINGE PPIMAGE.SHUTTER INPUT @DETDB CHIP IMAGE +PPIMAGE.LINEARITY INPUT @DETDB CHIP LINEARITY ## Files used by ppMerge Index: trunk/ippconfig/recipes/filerules-simple.mdc =================================================================== --- trunk/ippconfig/recipes/filerules-simple.mdc (revision 29831) +++ trunk/ippconfig/recipes/filerules-simple.mdc (revision 29833) @@ -28,4 +28,5 @@ PPIMAGE.FRINGE INPUT @DETDB CHIP FRINGE PPIMAGE.SHUTTER INPUT @DETDB CHIP IMAGE +PPIMAGE.LINEARITY INPUT @DETDB CHIP LINEARITY ## Files used by ppMerge Index: trunk/ippconfig/recipes/filerules-split.mdc =================================================================== --- trunk/ippconfig/recipes/filerules-split.mdc (revision 29831) +++ trunk/ippconfig/recipes/filerules-split.mdc (revision 29833) @@ -46,4 +46,5 @@ PPIMAGE.FRINGE INPUT @DETDB CHIP FRINGE PPIMAGE.SHUTTER INPUT @DETDB CHIP IMAGE +PPIMAGE.LINEARITY INPUT @DETDB CHIP LINEARITY ## Files used by ppMerge Index: trunk/ppImage/src/ppImage.h =================================================================== --- trunk/ppImage/src/ppImage.h (revision 29831) +++ trunk/ppImage/src/ppImage.h (revision 29833) @@ -155,5 +155,6 @@ 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); Index: trunk/ppImage/src/ppImageArguments.c =================================================================== --- trunk/ppImage/src/ppImageArguments.c (revision 29831) +++ trunk/ppImage/src/ppImageArguments.c (revision 29833) @@ -24,4 +24,5 @@ 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,4 +121,5 @@ 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 Index: trunk/ppImage/src/ppImageDefineFile.c =================================================================== --- trunk/ppImage/src/ppImageDefineFile.c (revision 29831) +++ trunk/ppImage/src/ppImageDefineFile.c (revision 29833) @@ -15,5 +15,5 @@ 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,5 +23,5 @@ 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,5 +31,5 @@ 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,5 +39,5 @@ 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; } Index: trunk/ppImage/src/ppImageDetrendFree.c =================================================================== --- trunk/ppImage/src/ppImageDetrendFree.c (revision 29831) +++ trunk/ppImage/src/ppImageDetrendFree.c (revision 29833) @@ -13,4 +13,5 @@ "PPIMAGE.FLAT", "PPIMAGE.SHUTTER", + "PPIMAGE.LINEARITY", NULL }; @@ -23,4 +24,5 @@ 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 Index: trunk/ppImage/src/ppImageDetrendNonLinear.c =================================================================== --- trunk/ppImage/src/ppImageDetrendNonLinear.c (revision 29831) +++ trunk/ppImage/src/ppImageDetrendNonLinear.c (revision 29833) @@ -44,5 +44,42 @@ } -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,6 +96,4 @@ 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); Index: trunk/ppImage/src/ppImageDetrendReadout.c =================================================================== --- trunk/ppImage/src/ppImageDetrendReadout.c (revision 29831) +++ trunk/ppImage/src/ppImageDetrendReadout.c (revision 29833) @@ -51,11 +51,23 @@ } - -# 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,7 +89,7 @@ } - // 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,5 +97,5 @@ } } - + // Weight on the basis of pixel value needs to be done after the overscan has been subtracted if (options->doVarianceBuild) { Index: trunk/ppImage/src/ppImageLoop.c =================================================================== --- trunk/ppImage/src/ppImageLoop.c (revision 29831) +++ trunk/ppImage/src/ppImageLoop.c (revision 29833) @@ -152,5 +152,5 @@ ESCAPE("Unable to free detrend images"); } - + // Apply the fringe correction if (options->doFringe) { Index: trunk/ppImage/src/ppImageOptions.c =================================================================== --- trunk/ppImage/src/ppImageOptions.c (revision 29831) +++ trunk/ppImage/src/ppImageOptions.c (revision 29833) @@ -8,6 +8,6 @@ { psFree(options->overscan); - psFree(options->nonLinearData); - psFree(options->nonLinearSource); + // psFree(options->nonLinearData); + // psFree(options->nonLinearSource); } @@ -130,6 +130,5 @@ 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,17 +146,12 @@ // 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); } Index: trunk/ppImage/src/ppImageParseCamera.c =================================================================== --- trunk/ppImage/src/ppImageParseCamera.c (revision 29831) +++ trunk/ppImage/src/ppImageParseCamera.c (revision 29833) @@ -73,4 +73,13 @@ 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; + } } Index: trunk/psLib/src/types/psMetadata.c =================================================================== --- trunk/psLib/src/types/psMetadata.c (revision 29831) +++ trunk/psLib/src/types/psMetadata.c (revision 29833) @@ -109,5 +109,5 @@ type = metadataItem->type; - + // fprintf(stderr,"%s %s %d\n",metadataItem->name,metadataItem->comment,type); psMemDecrRefCounter(metadataItem->name); psMemDecrRefCounter(metadataItem->comment); Index: trunk/psModules/src/camera/pmFPAfile.c =================================================================== --- trunk/psModules/src/camera/pmFPAfile.c (revision 29831) +++ trunk/psModules/src/camera/pmFPAfile.c (revision 29833) @@ -36,5 +36,5 @@ 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,4 +523,7 @@ if (!strcasecmp(type, "HEADER")) { return PM_FPA_FILE_HEADER; + } + if (!strcasecmp(type, "LINEARITY")) { + return PM_FPA_FILE_LINEARITY; } // deprecate this? Index: trunk/psModules/src/camera/pmFPAfile.h =================================================================== --- trunk/psModules/src/camera/pmFPAfile.h (revision 29831) +++ trunk/psModules/src/camera/pmFPAfile.h (revision 29833) @@ -50,4 +50,5 @@ PM_FPA_FILE_SRCTEXT, PM_FPA_FILE_PATTERN, + PM_FPA_FILE_LINEARITY, } pmFPAfileType; Index: trunk/psModules/src/camera/pmFPAfileFitsIO.c =================================================================== --- trunk/psModules/src/camera/pmFPAfileFitsIO.c (revision 29831) +++ trunk/psModules/src/camera/pmFPAfileFitsIO.c (revision 29833) @@ -145,4 +145,5 @@ 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: Index: trunk/psModules/src/camera/pmFPAfileIO.c =================================================================== --- trunk/psModules/src/camera/pmFPAfileIO.c (revision 29831) +++ trunk/psModules/src/camera/pmFPAfileIO.c (revision 29833) @@ -233,4 +233,5 @@ case PM_FPA_FILE_JPEG: case PM_FPA_FILE_KAPA: + case PM_FPA_FILE_LINEARITY: break; default: @@ -560,4 +561,5 @@ 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,4 +577,5 @@ case PM_FPA_FILE_JPEG: case PM_FPA_FILE_KAPA: + break; @@ -619,4 +622,5 @@ case PM_FPA_FILE_FRINGE: case PM_FPA_FILE_DARK: + // status = pmFPAviewFreeData(view, file); break; @@ -636,4 +640,5 @@ 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,4 +796,5 @@ 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); Index: trunk/psModules/src/detrend/pmBias.c =================================================================== --- trunk/psModules/src/detrend/pmBias.c (revision 29831) +++ trunk/psModules/src/detrend/pmBias.c (revision 29833) @@ -163,5 +163,5 @@ } -bool pmBiasSubtract(pmReadout *in, pmOverscanOptions *overscanOpts, +bool pmBiasSubtract(pmReadout *in, pmReadout *bias, pmReadout *dark, const pmFPAview *view) { @@ -184,8 +184,4 @@ pmHDU *hdu = pmHDUFromReadout(in); // HDU of interest - - if (!pmOverscanSubtract (in, overscanOpts)) { - return false; - } // Bias frame subtraction @@ -247,5 +243,5 @@ 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, ""); Index: trunk/psModules/src/detrend/pmBias.h =================================================================== --- trunk/psModules/src/detrend/pmBias.h (revision 29831) +++ trunk/psModules/src/detrend/pmBias.h (revision 29833) @@ -21,5 +21,4 @@ /// 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 Index: trunk/psModules/src/detrend/pmDetrendDB.c =================================================================== --- trunk/psModules/src/detrend/pmDetrendDB.c (revision 29831) +++ trunk/psModules/src/detrend/pmDetrendDB.c (revision 29833) @@ -106,4 +106,5 @@ DETREND_STRING_CASE(ASTROM); DETREND_STRING_CASE(NOISEMAP); + DETREND_STRING_CASE(LINEARITY); default: return NULL; Index: trunk/psModules/src/detrend/pmDetrendDB.h =================================================================== --- trunk/psModules/src/detrend/pmDetrendDB.h (revision 29831) +++ trunk/psModules/src/detrend/pmDetrendDB.h (revision 29833) @@ -36,4 +36,5 @@ PM_DETREND_TYPE_ASTROM, PM_DETREND_TYPE_NOISEMAP, + PM_DETREND_TYPE_LINEARITY, } pmDetrendType; Index: trunk/psModules/src/detrend/pmNonLinear.c =================================================================== --- trunk/psModules/src/detrend/pmNonLinear.c (revision 29831) +++ trunk/psModules/src/detrend/pmNonLinear.c (revision 29833) @@ -10,4 +10,5 @@ #include "pmNonLinear.h" +psF32 pmNonLinearityMeasureNoisy(psF32 flux, psVector *correction_fluxes, psVector *correction_factors); pmReadout *pmNonLinearityPolynomial(pmReadout *inputReadout, const psPolynomial1D *input1DPoly) { @@ -27,21 +28,45 @@ // 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,2 +117,291 @@ 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); +} + + + Index: trunk/psModules/src/detrend/pmNonLinear.h =================================================================== --- trunk/psModules/src/detrend/pmNonLinear.h (revision 29831) +++ trunk/psModules/src/detrend/pmNonLinear.h (revision 29833) @@ -31,5 +31,6 @@ 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