/*** File caphot.c *** January 30, 2002 *** By Jessica Mink from Fortran code by Sam Conner (MIT, 1984) *** Copyright (C) 2002 *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA Correspondence concerning WCSTools should be addressed as follows: Internet email: jmink@cfa.harvard.edu Postal address: Jessica Mink Smithsonian Astrophysical Observatory 60 Garden St. Cambridge, MA 02138 USA */ #include #include "fitshead.h" #include "imio.h" static double apint(); static double imapfr(); /* CAPHOT -- Perform rigorous circular aperture photometry using the imapfr() * subroutine to determine the fraction of a square pixel enclosed * by a circular aperture. */ void caphot (imbuff,header,cx,cy,rad,sumw,wsum) char *imbuff; /* address of start of image buffer */ char *header; /* image FITS header */ double cx,cy; /* center x,y of aperture */ double rad; /* radius of aperture */ double sumw; /* sum of values of pixel weights (returned) */ double wsum; /* sum of weighted pixel fluxes (returned) */ { double x, y, factor, flux, bs, bz; int ix1, ix2, iy1, iy2, ix, iy, bitpix, nx, ny; sumw = 0.0; wsum = 0.0; hgeti4 (header, "BITPIX", &bitpix); hgeti4 (header, "NAXIS1", &nx); hgeti4 (header, "NAXIS2", &ny); hgetr8 (header, "BSCALE", &bs); hgeti4 (header, "BZERO", &bz); /* Find range of ys to check */ iy1 = (int) (cy - rad); iy2 = (int) (cy + rad + 0.99999); if (iy1 < 1) iy1 = 1; if (iy2 > ny) iy2 = ny; /* Find range of xumns to check */ ix1 = (int) (cx - rad); ix2 = (int) (cx + rad + 0.99999); if (ix1 < 1) ix1 = 1; if (ix2 > nx) ix2 = nx; for (iy = iy1; iy <= iy2; iy++) { y = (double) iy; for (ix = ix1; ix <= ix2; ix++) { x = (double) ix; factor = imapfr (x,y,cx,cy,rad); sumw = sumw + factor; flux = getpix1 (imbuff, bitpix, nx, ny, bz, bs, ix, iy); wsum = wsum + (factor * flux); /* fprintf (stderr, "IMAPSB: (%d,%d)= %f weight= %f\n", ix,iy,flux,factor); */ } } /* fprintf (stderr, "IMAPSB: sum of weights = %f\n", sumw); fprintf (stderr, "IMAPSB: sum of weighted intensity = %f\n", wsum); */ return; } /* IMAPFR -- Determine the fraction of a square pixel that is included * within the boundary of a circle of arbitrary center and radius */ static double imapfr (pcol,prow,ccol,crow,rad) double pcol,prow; /* column, row of pixel */ double ccol,crow; /* column, row of center of circle */ double rad; /* radius of circle in pixels */ { int icorn[4], ncorn, j, i; double dist, xdiff, y0, y1, ydiff, x0, x1, x, y, dc5, dr5; double cterm, ulim, llim, dc, dr, dx, dy; double frac; /* fraction of pixel within circle */ /* Compute distance to center of circle from most distant corner of pixel */ dc = fabs (pcol - ccol); dr = fabs (prow - crow); dc5 = dc + 0.5; dr5 = dr + 0.5; dist = sqrt ((dc5 * dc5) + (dr5 * dr5)); dc5 = dc - 0.5; dr5 = dr - 0.5; /* if the pixel is completely enclosed by the aperture, return 1 */ if (dist <= rad) return (1.0); /* if the pixel is not completely included, compute distance from the center of the circle to the nearest point on the periphery of the pixel*/ if ((pcol-0.5) < ccol && ccol < (pcol+.5)) { dist = dr5; /* the pixel is completely excluded from the aperture */ if (dist >= rad) return (0.0); else if ((prow - 0.5) < crow && crow < (prow + 0.5)) { dist = dc5; /* the pixel is completely excluded from the aperture */ if (dist >= rad) return (0.0); } else dist = sqrt (dc5*dc5 + dr5*dr5); /* the pixel is completely excluded from the aperture */ if (dist >= rad) return (0.0); } /* If the pixel is partially included in the aperture, determine how many corners are enclosed */ ncorn = 1; y = prow - 1.5; for (i = 0; i < 2; i++) { y = y + 1.0; x = pcol - 1.50; for (j = 0; j < 2; j++) { icorn[j+(i*2)] = 0; x = x + 1.0; dx = x - ccol; dy = y - crow; dist = sqrt (dx*dx + dy*dy); if (dist < rad) { ncorn = ncorn + 1; icorn[j+(i*2)] = 1; } } } /* Depending on number of corners enclosed, branch to appropriate computation */ /* If no corners are enclosed, determine whether the slice is vertical or horizontal */ if (ncorn < 1) { /* If the slice is vertical (at the sides of the aperture, at extreme x), determine the limits of integration */ if ((pcol - 0.5) >= ccol || ccol >= (pcol + 0.5)) { xdiff = fabs (pcol - ccol) - 0.5; dx = sqrt (rad*rad - xdiff*xdiff); y0 = crow - dx; y1 = crow + dx; if (ccol >= pcol) x = pcol + 0.5; else x = pcol - 0.5; /* Evaluate integral */ frac = -xdiff * (y1 - y0) + apint ((y1 - crow), rad) - apint ((y0 - crow), rad); return (frac); } /* If the slice is horizontal (at top or bottom of the aperture, at extreme y), determine the limits of integration. */ else { ydiff = fabs (prow - crow) - 0.5; dy = sqrt (rad*rad - ydiff*ydiff); x0 = ccol - dy; x1 = ccol + dy; if (crow >= prow) y = prow + 0.5; else y = prow - 0.5; /* Evaluate integral */ frac = -ydiff*(x1 - x0) + apint ((x1 - ccol),rad) - apint ((x0 - ccol),rad); return (frac); } } /* if one corner is enclosed, find direction (ne,nw,se,sw) from the center*/ else if (ncorn < 2) { if (pcol < ccol) { /* Determine the limits of integration */ ydiff = fabs (prow - crow) - 0.5;; x0 = ccol - sqrt (rad*rad - ydiff*ydiff); if (crow >= prow) y = prow + 0.5; else y = prow - 0.5; /* Evaluate integral */ frac = -ydiff*(pcol + 0.5 - x0) + apint ((pcol + 0.5 - ccol),rad) - apint ((x0 - ccol),rad); return (frac); } else { /* Determine the limits of integration */ ydiff = fabs (crow - prow) - 0.5; x1 = ccol + sqrt (rad*rad - ydiff*ydiff); if (crow >= prow) y = prow + 0.5; else y = prow - 0.5; /* Evaluate integral */ frac = -ydiff*(x1 - pcol + 0.5) + apint ((x1 - ccol),rad) - apint ((pcol - .5 - ccol),rad); return (frac); } } /* If two corners are enclosed, this problem has two cases: one in which the aperture boundary intersects the pixel boundary twice, and another in which it intersects the pixel boundaries 4 times */ /* determine whether this may be a 4-intersection configuration */ else if (ncorn < 3) { if (((pcol - 0.5) < ccol && ccol < (pcol + .5) && (fabs(crow - prow) + 0.5) < rad) || ((prow - 0.5) < crow && crow < (prow + 0.5) && (fabs(ccol - pcol) + 0.5) < rad)) { /* if the pixel boundary intersects the aperture boundary in four places, determine whether the pixel is east-west or north-south. */ /* pixel is east-west */ if ((pcol - 0.5) >= ccol || (pcol + 0.5) <= ccol) { xdiff = dc + 0.5; /* Determine the limits of integration */ dx = sqrt (rad*rad - xdiff*xdiff); y0 = crow - dx; y1 = crow + dx; if (pcol >= ccol) x = pcol + .5; else x = pcol - .5; frac = 1. - xdiff*(y0 - y1 + 1.0) + apint ((y0 - crow), rad) - apint ((prow - 0.5 - crow), rad) + apint ((prow + 0.5 - crow), rad) - apint ((y1 - crow), rad); return (frac); } /* Pixel is north-south */ else { ydiff = dr + 0.5; /* determine the limits of integration */ dy = sqrt (rad*rad - ydiff*ydiff); /* x0 is the x-coordinate of the small x intercept */ x0 = ccol - dy; /* x1 is the x-coordinate of the large x intercept */ x1 = ccol + dy; if (prow >= crow) y = prow + 0.5; else y = prow - 0.5; frac = 1. - ydiff * (x0 - x1 + 1.0) + apint ((x0 - ccol), rad) - apint ((pcol - 0.5 - ccol), rad) + apint ((pcol + 0.5 - ccol), rad) - apint ((x1 - ccol), rad); return (frac); } } /* in the two-intersection case, determine which corners are included */ else { /* the pixel is north or south of the center of the aperture */ if (icorn[0]*icorn[1] == 1 || icorn[2]*icorn[3] == 1) { cterm = 0.5 - dr; ulim = apint ((pcol + 0.5 - ccol) ,rad); llim = apint ((pcol - 0.5 - ccol) ,rad); frac = cterm + ulim - llim; return (frac); } /* the pixel is east or west of the center */ else { cterm = 0.5 - dc; ulim = apint ((prow + 0.5 - crow), rad); llim = apint ((prow - 0.5 - crow), rad); frac = cterm + ulim - llim; return (frac); } } } /* if three corners are enclosed, determine whether the corner in question is east or west of the center of the aperture */ else { if (pcol > ccol) { ydiff = dr + 0.5; if (prow >= crow) y = prow + 0.5; else y = prow - 0.5; x0 = ccol + sqrt (rad*rad - ydiff*ydiff); frac = 1.0 - (ydiff * (pcol + 0.5 - x0)) + apint ((pcol + 0.5 - ccol), rad) - apint ((x0 - ccol), rad); return (frac); } else { ydiff = dr + 0.5; if (prow >= crow) y = prow + 0.5; else y = prow - 0.5; x0 = ccol - sqrt (rad*rad - ydiff*ydiff); frac = 1. - ydiff*(x0 - pcol + 0.5) + apint ((x0 - ccol),rad) - apint ((pcol - 0.5 - ccol),rad); return (frac); } } } /* APINT -- Evaluate the integral of sqrt (rad**2 - x**2) dx at one limit */ static double apint (x, rad) double x; double rad; { double x2,rad2,arg,arg2,arcsin,pi; pi = 3.141592654; arg = x / rad; x2 = x * x; rad2 = rad * rad; arg2 = x2 / rad2; arcsin = atan2 (arg, sqrt (1.0 - arg2)); if ((1. - fabs (arg)) < 0.000001) { if (arg >= 0) arcsin = pi / 2.0; else arcsin = -pi / 2.0; } return (0.5 * (x * sqrt (rad2 - x2) + rad2 * arcsin)); } /* Jul 4 1984 Original Fortran program by Sam Conner at MIT * Mar 2 1987 Original VAX Unix version * Mar 2 1987 declare undecldared variables x, y, and j * * Jan 8 1993 Modify to handle arbitrary byte-per-pixel images * Apr 16 1993 Declare undeclared variables iy, ic0, and ir0 * * Jan 30 2002 Translate from Fortran to C */