/* This file contains routines that handle background processing: CalcBack X1DBack Revision history: ---------------- 20 Feb 97 - Implemented (I.Busko) 11 Apr 97 - OPR 33790: flag background values that had more than 30% of their pixels affected by serious DQ flags (IB). 14 Apr 97 - OPR 33792: compute background errors using fit scatter instead of input pixel errors (IB). 15 Apr 97 - OPR 33789: perform sigma-clipping in background fit (IB). 28 Apr 97 - Fixed error in background box coordinate computation (IB). 28 Jul 97 - Print 1-indexed debug coordinates (IB). 10 Apr 98 - Remove debug file (IB). 14 Oct 98 - Add X1D_BAD_BACKGROUND macro (IB) 16 Dec 98 - Output background error: rms around fit / sqrt(n) 08 Sep 99 - Fixed sigma-clip test (IB) 13 Jan 00 - Scattered light correction (IB) */ # include # include # include # include # include # include "../stis.h" # include "../stiserr.h" # include "../stisdq.h" # include "calstis6.h" # define MAX_CLIP 5 /* max. number of sigma-clip iterations */ /* Compute background coefficients. Zero, first and nth order polynomials are supported. The zero-order background is computed as the unweighted average of all eligible pixels in both background extraction boxes. The first order polynomial is fitted to eligible pixels in each background extraction box. The independent variable is the pixel index in the A2 direction (in physical pixel array coordinates, not reference coordinates) and the dependent variable is the pixel content. 7th order is used exclusively by the scattered light correction algorithm. It was implemented later and the code had to be tweaked in order to support it. Note that sigma clipping is effectively turned off when scattered light correction is taking place. The background DQ flag is set according to the total fraction of flagged pixels in the background extraction boxes. */ int CalcBack (StisInfo6 *sts, XtractInfo *xtr, SingleGroup *in, FloatHdrData *ssgx, FloatHdrData *ssgy, int ipix, double ysp, int debug) { /* arguments: StisInfo6 *sts i: calibration switches and info XtractInfo *xtr i: extraction parameters SingleGroup *in i: input image FloatHdrData ssgx; i: small-scale distortion in X (not used) FloatHdrData ssgy; i: small-scale distortion in Y (not used) int ipix; i: index of image pixel in the A1 direction double ysp; i: center of spectrum extraction box int debug; i: debug control */ double rpix; /* index of reference pixel in A1 direction */ double xcent[2]; /* center of background boxes */ double ycent[2]; double bksize[2]; /* sizes of background boxes */ double bkoffset[2]; /* offsets of background boxes */ double npts; /* accumulators used for zero and first */ double sumx; /* order only */ double sumy; double sumx2; double sumy2; double sumxy; double sumvar; double delta; double sigma; /* scatter around current fit */ double new_sigma; /* updated scatter around current fit */ int nbck; /* total number of pixels in backgr. boxes */ int nfbck; /* discarded pixels due to flagging */ int nsbck; /* discarded pixels due to sigma-clipping */ double *yval; /* data for polynomial fit */ double *xval; double *wval; int ndata; /* # of data points in polynomial fit */ double *hold; int backord, iclip, maxclip, i; void X1DBack (StisInfo6 *, SingleGroup *, FloatHdrData *, FloatHdrData *, double, double, double, double *, double *, double *, double *, double *, double *, double *, int *, int *, int *, double, double *, double *, double *, double *, int *, int); int FitPoly (double *, double *, double *, int, int, double *); /* See if there are command-line overrides. Overrides are disabled in case the scatterd light correction algorithm is active. */ if (!(sts->scatter)) { for (i = 0; i < 2; i++) { if (sts->bksize[i] == NO_SIZE) bksize[i] = xtr->bksize[i]; else bksize[i] = sts->bksize[i]; if (sts->bkoffset[i] == NO_SIZE) bkoffset[i] = xtr->bkoffset[i]; else bkoffset[i] = sts->bkoffset[i]; } if (sts->bkord == NO_ORDER) backord = xtr->backord; else backord = sts->bkord; } else { /* Size and offset values are supposed to be re-defined by previous call to DefineBackRegions function. */ for (i = 0; i < 2; i++) { bksize[i] = sts->bksize[i]; bkoffset[i] = sts->bkoffset[i]; } backord = BACKP; } /* Translate image pixel index into reference pixel index. */ rpix = (ipix - sts->ltv[0]) / sts->ltm[0]; /* Compute centers of background boxes. */ xcent[0] = rpix + bkoffset[0] * sts->sin_bktilt; ycent[0] = ysp + bkoffset[0] * sts->cos_bktilt; xcent[1] = rpix + bkoffset[1] * sts->sin_bktilt; ycent[1] = ysp + bkoffset[1] * sts->cos_bktilt; sigma = DBL_MAX / 3.1; /*make sure first 3*sigma computation succeeds*/ sts->ebck = 0.0; sts->bck[0] = 0.0; sts->bck[1] = 0.0;; /* Main sigma-clip loop. */ maxclip = (sts->scatter) ? 1 : MAX_CLIP; for (iclip = 1; iclip <= maxclip; iclip++) { /* Clear all accumulators. */ npts = 0.0; sumx = 0.0; sumy = 0.0; sumx2 = 0.0; sumy2 = 0.0; sumxy = 0.0; sumvar = 0.0; new_sigma = 0.0; nbck = 0; nfbck = 0; nsbck = 0; /* Alloc memory used in 7th order background fit. This must be done even when the scattered light correction algorithm is not active, the reason being to avoid a possible addressing error when the calling sequence to X1DBack gets executed. Note that in this case the sigma-clip loop is executed only once. */ if ((xval = (double *) calloc (200, sizeof (double))) == NULL) return (OUT_OF_MEMORY); if ((yval = (double *) calloc (200, sizeof (double))) == NULL) return (OUT_OF_MEMORY); if ((wval = (double *) calloc (200, sizeof (double))) == NULL) return (OUT_OF_MEMORY); ndata = 0; /* Extract background in each box. */ X1DBack (sts, in, ssgx, ssgy, xcent[0], ycent[0], bksize[0], &npts, &sumx, &sumy, &sumx2, &sumy2, &sumxy, &sumvar, &nbck, &nfbck, &nsbck, sigma, &new_sigma, xval, yval, wval, &ndata, debug); X1DBack (sts, in, ssgx, ssgy, xcent[1], ycent[1], bksize[1], &npts, &sumx, &sumy, &sumx2, &sumy2, &sumxy, &sumvar, &nbck, &nfbck, &nsbck, sigma, &new_sigma, xval, yval, wval, &ndata, debug); /* If more than 30% of the background pixels were rejected by sigma-clip, or flagged, set background's 11th flag bit. */ if ((float)(nfbck+nsbck) / (float)nbck > 0.33333F) sts->dqbck = X1D_BAD_BACKGROUND; else sts->dqbck = 0; /* Compute background coefficients. Only unweighted fits are implemented in this version, so sumvar is never used. */ switch (backord) { case 0: if (npts > 0.0) { sts->bck[0] = sumy / npts; if (npts > 1.0) sts->vbck[0] = (sumy2 - sumy*sumy/npts) / (npts - 1.0); else sts->vbck[0] = 0.0; } else { sts->bck[0] = 0.0; sts->vbck[0] = 0.0; } sts->bck[1] = 0.0; sts->vbck[1] = 0.0; break; case 1: delta = npts * sumx2 - sumx * sumx; if (fabs(delta) > 0.0) { sts->bck[0] = (sumx2 * sumy - sumx * sumxy) / delta; sts->bck[1] = (sumxy * npts - sumx * sumy) / delta; sts->vbck[0] = sumx2 / delta; sts->vbck[1] = npts / delta; } else { sts->bck[0] = 0.0; sts->vbck[0] = 0.0; sts->bck[1] = 0.0; sts->vbck[1] = 0.0; } break; case BACKP: /* This code does NOT compute the coefficient variances ! */ if (ndata > BACKP) { if (!FitPoly (xval, yval, wval, ndata, BACKP, sts->bck)) { for (i = 0; i < BACKP+3; i++) sts->vbck[i] = 0.0; } else { for (i = 0; i < BACKP+3; i++) { sts->bck[i] = 0.0; sts->vbck[i] = 0.0; } } } else { for (i = 0; i < BACKP+3; i++) { sts->bck[i] = 0.0; sts->vbck[i] = 0.0; } } break; default: sts->bck[0] = 0.0; sts->vbck[0] = 0.0; sts->bck[1] = 0.0; sts->vbck[1] = 0.0; break; } /* Free memory used in 7th order background fit. */ free (wval); free (yval); free (xval); /* Update sigma for next iteration. */ if (npts > 1.0) { sigma = sqrt (new_sigma / (npts - 1.0)); sts->ebck = sigma; } else break; /* Exit the loop if no pixels where discarded by sigma-clip. */ /* test, debug (PEH, 1998 Dec 31) ... if (nsbck == 0) break; ... */ if (nsbck == 0 && iclip > 1) break; } return (0); } /* Extract background pixels and update fitting accumulators in a given background extraction box. The interpolation routine does *not* take into account that interpolated pixels may have a tilt respect to the input pixel array. */ void X1DBack (StisInfo6 *sts, SingleGroup *in, FloatHdrData *ssgx, FloatHdrData *ssgy, double xcenter, double ycenter, double size, double *npts, double *sumx, double *sumy, double *sumx2, double *sumy2, double *sumxy, double *sumvar, int *nbck, int *nfbck, int *nsbck, double sigma, double *new_sigma, double *xval, double *yval, double *wval, int *ndata, int debug) { /* arguments: StisInfo6 *sts i: calibration switches and info SingleGroup *in i: input image FloatHdrData ssgx; i: small-scale distortion in X (not used) FloatHdrData ssgy; i: small-scale distortion in Y (not used) double xcenter; i: center of spectrum extraction box in A1 direction double ycenter; i: center of spectrum extraction box in A2 direction double size; i: size of box double *npts... o: accumulators int *nbck; o: total number of accumulated pixels int *nfbck; o: number of flagged pixels int *nsbck; o: number of discarded pixels in this s-clip iteration double sigma; i: scatter around current fit double *new_sigma; io: updated scatter double *yval; io: data for polynomial fit double *xval; double *wval; int ndata; io: # of data points in polynomial fit int debug; i: debug control */ double x1, y1, y2; /* end points of background extraction box */ double xx, yy; /* coordinates of current pixel */ float oSci, oErr; /* interpolated values from image array */ short oDQ; double residual; /* residual background */ double tcent; /* center of trace in image coordinates */ double hold; void Interp2D (SingleGroup *, short, double, double, double, float *, float *, short *); /* Compute endpoints of background extraction box. Notice that the 0.5 pixel offsets below, used to make the endpoint coordinates point to the pixel center, only make sense for small tilt angles. In this case most of the 0.5 pixel correction goes into the A2 direction. A more sophisticated correction is needed for large angles. */ x1 = xcenter - (size / 2.0 * sts->sin_bktilt); y1 = ycenter - (size / 2.0 * sts->cos_bktilt) + 0.5; y2 = ycenter + (size / 2.0 * sts->cos_bktilt) - 0.5; /* Translate reference pixel coordinates into physical indices. */ x1 = x1 * sts->ltm[0] + sts->ltv[0]; y1 = y1 * sts->ltm[1] + sts->ltv[1]; y2 = y2 * sts->ltm[1] + sts->ltv[1]; tcent = ycenter * sts->ltm[1] + sts->ltv[1]; /* Scan background box, extract pixels and update accumulators. The update in A1 direction is clearly by steps of sin(bktilt) size, but in the A2 direction we have the choice of using steps of size 1.0 or cos(bktilt). Which one is better ? */ for (xx = x1, yy = y1; yy <= y2; xx += sts->sin_bktilt, yy += sts->cos_bktilt) { /* Interpolate, checking for out of bounds. */ Interp2D (in, sts->sdqflags, xx, yy, 1.0, &oSci, &oErr, &oDQ); /* Update accumulators, discarding flagged data and outliers. Notice that background is computed as a function of physical pixel coordinates. Two separate sets of tests are used because each background processing algorithm has a different set of requirements. */ (*nbck)++; if (!(sts->scatter)) { /* for zero-th and 1st order */ if (!(oDQ & sts->sdqflags)) { residual = fabs (oSci - (sts->bck[0] + sts->bck[1] * yy)); if (residual <= (3.0 * sigma)) { *npts += 1.0; *sumvar += oErr * oErr; *sumx += yy; *sumy += oSci; *sumx2 += yy * yy; *sumxy += yy * oSci; *sumy2 += oSci * oSci; *new_sigma += residual * residual; } else (*nsbck)++; } else (*nfbck)++; } else { /* for 7th order. Discard data points affected by out-of-image pixels as flagged by Interp2D. */ if (!(oDQ & DETECTORPROB & sts->sdqflags)) { xval[*ndata] = yy; yval[*ndata] = (double)oSci; hold = yy - tcent; if (hold < 0) hold = -hold; if (hold > 0) hold = 1.0 / sqrt (hold); else hold = 1.0; wval[*ndata] = hold; (*ndata)++; } else (*nfbck)++; } } }