/************************************************************************** * * * 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. * * * **************************************************************************/ /* given a list of star centroids for a CCD frame, do primitive aperture photometry for each star. PHOT image [INFILE= ] [OUTFILE= ] [BIN= ] [MCONST= ] [ADU=] [READNOISE=] [SKYINNER= ] [SKYOUTER=] [APER= ] ... [VERBOSE] [MEAN] [MODE] [SKY= ] [SKYSIG= ] [SATURATE= ] [FLUXZERO= ] [MASKFILE= ] uses a LOCAL value for the sky background for each star. 9/17/91 - fixed bug in 'get_sky()' routine by adding check to see that bin number is reasonable. MWR 4/9/92 - corrected calculation of estimate uncertainty. MWR 4/16/92 - made saturated stars' uncertanties print out as negative numbers. MWR 8/20/92 - added external sky value -rrt 10/29/92- changed from 'long' to 'double' some variables in the get_sky() routine. overflow was causing errors. also, added command-line option for SATURATE. also, changed 'get_sky()' to switch from using mode to using mean (even if not given MEAN flag) whenever the sky histogram has peak in top or bottom bin. MWR 11/19/92- changed the way in which pixels which fall PARTLY inside the perfectly circular aperture are handled. old version took whole value of pixel whose center was a distance APER (or less) away from center of pixel nearest to centroid. new version takes a fraction of each pixel which falls only partly within the aperture, now centered on the actual stellar centroid, rather than the nearest pixel. the fraction taken is very close the fraction of the area which falls within the circular aperture. MWR 12/11/92- fixed a place near line 411 where it was possible to take sqrt of negative number. Not any more! MWR 12//17/92 - fixed another place near the above where log is taken rrt 8/15/94 - add exit(0); 7/18/95 - increase APER_SIZE 30--> 100 I don't see what it does -rrt [sets limit on size of user-supplied apertures. MWR] 1/7/96 - added the "help" option, which now prints the default values. MWR 2/23/96 - changed the way we calculate sky values default is now to use (binned) median value "mode" keyword forces use of the mode (old default) "mean" keyword forces use of interquartile mean use standard deviation from interquartile mean as the reported "stdev" for all sky values MWR 3/16/97 - added "fluxzero=" option, where value is the number of ADU corresponding to no flux. For TASS data, reduced by me, this is -32768; for most people, one would expect 0. Make default=0. MWR 3/8/2001 - fixed calculations of noise due to sky and readnoise in 'do_phot()'. We had been using 'num_star_pix[]' as the number of pixels within an aperture, but that counts each fractional pixel as a whole one, so it was an overestimate. We now calculate an accurate number of pixels within the aperture as local array 'true_pix_in_aper[]'. This should decrease the improperly large noise values, and yield improved (and smaller) uncertainties for each magnitude measurement. MWR 5/20/2001 - No use a more accurate calculation of the fraction of each pixel within the synthetic aperture. The new 'frac_pixel_inside_aperture' routine uses a brute-force approach, dividing each pixel on the boundary into sub-pixels. Is slower than old method, but more accurate. MWR 5/27/2001 - Added new 'maskfile=' option, to specify a file containing list of bad regions. The 'phot' program will now print output with a NEW COLUMN: this column is a flag which describes each stellar measurement; it can notify that star is in a bad region, or that star is saturated, or other information. */ #include #include #include #include #include #include "pcvista.h" #include "fits.h" #undef DEBUG /* give lots of information about program's progress */ #ifndef M_PI #define M_PI 3.14159265359 #endif #define YES 1 #define NO 0 #define NUM_APER 10 /* max number of apertures which can be used at once */ #define APER_SIZE 100 /* maximum radius of any star or sky aperture */ #define SKY_INNER 10 /* default val for inner radius of sky annulus */ #define SKY_OUTER 20 /* default val for outer radius of sky annulus */ #define STAR_RAD 3 /* default val for outer radius of star aperture */ #define ADU 1 /* default # of electrons per adu */ #define READNOISE 1 /* default val for readout noise of chip (electrons) */ #define SATURATE 32767 /* if max value >= this, star is saturated */ #define MAG_CONST 25 /* arbitrary constant added to reported star mags */ /* (since are all negative, mags above sky) */ #define FLUXZERO 0 /* this many ADU in a pixel means "no flux" */ #define BAD_STAR_MAG 99 /* magnitude reported for a 'negative' star, or */ /* one with any sort of problem */ #define BIN_SIZE 1 /* default binsize used to find mode of the sky */ #define NUM_BIN 65535 /* number of bins used to find mode of sky */ #define NEG_VALUE 32768 /* add this to all pixel values when getting */ /* the mode, so that negative values don't all */ /* add up in the 'value=0' bin. in other words, */ /* any neg pixel values should be > NEG_VALUE */ #define SKY_MEDIAN 0 /* use median in annlus for sky */ #define SKY_MEAN 1 /* use interquartile mean in annulus for sky */ #define SKY_MODE 2 /* use mode in annlus for sky */ #define LINELEN 100 /* max length of lines in mask file */ #define FLAG_OK 0 /* no problems with the stellar measurement */ #define FLAG_SATURATED 1 /* star may be saturated */ #define FLAG_MASK 2 /* star falls close to masked region */ #define FLAG_EDGE 4 /* star falls close to edge of image */ /* this structure holds information about regions of "bad" pixels */ typedef struct region_struct { int region_id; /* ID number for internal use */ int numpix; /* number of connected bad pixels in region */ int min_row; /* minimum row of bounding box */ int max_row; /* maximum row of bounding box */ int min_col; /* minimum col of bounding box */ int max_col; /* maximum col of bounding box */ struct region_struct *next; } REGION_STRUCT; int main(int, char **); static void print_usage(FILE *fp); void make_masks(void); int do_phot(void); int get_sky(int, int, int, double *); void get_star(double, double, long); static double frac_pixel_inside_aperture(double aper_row_center, double aper_col_center, double aper_radius, int pixel_row, int pixel_col); static REGION_STRUCT * read_regions(FILE *maskfp); static void free_regions(REGION_STRUCT *head); static int is_close_to_bad_region(double row, double col, int largest_aper_radius); static int is_close_to_edge(double row, double col, int largest_aper_radius); /* the various arrays we have to deal with */ static int *sky_r, *sky_c; /* sky aperture mask arrays */ static int *star_r[NUM_APER], *star_c[NUM_APER]; /* star aperture mask arrays */ static int *bin; /* bins for calculating sky mode value */ static int16 data[NMAX]; /* array used to read in data from file */ static long star_counts[NUM_APER]; /* counts w/in each aperture */ /* some other global variables */ static FITS_HANDLE fh; /* image data file */ static int nrow, ncol; /* number of rows and cols in the image */ static double sky_inner; /* sky annulus inner radius */ static double sky_outer; /* sky annulus outer radius */ static double star_rad[NUM_APER]; /* star aperture radii */ static int num_aper; /* number of specified star apertures to use */ static int max_aper; /* largest radius aperture specified */ static int num_sky_pix; /* number of pixels in sky mask */ static int num_star_pix[NUM_APER]; /* number of pixels in each star mask */ static int num_bin; /* no. of bins in sky mode-finding */ static int bin_size; /* size of each bin */ static double mag_const; /* arbitrary const added to mags */ static double adu; /* # of electrons/adu for the chip */ static double readnoise; /* readout noise of chip */ static double saturate; /* saturation point of chip, in ADU */ static double fluxzero; /* this many ADU in a pixel means "no flux" */ static int verbose; /* if verbose=1, print out info on star candidates */ static int skyflag; /* indicates algorithm to use in sky calc */ FILE *infp; /* file from which input is read (or stdin, if none) */ FILE *outfp; /* file into which output goes (or stdout, if none) */ FILE *maskfp; /* file with list of bad regions in image */ static REGION_STRUCT *region_list = NULL; static int external_sky_flag; /* whether to use user supplied sky */ static double external_sky_value; static double external_skysigma=1.0; int main(argc,argv) int argc; char *argv[]; { static int i; static int nr, nc, boxno; char *gotit; static char fname[NBUF], buf[NBUF]; boxno = 0; outfp = stdout; infp = stdin; maskfp = NULL; sky_inner = (double) SKY_INNER; sky_outer = (double) SKY_OUTER; mag_const = (double) MAG_CONST; adu = (double) ADU; readnoise = (double) READNOISE; num_aper = 0; max_aper = (int) (STAR_RAD + 0.5); bin_size = BIN_SIZE; num_bin = NUM_BIN; skyflag = SKY_MEDIAN; verbose = 0; saturate = SATURATE; fluxzero = FLUXZERO; if (argc < 2) { print_usage(stderr); exit(-1); } if (strcmp(argv[1], "help") == 0) { print_usage(stdout); exit(0); } strcpy(fname, argv[1]); fh = fits_open(fname, "r", &nrow, &ncol); nr = nrow; nc = ncol; /* take care of the command-line arguments */ for (i = 2; i < argc; i++) { if (keyword("help", argv[i])) { print_usage(stdout); exit(0); } if ((gotit = find("infile", argv[i])) != NULL) { if ((infp = fopen(gotit, "r")) == NULL) error(-1, "phot: can't open input file"); } if ((gotit = find("outfile", argv[i])) != NULL) { if ((outfp = fopen(gotit, "w")) == NULL) error(-1, "phot: can't open output file"); } if ((gotit = find("maskfile", argv[i])) != NULL) { if ((maskfp = fopen(gotit, "r")) == NULL) error(-1, "phot: can't open mask file"); } if ((gotit = find("adu", argv[i])) != NULL) { adu = evaluate(gotit); } if ((gotit = find("readnoise", argv[i])) != NULL) { readnoise = evaluate(gotit); } if ((gotit = find("aper", argv[i])) != NULL) { if (++num_aper > NUM_APER) { sprintf(fname, "phot: too many apertures - can only handle %d\n", NUM_APER); error(-1, fname); } star_rad[num_aper - 1] = evaluate(gotit); if ((int) (star_rad[num_aper - 1] + 0.5) > max_aper) max_aper = (int) (star_rad[num_aper - 1] + 0.5); } if ((gotit = find("skyinner", argv[i])) != NULL) { sky_inner = evaluate(gotit); continue; } if ((gotit = find("skyouter", argv[i])) != NULL) { sky_outer = evaluate(gotit); continue; } if ((gotit = find("sky", argv[i])) != NULL) { external_sky_value = evaluate(gotit); external_sky_flag=1; continue; } if ((gotit = find("skysig", argv[i])) != NULL) { external_skysigma = evaluate(gotit); continue; } if ((gotit = find("saturate", argv[i])) != NULL) { saturate = evaluate(gotit); continue; } if ((gotit = find("fluxzero", argv[i])) != NULL) { fluxzero = evaluate(gotit); continue; } if ((gotit = find("bin", argv[i])) != NULL) { bin_size = (int) (evaluate(gotit) + 0.5); } if ((gotit = find("mconst", argv[i])) != NULL) { mag_const = evaluate(gotit); } if (keyword("mean", argv[i])) { skyflag = SKY_MEAN; } if (keyword("mode", argv[i])) { skyflag = SKY_MODE; } if (keyword("verbose", argv[i])) { verbose = 1; } } if (num_aper == 0) { num_aper = 1; star_rad[0] = (double) STAR_RAD; } if ((sky_inner < 0.0) || (sky_inner >= (double) APER_SIZE)) error(-1, "phot: invalid value for SKYINNER parameter"); if ((sky_outer < 0.0) || (sky_outer >= (double) APER_SIZE)) error(-1, "phot: invalid value for SKYOUTER parameter"); for (i = 0; i < num_aper; i++) if ((star_rad[i] < 0.0) || (star_rad[i] >= (double) APER_SIZE)) error(-1, "phot: invalid value for APER parameter"); if (adu <= 0.0) error(-1, "phot: invalid adu input - must be greater than zero"); if (readnoise < 0.0) error(-1, "phot: invalid readnoise input - must be at least zero"); /* try to allocate enough memory for the mode-finding bins */ if ((bin = (int *) malloc(num_bin*sizeof(int))) == NULL) error(-1, "phot: can't allocate space for enough bins - try larger binsize"); /* give a warning message if there may not be enough bins to find the sky value correctly */ if (((long) bin_size)*num_bin < ((long) MAX_DATA_VAL) + NEG_VALUE) { if (verbose) { sprintf(buf, "warning: bin size may be too small if sky values near %ld", ((long) bin_size)*num_bin); error(0, buf); } } /* create the sky annulus mask and the star aperture masks */ make_masks(); /* * if the user specified a file with list of bad regions, * read that file and create a linked list of regions */ if (maskfp != NULL) { region_list = read_regions(maskfp); } /* and now go to the photometry */ do_phot(); if (infp != stdin) fclose(infp); if (outfp != stdout) fclose(outfp); if (maskfp != NULL) { free_regions(region_list); fclose(maskfp); } exit(0); } /* * print out the list of arguments and their default values */ static void print_usage(fp) FILE *fp; { int i; fprintf(fp, "usage: phot image \n"); fprintf(fp, " [infile=] [outfile=] [maskfile=] \n"); fprintf(fp, " [bin=%d]\n", bin_size); fprintf(fp, " [mconst=%.2f]\n", mag_const); fprintf(fp, " [adu=%.2f]\n", adu); fprintf(fp, " [readnoise=%.2f]\n", adu); fprintf(fp, " [skyinner=%.2f] [skyouter=%.2f]\n", sky_inner, sky_outer); if (num_aper == 0) { fprintf(fp, " [aper=] ... \n"); } else { fprintf(fp, " "); for (i = 0; i < num_aper; i++) { fprintf(fp, "aper=%.2f ", star_rad[i]); } fprintf(fp, "\n"); } fprintf(fp, " [sky=%.2f] [skysig=%.2f]\n", external_sky_value, external_skysigma); fprintf(fp, " [saturate=%.2f]\n", saturate); fprintf(fp, " [fluxzero=%.2f]\n", fluxzero); fprintf(fp, " [mean] [mode] [verbose]\n"); } /* create lists of row- and col- offsets, in pixels, from the center of a star, that we should use to calculate 1. the sky value (sky mask) 2. the star brightness (star mask(s)) the n'th element of each array specfies a pixel offset (in row and col) that we use to add to the sky or star sum. */ void make_masks() { static int i, j, k, n, i2; static double d, inner2, outer2; /* first, allocate space for the sky arrays */ n = 4 * (int) (sky_outer*sky_outer); /* just to be safe */ if ((sky_r = (int *) malloc(n*sizeof(int))) == NULL) error(-1, "can't malloc for sky r array"); if ((sky_c = (int *) malloc(n*sizeof(int))) == NULL) error(-1, "can't malloc for sky c array"); /* next, allocate space for all of the star arrays */ for (i = 0; i < num_aper; i++) { if (star_rad[i] < 2.0) n = 16; else n = 4 * (int) ((star_rad[i] + 1)*(star_rad[i] + 1)); if ((star_r[i] = (int *)malloc(n*sizeof(int))) == NULL) error(-1,"can't malloc for star mask array"); if ((star_c[i] = (int *)malloc(n*sizeof(int))) == NULL) error(-1,"can't malloc for star mask array"); } /* now, fill in the sky mask arrays */ outer2 = sky_outer * sky_outer; inner2 = sky_inner * sky_inner; num_sky_pix = 0; for (i = (int)(-(sky_outer + 1)); i <= (int)(sky_outer + 1); i++) { i2 = i*i; for (j = (int)(-(sky_outer + 1)); j <= (int)(sky_outer + 1); j++) { d = (double) (i2 + j*j); if ((d <= outer2) && (d >= inner2)) { #ifdef DEBUG fprintf(stderr, "sky mask: add offset (%5d, %5d) \n", i, j); #endif sky_r[num_sky_pix] = i; sky_c[num_sky_pix] = j; num_sky_pix++; } } } /* and now fill in the star mask arrays */ for (k = 0; k < num_aper; k++) { outer2 = (star_rad[k] + 1)*(star_rad[k] + 1); num_star_pix[k] = 0; for (i = (int)(-(star_rad[k] + 1)); i <= (int)(star_rad[k] + 1); i++) { i2 = i*i; for (j = (int)(-(star_rad[k] + 1)); j <= (int)(star_rad[k] + 1); j++) { if ((double) (i2 + j*j) <= outer2) { #ifdef DEBUG fprintf(stderr, "star mask %d: add offset (%5d, %5d) \n", k, i, j); #endif star_r[k][num_star_pix[k]] = i; star_c[k][num_star_pix[k]] = j; num_star_pix[k]++; } } } } } /* do aperture photometry around each star centroid given in the input (assumed to be the output of FIND). returns 0 when it reaches end-of-file, or 1 if an invalid line of input is read */ int do_phot() { static int i, j, row, col, star_num; static int quality_flag; static int largest_aper_radius; static long sky, peakval; static char buf[NBUF]; static double star_mag[NUM_APER], star_err[NUM_APER], temp; static double row_cen, col_cen, sky_sig; static double true_pix_in_aper[NUM_APER]; /* * find the largest radius of all the apertures * (round to next-highest integer if necessary). We use this * to determine if the star is "close" to the edge of the image. */ largest_aper_radius = star_rad[0]; for (i = 1; i < num_aper; i++) { if (star_rad[i] > largest_aper_radius) { largest_aper_radius = star_rad[i]; } } /* * Pre-calculate the true number of pixels within each aperture, * including only the fraction of each pixel which falls within * the ideal circular aperture. We'll use this when calculating * the contribution to noise from sky and readout noise. * * This isn't the same as 'num_star_pix' -- because 'num_star_pix' * is a count of the number of pixels which fall within the aperture, * even by only a tiny fraction of a pixel. Thus, we should * find * true_pix_in_aper[i] < num_star_pix[i] * if we do the job right. * * -- MWR 3/8/2001 */ for (i = 0; i < num_aper; i++) { true_pix_in_aper[i] = M_PI*star_rad[i]*star_rad[i]; if (true_pix_in_aper[i] >= num_star_pix[i]) { error(1, "do_phot: true_pix_in_aper > num_star_pix?"); } } for (i = 1; ; i++) { /* first, get the next input line */ if (fgets(buf, NBUF, infp) == NULL) return(0); if (sscanf(buf, "%d %lf %lf %ld", &star_num, &row_cen, &col_cen, &peakval) != 4) { sprintf(buf, "invalid input in line %d\n", i); error(0, buf); return(1); } if (verbose) fprintf(stderr, "starting on star %5d, row=%6.3f col=%6.3f\n", star_num, row_cen, col_cen); /* start this star with no problems */ quality_flag = FLAG_OK; row = (int) (row_cen + 0.5); col = (int) (col_cen + 0.5); /* get the sky value for this star */ if(external_sky_flag){ sky=external_sky_value; sky_sig=external_skysigma; }else{ sky = get_sky(row, col, skyflag, &sky_sig); } peakval += sky; if (verbose) fprintf(stderr, " sky for this star is %5ld, sky_sig = %f\n", sky, sky_sig); /* now, for each aperture, add up the signal inside the aperture, subtracting the given sky value as we go. store the results inside the 'star_counts[]' array */ get_star(row_cen, col_cen, sky); /* for each aperture, convert the result to magnitudes */ for (j = 0; j < num_aper; j++) { if (verbose) { fprintf(stderr, " (star-sky) counts for aperture %d = %8ld\n", j, star_counts[j]); } if (star_counts[j] <= 0) { star_mag[j] = (double) BAD_STAR_MAG; star_err[j] = (double) BAD_STAR_MAG; } else { star_mag[j] = mag_const - 2.5*log10((double) star_counts[j]); /* new version does noise calculations in terms of electrons, not ADU */ /* * note that the equation below ASSUMES that the "sky" value * has a zero point at the level of no flux. This is appropriate * in the usual case, when one has subtracted bias and dark * current; but the user may have chosen to keep the full * 16-bit dynamic range in XVista by setting the zero-flux * level to -32768. The default is "fluxzero=0", but * the user may have changed it on the command line. * In any case, subtract "fluxzero" from sky value when * calculating the magnitude uncertainty. * MWR 3/16/97 */ temp = ((double) true_pix_in_aper[j]*(sky - fluxzero)*adu) + ((double) true_pix_in_aper[j]*readnoise*readnoise) + ((double) star_counts[j]*adu); if (temp > 0.){ temp = sqrt(temp); star_err[j] = 2.5*log10(1.0 + (temp/(star_counts[j]*adu))); if (verbose) { fprintf(stderr, " noise for aperture %d = %8.1f \n", j, temp); } }else { star_err[j] = (double)(BAD_STAR_MAG); if (verbose) { fprintf(stderr, " noise for aperture %d = %8.1f \n", j, (double) BAD_STAR_MAG); } } if (peakval >= saturate) { quality_flag |= FLAG_SATURATED; } } if (verbose) { fprintf(stderr, " mag for aperture %d = %6.2f +/- %6.3f\n", j, star_mag[j], star_err[j]); } } /* check to see if this star lies close to a bad region */ if (is_close_to_bad_region(row_cen, col_cen, largest_aper_radius) == 1) { quality_flag |= FLAG_MASK; } /* check to see if this star lies close to an edge of the image */ if (is_close_to_edge(row_cen, col_cen, largest_aper_radius) == 1) { quality_flag |= FLAG_EDGE; } /* now spit out the results */ sprintf(buf, "%5d %7.2f %7.2f %6ld %6.2f ", star_num, row_cen, col_cen, sky, sky_sig); for (j = 0; j < num_aper; j++) sprintf(buf, "%s %6.3f %6.3f ", buf, star_mag[j], star_err[j]); sprintf(buf, "%s %3d ", buf, quality_flag); fprintf(outfp, "%s\n", buf); } } /* return the sky value inside the mask around the given star - there are several choices for algorithm, depending on the value of "skyflag": skyflag = SKY_MEDIAN use binned median (this is default) skyflag = SKY_MEAN use interquartile mean skyflag = SKY_MODE use the mode Always place the standard deviation of interquartile mean into the 'std_dev' parameter. */ int get_sky(row, col, skyflag, std_dev) int row, col, skyflag; double *std_dev; { static int i, r, c, max, mode, nb; static int sr, er, sc, ec, nc; static int min_bin, max_bin; static int first_flag = 0; static int q25_bin, q50_bin, q75_bin, qmedian; static double x; static double sum, sum2, num; static double mean, stdev; static double q25, q50, q75, qmean, qstdev; num = 0.0; sum = 0.0; sum2 = 0.0; if (first_flag == 0) { for (i = 0; i < num_bin; i++) bin[i] = 0; } else { for (i = min_bin; i <= max_bin; i++) { bin[i] = 0; } } min_bin = NUM_BIN; max_bin = 0; /* the following assumes that sky_r[0] has the maximum size of the sky annulus */ if ((sr = row + sky_r[0]) < 0) sr = 0; if ((er = row - sky_r[0]) >= nrow) er = nrow; if ((sc = col + sky_r[0]) < 0) sc = 0; if ((ec = col - sky_r[0]) >= ncol) ec = ncol - 1; nc = (ec - sc) + 1; i = 0; for (r = sr; r < er; r++) { if (fits_get_data(fh, r, sc, data, nc) == FITS_FAIL) error(-1, "error reading data from file"); while ((i < num_sky_pix) && (sky_r[i] < r - row)) i++; while ((i < num_sky_pix) && (sky_r[i] == r - row)) { if (((c = col + sky_c[i] - sc) >= 0) && (c < nc)) { #ifdef DEBUG2 printf(" r %6d c %6d %6d\n", r, c, data[c]); #endif x = data[c]; sum += x; sum2 += x*x; num++; if ((nb = (data[c] + (int) NEG_VALUE)/bin_size) < 0) nb = 0; else if (nb >= NUM_BIN) nb = NUM_BIN - 1; if (nb < min_bin) { min_bin = nb; } else if (nb > max_bin) { max_bin = nb; } bin[nb]++; } i++; } } #ifdef DEBUG printf(" min_bin %6d max_bin %6d\n", min_bin, max_bin); #endif /* * here we calculate the mean and standard deviation, using all * bins. someone might want to know these numbers sometime, * but at the present, we never return them. * * MWR 2/23/97 */ if (num > (long) 0) mean = ((double)sum)/num; else mean = 0.0; if (num > (long) 1) stdev = sqrt(((double)sum2 - num*mean*mean)/((double)(num - 1))); else stdev = 0.0; /* * Now calculate the mode */ max = 0; mode = 0; for (i = min_bin; i <= max_bin; i++) { #ifdef DEBUG2 if (bin[i] > 0) { printf(" %6d %8d\n", i, bin[i]); } #endif if (bin[i] > max) { max = bin[i]; mode = i; #ifdef DEBUG printf(" new mode is %6d has %5d pix\n", mode, max); #endif } } /* if the mode is the 0'th bin, or the top bin, it's likely that some kind of bad data has gotten into the frame. In such a case, use the mean as the sky value */ if ((mode == 0) || (mode == (num_bin - 1))) { mode = (int) mean; } else { mode = (mode*bin_size) + (bin_size/2) - NEG_VALUE; } /* * find the 25'th and 75'th quartiles in the sky histogram, * and calc the mean and stdev of value between these quartiles only. */ q25 = 0.25*num; q50 = 0.50*num; q75 = 0.75*num; q25_bin = -1; q50_bin = -1; q75_bin = -1; qmean = -1; qstdev = -1; sum = 0; for (i = min_bin; i <= max_bin; i++) { sum += bin[i]; if ((q25_bin == -1) && (sum >= q25)) { q25_bin = i; } if ((q50_bin == -1) && (sum >= q50)) { q50_bin = i; } if ((q75_bin == -1) && (sum >= q75)) { q75_bin = i; break; } } #ifdef DEBUG printf(" q25 %6.0f at bin %5d q50 %6.0f at bin %5d q75 %6.0f at bin %5d\n", q25, q25_bin, q50, q50_bin, q75, q75_bin); #endif if (q75 > q25) { num = 0; sum = 0; sum2 = 0; for (i = q25_bin; i <= q75_bin; i++) { x = bin_size*i; sum += x*(bin[i]); sum2 += x*x*(bin[i]); num += bin[i]; } if (num > (long) 0) { qmean = ((double)sum)/num; } else { qmean = 0.0; } if (num > (long) 1) { qstdev = sqrt((sum2 - num*qmean*qmean)/((double)(num - 1))); } else { qstdev = 0.0; } qmean -= NEG_VALUE; qmedian = q50_bin - NEG_VALUE; } #ifdef DEBUG printf(" mean %6.0f qmn %6.0f qstd %5.0f qmedian %6d mode %6d ", mean, qmean, qstdev, q50_bin - NEG_VALUE, mode); #endif /* * now, decide which value to return, and set the "std_dev" output * argument to the "qstdev", if it has been calculated, or * to the overall stdev otherwise. */ if (qstdev != -1.0) { *std_dev = qstdev; } else { *std_dev = stdev; } if (skyflag == SKY_MEDIAN) { return(qmedian); } else if (skyflag == SKY_MEAN) { return(qmean); } else if (skyflag == SKY_MODE) { return(mode); } else { error(1, "phot: get_sky has invalid value of skyflag?"); } /* we never get here, but this placates compiler */ return(qmedian); } /* put the sum of all pixels within the each star aperture around the given row and col into the 'star_counts[]' array. note how much we have to go through because we can't fit the whole image into memory at once! */ /* modified to use a fractional-pixel weighting scheme where necessary, based on the true row and col centroid. The fractional-pixel weighting uses a linear approximation to the true fractional area included as a function of radius, but it's pretty good, especially for large apertures (radius >= 3 pixels). MWR 11/18/92. */ /* also modified so that the SKY is subtracted from each fractional pixel, as well. MWR 11/18/92 */ #define ROOT2OVER2 (sqrt(2)/2.0) void get_star(row_cen, col_cen, sky) double row_cen, col_cen; long sky; { static int i, r, c, sr, er, sc, ec, nc, aper, counted[NUM_APER]; static double min_dist_sq[NUM_APER], max_dist_sq[NUM_APER]; int aper_row_index, aper_col_index; int addnum; double dr, dc, dist_sq; double fraction; /* * the exact center of the aperture for this star lies somewhere * within the pixel with these integer coords */ aper_row_index = floor(row_cen); aper_col_index = floor(col_cen); #ifdef DEBUG fprintf(stderr, "into get_star for row=%f (%d), col=%f (%d)\n", row_cen, aper_row_index, col_cen, aper_col_index); #endif /* max_aper has the largest aperture size */ if ((sr = aper_row_index - (max_aper + 1)) < 0) sr = 0; if ((er = aper_row_index + (max_aper + 1)) >= nrow) er = (nrow - 1); if ((sc = aper_col_index - (max_aper + 1)) < 0) sc = 0; if ((ec = aper_col_index + (max_aper + 1)) >= ncol) ec = (ncol - 1); nc = (ec - sc) + 1; for (i = 0; i < NUM_APER; i++) { counted[i] = 0; star_counts[i] = (long) 0; } /* * calculate the critical distances we need to use when deciding * whether to include a pixel in the aperture; and, if we include * it, how much to include. */ for (aper = 0; aper < num_aper; aper++) { min_dist_sq[aper] = (star_rad[aper] - ROOT2OVER2)* (star_rad[aper] - ROOT2OVER2); max_dist_sq[aper] = (star_rad[aper] + ROOT2OVER2)* (star_rad[aper] + ROOT2OVER2); } for (r = sr; r <= er; r++) { dr = ((r + 0.5) - row_cen)*((r + 0.5) - row_cen); if (fits_get_data(fh, r, sc, data, nc) == FITS_FAIL) error(-1, "error reading data from file"); for (aper = 0; aper < num_aper; aper++) { if (r - aper_row_index < star_r[aper][counted[aper]]) continue; while ((counted[aper] < num_star_pix[aper]) && (star_r[aper][counted[aper]] < r - aper_row_index)) counted[aper]++; while ((counted[aper] < num_star_pix[aper]) && (star_r[aper][counted[aper]] == r - aper_row_index)) { if (((c = aper_col_index + star_c[aper][counted[aper]] - sc) >= 0) && (c < nc)) { dc = ((c + sc + 0.5) - col_cen)*((c + sc + 0.5) - col_cen); dist_sq = (dr + dc); if (dist_sq < min_dist_sq[aper]) { /* the entire pixel falls within the aperture */ addnum = data[c] - sky; star_counts[aper] += (long) addnum; #ifdef DEBUG fprintf(stderr, " aper %d: try pix %5d: offset (%5d, %5d)\n", aper, counted[aper], star_r[aper][counted[aper]], star_c[aper][counted[aper]]); fprintf(stderr, " dist_sq = %f about to add whole value %d\n", dist_sq, addnum); fprintf(stderr, " aper %d: added %5d in, now sum is %8ld\n", aper, addnum, star_counts[aper]); #endif } else if (dist_sq < max_dist_sq[aper]) { /* pixel falls partially within aperture */ fraction = frac_pixel_inside_aperture(row_cen, col_cen, star_rad[aper], r, c + sc); addnum = (int) ((data[c] - sky)*fraction); star_counts[aper] += (long) addnum; #ifdef DEBUG fprintf(stderr, " aper %d: try pix %5d: offset (%5d, %5d)\n", aper, counted[aper], star_r[aper][counted[aper]], star_c[aper][counted[aper]]); fprintf(stderr, " dist_sq = %f about to add %4.2f*%d = %d\n", dist_sq, fraction, data[c], addnum); fprintf(stderr, " aper %d: added %5d in, now sum is %8ld\n", aper, addnum, star_counts[aper]); #endif } else { /* pixel falls entirely outside aperture */ #ifdef DEBUG fprintf(stderr, " aper %d: reject pix %5d: offset (%5d, %5d)\n", aper, counted[aper], star_r[aper][counted[aper]], star_c[aper][counted[aper]]); #endif } } counted[aper]++; } } } } /*************************************************************************** * PROCEDURE: frac_pixel_inside_aperture * * DESCRIPTION: calculate the fraction of the given pixel which falls * within the photometric aperture. * * We use a brute force method: break the pixel into * NUM_SUB_PIXEL sub-pixels, and check to see if each one * falls within the aperture. * * RETURNS: * fraction of pixel within aperture (between 0 and 1) */ #define NUM_SUB_PIXEL 100 static double frac_pixel_inside_aperture ( double aper_row_center, /* I: central row of aperture */ double aper_col_center, /* I: central col of aperture */ double aper_radius, /* I: radius of aperture (pixels) */ int pixel_row, /* I: row of this pixel */ int pixel_col /* I: col of this pixel */ ) { int row_index, col_index; double aper_radius_sq; double sub_pix_size, sub_pix_area; double sub_row_center, sub_col_center; double drow_sq, dcol_sq, dist_sq; double sum; aper_radius_sq = aper_radius*aper_radius; sub_pix_size = 1.0/NUM_SUB_PIXEL; sub_pix_area = sub_pix_size*sub_pix_size; sum = 0.0; for (row_index = 0; row_index < NUM_SUB_PIXEL; row_index++) { sub_row_center = pixel_row + (row_index + 0.5)*sub_pix_size; drow_sq = (sub_row_center - aper_row_center)* (sub_row_center - aper_row_center); if (drow_sq > aper_radius_sq) { continue; } for (col_index = 0; col_index < NUM_SUB_PIXEL; col_index++) { sub_col_center = pixel_col + (col_index + 0.5)*sub_pix_size; dcol_sq = (sub_col_center - aper_col_center)* (sub_col_center - aper_col_center); dist_sq = drow_sq + dcol_sq; if (dist_sq < aper_radius_sq) { sum += sub_pix_area; } } } return(sum); } /************************************************************************** * PROCEDURE: read_regions * * DESCRIPTION: given an open FILE pointer to a file containing a list * of bad regions in the image, read the file, one line * at a time, and create a REGION_STRUCT for each one. * Place the REGION_STRUCTs into a linked list, and * return a pointer to the linked list. * * RETURNS: * linked list of REGION_STRUCTs */ static REGION_STRUCT * read_regions ( FILE *maskfp /* I: file containing region information */ ) { char line[LINELEN + 1]; int id, numpix, min_row, max_row, min_col, max_col; REGION_STRUCT *new; REGION_STRUCT *head = NULL, *tail = NULL; if (maskfp == NULL) { return(NULL); } while (fgets(line, LINELEN, maskfp) != NULL) { if (sscanf(line, "%d %d %d %d %d %d", &id, &numpix, &min_row, &max_row, &min_col, &max_col) != 6) { error(0, "phot: read_regions: skipping bad input line"); continue; } if ((new = (REGION_STRUCT *) malloc(sizeof(REGION_STRUCT))) == NULL) { error(1, "phot: read_regions: can't malloc for region "); } new->region_id = id; new->numpix = numpix; new->min_row = min_row; new->max_row = max_row; new->min_col = min_col; new->max_col = max_col; new->next = NULL; #ifdef DEBUG printf(" reg %6d %5d pix %5d %5d %5d %5d \n", id, numpix, min_row, max_row, min_col, max_col); #endif if (tail == NULL) { head = new; tail = new; } else { tail->next = new; tail = tail->next; } } return(head); } /*************************************************************************** * PROCEDURE: free_regions * * DESCRIPTION: de-allocate memory for all the REGION_STRUCTs in the * given linked list. * * RETURNS: * nothing * */ static void free_regions ( REGION_STRUCT *head /* free all structs on this linked list */ ) { REGION_STRUCT *removed; while (head != NULL) { removed = head; head = head->next; free(removed); } return; } /*************************************************************************** * PROCEDURE: is_close_to_edge * * DESCRIPTION: given the location of a particular star, and given the * largest aperture radius, determine if the star lies * within this radius of any image edge. If so, return 1; * otherwise, return 0. * * RETURNS: * 0 if star is not close to an edge * 1 if it is */ static int is_close_to_edge ( double row, /* I: row centroid of star */ double col, /* I: col centroid of star */ int largest_aper_radius /* I: critical distance from edge */ ) { if (row < (double) largest_aper_radius) { return(1); } if (row > (double) (nrow - largest_aper_radius)) { return(1); } if (col < (double) largest_aper_radius) { return(1); } if (col > (double) (ncol - largest_aper_radius)) { return(1); } /* no, this star is far from any edge */ return(0); } /*************************************************************************** * PROCEDURE: is_close_to_bad_region * * DESCRIPTION: given the location of a particular star, and given the * largest aperture radius, determine if the star lies * within this radius of any bad region in a linked list * of regions. * * If it is close to a bad region, return 1. Otherwise, * return 0. * * RETURNS: * 0 if star is not close to a bad region * 1 if it is */ static int is_close_to_bad_region ( double row, /* I: row centroid of star */ double col, /* I: col centroid of star */ int largest_aper_radius /* I: critical distance from bad_region */ ) { REGION_STRUCT *reg; reg = region_list; while (reg != NULL) { if (row > (double) (reg->min_row - largest_aper_radius)) { if (row < (double) (reg->max_row + largest_aper_radius)) { if (col > (double) (reg->min_col - largest_aper_radius)) { if (col < (double) (reg->max_col + largest_aper_radius)) { return(1); } } } } reg = reg->next; } /* no, this star is far from any bad region */ return(0); }