Index: /branches/eam_branches/ipp-20101103/psModules/src/camera/pmFPAfile.c =================================================================== --- /branches/eam_branches/ipp-20101103/psModules/src/camera/pmFPAfile.c (revision 29851) +++ /branches/eam_branches/ipp-20101103/psModules/src/camera/pmFPAfile.c (revision 29852) @@ -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: /branches/eam_branches/ipp-20101103/psModules/src/camera/pmFPAfile.h =================================================================== --- /branches/eam_branches/ipp-20101103/psModules/src/camera/pmFPAfile.h (revision 29851) +++ /branches/eam_branches/ipp-20101103/psModules/src/camera/pmFPAfile.h (revision 29852) @@ -50,4 +50,5 @@ PM_FPA_FILE_SRCTEXT, PM_FPA_FILE_PATTERN, + PM_FPA_FILE_LINEARITY, } pmFPAfileType; Index: /branches/eam_branches/ipp-20101103/psModules/src/camera/pmFPAfileFitsIO.c =================================================================== --- /branches/eam_branches/ipp-20101103/psModules/src/camera/pmFPAfileFitsIO.c (revision 29851) +++ /branches/eam_branches/ipp-20101103/psModules/src/camera/pmFPAfileFitsIO.c (revision 29852) @@ -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: /branches/eam_branches/ipp-20101103/psModules/src/camera/pmFPAfileIO.c =================================================================== --- /branches/eam_branches/ipp-20101103/psModules/src/camera/pmFPAfileIO.c (revision 29851) +++ /branches/eam_branches/ipp-20101103/psModules/src/camera/pmFPAfileIO.c (revision 29852) @@ -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: /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmBias.c =================================================================== --- /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmBias.c (revision 29851) +++ /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmBias.c (revision 29852) @@ -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: /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmBias.h =================================================================== --- /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmBias.h (revision 29851) +++ /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmBias.h (revision 29852) @@ -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: /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmDetrendDB.c =================================================================== --- /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmDetrendDB.c (revision 29851) +++ /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmDetrendDB.c (revision 29852) @@ -106,4 +106,5 @@ DETREND_STRING_CASE(ASTROM); DETREND_STRING_CASE(NOISEMAP); + DETREND_STRING_CASE(LINEARITY); default: return NULL; Index: /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmDetrendDB.h =================================================================== --- /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmDetrendDB.h (revision 29851) +++ /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmDetrendDB.h (revision 29852) @@ -36,4 +36,5 @@ PM_DETREND_TYPE_ASTROM, PM_DETREND_TYPE_NOISEMAP, + PM_DETREND_TYPE_LINEARITY, } pmDetrendType; Index: /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmNonLinear.c =================================================================== --- /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmNonLinear.c (revision 29851) +++ /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmNonLinear.c (revision 29852) @@ -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,267 @@ 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); + } + + 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; + } + } + + if (!isfinite(result)) { + result = 0.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); + } + + 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; + } + } + + if (!isfinite(result)) { + result = 0.0; + } + return(result); +} + + + Index: /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmNonLinear.h =================================================================== --- /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmNonLinear.h (revision 29851) +++ /branches/eam_branches/ipp-20101103/psModules/src/detrend/pmNonLinear.h (revision 29852) @@ -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 Index: /branches/eam_branches/ipp-20101103/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/eam_branches/ipp-20101103/psModules/src/imcombine/pmSubtraction.c (revision 29851) +++ /branches/eam_branches/ipp-20101103/psModules/src/imcombine/pmSubtraction.c (revision 29852) @@ -1391,4 +1391,8 @@ out1->covariance = psImageCovarianceAverage(covars); psFree(covars); + if (!out1->covariance) { + psError(PM_ERR_UNKNOWN, false, "psImageCovarianceAverage returned NULL for out1."); + return false; + } // Remove covariance factor from covariance, since we've put it in the variance map already float factor = psImageCovarianceFactor(out1->covariance); @@ -1408,4 +1412,8 @@ out2->covariance = psImageCovarianceAverage(covars); psFree(covars); + if (!out2->covariance) { + psError(PM_ERR_UNKNOWN, false, "psImageCovarianceAverage returned NULL for out2."); + return false; + } // Remove covariance factor from covariance, since we've put it in the variance map already float factor = psImageCovarianceFactor(out2->covariance);