From 3d23c95ab71b0ea77b3073a0a84c4295939be823 Mon Sep 17 00:00:00 2001 From: mdlpstsci Date: Wed, 19 Apr 2023 17:59:09 -0400 Subject: [PATCH] Basic infrastructure changes needed to support saturation map. (#563) * Basic infrastructure changes needed to support saturation map. * Implementation of the use of a full image to flag full-well saturated pixels. The new calibration reference file is associated with the SATUFILE FITS keyword. Full-welll saturated pixels are no longer flagged in the DoDQI step. The new function, doFullWellSat(), which controls the flagging is executed after the BLEVCORR and BIASCORR steps have been done. * Added parentheses to decision equations to make the operator precedence explicit. * Replace some hardcoded values in dofwsat.c with defined variables, and update the date of changes to correspond to a time close to the actual release date. * Forget to update the date in wf3version.h. * Very minor changes to a comment and added a space for pound include. --- pkg/wfc3/calwf3/Dates | 7 + pkg/wfc3/calwf3/History | 7 + pkg/wfc3/calwf3/Updates | 7 + pkg/wfc3/calwf3/include/wf3dq.h | 2 +- pkg/wfc3/calwf3/include/wf3info.h | 6 +- pkg/wfc3/calwf3/include/wf3version.h | 4 +- pkg/wfc3/calwf3/lib/dodqi.c | 28 +- pkg/wfc3/calwf3/lib/getccdtab.c | 19 +- pkg/wfc3/calwf3/lib/wf3info.c | 7 +- pkg/wfc3/calwf3/lib/wscript | 2 +- pkg/wfc3/calwf3/wf3ccd/doccd.c | 46 ++- pkg/wfc3/calwf3/wf3ccd/doccd.h | 2 + pkg/wfc3/calwf3/wf3ccd/dofwsat.c | 433 +++++++++++++++++++++++++++ pkg/wfc3/calwf3/wf3ccd/getflags.c | 50 +++- 14 files changed, 571 insertions(+), 49 deletions(-) create mode 100644 pkg/wfc3/calwf3/wf3ccd/dofwsat.c diff --git a/pkg/wfc3/calwf3/Dates b/pkg/wfc3/calwf3/Dates index 475111f28..ab9038d0f 100644 --- a/pkg/wfc3/calwf3/Dates +++ b/pkg/wfc3/calwf3/Dates @@ -1,3 +1,10 @@ +17-Apr-2023 - MDD - Version 3.7.0 + - Implementation to use an image to detect and flag full-well saturation versus a simple scalar. + The new routine, dofwsat.c, is similar to what has been done for calacs. The WFC3 implementation + is more complicated in that there are serial virtual overscan columns, as well as binned images, + to accommodate. The detection/flagging occurs after blev and bias correction while the output is + still in counts. + 27-May-2021 - MDD - Version 3.6.2 - Bug fix to address calwf3.e crashing (Abort trap: 6) when taking an existing *_ima.fits (IR) file (produced with FLATCORR=PERFORM and CRCORR=OMIT) with CRCORR set to PERFORM. This reentrant diff --git a/pkg/wfc3/calwf3/History b/pkg/wfc3/calwf3/History index ddb68fed7..15b04ab39 100644 --- a/pkg/wfc3/calwf3/History +++ b/pkg/wfc3/calwf3/History @@ -1,3 +1,10 @@ +### 17-Apr-2023 - MDD - Version 3.7.0 + - Implementation to use an image to detect and flag full-well saturation versus a simple scalar. + The new routine, dofwsat.c, is similar to what has been done for calacs. The WFC3 implementation + is more complicated in that there are serial virtual overscan columns, as well as binned images, + to accommodate. The detection/flagging occurs after blev and bias correction while the output is + still in counts. + ### 27-May-2021 - MDD - Version 3.6.2 - Bug fix to address calwf3.e crashing (Abort trap: 6) when taking an existing *_ima.fits (IR) file (produced with FLATCORR=PERFORM and CRCORR=OMIT) with CRCORR set to PERFORM. This reentrant diff --git a/pkg/wfc3/calwf3/Updates b/pkg/wfc3/calwf3/Updates index ba57f8615..646f45842 100644 --- a/pkg/wfc3/calwf3/Updates +++ b/pkg/wfc3/calwf3/Updates @@ -1,3 +1,10 @@ +Updates for Version 3.7.0 17-Apr-2023 - MDD + - Implementation to use an image to detect and flag full-well saturation versus a simple scalar. + The new routine, dofwsat.c, is similar to what has been done for calacs. The WFC3 implementation + is more complicated in that there are serial virtual overscan columns, as well as binned images, + to accommodate. The detection/flagging occurs after blev and bias correction while the output is + still in counts. + Updates for Version 3.6.2 27-May-2021 - MDD - Bug fix to address calwf3.e crashing (Abort trap: 6) when taking an existing *_ima.fits (IR) file (produced with FLATCORR=PERFORM and CRCORR=OMIT) with CRCORR set to PERFORM. This reentrant diff --git a/pkg/wfc3/calwf3/include/wf3dq.h b/pkg/wfc3/calwf3/include/wf3dq.h index 478132f46..b8d72909c 100644 --- a/pkg/wfc3/calwf3/include/wf3dq.h +++ b/pkg/wfc3/calwf3/include/wf3dq.h @@ -17,7 +17,7 @@ # define UNSTABLE 32 /* IR unstable pixel */ # define WARMPIX 64 /* warm pixel */ # define BADBIAS 128 /* bad bias value */ -# define SATPIXEL 256 /* full-well or a-to-d saturated pixel */ +# define SATPIXEL 256 /* full-well saturated pixel */ # define BADFLAT 512 /* bad flatfield value */ # define TRAP 1024 /* UVIS charge trap, SINK pixel */ # define SPIKE 1024 /* CR spike detected during cridcalc IR */ diff --git a/pkg/wfc3/calwf3/include/wf3info.h b/pkg/wfc3/calwf3/include/wf3info.h index 382b01289..1592f2c45 100644 --- a/pkg/wfc3/calwf3/include/wf3info.h +++ b/pkg/wfc3/calwf3/include/wf3info.h @@ -86,6 +86,10 @@ M. De La Pena 2020 March: Read the PCTERNOI value from the raw image header for possible use in the CTE reduction. + + M. De La Pena 2022 February: + Added the satmap variable which rendered the "saturate" variable + obsolete. Removed "saturate" as a cleanup operation. */ @@ -153,7 +157,6 @@ typedef struct { float mean_gain; /* mean actual gain of all amps */ int ampx; /* first column affected by amps on 2/4amp readout*/ int ampy; /* first row affected by amps on 2/4amp readout*/ - float saturate; /* CCD saturation level */ int trimx[4]; /* Width of overscan to trim off ends of each line */ int trimy[2]; /* Amount of overscan to trim off ends of each col */ int vx[4]; @@ -227,6 +230,7 @@ typedef struct { RefImage nlin; /* non-linearity coefficients image */ RefImage zoff; /* super zero */ RefTab pctetab; /*uvis CTE parameter table*/ + RefImage satmap; /* full-well saturation image */ } WF3Info; diff --git a/pkg/wfc3/calwf3/include/wf3version.h b/pkg/wfc3/calwf3/include/wf3version.h index f9da1db0d..7896267eb 100644 --- a/pkg/wfc3/calwf3/include/wf3version.h +++ b/pkg/wfc3/calwf3/include/wf3version.h @@ -4,7 +4,7 @@ /* This string is written to the output primary header as CAL_VER. */ -# define WF3_CAL_VER "3.6.2 (May-27-2021)" -# define WF3_CAL_VER_NUM "3.6.2" +# define WF3_CAL_VER "3.7.0 (Apr-17-2023)" +# define WF3_CAL_VER_NUM "3.7.0" #endif /* INCL_WF3VERSION_H */ diff --git a/pkg/wfc3/calwf3/lib/dodqi.c b/pkg/wfc3/calwf3/lib/dodqi.c index 2db82ceb6..3e438c302 100644 --- a/pkg/wfc3/calwf3/lib/dodqi.c +++ b/pkg/wfc3/calwf3/lib/dodqi.c @@ -121,6 +121,12 @@ static void FirstLast (double *, double *, int *, int *, int *, int *, H.Bushouse, 2007 Dec 21: Updated to use new rewrite of FirstLast routine provided by P. Hodge from calstis. + + M. De La Pena, 2022 February + Removed flagging of full-well saturated pixels based upon a science + pixel value being greater than a defined scalar value. Use of a + new full-well saturation image supersedes the functionality previously + done in this routine. */ int doDQI (WF3Info *wf3, SingleGroup *x) { @@ -154,7 +160,6 @@ SingleGroup *x io: image to be calibrated; DQ array written to in-place int i, j, i0, j0; /* indexes for scratch array ydq */ int m, n; /* indexes for data quality array in x */ short sum_dq; /* for binning data quality array */ - float sat; int row; /* loop index for row number */ int dimx, dimy; @@ -178,25 +183,18 @@ SingleGroup *x io: image to be calibrated; DQ array written to in-place return (status); /* For the CCD, check for and flag saturation. */ - sat = wf3->saturate; if (wf3->detector != IR_DETECTOR) { dimx = x->sci.data.nx; dimy = x->sci.data.ny; for (j = 0; j < dimy; j++) { - for (i = 0; i < dimx; i++) { - /* Flag a-to-d saturated pixels with 2048 dq bit */ - if (Pix (x->sci.data, i, j) > ATOD_SATURATE) { - sum_dq = DQPix (x->dq.data, i, j) | ATODSAT; - DQSetPix (x->dq.data, i, j, sum_dq); /* atod sat */ - } - /* Flag full-well or atod saturated pixels with 256 bit */ - if (Pix (x->sci.data, i, j) > sat || - Pix (x->sci.data, i, j) > ATOD_SATURATE) { - sum_dq = DQPix (x->dq.data, i, j) | SATPIXEL; - DQSetPix (x->dq.data, i, j, sum_dq); /* saturated */ - } - } + for (i = 0; i < dimx; i++) { + /* Flag a-to-d saturated pixels with 2048 dq bit */ + if (Pix (x->sci.data, i, j) > ATOD_SATURATE) { + sum_dq = DQPix (x->dq.data, i, j) | ATODSAT; + DQSetPix (x->dq.data, i, j, sum_dq); /* atod sat */ + } + } } } diff --git a/pkg/wfc3/calwf3/lib/getccdtab.c b/pkg/wfc3/calwf3/lib/getccdtab.c index 5582dbc75..47d0a7f2b 100644 --- a/pkg/wfc3/calwf3/lib/getccdtab.c +++ b/pkg/wfc3/calwf3/lib/getccdtab.c @@ -23,7 +23,6 @@ typedef struct { IRAFPointer cp_readnoise[NAMPS]; IRAFPointer cp_ampx; IRAFPointer cp_ampy; - IRAFPointer cp_saturate; IRAFPointer cp_pedigree; IRAFPointer cp_descrip; int nrows; /* number of rows in table */ @@ -40,7 +39,6 @@ typedef struct { float readnoise[NAMPS]; int ampx; int ampy; - float saturate; } TblRow; static int OpenCCDTab (char *, TblInfo *); @@ -112,6 +110,10 @@ static int CloseCCDTab (TblInfo *); 21 Oct 2009: H.Bushouse: Added computation of mean_gain to GetCCDTab. mean_gain is now used in flatcorr steps for doing gain conversion. + + 16 Feb 2022: M. De La Pena: + The "saturate" variable became obsolete once the full-well saturation + map was implemented. Removed "saturate" as a cleanup operation. */ int GetCCDTab (WF3Info *wf3, int dimx, int dimy) { @@ -188,15 +190,14 @@ int dimy i: number of lines in exposure } wf3->mean_gain /= NAMPS; - /* For WFC3/UVIS exposures, correct ampx to match the actual - size of the image for more seamless processing of subarrays. - Leave ampx as is for WFC3/IR exposures. */ + /* For WFC3/UVIS exposures, correct ampx to match the actual + size of the image for more seamless processing of subarrays. + Leave ampx as is for WFC3/IR exposures. */ if (wf3->detector == CCD_DETECTOR) wf3->ampx = (tabrow.ampx > dimx) ? dimx : tabrow.ampx; else wf3->ampx = tabrow.ampx; wf3->ampy = tabrow.ampy; - wf3->saturate = tabrow.saturate; break; } @@ -276,7 +277,6 @@ static int OpenCCDTab (char *tname, TblInfo *tabinfo) { c_tbcfnd1 (tabinfo->tp, "READNSED", &tabinfo->cp_readnoise[3]); c_tbcfnd1 (tabinfo->tp, "AMPX", &tabinfo->cp_ampx); c_tbcfnd1 (tabinfo->tp, "AMPY", &tabinfo->cp_ampy); - c_tbcfnd1 (tabinfo->tp, "SATURATE", &tabinfo->cp_saturate); /* Initialize counters here... */ missing = 0; @@ -308,7 +308,6 @@ static int OpenCCDTab (char *tname, TblInfo *tabinfo) { if (tabinfo->cp_readnoise[3] == 0 ) { missing++; nocol[i] = YES;} i++; if (tabinfo->cp_ampx == 0 ) { missing++; nocol[i] = YES;} i++; if (tabinfo->cp_ampy == 0 ) { missing++; nocol[i] = YES;} i++; - if (tabinfo->cp_saturate == 0) { missing++; nocol[i] = YES;} i++; if (PrintMissingCols (missing, NUMCOLS, nocol, colnames, "CCDTAB", tabinfo->tp) ) return(status); @@ -437,10 +436,6 @@ static int ReadCCDTab (TblInfo *tabinfo, int row, TblRow *tabrow) { if (c_iraferr()) return (status = TABLE_ERROR); - c_tbegtr (tabinfo->tp, tabinfo->cp_saturate, row, &tabrow->saturate); - if (c_iraferr()) - return (status = TABLE_ERROR); - return (status); } diff --git a/pkg/wfc3/calwf3/lib/wf3info.c b/pkg/wfc3/calwf3/lib/wf3info.c index 0f0837ca8..90e8edb99 100644 --- a/pkg/wfc3/calwf3/lib/wf3info.c +++ b/pkg/wfc3/calwf3/lib/wf3info.c @@ -71,6 +71,11 @@ M. De La Pena 2020 March: Added reading of PCTERNOI keyword from the raw primary header for possible use in the CTE reduction + + M. De La Pena 2022 February: + Added variable, satmap - the reference image for full-well + saturation. Use of this image rendered wf3->saturate variable + obsolete. Removed "wf3->saturate" as part of the cleanup operation. */ void WF3Init (WF3Info *wf3) { @@ -128,7 +133,6 @@ void WF3Init (WF3Info *wf3) { } wf3->ampx = 0; wf3->ampy = 0; - wf3->saturate = 0.; wf3->trimx[0] = 0; wf3->trimx[1] = 0; wf3->trimx[2] = 0; @@ -185,6 +189,7 @@ void WF3Init (WF3Info *wf3) { InitRefTab (&(wf3->oscn)); InitRefTab (&(wf3->atod)); InitRefTab (&(wf3->pctetab)); + InitRefImg (&(wf3->satmap)); /* Initialize reference images and tables for WF32D */ InitRefImg (&(wf3->dark)); diff --git a/pkg/wfc3/calwf3/lib/wscript b/pkg/wfc3/calwf3/lib/wscript index 45ab329d9..3243ea667 100644 --- a/pkg/wfc3/calwf3/lib/wscript +++ b/pkg/wfc3/calwf3/lib/wscript @@ -30,7 +30,7 @@ def build(bld): ../wf3ccd/doatod.c ../wf3ccd/dobias.c ../wf3ccd/doblev.c ../wf3ccd/doccd.c ../wf3ccd/doflash.c ../wf3ccd/findblev.c ../wf3ccd/findover.c ../wf3ccd/getflags.c - ../wf3ccd/getccdsw.c ../wf3ccd/sink.c + ../wf3ccd/getccdsw.c ../wf3ccd/sink.c ../wf3ccd/dofwsat.c ../wf3ir/blevcorr.c ../wf3ir/cridcalc.c ../wf3ir/darkcorr.c ../wf3ir/doir.c ../wf3ir/dqicorr.c ../wf3ir/flatcorr.c diff --git a/pkg/wfc3/calwf3/wf3ccd/doccd.c b/pkg/wfc3/calwf3/wf3ccd/doccd.c index 45816b771..22ac42ecb 100644 --- a/pkg/wfc3/calwf3/wf3ccd/doccd.c +++ b/pkg/wfc3/calwf3/wf3ccd/doccd.c @@ -55,11 +55,15 @@ Moved the subarray codes to a general function and changed the calculation of where the subarray is. Also moved the DQ bit to the singlegroup extension + M. De La Pena February 2022 + Modified to apply the full-well saturation flags stored as an image + to the data in doFullWellSat() instead of in the doDQI step. + */ # include # include -#include "hstcal.h" +# include "hstcal.h" # include "hstio.h" # include "wf3.h" # include "wf3info.h" @@ -235,7 +239,6 @@ int DoCCD (WF3Info *wf3, int extver) { if (dqiHistory (wf3, x.globalhdr)) return (status); - /* ANALOG TO DIGITAL CORRECTION. */ AtoDMsg (wf3, extver); if (wf3->atodcorr == PERFORM) { @@ -313,6 +316,27 @@ int DoCCD (WF3Info *wf3, int extver) { if (biasHistory (wf3, x.globalhdr)) return (status); + /* Apply the saturation image. + Strictly speaking, the application of the full-well saturation image is + not a calibration step (i.e., there is no SATCORR), but the application + of a 2D image to flag pixels versus using a single scalar to flag + saturated pixels as previously done in DQICORR will be done in doFullWellSat() + after BLEVCORR and BIASCORR. This correction should only be done if both + BLEVCORR and BIASCORR have been performed. This flagging is only applicable + for the UVIS. */ + + if (wf3->biascorr == PERFORM && wf3->blevcorr == PERFORM) { + SatMsg (wf3, extver); + sprintf(MsgText, "\nFull-well saturation flagging being performed."); + trlmessage(MsgText); + if (doFullWellSat(wf3, &x)) { + return (status); + } + } else { + sprintf(MsgText, "\nNo Full-well saturation flagging being performed.\n"); + trlwarn(MsgText); + } + /*UPDATE THE SINK PIXELS IN THE DQ MASK OF BOTH SCIENCE IMAGE SETS IT'S DONE HERE WITH ONE CALL TO THE FILE BECAUSE THEY NEED TO BE PROCESSED IN THE RAZ FORMAT JAY USES, THOUGH ONLY ONE CHIP DONE HERE @@ -340,8 +364,8 @@ int DoCCD (WF3Info *wf3, int extver) { */ /* Make temporary full group - use the sinkfile as the reference - sub2full will reset the pixel values + use the sinkfile as the reference + sub2full will reset the pixel values */ initSingleGroup(&fullarray); allocSingleGroup(&fullarray,RAZ_COLS/2,RAZ_ROWS, True); @@ -551,6 +575,20 @@ static void BiasMsg (WF3Info *wf3, int extver) { } } +static void SatMsg (WF3Info *wf3, int extver) { + + int OmitStep (int); + void PrSwitch (char *, int); + void PrRefInfo (char *, char *, char *, char *, char *); + + trlmessage (""); + if (extver == 1 && !OmitStep (wf3->biascorr)) { + + PrRefInfo ("satufile", wf3->satmap.name, wf3->satmap.pedigree, + wf3->satmap.descrip, ""); + } +} + static void FlashMsg (WF3Info *wf3, int extver) { int OmitStep (int); diff --git a/pkg/wfc3/calwf3/wf3ccd/doccd.h b/pkg/wfc3/calwf3/wf3ccd/doccd.h index c07ea8cdd..f2538a38d 100644 --- a/pkg/wfc3/calwf3/wf3ccd/doccd.h +++ b/pkg/wfc3/calwf3/wf3ccd/doccd.h @@ -4,6 +4,7 @@ /*USEFUL LIB FUNCTIONS*/ static void AtoDMsg (WF3Info *, int); static void BiasMsg (WF3Info *, int); +static void SatMsg (WF3Info *, int); static void FlashMsg (WF3Info *, int); static void BlevMsg (WF3Info *, int); static void dqiMsg (WF3Info *, int); @@ -13,6 +14,7 @@ int GetCorner (Hdr *, int , int *, int *); int doAtoD (WF3Info *, SingleGroup *); int atodHistory (WF3Info *, Hdr *); int doBias (WF3Info *, SingleGroup *); +int doFullWellSat(WF3Info *, SingleGroup *); int biasHistory (WF3Info *, Hdr *); int doFlash (WF3Info *, SingleGroup *, float *); int flashHistory (WF3Info *, Hdr *); diff --git a/pkg/wfc3/calwf3/wf3ccd/dofwsat.c b/pkg/wfc3/calwf3/wf3ccd/dofwsat.c new file mode 100644 index 000000000..6fa6413af --- /dev/null +++ b/pkg/wfc3/calwf3/wf3ccd/dofwsat.c @@ -0,0 +1,433 @@ +# include +# include + +# include "hstcal.h" +# include "hstio.h" +# include "wf3.h" +# include "wf3dq.h" +# include "wf3info.h" +# include "hstcalerr.h" + +/* + This routine compares the full-well saturation image pixels against + those of the SingleGroup x science array pixels. If the science array + pixel value >= the value in the saturation image, the corresponding + dq array pixel of the SingleGroup x is ORed with the full-well saturation + flag, SATPIXEL, and the dq array is modified in-place. + Note the saturation image data is only valid in the region of a + trimmed full-frame image. In order to apply the saturation image to the + science data (full-frame or subarray), the starting X and Y pixels, as well + as the X and Y sizes need to be obtained, and the saturation image properly + overlaid on the science data. Note that formerly the full-well saturation + flags were applied during doDQI. + + Michele De La Pena, 2022 March 14 + Initial implementation of the full-well saturation flagging done + via an image. Based on a similar routine written for ACS, as well as + WFC3 doflash to accommodate binned or subarray data. + + Michele De La Pena, 2023 April 17 + Replace hardcoded values with variables isolated to this module. + + */ + +/* SIZE_SV_OVERSCAN is the size of the serial virtual overscan region + (in pixels) between amps on the same chip. + END_PIX_AC_AMP is the last pixel value (0-based system) of amp A + or C (CCD area plus serial phyical overscan */ +# define SIZE_SV_OVERSCAN 60 +# define END_PIX_AC_AMP 2072 + +int doFullWellSat(WF3Info *wf3, SingleGroup *x) { + + /* + WF3Info *wf3 i: calibration switches, etc. + SingleGroup *x io: image to be modified; dq array is updated in-place + */ + + extern int status; + + SingleGroupLine y; /* scratch space */ + SingleGroup satimage; /* storage for entire saturation image */ + int extver = 1; /* get this imset from bias image */ + int rx, ry; /* for binning bias image down to size of x */ + int x0, y0; /* offsets of sci image */ + int same_size; /* true if no binning of ref image required */ + int xdim; /* number of columns in science image */ + int ydim; /* number of lines in science image */ + short sum_dq; /* total of the DQ bits for each pixel */ + int is_subarray = 0; /* identification of data as a subarray */ + int straddle = 0; /* subarray starts in A or C and straddles the virtual overscan in the reference image */ + int overstart = -1; /* location where the overscan starts in the cut science image */ + int xbeg[2], ybeg[2]; + int xend[2], yend[2]; + + int FindLine (SingleGroup *, SingleGroupLine *, int *, int *, int *, int *, int *); + int DetCCDChip (char *, int, int, int *); + void computeLimits(WF3Info *, int, int, int *, int *, int *, int *); + + {unsigned int i; + for (i = 0; i < 2; i++) { + xbeg[i] = -1; + xend[i] = -1; + ybeg[i] = -1; + yend[i] = -1; + }} + + /* Initialize local variables */ + rx = 1; + ry = 1; + x0 = 0; + y0 = 0; + + /* + Compute correct extension version number to extract from + reference image to correspond to CHIP in science data. + Note: wf3->nimsets is no longer used in DetCCDChip() routine. + */ + if (DetCCDChip (wf3->satmap.name, wf3->chip, wf3->nimsets, &extver)) + return (status); + + /* Get the first line of saturation image data */ + initSingleGroupLine (&y); + openSingleGroupLine (wf3->satmap.name, extver, &y); + if (hstio_err()) + return (status = OPEN_FAILED); + + /* + Reference image should already be selected to have the + same binning factor as the science image. All we need to + make sure of is whether the science array is a subarray. + + x0,y0 is the location of the start of the + subimage in the reference image. + */ + if (FindLine (x, &y, &same_size, &rx, &ry, &x0, &y0)) + return (status); + sprintf(MsgText,"Image has starting location of %d,%d in the reference image", x0, y0); + trlmessage(MsgText); + + /* Clean up the SingleGroupLine object */ + closeSingleGroupLine (&y); + freeSingleGroupLine (&y); + + /* If reference data not binned same as input, return with an error */ + if (rx != 1 || ry != 1) { + sprintf (MsgText, + "Saturation image and input are not binned to the same pixel size!"); + trlerror (MsgText); + return (status = SIZE_MISMATCH); + } + + /* If the science and reference images are not the same size, check to + see if we are calibrating a subarray image against the untrimmed, + full-frame, saturation image. The full-frame saturation image will contain + serial virtual overscan columns in the middle of it, which must be accounted + for in the value of x0 if the science image subarray is located to the + right (higher column numbers) of the virtual overscan. This will be the + case for science image subarrays located in the amp B and D quadrants. + + For subarrays which START in B or D we can just move the x0 over 60 pixels + to avoid the serial virtual overscan in the reference image, which is always + an untrimmed, full-frame image. + + Otherwise, only part of the subarray overlaps the serial virtual overscan + and special measures must be taken to avoid it. + */ + + /* Science image dimensions */ + xdim = x->sci.data.nx; + ydim = x->sci.data.ny; + + /* Full-frame */ + if (same_size) { + sprintf (MsgText,"Saturation image and input are the same size."); + trlmessage (MsgText); + /* Subarray */ + } else { + is_subarray = 1; + sprintf (MsgText,"Saturation image and input are NOT the same size - SUBARRAY found, amp %s", wf3->ccdamp); + trlmessage (MsgText); + + /* + ONLY 1 AMP is used to read subarrays, so ampx and ampy should be + set to the size of the image for all cases + */ + wf3->ampx = xdim; + wf3->ampy = ydim; + + /* The image starts in the B or D regions so we can just shift the starting pixel */ + if (x0 > END_PIX_AC_AMP) { + if (wf3->verbose) { + sprintf (MsgText,"Subarray starts in B or D region, moved from (%d,%d) to ", x0, y0); + trlmessage (MsgText); + } + x0 += SIZE_SV_OVERSCAN; + if (wf3->verbose) { + sprintf (MsgText,"(%d,%d) to avoid virtual overscan in reference image", x0, y0); + trlmessage (MsgText); + } + + /* The subarray starts somewhere in A or C and might straddle the virtual overscan region */ + } else { + if ((x0 + xdim) > END_PIX_AC_AMP) { + straddle = 1; + overstart = (END_PIX_AC_AMP + 1) - x0; + } + } + } + + if (wf3->verbose) { + sprintf(MsgText,"ccdamp = %s, straddle = %d, (xdim,ydim) = (%d,%d), (ampx,ampy) = (%d,%d), (x0,y0) = (%d,%d)", + wf3->ccdamp, straddle, xdim, ydim, wf3->ampx, wf3->ampy, x0, y0); + trlmessage(MsgText); + } + + /* Get the full saturation image */ + initSingleGroup(&satimage); + getSingleGroup(wf3->satmap.name, extver, &satimage); + if (hstio_err()) { + freeSingleGroup(&satimage); + return (status = OPEN_FAILED); + } + + /* If the science image is binned, it will be assumed to have + same size as the reference image, since reading subarrays + of a binned chip is not supported in current flight software. + + Need to ignore the columns which represent the serial virtual + overscan between Amps A and B and Amps C and D as there is no + saturation information for this data. The values in the + saturation image are set to zero in this location. + Bin size: 1x1 ==> overscan = 60 pixels, pixels (2073 - 2132) + Bin size: 2x2 ==> overscan = 28 pixels, pixels (1037 - 1064) + Bin size: 3x3 ==> overscan = 20 pixels, pixels (691 - 710) + The pixel ranges documented here (begin - end) are zero-based. + */ + + /* Determine the limits of real data to process - note the end point + in the loops is "<" and NOT "<=". + */ + computeLimits(wf3, xdim, ydim, xbeg, ybeg, xend, yend); + if (wf3->verbose) { + sprintf(MsgText,"xbeg[0]: %d xend[0]: %d xbeg[1]: %d xend[1]: %d", xbeg[0], xend[0], xbeg[1], xend[1]); + trlmessage(MsgText); + sprintf(MsgText,"ybeg[0]: %d yend[0]: %d ybeg[1]: %d yend[1]: %d", ybeg[0], yend[0], ybeg[1], yend[1]); + trlmessage(MsgText); + } + + /* The saturation image already has the gain applied, but the science data + at this stage in the WFC3 pipeline is still in counts. It is + necessary to divide out the gain from the saturation data before + any comparison is done. + */ + + if (wf3->verbose) { + sprintf(MsgText, "Mean gain: %f", wf3->mean_gain); + trlmessage(MsgText); + } + + /* Full-frame */ + if (same_size) { + + /* Loop over the lines in the science image */ + {unsigned int j; + for (j=ybeg[0]; j < yend[0]; j++) { + + /* Loop over the indices in the line in the science image */ + {unsigned int i; + for (i = xbeg[0]; i < xend[0]; i++) { + /* Flag full-well saturated pixels with 256 dq bit*/ + if (Pix(x->sci.data, i, j) > (Pix(satimage.sci.data, i, j) / wf3->mean_gain)) { + sum_dq = DQPix(x->dq.data, i, j) | SATPIXEL; + DQSetPix(x->dq.data, i, j, sum_dq); + } + }} + + /* If there is a second Amp in play, complete the processing of the line */ + {unsigned int i; + for (i = xbeg[1]; i < xend[1]; i++) { + /* Flag full-well saturated pixels with 256 dq bit*/ + if (Pix(x->sci.data, i, j) > (Pix(satimage.sci.data, i, j) / wf3->mean_gain)) { + sum_dq = DQPix(x->dq.data, i, j) | SATPIXEL; + DQSetPix(x->dq.data, i, j, sum_dq); + } + }} + }} + sprintf(MsgText, "Full-frame full-well saturation image flagging step done."); + trlmessage(MsgText); + + /* Subarray */ + } else { + + /* Loop over all the indices and lines in the science array, and + match the data to the appropriate indices and lines in the reference image. + + i - index in science image + k - index in reference image + j - line in the science image + l - line in reference image + */ + + {unsigned int j, l; + for (j = 0, l = y0; j < ydim; j++, l++) { + + /* Working with a subarray so need to apply the proper + section from the reference image to the science image. + */ + + {unsigned int i, k; + for (i = 0, k = x0; i < xdim; i++, k++) { + + /* Increase the value of l to jump over the virtual overscan */ + if (i == overstart) + k += SIZE_SV_OVERSCAN; + + /* Flag full-well saturated pixels with 256 dq bit*/ + if (Pix(x->sci.data, i, j) > (Pix(satimage.sci.data, k, l) / wf3->mean_gain)) { + sum_dq = DQPix(x->dq.data, i, j) | SATPIXEL; + DQSetPix(x->dq.data, i, j, sum_dq); + } + }} + }} + sprintf(MsgText, "Subarray full-well saturation image flagging step done.\n"); + trlmessage(MsgText); + } + + freeSingleGroup(&satimage); + + return (status); +} + + +/* Because the pixels representing the serial virtual overscan are present + in the saturation image, it is necessary to skip over these pixels when + determining the flagging for full-well saturation. + + This routine is based upon code in doblev.c. +*/ +void +computeLimits(WF3Info *wf3, int xdim, int ydim, int begx[2], int begy[2], int endx[2], int endy[2]) { + int amp; /* Counter for amps used */ + int numamps; /* Number of AMPS used to readout chip */ + char *ccdamp; /* Amps which are used for this chip */ + char *amploc; + int trimx1, trimx2; /* width to trim off ends of each line */ + int trimx3, trimx4; /* width to trim off middle of each line */ + int trimy1, trimy2; /* amount to trim off ends of each column */ + int bias_loc, amp_indx; + int bias_ampx, bias_ampy; + int biassect[4]; /* section(s) to use */ + int bias_orderx[4] = {0,1,0,1}; + int bias_ordery[4] = {0,0,1,1}; + + void parseWFCamps (char *, int, char *); + + /* Copy out overscan size information for ease of reference in this function */ + trimx1 = wf3->trimx[0]; /* Left serial physical */ + trimx2 = wf3->trimx[1]; /* Right serial physical */ + trimx3 = wf3->trimx[2]; /* Left serial virtual */ + trimx4 = wf3->trimx[3]; /* Right serial virtual */ + trimy1 = wf3->trimy[0]; /* Bottom (Chip 1) parallel virtual */ + trimy2 = wf3->trimy[1]; /* Top (Chip 2) parallel virtual */ + + /* Establish which amps from ccdamp string are appropriate + ** for this chip */ + ccdamp = (char *) calloc (NAMPS+1, sizeof(char)); + ccdamp[0] = '\0'; + + /* Set up the 2 AMPS used per line */ + if (wf3->detector == CCD_DETECTOR) { + parseWFCamps (wf3->ccdamp, wf3->chip, ccdamp); + } + sprintf(MsgText,"Amp: %s chip: %d ccdamp: %s ", wf3->ccdamp, wf3->chip, ccdamp); + trlmessage(MsgText); + + /* Set up the 2 AMPS used per line */ + + /* How many amps are used for this chip */ + numamps = strlen(ccdamp); + sprintf(MsgText,"Number of amps %d", numamps); + trlmessage(MsgText); + + for (amp = 0; amp < numamps; amp++) { + bias_loc = 0; + + /* determine which section goes with the amp */ + /* bias_amp = 0 for BIASSECTA, = 1 for BIASSECTB */ + amploc = strchr (AMPSORDER, ccdamp[amp]); + bias_loc = *amploc - *ccdamp; + amp_indx = *amploc - *AMPSORDER; + bias_ampx = bias_orderx[bias_loc]; + bias_ampy = bias_ordery[bias_loc]; + + /* Requirement: at least 1 section should be specified! */ + /* If both bias sections are specified,... */ + if (wf3->biassectc[1] > 0 && wf3->biassectd[1] > 0) { + + /* select section nearest the amp based on bias_amp */ + biassect[0] = (bias_ampx == 0) ? wf3->biassectc[0] : + wf3->biassectd[0]; + biassect[1] = (bias_ampx == 0) ? wf3->biassectc[1] : + wf3->biassectd[1]; + + /* If only 1 AMP is used for each line, use both sections */ + if (numamps == 1) { + biassect[2] = (bias_ampx == 0) ? wf3->biassectd[0] : + wf3->biassectc[0]; + biassect[3] = (bias_ampx == 0) ? wf3->biassectd[1] : + wf3->biassectc[1]; + } else { + biassect[2] = 0; + biassect[3] = 0; + } + + /* if neither of the virtual overscan regions are specified + ** use the physical overscan instead */ + } else if (wf3->biassectc[1] <= 0 && wf3->biassectd[1] <= 0) { + + biassect[0] = (wf3->biassecta[1] == 0) ? wf3->biassectb[0] : + wf3->biassecta[0]; + biassect[1] = (wf3->biassecta[1] == 0) ? wf3->biassectb[1] : + wf3->biassecta[1]; + /* then set trailing section to zero indicating it's not used*/ + biassect[2] = 0; + biassect[3] = 0; + + } else { + + /* otherwise, select the non-zero bias section */ + biassect[0] = (wf3->biassectc[1] == 0) ? wf3->biassectd[0] : + wf3->biassectc[0]; + biassect[1] = (wf3->biassectc[1] == 0) ? wf3->biassectd[1] : + wf3->biassectc[1]; + /* then set trailing section to zero indicating it's not used*/ + biassect[2] = 0; + biassect[3] = 0; + } + if (wf3->verbose) { + sprintf (MsgText, "Using overscan columns %d to %d", + biassect[0]+1, biassect[1]+1); + trlmessage (MsgText); + if (biassect[2] != 0) { + sprintf (MsgText, " and columns %d to %d", + biassect[2]+1, biassect[3]+1); + trlmessage (MsgText); } + } + + /* Compute range of pixels affected by each amp */ + begx[amp] = trimx1 + (wf3->ampx + trimx3 + trimx4) * bias_ampx; + endx[amp] = (bias_ampx == 0 && wf3->ampx != 0) ? wf3->ampx + trimx1 : + xdim - trimx2; + begy[amp] = trimy1 + wf3->ampy* bias_ampy; + endy[amp] = (bias_ampy == 0 && wf3->ampy != 0) ? wf3->ampy +trimy1 : + ydim - trimy2; + + /* Make sure that endx and endy do not extend beyond the bounds of + ** the image... HAB 7 May 2001 (same as calacs change). */ + if (endx[amp] > xdim) endx[amp] = xdim; + if (endy[amp] > ydim) endy[amp] = ydim; + } + + free(ccdamp); +} diff --git a/pkg/wfc3/calwf3/wf3ccd/getflags.c b/pkg/wfc3/calwf3/wf3ccd/getflags.c index 4e8842266..0c876874e 100644 --- a/pkg/wfc3/calwf3/wf3ccd/getflags.c +++ b/pkg/wfc3/calwf3/wf3ccd/getflags.c @@ -34,6 +34,9 @@ int GetImageRef (RefFileInfo *, Hdr *, char *, RefImage *, int *); for each reference file, as well as verifying correct selection criteria such as DETECTOR, FILTER, and CCDGAIN. + M. De La Pena, 2022 February + Added new SATUFILE: Full-well saturation image. + */ int GetFlags (WF3Info *wf3, Hdr *phdr) { @@ -148,6 +151,8 @@ int *nsteps io: incremented if this step can be performed extern int status; + int saveBiasCorr = GOOD_PEDIGREE; + int calswitch; int GetSwitch (Hdr *, char *, int *); void MissingFile (char *, char *, int *); @@ -164,24 +169,45 @@ int *nsteps io: incremented if this step can be performed return (status); } - if (GetImageRef (wf3->refnames, phdr, "BIASFILE", &wf3->bias, - &wf3->biascorr)) - return (status); - if (wf3->bias.exists != EXISTS_YES) { - MissingFile ("BIASFILE", wf3->bias.name, missing); + if (GetImageRef (wf3->refnames, phdr, "BIASFILE", &wf3->bias, &wf3->biascorr)) + return (status); - } else { + if (wf3->bias.exists != EXISTS_YES) { + MissingFile ("BIASFILE", wf3->bias.name, missing); + } else { - /* Is the FILETYPE appropriate for a BIAS file? */ - CheckImgType (&wf3->bias, "BIAS", "BIASFILE", missing); + /* Is the FILETYPE appropriate for a BIAS file? */ + CheckImgType (&wf3->bias, "BIAS", "BIASFILE", missing); - /* Does it have the correct GAIN value? */ - if (CheckGain(wf3->bias.name, wf3->ccdgain, "CCDGAIN", missing)) - return (status); - } + /* Does it have the correct GAIN value? */ + if (CheckGain(wf3->bias.name, wf3->ccdgain, "CCDGAIN", missing)) + return (status); + } if (wf3->biascorr == PERFORM) (*nsteps)++; + + /* Save the value for recovery */ + saveBiasCorr = wf3->biascorr; + + /* + Also check for the new full-well saturation image which is + applied after BLEVCORR and BIASCORR are done. Since the reference + file is not associated with its own "calibration step keyword" + (e.g., SATUCORR), just using the BIASCORR key as a standin here - + make sure the BIASCORR retains its value as set in the above code. + + This is a kludge. + */ + if (GetImageRef (wf3->refnames, phdr, + "SATUFILE", &wf3->satmap, &wf3->biascorr)) + return (status); + + /* Recover the biascorr setting */ + wf3->biascorr = saveBiasCorr; + + if (wf3->satmap.exists != EXISTS_YES) + MissingFile ("SATUFILE", wf3->satmap.name, missing); } return (status);