# include # include # include # include "pedsub.h" extern float A[11]; extern int MA; # define BADVAL -1.0e30 /* Algorithm for adjusting the DC offset levels that may exist between the four quadrants of a NICMOS image. The basic idea behind the algorithm is to emulate the process used by the human eye in detecting the mismatches between quadrants. The eye is a very good detector of edges, even at very low contrast. Thus the algorithm basically builds an estimate of the high-frequency Fourier power along the edges between quadrants, and uses a minimization code to get the best set of DC offsets that minimize that power. The Fourier power is proportional to the squared differences between intensities extracted at each side of an edge. Intensities are extracted from within a band with given width. The exact position where the edge sits (128) is ignored. For the vertical edge, at a single image line, pixels sitting at the left side of the edge, within the band, are averaged into a1. Pixels sitting at the right side are averaged into a2. The Fourier power contribution from that image line is just (a2 - a1)^2, and contributions from all 128 lines along the vertical edge are added together. The same mechanism is applied to the horizontal edge, resulting in the end in a single scalar quantity that measures the global Fourier power along the edges. By ignoring pixels across most of the image, the power estimator is rather insensitive to the presence of bright sources and background spatial gradients. The code also looks at the Data Quality information for each pixel and can exclude pixels with selected anomalies. The code can also exclude pixels above a certain threshold, and can exclude the largest differences when computing the power estimator. At this point no information from the error array is used. The code uses the "amoeba" (downhill simplex) minimization method to search for the set of DC offsets that minimize the power estimator. The upper-left quadrant is used as a reference, and the three others are varied until the method converges. The resultant image is by definition the one which has the minimum possible Fourier power at the edges, thus is the one most pleasing to the human eye, as far as the edge visual appearance goes. C function that implements the algorithm: ---------------------------------------- The algorithm is activated by simply calling function dcquad. The main parameter is the pixel array itself, which on output gets corrected by the proper offsets. There are some fine-tuning parameters which can for the most part be set to "reasonable" values that result in good fits in any data set the algorithm was tested on. void dcquad (float *sci, short *dq, short mask, int b1, int b2, float rejfrac, float threshold, float *of1, float *of2, float *of3); sci - pointer to the 256X256 buffer that stores the input image data. On output, the buffer will be modified with the DC offsets applied to the three respective quadrants. dq - pointer to the 256X256 buffer that stores the DQ image. mask - used to selectively control pixel rejection with DQ info. To reject all anomalies, set mask to 32767. To include all pixels, set mask to zero. b1, b2 - beginning and ending pixel indexes for band that is sampled at the edge of each quadrant. rejfrac - fraction of pixels/differences to reject. This is used both for rejecting bright pixels when computing the first guesses for amoeba, and for rejecting the largest differences when computing the power estimator. A value of 10-20% seems to be appropriate. threshold - pixels above this are rejected from all computations. This could usually be set to a very large number in order to accept all pixels. of1 - output values of2 of3 External dependencies: --------------------- amoeba.c - Numerical Recipes amoeba function. amotry.c - Used by the above. Revision history: ---------------- I. Busko 20-May-1998 Implemented for NICMOS "dccorr" task. H. Bushouse 25-May-2000 Final revisions for use in pedestal removal tasks. H. Bushouse 21-Jun-2000 Replaced dshellsort routine with shell routine in Num. Rec. library. */ /* This structure stores data used by the amoeba code. */ typedef struct { float *sci; short *dq; short mask; int b1; int b2; float threshold; float rejfrac; int eqorder; } FitData; extern FitData fitdata; FitData fitdata; /* These are fixed amoeba parameters. */ # define NDIM 3 # define DTOL 1.0e-6 /* Macros for accessing 2-D arrays. */ # define PIX(i,j) *(sci+(j)*NIC_XSIZE+(i)) # define DQ(i,j) *(dq +(j)*NIC_XSIZE+(i)) /* AMISMATCH -- This function exists only to enable the amoeba * code to access the extra parameters needed for * the mismatch computation, such as the data * arrays, pixel mask and band limits. It must be * scanned by the compiler *before* the function that * actually calls amoeba. */ float amismatch (float x[]) { float mismatch (float *, short *, short, int, int, int, float, float, float); /* 1-indexed according to Numerical Recipes standard. */ return (mismatch (fitdata.sci, fitdata.dq, fitdata.mask, fitdata.b1, fitdata.b2, fitdata.eqorder, x[1], x[2], x[3])); } /* DCQUAD -- Main algorithm. * */ void dcquad (float *sci, short *dq, short mask, int b1, int b2, float rejfrac, float threshold, int eqorder, float *of1, float *of2, float *of3) { /* Local variables */ float *x,*y,**p; /* amoeba storage */ int nit; /* # of iterations from amoeba */ int i; float hold; /* Function declarations */ float *vector (int nl, int nh); void free_vector (float *v, int nl, int nh); float **matrix (int nrl, int nrh, int ncl, int nch); void free_matrix (float **m, int nrl, int nrh, int ncl, int nch); int initoffs (float *, short *, short, int, int, int, float *, float *, float *, float *); void amoeba (float **p, float y[], int ndim, float dtol, float (*funk)(float []), int *nfunk); /* Decrement band limit indexes to make them zero-indexed */ b1 -= 1; b2 -= 1; /* Store parameters in amoeba data structure. */ fitdata.sci = sci; fitdata.dq = dq; fitdata.mask = mask; fitdata.b1 = b1; fitdata.b2 = b2; fitdata.threshold = threshold; fitdata.rejfrac = rejfrac; fitdata.eqorder = eqorder; /* Starting guesses for the three variable parameters. ** The starting guesses are obtained by getting the clipped ** median value of the offsets between each quadrant and increasing ** them by 50%. This increase is necessary to give amoeba enough ** leverage to properly scan solution space. */ if (initoffs (sci, dq, mask, b1, b2, eqorder, of1, of2, of3, &hold)) { *of1 = *of2 = *of3 = 0.0; return; } *of2 = ((*of1 + *of2) + (*of3 - hold))/2.0; *of1 *= 1.5; *of2 *= 1.5; *of3 *= 1.5; /* Prepare data in format accessible for Numerical Recipes amoeba. ** This uses the vector and matrix 1-indexed data formats defined ** in the NICMOS lib numrec.c library. Simplex is built with steps ** of 10, 20 30 and 40%. */ x = vector (1, NDIM); y = vector (1, NDIM+1); p = matrix (1, NDIM+1, 1, NDIM); for (i = 1; i <= NDIM+1; i++) { x[1] = p[i][1] = *of1 * (1.0 + (float)i/10); x[2] = p[i][2] = *of2 * (1.0 + (float)i/10); x[3] = p[i][3] = *of3 * (1.0 + (float)i/10); y[i] = amismatch (x); } /* Call the amoeba routine to find best set of offset values. */ amoeba (p, y, NDIM, DTOL, amismatch, &nit); /* Get results from amoeba data structure. */ *of1 = p[1][1]; *of2 = p[1][2]; *of3 = p[1][3]; /* Free amoeba resources. */ free_matrix (p, 1, NDIM+1, 1, NDIM); free_vector (y, 1, NDIM+1); free_vector (x, 1, NDIM); } /* MISMATCH -- Compute mismatch between quadrants. ** ** Quadrants are numbered from zero to 3, clockwise, starting with ** the upper left quadrant. */ float mismatch (float *sci, short *dq, short mask, int b1, int b2, int eqorder, float of1, float of2, float of3) { /* Local variables */ int i; /* loop index */ int nb1, nb2; /* difference counters */ double mism; /* final result */ float *buff1, *buff2; /* storage buffers */ /* Function declarations */ int doVert (float *, short *, short, int, int, int, int, int, float *, int *); int doHori (float *, short *, short, int, int, int, int, int, float *, int *); void shell (float *, int); /* Alloc buffers for storing differences. */ buff1 = (float *) malloc ((4*NIC_QSIZE + 10) * sizeof(float)); buff2 = (float *) malloc ((NIC_QSIZE + 10) * sizeof(float)); nb1 = 0; /* Collect differences from vertical band between quadrants 0 and 1 */ if (doVert (sci, dq, mask, b1, b2, eqorder, 2, 3, buff2, &nb2)) return (1); /* Load square of differences into buffer */ for (i = 0; i < nb2; i++) { buff1[nb1] = buff2[i] + of1; buff1[nb1] *= buff1[nb1]; nb1++; } /* Collect differences from vertical band between quadrants 2 and 3 */ if (doVert (sci, dq, mask, b1, b2, eqorder, 0, 1, buff2, &nb2)) return (1); /* Load square of differences into buffer */ for (i = 0; i < nb2; i++) { buff1[nb1] = buff2[i] + of2 - of3; buff1[nb1] *= buff1[nb1]; nb1++; } /* Collect differences from horiz. band between quadrants 0 and 3 */ if (doHori (sci, dq, mask, b1, b2, eqorder, 0, 2, buff2, &nb2)) return (1); /* Load square of differences into buffer */ for (i = 0; i < nb2; i++) { buff1[nb1] = buff2[i] - of3; buff1[nb1] *= buff1[nb1]; nb1++; } /* Collect differences from horiz. band between quadrants 1 and 2 */ if (doHori (sci, dq, mask, b1, b2, eqorder, 1, 3, buff2, &nb2)) return (1); /* Load square of differences into buffer */ for (i = 0; i < nb2; i++) { buff1[nb1] = buff2[i] + of1 - of2; buff1[nb1] *= buff1[nb1]; nb1++; } /* Sort the entire list of differences. */ shell (buff1-1, nb1); /* Discard largest ones. */ nb1 = (int)(nb1 * (1.0 - fitdata.rejfrac/100.)); /* Compute total mismatch. */ for (i = 0, mism = 0.0; i < nb1; mism += buff1[i++]); /* Free work buffers */ free (buff1); free (buff2); /* Return the mismatch value */ return ((float)mism); } /* * DOFIT - Do polynomial fit. * * Fits a polynomial of the requested order to the input vectors * xval vs. yval, containing npts data points. The resulting * fit is evaluated at x pixel coordinate 0.5 and returned. * */ int doFit (float *xval, float *yval, int npts, int order, float *pred) { /* Local variables */ int i; /* loop index */ float chisq; /* chi-square of fit */ float *sval; /* statistical weight vector */ float **u, **v, *w; /* fit work arrays */ /* Function declarations */ void svdfit (float x[], float y[], float sig[], int ndata, float a[], int ma, float **u, float **v, float *w, float *chisq, void (*funcs)(float, float[], int)); void fpoly (float x, float p[], int np); float spoly (float x); int sort (float *arr, int npts); float findMode (float *arr, int npts); float findMedian (float *arr, int npts); float *vector (int nl, int nh); void free_vector (float *v, int nl, int nh); float **matrix (int nrl, int nrh, int ncl, int nch); void free_matrix (float **m, int nrl, int nrh, int ncl, int nch); if (npts <= 0) { *pred = BADVAL; return (0); } /* Allocate memory for fit arrays */ if (order >= 0) { MA = order + 1; u = matrix (1, npts, 1, MA); v = matrix (1, MA, 1, MA); w = vector (1, MA); sval = vector (1, npts); for (i=1; i<=npts; i++) sval[i] = 1.0; } /* Determine the predicted boundary values */ if (order >= 0) { if (npts >= MA+1) { svdfit (xval, yval, sval, npts, A, MA, u, v, w, &chisq, fpoly); *pred = spoly(0.5); } else { *pred = BADVAL; } } else { if (npts >= 1) { if (sort (yval, npts)) return (1); if (order == -1) *pred = findMedian (yval+1, npts); else *pred = findMode (yval+1, npts); } else { *pred = BADVAL; } } /* Free local arrays */ if (order >= 0) { free_matrix (u, 1, npts, 1, MA); free_matrix (v, 1, MA, 1, MA); free_vector (w, 1, MA); free_vector (sval, 1, npts); } return (0); } /* * INITOFFS - Initialize quadrant offset values. * * The initial quadrant offset values are computed by determining the * difference in signal level across each quadrant boundary for each * row or column in each quadrant and returning the clipped median of * the difference values for each quadrant. The difference at each row * or column is computed by collecting the pixel values within a band * on either side of the boundary and fitting a polynomial to those * values or simply computing the median or mode of the values. * */ int initoffs (float *sci, short *dq, short mask, int b1, int b2, int eqorder, float *of1, float *of2, float *of3, float *of4) { /* Local variables */ int i; /* loop index */ int nb; /* difference counter */ float *buffer; /* storage buffer */ /* Function declarations */ int doVert (float *, short *, short, int, int, int, int, int, float *, int *); int doHori (float *, short *, short, int, int, int, int, int, float *, int *); void shell (float *, int); /* Alloc buffers for storing differences */ buffer = (float *) malloc ((NIC_QSIZE + 10) * sizeof(float)); /* Collect differences from vertical band between quadrants 0 and 1 */ if (doVert (sci, dq, mask, b1, b2, eqorder, 2, 3, buffer, &nb)) return (1); /* Flip the sign of the offsets for this quadrant */ for (i = 0; i < nb; i++) buffer[i] = -1.0 * buffer[i]; /* Find the clipped median of the offset values for this quadrant */ if (nb > 0) { shell (buffer-1, nb); nb = (int)(nb * (1.0 - fitdata.rejfrac/100.)); *of1 = buffer[nb/2]; } else { sprintf (MsgText, "No valid pixels for fitting"); n_error (MsgText); return (1); } /* Collect differences from vertical band between quadrants 2 and 3 */ if (doVert (sci, dq, mask, b1, b2, eqorder, 0, 1, buffer, &nb)) return (1); /* Find the clipped median of the offset values for this quadrant */ if (nb > 0) { shell (buffer-1, nb); nb = (int)(nb * (1.0 - fitdata.rejfrac/100.)); *of4 = buffer[nb/2]; } else { sprintf (MsgText, "No valid pixels for fitting"); n_error (MsgText); return (1); } /* Collect differences from horiz. band between quadrants 0 and 3 */ if (doHori (sci, dq, mask, b1, b2, eqorder, 0, 2, buffer, &nb)) return (1); /* Find the clipped median of the offset values for this quadrant */ if (nb > 0) { shell (buffer-1, nb); nb = (int)(nb * (1.0 - fitdata.rejfrac/100.)); *of3 = buffer[nb/2]; } else { sprintf (MsgText, "No valid pixels for fitting"); n_error (MsgText); return (1); } /* Collect differences from horiz. band between quadrants 1 and 2 */ if (doHori (sci, dq, mask, b1, b2, eqorder, 1, 3, buffer, &nb)) return (1); /* Find the clipped median of the offset values for this quadrant */ if (nb > 0) { shell (buffer-1, nb); nb = (int)(nb * (1.0 - fitdata.rejfrac/100.)); *of2 = buffer[nb/2]; } else { sprintf (MsgText, "No valid pixels for fitting"); n_error (MsgText); return (1); } free (buffer); return (0); } /* * DOVERT - Do a vertical boundary. * * Collect quadrant boundary differences for a vertical boundary. * */ int doVert (float *sci, short *dq, short mask, int b1, int b2, int eqorder, int q1, int q2, float *buffer, int *nb) { int i, j; int n1, n2; float fit1, fit2; float *xval, *yval; float *vector (int nl, int nh); void free_vector (float *v, int nl, int nh); int doFit (float *xval, float *yval, int npts, int order, float *pred); /* Allocate memory for the x/y vectors */ xval = vector (1, b2-b1+1); yval = vector (1, b2-b1+1); *nb = 0; /* Loop from bottom to top of vertical boundary */ for (j = QYI[q1]; j <= QYF[q1]; j++) { /* Pixels at edge's left side */ fit1 = 0.0; n1 = 0; for (i = QXF[q1]-b1; i >= QXF[q1]-b2; i--) { if (PIX(i,j) < fitdata.threshold && !(DQ(i,j) & mask)) { n1++; xval[n1] = (float)(QXF[q1]-i+1); yval[n1] = PIX(i,j); } } if (doFit (xval, yval, n1, eqorder, &fit1)) return (1); /* Pixels at edge's right side */ fit2 = 0.0; n2 = 0; for (i = QXI[q2]+b1; i <= QXI[q2]+b2; i++) { if (PIX(i,j) < fitdata.threshold && !(DQ(i,j) & mask)) { n2++; xval[n2] = (float)(i-QXI[q2]+1); yval[n2] = PIX(i,j); } } if (doFit (xval, yval, n2, eqorder, &fit2)) return (1); if (n1 > 0 && n2 > 0 && fit1 > BADVAL && fit2 > BADVAL) { buffer[(*nb)] = fit2 - fit1; (*nb)++; } } free_vector (xval, 1, b2-b1+1); free_vector (yval, 1, b2-b1+1); return (0); } /* * DOHORI - Do a horizontal boundary. * * Collect quadrant boundary differences for a horizontal boundary. * */ int doHori (float *sci, short *dq, short mask, int b1, int b2, int eqorder, int q1, int q2, float *buffer, int *nb) { int i, j; int n1, n2; float fit1, fit2; float *xval, *yval; float *vector (int nl, int nh); void free_vector (float *v, int nl, int nh); int doFit (float *xval, float *yval, int npts, int order, float *pred); xval = vector (1, b2-b1+1); yval = vector (1, b2-b1+1); *nb = 0; /* Loop from left to right along horizontal boundary */ for (i = QXI[q1]; i <= QXF[q1]; i++) { /* Pixels at edge's bottom side */ fit1 = 0.0; n1 = 0; for (j = QYF[q1]-b1; j >= QYF[q1]-b2; j--) { if (PIX(i,j) < fitdata.threshold && !(DQ(i,j) & mask)) { n1++; xval[n1] = (float)(QYF[q1]-j+1); yval[n1] = PIX(i,j); } } if (doFit (xval, yval, n1, eqorder, &fit1)) return (1); /* Pixels at edge's top side */ fit2 = 0.0; n2 = 0; for (j = QYI[q2]+b1; j <= QYI[q2]+b2; j++) { if (PIX(i,j) < fitdata.threshold && !(DQ(i,j) & mask)) { n2++; xval[n2] = (float)(j-QYI[q2]+1); yval[n2] = PIX(i,j); } } if (doFit (xval, yval, n2, eqorder, &fit2)) return (1); if (n1 > 0 && n2 > 0 && fit1 > BADVAL && fit2 > BADVAL) { buffer[(*nb)] = fit2 - fit1; (*nb)++; } } free_vector (xval, 1, b2-b1+1); free_vector (yval, 1, b2-b1+1); return (0); }