define SZ_KEYW 11 # size of buffer for header keyword include # encounts -- encircled counts # This routine adds up the counts encircled within a specified radial # distance of each star in an image; the image is of the form created # by fextract. The number of counts is normalized by dividing by the # exposure time, as gotten from the input image header. # If a box size is specified (cboxsize >= 3) then mpc_cntr is called to # find an improved value for the position. # The star position and value of the normalized counts for each star are # written to the output text file. These are preceded by "#" to comment # them out. For each refoc position, the sum of the values for all stars # for that refoc position is also written, preceded by the refoc value. # # Phil Hodge, 22-Mar-1991 task created # Phil Hodge, 24-Jun-1991 include option to call mpc_cntr; only one radius # Phil Hodge, 9-Jan-1992 read fextract output procedure encounts() pointer input # scr for name of input image pointer output # scr for name of output ascii file real radius # include counts within this radius (pixels) int cboxsize # size of box for determining position #-- pointer sp pointer rownum # scratch for list of row numbers pointer im # pointer to imhdr struct pointer v # pointer to a section of data from the image char keyword[SZ_KEYW] # keyword for getting refoc position real refoc # refoc mirror position (just copied to output) real minval, maxval # limits of values in image section real counts # counts within radius real sum # sum of counts for each star real x, y # location of a point real xp, yp # x & y within image section int fd # fd for output text file int box1 # separation betw different stars on a row int box2 # separation betw stars at diff focus positions int x1, x2, y1, y2 # limits of image section int j, k # loop indexes pointer immap(), imgs2r() real clgetr(), imgetr() int open(), imgeti(), imaccf() int clgeti() int findrow() begin call smark (sp) call salloc (input, SZ_FNAME, TY_CHAR) call salloc (output, SZ_FNAME, TY_CHAR) call clgstr ("input", Memc[input], SZ_FNAME) call clgstr ("output", Memc[output], SZ_FNAME) radius = clgetr ("radius") cboxsize = clgeti ("cboxsize") if (cboxsize < 3) cboxsize = 0 im = immap (Memc[input], READ_ONLY, NULL) fd = open (Memc[output], NEW_FILE, TEXT_FILE) # Get the list of numbers of rows that have been filled in. call salloc (rownum, SZ_LINE, TY_CHAR) call imgstr (im, "rownum", Memc[rownum], SZ_LINE) # Get box size from image header. box1 = imgeti (im, "box1") box2 = imgeti (im, "box2") # Print header giving radius and image name. call fprintf (fd, "# r =%10.2f\n") call pargr (radius) call fprintf (fd, "# [%s]\n") call pargstr (Memc[input]) # Do for each refoc position. do k = 1, IM_LEN(im,2) / box2 { # Has this row been filled in? If not, skip it. if (findrow (Memc[rownum], k) == NO) next # Get the refoc mirror position from the input header. call sprintf (keyword, SZ_KEYW, "refoc%d") call pargi (k) if (imaccf (im, keyword) == YES) refoc = imgetr (im, keyword) else refoc = 0. # We're going to add up the counts for all stars at a given # refoc position, but we'll print individual values as well. sum = 0. # initial value # Do for each star along a row at a given refoc position. do j = 1, IM_LEN(im,1) / box1 { x1 = (j - 1) * box1 + 1 x2 = j * box1 y1 = (k - 1) * box2 + 1 y2 = k * box2 x = real (x1 + x2) / 2. # nominal position of star y = real (y1 + y2) / 2. # Get image section for current star. v = imgs2r (im, x1, x2, y1, y2) # If there's no data for this star, skip it. call alimr (Memr[v], box1*box2, minval, maxval) if (minval == maxval) next # Get an improved position; x & y are modified in-place. if (cboxsize > 0) call mpc_cntr (im, x, y, cboxsize, x, y) xp = x - x1 + 1. # relative to image section yp = y - y1 + 1. # Get section again (because mpc_cntr did so also). v = imgs2r (im, x1, x2, y1, y2) # Compute encircled counts. call en_calc (Memr[v], box1, box2, xp, yp, radius, counts) call fprintf (fd, "#%8.2f %8.2f") call pargr (x) call pargr (y) # Accumulate sum of counts and print. if (IS_INDEF(counts)) { call fprintf (fd, " * \n") } else { sum = sum + counts call fprintf (fd, "%15.6g\n") call pargr (counts) } } # This is the result for current refoc position. call fprintf (fd, "%8.5g %15.6g\n") call pargr (refoc) call pargr (sum) } call close (fd) call imunmap (im) call sfree (sp) end # en_calc -- calculate encircled counts procedure en_calc (v, box1, box2, x, y, radius, counts) real v[box1,box2] # i: image section int box1, box2 # i: size of image section real x, y # i: pixel coordinates of star within section real radius # i: sum counts within this number of pixels real counts # o: sum of counts within radius #-- real sum # for incrementing counts real fraction # fraction of a pixel included real dx, dy # deltas from x,y to a pixel real r # distance from x,y to a pixel real d # radius - r int i, j begin # Check whether radius extends beyond image section if ((x + radius > box1) || (x - radius < 1) || (y + radius > box2) || (y - radius < 1)) { counts = INDEF } else { sum = 0 do j = 1, box2 { do i = 1, box1 { dx = i - x dy = j - y r = sqrt (dx * dx + dy * dy) d = radius - r if (d > 0.5) { sum = sum + v[i,j] } else if (d > -0.5) { fraction = d + 0.5 sum = sum + v[i,j] * fraction } } } counts = sum } end # MPC_CNTR -- Compute star center using MPC algorithm # as outlined in Stellar Magnitudes from Digital Images. procedure mpc_cntr (im, xstart, ystart, boxsize, xcntr, ycntr) pointer im real xstart, ystart int boxsize real xcntr, ycntr int x1, x2, y1, y2, half_box int ncols, nrows, nx, ny, try real xinit, yinit pointer bufptr, sp, x_vect, y_vect int imgs2r() begin half_box = (boxsize - 1) / 2 xinit = xstart yinit = ystart # Mark region to extract - use box size ncols = IM_LEN (im, 1) nrows = IM_LEN (im, 2) try = 0 repeat { x1 = amax1 (xinit - half_box, 1.0) +0.5 x2 = amin1 (xinit + half_box, real(ncols)) +0.5 y1 = amax1 (yinit - half_box, 1.0) +0.5 y2 = amin1 (yinit + half_box, real(nrows)) +0.5 nx = x2 - x1 + 1 ny = y2 - y1 + 1 # Extract region around center bufptr = imgs2r (im, x1, x2, y1, y2) # Collapse to two 1-D arrays call smark (sp) call salloc (x_vect, nx, TY_REAL) call salloc (y_vect, ny, TY_REAL) call aclrr (Memr[x_vect], nx) call aclrr (Memr[y_vect], ny) # Sum all rows call mpc_rowsum (Memr[bufptr], Memr[x_vect], nx, ny) # Sum all columns call mpc_colsum (Memr[bufptr], Memr[y_vect], nx, ny) # Find centers call mpc_getcenter (Memr[x_vect], nx, xcntr) call mpc_getcenter (Memr[y_vect], ny, ycntr) # Add in offsets xcntr = xcntr + x1 ycntr = ycntr + y1 call sfree (sp) try = try + 1 if (try == 1) { if ((abs(xcntr-xinit) > 1.0) || (abs(ycntr-yinit) > 1.0)) { xinit = xcntr yinit = ycntr } } else break } end # ROWSUM -- Sum all rows in a raster procedure mpc_rowsum (v, row, nx, ny) int nx, ny real v[nx,ny] real row[ARB] int i, j begin do i = 1, ny do j = 1, nx row[j] = row[j] + v[j,i] end # COLSUM -- Sum all columns in a raster. procedure mpc_colsum (v, col, nx, ny) int nx, ny real v[nx,ny] real col[ARB] int i, j begin do i = 1, ny do j = 1, nx col[j] = col[j] + v[i,j] end # GETCENTER -- Compute center of gravity of array. procedure mpc_getcenter (v, nv, vc) real v[ARB] int nv real vc int i real sum1, sum2, sigma, cont begin # Assume continuum level is at endpoints # Compute first moment sum1 = 0.0 sum2 = 0.0 call aavgr (v, nv, cont, sigma) do i = 1, nv if (v[i] > cont) { sum1 = sum1 + (i-1) * (v[i] - cont) sum2 = sum2 + (v[i] - cont) } # Determine center vc = sum1 / sum2 end