/************************************************************************** * * * Copyright (c) 1996 Michael Richmond and Richard Treffers * * * * This software may be copied and distributed for * * educational, research and not for profit services * * provided that this copyright and statement are * * included in all such copies. * * * **************************************************************************/ /* convolve an image with a gaussian of the given size (plus a little more computation) to create an image which may be better searched for real stars. */ /* * change "apply_gaussian" so that any negative values in the * convolved image are replaced by zero, if the * #define CLIPNEGATIVE is set to 1. * MWR 3/16/1997 */ #undef CLIPNEGATIVE /* if defined, replace negative values in */ /* convolved image with zero. */ #include #include #include #include #include "pcvista.h" #include "fits.h" #include "stars.h" #define CLIPSIG 3.0 /* clip values more than this many stdev */ /* from 0.0 when calculating sky-sig in */ /* sg_skysig routine. */ extern int16 *image[NMAX]; extern int16 *conv_image[NMAX]; double **gauss; double gsum; /* sum of terms in gaussian */ double gsum2; /* sum of squares in terms in gaussian */ double gnum; /* denominator of expression for H */ /* convolve 'image' array with a gaussian which has the given fwhm, and place the convolved version into the 'conv_image' array. return the sigma of the convolved image in 'conv_sig' variable */ static void sg_makegauss( double, int *); static void apply_gaussian (int sr, int sc, int nr, int nc, int ngauss, int signif, int verbose, char *filename); static int16 do_gauss(int, int , int, int); void star_convolve(sr, sc, nr, nc, fwhm, signif, conv_sig, verbose, imagename) int sr, sc, nr, nc, verbose; int signif; double fwhm, *conv_sig; char *imagename; { int ngauss; sg_makegauss(fwhm, &ngauss); apply_gaussian(sr, sc, nr, nc, ngauss, signif, verbose, imagename); *conv_sig = sg_skysig(sr, sc, nr, nc, *conv_sig); if (verbose) { printf("sigma of convolved image is %f\n", *conv_sig); } } /* create a gaussian of the appropriate size for the given FWHM and return in the given parameter its size (N x N). The gaussian is stored as an array of doubles, and is normalized so its sum is zero when applied to a "flat" surface. */ static void sg_makegauss(fwhm, n) double fwhm; int *n; { int i, j; double sigma, mid, ss, x, y, z, num, sum, sum2; sigma = fwhm/2.355; /* allocate space for the array */ if (((i = (int) ((fwhm*2.0) + 0.5)) % 2) == 0) i++; *n = i; if ((gauss = (double **) malloc((*n)*sizeof(double *))) == NULL) { error(1, "sg_makegauss: can't malloc for gauss[] ?! "); } for (i = 0; i < *n; i++) { if ((gauss[i] = (double *) malloc((*n)*sizeof(double))) == NULL) { error(1, "sg_makegauss: can't malloc for gaussian ?!"); } } /* and now fill the array with appropriate numbers */ mid = (*n - 1)/2; ss = 2.0*sigma*sigma; for (i = 0; i < *n; i++) { x = (mid - i); for (j = 0; j < *n; j++) { y = (mid - j); z = (x*x + y*y)/ss; gauss[i][j] = exp(-z); } } /* and now, calculate quantities we'll need later */ num = (*n)*(*n); sum = 0; sum2 = 0; for (i = 0; i < *n; i++) { for (j = 0; j < *n; j++) { sum += gauss[i][j]; sum2 += gauss[i][j]*gauss[i][j]; } } gsum = sum; gsum2 = sum2; gnum = (sum2 - (sum*sum/num)); gsum /= num; } /* apply the gaussian to the given sub-image, which starts at the given (sr, sc) and has a size of nr by nc pixels. Just ignore the outermost pixels, those which are closer than 'ngaussian/2' to an edge. If verbose == 1, then print out messages telling how far along the convolution is. If filename != "", then write the result of the convolution into a new FITS file on disk with the given name. */ static void apply_gaussian(sr, sc, nr, nc, ngauss, signif, verbose, filename) int sr, sc, nr, nc, ngauss, verbose; char *filename; int signif; { int i, j, k, offset, flag; int16 data[NMAX]; FITS_HANDLE fout; /* placate compiler */ fout = (FITS_HANDLE) 0; if (strcmp(filename, "") != 0) { flag = 1; } else { flag = 0; } /* first, open the new FITS file for writing */ if (flag == 1) { fout = fits_open(filename, "w", &nr, &nc); } offset = ngauss/2; /* first, write zeros in all rows too close to the near edge */ if (flag == 1) { for (j = 0; j < nc; j++) data[j] = 0; for (i = 0; i < offset; i++) fits_put_data(fout, i, 0, data, nc); } /* now, write real numbers to all of the image that we can */ for (i = sr + offset; i < (sr + nr) - offset; i++) { if ((verbose == 1) && (i % 20 == 0)) { printf("\r convolution up to row %4d of %4d ", i, (sr + nr) - offset); fflush(stdout); } for (j = sc + offset; j < (sc + nc) - offset; j++) { k = j - sc; data[k] = do_gauss(i, j, ngauss, signif); /* * tests on TASS images show that the negative * "moats" around star centers in the convolved * image can cause the centroiding routine to go * haywire. Let's replace negative values with zero * in the convolved image, to avoid this. * MWR 3/16/1997 */ #ifdef CLIPNEGATIVE if (data[k] < 0) { data[k] = 0; } #endif conv_image[i][j] = data[k]; } if (flag == 1) fits_put_data(fout, i - sr, 0, data, nc - offset); } if (verbose) printf("\n"); /* finally, write zeros in all rows too close to the far edge */ if (flag == 1) { for (j = 0; j < nc; j++) data[j] = 0; for (i = nr - offset; i < nr; i++) fits_put_data(fout, i, 0, data, nc); fits_close(fout); } } /* return the value of the convolution of the image and the gaussian at the given row and col. actually, do a little more than just the convolution; in effect, fit the gaussian to the image as described in Stetson's DAOPHOT paper. return an int16, rather than a double. */ static int16 do_gauss(row, col, ngauss, signif) int row, col, ngauss; int16 signif; { int i, j, sr, sc; int16 *q; double sum, dsum, *dq; sr = row - ngauss/2; sc = col - ngauss/2; sum = 0; dsum = 0; for (i = 0; i < ngauss; i++) { q = image[sr + i] + sc; dq = gauss[i]; for (j = 0; j < ngauss; j++) { sum += ((double) (*q))*(*dq++); dsum += *q++; } } sum = (sum - dsum*gsum)/gnum; /* this is the "extra" bit */ if (sum > MAX_DATA_VAL) { return((int16) MAX_DATA_VAL); } else if (sum < MIN_DATA_VAL) { return((int16) MIN_DATA_VAL); } else { return((int16) sum); } } /* calculate the standard deviation of pixels from the mean in the convolved image. We are a little tricky here: - we know the mean of the convolved image will be very close to 0.0 - we can use the "rough_sig" argument (which is based on the standard deviation of sky in the UN-convolved image) as a guess at the width of the sky in the convolved image So, what we do is calculate a CLIPPED mean and standard deviation, using only pixels which have values between -CLIPSIG*rough_sig < value < CLIPSIG*rough_sig and we ignore all others. Return a value of 1.0 if we can't calculate the standard deviation formally. */ double sg_skysig(sr, sc, nr, nc, rough_sig) int sr, sc, nr, nc; double rough_sig; { int row, col; int16 number; double mean, diff, stdev, total, pixnum, sumsq; double minval, maxval; total = 0.0; pixnum = 0.0; sumsq = 0.0; minval = -CLIPSIG*rough_sig; maxval = CLIPSIG*rough_sig; for (row = sr; row < (sr + nr); row++) { for (col = sc; col < (sc + nc); col++) { number = conv_image[row][col]; if ((number >= minval) && (number <= maxval)) { total += (double)number; sumsq += (double)number*(double)number; pixnum++; } } } if (pixnum < 2) { return(1.0); } mean = total/pixnum; diff = (sumsq - pixnum*mean*mean); if (diff > 0.0) { stdev = sqrt(diff/(pixnum - 1.0)); } else { stdev = 1.0; } return(stdev); }