include include include include "../lib/xphot.h" include "../lib/impars.h" include "../lib/find.h" include "../lib/objects.h" include "../lib/center.h" include "../lib/fitsky.h" # XP_BPHOT -- Compute centers, sky values, and photometry for a list of # objects. int procedure xp_bphot (im, ol, rl, xp, verbose) pointer im #I the input image descriptor int ol #I the input objects file descriptor int rl #I the output results file descriptor pointer xp #I the pointer to the main xapphot structure int verbose #I print summary of output ? int cier, sier, pier, loseqno, seqno pointer sp, imname, olname, symbol real xc, yc, xver, yver int xp_robjects(), xp_fitcenter(), xp_ofitsky(), xp_magp(), xp_stati() pointer sthead(), xp_statp() real xp_statr() begin if (ol == NULL) return call smark (sp) call salloc (imname, SZ_FNAME, TY_CHAR) call salloc (olname, SZ_FNAME, TY_CHAR) call xp_stats (xp, IMAGE, Memc[imname], SZ_FNAME) call xp_stats (xp, OBJECTS, Memc[olname], SZ_FNAME) loseqno = 0 seqno = SEQNO(xp_statp(xp,PSTATUS)) while (xp_robjects (ol, xp, RLIST_TEMP) > 0) { #psymbolo = sthead (xp_statp (xp, POLYGONLIST)) #psymbols = stnext (xp_statp (xp, POLYGONLIST), psymbolo) #if (psymbolo != NULL && psymbols != NULL) { #if (XP_ONVERTICES(psymbols) > 0) { #symbol = psymbolo #psymbolo = psymbols #psymbols = symbol #} #} else if (psymbolo != NULL) { #if (XP_SNVERTICES(psymbolo) > 0) { #psymbols = psymbolo #psymbolo = NULL #} else #psymbols = NULL #} else if (psymbols != NULL) { #if (XP_ONVERTICES(psymbols) > 0) { #psymbolo = psymbols #psymbols = NULL #} else #psymbolo = NULL #} symbol = sthead (xp_statp (xp, OBJLIST)) if (symbol == NULL) break loseqno = loseqno + 1 xc = XP_OXINIT(symbol) yc = XP_OYINIT(symbol) cier = xp_fitcenter (xp, im, xc, yc) sier = xp_ofitsky (xp, im, symbol, xp_statr (xp, XSHIFT), xp_statr (xp, YSHIFT), xver, yver, 0, NULL, NULL) pier = xp_magp (xp, im, symbol, xp_statr (xp, XSHIFT), xp_statr (xp, YSHIFT), xver, yver, 0, xp_statr(xp, SKY_MODE), xp_statr(xp, SKY_STDEV), xp_stati (xp, NSKY)) if (rl == NULL) { if (verbose == YES) call xp_pqprint (xp, Memc[imname], cier, sier, pier, NO) next } else if (verbose == YES) call xp_pqprint (xp, Memc[imname], cier, sier, pier, YES) seqno = seqno + 1 call xp_pwrite (xp, rl, seqno, Memc[olname], loseqno, cier, sier, pier) } SEQNO(xp_statp(xp,PSTATUS)) = seqno call sfree (sp) return (loseqno) end # XP_APHOT -- Compute magnitudes for a list of objects. int procedure xp_aphot (gd, im, rl, xp, verbose) pointer gd #I the graphics stream descriptor pointer im #I the input image descriptor int rl #I the output results file descriptor pointer xp #I the pointer to the main xapphot structure int verbose #I print summary of output ? int i, j, k, fwidth, swidth, norm, nxb, nyb, nstars, l1, l2, nlines, ntotal int c1, c2, ncols, cier, sier, pier, seqno pointer gker2d, ngker2d, skip, imbuf, denbuf, x, y real sigma, nsigma, a, b, c, f, dmin, dmax, relerr, xver, yver real gsums[LEN_GAUSS] int xp_staget(), xp_fitcenter(), xp_fitsky(), xp_mag(), xp_stati() pointer xp_statp() real xp_statr(), xp_egkernel() begin # Compute the parameters of the Gaussian kernel. sigma = HWHM_TO_SIGMA * xp_statr (xp, IHWHMPSF) * xp_statr (xp, ISCALE) nsigma = xp_statr (xp, FRADIUS) * xp_statr (xp, ISCALE) / HWHM_TO_SIGMA call xp_egparams (sigma, 1.0, 0.0, nsigma, a, b, c, f, fwidth, fwidth) # Compute the separation parameter swidth = max (2, int (xp_statr (xp, FSEPMIN) * xp_statr (xp, IHWHMPSF) * xp_statr (xp, ISCALE) + 0.5)) # Compute the minimum and maximum pixel values. if (IS_INDEFR(xp_statr(xp,IMINDATA)) && IS_INDEFR(xp_statr(xp,IMAXDATA))) { norm = YES dmin = -MAX_REAL dmax = MAX_REAL } else { norm = NO if (IS_INDEFR(xp_statr(xp,IMINDATA))) dmin = -MAX_REAL else dmin = xp_statr (xp,IMINDATA) if (IS_INDEFR(xp_statr(xp,IMAXDATA))) dmax = MAX_REAL else dmax = xp_statr(xp,IMAXDATA) } # Set up the image boundary extension characteristics. call imseti (im, IM_TYBNDRY, BT_NEAREST) call imseti (im, IM_NBNDRYPIX, 1 + fwidth / 2 + swidth) # Set up the blocking factor. nxb = IM_LEN(im,1) nyb = DEF_NYBLOCK # Process the image block by block. ntotal = 0 seqno = SEQNO(xp_statp(xp,PSTATUS)) do j = 1, IM_LEN(im,2), nyb { l1 = j l2 = min (IM_LEN(im,2), j + nyb - 1) nlines = l2 - l1 + 1 + 2 * (fwidth / 2 + swidth) do i = 1, IM_LEN(im,1), nxb { # Allocate space for the convolution kernel. call malloc (gker2d, fwidth * fwidth, TY_REAL) call malloc (ngker2d, fwidth * fwidth, TY_REAL) call malloc (skip, fwidth * fwidth, TY_INT) # Compute the convolution kernels. relerr = xp_egkernel (Memr[gker2d], Memr[ngker2d], Memi[skip], fwidth, fwidth, gsums, a, b, c, f) # Allocate space for the data and the convolution. c1 = i c2 = min (IM_LEN(im,1), i + nxb - 1) ncols = c2 - c1 + 1 + 2 * (fwidth / 2 + swidth) call malloc (imbuf, ncols * nlines, TY_REAL) call malloc (denbuf, ncols * nlines, TY_REAL) call malloc (x, ncols * nlines, TY_REAL) call malloc (y, ncols * nlines, TY_REAL) # Do the convolution. if (norm == YES) call xp_fconvolve (im, c1, c2, l1, l2, swidth, Memr[imbuf], Memr[denbuf], ncols, nlines, Memr[ngker2d], Memi[skip], fwidth, fwidth) else call xp_gconvolve (im, c1, c2, l1, l2, swidth, Memr[imbuf], Memr[denbuf], ncols, nlines, Memr[gker2d], Memi[skip], fwidth, fwidth, gsums, dmin, dmax) # Find, the stars. nstars = xp_staget (xp, Memr[imbuf], Memr[denbuf], Memr[x], Memr[y], ncols, nlines, c1, c2, l1, l2, swidth, Memi[skip], fwidth, fwidth, dmin, dmax) # Center and record the stars do k = 1, nstars { cier = xp_fitcenter (xp, im, Memr[x+k-1], Memr[y+k-1]) Memr[x+k-1] = xp_statr (xp, XCENTER) Memr[y+k-1] = xp_statr (xp, YCENTER) sier = xp_fitsky (xp, im, Memr[x+k-1], Memr[y+k-1], xver, yver, 0, NULL, gd) pier = xp_mag (xp, im, Memr[x+k-1], Memr[y+k-1], xver, yver, 0, xp_statr(xp,SKY_MODE), xp_statr (xp, SKY_STDEV), xp_stati(xp, NSKY)) if (gd != NULL) call xp_amark (gd, xp, xver, yver, 0, xver, yver, 0, 1, 1) if (rl == NULL) { if (verbose == YES) call xp_pqprint (xp, IM_HDRFILE(im), cier, sier, pier, NO) next } else if (verbose == YES) call xp_pqprint (xp, IM_HDRFILE(im), cier, sier, pier, YES) call xp_pwrite (xp, rl, seqno + k, "none", 0, cier, sier, pier) } # Mark the stars using the computed centers and the # object marking routine. if (gd != NULL) call xp_mkfobjects (gd, xp, Memr[x], Memr[y], nstars, ntotal + seqno, 1, 1) # Increment the star counter and sequence number. ntotal = ntotal + nstars if (rl != NULL) seqno = seqno + nstars # Free the memory. call mfree (gker2d, TY_REAL) call mfree (ngker2d, TY_REAL) call mfree (skip, TY_INT) call mfree (imbuf, TY_REAL) call mfree (denbuf, TY_REAL) call mfree (x, TY_REAL) call mfree (y, TY_REAL) } } SEQNO(xp_statp(xp,PSTATUS)) = seqno return (ntotal) end # XP_UAPHOT -- Compute centers for a list of objects. int procedure xp_uaphot (gd, ui, im, rl, xp, verbose) pointer gd #I the graphics stream descriptor pointer ui #I the pointer to the user interface pointer im #I the input image descriptor int rl #I the output results file descriptor pointer xp #I the pointer to the main xapphot structure int verbose #I print summary of output ? int i, j, k, fwidth, swidth, norm, nxb, nyb, nstars, l1, l2, nlines, ntotal int c1, c2, ncols, cier, sier, pier, seqno pointer gker2d, ngker2d, skip, imbuf, denbuf, x, y real sigma, nsigma, a, b, c, f, dmin, dmax, relerr, xver, yver real gsums[LEN_GAUSS] int xp_staget(), xp_fitcenter(), xp_fitsky(), xp_mag(), xp_stati() pointer xp_statp() real xp_statr(), xp_egkernel() begin # Compute the parameters of the Gaussian kernel. sigma = HWHM_TO_SIGMA * xp_statr (xp, IHWHMPSF) * xp_statr (xp, ISCALE) nsigma = xp_statr (xp, FRADIUS) * xp_statr (xp, ISCALE) / HWHM_TO_SIGMA call xp_egparams (sigma, 1.0, 0.0, nsigma, a, b, c, f, fwidth, fwidth) # Compute the separation parameter swidth = max (2, int (xp_statr (xp, FSEPMIN) * xp_statr (xp, IHWHMPSF) * xp_statr (xp, ISCALE) + 0.5)) # Compute the minimum and maximum pixel values. if (IS_INDEFR(xp_statr(xp,IMINDATA)) && IS_INDEFR(xp_statr(xp,IMAXDATA))) { norm = YES dmin = -MAX_REAL dmax = MAX_REAL } else { norm = NO if (IS_INDEFR(xp_statr(xp,IMINDATA))) dmin = -MAX_REAL else dmin = xp_statr (xp,IMINDATA) if (IS_INDEFR(xp_statr(xp,IMAXDATA))) dmax = MAX_REAL else dmax = xp_statr(xp,IMAXDATA) } # Set up the image boundary extension characteristics. call imseti (im, IM_TYBNDRY, BT_NEAREST) call imseti (im, IM_NBNDRYPIX, 1 + fwidth / 2 + swidth) # Set up the blocking factor. nxb = IM_LEN(im,1) nyb = DEF_NYBLOCK # Process the image block by block. ntotal = 0 seqno = SEQNO(xp_statp(xp,PSTATUS)) do j = 1, IM_LEN(im,2), nyb { l1 = j l2 = min (IM_LEN(im,2), j + nyb - 1) nlines = l2 - l1 + 1 + 2 * (fwidth / 2 + swidth) do i = 1, IM_LEN(im,1), nxb { # Allocate space for the convolution kernel. call malloc (gker2d, fwidth * fwidth, TY_REAL) call malloc (ngker2d, fwidth * fwidth, TY_REAL) call malloc (skip, fwidth * fwidth, TY_INT) # Compute the convolution kernels. relerr = xp_egkernel (Memr[gker2d], Memr[ngker2d], Memi[skip], fwidth, fwidth, gsums, a, b, c, f) # Allocate space for the data and the convolution. c1 = i c2 = min (IM_LEN(im,1), i + nxb - 1) ncols = c2 - c1 + 1 + 2 * (fwidth / 2 + swidth) call malloc (imbuf, ncols * nlines, TY_REAL) call malloc (denbuf, ncols * nlines, TY_REAL) call malloc (x, ncols * nlines, TY_REAL) call malloc (y, ncols * nlines, TY_REAL) # Do the convolution. if (norm == YES) call xp_fconvolve (im, c1, c2, l1, l2, swidth, Memr[imbuf], Memr[denbuf], ncols, nlines, Memr[ngker2d], Memi[skip], fwidth, fwidth) else call xp_gconvolve (im, c1, c2, l1, l2, swidth, Memr[imbuf], Memr[denbuf], ncols, nlines, Memr[gker2d], Memi[skip], fwidth, fwidth, gsums, dmin, dmax) # Find, the stars. nstars = xp_staget (xp, Memr[imbuf], Memr[denbuf], Memr[x], Memr[y], ncols, nlines, c1, c2, l1, l2, swidth, Memi[skip], fwidth, fwidth, dmin, dmax) # Center and record the stars do k = 1, nstars { cier = xp_fitcenter (xp, im, Memr[x+k-1], Memr[y+k-1]) Memr[x+k-1] = xp_statr (xp, XCENTER) Memr[y+k-1] = xp_statr (xp, YCENTER) sier = xp_fitsky (xp, im, Memr[x+k-1], Memr[y+k-1], xver, yver, 0, NULL, gd) pier = xp_mag (xp, im, Memr[x+k-1], Memr[y+k-1], xver, yver, 0, xp_statr(xp,SKY_MODE), xp_statr (xp, SKY_STDEV), xp_stati(xp, NSKY)) call xp_tmkresults (gd, ui, xp) if (gd != NULL) call xp_amark (gd, xp, xver, yver, 0, xver, yver, 0, 1, 1) if (rl == NULL) next call xp_pwrite (xp, rl, seqno + k, "none", 0, cier, sier, pier) } # Mark the stars using the computed centers and the # object marking routine. if (gd != NULL) call xp_mkfobjects (gd, xp, Memr[x], Memr[y], nstars, ntotal + seqno, 1, 1) # Increment the star counter and sequence number. ntotal = ntotal + nstars if (rl != NULL) seqno = seqno + nstars # Free the memory. call mfree (gker2d, TY_REAL) call mfree (ngker2d, TY_REAL) call mfree (skip, TY_INT) call mfree (imbuf, TY_REAL) call mfree (denbuf, TY_REAL) call mfree (x, TY_REAL) call mfree (y, TY_REAL) } } SEQNO(xp_statp(xp,PSTATUS)) = seqno return (ntotal) end