include include # t_imagfe -- compute the necessary parameters to register an image by # using common features of image data # # R.L. Williamson II, 29-Apr-92 Task derived from REGIST_SX in old SDAS # # Description: # ------------ # This routine determines the amount of translation, rotation, and/or # magnification needed to make two or more images congruent. # The determination involves the use of reference points (specified # by cursor input or by time, wavelength, or position in a standard # coordinate system). # # Usage Notes: # ------------ # Accuracy of input reference point positions can be improved by using # the task RIMCURSOR's output linked to this routine in a command # language procedure. (In fact, this will probably be the standard way of # using this module). # # A resampled image may be produced by linking the output of Register # Files to Resample in a command language procedure. procedure t_imagfe() char reftable[SZ_FNAME] # reference input table name char sectable[SZ_FNAME] # second input table name char outtable[SZ_FNAME] # output table name char secfile[SZ_FNAME] # secondary image file char temp[SZ_FNAME] # place holder char colname1[SZ_COLNAME] # input column names char colname2[SZ_COLNAME] char ctype[SZ_CTYPE,2] # Ctype matrix double crval[2] # world coordinates at reference pixel double crpix[2] # Reference pixels double cd[2,2] # CD Matrix double xcentr, ycentr # x,y of center of secondary image double rad_to_deg # convert from radians to degrees double parmval[11] # output parameters of registration int ifit # which fit to use? int rnfeat, snfeat,mxfea # number of features in each table, max int nfeat # number of features comparing int maxit # maximum nmber of iterations int nreject # number of points rejected in fit int naxis # number of axes in secondary image int i # Counter int cnt # Counter real gamma, epsilon # input tolerances used by magrotran real deltax, deltay, xmag, ymag, theta # output fit parameters pointer sp pointer rp # reference table descriptor pointer secp # secondary table descriptor pointer sip # secondary image descriptor pointer otp # output table pointer pointer rcolptr[2], scolptr[2] # pointers to columns of input tables pointer nullflag1 # null pointer pointer nullflag2 # null pointer pointer nullflag3 # null pointer pointer nullflag4 # null pointer pointer tempx1,tempx2,tempy1,tempy2 # pointers to temporary arrays pointer refx,secx,refy,secy # pointers to arrays of input coordinates pointer coeff # scratch for fit coefficients pointer out[11] # pointer for columns of output table int tbpsta() #Functions int clgeti() int clgwrd() real clgetr() pointer tbtopn() pointer immap() data rad_to_deg/57.2957795130823/ begin call smark (sp) #Read the method of registration to use ifit = clgwrd ("method", temp, SZ_LINE, "|translat|rotrans|constrained|general|") if (ifit == 0) { call error (1, "Illegal method") } #Get size of the reference group call clgstr ("reftable", reftable, SZ_FNAME) rp = tbtopn (reftable, READ_ONLY, 0) rnfeat = tbpsta(rp, TBL_NROWS) #Now get the secondary call clgstr ("sectable", sectable, SZ_FNAME) secp = tbtopn (sectable, READ_ONLY, 0) snfeat = tbpsta(secp, TBL_NROWS) nfeat = min (rnfeat, snfeat) mxfea = max (rnfeat,snfeat) # Read the names of the table columns from the par file call clgstr ("column1", colname1, SZ_COLNAME) call clgstr ("column2", colname2, SZ_COLNAME) # Read the data from the tables call salloc (nullflag1, mxfea, TY_BOOL) call salloc (nullflag2, mxfea, TY_BOOL) call salloc (nullflag3, mxfea, TY_BOOL) call salloc (nullflag4, mxfea, TY_BOOL) call salloc (tempx1, rnfeat, TY_DOUBLE) call salloc (tempy1, rnfeat, TY_DOUBLE) call salloc (tempx2, snfeat, TY_DOUBLE) call salloc (tempy2, snfeat, TY_DOUBLE) call salloc (refx, rnfeat, TY_DOUBLE) call salloc (refy, rnfeat, TY_DOUBLE) call salloc (secx, snfeat, TY_DOUBLE) call salloc (secy, snfeat, TY_DOUBLE) call tbcfnd (rp, colname1, rcolptr[1], 1) call tbcfnd (rp, colname2, rcolptr[2], 1) call tbcgtd (rp, rcolptr[1], Memd[tempx1], Memb[nullflag1], 1, rnfeat) call tbcgtd (rp, rcolptr[2], Memd[tempy1], Memb[nullflag2], 1, rnfeat) # Weed out nulls cnt=0 do i=0,rnfeat-1 { if(!Memb[nullflag1+i] || !Memb[nullflag2+i]) { Memd[refx+cnt]=Memd[tempx1+i] Memd[refy+cnt]=Memd[tempy1+i] cnt=cnt+1 } } rnfeat=cnt call tbcfnd (secp, colname1, scolptr[1], 1) call tbcfnd (secp, colname2, scolptr[2], 1) call tbcgtd (secp, scolptr[1], Memd[tempx2], Memb[nullflag3], 1, snfeat) call tbcgtd (secp, scolptr[2], Memd[tempy2], Memb[nullflag4], 1, snfeat) # Weed out nulls cnt=0 do i=0,snfeat-1 { if(!Memb[nullflag1+i] || !Memb[nullflag2+i]) { Memd[secx+cnt]=Memd[tempx2+i] Memd[secy+cnt]=Memd[tempy2+i] cnt=cnt+1 } } snfeat=cnt # Close input tables call tbtclo (rp) call tbtclo (secp) # Allocate some memory for fit coefficients call salloc (coeff,6,TY_DOUBLE) # Get some parameters from the cl needed by MAGROTRAN maxit = clgeti ("maxit") gamma = clgetr ("gamma") epsilon = clgetr ("fittoler") # Call the MAGROTRAN to do the calculations call magrotran ( Memd[secx], Memd[secy], Memd[refx], Memd[refy], nfeat, ifit, maxit, gamma, epsilon, Memd[coeff], nreject) # Compute user parameters # Need to get info from secondary image call clgstr ("secfile", secfile, SZ_FNAME) sip = immap (secfile, READ_ONLY, 0) # Get number of dimensions naxis = IM_NDIM(sip) # Get MWCS stuff from the secondary file call stx_getcoord (sip, crpix, crval, cd, naxis, ctype, SZ_CTYPE) call imunmap (sip) if (naxis == 1) { deltax = Memd[coeff] xmag = Memd[coeff+1] } else { # Assign the values xcentr = crpix[1] ycentr = crpix[2] # Watch out for reflection in coordinate transform! if (Memd[coeff+1]*Memd[coeff+5] < 0.0) { if (Memd[coeff+1] < 0.0 ) { theta = atan2 (-Memd[coeff+2], Memd[coeff+5] ) * rad_to_deg } else { theta = atan2 ( Memd[coeff+4], Memd[coeff+1] ) * rad_to_deg } } else { theta = atan2 ( Memd[coeff+4], Memd[coeff+1] ) * rad_to_deg } xmag = sqrt (Memd[coeff+1]**2 + Memd[coeff+4]**2) ymag = sqrt (Memd[coeff+2]**2 + Memd[coeff+5]**2) deltax = Memd[coeff] + Memd[coeff+1]*xcentr + Memd[coeff+2]*ycentr deltay = Memd[coeff+3] + Memd[coeff+4]*xcentr + Memd[coeff+5]*ycentr deltax = deltax - xcentr # Just give the offset between deltay = deltay - ycentr # reference pixels } # Place Values into array parmval[1] = Memd[coeff] parmval[2] = Memd[coeff+1] parmval[3] = Memd[coeff+2] parmval[4] = Memd[coeff+3] parmval[5] = Memd[coeff+4] parmval[6] = Memd[coeff+5] parmval[7] = deltax parmval[8] = deltay parmval[9] = xmag parmval[10] = ymag parmval[11] = theta # Write output from magrotran into table call clgstr ("outtable", outtable, SZ_FNAME) otp = tbtopn (outtable, NEW_FILE, 0) call salloc (out,11,TY_POINTER) call tbcdef (otp, out[1], "coeff1", "", "", TY_DOUBLE, 1, 1) call tbcdef (otp, out[2], "coeff2", "", "", TY_DOUBLE, 1, 1) call tbcdef (otp, out[3], "coeff3", "", "", TY_DOUBLE, 1, 1) call tbcdef (otp, out[4], "coeff4", "", "", TY_DOUBLE, 1, 1) call tbcdef (otp, out[5], "coeff5", "", "", TY_DOUBLE, 1, 1) call tbcdef (otp, out[6], "coeff6", "", "", TY_DOUBLE, 1, 1) call tbcdef (otp, out[7], "deltax", "", "", TY_DOUBLE, 1, 1) call tbcdef (otp, out[8], "deltay", "", "", TY_DOUBLE, 1, 1) call tbcdef (otp, out[9], "xmag", "", "", TY_DOUBLE, 1, 1) call tbcdef (otp, out[10], "ymag", "", "", TY_DOUBLE, 1, 1) call tbcdef (otp, out[11], "rotang", "", "", TY_DOUBLE, 1, 1) call tbtcre (otp) # Dump the parameters to the table call tbrptd (otp, out, parmval, 11,1) call tbtclo (otp) call sfree (sp) end