From d54fe7c1f704a63824c5bfa0ece65245572e9b27 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 4 Mar 2015 21:21:30 -0500 Subject: Initial commit --- src/libcf/cf_make_mask.c | 175 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 175 insertions(+) create mode 100644 src/libcf/cf_make_mask.c (limited to 'src/libcf/cf_make_mask.c') diff --git a/src/libcf/cf_make_mask.c b/src/libcf/cf_make_mask.c new file mode 100644 index 0000000..090222c --- /dev/null +++ b/src/libcf/cf_make_mask.c @@ -0,0 +1,175 @@ +/*********************************************************************** + * + * CF_MAKE_MASK + * + * Procedure to create a mask of the bad pixels from the pixel list + * generated by cf_bad_pixels + * + * CALL: cf_make_mask(header, ap, nout, wave_out, bny, bymin, bpmask) + * + * Arguments: fitsfile header name of the input IDF + * int ap Aperture designation + * long nout number of points in wavelength array + * float wave_out wavelength array + * int bny Y dimension of the output image - + * should be same as background + * int bymin Y coordinate of first element in + * output array. Should match background + * float bpmask Output bad pixel mask - dimensions + * are nout * bny + * + * History: 04/29/03 rdr v1.0 Begin work + * 04/30/03 wvd v1.1 Install + * 09/30/03 wvd v1.2 Use cf_read_col to read BPM_CAL + * file. + * 11/11/03 wvd v1.3 Fix bug in normalization of + * pothole map. + * 07/21/04 wvd v1.4 Fix bug in construction of + * pothole mask: change & to && + * 11/25/05 wvd v1.5 Fill gaps in the mask that + * open when the bad-pixel map + * is converted from pixels to + * wavelengths. + * + ***********************************************************************/ + +#include "calfuse.h" + +int cf_make_mask(fitsfile *header, int ap, long nout, float *wave_out, +int bny, int bymin, float **bpmask){ + + char CF_PRGM_ID[] = "cf_make_mask"; + char CF_VER_NUM[] = "1.5"; + + char bpm_file[FLEN_VALUE], *chan; + int *bphist, nhist, status=0; + long bpsum, i, j, k, ndx, nsam, npix, npts; + float *bpmaskt, w0, wpc, pnorm=1; + float *x, *y, *wt, *lam ; + float bpmax ; + fitsfile *bpmfits ; + + cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr); + cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing"); + + cf_verbose(3, "Entering routine cf_make_mask") ; + + /* Generate an array (nout * bny) to contain the map. */ + npix = nout * bny; + bpmaskt = cf_calloc(npix, sizeof(float)) ; + + /* Read wavelength information from input file header. */ + FITS_movabs_hdu(header, 1, NULL, &status); + FITS_read_key(header, TFLOAT, "WPC", &wpc, NULL, &status); + w0 = wave_out[0]; + cf_verbose(3, "Wavelengths: w0=%g, dw=%g", w0, wpc); + + /* Open the bad-pixel file. */ + if (fits_read_key(header, TSTRING, "BPM_CAL", bpm_file, NULL, &status)) { + cf_verbose(0, + "No bad pixel (BPM_CAL) file specified. Assuming no potholes.") ; + for (i=0; i 0 && ndx < npix) bpmaskt[ndx] += wt[i] ; + } + } + } + cf_verbose(3, "Channel %d contains %d bad pixels.", ap, nsam) ; + + /* If any potholes are present, generate the pothole map. */ + if(nsam > 0) { + + /* Determine a normalization for the potholes. This is done by + creating a histogram of the weights between 0 and the maximum + and then selecting the 95% level as the normalization factor. + We use 95%, rather than the maximum, to avoid the possibility + of having a few spurious points define the normalization. */ + + /* Determine the maximum value. */ + bpmax = 0.; + for (i=0; i bpmax) bpmax = bpmaskt[i] ; + + /* Generate the histogram. */ + nhist = (int) (bpmax * 10.) + 1 ; + bphist = (int *) cf_calloc(nhist, sizeof(int)) ; + bpsum = 0; + for (i=0; i < npix; i++) { + ndx = (int) (bpmaskt[i] * 10.) ; + if (ndx > 0 && ndx < nhist) { + bphist[ndx]++; + bpsum++; + } + } + + /* Now select the 95% level. */ + nsam = 0; + i = nhist-1 ; + while (nsam < bpsum/20 && i >= 0) { + nsam += bphist[i] ; + pnorm = bpmax - (nhist-1-i) / 10. ; + i--; + } + + cf_verbose(3, "bpmax = %.1f, normalizing factor = %.2f", bpmax, pnorm) ; + + /* Normalize and invert the mask so that the center of the + pothole is zero and the region outside is 1. */ + + for (i=0; i bpmaskt[ndx+1]) + bpmaskt[ndx] = (bpmaskt[ndx-1] + bpmaskt[ndx+1]) / 2.; + } + } + } + + /* If there are no potholes, generate an array with all ones. */ + else + for (i=0; i