include # for IM_LEN include "resfind.h" # gridextrap -- Reseau finding routine # # Dave Baxter, May-1993 Original Fortran version. # Phil Hodge, 14-Jun-1993 Convert to SPP; call rfind subroutines. # Phil Hodge, 14-Feb-1994 Remove pc_coords from calling sequence. # # # Method: # This routine will attempt to find the approximate (i.e. # nearest pixel) reseau positions by the relatively simple method # of `Grid Extrapolation'. Basically, given an estimate of the # offsets to the next reseaux in the row (Dxr & Dyr) and to the # next reseaux in the column (Dxc & Dyc), the procedure will com- # plete the small 3x3 set of reseaux surrounding a reseau with a # previously determined position (Px,Py). # # X # # X # # X # # X # # # # X # # X # Dxc (8) # X--- # (7) | # X | # (6) | # | # | Dyc # | X # | # | X # Dxr | (5) # | ------------- X # Dyr| (Px,Py) # X # (4) # # # X # # X # (3) # X # (2) # X # (1) # # In order for the procedure to start, it requires intial # estimates for the central reseau (Xo,Yo) and the column and row # offsets (Dxc,Dyc,Dxr,Dyr). The routine first refines the position # for the central reseau, and then completes the 3x3 region around it # (in the order shown by the number scheme above). The initial values # are then discarded and new values calculated for each subsequent row # based on the offsets between reseaux which are already found. # # As they are found, the reseau positions are written into a # 3-Dimensional array, IRp(17,17,3), as follows: # # IRp(I,J,1) = IRp(Column,Row,Xpos) # IRp(I,J,2) = IRp(Column,Row,Ypos) # IRp(I,J,3) = YES if position was found, NO if not procedure gridextrap (im, rm, Xo, Yo, Dxr, Dyr, Dxc, Dyc, xoff, yoff, searchx, searchy, depth, verbose, noindefs, IRp, nfound) pointer im # i: imhdr struct for input image pointer rm # i: reseau model data real Xo, Yo # i: estimated position of central reseau real Dxr # i: X-offset along row \ Depends on position real Dyr # i: Y-offset along row | of (Irx,Iry) relative real Dxc # i: X-offset along column | to reseau (Px,Py) real Dyc # i: Y-offset along column / real xoff # o: X offset of image from beginning of photocathode real yoff # o: Y offset of image from beginning of photocathode int searchx, searchy # i: size of search region real depth # i: correlation depth required bool verbose # i: print info about positions as they are found? bool noindefs # i: save approximate position if reseau not found? real IRp[NRX,NRY,3] # o: reseau positions int nfound # o: number of reseaux actually found #-- pointer sp pointer corr # work area int mx, my # size of largest reseau model int csx, csy # size of search region for finding central reseau int i int xcent, ycent # indexes (in IRp) of central reseau int n_l, n_r, n_t, n_b # for finding maxPres int maxPres # upper limit to loop index Pres int Pres # Row/Column offset to current perimeter. int Px # Column number of current 3x3 grid centre int Py # Row number of current 3x3 grid centre. int Irx # Column number of current reseau int Iry # Row number of current reseau real Ix # Estimated X-coord of (Irx,Iry) real Iy # Estimated Y-coord of (Irx,Iry) bool notfound # true if reseau was not found begin call smark (sp) csx = searchx * 2 + 1 csy = searchy * 2 + 1 call salloc (corr, csx * csy, TY_REAL) # These are the indexes in IRp of the central reseau mark. xcent = NRX / 2 + 1 ycent = NRY / 2 + 1 # Initialize array of flags to say we haven't looked for any yet. do iry = 1, NRY do irx = 1, NRX IRp(irx,iry,3) = INDEFR # Get the size of the largest model array. mx = 0 my = 0 do i = 1, MODEL_NUM(rm) { mx = max (MODEL_NAXIS1(rm,i), mx) my = max (MODEL_NAXIS2(rm,i), my) } # Take care of CENTRAL reseau. Irx = xcent Iry = ycent Ix = Xo Iy = Yo # Find central reseau mark; update IRp. # Note that we use a larger search box for this one. call xresfind1 (im, Irx, Iry, Ix, Iy, rm, csx, csy, depth, mx, my, verbose, noindefs, Memr[corr], IRp, notfound) if (notfound) call error (1, "can't find central reseau") # Find out how far away from the central reseau we must go. n_l = (IRp[Irx,Iry,1] - 1) / RES_SPACING_X n_r = (IM_LEN(im,1) - IRp[Irx,Iry,1]) / RES_SPACING_X n_b = (IRp[Irx,Iry,2] - 1) / RES_SPACING_Y n_t = (IM_LEN(im,2) - IRp[Irx,Iry,2]) / RES_SPACING_Y maxPres = max (n_l, n_r, n_b, n_t) maxPres = min (maxPres, NRESX, NRESY) maxPres = maxPres - 1 # - 1 because of 3 x 3 loop # Run rings around the central reseau. nfound = 1 # number found so far Do Pres = 0, maxPres { # Loop around the known perimeter. Do Py = ycent-Pres, ycent+Pres { Do Px = xcent-Pres, xcent+Pres { # Loop around a 3x3 subset. Do Iry = Py-1, Py+1 { Do Irx = Px-1, Px+1 { if (Irx < 1 || Irx > NRX || Iry < 1 || Iry > NRY) next # out of limits # Goto next if position already known or if we have # looked unsuccessfully for it. If (!IS_INDEF (IRp(Irx,Iry,3))) next If (Irx < Px) { If (Iry < Py) { If (Pres > 1) { Dxr = Abs (IRp(Px,Py,1) - IRp(Px+1,Py,1)) Dyr = Abs (IRp(Px,Py,2) - IRp(Px+1,Py,2)) Dxc = Abs (IRp(Px,Py,1) - IRp(Px,Py+1,1)) Dyc = Abs (IRp(Px,Py,2) - IRp(Px,Py+1,2)) } Ix = IRp(Px,Py,1) - Dxr + Dxc Iy = IRp(Px,Py,2) - Dyr - Dyc } Else if (Iry > Py) { If (Pres > 1) { Dxr = Abs (IRp(Px,Py,1) - IRp(Px+1,Py,1)) Dyr = Abs (IRp(Px,Py,2) - IRp(Px+1,Py,2)) Dxc = Abs (IRp(Px,Py,1) - IRp(Px,Py-1,1)) Dyc = Abs (IRp(Px,Py,2) - IRp(Px,Py-1,2)) } Ix = IRp(Px,Py,1) - Dxr - Dxc Iy = IRp(Px,Py,2) - Dyr + Dyc } Else { If (Pres > 1) { Dxr = Abs (IRp(Px,Py,1) - IRp(Px+1,Py,1)) Dyr = Abs (IRp(Px,Py,2) - IRp(Px+1,Py,2)) } Ix = IRp(Px,Py,1) - Dxr Iy = IRp(Px,Py,2) - Dyr } } If (Irx > Px) { If (Iry < Py) { If (Pres > 1) { Dxr = Abs (IRp(Px,Py-1,1) - IRp(Px-1,Py-1,1)) Dyr = Abs (IRp(Px,Py-1,2) - IRp(Px-1,Py-1,2)) Dxc = Abs (IRp(Px,Py,1) - IRp(Px,Py+1,1)) Dyc = Abs (IRp(Px,Py,2) - IRp(Px,Py+1,2)) } Ix = IRp(Px,Py,1) + Dxr + Dxc Iy = IRp(Px,Py,2) + Dyr - Dyc } Else if (Iry > Py) { If (Pres > 1) { Dxr = Abs (IRp(Px-1,Py,1) - IRp(Px,Py,1)) Dyr = Abs (IRp(Px-1,Py,2) - IRp(Px,Py,2)) Dxc = Abs (IRp(Px,Py,1) - IRp(Px,Py-1,1)) Dyc = Abs (IRp(Px,Py,2) - IRp(Px,Py-1,2)) } Ix = IRp(Px,Py,1) + Dxr - Dxc Iy = IRp(Px,Py,2) + Dyr + Dyc } Else { If (Pres > 1) { Dxr = Abs (IRp(Px-1,Py,1) - IRp(Px,Py,1)) Dyr = Abs (IRp(Px-1,Py,2) - IRp(Px,Py,2)) } Ix = IRp(Px,Py,1) + Dxr Iy = IRp(Px,Py,2) + Dyr } } If (Irx == Px) { If (Iry < Py) { If (Pres > 1) { Dxc = Abs(IRp(Px,Py,1) - IRp(Px,Py+1,1)) Dyc = Abs(IRp(Px,Py,2) - IRp(Px,Py+1,2)) } Ix = IRp(Px,Py,1) + Dxc Iy = IRp(Px,Py,2) - Dyc } Else { If (Pres > 1) { Dxc = Abs(IRp(Px,Py,1) - IRp(Px,Py-1,1)) Dyc = Abs(IRp(Px,Py,2) - IRp(Px,Py-1,2)) } Ix = IRp(Px,Py,1) - Dxc Iy = IRp(Px,Py,2) + Dyc } } # Find one reseau mark; update IRp. call xresfind1 (im, Irx, Iry, Ix, Iy, rm, searchx, searchy, depth, mx, my, verbose, noindefs, Memr[corr], IRp, notfound) if (!notfound) nfound = nfound + 1 } } } } } # Add sampoff and lineoff (actually xoff & yoff) to positions # to convert to (flipped) photocathode coordinates. do Iry = 1, NRY { do Irx = 1, NRX { If (IS_INDEF (IRp(Irx,Iry,3))) { IRp(Irx,Iry,3) = NO # not found } else { IRp(Irx,Iry,1) = IRp(Irx,Iry,1) + xoff IRp(Irx,Iry,2) = IRp(Irx,Iry,2) + yoff } } } call sfree (sp) end