include include include include include include include "../dbfit.h" define SZ_KEYVAL 11 # size of string for a keyword value define MAXITER 40 # Maximum iterations define OKERROR 0.005 # Acceptable error in Y position # T_REVALFITX -- Evaluates a set of transformation on a given set of # grid points # # This STSDAS application transforms a set of grid positions held # in a reseau file. The transform may be one or two dimensional, the # former usually coming from longslit distortion mapping using the # noao.longslit package and the latter from mapping image tube # distortion using reseau makes. The transforms may be used 'backwards' # i.e. the program will interatively search for points that would be # transformed to the given points by the transformation. This operation # is only likely to work for reasonably well-behaved functions as it is # rather primitive. # # The dimensionality of the transform is deduced from the presence or # otherwise of the "axis" entry in the transform header. # # In the 1d case the direction (i.e. whether X or Y) is determined by # the value of the 'axis' parameter in the transform header. # # This program was originally written by Dave Giaretta of STScI for the # 2d 'forward only' case. It has been generalised to handle the extra # options (largely to handle longslit images) by Richard Hook at ST-ECF. # # Richard Hook, ST-ECF, Jan 88 # Phil Hodge, 3-Mar-1991 Remove print statements. # Phil Hodge, 5-Apr-1991 Remove db_readr, etc from this file. # Remove TYPE from argument list of revalfitx. procedure t_revalfitx() char input1[SZ_FNAME] # input transformation char input2[SZ_FNAME] # ref reseau file with grid char output[SZ_FNAME] # output reseau file char entry[SZ_R_PATTERN] # transform entries to use char refentry[SZ_R_PATTERN] # reference reseau entry char dir[SZ_KEYVAL] # direction as string #-- pointer sp pointer fit, dt pointer dtmap() pointer rs_open() pointer rpout, respt pointer rpref, refpt, xrefpt, yrefpt int dirn # direction of transform int axis int npts int rec int nowhite(), clgwrd() int rs_readr(), db_readr() bool streq() begin call smark (sp) # space for fit pointers call salloc( fit, LEN_FIT, TY_LONG) # get infile1 and check for blanks call clgstr( "input1", input1, SZ_FNAME) if ( nowhite( input1, input1, SZ_FNAME) == 0 ) call error( 0, "input file must be given") dt = dtmap ( input1, READ_ONLY) # input2 call clgstr( "input2", input2, SZ_FNAME) if ( nowhite( input2, input2, SZ_FNAME) == 0) call error( 0, "input2 reseau file must be given") # get entry defs call clgstr( "entry", entry, SZ_R_PATTERN) call clgstr( "refentry", refentry, SZ_R_PATTERN) # open output reseau file call clgstr ("output", output, SZ_FNAME) if ( nowhite( output, output, SZ_FNAME) == 0 ) call error( 0, "output reseau file must be given") if ( streq( input2, output) ) { rpref = rs_open( input2, READ_WRITE, 0) rpout = rpref } else { rpref = rs_open( input2, READ_ONLY, 0) rpout = rs_open( output, NEW_COPY, rpref) } # allocate space for reference & reseau call rs_alloc( rpref, TY_REAL, refpt) call rs_alloc( rpout, TY_REAL, respt) # read reference if (rs_readr (rpref, refpt, refentry, 0) == RES_F_NONEXTENTRY) call error (0, "cannot read reference entry") npts = RES_NROWS(rpref) * RES_NCOLS(rpref) # save the xref, yref pointers, and swap with o/p # in order to transfer the header info xrefpt = RES_XPT(refpt) yrefpt = RES_YPT(refpt) RES_XPT(refpt) = RES_XPT(respt) RES_YPT(refpt) = RES_YPT(respt) # get the direction of the transform dirn = clgwrd ("direction", dir, SZ_KEYVAL, "|forwards|backwards|") if (db_readr (dt, dirn, entry, rec, axis, fit) == RES_F_NONEXTENTRY) call error (1, "entry not found") # set entry name in structure call strcpy (FIT_NAME(fit), RES_ENT_ENTRY(refpt), SZ_FNAME) call revalfitx (fit, axis, npts, Memr[xrefpt], Memr[yrefpt], Memb[RES_FPT(refpt)], Memr[RES_XPT(refpt)], Memr[RES_YPT(refpt)]) call rs_write( rpout, refpt, 0) # close up call dtunmap (dt ) call rs_close ( rpref ) call rs_close ( rpout ) call sfree (sp) end # REVALFITX -- evaluate the coordinate transformations # # Modified, Richard Hook, ST-ECF, August 1988 to handle # longslit mode. # # Further modifications, Richard Hook, ST-ECF, Jan 89 # additional 'axis' argument to allow for both X and Y distortions. procedure revalfitx (fit, axis, npts, x, y, flag, xout, yout) pointer fit # i: pointer to fit parameters int axis # i: Axis number of 1d distortion int npts # i: number of points real x[ARB] # i: x coords to eval fit real y[ARB] # i: y coords to eval fit bool flag[ARB] # i: flag values real xout[ARB] # o: evaluated x positions real yout[ARB] # o: evaluated y positions #-- pointer sx1, sy1, sx2, sy2 real gseval() int i real lamtoy(),lamtox() errchk gseval begin # First consider the original 'simple' case switch (FIT_TYPE(fit)) { case FOR2D: sx1 = FIT_SX1(fit) sy1 = FIT_SY1(fit) sx2 = FIT_SX2(fit) sy2 = FIT_SY2(fit) do i = 1, npts { if (flag[i]) { xout[i] = INDEFR yout[i] = INDEFR } else { xout[i] = gseval( sx1, x[i], y[i] ) yout[i] = gseval( sy1, x[i], y[i] ) if (sx2 != NULL) xout[i] = xout[i] + gseval (sx2, x[i], y[i] ) if (sy2 != NULL) yout[i] = yout[i] + gseval (sy2, x[i], y[i] ) } } case BACK2D: sx1 = FIT_SX1(fit) sy1 = FIT_SY1(fit) sx2 = FIT_SX2(fit) sy2 = FIT_SY2(fit) do i = 1, npts { if (flag[i]) { xout[i] = INDEFR yout[i] = INDEFR } else { call findxy (sx1, sx2, sy1, sy2, x[i], y[i], xout[i], yout[i]) } } case FOR1D: # FOR1D case - the transformation is wavelength as a function # of both x and y. The output values are x,lambda - x is unchanged # from the input sx1 = FIT_SX1(fit) do i = 1, npts { if (flag[i]) { xout[i] = INDEFR yout[i] = INDEFR } else if (axis == 2) { xout[i] = x[i] # X is unchanged yout[i] = gseval( sx1, x[i], y[i] ) } else { yout[i] = y[i] # Y is unchanged xout[i] = gseval( sx1, x[i], y[i] ) } } case BACK1D: # BACK1D case - the transformation is from Lambda to Y at known # X, this is done iteratively sx1 = FIT_SX1(fit) do i = 1, npts { if (flag[i]) { xout[i] = INDEFR yout[i] = INDEFR } else if (axis == 2) { xout[i] = x[i] # X is unchanged yout[i] = lamtoy( sx1, x[i], y[i] ) } else { yout[i] = y[i] # Y is unchanged xout[i] = lamtox( sx1, x[i], y[i] ) } } } end # LAMTOY - find Y such that f(x,y) = Lambda where X and Lambda are known. # # Richard Hook, ST-ECF, 9/8/88 real procedure lamtoy (sf, x, lambda) pointer sf # GSURFIT pointer real x # Known X value real lambda # Known Lambda value #-- real gseval() # Function evaluator real gsgetr() # Get a fit parameter real c # Temp value, lin. coeffs. int i # Loop counter real guess # Guess value real newguess # Improved guess real ymin,ymax # Range of fits real lammin,lammax begin # Get the fits range ymin=gsgetr(sf, GSYMIN) ymax=gsgetr(sf, GSYMAX) # Get a first guess and the slope (c) # by assuming approx linear shape lammin = gseval(sf, x, ymin) lammax = gseval(sf, x, ymax) c=(lammax-lammin)/(ymax-ymin) guess = ((ymax-ymin)*(lambda-lammin) / (lammax-lammin)) + ymin if (c==0.0 || lammin==lammax) return (INDEFR) # Iterate until too many tries or error acceptable do i = 1, MAXITER { newguess = guess + (lambda - gseval(sf, x, guess))/c if (abs(newguess-guess) < OKERROR) return (newguess) guess = newguess } return (newguess) end # LAMTOX - find X such that f(x,y) = Lambda where Y and Lambda are # known # # Richard Hook, ST-ECF, 3/1/89 (from LAMTOY) real procedure lamtox (sf, lambda, y) pointer sf # GSURFIT pointer real y # Known Y value real lambda # Known Lambda value #-- real gseval() # Function evaluator real gsgetr() # Get a fit parameter real c # Temp value, lin. coeffs. int i # Loop counter real guess # Guess value real newguess # Improved guess real xmin,xmax # Range of fits real lammin,lammax begin # Get the fits range xmin=gsgetr(sf, GSXMIN) xmax=gsgetr(sf, GSXMAX) # Get a first guess and the slope (c) # by assuming approx linear shape lammin = gseval(sf, xmin, y) lammax = gseval(sf, xmax, y) c=(lammax-lammin)/(xmax-xmin) guess = ((xmax-xmin)*(lambda-lammin) / (lammax-lammin)) + xmin if (c==0.0 || lammin==lammax) return (INDEFR) # Iterate until too many tries or error acceptable do i = 1, MAXITER { newguess = guess + (lambda - gseval(sf, guess, y)) / c if (abs(newguess-guess) < OKERROR) return (newguess) guess = newguess } return (newguess) end # FINDXY - find X,Y such that f(X,Y) = x,y where x and y and f are known. # # This routine makes rather a lot of simple assumptions and will # not work for unusual cases. # # Richard Hook, ST-ECF, 3/1/89 (from LAMTOY etc) procedure findxy (sx1, sx2, sy1, sy2, x, y, xin, yin) pointer sx1, sx2 # i: GSURFIT pointers (X) pointer sy1, sy2 # i: GSURFIT pointers (Y) real x, y # i: Known output X,Y values real xin, yin # o: Input x,y values to be found #-- # Functions real gseval() # Function evaluator real gsgetr() # Get a fit parameter real cx,cy # Temp values, lin. coeffs. - slopes in X and Y real ex,ey # More temporary values int i # Loop counter real gx,gy # Guess values real newgx,newgy # Improved guesses real xmin,xmax # Range of fits (X) real ymin,ymax # Range of fits (Y) real lmin,lmax # Temporary values begin # Get the fits range xmin=gsgetr(sx1, GSXMIN) xmax=gsgetr(sx1, GSXMAX) ymin=gsgetr(sy1, GSYMIN) ymax=gsgetr(sy1, GSYMAX) # Get first guesses and the slopes # by assuming approx linear shape lmin = gseval(sx1, xmin, ymin) lmax = gseval(sx1, xmax, ymin) cx=(lmax-lmin)/(xmax-xmin) gx = ((xmax-xmin)*(x-lmin) / (lmax-lmin)) + xmin if(cx<1.0e-10 || lmin==lmax) { xin=INDEFR yin=INDEFR return } # Also in Y... lmin = gseval(sy1, xmin, ymin) lmax = gseval(sy1, xmin, ymax) cy=(lmax-lmin)/(ymax-ymin) gy = ((ymax-ymin)*(y-lmin) / (lmax-lmin)) + ymin # In pathological cases set to undefined if(cy<1.0e-10 || lmin==lmax) { xin=INDEFR yin=INDEFR return } # Iterate until too many tries or error acceptable do i = 1, MAXITER { ex=gseval(sx1, gx, gy) if(sx2 != NULL) ex=ex+gseval(sx2,gx,gy) newgx = gx + (x - ex)/cx ey=gseval(sy1, gx, gy) if(sy2 != NULL) ey=ey+gseval(sy2,gx,gy) newgy = gy + (y - ey)/cy if ((abs(newgx-gx) < OKERROR) && (abs(newgy-gy) < OKERROR)) { xin=gx yin=gy return } gx=newgx gy=newgy } end