Index: /branches/2dbias/ippconfig/gpc1/format_20100723.config =================================================================== --- /branches/2dbias/ippconfig/gpc1/format_20100723.config (revision 42678) +++ /branches/2dbias/ippconfig/gpc1/format_20100723.config (revision 42679) @@ -166,5 +166,5 @@ # 2D bias subtraction requires ppImage.config OVERSCAN.2D to be TRUE CELL.BIASSEC.SOURCE STR VALUE - CELL.BIASSEC STR [591:624,1:598],[1:590,598:608] + CELL.BIASSEC STR [591:624,1:598],[1:624,599:608] # # a single region will be used for constant and 1D bias subtraction @@ -174,4 +174,70 @@ END +OVERSCAN.2D.SUBSET METADATA + XY00 BOOL TRUE + XY01 BOOL TRUE + XY02 BOOL TRUE + XY03 BOOL TRUE + XY04 BOOL TRUE + XY05 BOOL TRUE + XY06 BOOL TRUE + XY07 BOOL TRUE + XY10 BOOL TRUE + XY11 BOOL TRUE + XY12 BOOL TRUE + XY13 BOOL TRUE + XY14 BOOL TRUE + XY15 BOOL TRUE + XY16 BOOL TRUE + XY17 BOOL TRUE + XY20 BOOL TRUE + XY21 BOOL TRUE + XY22 BOOL TRUE + XY23 BOOL TRUE + XY24 BOOL TRUE + XY25 BOOL TRUE + XY26 BOOL TRUE + XY27 BOOL TRUE + XY30 STR 1111100011111000111100001111000011110000111110001111000011110000 + XY31 BOOL TRUE + XY32 BOOL TRUE + XY33 BOOL TRUE + XY34 BOOL TRUE + XY35 BOOL TRUE + XY36 BOOL TRUE + XY37 BOOL TRUE + XY40 BOOL TRUE + XY41 BOOL TRUE + XY42 BOOL TRUE + XY43 BOOL TRUE + XY44 BOOL TRUE + XY45 BOOL TRUE + XY46 BOOL TRUE + XY47 BOOL TRUE + XY50 BOOL TRUE + XY51 BOOL TRUE + XY52 BOOL TRUE + XY53 BOOL TRUE + XY54 BOOL TRUE + XY55 BOOL TRUE + XY56 BOOL TRUE + XY57 BOOL TRUE + XY60 BOOL TRUE + XY61 BOOL TRUE + XY62 BOOL TRUE + XY63 BOOL TRUE + XY64 BOOL TRUE + XY65 BOOL TRUE + XY66 BOOL TRUE + XY67 BOOL TRUE + XY70 BOOL TRUE + XY71 BOOL TRUE + XY72 BOOL TRUE + XY73 BOOL TRUE + XY74 BOOL TRUE + XY75 BOOL TRUE + XY76 BOOL TRUE + XY77 BOOL TRUE +END # How to translate PS concepts into FITS headers Index: /branches/2dbias/ippconfig/gpc1/ppImage.config =================================================================== --- /branches/2dbias/ippconfig/gpc1/ppImage.config (revision 42678) +++ /branches/2dbias/ippconfig/gpc1/ppImage.config (revision 42679) @@ -9,4 +9,71 @@ OVERSCAN.STAT STR MEDIAN OVERSCAN.BOXCAR S32 3 + +OVERSCAN.2D.SUBSET METADATA + XY00 BOOL TRUE + XY01 BOOL TRUE + XY02 BOOL TRUE + XY03 BOOL TRUE + XY04 BOOL TRUE + XY05 BOOL TRUE + XY06 BOOL TRUE + XY07 BOOL TRUE + XY10 BOOL TRUE + XY11 BOOL TRUE + XY12 BOOL TRUE + XY13 BOOL TRUE + XY14 BOOL TRUE + XY15 BOOL TRUE + XY16 BOOL TRUE + XY17 BOOL TRUE + XY20 BOOL TRUE + XY21 BOOL TRUE + XY22 BOOL TRUE + XY23 BOOL TRUE + XY24 BOOL TRUE + XY25 BOOL TRUE + XY26 BOOL TRUE + XY27 BOOL TRUE + XY30 BOOL TRUE + XY31 BOOL TRUE + XY32 BOOL TRUE + XY33 BOOL TRUE + XY34 BOOL TRUE + XY35 BOOL TRUE + XY36 BOOL TRUE + XY37 BOOL TRUE + XY40 BOOL TRUE + XY41 BOOL TRUE + XY42 BOOL TRUE + XY43 BOOL TRUE + XY44 BOOL TRUE + XY45 BOOL TRUE + XY46 BOOL TRUE + XY47 BOOL TRUE + XY50 BOOL TRUE + XY51 BOOL TRUE + XY52 BOOL TRUE + XY53 BOOL TRUE + XY54 BOOL TRUE + XY55 BOOL TRUE + XY56 BOOL TRUE + XY57 BOOL TRUE + XY60 BOOL TRUE + XY61 BOOL TRUE + XY62 BOOL TRUE + XY63 BOOL TRUE + XY64 BOOL TRUE + XY65 BOOL TRUE + XY66 BOOL TRUE + XY67 BOOL TRUE + XY70 BOOL TRUE + XY71 BOOL TRUE + XY72 BOOL TRUE + XY73 BOOL TRUE + XY74 BOOL TRUE + XY75 BOOL TRUE + XY76 BOOL TRUE + XY77 BOOL TRUE +END # Normalization class Index: /branches/2dbias/ippconfig/recipes/ppImage.config =================================================================== --- /branches/2dbias/ippconfig/recipes/ppImage.config (revision 42678) +++ /branches/2dbias/ippconfig/recipes/ppImage.config (revision 42679) @@ -110,4 +110,6 @@ OVERSCAN.2D BOOL FALSE # use a 2D model for the overscan subtraction? +OVERSCAN.2D.SUBSET METADATA +END # statistics to use to convert 2D top-overscan region to a vector Index: /branches/2dbias/ppImage/src/ppImage.h =================================================================== --- /branches/2dbias/ppImage/src/ppImage.h (revision 42678) +++ /branches/2dbias/ppImage/src/ppImage.h (revision 42679) @@ -188,4 +188,5 @@ bool ppImageDetrendFree(pmConfig *config, pmFPAview *view); +bool ppImageDoPatternForView(bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipeName, const char *recipeValue); bool ppImageFringeFree(pmConfig *config, pmFPAview *view); Index: /branches/2dbias/ppImage/src/ppImageDetrendPattern.c =================================================================== --- /branches/2dbias/ppImage/src/ppImageDetrendPattern.c (revision 42678) +++ /branches/2dbias/ppImage/src/ppImageDetrendPattern.c (revision 42679) @@ -5,11 +5,12 @@ #include "ppImage.h" -#define ESCAPE(STATUS,...) { \ - psError(PS_ERR_UNKNOWN, STATUS, __VA_ARGS__); \ - psFree(view); \ - return false; \ - } - -static bool doPatternForView (bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipename, const char *recipevalue); +#define ESCAPE(STATUS, ...) \ + { \ + psError(PS_ERR_UNKNOWN, STATUS, __VA_ARGS__); \ + psFree(view); \ + return false; \ + } + +bool ppImageDoPatternForView(bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipename, const char *recipevalue); bool ppImageDetrendPatternRowApply(pmConfig *config, pmChip *chip, const pmFPAview *inputView, ppImageOptions *options); @@ -26,17 +27,21 @@ assert(inputView->readout == -1); - if (!ppImageDetrendPatternRowApply (config, chip, inputView, options)) { - return false; - } - if (!ppImageDetrendPatternContinuityApply (config, chip, inputView, options)) { - return false; - } - if (!ppImageDetrendPatternCellApply (config, chip, inputView, options)) { - return false; - } - if (!ppImageDetrendPatternDeadCellsApply (config, chip, inputView, options)) { - return false; - } - return(true); + if (!ppImageDetrendPatternRowApply(config, chip, inputView, options)) + { + return false; + } + if (!ppImageDetrendPatternContinuityApply(config, chip, inputView, options)) + { + return false; + } + if (!ppImageDetrendPatternCellApply(config, chip, inputView, options)) + { + return false; + } + if (!ppImageDetrendPatternDeadCellsApply(config, chip, inputView, options)) + { + return false; + } + return (true); } @@ -45,6 +50,7 @@ bool status; - if (!options->doPatternRow) return true; - + if (!options->doPatternRow) + return true; + assert(inputView->chip != -1); assert(inputView->cell == -1); @@ -73,7 +79,8 @@ pmHDU *hdu = pmHDUFromChip(chip); - if (psMetadataLookupBool(&status, hdu->header, "PTRN_ROW")) { - psLogMsg("ppImage", PS_LOG_INFO, "Not performing row pattern correction as it has already been done."); - return true; + if (psMetadataLookupBool(&status, hdu->header, "PTRN_ROW")) + { + psLogMsg("ppImage", PS_LOG_INFO, "Not performing row pattern correction as it has already been done."); + return true; } @@ -85,92 +92,108 @@ // grab the PATTERN.ROW.AMP file - pmFPAfile *PRAfile = psMetadataLookupPtr (&status, config->files, "PPIMAGE.PATTERN.ROW.AMP"); - if (!PRAfile) { - psLogMsg("ppImage", PS_LOG_INFO, "Pattern Row Amplitude file not found, applying to all (with limits from ROW.SUBSET) "); - } - + pmFPAfile *PRAfile = psMetadataLookupPtr(&status, config->files, "PPIMAGE.PATTERN.ROW.AMP"); + if (!PRAfile) + { + psLogMsg("ppImage", PS_LOG_INFO, "Pattern Row Amplitude file not found, applying to all (with limits from ROW.SUBSET) "); + } + pmFPAview *view = pmFPAviewAlloc(0); // View for local processing *view = *inputView; pmCell *cell = NULL; - while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) { - // const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); - // const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); - // psWarning ("Looping through %s, %s\n", chipName, cellName); - - if (!cell->process || !cell->file_exists) { - continue; - } - - // this forces pmFPAfileRead of the PATTERN.ROW.AMP file (XXX but is this needed?) - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - ESCAPE(false, "load failure for Cell"); - } - - if (!cell->data_exists) { - continue; - } - - if (cell->readouts->n > 1) { - psWarning ("Skipping Video Cell for ppImageDetrendPatternApply"); - continue; - } - - // grab the corresponding cell - if (PRAfile) { - pmCell *PRAcell = pmFPAviewThisCell (view, PRAfile->fpa); - psAssert (PRAcell, "found Pattern Row Amplitude file, but not cell?"); - - // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file) - float amplitude = psMetadataLookupF32 (&status, PRAcell->analysis, "PTN.ROW.AMP"); - if (!status) amplitude = NAN; - - // put the value on the science cell - psMetadataAddF32 (cell->analysis, PS_LIST_TAIL, "PTN.ROW.AMP", PS_META_REPLACE, "", amplitude); - } - - bool doPattern = false; - if (!doPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.ROW.SUBSET")) { - ESCAPE(false, "Unable to determine whether row pattern matching should be applied."); - } - if (!doPattern) continue; - - // switch to test threaded version - if (true) { - // I need to allocate a view here to be freed by the - // called function below. - pmFPAview *myView = pmFPAviewAlloc(0); // View for local processing - *myView = *view; - - // allocate a job, construct the arguments for this job - psThreadJob *job = psThreadJobAlloc("PPIMAGE_PATTERN_ROW_CELL"); - psArrayAdd(job->args, 1, config); - psArrayAdd(job->args, 1, input->fpa); - psArrayAdd(job->args, 1, chip); - psArrayAdd(job->args, 1, cell); - psArrayAdd(job->args, 1, myView); - psArrayAdd(job->args, 1, options); - if (!psThreadJobAddPending(job)) { - return false; - } - } else { - // bump the counter since it must be freed by the function below. - psMemIncrRefCounter (view); // View for local processing - if (!ppImageDetrendPatternApplyCell (config, input->fpa, chip, cell, view, options)) { - return false; - } - } + while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) + { + // const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); + // const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); + // psWarning ("Looping through %s, %s\n", chipName, cellName); + + if (!cell->process || !cell->file_exists) + { + continue; + } + + // this forces pmFPAfileRead of the PATTERN.ROW.AMP file (XXX but is this needed?) + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) + { + ESCAPE(false, "load failure for Cell"); + } + + if (!cell->data_exists) + { + continue; + } + + if (cell->readouts->n > 1) + { + psWarning("Skipping Video Cell for ppImageDetrendPatternApply"); + continue; + } + + // grab the corresponding cell + if (PRAfile) + { + pmCell *PRAcell = pmFPAviewThisCell(view, PRAfile->fpa); + psAssert(PRAcell, "found Pattern Row Amplitude file, but not cell?"); + + // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file) + float amplitude = psMetadataLookupF32(&status, PRAcell->analysis, "PTN.ROW.AMP"); + if (!status) + amplitude = NAN; + + // put the value on the science cell + psMetadataAddF32(cell->analysis, PS_LIST_TAIL, "PTN.ROW.AMP", PS_META_REPLACE, "", amplitude); + } + + bool doPattern = false; + if (!ppImageDoPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.ROW.SUBSET")) + { + ESCAPE(false, "Unable to determine whether row pattern matching should be applied."); + } + if (!doPattern) + continue; + + // switch to test threaded version + if (true) + { + // I need to allocate a view here to be freed by the + // called function below. + pmFPAview *myView = pmFPAviewAlloc(0); // View for local processing + *myView = *view; + + // allocate a job, construct the arguments for this job + psThreadJob *job = psThreadJobAlloc("PPIMAGE_PATTERN_ROW_CELL"); + psArrayAdd(job->args, 1, config); + psArrayAdd(job->args, 1, input->fpa); + psArrayAdd(job->args, 1, chip); + psArrayAdd(job->args, 1, cell); + psArrayAdd(job->args, 1, myView); + psArrayAdd(job->args, 1, options); + if (!psThreadJobAddPending(job)) + { + return false; + } + } + else + { + // bump the counter since it must be freed by the function below. + psMemIncrRefCounter(view); // View for local processing + if (!ppImageDetrendPatternApplyCell(config, input->fpa, chip, cell, view, options)) + { + return false; + } + } } // wait here for the threaded jobs to finish - // if no threads are allocated, this - if (!psThreadPoolWait(true, true)) { - psError(PS_ERR_UNKNOWN, false, "Unable to apply bias correction."); - return false; - } - - psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_ROW",PS_META_REPLACE,"PATTERN.ROW correction applied",true); + // if no threads are allocated, this + if (!psThreadPoolWait(true, true)) + { + psError(PS_ERR_UNKNOWN, false, "Unable to apply bias correction."); + return false; + } + + psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_ROW", PS_META_REPLACE, "PATTERN.ROW correction applied", true); psFree(view); - + return true; } @@ -180,5 +203,6 @@ bool status; - if (!options->doPatternContinuity) return true; + if (!options->doPatternContinuity) + return true; assert(inputView->chip != -1); @@ -188,13 +212,14 @@ // see the comment for PATTERN.ROW; the same rules apply for PATTERN.CELL - int numCells = chip->cells->n; // Number of cells + int numCells = chip->cells->n; // Number of cells psVector *tweak = psVectorAlloc(numCells, PS_TYPE_U8); // Tweak cell? - pmFPAview *view = pmFPAviewAlloc(0); // View for local processing + pmFPAview *view = pmFPAviewAlloc(0); // View for local processing *view = *inputView; pmHDU *hdu = pmHDUFromChip(chip); - if (psMetadataLookupBool(&status, hdu->header, "PTRN_CON")) { - psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell continuity correction as it has already been done."); - return true; + if (psMetadataLookupBool(&status, hdu->header, "PTRN_CON")) + { + psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell continuity correction as it has already been done."); + return true; } @@ -202,35 +227,40 @@ psLogMsg("ppImage", PS_LOG_INFO, "Performing pattern continuity correction for %s\n", chipName); - for (int i = 0; i < chip->cells->n; i++) { - view->cell = i; - - pmCell *cell = chip->cells->data[i]; // Cell of interest - - if (cell->readouts->n > 1) { - psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell continuity correction on video cell."); - continue; - } - - bool doPattern = false; - if (!doPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.CONTINUITY.SUBSET")) { - ESCAPE(false, "Unable to determine whether row pattern matching should be applied."); - } - - if (doPattern) { - tweak->data.U8[i] = 0xFF; - } + for (int i = 0; i < chip->cells->n; i++) + { + view->cell = i; + + pmCell *cell = chip->cells->data[i]; // Cell of interest + + if (cell->readouts->n > 1) + { + psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell continuity correction on video cell."); + continue; + } + + bool doPattern = false; + if (!ppImageDoPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.CONTINUITY.SUBSET")) + { + ESCAPE(false, "Unable to determine whether row pattern matching should be applied."); + } + + if (doPattern) + { + tweak->data.U8[i] = 0xFF; + } } // Tweak the cells if (!pmPatternContinuity(chip, tweak, options->patternCellBG, options->patternCellMean, - options->maskValue, options->darkMask,options->patternContinuityEdgeWidth)) { - psFree(tweak); - psFree(view); - return false; + options->maskValue, options->darkMask, options->patternContinuityEdgeWidth)) + { + psFree(tweak); + psFree(view); + return false; } psFree(tweak); psFree(view); - psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_CON",PS_META_REPLACE,"PATTERN.CONTINUITY correction applied",true); + psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_CON", PS_META_REPLACE, "PATTERN.CONTINUITY correction applied", true); return true; } @@ -240,5 +270,6 @@ bool status; - if (!options->doPatternCell) return true; + if (!options->doPatternCell) + return true; assert(inputView->chip != -1); @@ -246,13 +277,14 @@ assert(inputView->readout == -1); - int numCells = chip->cells->n; // Number of cells + int numCells = chip->cells->n; // Number of cells psVector *tweak = psVectorAlloc(numCells, PS_TYPE_U8); // Tweak cell? - pmFPAview *view = pmFPAviewAlloc(0); // View for local processing + pmFPAview *view = pmFPAviewAlloc(0); // View for local processing *view = *inputView; pmHDU *hdu = pmHDUFromChip(chip); - if (psMetadataLookupBool(&status, hdu->header, "PTRN_CEL")) { - psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell pattern correction as it has already been done."); - return true; + if (psMetadataLookupBool(&status, hdu->header, "PTRN_CEL")) + { + psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell pattern correction as it has already been done."); + return true; } @@ -260,115 +292,135 @@ psLogMsg("ppImage", PS_LOG_INFO, "Performing cell pattern correction for %s\n", chipName); - for (int i = 0; i < chip->cells->n; i++) { - view->cell = i; - - pmCell *cell = chip->cells->data[i]; // Cell of interest - - if (cell->readouts->n > 1) { - psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell pattern correction on video cell."); - continue; - } - - bool doPattern = false; - if (!doPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.CELL.SUBSET")) { - ESCAPE(false, "Unable to determine whether row pattern matching should be applied."); - } - if (doPattern) { - tweak->data.U8[i] = 0xFF; - } + for (int i = 0; i < chip->cells->n; i++) + { + view->cell = i; + + pmCell *cell = chip->cells->data[i]; // Cell of interest + + if (cell->readouts->n > 1) + { + psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell pattern correction on video cell."); + continue; + } + + bool doPattern = false; + if (!ppImageDoPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.CELL.SUBSET")) + { + ESCAPE(false, "Unable to determine whether row pattern matching should be applied."); + } + if (doPattern) + { + tweak->data.U8[i] = 0xFF; + } } // Tweak the cells if (!pmPatternCell(chip, tweak, options->patternCellBG, options->patternCellMean, - options->maskValue, options->darkMask)) { - psFree(tweak); - psFree(view); - return false; + options->maskValue, options->darkMask)) + { + psFree(tweak); + psFree(view); + return false; } psFree(tweak); psFree(view); - psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_CEL",PS_META_REPLACE,"PATTERN.CELL correction applied",true); + psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_CEL", PS_META_REPLACE, "PATTERN.CELL correction applied", true); return true; } -bool ppImageDetrendPatternDeadCellsMask (pmChip *chip, psImageMaskType maskVal) { - - int numCells = chip->cells->n; // Number of cells +bool ppImageDetrendPatternDeadCellsMask(pmChip *chip, psImageMaskType maskVal) +{ + + int numCells = chip->cells->n; // Number of cells // now mask bad cells - for (int i = 0; i < numCells; i++) { - pmCell *cell = chip->cells->data[i]; // Cell of interest + for (int i = 0; i < numCells; i++) + { + pmCell *cell = chip->cells->data[i]; // Cell of interest pmReadout *ro = cell->readouts->data[0]; // Readout of interest - psImage *mask = ro->mask; // mask of interest - int numCols = mask->numCols, numRows = mask->numRows; // Size of image - - for (int y = 0; y < numRows; y++) { - for (int x = 0; x < numCols; x++) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskVal; - } - } + psImage *mask = ro->mask; // mask of interest + int numCols = mask->numCols, numRows = mask->numRows; // Size of image + + for (int y = 0; y < numRows; y++) + { + for (int x = 0; x < numCols; x++) + { + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskVal; + } + } } return true; } -bool ppImageDetrendPatternDeadCellsXY42 (pmChip *chip, psImageMaskType maskVal) { +bool ppImageDetrendPatternDeadCellsXY42(pmChip *chip, psImageMaskType maskVal) +{ bool status; // extract the MEDIAN_CELL_BACKGROUND values to check for problems - int numCells = chip->cells->n; // Number of cells - psVector *cellBackground = psVectorAllocEmpty (numCells, PS_DATA_F32); - - for (int i = 0; i < numCells; i++) { - pmCell *cell = chip->cells->data[i];; - - if (!cell->process || !cell->file_exists || !cell->data_exists) continue; - - // select the HDU for this cell - pmHDU *cellHDU = pmHDUFromCell(cell); // HDU of interest - - psF32 BackValue = psMetadataLookupF32(&status, cellHDU->header, "BACK_VAL"); - if (!status) continue; - - psF32 BackError = psMetadataLookupF32(&status, cellHDU->header, "BACK_ERR"); - if (!status) continue; - - if (BackValue == 0.0) continue; - if (BackError == 0.0) continue; - - psVectorAppend (cellBackground, BackValue / BackError); - } - - if (cellBackground->n < 0.375*numCells) { - // Chip is bad, mask the whole thing - ppImageDetrendPatternDeadCellsMask (chip, maskVal); - psFree (cellBackground); - return true; + int numCells = chip->cells->n; // Number of cells + psVector *cellBackground = psVectorAllocEmpty(numCells, PS_DATA_F32); + + for (int i = 0; i < numCells; i++) + { + pmCell *cell = chip->cells->data[i]; + ; + + if (!cell->process || !cell->file_exists || !cell->data_exists) + continue; + + // select the HDU for this cell + pmHDU *cellHDU = pmHDUFromCell(cell); // HDU of interest + + psF32 BackValue = psMetadataLookupF32(&status, cellHDU->header, "BACK_VAL"); + if (!status) + continue; + + psF32 BackError = psMetadataLookupF32(&status, cellHDU->header, "BACK_ERR"); + if (!status) + continue; + + if (BackValue == 0.0) + continue; + if (BackError == 0.0) + continue; + + psVectorAppend(cellBackground, BackValue / BackError); + } + + if (cellBackground->n < 0.375 * numCells) + { + // Chip is bad, mask the whole thing + ppImageDetrendPatternDeadCellsMask(chip, maskVal); + psFree(cellBackground); + return true; } // Second, calculate the median - psVectorSortInPlace (cellBackground); + psVectorSortInPlace(cellBackground); int midPt = cellBackground->n / 2.0; - float median = cellBackground->n % 2 ? cellBackground->data.F32[midPt] : 0.5*(cellBackground->data.F32[midPt] + cellBackground->data.F32[midPt-1]); - psFree (cellBackground); + float median = cellBackground->n % 2 ? cellBackground->data.F32[midPt] : 0.5 * (cellBackground->data.F32[midPt] + cellBackground->data.F32[midPt - 1]); + psFree(cellBackground); // XXX hardwired number - if (median < 3.0) { - // Chip is bad, mask the whole thing - ppImageDetrendPatternDeadCellsMask (chip, maskVal); - } - psFree (cellBackground); + if (median < 3.0) + { + // Chip is bad, mask the whole thing + ppImageDetrendPatternDeadCellsMask(chip, maskVal); + } + psFree(cellBackground); return true; } // to apply the dead cell pattern, we need to transfer the pattern for this chip from the -// pmFPAfile for the pattern to the one for the image being processed. +// pmFPAfile for the pattern to the one for the image being processed. bool ppImageDetrendPatternDeadCellsApply(pmConfig *config, pmChip *chip, const pmFPAview *inputView, ppImageOptions *options) { bool status; - if (!options->doPatternDeadCells) return true; + if (!options->doPatternDeadCells) + return true; assert(inputView->chip != -1); @@ -380,20 +432,23 @@ pmHDU *hdu = pmHDUFromChip(chip); - if (psMetadataLookupBool(&status, hdu->header, "PTRN_DED")) { - psLogMsg("ppImage", PS_LOG_INFO, "Not performing dead cell pattern correction as it has already been done."); - return true; - } - - if (!strcmp (chipName, "XY40") || !strcmp (chipName, "XY42")) { - // special case : check for BACK_VAL / BACK_ERR > 3 or < 3 - psLogMsg("ppImage", PS_LOG_INFO, "Using special case for XY40 and XY42"); - ppImageDetrendPatternDeadCellsXY42 (chip, options->blankMask); - return true; + if (psMetadataLookupBool(&status, hdu->header, "PTRN_DED")) + { + psLogMsg("ppImage", PS_LOG_INFO, "Not performing dead cell pattern correction as it has already been done."); + return true; + } + + if (!strcmp(chipName, "XY40") || !strcmp(chipName, "XY42")) + { + // special case : check for BACK_VAL / BACK_ERR > 3 or < 3 + psLogMsg("ppImage", PS_LOG_INFO, "Using special case for XY40 and XY42"); + ppImageDetrendPatternDeadCellsXY42(chip, options->blankMask); + return true; } pmFPAfile *pattern = psMetadataLookupPtr(&status, config->files, "PPIMAGE.PATTERN.DEAD.CELLS"); - if (!pattern) { - psLogMsg("ppImage", PS_LOG_INFO, "Pattern Dead Cells file not found, skipping"); - return true; + if (!pattern) + { + psLogMsg("ppImage", PS_LOG_INFO, "Pattern Dead Cells file not found, skipping"); + return true; } @@ -401,44 +456,49 @@ *view = *inputView; - pmChip *patternChip = pmFPAviewThisChip (view, pattern->fpa); + pmChip *patternChip = pmFPAviewThisChip(view, pattern->fpa); // grab the pattern from the input pattern file chip - psImage *deadCellPattern = (psImage *) psMetadataLookupPtr (&status, patternChip->analysis, "PTN.DEAD.CELL"); - if (!deadCellPattern) { + psImage *deadCellPattern = (psImage *)psMetadataLookupPtr(&status, patternChip->analysis, "PTN.DEAD.CELL"); + if (!deadCellPattern) + { psLogMsg("psModules.detrend", PS_LOG_DETAIL, "No DEAD CELL pattern for chip, skipping\n"); - psFree (view); - return true; + psFree(view); + return true; } // copy the pattern pointer to the chip of the image being processed - psMetadataAddImage (chip->analysis, PS_LIST_TAIL, "PTN.DEAD.CELL", PS_META_REPLACE, "", deadCellPattern); + psMetadataAddImage(chip->analysis, PS_LIST_TAIL, "PTN.DEAD.CELL", PS_META_REPLACE, "", deadCellPattern); // extract the MEDIAN_CELL_BACKGROUND values to check for problems - int numCells = chip->cells->n; // Number of cells - psVector *cellBackground = psVectorAllocEmpty (numCells, PS_DATA_F32); + int numCells = chip->cells->n; // Number of cells + psVector *cellBackground = psVectorAllocEmpty(numCells, PS_DATA_F32); pmCell *cell = NULL; - while ((cell = pmFPAviewNextCell(view, chip->parent, 1)) != NULL) { - if (!cell->process || !cell->file_exists || !cell->data_exists) { - psVectorAppend (cellBackground, NAN); - continue; - } - - // select the HDU for this cell - pmHDU *cellHDU = pmHDUFromCell(cell); // HDU of interest - - psF32 value = psMetadataLookupF32(&status, cellHDU->header, "BACK_VAL"); - if (!status) { - psVectorAppend (cellBackground, NAN); - continue; - } - psVectorAppend (cellBackground, value); + while ((cell = pmFPAviewNextCell(view, chip->parent, 1)) != NULL) + { + if (!cell->process || !cell->file_exists || !cell->data_exists) + { + psVectorAppend(cellBackground, NAN); + continue; + } + + // select the HDU for this cell + pmHDU *cellHDU = pmHDUFromCell(cell); // HDU of interest + + psF32 value = psMetadataLookupF32(&status, cellHDU->header, "BACK_VAL"); + if (!status) + { + psVectorAppend(cellBackground, NAN); + continue; + } + psVectorAppend(cellBackground, value); } // match cellBackground pattern to registered patterns and mask as needed - if (!pmPatternDeadCells(chip, cellBackground, options->blankMask)) { - psFree(cellBackground); - psFree(view); - return false; + if (!pmPatternDeadCells(chip, cellBackground, options->blankMask)) + { + psFree(cellBackground); + psFree(view); + return false; } psFree(cellBackground); @@ -449,5 +509,6 @@ } -static bool doPatternForView (bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipeName, const char *recipeValue) { +bool ppImageDoPatternForView(bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipeName, const char *recipeValue) +{ *doit = false; @@ -456,36 +517,41 @@ doPattern = pmConfigRecipeValueByView(config, recipeName, recipeValue, chip->parent, view); - - if (!doPattern) { + + if (!doPattern) + { psError(PS_ERR_UNKNOWN, false, "Unable to determine whether row pattern matching should be applied."); return false; } - if (doPattern->type == PS_DATA_BOOL) { + if (doPattern->type == PS_DATA_BOOL) + { *doit = doPattern->data.B; return true; } - if (doPattern->type == PS_DATA_STRING) { + if (doPattern->type == PS_DATA_STRING) + { // expect a string of the form "000110001001" with at least view->cell entries char *string = doPattern->data.str; - if (strlen(string) < view->cell) { - psError(PS_ERR_UNKNOWN, true, "error in PATTERN.ROW.SUBSET chip string (too few elements %d)", (int) strlen(string)); + if (strlen(string) < view->cell) + { + psError(PS_ERR_UNKNOWN, true, "error in PATTERN.ROW.SUBSET chip string (too few elements %d)", (int)strlen(string)); return false; } - switch (string[view->cell]) { - case '0': - case 'f': - case 'F': - case 'n': - case 'N': + switch (string[view->cell]) + { + case '0': + case 'f': + case 'F': + case 'n': + case 'N': *doit = false; return true; - case '1': - case 't': - case 'T': - case 'y': - case 'Y': + case '1': + case 't': + case 'T': + case 'y': + case 'Y': *doit = true; return true; - default: + default: psError(PS_ERR_UNKNOWN, true, "error in PATTERN.ROW.SUBSET chip string %s (unknown value %c))", string, string[view->cell]); return false; @@ -497,47 +563,54 @@ } -// thread safety : -bool ppImageDetrendPatternApplyCell (pmConfig *config, pmFPA *fpa, pmChip *chip, pmCell *cell, pmFPAview *view, ppImageOptions *options) { - +// thread safety : +bool ppImageDetrendPatternApplyCell(pmConfig *config, pmFPA *fpa, pmChip *chip, pmCell *cell, pmFPAview *view, ppImageOptions *options) +{ + // process each of the readouts - pmReadout *readout; // Readout from cell - while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { - if (!readout->data_exists) { - continue; - } - - // Perform pattern correction - if (!pmPatternRow(readout, options->patternRowOrder, options->patternRowIter, - options->patternRowRej, options->patternRowThresh, options->patternRowMean, - options->patternRowStdev, options->maskValue, options->darkMask)) { - psFree (view); - return false; - } - } - psFree (view); // supplied view must be freeable + pmReadout *readout; // Readout from cell + while ((readout = pmFPAviewNextReadout(view, fpa, 1)) != NULL) + { + if (!readout->data_exists) + { + continue; + } + + // Perform pattern correction + if (!pmPatternRow(readout, options->patternRowOrder, options->patternRowIter, + options->patternRowRej, options->patternRowThresh, options->patternRowMean, + options->patternRowStdev, options->maskValue, options->darkMask)) + { + psFree(view); + return false; + } + } + psFree(view); // supplied view must be freeable return true; } -psVector *ppImageDetrendPatternCellFailures(pmFPAfile *input, const pmFPAview *inputView) { - +psVector *ppImageDetrendPatternCellFailures(pmFPAfile *input, const pmFPAview *inputView) +{ + bool status = false; pmCell *cell = NULL; - + pmFPAview *view = pmFPAviewAlloc(0); // View for local processing *view = *inputView; // extract the MEDIAN_CELL_BACKGROUND values to check for problems - psVector *cellBackground = psVectorAllocEmpty (64, PS_DATA_F32); - while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) { - if (!cell->process || !cell->file_exists || !cell->data_exists) { - psVectorAppend (cellBackground, NAN); - continue; - } - - // select the HDU for this cell - pmHDU *hdu = pmHDUFromCell(cell); // HDU of interest - - psF32 value = psMetadataLookupF32(&status, hdu->header, "BACK_VAL"); - psVectorAppend (cellBackground, value); + psVector *cellBackground = psVectorAllocEmpty(64, PS_DATA_F32); + while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) + { + if (!cell->process || !cell->file_exists || !cell->data_exists) + { + psVectorAppend(cellBackground, NAN); + continue; + } + + // select the HDU for this cell + pmHDU *hdu = pmHDUFromCell(cell); // HDU of interest + + psF32 value = psMetadataLookupF32(&status, hdu->header, "BACK_VAL"); + psVectorAppend(cellBackground, value); } return cellBackground; Index: /branches/2dbias/ppImage/src/ppImageDetrendReadout.c =================================================================== --- /branches/2dbias/ppImage/src/ppImageDetrendReadout.c (revision 42678) +++ /branches/2dbias/ppImage/src/ppImageDetrendReadout.c (revision 42679) @@ -5,7 +5,14 @@ #include "ppImage.h" +#define ESCAPE(STATUS, ...) \ + { \ + psError(PS_ERR_UNKNOWN, STATUS, __VA_ARGS__); \ + psFree(view); \ + return false; \ + } + bool ppImageDetrendReadout(pmConfig *config, ppImageOptions *options, pmFPAview *view) { - // psTimerStart("detrend.readout"); + // psTimerStart("detrend.readout"); // construct a view for the detrend images (which have only one readout) @@ -18,22 +25,26 @@ // Check that the gain is set (this is used by both pmReadoutGenerateMask and pmReadoutGenerateVariance) - { - float gain = psMetadataLookupF32(NULL, input->parent->concepts, "CELL.GAIN"); // Gain for cell - if (!isfinite(gain)) { - psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, RECIPE_NAME); // Recipe - psAssert(recipe, "Should be there!"); - bool override = psMetadataLookupBool(NULL, recipe, "GAIN.OVERRIDE"); // Override the bad gain? - if (override) { - psWarning("CELL.GAIN is not set for readout (%d,%d,%d) --- setting to unity.", view->chip, view->cell, view->readout); - psMetadataItem *item = psMetadataLookup(input->parent->concepts, "CELL.GAIN"); // Gain item - psAssert(item, "Should be there!"); - item->data.F32 = 1.0; - - // for unity gain, there is no modification for the readnoise, note that it has (effectively) been updated - psMetadataRemoveKey(input->parent->concepts, "CELL.READNOISE.UPDATE"); - } else { - psWarning("CELL.GAIN is NAN for readout (%d,%d,%d), image will be masked.", view->chip, view->cell, view->readout); - } - } + { + float gain = psMetadataLookupF32(NULL, input->parent->concepts, "CELL.GAIN"); // Gain for cell + if (!isfinite(gain)) + { + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, RECIPE_NAME); // Recipe + psAssert(recipe, "Should be there!"); + bool override = psMetadataLookupBool(NULL, recipe, "GAIN.OVERRIDE"); // Override the bad gain? + if (override) + { + psWarning("CELL.GAIN is not set for readout (%d,%d,%d) --- setting to unity.", view->chip, view->cell, view->readout); + psMetadataItem *item = psMetadataLookup(input->parent->concepts, "CELL.GAIN"); // Gain item + psAssert(item, "Should be there!"); + item->data.F32 = 1.0; + + // for unity gain, there is no modification for the readnoise, note that it has (effectively) been updated + psMetadataRemoveKey(input->parent->concepts, "CELL.READNOISE.UPDATE"); + } + else + { + psWarning("CELL.GAIN is NAN for readout (%d,%d,%d), image will be masked.", view->chip, view->cell, view->readout); + } + } } // psLogMsg ("ppImage", 6, "check gain: %f sec\n", psTimerMark ("detrend.readout")); @@ -43,22 +54,29 @@ bool hasVideo = false; { - // XXX test for GPC1? or CAN_HAVE_VIDEO in camera config? - bool status = false; - if (!input) goto done_video_check; - if (!input->parent) goto done_video_check; - if (!input->parent->parent) goto done_video_check; - if (!input->parent->parent->hdu) goto done_video_check; - if (!input->parent->parent->hdu->header) goto done_video_check; - char *ptr = psMetadataLookupStr(&status,input->parent->parent->hdu->header,"CELLMODE"); - if (status) { - psLogMsg ("ppImage.detrend", PS_LOG_DETAIL, "VIDEO: %d %d %d\n",(int) options->hasVideo,(int) options->useVideoDark, (int) options->useVideoMask); - char *Vptr = strchr(ptr,'V'); - if (Vptr) { - hasVideo = options->hasVideo = true; - psLogMsg ("ppImage.detrend", PS_LOG_INFO, "VIDEO: %d %d %d\n",(int) options->hasVideo,(int) options->useVideoDark, (int) options->useVideoMask); - } - } - } - done_video_check: + // XXX test for GPC1? or CAN_HAVE_VIDEO in camera config? + bool status = false; + if (!input) + goto done_video_check; + if (!input->parent) + goto done_video_check; + if (!input->parent->parent) + goto done_video_check; + if (!input->parent->parent->hdu) + goto done_video_check; + if (!input->parent->parent->hdu->header) + goto done_video_check; + char *ptr = psMetadataLookupStr(&status, input->parent->parent->hdu->header, "CELLMODE"); + if (status) + { + psLogMsg("ppImage.detrend", PS_LOG_DETAIL, "VIDEO: %d %d %d\n", (int)options->hasVideo, (int)options->useVideoDark, (int)options->useVideoMask); + char *Vptr = strchr(ptr, 'V'); + if (Vptr) + { + hasVideo = options->hasVideo = true; + psLogMsg("ppImage.detrend", PS_LOG_INFO, "VIDEO: %d %d %d\n", (int)options->hasVideo, (int)options->useVideoDark, (int)options->useVideoMask); + } + } + } +done_video_check: // psLogMsg ("ppImage", 6, "check video: %f sec\n", psTimerMark ("detrend.readout")); @@ -74,25 +92,30 @@ // Masking on the basis of pixel value needs to be done before anything else, so the values are pristine. - if (options->doMaskBuild) { + if (options->doMaskBuild) + { psImageMaskType satMask = options->doMaskSat ? options->satMask : 0; psImageMaskType lowMask = options->doMaskLow ? options->lowMask : 0; - if (!pmReadoutGenerateMask(input, satMask, lowMask)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate a mask."); - psFree(detview); - return false; - } - // psLogMsg ("ppImage", 6, "generate mask: %f sec\n", psTimerMark ("detrend.readout")); + if (!pmReadoutGenerateMask(input, satMask, lowMask)) + { + psError(PS_ERR_UNKNOWN, false, "Unable to generate a mask."); + psFree(detview); + return false; + } + // psLogMsg ("ppImage", 6, "generate mask: %f sec\n", psTimerMark ("detrend.readout")); } // apply the externally supplied mask to the input->mask pixels - if (options->doMask) { - pmReadout *mask; - if ((options->useVideoMask)&&(hasVideo)) { - mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.VIDEOMASK"); - } - else { - mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.MASK"); - } - pmMaskBadPixels(input, mask, options->maskValue); - // psLogMsg ("ppImage", 6, "apply mask: %f sec\n", psTimerMark ("detrend.readout")); + if (options->doMask) + { + pmReadout *mask; + if ((options->useVideoMask) && (hasVideo)) + { + mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.VIDEOMASK"); + } + else + { + mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.MASK"); + } + pmMaskBadPixels(input, mask, options->maskValue); + // psLogMsg ("ppImage", 6, "apply mask: %f sec\n", psTimerMark ("detrend.readout")); } @@ -101,116 +124,150 @@ // extern bool ppImageBurntoolApply(pmConfig *, ppImageOptions *, pmFPAview *, pmReadout *); ppImageBurntoolApply(config, options, view, input); - // psLogMsg ("ppImage", 6, "apply burntool: %f sec\n", psTimerMark ("detrend.readout")); - } - - if (options->doMaskBurntool) { - if (options->doApplyBurntool) { + // psLogMsg ("ppImage", 6, "apply burntool: %f sec\n", psTimerMark ("detrend.readout")); + } + + if (options->doMaskBurntool) + { + if (options->doApplyBurntool) + { // build burntool mask from data input burntool table - ppImageBurntoolMaskFromTable(config,options,view,input); - } else { + ppImageBurntoolMaskFromTable(config, options, view, input); + } + else + { // build burntool mask from data in the input image's fits extension - ppImageBurntoolMask(config,options,view,input); - } - // psLogMsg ("ppImage", 6, "apply burntool mask: %f sec\n", psTimerMark ("detrend.readout")); + ppImageBurntoolMask(config, options, view, input); + } + // psLogMsg ("ppImage", 6, "apply burntool mask: %f sec\n", psTimerMark ("detrend.readout")); } // Subtract the overscan - if (options->doOverscan) { - if (!pmOverscanSubtract (input, options->overscan)) { - psError(PS_ERR_UNKNOWN, false, "Unable to subtract overscan."); - psFree(detview); - return false; - } - // psLogMsg ("ppImage", 6, "subtract overscan: %f sec\n", psTimerMark ("detrend.readout")); + if (options->doOverscan) + { + bool doTwoDOverscan = false; + if (options->overscan->TwoD) + { + pmChip *chip = pmFPAfileThisChip(config->files, view, "PPIMAGE.INPUT"); + if (!ppImageDoPatternForView(&doTwoDOverscan, config, chip, view, RECIPE_NAME, "OVERSCAN.2D.SUBSET")) + { + ESCAPE(false, "Unable to determine whether 2D Overscan subtraction should be applied."); + } + } + if (!pmOverscanSubtract(input, options->overscan, doTwoDOverscan)) + { + psError(PS_ERR_UNKNOWN, false, "Unable to subtract overscan."); + psFree(detview); + return false; + } + // psLogMsg ("ppImage", 6, "subtract overscan: %f sec\n", psTimerMark ("detrend.readout")); } // measure the overscan-subtracted readoutBackground here (or subtract the overscan value?) // XXX this is the measurements and should be independent of the pattern masking - if (options->doPatternDeadCells) { - psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); - psImageBackground (stats, NULL, input->image, NULL, 0xffff, rng); - - // save this value somewhere - pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "BACK_VAL", PS_META_REPLACE, "Median cell background", stats->robustMedian); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "BACK_ERR", PS_META_REPLACE, "Stdev of cell background", stats->robustStdev); - psFree (stats); - psFree (rng); + if (options->doPatternDeadCells) + { + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + psImageBackground(stats, NULL, input->image, NULL, 0xffff, rng); + + // save this value somewhere + pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "BACK_VAL", PS_META_REPLACE, "Median cell background", stats->robustMedian); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "BACK_ERR", PS_META_REPLACE, "Stdev of cell background", stats->robustStdev); + psFree(stats); + psFree(rng); } // Non-linearity correction - if (options->doNonLin) { - if (!ppImageDetrendNonLinear(input,detview,config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to correct NonLinearity"); - psFree(detview); - return(false); - } - // psLogMsg ("ppImage", 6, "nonlinear correction: %f sec\n", psTimerMark ("detrend.readout")); + if (options->doNonLin) + { + if (!ppImageDetrendNonLinear(input, detview, config)) + { + psError(PS_ERR_UNKNOWN, false, "Unable to correct NonLinearity"); + psFree(detview); + return (false); + } + // psLogMsg ("ppImage", 6, "nonlinear correction: %f sec\n", psTimerMark ("detrend.readout")); } // New Non-linearity correction (exclusive of the above) - if (options->doNewNonLin) { - if (!ppImageDetrendNewNonLinear(input, detview, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to correct Non-Linearity with new version (2023)"); - psFree(detview); - return(false); - } - // psLogMsg ("ppImage", 6, "nonlinear correction: %f sec\n", psTimerMark ("detrend.readout")); + if (options->doNewNonLin) + { + if (!ppImageDetrendNewNonLinear(input, detview, config)) + { + psError(PS_ERR_UNKNOWN, false, "Unable to correct Non-Linearity with new version (2023)"); + psFree(detview); + return (false); + } + // psLogMsg ("ppImage", 6, "nonlinear correction: %f sec\n", psTimerMark ("detrend.readout")); } // set up the dark and bias - pmCell *dark = NULL; // Multi-dark - pmReadout *oldDark = NULL; // Old-fashioned dark - pmReadout *bias = NULL; - if (options->doBias) { + pmCell *dark = NULL; // Multi-dark + pmReadout *oldDark = NULL; // Old-fashioned dark + pmReadout *bias = NULL; + if (options->doBias) + { bias = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.BIAS"); } - if (options->doDark) { - bool mdok; // Status of MD lookup - psMetadata *recipe = psMetadataLookupPtr (&mdok, config->recipes, RECIPE_NAME); + if (options->doDark) + { + bool mdok; // Status of MD lookup + psMetadata *recipe = psMetadataLookupPtr(&mdok, config->recipes, RECIPE_NAME); assert(mdok && recipe); - if ((options->useVideoDark)&&(hasVideo)) { - if (psMetadataLookupBool(&mdok, recipe, "OLDDARK")) { - oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.VIDEODARK"); - } else { - dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.VIDEODARK"); - } - } - else { - if (psMetadataLookupBool(&mdok, recipe, "OLDDARK")) { - oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.DARK"); - } else { - dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.DARK"); - } - } + if ((options->useVideoDark) && (hasVideo)) + { + if (psMetadataLookupBool(&mdok, recipe, "OLDDARK")) + { + oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.VIDEODARK"); + } + else + { + dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.VIDEODARK"); + } + } + else + { + if (psMetadataLookupBool(&mdok, recipe, "OLDDARK")) + { + oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.DARK"); + } + else + { + dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.DARK"); + } + } } // Bias and temperature-independent-dark subtraction are merged. - if (options->doBias) { - if (!pmBiasSubtract(input, bias, oldDark, view)) { + if (options->doBias) + { + if (!pmBiasSubtract(input, bias, oldDark, view)) + { psError(PS_ERR_UNKNOWN, false, "Unable to subtract bias."); psFree(detview); return false; } - // psLogMsg ("ppImage", 6, "apply bias: %f sec\n", psTimerMark ("detrend.readout")); - } - + // psLogMsg ("ppImage", 6, "apply bias: %f sec\n", psTimerMark ("detrend.readout")); + } + // Weight on the basis of pixel value needs to be done after the overscan has been subtracted - if (options->doVarianceBuild) { + if (options->doVarianceBuild) + { // create the target mask and variance images psImage *noiseImage = NULL; - if (options->doNoiseMap) { + if (options->doNoiseMap) + { // XXX convert the noiseMap image to a binned image pmReadout *noiseMap = NULL; noiseMap = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.NOISEMAP"); - noiseImage = psImageCopy (NULL, input->image, PS_TYPE_F32); - psImageInit (noiseImage, 0.0); + noiseImage = psImageCopy(NULL, input->image, PS_TYPE_F32); + psImageInit(noiseImage, 0.0); // XXX this works, but is not really quite right: the model shoud include the // offset information, we are not really getting exactly the right mapping from the // original file. - // CZW 2012-03-21: I do not believe this is true anymore. In any case, this seems to be what - // ppImageReplaceBackground does to do sky subtraction. + // CZW 2012-03-21: I do not believe this is true anymore. In any case, this seems to be what + // ppImageReplaceBackground does to do sky subtraction. psImageBinning *binning = psImageBinningAlloc(); binning->nXruff = noiseMap->image->numCols; @@ -220,72 +277,82 @@ psImageBinningSetScale(binning, PS_IMAGE_BINNING_LEFT); - psImageUnbin (noiseImage, noiseMap->image, binning); - psFree (binning); - } - - if (!pmReadoutGenerateVariance(input, noiseImage, true)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate a variance image."); - psFree(detview); - return false; - } - psFree (noiseImage); - // psLogMsg ("ppImage", 6, "generate variance: %f sec\n", psTimerMark ("detrend.readout")); - } - - if (options->doDark && dark) { - if (!pmDarkApply(input, dark, options->darkMask)) { + psImageUnbin(noiseImage, noiseMap->image, binning); + psFree(binning); + } + + if (!pmReadoutGenerateVariance(input, noiseImage, true)) + { + psError(PS_ERR_UNKNOWN, false, "Unable to generate a variance image."); + psFree(detview); + return false; + } + psFree(noiseImage); + // psLogMsg ("ppImage", 6, "generate variance: %f sec\n", psTimerMark ("detrend.readout")); + } + + if (options->doDark && dark) + { + if (!pmDarkApply(input, dark, options->darkMask)) + { psError(PS_ERR_UNKNOWN, false, "Unable to subtract dark."); psFree(detview); return false; } - // psLogMsg ("ppImage", 6, "apply dark: %f sec\n", psTimerMark ("detrend.readout")); - } - - if (options->doRemnance) { + // psLogMsg ("ppImage", 6, "apply dark: %f sec\n", psTimerMark ("detrend.readout")); + } + + if (options->doRemnance) + { if (!pmRemnance(input, options->maskValue, options->lowMask, - options->remnanceSize, options->remnanceThresh)) { + options->remnanceSize, options->remnanceThresh)) + { psError(PS_ERR_UNKNOWN, false, "Unable to mask remnance."); psFree(detview); return false; } - // psLogMsg ("ppImage", 6, "mask remnance: %f sec\n", psTimerMark ("detrend.readout")); + // psLogMsg ("ppImage", 6, "mask remnance: %f sec\n", psTimerMark ("detrend.readout")); } // Shutter correction - if (options->doShutter) { + if (options->doShutter) + { pmReadout *shutter = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.SHUTTER"); - if (!pmShutterCorrectionApply(input, shutter, pmConfigMaskGet("FLAT", config))) { - psFree(detview); - return false; - } - // psLogMsg ("ppImage", 6, "shutter correction: %f sec\n", psTimerMark ("detrend.readout")); + if (!pmShutterCorrectionApply(input, shutter, pmConfigMaskGet("FLAT", config))) + { + psFree(detview); + return false; + } + // psLogMsg ("ppImage", 6, "shutter correction: %f sec\n", psTimerMark ("detrend.readout")); } // Flat-field correction (no options used?) - if (options->doFlat) { + if (options->doFlat) + { pmReadout *flat = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.FLAT"); - if (!pmFlatField(input, flat, options->flatMask)) { - psFree(detview); - return false; - } - // psLogMsg ("ppImage", 6, "apply flat: %f sec\n", psTimerMark ("detrend.readout")); - } - -/* // Pattern noise correction */ -/* if (options->doPattern) { */ -/* if (!pmPatternRow(input, options->patternOrder, options->patternIter, options->patternRej, */ -/* options->patternThresh, options->patternMean, options->patternStdev, */ -/* options->maskValue, options->darkMask)) { */ -/* psFree(detview); */ -/* return false; */ -/* } */ -/* } */ + if (!pmFlatField(input, flat, options->flatMask)) + { + psFree(detview); + return false; + } + // psLogMsg ("ppImage", 6, "apply flat: %f sec\n", psTimerMark ("detrend.readout")); + } + + /* // Pattern noise correction */ + /* if (options->doPattern) { */ + /* if (!pmPatternRow(input, options->patternOrder, options->patternIter, options->patternRej, */ + /* options->patternThresh, options->patternMean, options->patternStdev, */ + /* options->maskValue, options->darkMask)) { */ + /* psFree(detview); */ + /* return false; */ + /* } */ + /* } */ // Normalization by a single (known) constant - bool mdok; // Status of MD lookup + bool mdok; // Status of MD lookup float norm = psMetadataLookupF32(&mdok, config->arguments, "NORMALIZATION"); - if (mdok && isfinite(norm) && norm != 1.0) { + if (mdok && isfinite(norm) && norm != 1.0) + { pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest - psString comment = NULL; // Comment to add + psString comment = NULL; // Comment to add psStringAppend(&comment, "Normalization: %f", norm); psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); @@ -293,16 +360,18 @@ psBinaryOp(input->image, input->image, "*", psScalarAlloc(norm, PS_TYPE_F32)); - // psLogMsg ("ppImage", 6, "renormalize: %f sec\n", psTimerMark ("detrend.readout")); - } - -# if (1) + // psLogMsg ("ppImage", 6, "renormalize: %f sec\n", psTimerMark ("detrend.readout")); + } + +#if (1) // Normalization by per-class values psMetadata *normlist = psMetadataLookupMetadata(&mdok, config->arguments, "NORMALIZATION.TABLE"); - if (normlist) { + if (normlist) + { pmFPAfile *inputFile = psMetadataLookupPtr(&mdok, config->files, "PPIMAGE.INPUT"); // get the menu of class IDs psMetadata *menu = psMetadataLookupMetadata(&mdok, inputFile->camera, "CLASSID"); - if (!menu) { + if (!menu) + { psError(PS_ERR_IO, false, "Unable to find CLASSID metadata in camera configuration"); psFree(detview); @@ -311,5 +380,6 @@ // get the rule for class_id for the desired class const char *rule = psMetadataLookupStr(&mdok, menu, options->normClass); - if (!rule) { + if (!rule) + { psError(PS_ERR_IO, false, "Unable to find NORM.CLASS value %s in CLASSID in camera configuration", options->normClass); psFree(detview); @@ -318,5 +388,6 @@ // get the class_id from the rule char *classID = pmFPAfileNameFromRule(rule, inputFile, view); - if (!classID) { + if (!classID) + { psError(PS_ERR_IO, false, "error converting CLASSID rule %s to name\n", rule); psFree(detview); @@ -325,6 +396,7 @@ // get normalization from the class_id - float norm = psMetadataLookupF32 (&mdok, normlist, classID); - if (!mdok) { + float norm = psMetadataLookupF32(&mdok, normlist, classID); + if (!mdok) + { psError(PS_ERR_IO, false, "failed to find class ID %s in normalization table\n", classID); psFree(detview); @@ -333,5 +405,5 @@ pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest - psString comment = NULL; // Comment to add + psString comment = NULL; // Comment to add psStringAppend(&comment, "Normalization: %f", norm); psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); @@ -341,17 +413,19 @@ psBinaryOp(input->image, input->image, "*", psScalarAlloc(norm, PS_TYPE_F32)); - psFree (classID); - - // psLogMsg ("ppImage", 6, "renormalize by class: %f sec\n", psTimerMark ("detrend.readout")); - } -# endif - - if (options->doFringe) { + psFree(classID); + + // psLogMsg ("ppImage", 6, "renormalize by class: %f sec\n", psTimerMark ("detrend.readout")); + } +#endif + + if (options->doFringe) + { pmCell *fringe = pmFPAfileThisCell(config->files, detview, "PPIMAGE.FRINGE"); - if (!ppImageDetrendFringeMeasure(input, fringe, false, options)) { - psFree(detview); - return false; - } - // psLogMsg ("ppImage", 6, "measure fringe: %f sec\n", psTimerMark ("detrend.readout")); + if (!ppImageDetrendFringeMeasure(input, fringe, false, options)) + { + psFree(detview); + return false; + } + // psLogMsg ("ppImage", 6, "measure fringe: %f sec\n", psTimerMark ("detrend.readout")); } Index: /branches/2dbias/psModules/src/detrend/pmOverscan.c =================================================================== --- /branches/2dbias/psModules/src/detrend/pmOverscan.c (revision 42678) +++ /branches/2dbias/psModules/src/detrend/pmOverscan.c (revision 42679) @@ -17,419 +17,482 @@ #include "pmOverscan.h" -#define SMOOTH_NSIGMA 4.0 // Number of Gaussian sigma the smoothing kernel extends +#define SMOOTH_NSIGMA 4.0 // Number of Gaussian sigma the smoothing kernel extends static void pmOverscanOptionsFree(pmOverscanOptions *options); static void pmOverscanStatOptionsFree(pmOverscanStatOptions *options); -bool pmOverscanUpdateHeaderVector (pmReadout *input, pmHDU *hdu, pmOverscanOptions *overscanOpts, psVector *reduced); - -bool pmOverscanSubtract (pmReadout *input, pmOverscanOptions *overscanOpts) { - - assert (input); - - if (overscanOpts == NULL) return true; // no overscan subtraction requested - - pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest - psImage *image = input->image; - - // check for 'soft bias' (simple, fixed offset to be subtracted) - if (overscanOpts->constant) { - // write metadata header value - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", overscanOpts->value); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN); - - // NOTE psBinaryOp frees arg2 if it is a scalar - (void)psBinaryOp(input->image, input->image, "-", psScalarAlloc((float)overscanOpts->value, PS_TYPE_F32)); - +bool pmOverscanUpdateHeaderVector(pmReadout *input, pmHDU *hdu, pmOverscanOptions *overscanOpts, psVector *reduced); + +bool pmOverscanSubtract(pmReadout *input, pmOverscanOptions *overscanOpts, bool doTwoDOverscan) +{ + + assert(input); + + if (overscanOpts == NULL) + return true; // no overscan subtraction requested + + pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest + psImage *image = input->image; + + // check for 'soft bias' (simple, fixed offset to be subtracted) + if (overscanOpts->constant) + { + // write metadata header value + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", overscanOpts->value); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN); + + // NOTE psBinaryOp frees arg2 if it is a scalar + (void)psBinaryOp(input->image, input->image, "-", psScalarAlloc((float)overscanOpts->value, PS_TYPE_F32)); + + return true; + } + + // we are performing a statitical analysis of the overscan region + + // Check for an unallowable pmFit. + if (overscanOpts->primary) + { + if (overscanOpts->primary->fitType != PM_FIT_NONE && + overscanOpts->primary->fitType != PM_FIT_POLY_ORD && + overscanOpts->primary->fitType != PM_FIT_POLY_CHEBY && + overscanOpts->primary->fitType != PM_FIT_SPLINE) + { + psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d). Returning original image.\n", + overscanOpts->primary->fitType); + return false; + } + } + if (overscanOpts->secondary) + { + if (overscanOpts->secondary->fitType != PM_FIT_NONE && + overscanOpts->secondary->fitType != PM_FIT_POLY_ORD && + overscanOpts->secondary->fitType != PM_FIT_POLY_CHEBY && + overscanOpts->secondary->fitType != PM_FIT_SPLINE) + { + psError(PS_ERR_UNKNOWN, true, "Invalid fit type (2D) (%d). Returning original image.\n", + overscanOpts->secondary->fitType); + return false; + } + } + + psList *overscans = input->bias; // List of the overscan images + + // Reduce all overscan pixels to a single value + if (overscanOpts->single) + { + + // extract overscan pixels to a single vector + psVector *pixels = psVectorAlloc(0, PS_TYPE_F32); + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) + { + int index = pixels->n; // Index + pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols); + pixels->n += overscan->numRows * overscan->numCols; + for (int i = 0; i < overscan->numRows; i++) + { + memcpy(&pixels->data.F32[index], overscan->data.F32[i], + overscan->numCols * sizeof(psF32)); + index += overscan->numCols; + } + } + psFree(iter); + + // statistic to be calculated + psStatsOptions statistic = psStatsSingleOption(overscanOpts->primary->stat->options); // Statistic to use + if (!statistic) + { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n", + overscanOpts->primary->stat); + return false; + } + psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original + + if (!psVectorStats(stats, pixels, NULL, NULL, 0)) + { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } + psFree(pixels); + double reduced = psStatsGetValue(stats, statistic); // Result of statistics + + psString comment = NULL; // Comment to add + psStringAppend(&comment, "Overscan value: %f", reduced); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + + // write metadata header value + // XXX EAM : this could / should write the stdev of the overscan region + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", reduced); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN); + + psScalar *reducedScalar = psScalarAlloc(reduced, PS_TYPE_F32); + psBinaryOp(image, image, "-", psMemIncrRefCounter(reducedScalar)); // NOTE: psBinaryOp frees arg2 if it a scalar, so we need to bump to re-use + + // subtract the measured value from each overscan region as well + iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) + { + psBinaryOp(overscan, overscan, "-", psMemIncrRefCounter(reducedScalar)); // NOTE: psBinaryOp frees arg2 if it a scalar, so we need to bump to re-use + } + psFree(iter); + psFree(reducedScalar); + + // EAM 2022.03.29 : if the calculated overscan value is below the threshold, + // declare the readout dead and mask + + if ((reduced < overscanOpts->minValid) || (reduced > overscanOpts->maxValid)) + { + fprintf(stderr, "bad overscan (1) %f, masking readout\n", reduced); + psImage *mask = input->mask; + for (int y = 0; y < mask->numRows; y++) + { + for (int x = 0; x < mask->numCols; x++) + { + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= overscanOpts->maskVal; + } + } + } + + psFree(stats); + return true; + } + + bool mdok = false; + + // We are performing a row-by-row overscan subtraction + int cellreaddir = psMetadataLookupS32(&mdok, input->parent->concepts, "CELL.READDIR"); // Read direction + if ((cellreaddir != 1) && (cellreaddir != 2)) + { + psError(PS_ERR_UNKNOWN, true, "CELL.READDIR must be 1 (rows) or 2 (cols)\n"); + return false; + } + + float chi2 = NAN; // chi^2 from fit + + // adjust operation depending on the read direction : need to re-org pixels for columns + if (!doTwoDOverscan && (cellreaddir == 1)) + { + // The read direction is rows + psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels + for (int i = 0; i < pixels->n; i++) + { + pixels->data[i] = psVectorAlloc(0, PS_TYPE_F32); + } + + // Pull the pixels out into the vectors + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) + { + // the overscan and image might not be aligned. pixels->data represents + // the image row pixels. + int diff = overscan->row0 - image->row0; // Offset between the two regions + for (int i = PS_MAX(0, diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++) + { + int j = i - diff; + // i is row on image + // j is row on overscan + psVector *values = pixels->data[i]; + int index = values->n; // Index in the vector + values = psVectorRealloc(values, values->n + overscan->numCols); + values->n += overscan->numCols; + // XXX double-check the range of values->n here + memcpy(&values->data.F32[index], overscan->data.F32[j], + overscan->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + index += overscan->numCols; + pixels->data[i] = values; // Update the pointer in case it's moved + } + } + psFree(iter); + + // Reduce the overscans + psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels); + psFree(pixels); + if (!reduced) + { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n"); + return false; + } + + if (!pmOverscanUpdateHeaderVector(input, hdu, overscanOpts, reduced)) + { + psError(PS_ERR_UNKNOWN, false, "failure to update header"); + return false; + } + + // Subtract row by row + for (int i = 0; i < image->numRows; i++) + { + for (int j = 0; j < image->numCols; j++) + { + image->data.F32[i][j] -= reduced->data.F32[i]; + } + } + + // subtract from the overscan regions + { + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) + { + // the overscan and image might not be aligned. + int diff = overscan->row0 - image->row0; // Offset between the two regions + for (int i = PS_MAX(0, diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++) + { + int j = i - diff; + // i is row on image + // j is row on overscan + for (int k = 0; k < overscan->numCols; k++) + { + overscan->data.F32[j][k] -= reduced->data.F32[j]; + } + } + } + psFree(iter); + } + psFree(reduced); + } + + if (!doTwoDOverscan && (cellreaddir == 2)) + { + // The read direction is columns + psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels + for (int i = 0; i < pixels->n; i++) + { + psVector *values = psVectorAlloc(0, PS_TYPE_F32); + pixels->data[i] = values; + } + + // Pull the pixels out into the vectors + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) + { + // the overscan and image might not be aligned. pixels->data represents + // the image row pixels. + int diff = overscan->col0 - image->col0; // Offset between the two regions + for (int i = PS_MAX(0, diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++) + { + int iFixed = i - diff; + // i is column on image + // iFixed is column on overscan + psVector *values = pixels->data[i]; + int index = values->n; // Index in the vector + values = psVectorRealloc(values, values->n + overscan->numRows); + for (int j = 0; j < overscan->numRows; j++) + { + values->data.F32[index++] = overscan->data.F32[j][iFixed]; + } + values->n += overscan->numRows; + pixels->data[i] = values; // Update the pointer in case it's moved + } + } + psFree(iter); + + // Reduce the overscans + psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels); + psFree(pixels); + if (!reduced) + { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n"); + return false; + } + + if (!pmOverscanUpdateHeaderVector(input, hdu, overscanOpts, reduced)) + { + psError(PS_ERR_UNKNOWN, false, "failure to update header"); + return false; + } + + // Subtract column by column + for (int j = 0; j < image->numRows; j++) + { + for (int i = 0; i < image->numCols; i++) + { + image->data.F32[j][i] -= reduced->data.F32[i]; + } + } + + // subtract from the overscan regions + { + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) + { + // the overscan and image might not be aligned. + int diff = overscan->col0 - image->col0; // Offset between the two regions + for (int i = PS_MAX(0, diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++) + { + int j = i - diff; + // i is col on image + // j is col on overscan + for (int k = 0; i < overscan->numRows; k++) + { + overscan->data.F32[k][j] -= reduced->data.F32[j]; + } + } + } + psFree(iter); + } + + psFree(reduced); + } + + // 2D bias subtraction with x-dir readout direction: the + // bias is constructed by combining a 1D pattern in the + // readout direction from the top overscan region and a second + // 1D pattern in the cross direction from the overscan + if (doTwoDOverscan && (cellreaddir == 1)) + { + psAssert(overscanOpts->secondary, "2D overscan subtraction requires OVERSCAN.2D parameters"); + // we require 2 overscan regions, and they must match these directions: + if (overscans->n != 2) + { + psLogMsg(__func__, PS_LOG_ERROR, "OVERSCAN.2D requires 2 overscan regions but %d supplied", (int)overscans->n); + psLogMsg(__func__, PS_LOG_ERROR, "e.g.: CELL.BIASSEC STR [591:624,1:608],[1:624,599:608]"); + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to generate overscan vector.\n"); + return false; + } + + // the serial (fast readout) direction is columns (x-dir) + psImage *yscan = psListGet(overscans, 0); // overscan region spanning all rows + psImage *xscan = psListGet(overscans, 1); // overscan region spanning all columns + + // Extract the y-dir overscan vector. The overscan and image might not be aligned: + // diff represents the offset between the rows in the image data and the overscan. + // pixels->data represents the image row pixels. For example, the image region may be + // inset in the y-direction but the overscan could cover the entire y-range + + // The read direction is rows + psArray *yscanPixels = psArrayAlloc(yscan->numRows); // Array of vectors containing pixels + for (int i = 0; i < yscanPixels->n; i++) + { + yscanPixels->data[i] = psVectorAlloc(0, PS_TYPE_F32); + } + + // XXX this code allows multiple yscans to be appended, but this does not + // match the concept of how they are assigned above: biassec[0] = yscan + // int yDiff = yscan->row0 - image->row0; // Offset between the two regions + for (int i = 0; i < yscanPixels->n; i++) + { + psVector *values = yscanPixels->data[i]; + int index = values->n; // Index in the vector + values = psVectorRealloc(values, values->n + yscan->numCols); + values->n += yscan->numCols; + // XXX double-check the range of values->n here + memcpy(&values->data.F32[index], yscan->data.F32[i], + yscan->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + yscanPixels->data[i] = values; // Update the pointer in case it's moved + } + + // Extract the x-dir overscan vector. The overscan and image might not be aligned: + // diff represents the offset between the rows in the image data and the overscan. + // pixels->data represents the image row pixels. For example, the image region may be + // inset in the x-direction but the overscan could cover the entire y-range + + // Extract the top region as a vector of the columns + psArray *xscanPixels = psArrayAlloc(xscan->numCols); // Array of vectors containing pixels + for (int i = 0; i < xscanPixels->n; i++) + { + xscanPixels->data[i] = psVectorAlloc(0, PS_TYPE_F32); + } + + // int xDiff = xscan->col0 - image->col0; // Offset between the two regions + for (int ix = 0; ix < xscanPixels->n; ix++) + { + psVector *values = xscanPixels->data[ix]; + values = psVectorRealloc(values, xscan->numRows); + values->n = xscan->numRows; + for (int iy = 0; iy < xscan->numRows; iy++) + { + values->data.F32[iy] = xscan->data.F32[iy][ix]; + } + xscanPixels->data[ix] = values; // Update the pointer in case it's moved + } + + // Reduce the overscans + // XXX need to save 2 different chi-square values + psVector *yReduced = pmOverscanVector(&chi2, overscanOpts->primary, yscanPixels); + psFree(yscanPixels); + if (!yReduced) + { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate y-dir overscan vector.\n"); + return false; + } + if (!pmOverscanUpdateHeaderVector(input, hdu, overscanOpts, yReduced)) + { + psError(PS_ERR_UNKNOWN, false, "failure to update header"); + return false; + } + + // Reduce the overscans + // XXX need to save 2 different chi-square values + psVector *xReduced = pmOverscanVector(&chi2, overscanOpts->secondary, xscanPixels); + psFree(xscanPixels); + if (!xReduced) + { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate x-dir overscan vector.\n"); + return false; + } + + // subtract the 2D bias from the image + if (yscan->col0 >= xscan->col0 && yscan->col0 + yscan->numCols <= xscan->col0 + xscan->numCols) + { + // define where to normalize the xReduced vector + int j0 = yscan->col0 - xscan->col0 + (int)(yscan->numCols / 2) - 1; + // XXX apply the 2D bias correction here + int yDiff = yscan->row0 - image->row0; // y offset between the science and the yoverscan region + int xDiff = xscan->col0 - image->col0; // x offset between the science and the xoverscan region + for (int i = 0; i < image->numRows; i++) + { + for (int j = 0; j < image->numCols; j++) + { + int iy = i + yDiff; + int jx = j + xDiff; + image->data.F32[i][j] -= yReduced->data.F32[iy] - xReduced->data.F32[j0] + xReduced->data.F32[jx]; + } + } + } + else + { + psError(PS_ERR_UNKNOWN, true, "x dimension of yscan is not fully contained by xscan\n"); + return false; + } + + // subtract the y-dir vector from the y-dir overscan regions (why?) + { + // the overscan and image might not be aligned. + // int diff = yscan->row0 - image->row0; // Offset between the two regions + for (int i = 0; i < yscan->numRows; i++) + { + for (int j = 0; j < yscan->numCols; j++) + { + yscan->data.F32[i][j] -= yReduced->data.F32[i]; + } + } + } + + // subtract the x-dir vector from the x-dir overscan regions (why?) + { + // the overscan and image might not be aligned. + // int diff = xscan->col0 - image->col0; // Offset between the two regions + for (int i = 0; i < xscan->numCols; i++) + { + // int j = i - diff; + // i is column on image (aligned with xReduced) + // j is column on xscan + for (int j = 0; j < xscan->numRows; j++) + { + xscan->data.F32[j][i] -= xReduced->data.F32[i]; + } + } + } + psFree(xReduced); + psFree(yReduced); + } + // pmOverscanUpdateHeader (hdu, overscanOpts, chi2); return true; - } - - // we are performing a statitical analysis of the overscan region - - // Check for an unallowable pmFit. - if (overscanOpts->primary) { - if (overscanOpts->primary->fitType != PM_FIT_NONE && - overscanOpts->primary->fitType != PM_FIT_POLY_ORD && - overscanOpts->primary->fitType != PM_FIT_POLY_CHEBY && - overscanOpts->primary->fitType != PM_FIT_SPLINE) { - psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d). Returning original image.\n", - overscanOpts->primary->fitType); - return false; - } - } - if (overscanOpts->secondary) { - if (overscanOpts->secondary->fitType != PM_FIT_NONE && - overscanOpts->secondary->fitType != PM_FIT_POLY_ORD && - overscanOpts->secondary->fitType != PM_FIT_POLY_CHEBY && - overscanOpts->secondary->fitType != PM_FIT_SPLINE) { - psError(PS_ERR_UNKNOWN, true, "Invalid fit type (2D) (%d). Returning original image.\n", - overscanOpts->secondary->fitType); - return false; - } - } - - psList *overscans = input->bias; // List of the overscan images - - // Reduce all overscan pixels to a single value - if (overscanOpts->single) { - - // extract overscan pixels to a single vector - psVector *pixels = psVectorAlloc(0, PS_TYPE_F32); - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - int index = pixels->n; // Index - pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols); - pixels->n += overscan->numRows * overscan->numCols; - for (int i = 0; i < overscan->numRows; i++) { - memcpy(&pixels->data.F32[index], overscan->data.F32[i], - overscan->numCols * sizeof(psF32)); - index += overscan->numCols; - } - } - psFree(iter); - - // statistic to be calculated - psStatsOptions statistic = psStatsSingleOption(overscanOpts->primary->stat->options); // Statistic to use - if (!statistic) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n", - overscanOpts->primary->stat); - return false; - } - psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original - - if (!psVectorStats(stats, pixels, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - psFree(pixels); - double reduced = psStatsGetValue(stats, statistic); // Result of statistics - - psString comment = NULL; // Comment to add - psStringAppend(&comment, "Overscan value: %f", reduced); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - - // write metadata header value - // XXX EAM : this could / should write the stdev of the overscan region - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", reduced); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN); - - psScalar *reducedScalar = psScalarAlloc(reduced, PS_TYPE_F32); - psBinaryOp (image, image, "-", psMemIncrRefCounter(reducedScalar)); // NOTE: psBinaryOp frees arg2 if it a scalar, so we need to bump to re-use - - // subtract the measured value from each overscan region as well - iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - psBinaryOp(overscan, overscan, "-", psMemIncrRefCounter(reducedScalar)); // NOTE: psBinaryOp frees arg2 if it a scalar, so we need to bump to re-use - } - psFree(iter); - psFree(reducedScalar); - - // EAM 2022.03.29 : if the calculated overscan value is below the threshold, - // declare the readout dead and mask - - if ((reduced < overscanOpts->minValid) || (reduced > overscanOpts->maxValid)) { - fprintf (stderr, "bad overscan (1) %f, masking readout\n", reduced); - psImage *mask = input->mask; - for (int y = 0; y < mask->numRows; y++) { - for (int x = 0; x < mask->numCols; x++) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= overscanOpts->maskVal; - } - } - } - - psFree(stats); - return true; - } - - bool mdok = false; - - // We are performing a row-by-row overscan subtraction - int cellreaddir = psMetadataLookupS32(&mdok, input->parent->concepts, "CELL.READDIR"); // Read direction - if ((cellreaddir != 1) && (cellreaddir != 2)) { - psError(PS_ERR_UNKNOWN, true, "CELL.READDIR must be 1 (rows) or 2 (cols)\n"); - return false; - } - - float chi2 = NAN; // chi^2 from fit - - // adjust operation depending on the read direction : need to re-org pixels for columns - if (!overscanOpts->TwoD && (cellreaddir == 1)) { - // The read direction is rows - psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels - for (int i = 0; i < pixels->n; i++) { - pixels->data[i] = psVectorAlloc(0, PS_TYPE_F32); - } - - // Pull the pixels out into the vectors - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - // the overscan and image might not be aligned. pixels->data represents - // the image row pixels. - int diff = overscan->row0 - image->row0; // Offset between the two regions - for (int i = PS_MAX(0,diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++) { - int j = i - diff; - // i is row on image - // j is row on overscan - psVector *values = pixels->data[i]; - int index = values->n; // Index in the vector - values = psVectorRealloc(values, values->n + overscan->numCols); - values->n += overscan->numCols; - // XXX double-check the range of values->n here - memcpy(&values->data.F32[index], overscan->data.F32[j], - overscan->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); - index += overscan->numCols; - pixels->data[i] = values; // Update the pointer in case it's moved - } - } - psFree(iter); - - // Reduce the overscans - psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels); - psFree(pixels); - if (! reduced) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n"); - return false; - } - - if (!pmOverscanUpdateHeaderVector (input, hdu, overscanOpts, reduced)) { - psError(PS_ERR_UNKNOWN, false, "failure to update header"); - return false; - } - - // Subtract row by row - for (int i = 0; i < image->numRows; i++) { - for (int j = 0; j < image->numCols; j++) { - image->data.F32[i][j] -= reduced->data.F32[i]; - } - } - - // subtract from the overscan regions - { - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - // the overscan and image might not be aligned. - int diff = overscan->row0 - image->row0; // Offset between the two regions - for (int i = PS_MAX(0,diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++) { - int j = i - diff; - // i is row on image - // j is row on overscan - for (int k = 0; k < overscan->numCols; k++) { - overscan->data.F32[j][k] -= reduced->data.F32[j]; - } - } - } - psFree(iter); - } - psFree(reduced); - } - - if (!overscanOpts->TwoD && (cellreaddir == 2)) { - // The read direction is columns - psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels - for (int i = 0; i < pixels->n; i++) { - psVector *values = psVectorAlloc(0, PS_TYPE_F32); - pixels->data[i] = values; - } - - // Pull the pixels out into the vectors - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - // the overscan and image might not be aligned. pixels->data represents - // the image row pixels. - int diff = overscan->col0 - image->col0; // Offset between the two regions - for (int i = PS_MAX(0,diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++) { - int iFixed = i - diff; - // i is column on image - // iFixed is column on overscan - psVector *values = pixels->data[i]; - int index = values->n; // Index in the vector - values = psVectorRealloc(values, values->n + overscan->numRows); - for (int j = 0; j < overscan->numRows; j++) { - values->data.F32[index++] = overscan->data.F32[j][iFixed]; - } - values->n += overscan->numRows; - pixels->data[i] = values; // Update the pointer in case it's moved - } - } - psFree(iter); - - // Reduce the overscans - psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels); - psFree(pixels); - if (! reduced) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n"); - return false; - } - - if (!pmOverscanUpdateHeaderVector (input, hdu, overscanOpts, reduced)) { - psError(PS_ERR_UNKNOWN, false, "failure to update header"); - return false; - } - - // Subtract column by column - for (int j = 0; j < image->numRows; j++) { - for (int i = 0; i < image->numCols; i++) { - image->data.F32[j][i] -= reduced->data.F32[i]; - } - } - - // subtract from the overscan regions - { - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - // the overscan and image might not be aligned. - int diff = overscan->col0 - image->col0; // Offset between the two regions - for (int i = PS_MAX(0,diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++) { - int j = i - diff; - // i is col on image - // j is col on overscan - for (int k = 0; i < overscan->numRows; k++) { - overscan->data.F32[k][j] -= reduced->data.F32[j]; - } - } - } - psFree(iter); - } - - psFree(reduced); - } - - // 2D bias subtraction with x-dir readout direction: the - // bias is constructed by combining a 1D pattern in the - // readout direction from the top overscan region and a second - // 1D pattern in the cross direction from the overscan - if (overscanOpts->TwoD && (cellreaddir == 1)) { - psAssert (overscanOpts->secondary, "2D overscan subtraction requires OVERSCAN.2D parameters"); - // we require 2 overscan regions, and they must match these directions: - if (overscans->n != 2) { - psLogMsg (__func__, PS_LOG_ERROR, "OVERSCAN.2D requires 2 overscan regions but %d supplied", (int) overscans->n); - psLogMsg (__func__, PS_LOG_ERROR, "e.g.: CELL.BIASSEC STR [591:624,1:608],[1:590,598:608]"); - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to generate overscan vector.\n"); - return false; - } - - // the serial (fast readout) direction is columns (x-dir) - psImage *yscan = psListGet (overscans, 0); // overscan region spanning all rows - psImage *xscan = psListGet (overscans, 1); // overscan region spanning all columns - - // Extract the y-dir overscan vector. The overscan and image might not be aligned: - // diff represents the offset between the rows in the image data and the overscan. - // pixels->data represents the image row pixels. For example, the image region may be - // inset in the y-direction but the overscan could cover the entire y-range - - // The read direction is rows - psArray *yscanPixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels - for (int i = 0; i < yscanPixels->n; i++) { - yscanPixels->data[i] = psVectorAlloc(0, PS_TYPE_F32); - } - - // XXX this code allows multiple yscans to be appended, but this does not - // match the concept of how they are assigned above: biassec[0] = yscan - int yDiff = yscan->row0 - image->row0; // Offset between the two regions - for (int i = PS_MAX(0,yDiff); i < PS_MIN(image->numRows, yscan->numRows + yDiff); i++) { - int j = i - yDiff; - // i is row on image, j is row on yscan - psVector *values = yscanPixels->data[i]; - int index = values->n; // Index in the vector - values = psVectorRealloc(values, values->n + yscan->numCols); - values->n += yscan->numCols; - // XXX double-check the range of values->n here - memcpy(&values->data.F32[index], yscan->data.F32[j], - yscan->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); - index += yscan->numCols; - yscanPixels->data[i] = values; // Update the pointer in case it's moved - } - - // Extract the x-dir overscan vector. The overscan and image might not be aligned: - // diff represents the offset between the rows in the image data and the overscan. - // pixels->data represents the image row pixels. For example, the image region may be - // inset in the x-direction but the overscan could cover the entire y-range - - // Extract the top region as a vector of the columns - psArray *xscanPixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels - for (int i = 0; i < xscanPixels->n; i++) { - xscanPixels->data[i] = psVectorAlloc(0, PS_TYPE_F32); - } - - int xDiff = xscan->col0 - image->col0; // Offset between the two regions - for (int ix = PS_MAX(0,xDiff); ix < PS_MIN(image->numCols, xscan->numCols + xDiff); ix++) { - int jx = ix - xDiff; - // ix is row on image, jx is column on xscan - psVector *values = xscanPixels->data[ix]; - values = psVectorRealloc(values, xscan->numRows); - values->n = xscan->numRows; - for (int iy = 0; iy < xscan->numRows; iy++) { - values->data.F32[iy] = xscan->data.F32[iy][jx]; - } - xscanPixels->data[ix] = values; // Update the pointer in case it's moved - } - - // Reduce the overscans - // XXX need to save 2 different chi-square values - psVector *yReduced = pmOverscanVector(&chi2, overscanOpts->primary, yscanPixels); - psFree(yscanPixels); - if (!yReduced) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate y-dir overscan vector.\n"); - return false; - } - if (!pmOverscanUpdateHeaderVector (input, hdu, overscanOpts, yReduced)) { - psError(PS_ERR_UNKNOWN, false, "failure to update header"); - return false; - } - - // Reduce the overscans - // XXX need to save 2 different chi-square values - psVector *xReduced = pmOverscanVector(&chi2, overscanOpts->secondary, xscanPixels); - psFree(xscanPixels); - if (!xReduced) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate x-dir overscan vector.\n"); - return false; - } - - // XXX apply the 2D bias correction here - for (int i = 0; i < image->numRows; i++) { - for (int j = 0; j < image->numCols; j++) { - image->data.F32[i][j] -= yReduced->data.F32[i]; - } - } - - // subtract the y-dir vector from the y-dir overscan regions (why?) - { - // the overscan and image might not be aligned. - int diff = yscan->row0 - image->row0; // Offset between the two regions - for (int i = PS_MAX(0,diff); i < PS_MIN(image->numRows, yscan->numRows + diff); i++) { - int j = i - diff; - // i is row on image (aligned with yReduced) - // j is row on yscan - for (int k = 0; k < yscan->numCols; k++) { - yscan->data.F32[j][k] -= yReduced->data.F32[i]; - } - } - } - - // subtract the x-dir vector from the x-dir overscan regions (why?) - { - // the overscan and image might not be aligned. - int diff = xscan->col0 - image->col0; // Offset between the two regions - for (int i = PS_MAX(0,diff); i < PS_MIN(image->numCols, xscan->numCols + diff); i++) { - int j = i - diff; - // i is column on image (aligned with xReduced) - // j is column on xscan - for (int k = 0; k < xscan->numRows; k++) { - xscan->data.F32[k][j] -= xReduced->data.F32[i]; - } - } - } - psFree(xReduced); - psFree(yReduced); - } - // pmOverscanUpdateHeader (hdu, overscanOpts, chi2); - return true; } // End of overscan subtraction @@ -437,13 +500,13 @@ static void pmOverscanOptionsFree(pmOverscanOptions *options) { - psFree(options->primary); - psFree(options->secondary); + psFree(options->primary); + psFree(options->secondary); } static void pmOverscanStatOptionsFree(pmOverscanStatOptions *options) { - psFree(options->stat); - psFree(options->poly); - psFree(options->spline); + psFree(options->stat); + psFree(options->poly); + psFree(options->spline); } @@ -451,21 +514,21 @@ pmOverscanStatOptions *pmOverscanStatOptionsAlloc(void) { - pmOverscanStatOptions *opts = psAlloc(sizeof(pmOverscanStatOptions)); - psMemSetDeallocator(opts, (psFreeFunc)pmOverscanStatOptionsFree); - - // Inputs - opts->fitType = PM_FIT_NONE; - opts->order = 0; - opts->stat = NULL; - - // Smoothing - opts->boxcar = 0; - opts->gauss = 0.0; - - // Outputs - opts->poly = NULL; - opts->spline = NULL; - - return opts; + pmOverscanStatOptions *opts = psAlloc(sizeof(pmOverscanStatOptions)); + psMemSetDeallocator(opts, (psFreeFunc)pmOverscanStatOptionsFree); + + // Inputs + opts->fitType = PM_FIT_NONE; + opts->order = 0; + opts->stat = NULL; + + // Smoothing + opts->boxcar = 0; + opts->gauss = 0.0; + + // Outputs + opts->poly = NULL; + opts->spline = NULL; + + return opts; } @@ -473,147 +536,166 @@ pmOverscanOptions *pmOverscanOptionsAlloc(void) { - pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions)); - psMemSetDeallocator(opts, (psFreeFunc)pmOverscanOptionsFree); - - // Inputs - opts->single = false; - opts->constant = false; - opts->TwoD = false; - - opts->value = 0.0; - - opts->minValid = 0.0; // default value if not defined - opts->maxValid = (float) 0x10000; // default value if not defined - opts->maskVal = 0x0001; // default value if not defined - - // stat options - opts->primary = NULL; - opts->secondary = NULL; - - return opts; + pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions)); + psMemSetDeallocator(opts, (psFreeFunc)pmOverscanOptionsFree); + + // Inputs + opts->single = false; + opts->constant = false; + opts->TwoD = false; + + opts->value = 0.0; + + opts->minValid = 0.0; // default value if not defined + opts->maxValid = (float)0x10000; // default value if not defined + opts->maskVal = 0x0001; // default value if not defined + + // stat options + opts->primary = NULL; + opts->secondary = NULL; + + return opts; } // Produce an overscan vector from an array of pixels -psVector *pmOverscanVector(float *chi2, // chi^2 from fit - pmOverscanStatOptions *overscanOpts, // Overscan statistic options - const psArray *pixels // Array of vectors containing the pixel values - ) +psVector *pmOverscanVector(float *chi2, // chi^2 from fit + pmOverscanStatOptions *overscanOpts, // Overscan statistic options + const psArray *pixels // Array of vectors containing the pixel values +) { - assert(overscanOpts); - assert(pixels); - - // statisctic to be calculated - psStatsOptions statistic = psStatsSingleOption(overscanOpts->stat->options); // Statistic to use - if (!statistic) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n", overscanOpts->stat); - return false; - } - psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original - - // Reduce the overscans - psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row - psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate - psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_VECTOR_MASK); // Mask for fitting - - for (int i = 0; i < pixels->n; i++) { - psVector *values = pixels->data[i]; // Vector with overscan values - if (values->n > 0) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0; - ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1] - if (!psVectorStats(stats, values, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + assert(overscanOpts); + assert(pixels); + + // statisctic to be calculated + psStatsOptions statistic = psStatsSingleOption(overscanOpts->stat->options); // Statistic to use + if (!statistic) + { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n", overscanOpts->stat); + return false; + } + psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original + + // Reduce the overscans + psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row + psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate + psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_VECTOR_MASK); // Mask for fitting + + for (int i = 0; i < pixels->n; i++) + { + psVector *values = pixels->data[i]; // Vector with overscan values + if (values->n > 0) + { + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0; + ordinate->data.F32[i] = 2.0 * (float)i / (float)pixels->n - 1.0; // Scale to [-1,1] + if (!psVectorStats(stats, values, NULL, NULL, 0)) + { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + goto escape; + } + reduced->data.F32[i] = psStatsGetValue(stats, statistic); + } + else + { + if (overscanOpts->fitType == PM_FIT_NONE) + { + psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the image, and no fit is requested.\n"); + goto escape; + } + else + { + // We'll fit this one out + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1; + } + } + } + // Smooth the reduced vector + if (overscanOpts->boxcar > 0) + { + psVector *smoothed = psVectorBoxcar(NULL, reduced, overscanOpts->boxcar); // Smoothed vector + psFree(reduced); + reduced = smoothed; + } + if (isfinite(overscanOpts->gauss) && overscanOpts->gauss > 0) + { + if (overscanOpts->boxcar > 0) + { + psWarning("Gaussian smoothing the boxcar smoothed overscan --- you asked for it."); + } + psVector *smoothed = psVectorSmooth(NULL, reduced, overscanOpts->gauss, SMOOTH_NSIGMA); + psFree(reduced); + reduced = smoothed; + } + + // Fit the overscan, if required + psVector *fitted = NULL; // Fitted overscan values + switch (overscanOpts->fitType) + { + case PM_FIT_NONE: + // No fitting --- that's easy. + fitted = psMemIncrRefCounter(reduced); + break; + case PM_FIT_POLY_ORD: + if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || + overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) + { + psFree(overscanOpts->poly); + overscanOpts->poly = NULL; + } + if (!overscanOpts->poly) + { + overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order); + } + psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); + fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); + break; + case PM_FIT_POLY_CHEBY: + if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || + overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) + { + psFree(overscanOpts->poly); + overscanOpts->poly = NULL; + } + if (!overscanOpts->poly) + { + overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order); + } + psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); + fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); + break; + case PM_FIT_SPLINE: + + // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and it assumes + // a knot for every input point. it needs an argument like 'number of knots' for the + // output spline. EAM: still true 2023.01.22 + + // overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate); + // fitted = psSpline1DEvalVector(overscanOpts->spline, ordinate); + psError(PS_ERR_UNKNOWN, true, "Spline overscan fitting is broken\n"); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType); goto escape; - } - reduced->data.F32[i] = psStatsGetValue(stats, statistic); - } else { - if (overscanOpts->fitType == PM_FIT_NONE) { - psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the image, and no fit is requested.\n"); - goto escape; - } else { - // We'll fit this one out - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1; - } - } - } - // Smooth the reduced vector - if (overscanOpts->boxcar > 0) { - psVector *smoothed = psVectorBoxcar(NULL, reduced, overscanOpts->boxcar); // Smoothed vector + } + + if (chi2) + { + *chi2 = 0.0; // chi^2 (sort of) + for (int i = 0; i < reduced->n; i++) + { + *chi2 += PS_SQR(fitted->data.F32[i] - reduced->data.F32[i]); + } + } + psFree(reduced); - reduced = smoothed; - } - if (isfinite(overscanOpts->gauss) && overscanOpts->gauss > 0) { - if (overscanOpts->boxcar > 0) { - psWarning("Gaussian smoothing the boxcar smoothed overscan --- you asked for it."); - } - psVector *smoothed = psVectorSmooth(NULL, reduced, overscanOpts->gauss, SMOOTH_NSIGMA); + psFree(ordinate); + psFree(mask); + psFree(stats); + return fitted; + +escape: psFree(reduced); - reduced = smoothed; - } - - // Fit the overscan, if required - psVector *fitted = NULL; // Fitted overscan values - switch (overscanOpts->fitType) { - case PM_FIT_NONE: - // No fitting --- that's easy. - fitted = psMemIncrRefCounter(reduced); - break; - case PM_FIT_POLY_ORD: - if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || - overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) { - psFree(overscanOpts->poly); - overscanOpts->poly = NULL; - } - if (! overscanOpts->poly) { - overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order); - } - psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); - fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); - break; - case PM_FIT_POLY_CHEBY: - if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || - overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) { - psFree(overscanOpts->poly); - overscanOpts->poly = NULL; - } - if (! overscanOpts->poly) { - overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order); - } - psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); - fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); - break; - case PM_FIT_SPLINE: - - // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and it assumes - // a knot for every input point. it needs an argument like 'number of knots' for the - // output spline. EAM: still true 2023.01.22 - - // overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate); - // fitted = psSpline1DEvalVector(overscanOpts->spline, ordinate); - psError(PS_ERR_UNKNOWN, true, "Spline overscan fitting is broken\n"); - break; - default: - psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType); - goto escape; - } - - if (chi2) { - *chi2 = 0.0; // chi^2 (sort of) - for (int i = 0; i < reduced->n; i++) { - *chi2 += PS_SQR(fitted->data.F32[i] - reduced->data.F32[i]); - } - } - - psFree(reduced); - psFree(ordinate); - psFree(mask); - psFree(stats); - return fitted; - -escape: - psFree(reduced); - psFree(ordinate); - psFree(mask); - psFree(stats); - return NULL; + psFree(ordinate); + psFree(mask); + psFree(stats); + return NULL; } @@ -621,88 +703,97 @@ // the reduced overscan vector with the 0-order element from the polynomial fit. Not clear is // adding any useful information -bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanStatOptions *overscanOpts, float chi2) { - - psString comment = NULL; // Comment to add - - switch (overscanOpts->fitType) { - case PM_FIT_POLY_ORD: - case PM_FIT_POLY_CHEBY: { - psStringAppend(&comment, "Overscan fit (chi2: %.2f): ", chi2); - psPolynomial1D *poly = overscanOpts->poly; // The polynomial - for (int i = 0; i < poly->nX; i++) { - psStringAppend(&comment, "%.1f ", poly->coeff[i]); - } - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - comment = NULL; - - // write metadata header value - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", poly->coeff[0]); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", poly->coeffErr[0]); - break; - } - case PM_FIT_SPLINE: { - /* - psSpline1D *spline = overscanOpts->spline; // The spline - for (int i = 0; i < spline->n; i++) { - psStringAppend(&comment, "Overscan fit (chi2: %.2f) %d:", chi2, i); - psPolynomial1D *poly = spline->spline[i]; // i-th polynomial - for (int j = 0; j < poly->nX; j++) { - psStringAppend(&comment, "%.1f ", poly->coeff[i]); - } - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, - comment, ""); - psFree(comment); - comment = NULL; - } - */ - // write metadata header value - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, - "Overscan value", NAN); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, - "Overscan stdev", NAN); - break; - } - case PM_FIT_NONE: - break; - default: - psAbort("Should never get here!!!\n"); - } - return true; +bool pmOverscanUpdateHeader(pmHDU *hdu, pmOverscanStatOptions *overscanOpts, float chi2) +{ + + psString comment = NULL; // Comment to add + + switch (overscanOpts->fitType) + { + case PM_FIT_POLY_ORD: + case PM_FIT_POLY_CHEBY: + { + psStringAppend(&comment, "Overscan fit (chi2: %.2f): ", chi2); + psPolynomial1D *poly = overscanOpts->poly; // The polynomial + for (int i = 0; i < poly->nX; i++) + { + psStringAppend(&comment, "%.1f ", poly->coeff[i]); + } + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + comment = NULL; + + // write metadata header value + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", poly->coeff[0]); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", poly->coeffErr[0]); + break; + } + case PM_FIT_SPLINE: + { + /* + psSpline1D *spline = overscanOpts->spline; // The spline + for (int i = 0; i < spline->n; i++) { + psStringAppend(&comment, "Overscan fit (chi2: %.2f) %d:", chi2, i); + psPolynomial1D *poly = spline->spline[i]; // i-th polynomial + for (int j = 0; j < poly->nX; j++) { + psStringAppend(&comment, "%.1f ", poly->coeff[i]); + } + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, + comment, ""); + psFree(comment); + comment = NULL; + } + */ + // write metadata header value + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, + "Overscan value", NAN); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, + "Overscan stdev", NAN); + break; + } + case PM_FIT_NONE: + break; + default: + psAbort("Should never get here!!!\n"); + } + return true; } // generate stats of overscan vector for header // reduced: 1D vector with overscan stats -bool pmOverscanUpdateHeaderVector (pmReadout *input, pmHDU *hdu, pmOverscanOptions *overscanOpts, psVector *reduced) { - - psString comment = NULL; // Comment to add - psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); - if (!psVectorStats (vectorStats, reduced, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - - // write metadata header value - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan mean", vectorStats->sampleMean); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", vectorStats->sampleStdev); - - // EAM 2022.03.29 : if the calculated overscan value is below the threshold, - // declare the readout dead and mask - - if ((vectorStats->sampleMean < overscanOpts->minValid) || (vectorStats->sampleMean > overscanOpts->maxValid)) { - fprintf (stderr, "bad overscan (2) %f, masking readout\n", vectorStats->sampleMean); - psImage *mask = input->mask; - for (int y = 0; y < mask->numRows; y++) { - for (int x = 0; x < mask->numCols; x++) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= overscanOpts->maskVal; - } - } - } - - psFree (vectorStats); - return true; +bool pmOverscanUpdateHeaderVector(pmReadout *input, pmHDU *hdu, pmOverscanOptions *overscanOpts, psVector *reduced) +{ + + psString comment = NULL; // Comment to add + psStats *vectorStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + if (!psVectorStats(vectorStats, reduced, NULL, NULL, 0)) + { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } + psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + + // write metadata header value + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan mean", vectorStats->sampleMean); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", vectorStats->sampleStdev); + + // EAM 2022.03.29 : if the calculated overscan value is below the threshold, + // declare the readout dead and mask + + if ((vectorStats->sampleMean < overscanOpts->minValid) || (vectorStats->sampleMean > overscanOpts->maxValid)) + { + fprintf(stderr, "bad overscan (2) %f, masking readout\n", vectorStats->sampleMean); + psImage *mask = input->mask; + for (int y = 0; y < mask->numRows; y++) + { + for (int x = 0; x < mask->numCols; x++) + { + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= overscanOpts->maskVal; + } + } + } + + psFree(vectorStats); + return true; } - Index: /branches/2dbias/psModules/src/detrend/pmOverscan.h =================================================================== --- /branches/2dbias/psModules/src/detrend/pmOverscan.h (revision 42678) +++ /branches/2dbias/psModules/src/detrend/pmOverscan.h (revision 42679) @@ -81,5 +81,5 @@ // bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanOptions *overscanOpts, float chi2); -bool pmOverscanSubtract (pmReadout *input, pmOverscanOptions *overscanOpts); +bool pmOverscanSubtract (pmReadout *input, pmOverscanOptions *overscanOpts, bool doTwoDOverscan); /// @}