Index: branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfile.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfile.c (revision 41855) +++ branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfile.c (revision 41856) @@ -551,4 +551,7 @@ if (!strcasecmp(type, "KH.CORRECT")) { return PM_FPA_FILE_KH_CORRECT; + } + if (!strcasecmp(type, "PATTERN.ROW.AMP")) { + return PM_FPA_FILE_PATTERN_ROW_AMP; } if (!strcasecmp(type, "SUBKERNEL")) { @@ -606,4 +609,6 @@ case PM_FPA_FILE_KH_CORRECT: return ("KH.CORRECT"); + case PM_FPA_FILE_PATTERN_ROW_AMP: + return ("PATTERN.ROW.AMP"); case PM_FPA_FILE_SUBKERNEL: return ("SUBKERNEL"); Index: branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfile.h =================================================================== --- branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfile.h (revision 41855) +++ branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfile.h (revision 41856) @@ -52,4 +52,5 @@ PM_FPA_FILE_SRCTEXT, PM_FPA_FILE_PATTERN, + PM_FPA_FILE_PATTERN_ROW_AMP, PM_FPA_FILE_LINEARITY, PM_FPA_FILE_EXPNUM, Index: branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfileFitsIO.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfileFitsIO.c (revision 41855) +++ branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfileFitsIO.c (revision 41856) @@ -152,4 +152,5 @@ case PM_FPA_FILE_ASTROM_REFSTARS: case PM_FPA_FILE_KH_CORRECT: + case PM_FPA_FILE_PATTERN_ROW_AMP: { pmHDU *hdu = pmFPAviewThisHDU(view, fpa); Index: branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfileIO.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfileIO.c (revision 41855) +++ branches/eam_branches/ipp-dev-20210817/psModules/src/camera/pmFPAfileIO.c (revision 41856) @@ -253,4 +253,7 @@ status = pmPatternRead(view, file, config); break; + case PM_FPA_FILE_PATTERN_ROW_AMP: + status = pmPatternRowAmpRead(view, file, config); + break; case PM_FPA_FILE_SX: case PM_FPA_FILE_RAW: @@ -363,4 +366,5 @@ case PM_FPA_FILE_ASTROM_REFSTARS: case PM_FPA_FILE_KH_CORRECT: + case PM_FPA_FILE_PATTERN_ROW_AMP: case PM_FPA_FILE_JPEG: case PM_FPA_FILE_KAPA: @@ -451,4 +455,8 @@ return true; } + if (file->type == PM_FPA_FILE_PATTERN_ROW_AMP) { + psTrace("psModules.camera", 6, "skip write for %s, no write function defined", file->name); + return true; + } // open the file if not yet opened @@ -547,4 +555,8 @@ case PM_FPA_FILE_KH_CORRECT: psError(PS_ERR_IO, true, "cannot write type KH.CORRECT (%s)", file->name); + break; + + case PM_FPA_FILE_PATTERN_ROW_AMP: + psError(PS_ERR_IO, true, "cannot write type PATTERN.ROW.AMP (%s)", file->name); break; @@ -624,4 +636,5 @@ case PM_FPA_FILE_ASTROM_REFSTARS: case PM_FPA_FILE_KH_CORRECT: + case PM_FPA_FILE_PATTERN_ROW_AMP: case PM_FPA_FILE_LINEARITY: case PM_FPA_FILE_EXPNUM: @@ -702,4 +715,5 @@ case PM_FPA_FILE_ASTROM_REFSTARS: case PM_FPA_FILE_KH_CORRECT: + case PM_FPA_FILE_PATTERN_ROW_AMP: case PM_FPA_FILE_EXPNUM: psTrace ("psModules.camera", 6, "NOT freeing %s (%s) : save for further analysis\n", file->filename, file->name); @@ -865,4 +879,5 @@ case PM_FPA_FILE_ASTROM_REFSTARS: case PM_FPA_FILE_KH_CORRECT: + case PM_FPA_FILE_PATTERN_ROW_AMP: case PM_FPA_FILE_LINEARITY: case PM_FPA_FILE_EXPNUM: @@ -1069,4 +1084,5 @@ case PM_FPA_FILE_ASTROM_MODEL: case PM_FPA_FILE_KH_CORRECT: + case PM_FPA_FILE_PATTERN_ROW_AMP: case PM_FPA_FILE_SX: case PM_FPA_FILE_RAW: Index: branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmDetrendDB.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmDetrendDB.c (revision 41855) +++ branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmDetrendDB.c (revision 41856) @@ -112,4 +112,5 @@ DETREND_STRING_CASE(AUXMASK); DETREND_STRING_CASE(KH_CORRECT); + DETREND_STRING_CASE(PATTERN_ROW_AMP); default: return NULL; Index: branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmDetrendDB.h =================================================================== --- branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmDetrendDB.h (revision 41855) +++ branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmDetrendDB.h (revision 41856) @@ -41,4 +41,5 @@ PM_DETREND_TYPE_AUXMASK, PM_DETREND_TYPE_KH_CORRECT, + PM_DETREND_TYPE_PATTERN_ROW_AMP, } pmDetrendType; Index: branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPattern.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPattern.c (revision 41855) +++ branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPattern.c (revision 41856) @@ -160,15 +160,20 @@ +// XXX allow user choice of binned vs unbinned analysis? bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, float thresh, psStatsOptions clipMean, psStatsOptions clipStdev, psImageMaskType maskVal, psImageMaskType maskBad) { - bool status = pmPatternRowBinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad); - // bool status = pmPatternRowUnbinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad); - return status; + bool status = false; + if (true) { + status = pmPatternRowBinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad); + } else { + status = pmPatternRowUnbinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad); + } + return status; } -// USE_BACKGROUND_STDEV: if true, the analysis will use the measured robust stdev to clip the out-of-range pixles -// if false, the stdev will be estimated based on the Poisson statistics of the median (sky level). +// USE_BACKGROUND_STDEV: if TRUE, the analysis will use the measured robust stdev to clip the out-of-range pixles +// if FALSE, the stdev will be estimated based on the Poisson statistics of the median (sky level). # define USE_BACKGROUND_STDEV 0 @@ -184,4 +189,11 @@ PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false); + bool mdok; // Status of MD lookup + + pmCell *cell = ro->parent; + pmChip *chip = cell->parent; + const char *chipName = psMetadataLookupStr(&mdok, chip->concepts, "CHIP.NAME"); // Name of chip + const char *cellName = psMetadataLookupStr(&mdok, cell->concepts, "CELL.NAME"); // Name of cell + psImage *image = ro->image; // Image to correct psImage *mask = ro->mask; // Mask for image @@ -191,15 +203,8 @@ psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskVal, rng)) { - psWarning("Unable to calculate statistics on readout; masking entire readout."); + psWarning("Unable to calculate statistics on readout; skipping pattern correction for %s, %s.", chipName, cellName); psErrorClear(); psFree(stats); psFree(rng); - psImageInit(image, NAN); - if (mask) { - psBinaryOp(mask, mask, "|", psScalarAlloc(maskBad, PS_TYPE_IMAGE_MASK)); - } - if (ro->variance) { - psImageInit(image, NAN); - } return true; } @@ -226,19 +231,29 @@ // so smaller by a factor of 1.4 than what we predict, which is not very large - # define READNOISE 10 + // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file) + float nominalAmplitude = psMetadataLookupF32 (&mdok, cell->analysis, "PTN.ROW.AMP"); + if (!mdok) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number + // If we cannot determine the nominal amplitude, we fall-back on the worst case + + // XXX retrieve noise and gain from 'concepts' if possible +# define READNOISE 10 /* arbitrary number */ float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE)); - float lower = stats->robustMedian - thresh * sigma; // Lower bound for data - float upper = stats->robustMedian + thresh * sigma; // Upper bound for data + float delta = PS_MIN (thresh * sigma, 2*nominalAmplitude); + float lower = stats->robustMedian - delta; // Lower bound for data + float upper = stats->robustMedian + delta; // Upper bound for data float background = stats->robustMedian; + + // if the noise from the background is too large + float significance = nominalAmplitude / sigma; + + // XXX EAM : arbitrary number + if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) { + psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); + return true; + } + fprintf (stderr, "correcting pattern row background %s: %f - %f - %f : %f %f %f\n", cellName, lower, background, upper, sigma, nominalAmplitude, significance); + psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); # endif - pmCell *cell = ro->parent; - const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell - fprintf (stderr, "pattern row background %s: %f - %f - %f\n", cellName,lower, background, upper); - - - // XXX add code to skip fit if background stdev is large compared to expected signal - // XXX add a constraint to the fit so the |amplitude| < x (x ~ 20 - 25) - psFree(stats); psFree(rng); @@ -458,5 +473,5 @@ } -# define NPIX 31 +# define NPIX 15 // bin by NPIX in the x-direction to reduce the number of calculations needed to measure @@ -473,15 +488,19 @@ PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false); + bool mdok; // Status of MD lookup + + pmCell *cell = ro->parent; + pmChip *chip = cell->parent; + const char *chipName = psMetadataLookupStr(&mdok, chip->concepts, "CHIP.NAME"); // Name of chip + const char *cellName = psMetadataLookupStr(&mdok, cell->concepts, "CELL.NAME"); // Name of cell + psImage *image = ro->image; // Image to correct psImage *mask = ro->mask; // Mask for image int numCols = image->numCols, numRows = image->numRows; // Size of image - // use the background and typical amplitude value to choose to correct or not - // pmCell *cell = ro->parent; - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskVal, rng)) { - psWarning("Unable to calculate statistics on readout; skipping pattern correction."); + psWarning("Unable to calculate statistics on readout; skipping pattern correction for %s, %s.", chipName, cellName); psErrorClear(); psFree(stats); @@ -525,22 +544,29 @@ // so smaller by a factor of 1.4 than what we predict, which is not very large + // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file) + float nominalAmplitude = psMetadataLookupF32 (&mdok, cell->analysis, "PTN.ROW.AMP"); + if (!mdok) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number + if (!isfinite(nominalAmplitude)) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number + // If we cannot determine the nominal amplitude, we fall-back on the worst case + // retrieve noise and gain from 'concepts' if possible # define READNOISE 10 /* arbitrary number */ float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE)); - float delta = PS_MIN (thresh * sigma, 40); // XXX EAM : arbitrary number + float delta = PS_MIN (thresh * sigma, 5*nominalAmplitude); float lower = stats->robustMedian - delta; // Lower bound for data float upper = stats->robustMedian + delta; // Upper bound for data float background = stats->robustMedian; - // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file) - float nominalAmplitude = 20; // XXX EAM : arbitrary number - // if the noise from the background is too large float significance = nominalAmplitude / sigma; - if (significance < 0.5) return true; // XXX EAM : arbitrary number - - pmCell *cell = ro->parent; - const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell - fprintf (stderr, "pattern row background %s: %f - %f - %f\n", cellName, lower, background, upper); + + // XXX EAM : arbitrary number + if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) { + psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); + psFree(stats); + psFree(rng); + return true; + } + psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); # endif @@ -643,9 +669,6 @@ } - // XXX how much time is spent in the fitting? to test: make this always true + // If not enough points are valid, skip if (num < validNmin) { - // Not enough points to fit - // EAM 20211011 : we used to mask rows which could not be fit, but this is excessive. just skip. - // patternMaskRow(ro, y, maskBad); // Ignore this row in our subsequent fits, because the fit failed. yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; @@ -656,6 +679,4 @@ psWarning("Unable to fit polynomial to row %d", y); psErrorClear(); - // EAM 20211011 : we used to mask rows which could not be fit, but this is excessive. just skip. - // patternMaskRow(ro, y, maskBad); // Ignore this row in our subsequent fits, because the fit failed. yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; @@ -672,6 +693,4 @@ psWarning("Unable to evaluate polynomial for row %d", y); psErrorClear(); - // EAM 20211011 : we used to mask rows which could not be fit, but this is excessive. just skip. - // patternMaskRow(ro, y, maskBad); yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; continue; Index: branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPatternIO.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPatternIO.c (revision 41855) +++ branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPatternIO.c (revision 41856) @@ -8,4 +8,5 @@ #include "pmFPAview.h" #include "pmFPAfile.h" +#include "pmFPAUtils.h" #include "pmFPAfileFitsIO.h" #include "pmFPAHeader.h" @@ -455,8 +456,5 @@ /**************** PatternRowAmp(litude) I/O *************************/ -# if (0) -/********************* Read Data functions *****************************/ - -bool pmPatternRowAmpReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) +bool pmPatternRowAmpRead (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) { // read the full model in one pass: require the level to be FPA @@ -499,4 +497,6 @@ } + int nGood = 0; + int nTotal = 0; while (haveData) { @@ -511,9 +511,10 @@ // determine the chip: char *extname = psMetadataLookupStr (&status, header, "EXTNAME"); - psLogMsg ("psModules.detrend", 4, "read %ld rows from Pattern Row Amplitude file, extname %s\n", model->n, extname); + psLogMsg ("psModules.detrend", 8, "read %ld rows from Pattern Row Amplitude file, extname %s\n", model->n, extname); + nTotal += model->n; // I expect to find a name of the form: chipName.ptn (eg, XY01.ptn) // where chipName like 'XY01' - psAssert (strlen(extname) == 10, "invalid extension %s", extname); + psAssert (strlen(extname) == 8, "invalid extension %s", extname); psAssert (extname[5] == 'p', "invalid extension %s", extname); psAssert (extname[6] == 't', "invalid extension %s", extname); @@ -530,12 +531,18 @@ for (int i = 0; i < model->n; i++) { psMetadata *row = model->data[i]; - char *cellName = psMetadataLookupStr(&status, row, "cellname"); - float *amplitude = psMetadataLookupF32(&status, row, "medianValue"); + + char *cellName = psMetadataLookupStr(&status, row, "CELL_NAME"); + if (!cellName) continue; + + float amplitude = psMetadataLookupF32(&status, row, "VALUE_MEDIAN"); int cellNumber = pmChipFindCell (chip, cellName); - pmCell *cell = chip->data[cellNumber]; + pmCell *cell = chip->cells->data[cellNumber]; + if (!cell) continue; psMetadataAddF32 (cell->analysis, PS_LIST_TAIL, "PTN.ROW.AMP", PS_META_REPLACE, "", amplitude); + nGood ++; } + psFree (model); psFree (header); @@ -544,6 +551,6 @@ haveData = psFitsMoveExtNum (file->fits, 1, true); } - - return true; -} -# endif + psLogMsg ("psModules.detrend", 4, "read %d of %d rows from Pattern Row Amplitude file\n", nGood, nTotal); + + return true; +} Index: branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPatternIO.h =================================================================== --- branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPatternIO.h (revision 41855) +++ branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPatternIO.h (revision 41856) @@ -35,3 +35,8 @@ ); + +bool pmPatternRowAmpRead (const pmFPAview *view, pmFPAfile *file, const pmConfig *config); +bool pmPatternRowAmpReadFPA (pmFPAfile *file); +bool pmPatternRowAmpReadChips (pmFPAfile *file); + #endif