diff options
Diffstat (limited to 'noao/digiphot/daophot/psf/dppsubrast.x')
| -rw-r--r-- | noao/digiphot/daophot/psf/dppsubrast.x | 172 |
1 files changed, 172 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/psf/dppsubrast.x b/noao/digiphot/daophot/psf/dppsubrast.x new file mode 100644 index 00000000..1fafb408 --- /dev/null +++ b/noao/digiphot/daophot/psf/dppsubrast.x @@ -0,0 +1,172 @@ +include <mach.h> +include <imhdr.h> +include "../lib/daophotdef.h" +include "../lib/psfdef.h" + +# DP_PSUBRAST -- Fetch the prospective PSF star data and check it for bad +# pixel values. + +pointer procedure dp_psubrast (dao, im, lowbad, highbad, x1, x2, y1, y2, + saturated) + +pointer dao # pointer to the daophot structure +pointer im # pointer to the input image +real lowbad, highbad # minimum and maximum good data values +int x1, x2 # output x limits of the extracted subraster +int y1, y2 # output y limits of the extracted subraster +int saturated # is the star saturated ? + +pointer psf, buf +real psfrad, fitrad +int dp_chksr() +pointer dp_subrast() + +begin + # Initialize. + psf = DP_PSF(dao) + buf = NULL + psfrad = DP_PSFRAD(dao) + fitrad = DP_FITRAD(dao) + + # Get the data. + buf = dp_subrast (im, DP_CUR_PSFX(psf), DP_CUR_PSFY(psf), psfrad, + x1, x2, y1, y2) + if (buf == NULL) + return (NULL) + + # Check for bad pixels in subraster, compute the min and max. + if (dp_chksr (DP_CUR_PSFX(psf), DP_CUR_PSFY(psf), Memr[buf], + x2 - x1 + 1, y2 - y1 + 1, x1, y1, psfrad, fitrad, + lowbad, highbad, saturated, DP_CUR_PSFMIN(psf), + DP_CUR_PSFMAX(psf), DP_CUR_PSFGMAX(psf)) == ERR) { + call mfree (buf, TY_REAL) + return (NULL) + } + + return (buf) +end + + +# DL_LSUBRAST -- Give a valid PSF star compute the limits of the data to +# be extracted around it. + +int procedure dp_lsubrast (im, xcen, ycen, radius, x1, x2, y1, y2) + +pointer im # input image descriptor +real xcen, ycen # center of subraster +real radius # radius of the box +int x1, y1, x2, y2 # boundaries of subraster + +begin + # Calculate start position of extraction box. + x1 = int (xcen - radius) - 2 + x2 = int (xcen + radius) + 3 + y1 = int (ycen - radius) - 2 + y2 = int (ycen + radius) + 3 + if (x1 > IM_LEN(im,1) || x2 < 1 || y1 > IM_LEN(im,2) || y2 < 1) + return (ERR) + + x1 = max (1, x1) + x2 = min (IM_LEN(im,1), x2) + y1 = max (1, y1) + y2 = min (IM_LEN(im,2), y2) + return (OK) +end + + +# DP_SUBRAST -- Given a valid PSF star extract the data around it. + +pointer procedure dp_subrast (im, xcen, ycen, radius, x1, x2, y1, y2) + +pointer im # input image descriptor +real xcen, ycen # center of subraster +real radius # radius of the box +int x1, y1, x2, y2 # boundaries of subraster + +int j, ncols +pointer buf, ptr, imbuf +pointer imgs2r() + +begin + # Calculate start position of extraction box. + x1 = int (xcen - radius) - 2 + x2 = int (xcen + radius) + 3 + y1 = int (ycen - radius) - 2 + y2 = int (ycen + radius) + 3 + if (x1 > IM_LEN(im,1) || x2 < 1 || y1 > IM_LEN(im,2) || y2 < 1) + return (NULL) + + x1 = max (1, x1) + x2 = min (IM_LEN(im,1), x2) + y1 = max (1, y1) + y2 = min (IM_LEN(im,2), y2) + call malloc (buf, (x2 - x1 + 1) * (y2 - y1 + 1), TY_REAL) + + ptr = buf + ncols = x2 - x1 + 1 + do j = y1, y2 { + imbuf = imgs2r (im, x1, x2, j, j) + call amovr (Memr[imbuf], Memr[ptr], ncols) + ptr = ptr + ncols + } + + return (buf) +end + + +# DP_CHKSR -- Check the input subraster for bad pixels. + +int procedure dp_chksr (x, y, sr, xdim, ydim, x1, y1, psfrad, fitrad, lowbad, + highbad, saturated, dmin, dmax, gmax) + +real x, y # position of the star +real sr[xdim,ydim] # the data subraster +int xdim, ydim # the dimensions of the subraster +int x1, y1 # the lower left hand coordinates of the array +real psfrad # the psf radius +real fitrad # the fitting radius +real lowbad, highbad # the good data limits +int saturated # is the star saturated +real dmin, dmax # output data limits +real gmax # maximum good data limit + +int i,j +real pradsq, fradsq, dy2, r2 + +begin + pradsq = psfrad * psfrad + fradsq = fitrad * fitrad + dmin = MAX_REAL + dmax = -MAX_REAL + gmax = -MAX_REAL + saturated = NO + + # Loop through the pixels looking for bad values. + do j = 1, ydim { + dy2 = (y - (y1 + j -1)) ** 2 + if (dy2 > pradsq) + next + do i = 1, xdim { + r2 = (x - (x1 + i - 1)) ** 2 + dy2 + if (r2 > pradsq) + next + if (sr[i,j] < dmin) + dmin = sr[i,j] + if (sr[i,j] > dmax) + dmax = sr[i,j] + if (sr[i,j] < lowbad || sr[i,j] > highbad) { + if (r2 <= fradsq) { + if (sr[i,j] < lowbad) + return (ERR) + else if (saturated == NO) + saturated = YES + } + } else { + if (sr[i,j] > gmax) + gmax = sr[i,j] + } + } + } + + return (OK) +end |