/************************************************************************** * * * 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. * * * **************************************************************************/ /* more routines which pertain to image windows 10/22/96 - added new function "cleanup_profile" to take care of temporary files left over after plotting 4/18/99 - fixed bug in "do_zoom" which calculated the zoom boundaries incorrectly. 7/16/2002 - change 'find_sky' so it calculates interquartile mean of sky pixels 2/23/2003 - change name of "find_sky" routine to "image_find_sky" so that it doesn't conflict with a similar library routine */ #include #include #include #include #include "pcvista.h" #include "pcvistax.h" #include "pcvimage.h" #include "xplotlib.h" #include "fits.h" #undef DEBUG #undef DEBUG2 /* lots of messages from "sum_counts" */ #define CROSSHAIR_SIZE 10 /* the size of a crosshair, in pixels */ typedef int (*PFI)(); extern Display *display; extern int screen; extern Window window; extern GC gc, copygc, rubbergc; extern char *progname; /* * this variable starts out at 0. If the user tries to make a * radial-profile plot with the 'r' key, then we set this to 1 * to note the fact. We check its value to see if its necessary * to call "xplot_init()" to start a plot, or to call "xplot_quit()" * when we are about to quit the TV process. */ static int made_xplot = 0; static int16 ** get_prof_data(char *fname, int sr, int sc, int er, int ec); static void free_prof_data(int16 **p, int nr); static double image_find_sky(int16 **p, int nrow, int ncol, double innersky, double outersky); static int compar(int16 *a, int16 *b); static int sum_counts(int16 **p, int nrow, int ncol, double cr, double cc, double radius, double sky, double *summed_counts, double *summed_area); static double frac_pixel_inside_aperture(double aper_row_center, double aper_col_center, double aper_radius, int pixel_row, int pixel_col); /* zoom in on an image, forking off a new tv process that will focus in on the current position of the pointer, but at a zoom factor higher by a factor of ZOOM_FACTOR (if keystr == ZOOM_IN_STR) or lower by a factor of ZOOM_FACTOR (if keystr == ZOOM_OUT_STR). */ void do_zoom(event, key) XEvent *event; KeySym *key; { int r, c, sr, sc, er, ec, bin, zoom, pid, nr, nc, zero, span, dither; int box_sr, box_sc, box_nr, box_nc; int16 pixval; Window root_win, cur_win; int root_x, root_y, cur_x, cur_y, xo, xe, yo, ye; unsigned int keys_buttons; char *fname, box_buf[NBUF], buf[15][NBUF]; double zf; /* if False, pointer is not on the same screen as our image window */ if (XQueryPointer(display, window, &root_win, &cur_win, &root_x, &root_y, &cur_x, &cur_y, &keys_buttons) == False) { return; } if (get_data(cur_y, cur_x, &r, &c, &pixval, 1, &fname, &bin, &zoom) != 1) { return; } /* get all the particulars of this window */ if (get_win_data(&fname, &xo, &xe, &yo, &ye, &zoom, &sr, &er, &sc, &ec, &bin, &zero, &span, &dither) != 0) { return; } nr = (er - sr) + 1; nc = (ec - sc) + 1; /* now, figure out how to zoom in; try to make the new window roughly ZOOM_SIZE by ZOOM_SIZE screen pixels, but zoom in (or out) by a factor of ZOOM_FACTOR. define a box of the correct size in the symbol table. try to make it centered on the cursor position, but allow it not to be if near the edge. */ if (*key == ZOOM_IN_KEY) zf = (double) zoom*ZOOM_FACTOR/bin; else zf = (double) zoom/(ZOOM_FACTOR*bin); if (zf > 0.5) zf = (double) (floor((zf + SMALL) + 0.5)); else zf = (double) 1.0/(floor((1.0/zf) + SMALL + 0.5)); box_nr = (int) ZOOM_SIZE/zf; box_nc = (int) ZOOM_SIZE/zf; if (box_nr > nr) box_nr = nr; if (box_nc > nc) box_nc = nc; if ((box_sr = r - (box_nr/2)) < 0) { box_sr = 0; } if ((box_sc = c - (box_nc/2)) < 0) { box_sc = 0; } if (box_sr + box_nr > er) box_nr = er - (box_sr + 1); if (box_sc + box_nc > ec) box_nc = ec - (box_sc + 1); sprintf(box_buf, "box%d=%d %d %d %d", ZOOM_BOX, box_sr, box_sc, box_nr, box_nc); putsym(box_buf); /* and try to move the new image slightly closer to the center of the screen than the old one. */ if (xo > DisplayWidth(display, screen)/2) xo -= ZOOM_OFFSET; else xo += ZOOM_OFFSET; if (yo > DisplayHeight(display, screen)/2) yo -= ZOOM_OFFSET; else yo += ZOOM_OFFSET; strcpy(buf[0], progname); strcpy(buf[1], progname); /* must duplicate zero'th arg */ strcpy(buf[2], fname); sprintf(buf[3], "z=%d", zero); sprintf(buf[4], "l=%d", span); sprintf(buf[5], "zoom=%f", zf); sprintf(buf[6], "box=%d", ZOOM_BOX); sprintf(buf[7], "xo=%d", xo); sprintf(buf[8], "yo=%d", yo); if (dither) { sprintf(buf[9], "dither"); } if ((pid = fork()) != 0) { return; } signal(SIGINT, SIG_IGN); signal(SIGQUIT, SIG_IGN); signal(SIGTERM, SIG_IGN); if (dither) { execlp(buf[0], buf[1], buf[2], buf[3], buf[4], buf[5], buf[6], buf[7], buf[8], buf[9], (char *) 0); } else { execlp(buf[0], buf[1], buf[2], buf[3], buf[4], buf[5], buf[6], buf[7], buf[8], (char *) 0); } } /* move the cursor by a single screen pixel because the user typed one of the cursor keys */ void do_cursor_key(event, key) XEvent *event; KeySym *key; { int r, c, bin, zoom; int inc_x, inc_y; int16 pixval; Window root_win, cur_win; int root_x, root_y, cur_x, cur_y; unsigned int keys_buttons; char *fname; /* if False, pointer is not on the same screen as our image window */ if (XQueryPointer(display, window, &root_win, &cur_win, &root_x, &root_y, &cur_x, &cur_y, &keys_buttons) == False) { return; } if (get_data(cur_y, cur_x, &r, &c, &pixval, 1, &fname, &bin, &zoom) != 1) { return; } inc_x = 0; inc_y = 0; switch (*key) { case LEFT_KEY: cur_x--; inc_x--; break; case RIGHT_KEY: cur_x++; inc_x++; break; case UP_KEY: cur_y--; inc_y--; break; case DOWN_KEY: cur_y++; inc_y++; break; default: break; } /* check to see if the new position would be out of the image - if so, just do nothing and return. */ if (get_data(cur_y, cur_x, &r, &c, &pixval, 1, &fname, &bin, &zoom) != 1) { return; } /* since the new position would be valid, move the pointer there. This will generate a MotionNotify Event for the window, and it will then take appropriate action (printing the new position and value) */ XWarpPointer(display, RootWindow(display, screen), None, 0, 0, 0, 0, inc_x, inc_y); } /* handle getting a "crosshair" state: all this means is "draw a crosshair at some position." the center of the crosshair should be placed at (cr, cc) = (data[2], data[3]). if the given position doesn't fit within this window, just do nothing. */ void do_crosshair() { int x1, y1, cr, cc; int data[PROP_NITEMS]; /* get the information on the box to be drawn from the RootWindow property */ read_property(data); cr = data[2]; /* center row */ cc = data[3]; /* center col */ if (xform_in_window(0, cr, cc, &x1, &y1) != 1) { return; } /* * originally used "rubbergc" here, but that means that two * crosshairs that overlap will erase each other. Try using * "copygc", which ought to always appear as white (or black, * if the user asks for reverse video). */ XDrawLine(display, window, copygc, x1, y1, x1 + CROSSHAIR_SIZE, y1); XDrawLine(display, window, copygc, x1, y1, x1 - CROSSHAIR_SIZE, y1); XDrawLine(display, window, copygc, x1, y1, x1, y1 + CROSSHAIR_SIZE); XDrawLine(display, window, copygc, x1, y1, x1, y1 - CROSSHAIR_SIZE); } /* * create a radial profile for a subsection of the image, centered * near the position of the cursor. We spawn an "xplot" process * to make the plot. */ void do_profile_key(event, key) XEvent *event; KeySym *key; { int r, c, sr, sc, er, ec, bin, zoom, zero, span, dither; int16 pixval; Window root_win, cur_win; int root_x, root_y, cur_x, cur_y, xo, xe, yo, ye; unsigned int keys_buttons; char *fname, *symp, buf[50]; int boxrad; double innersky, outersky; double fitrad; int row, col, start_row, start_col, end_row, end_col, nrow, ncol; int16 **p; double row_cen, col_cen, fwhm, peak, eccen, major_axis, minor_axis; double row_width, col_width; FILE *fp; double rdist, cdist, dist, z, fitval; double sky; char command[100]; /* if False, pointer is not on the same screen as our image window */ if (XQueryPointer(display, window, &root_win, &cur_win, &root_x, &root_y, &cur_x, &cur_y, &keys_buttons) == False) { return; } #ifdef DEBUG printf("cur_x is %d, cur_y is %d\n", cur_x, cur_y); #endif if (get_data(cur_y, cur_x, &r, &c, &pixval, 1, &fname, &bin, &zoom) != 1) { #ifdef DEBUG printf("get_data returns non-1, so quit\n"); #endif return; } /* get all the particulars of this window */ if (get_win_data(&fname, &xo, &xe, &yo, &ye, &zoom, &sr, &er, &sc, &ec, &bin, &zero, &span, &dither) != 0) { return; } #ifdef DEBUG printf("got past get_win_data, we have sr %d sc %d er %d ec %d\n", sr, sc, er, ec); #endif /* * look for "profile_boxrad" in the symbol table -- if present, * use it as the "radius" of the square box we'll extract. * If not, use PROFILE_BOXRAD. */ boxrad = PROFILE_BOXRAD; if ((symp = getsym("profile_boxrad")) != NULL) { if ((sscanf(symp, "%d", &boxrad) != 1) || (boxrad < 1)) { boxrad = PROFILE_BOXRAD; } } else { sprintf(buf, "profile_boxrad=%-d", boxrad); putsym(buf); } /* * look for "profile_fitrad" in the symbol table -- if present, * use it as the radius for fitting a star's width. * If not, use PROFILE_FITRAD. */ fitrad = PROFILE_FITRAD; if ((symp = getsym("profile_fitrad")) != NULL) { if ((sscanf(symp, "%lf", &fitrad) != 1) || (fitrad <= 0.0)) { fitrad = PROFILE_FITRAD; } } else { sprintf(buf, "profile_fitrad=%-.2f", fitrad); putsym(buf); } /* * decide on the position and size of a box centered near the cursor */ if ((start_row = r - boxrad) < sr) { start_row = sr; } if ((end_row = r + boxrad) >= er) { end_row = er; } if ((start_col = c - boxrad) < sc) { start_col = sc; } if ((end_col = c + boxrad) >= ec) { end_col = ec; } nrow = 1 + end_row - start_row; ncol = 1 + end_col - start_col; /* * get the data inside that box */ if ((p = get_prof_data(fname, start_row, start_col, end_row, end_col)) == NULL) { error(0, "do_profile_key: can't get box data"); return; } #ifdef DEBUG printf("got past get_prof_data\n"); #endif /* * figure out a "sky" value -- just a rough one, based on outermost * edges of the box */ innersky = boxrad - 2.0; outersky = 2.0*boxrad; sky = image_find_sky(p, nrow, ncol, innersky, outersky); /* * calculate the centroid of pixel values inside this box */ if (do_axes(p, 0, 0, nrow, ncol, sky, fitrad, &row_cen, &col_cen, &row_width, &col_width, &fwhm, &peak, &eccen, &major_axis, &minor_axis) != 0) { free_prof_data(p, nrow); error(0, "do_profile_key: do_axes returns with error"); return; } row_cen += start_row; col_cen += start_col; #ifdef DEBUG printf("got past do_axes\n"); #endif /* * create a file which will temporarily hold one command * for the XPLOT program, telling it to wait for the * _real_ command file to appear. Spawn an XPLOT * process, with the temporary file as arg. */ if ((fp = fopen(PROFILE_CMDFILE, "w+")) == NULL) { free_prof_data(p, nrow); error(0, "do_profile_key: can't open command file for writing"); return; } fprintf(fp, "wait %s\n", PROFILE_CMDFILE); fclose(fp); if (made_xplot == 0) { xplot_init(); made_xplot = 1; } /* * create a temporary data file, filling it with 2-column lines, * one line per pixel in the box: * * distance_from_center pixel_value * */ if ((fp = fopen(PROFILE_DATAFILE, "w+")) == NULL) { free_prof_data(p, nrow); error(0, "do_profile_key: can't open data file for writing"); return; } for (row = 0; row < nrow; row++) { rdist = (row + start_row) - row_cen; for (col = 0; col < ncol; col++) { cdist = (col + start_col) - col_cen; dist = sqrt(rdist*rdist + cdist*cdist); fprintf(fp, "%f %f\n", dist, (p[row][col] - sky)); } } fclose(fp); #ifdef DEBUG printf("finished writing data file\n"); #endif /* * create a temporary data file, filling it with 2-column lines, * which show the best-fitting gaussian as a function of radius * * distance_from_center fitted_profile_value * */ if ((fp = fopen(PROFILE_FITFILE, "w+")) == NULL) { free_prof_data(p, nrow); error(0, "do_profile_key: can't open fitted file for writing"); return; } z = (-2.35*2.35)/(2.0*fwhm*fwhm); for (dist = 0.0; dist < boxrad; dist += PROFILE_DELTAX) { fitval = peak*exp(dist*dist*z); fprintf(fp, "%f %f\n", dist, fitval); } fclose(fp); #ifdef DEBUG printf("finished writing fitted file\n"); #endif /* * run a set of XPLOT commands to make a radial-profile plot. */ xplot_command("erase"); sprintf(command, "data %s", PROFILE_DATAFILE); xplot_command(command); xplot_command("xc 1"); xplot_command("yc 2"); xplot_command("limits"); xplot_command("box"); sprintf(command, "title (%7.2f %7.2f) FWHM %-5.2f ec %-5.2f PA %-5.1f sky %-7.1f\n", row_cen, col_cen, fwhm, eccen, major_axis, sky); xplot_command(command); xplot_command("expand 0.5"); xplot_command("points x y"); xplot_command("expand 1.0"); sprintf(command, "data %s", PROFILE_FITFILE); xplot_command(command); xplot_command("ltype 1"); xplot_command("xc 1"); xplot_command("yc 2"); xplot_command("connect x y"); xplot_command("ltype 0"); /* * we're all done, so we can free up the little chunk of memory * we allocated to make the profile */ free_prof_data(p, nrow); } /* * Check to see if we have made any profile plots during the * lifetime of this process. If so, call "xplot_quit()" to * shut down the plotting window nicely. * * If we _don't_ call "xplot_quit()", we have a leftover * file in the /tmp directory! * * In addition, delete the temporary files we used for * holding the XPLOT commands and data (which live in the * current directory). */ void cleanup_profile(void) { int ret; if (made_xplot == 1) { xplot_quit(); ret = unlink(PROFILE_DATAFILE); ret = unlink(PROFILE_FITFILE); ret = unlink(PROFILE_CMDFILE); } } /* * allocate an area in memory big enough to hold the data we'll * use to calculate radial profile, and return a pointer to it. * Return NULL if we fail. */ static int16 ** get_prof_data(char *fname, int sr, int sc, int er, int ec) { FITS_HANDLE fh; int16 **p, **pp, *q; int nr, nc, nbytes; int row; fh = fits_open(fname, "r", &nr, &nc); nr = 1 + (er - sr); nc = 1 + (ec - sc); nbytes = nr*sizeof(int16 *); if ((p = (int16 **) malloc(nbytes)) == NULL) { error(0, "get_prof_data: can't alloc for pointers"); return(NULL); } pp = p; nbytes = nc*sizeof(int16); for (row = sr; row <= er; row++) { if ((q = (int16 *) malloc(nbytes)) == NULL) { fprintf(stderr, "get_prof_data: can't allocate for data array\n"); return(NULL); } *pp++ = q; fits_get_data(fh, row, sc, q, nc); } fits_close(fh); return(p); } /* * free the data allocated by "get_prof_data" */ static void free_prof_data(int16 **p, int nr) { int i; for (i = 0; i < nr; i++) { free(p[i]); } free(p); } /* * perform aperture photometry for a subsection of the image, centered * near the position of the cursor. */ void do_aperture_key(event, key) XEvent *event; KeySym *key; { int r, c, sr, sc, er, ec, bin, zoom, zero, span, dither; int16 pixval; Window root_win, cur_win; int root_x, root_y, cur_x, cur_y, xo, xe, yo, ye; unsigned int keys_buttons; char *fname, *symp, buf[50]; char errmsg[200]; int boxrad; double fitrad; int start_row, start_col, end_row, end_col, nrow, ncol; int16 **p; double row_cen, col_cen, fwhm, peak, eccen, major_axis, minor_axis; double row_width, col_width; double radius, magzero, sum, mag, area; double innersky, outersky, sky; /* if False, pointer is not on the same screen as our image window */ if (XQueryPointer(display, window, &root_win, &cur_win, &root_x, &root_y, &cur_x, &cur_y, &keys_buttons) == False) { return; } #ifdef DEBUG printf("cur_x is %d, cur_y is %d\n", cur_x, cur_y); #endif if (get_data(cur_y, cur_x, &r, &c, &pixval, 1, &fname, &bin, &zoom) != 1) { #ifdef DEBUG printf("get_data returns non-1, so quit\n"); #endif return; } /* get all the particulars of this window */ if (get_win_data(&fname, &xo, &xe, &yo, &ye, &zoom, &sr, &er, &sc, &ec, &bin, &zero, &span, &dither) != 0) { return; } #ifdef DEBUG printf("got past get_win_data, we have sr %d sc %d er %d ec %d\n", sr, sc, er, ec); #endif /* * look for "aperture_outersky" in the symbol table -- if present, * use it as the "radius" of the square box we'll extract. * If not, use APERTURE_OUTERSKY. */ outersky = APERTURE_OUTERSKY; if ((symp = getsym("aperture_outersky")) != NULL) { if ((sscanf(symp, "%lf", &outersky) != 1) || (outersky < 1)) { outersky = APERTURE_OUTERSKY; } } else { sprintf(buf, "aperture_outersky=%-.2f", outersky); putsym(buf); } boxrad = (int) (outersky + 0.5); /* * look for "aperture_innersky" in the symbol table -- if present, * set "innersky" to it. * If not, use APERTURE_INNERSKY. */ innersky = APERTURE_INNERSKY; if ((symp = getsym("aperture_innersky")) != NULL) { if ((sscanf(symp, "%lf", &innersky) != 1) || (innersky <= 0.0)) { innersky = APERTURE_INNERSKY; } } else { sprintf(buf, "aperture_innersky=%-.2f", innersky); putsym(buf); } /* * look for "aperture_radius" in the symbol table -- if present, * set "radius" to it. * If not, use APERTURE_RADIUS. */ radius = APERTURE_RADIUS; if ((symp = getsym("aperture_radius")) != NULL) { if ((sscanf(symp, "%lf", &radius) != 1) || (innersky <= 0.0)) { innersky = APERTURE_RADIUS; } } else { sprintf(buf, "aperture_radius=%-.2f", radius); putsym(buf); } /* * look for "aperture_magzero" in the symbol table -- if present, * set "magzero" to it. * If not, use APERTURE_MAGZERO. */ magzero = APERTURE_MAGZERO; if ((symp = getsym("aperture_magzero")) != NULL) { if (sscanf(symp, "%lf", &magzero) != 1) { magzero = APERTURE_MAGZERO; } } else { sprintf(buf, "aperture_magzero=%-.2f", magzero); putsym(buf); } /* * decide on the position and size of a box centered near the cursor */ if ((start_row = r - boxrad) < sr) { start_row = sr; } if ((end_row = r + boxrad) >= er) { end_row = er; } if ((start_col = c - boxrad) < sc) { start_col = sc; } if ((end_col = c + boxrad) >= ec) { end_col = ec; } nrow = 1 + end_row - start_row; ncol = 1 + end_col - start_col; /* * get the data inside that box */ #ifdef DEBUG printf("about to call get_prof_data with sr %4d sc %4d er %4d ec %4d \n", start_row, start_col, end_row, end_col); #endif if ((p = get_prof_data(fname, start_row, start_col, end_row, end_col)) == NULL) { error(0, "do_aperture_key: can't get box data"); return; } #ifdef DEBUG printf("got past get_prof_data\n"); #endif /* * figure out a "sky" value, using pixels at radii between * innersky and outersky */ sky = image_find_sky(p, nrow, ncol, innersky, outersky); /* * calculate the centroid of pixel values inside this box * If we can't calculate the centroid, just use the position * of the cursor as the centroid, and warn the user! */ #if 1 fitrad = radius; #else fitrad = 5.0; #endif if (do_axes(p, 0, 0, nrow, ncol, sky, fitrad, &row_cen, &col_cen, &row_width, &col_width, &fwhm, &peak, &eccen, &major_axis, &minor_axis) != 0) { error(0, "do_aperture_key: do_axes returns with error"); /* * so, use position of cursor as centroid. We have to add and * subtract some stuff to make it look like we just got * values "row_cen" and "col_cen" back from do_axes, which * calculated their position with respect to the little * box of pixels "p". */ row_cen = (r + 0.5) - start_row; col_cen = (c + 0.5) - start_col; sprintf(errmsg, "do_aperture_key: assume centroid at %6.1f,%6.1f, as given by cursor", row_cen + start_row, col_cen + start_col); error(0, errmsg); #if 0 /* we would do this if we didn't continue at cursor position */ free_prof_data(p, nrow); return; #endif } row_cen += start_row; col_cen += start_col; #ifdef DEBUG printf("got past do_axes\n"); #endif /* * now add up the flux inside a circular aperture of radius "radius" * pixels, centered on the star. */ if (sum_counts(p, nrow, ncol, row_cen - start_row, col_cen - start_col, radius, sky, &sum, &area) != 0) { free_prof_data(p, nrow); error(0, "do_aperture_key: sum_counts returns with error"); return; } /* * convert the result to magnitudes, and print it out */ if (sum > 0.0) { mag = magzero - 2.5*log10(sum); } else { mag = magzero; } printf("(%7.2f %7.2f) flux %8.1f npix %5.1f mag %6.3f sky %6.1f [%4.1f %2.0f %2.0f]\n", row_cen, col_cen, sum, area, mag, sky, radius, innersky, outersky); /* * we're all done, so we can free up the little chunk of memory * we allocated to make the profile */ free_prof_data(p, nrow); } /* * do a simple curve-of-growth analysis on the star under the cursor. * This is just a copy of the "do_aperture_key" function above, * except that it adds the flux within several different * circular apertures at the end. */ void do_curveofgrowth_key(event, key) XEvent *event; KeySym *key; { int r, c, sr, sc, er, ec, bin, zoom, zero, span, dither; int16 pixval; Window root_win, cur_win; int root_x, root_y, cur_x, cur_y, xo, xe, yo, ye; unsigned int keys_buttons; char *fname, *symp, buf[50]; char errmsg[200]; int boxrad; double fitrad; int start_row, start_col, end_row, end_col, nrow, ncol; int i, npoints; int16 **p; double row_cen, col_cen, fwhm, peak, eccen, major_axis, minor_axis; double row_width, col_width; double sum, area; double innersky, outersky, sky; double curve_radii[CURVEOFG_MAXPOINTS]; double curve_flux[CURVEOFG_MAXPOINTS]; double curve_value, curve_max; /* if False, pointer is not on the same screen as our image window */ if (XQueryPointer(display, window, &root_win, &cur_win, &root_x, &root_y, &cur_x, &cur_y, &keys_buttons) == False) { return; } #ifdef DEBUG printf("cur_x is %d, cur_y is %d\n", cur_x, cur_y); #endif if (get_data(cur_y, cur_x, &r, &c, &pixval, 1, &fname, &bin, &zoom) != 1) { #ifdef DEBUG printf("get_data returns non-1, so quit\n"); #endif return; } /* get all the particulars of this window */ if (get_win_data(&fname, &xo, &xe, &yo, &ye, &zoom, &sr, &er, &sc, &ec, &bin, &zero, &span, &dither) != 0) { return; } #ifdef DEBUG printf("got past get_win_data, we have sr %d sc %d er %d ec %d\n", sr, sc, er, ec); #endif /* * look for "aperture_outersky" in the symbol table -- if present, * use it as the "radius" of the square box we'll extract. * If not, use APERTURE_OUTERSKY. */ outersky = APERTURE_OUTERSKY; if ((symp = getsym("aperture_outersky")) != NULL) { if ((sscanf(symp, "%lf", &outersky) != 1) || (outersky < 1)) { outersky = APERTURE_OUTERSKY; } } else { sprintf(buf, "aperture_outersky=%-.2f", outersky); putsym(buf); } boxrad = (int) (outersky + 0.5); /* * look for "aperture_innersky" in the symbol table -- if present, * set "innersky" to it. * If not, use APERTURE_INNERSKY. */ innersky = APERTURE_INNERSKY; if ((symp = getsym("aperture_innersky")) != NULL) { if ((sscanf(symp, "%lf", &innersky) != 1) || (innersky <= 0.0)) { innersky = APERTURE_INNERSKY; } } else { sprintf(buf, "aperture_innersky=%-.2f", innersky); putsym(buf); } /* * look for "curveofg_npoints" in the symbol table -- if present, * set "npoints" to it. * If not, use CURVEOFG_NPOINTS */ npoints = CURVEOFG_NPOINTS; if ((symp = getsym("curveofg_npoints")) != NULL) { if ((sscanf(symp, "%d", &npoints) != 1) || (npoints <= 0)) { npoints = CURVEOFG_NPOINTS; } } else { #if 0 /* I don't think I want to do this */ sprintf(buf, "curveofg_npoints=%d", npoints); putsym(buf); #endif } #ifdef DEBUG printf("do_curveofgrowth_key: npoints = %d\n", npoints); #endif /* * sanity check */ if (npoints >= CURVEOFG_MAXPOINTS) { error(0, "do_curveofgrowth_key: too many points in curve"); return; } /* * look for the radii of each aperture to be used in curve of growth. * If we can't find it in the symbol table, use default values * for the first three points -- then stop. */ for (i = 0; i < npoints; i++) { sprintf(buf, "curveofg_r%d", i+1); if (((symp = getsym(buf)) == NULL) || (sscanf(symp, "%lf", &(curve_radii[i])) != 1) || (curve_radii[i] <= 0.0)) { #ifdef DEBUG printf("do_curveofgrowth_key: bad radius %d value %f\n", i, curve_radii[i]); #endif if (i == 0) { curve_radii[i] = CURVEOFG_R1; } else if (i == 1) { curve_radii[i] = CURVEOFG_R2; } else if (i == 2) { curve_radii[i] = CURVEOFG_R3; } else { /* make all other un-specified radii = CURVE_R3 */ curve_radii[i] = CURVEOFG_R3; } } #ifdef DEBUG printf("do_curveofgrowth_key: radius %d value %f\n", i, curve_radii[i]); #endif #if 0 /* I don't think I want to do this ... */ sprintf(buf, "curveofg_r%d=%-.1f", i, curve_radii[i]); putsym(buf); #endif } /* * decide on the position and size of a box centered near the cursor */ if ((start_row = r - boxrad) < sr) { start_row = sr; } if ((end_row = r + boxrad) >= er) { end_row = er; } if ((start_col = c - boxrad) < sc) { start_col = sc; } if ((end_col = c + boxrad) >= ec) { end_col = ec; } nrow = 1 + end_row - start_row; ncol = 1 + end_col - start_col; /* * get the data inside that box */ if ((p = get_prof_data(fname, start_row, start_col, end_row, end_col)) == NULL) { error(0, "do_aperture_key: can't get box data"); return; } #ifdef DEBUG printf("got past get_prof_data\n"); #endif /* * figure out a "sky" value, using pixels at radii between * innersky and outersky */ sky = image_find_sky(p, nrow, ncol, innersky, outersky); /* * calculate the centroid of pixel values inside this box * If we can't calculate the centroid, just use the position * of the cursor as the centroid, and warn the user! */ #if 0 fitrad = 2.0*radius; #else fitrad = 5.0; #endif if (do_axes(p, 0, 0, nrow, ncol, sky, fitrad, &row_cen, &col_cen, &row_width, &col_width, &fwhm, &peak, &eccen, &major_axis, &minor_axis) != 0) { error(0, "do_aperture_key: do_axes returns with error"); /* * so, use position of cursor as centroid. We have to add and * subtract some stuff to make it look like we just got * values "row_cen" and "col_cen" back from do_axes, which * calculated their position with respect to the little * box of pixels "p". */ row_cen = (r + 0.5) - start_row; col_cen = (c + 0.5) - start_col; sprintf(errmsg, "do_aperture_key: assume centroid at %6.1f,%6.1f, as given by cursor", row_cen + start_row, col_cen + start_col); error(0, errmsg); #if 0 /* we would do this if we didn't continue at cursor position */ free_prof_data(p, nrow); return; #endif } row_cen += start_row; col_cen += start_col; #ifdef DEBUG printf("got past do_axes\n"); #endif /* * Okay, we're ready to calculate the curve of growth. * What we do is to add up the flux inside all the apertures, * and then divide each one by the flux of the outermost aperture. * * now add up the flux inside a circular aperture of radius "radius" * pixels, centered on the star. */ for (i = 0; i < npoints; i++) { if (sum_counts(p, nrow, ncol, row_cen - start_row, col_cen - start_col, curve_radii[i], sky, &sum, &area) != 0) { free_prof_data(p, nrow); error(0, "do_curveofgrowth_key: sum_counts returns with error"); return; } curve_flux[i] = sum; } /* * convert each flux to a fraction of flux in outermost aperture, * and print out the list of values in a nice format. */ printf("%7.2f %7.2f ", row_cen, col_cen); curve_max = curve_flux[npoints - 1]; for (i = 0; i < npoints; i++) { if (curve_max > 0.0) { curve_value = curve_flux[i] / curve_max; } printf("%5.3f ", curve_value); } printf("\n"); /* * we're all done, so we can free up the little chunk of memory * we allocated to make the profile */ free_prof_data(p, nrow); } /**************************************************************************** * figure out a rough sky value, based on the pixels in the given box * which are between "innersky" and "outersky" pixels from the * center of the box. * * We use a median of pixels in the annulus to find the sky value. * * return the value of the sky we used. */ static double image_find_sky(int16 **p, int nrow, int ncol, double innersky, double outersky) { int16 *array; int narray; int col, row; int cr, cc; long sum, npix; long dr2, dc2, d2; double inner2, outer2; double sky; /* allocate space for the array of annulus pixel values */ narray = nrow*ncol; if ((array = (int16 *) malloc(narray*sizeof(int16))) == NULL) { error(1, "image_find_sky: can't allocate for array"); } /* square the inner and outer radii, for quicker comparison */ inner2 = innersky*innersky; outer2 = outersky*outersky; /* we calculate distances relative to center of box */ cr = nrow/2; cc = ncol/2; sky = 0.0; sum = 0; npix = 0; for (row = 0; row < nrow; row++) { dr2 = (cr - row)*(cr - row); if (dr2 > outer2) { continue; } for (col = 0; col < ncol; col++) { dc2 = (cc - col)*(cc - col); if (dc2 > outer2) { continue; } d2 = dr2 + dc2; if ((d2 < inner2) || (d2 > outer2)) { continue; } sum += p[row][col]; array[npix] = p[row][col]; npix++; } } #ifdef MEAN /* we calculate regular mean sky value here .... */ if (npix > 0) { sky = sum/(double) npix; } else { sky = 0.0; } #else /* or calculate interquartile mean sky value here ... */ if (npix > 4) { qsort((char *) array, npix, sizeof(int16), (PFI) compar); sum = 0; for (row = (npix/4); row < (3*npix/4); row++) { sum += array[row]; } sky = sum/((3*npix/4) - (npix/4)); } else { sky = 0.0; } free(array); #endif return(sky); } /*************************************************************************** * return 1 if the first value is greater, -1 if second is greater, * or 0 if the two are equal */ static int compar(int16 *a, int16 *b) { if (*a > *b) { return(1); } else if (*a < *b) { return(-1); } else { return(0); } } /**************************************************************************** * Add up all the counts inside a circular aperture of "radius" pixels * around (cr, cc). Subtract an appropriate amount of sky from the * sum, and return the sum in the last argument. * * Note that we deal with pixels which fall on the edge of an aperture * by counting a fraction of the total flux ... and by adding the * same fraction of that pixel to the sky contribution. The rule is * * calculate distance "D" from center of pixel to center of aperture * * if (D < radius - 0.5) fraction = 1.0 * if (radius - 0.5 < D < radius + 0.5) fraction = (radius + 0.5) - D * if (D > radius + 0.5) fraction = 0.0 * * Place the sum of counts into the "summed_counts" argument, and * place the sum of pixel area (including fractional pixels) into the * "summed_area" argument. * * Return 0 if all goes well, * 1 if there is an error. * */ static int sum_counts(int16 **p, int nrow, int ncol, double cr, double cc, double radius, double sky, double *summed_counts, double *summed_area) { int col, row; double sum, pixval, npix; double dr2, dc2, d2; double inner2, outer2; double fraction; if ((cr < 0) || (cr > (nrow - 1))) { error(0, "sum_counts: central row out of bounds"); return(1); } if ((cc < 0) || (cc > (ncol - 1))) { error(0, "sum_counts: central col out of bounds"); return(1); } /* make some temporary variables, for quicker comparison */ outer2 = (radius + 0.5)*(radius + 0.5); inner2 = (radius - 0.5)*(radius - 0.5); #ifdef DEBUG2 printf(" inner2=%-7.2f outer2=%-7.2f \n", inner2, outer2); #endif sum = 0; npix = 0; for (row = 0; row < nrow; row++) { dr2 = (cr - row)*(cr - row); if (dr2 > outer2) { continue; } for (col = 0; col < ncol; col++) { dc2 = (cc - col)*(cc - col); if (dc2 > outer2) { continue; } d2 = dr2 + dc2; #ifdef DEBUG2 printf(" d2 is %7.2f ", d2); #endif if (d2 < inner2) { /* we add entire pixel */ npix += 1.0; pixval = p[row][col]; pixval -= sky; sum += pixval; #ifdef DEBUG2 printf("adding entire pix at (%4d,%4d)\n", row, col); #endif } else if ((d2 >= inner2) && (d2 <= outer2)) { /* we add a fraction of the pixel */ #if 0 fraction = (radius + 0.5) - sqrt(d2); #else fraction = frac_pixel_inside_aperture(cr, cc, radius, row, col); #endif npix += fraction; pixval = p[row][col]; pixval -= sky; sum += fraction*pixval; #ifdef DEBUG2 printf("adding frac %5.2f at (%4d,%4d)\n", fraction, row, col); #endif } else { #ifdef DEBUG2 printf(" \n"); #endif } } } *summed_counts = sum; *summed_area = npix; return(0); } /*************************************************************************** * 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); }