include define SZ_KEYW 11 # size of buffer for header keyword # faradius -- compute a measure of the radius # This task uses autocorrelation to obtain a value indicating the spread # in the psf. # The output from fextract is the input to this routine. # # Phil Hodge, 2-Feb-1990 Task created. # Phil Hodge, 20-Apr-1990 Do autocorrelation separately for each star. procedure faradius() pointer input # scr for name of input image pointer output # scr for name of output ascii file pointer outim # scr for name of output image bool verbose # print refoc position? #-- pointer sp pointer rownum # scratch for list of row numbers pointer ac # scr for autocorrelation data pointer im # pointer to imhdr struct pointer oim # imhdr for output image pointer v # pointer to input data char keyword[SZ_KEYW] # keyword for getting refoc position real refoc # position of refoc mirror real radius # characteristic radius of the autocorrelation int fd # fd for output file int numimages # number of focus positions int naxis1 # width of input image int box1 # separation betw different stars in a slice int box2 # separation betw slices at diff focus positions int nout # size of autocorrelation (square) int y1, y2 # limits of image section on second axis int k # loop index for focus position pointer immap(), imgs2r() int open(), imgeti(), imaccf() int findrow() real imgetr() bool clgetb() begin call smark (sp) call salloc (input, SZ_FNAME, TY_CHAR) call salloc (output, SZ_FNAME, TY_CHAR) call salloc (outim, SZ_FNAME, TY_CHAR) call clgstr ("input", Memc[input], SZ_FNAME) call clgstr ("output", Memc[output], SZ_FNAME) call clgstr ("outim", Memc[outim], SZ_FNAME) verbose = clgetb ("verbose") 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) # fextract puts different stars along the horizontal direction # and different focus positions (images input to fextract) in # the vertical direction. The slice we are going to take from # the input image is of size naxis1 x box2. naxis1 = IM_LEN(im,1) box1 = imgeti (im, "box1") box2 = imgeti (im, "box2") numimages = IM_LEN(im,2) / box2 # should divide exactly # Allocate space for the result of the autocorrelation, which may be # smaller than the input slice. nout = min (box1, box2) call salloc (ac, nout*nout, TY_REAL) # Create an output image for the autocorrelations. The width will # be nout, for one AC, but the height will be numimages times larger, # one for the AC of each slice. oim = immap (Memc[outim], NEW_IMAGE, NULL) IM_NDIM(oim) = 2 IM_LEN(oim,1) = nout IM_LEN(oim,2) = nout * numimages # Process each focus position. do k = 1, numimages { # Has this row been filled in? If not, just fill the current # box of the autocorrelation with zero. if (findrow (Memc[rownum], k) == NO) { call amovkr (0., Memr[ac], nout*nout) call far_ac_save (oim, Memr[ac], nout, k) next } # Get the refoc mirror position from the input header, # and add it to the output header. call sprintf (keyword, SZ_KEYW, "refoc%d") call pargi (k) if (imaccf (im, keyword) == YES) refoc = imgetr (im, keyword) else refoc = 0. call imaddr (oim, keyword, refoc) if (verbose) { call printf ("refoc %g\n") call pargr (refoc) call flush (STDOUT) } # Get the data (a horizontal slice of size naxis1 x box2) # for all stars for a given focus position. y1 = (k - 1) * box2 + 1 y2 = k * box2 v = imgs2r (im, 1, naxis1, y1, y2) # Compute the autocorrelation of this image section. call far_auto (Memr[v], naxis1, box1, box2, Memr[ac], nout) # Compute the radius of the autocorrelation. call far_getrad (Memr[ac], nout, radius) # Write the refoc position and the radius to the file. call fprintf (fd, "%15.7g %15.7g\n") call pargr (refoc) call pargr (radius) call flush (fd) # Save the autocorrelation in the output image. call far_ac_save (oim, Memr[ac], nout, k) } call close (fd) call imunmap (oim) call imunmap (im) call sfree (sp) end # far_auto -- compute autocorrelation # This routine computes the autocorrelation of a 2-D array. Zero lag is # output in pixel [1,1] of ac. The autocorrelation is computed only for # a portion of the input, from lag = 0 to lag = nout-1, and similarly for # the second axis. The input array consists of nstars boxes of size # box1 x box2. The AC is computed separately for each such box and added # together. procedure far_auto (v, naxis1, box1, box2, ac, nout) real v[naxis1,box2] # i: input section int naxis1 # i: size of first axis of input image int box1, box2 # i: size of one box real ac[nout,nout] # o: autocorrelation of v int nout # i: size of ac #-- real sum # for accumulating sum of products int nstars # number of stars (boxes) at each focus position int star # loop index for star (box) number int off # offset in v for current star int j, k # loop indexes int m, n # offsets for autocorrelation begin nstars = naxis1 / box1 # should divide exactly do n = 1, nout do m = 1, nout ac[m,n] = 0. # Do for each box (star) in the input slice. do star = 1, nstars { # Offset along first axis of v to beginning of current box. off = (star - 1) * box1 # Loop over different values of lag for each axis. do n = 0, nout-1 { do m = 0, nout-1 { # Accumulate the sum by looping over the portion of the # current box that overlaps with itself for the current lag. sum = 0. do k = 1, box2-n do j = 1, box1-m sum = sum + v[j+off,k] * v[j+off+m,k+n] ac[m+1,n+1] = ac[m+1,n+1] + sum } } } # Divide the autocorrelation by the number of pixels. do n = 1, nout do m = 1, nout ac[m,n] = ac[m,n] / (nout*nout) end # far_getrad -- get the radius of the autocorrelation # This computes the average radius weighted by the value of the auto- # correlation. Points farther than nout from pixel [1,1] are ignored. procedure far_getrad (ac, nout, radius) real ac[nout,nout] # i: autocorrelation int nout # i: size of ac real radius # o: radius of autocorrelation #-- double sumv, sumrv # for accumulating sums double r # distance from [1,1] to a pixel double x, y # distance from [1,1] along x & y axes int j, k # loop indexes begin sumv = 0. sumrv = 0. do k = 1, nout { y = k - 1. do j = 1, nout { x = j - 1. r = sqrt (x*x + y*y) if (r <= nout) { sumv = sumv + ac[j,k] sumrv = sumrv + r * ac[j,k] } } } if (sumv > 0.) radius = sumrv / sumv else radius = 0. end # far_ac_save -- save the autocorrelation procedure far_ac_save (oim, ac, nout, k) pointer oim # i: imhdr pointer for output image real ac[nout,nout] # i: autocorrelation int nout # i: size of ac int k #-- pointer x int j # index in ac array int line # loop index for line of image int first, last # range of output lines pointer impl2r() begin first = (k - 1) * nout + 1 last = k * nout j = 1 do line = first, last { x = impl2r (oim, line) call amovr (ac[1,j], Memr[x], nout) j = j + 1 } end