Implement new (via image) and maintain old (scalar) saturation image …

…flagging (#585)

* Updates to not only implement the new saturation image, but to aslo to
be backwards compatible when the SATFILE FITS keyword is missing or
contains a bad value in the PHDU of the RAW file.

* Added a comment to doccd.c to indicate the original flagging of
full-well saturated data is done if the SATFILE keyword is missing
or has a bad value.

* Fixed a typo

* Changed spaces to tabs, lines 193-199. */

pkg/wfc3/calwf3/Dates

@@ -1,9 +1,12 @@
-17-Apr-2023 - MDD - Version 3.7.0
+08-June-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.
+       to accommodate.  The detection/flagging occurs after blev and bias corrections while the output is
+       still in counts.  If the SATUFILE keyword is missing from the FITS header, or the keyword does
+       not have a valid filename as a value, the code will revert to using the original method of flagging
+       full-well saturation.  The flagging will be done in doDQI and use a single value as the saturation 
+       threshold.
 
 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

pkg/wfc3/calwf3/History

@@ -1,9 +1,12 @@
-### 17-Apr-2023 - MDD - Version 3.7.0
+### 08-Jun-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.
+       to accommodate.  The detection/flagging occurs after blev and bias corrections while the output is
+       still in counts.  If the SATUFILE keyword is missing from the FITS header, or the keyword does
+       not have a valid filename as a value, the code will revert to using the original method of flagging
+       full-well saturation.  The flagging will be done in doDQI and use a single value as the saturation 
+       threshold.
 
 ### 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

pkg/wfc3/calwf3/Updates

@@ -1,9 +1,12 @@
-Updates for Version 3.7.0 17-Apr-2023 - MDD
+Updates for Version 3.7.0 08-Jun-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.
+       to accommodate.  The detection/flagging occurs after blev and bias corrections while the output is
+       still in counts.  If the SATUFILE keyword is missing from the FITS header, or the keyword does
+       not have a valid filename as a value, the code will revert to using the original method of flagging
+       full-well saturation.  The flagging will be done in doDQI and use a single value as the saturation 
+       threshold.
 
 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

pkg/wfc3/calwf3/include/wf3info.h

@@ -90,6 +90,14 @@
     M. De La Pena 2022 February:
     Added the satmap variable which rendered the "saturate" variable
     obsolete.  Removed "saturate" as a cleanup operation.
+
+    M. De La Pena 2023 June:
+    Resurrect the "saturate" variable as the full-well flagging can
+    be done using a saturation image, or in the case of a missing
+    SATUFILE FITS keyword, using a single threshold value which was
+    the original algorithm in doDQI. Added a flag, scalar_satflag, to
+    indicate which method should be used for the full-well saturation
+    flagging.
  
 */
 
@@ -157,6 +165,8 @@ 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;     /* full-well saturation threshold as a single value */
+    Bool scalar_satflag;/* indicator whether threshold value vs saturation image should be used */
     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];

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.7.0 (Apr-17-2023)"
+# define WF3_CAL_VER "3.7.0 (Jun-08-2023)"
 # define WF3_CAL_VER_NUM "3.7.0"
 
 #endif /* INCL_WF3VERSION_H */

pkg/wfc3/calwf3/lib/dodqi.c

@@ -127,6 +127,13 @@ static void FirstLast (double *, double *, int *, int *, int *, int *,
     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.
+
+   M. De La Pena, 2023 May
+    Resurrected the ability to flag full-well saturated pixels based upon
+    a science pixel value being greater than a defined scalar value.  This
+    was done at the request of the WFC3 team so that a user can have their
+    processed data flagged even if no saturation image is present. 
+
 */
 
 int doDQI (WF3Info *wf3, SingleGroup *x) {
@@ -160,6 +167,7 @@ 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;			/* saturation threshold */
 
 	int row;		/* loop index for row number */
 	int dimx, dimy;
@@ -182,7 +190,16 @@ SingleGroup *x    io: image to be calibrated; DQ array written to in-place
 	if (wf3->dqicorr != PERFORM && wf3->dqicorr != DUMMY)
 	    return (status);
 
+	/* Issue a message so it is clear that saturation flagging is 
+		being done here in DODQI using scalar *IF* wf3->scalar_satflag is True.
+	*/
+	if (wf3->scalar_satflag == True) {
+		sprintf (MsgText, "Full-well saturation flagging being applied during doDQI using a single threshold value.");
+		trlmessage (MsgText);
+	}
+
 	/* 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;
@@ -194,6 +211,14 @@ SingleGroup *x    io: image to be calibrated; DQ array written to in-place
 					sum_dq = DQPix (x->dq.data, i, j) | ATODSAT;
 					DQSetPix (x->dq.data, i, j, sum_dq); /* atod sat */
 				}
+
+				if (wf3->scalar_satflag == True) {
+					/* Flag full-well or a-to-d 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 */
+					}
+				}
 			}
 	    }
 	}

pkg/wfc3/calwf3/lib/getccdtab.c

@@ -23,6 +23,7 @@ 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 */
@@ -39,6 +40,7 @@ typedef struct {
 	float readnoise[NAMPS];
 	int ampx;
 	int ampy;
+    float saturate;
 } TblRow;
 
 static int OpenCCDTab (char *, TblInfo *);
@@ -114,6 +116,10 @@ static int CloseCCDTab (TblInfo *);
     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.
+
+    25 May 2023: M. De La Pena:
+    Resurrect the "saturate" variable so the original method of flagging
+    saturated pixels can be used if the saturation image is not available.
 */
 
 int GetCCDTab (WF3Info *wf3, int dimx, int dimy) {
@@ -198,6 +204,7 @@ int     dimy      i: number of lines in exposure
 		else
 		    wf3->ampx = tabrow.ampx;
 		wf3->ampy = tabrow.ampy;		
+        wf3->saturate = tabrow.saturate;
 
 		break;
 	    }
@@ -277,6 +284,7 @@ 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,6 +316,7 @@ 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);
@@ -435,6 +444,10 @@ static int ReadCCDTab (TblInfo *tabinfo, int row, TblRow *tabrow) {
 	c_tbegti (tabinfo->tp, tabinfo->cp_ampy, row, &tabrow->ampy);
 	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);
 }

pkg/wfc3/calwf3/lib/getrefname.c

@@ -40,7 +40,7 @@ char *refname     o: directory name and reference file name
 
 	if (!foundit) {
 	    if (GetKeyStr (phdr, keyword, USE_DEFAULT, "", refname, CHAR_LINE_LENGTH))
-		return (status);
+		    return (status);
 	}
 
 	return (status);

pkg/wfc3/calwf3/lib/wf3info.c

@@ -76,6 +76,10 @@
     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.
+
+  M. De La Pena 2023 May:
+    Resurrected the "wf3->saturate" variable so it can be used when a
+    saturation image is not available.
 */
 
 void WF3Init (WF3Info *wf3) {
@@ -133,6 +137,8 @@ void WF3Init (WF3Info *wf3) {
 	}
 	wf3->ampx = 0;
 	wf3->ampy = 0;
+	wf3->saturate = 0.;
+	wf3->scalar_satflag = False;
 	wf3->trimx[0] = 0;
 	wf3->trimx[1] = 0;
 	wf3->trimx[2] = 0;
@@ -295,7 +301,7 @@ int GetTabRef (RefFileInfo *refnames, Hdr *phdr,
 
 void MissingFile (char *keyword, char *filename, int *missing) {
 
-	sprintf (MsgText, "%s `%s' not found or can't open.", keyword, filename);
+	sprintf (MsgText, "%s '%s' not found or cannot be opened.", keyword, filename);
 	trlerror (MsgText);
 	(*missing)++;
 }

pkg/wfc3/calwf3/wf3ccd/doccd.c

@@ -59,6 +59,12 @@
     Modified to apply the full-well saturation flags stored as an image
     to the data in doFullWellSat() instead of in the doDQI step.
 
+    M. De La Pena June 2023
+    Only try to access the SATUFILE keyword if it is actually available in
+    the header.  If the keyword is missing or does not contain a filename,
+    the algorithm will indicate the original method of flagging saturated
+    pixels by using a single value threshold should be used.
+
  */
 
 # include <string.h>
@@ -325,17 +331,14 @@ int DoCCD (WF3Info *wf3, int extver) {
        BLEVCORR and BIASCORR have been performed.  This flagging is only applicable
        for the UVIS. */
 
-    if (wf3->biascorr == PERFORM && wf3->blevcorr == PERFORM) {
+    if (wf3->biascorr == PERFORM && wf3->blevcorr == PERFORM && wf3->scalar_satflag == False) {
         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

pkg/wfc3/calwf3/wf3ccd/dofwsat.c

@@ -254,7 +254,7 @@ int doFullWellSat(WF3Info *wf3, SingleGroup *x) {
                 }
             }}
         }}
-        sprintf(MsgText, "Full-frame full-well saturation image flagging step done.");
+        sprintf(MsgText, "Full-frame full-well saturation image flagging step done.\n");
         trlmessage(MsgText);
 
     /* Subarray */

pkg/wfc3/calwf3/wf3ccd/getflags.c

@@ -37,6 +37,12 @@ int GetImageRef (RefFileInfo *, Hdr *, char *, RefImage *, int *);
    M. De La Pena, 2022 February 
     Added new SATUFILE: Full-well saturation image.
 
+   M. De La Pena, 2023 May
+    Only try to access the SATUFILE keyword if it is actually available in
+    the header.  If the keyword is missing or does not contain a filename,
+    the algorithm will indicate the original method of flagging saturated
+    pixels by using a single value threshold should be used.
+
 */
 
 int GetFlags (WF3Info *wf3, Hdr *phdr) {
@@ -201,13 +207,26 @@ int *nsteps      io: incremented if this step can be performed
        */
         if (GetImageRef (wf3->refnames, phdr,
                          "SATUFILE", &wf3->satmap, &wf3->biascorr))
+        {
+            wf3->scalar_satflag = True;
+	        sprintf (MsgText, "SATUFILE not found or cannot be opened.");
+	        trlerror (MsgText);
+	        sprintf (MsgText, "A single threshold value will be used for full-well saturation flagging.");
+	        trlmessage(MsgText);
             return (status);
+        }
 
         /* Recover the biascorr setting */
         wf3->biascorr = saveBiasCorr;
 
-        if (wf3->satmap.exists != EXISTS_YES)
+        /* Accommodate a missing SATUFILE keyword or associated value */
+        if (wf3->satmap.exists != EXISTS_YES) {
+            wf3->scalar_satflag = True;
             MissingFile ("SATUFILE", wf3->satmap.name, missing);
+            *missing = 0;
+	        sprintf (MsgText, "A single threshold value will be used for full-well saturation flagging.");
+	        trlmessage(MsgText);
+        }
 	}
 
 	return (status);