include include include "register.h" # t_headers -- compute necessary registration parameters by # using FITS coordinate parameters # # R.L. Williamson II, 3-Jun-92 Task derived from HEADER_SX in old SDAS # # This task allows the user to compute the parameters needed to register # a "Secondary" file to the pixel coordinate scale of a "Reference" file # using FITS coordinate parameters (CRVAL, CRPIX, etc.) contained # in the group parameter block. The coordinate types (CTYPE) in the # Reference and Secondary files must be the same. procedure t_headers() bool rlinear, slinear # Flags to tell if coords are linear bool inpixels # Flag to tell if rectangular coord # are used(i.e. CTYPE ="PIXELS") char outtable[SZ_FNAME] # output table name char reffile[SZ_FNAME] # reference image file char secfile[SZ_FNAME] # secondary image file double rcrpix[2] # reference image reference pixel double scrpix[2] # secondary image reference pixel double rcrval[2] # world coordinates at reference pixel double scrval[2] # world coordinates at reference pixel char rctype[SZ_CTYPE,2] # Ctype matrix char sctype[SZ_CTYPE,2] # Ctype matrix double rcd[2,2], scd[2,2] # CD Matrix double det # Determinant of rcd double adj11, adj12, adj21, adj22 # Values of Adjoint of the rcd * scd 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 double ra double delalfa, deldec, crota, radcrota int rrotaxes[2], srotaxes[2] int rnaxis # number of axes in reference image int snaxis # number of axes in secondary image int naxis # number of axes int status # Status real deltax, deltay, xmag, ymag, theta # output fit parameters real sign pointer sp pointer rip # reference image descriptor pointer sip # secondary image descriptor pointer otp # output table pointer pointer coeff # scratch for fit coefficients pointer out[11] # pointer for columns of output table pointer tbtopn() # Functions pointer immap() bool strne() string typmism "Reference and Secondary files have mismatch in coord.types!" string rcoordtype "I don't know the coordinate type of the reference file!" string scoordtype "I don't know the coordinate type of the secondary file!" data rad_to_deg /57.2957795130823/ begin call smark (sp) # Open Reference File call clgstr ("reffile", reffile, SZ_FNAME) rip = immap (reffile, READ_ONLY, 0) # Get number of dimensions rnaxis = IM_NDIM(rip) # Get CTYPE from header # Get MWCS stuff from the reference file call stx_getcoord (rip, rcrpix, rcrval, rcd, 2, rctype, SZ_CTYPE) # Open Secondary File call clgstr ("secfile", secfile, SZ_FNAME) sip = immap (secfile, READ_ONLY, 0) # Get number of dimensions snaxis = IM_NDIM(sip) # Get MWCS stuff from the secondary file call stx_getcoord (sip, scrpix, scrval, scd, 2, sctype, SZ_CTYPE) # If the two images are not the same dimension, forget it! if (rnaxis != snaxis) { call error(1, "Two images have different dimensionality!") } naxis=rnaxis # Do the 1-dimensionsal case first if (naxis == 1) { if (rcd[1,1] == 0.0d0 || scd[1,1] == 0.0d0) { call error (1, "One of the scale parameters is zero!") } else if (strne(rctype[1,1],sctype[1,1])) { call error (1, typmism) } else { # all appears well, so calculate coefficients # Allocate some memory for fit coefficients call salloc (coeff, 2, TY_DOUBLE) Memd[coeff+1] = scd[1,1]/rcd[1,1] Memd[coeff] = rcrpix[1] + (scrval[1] - rcrval[1] - scrpix[1]*scd[1,1])/ rcd[1,1] } # Place Values into array parmval[1] = Memd[coeff] parmval[2] = Memd[coeff+1] parmval[3] = Memd[coeff] parmval[4] = Memd[coeff+1] # Write output table call clgstr ("outtable", outtable, SZ_FNAME) otp = tbtopn (outtable, NEW_FILE, 0) call salloc (out,4,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], "deltax", "", "", TY_DOUBLE, 1, 1) call tbcdef (otp, out[4], "xmag", "", "", TY_DOUBLE, 1, 1) call tbtcre (otp) # Dump the parameters to the table call tbrptd (otp, out, parmval, 4,1) call tbtclo (otp) } else { #2 dimensions inpixels = false call crdtyp (naxis,rctype,rrotaxes,rlinear, status) if (status != OK) { if (status == RECT) { inpixels = true } else { call error (1, rcoordtype) } } ra = rrotaxes[1] call crdtyp (naxis,sctype,srotaxes,slinear, status) if (status != OK) { if (status != RECT) { call error (1, scoordtype) } } # Find the transformation coefficients of the form: # a + bx + cy, d + ex + fy # Check for rectangular versus spherical if (inpixels) { delalfa = 0.0 deldec = 0.0 } else { if (ra == 1) { # Cosine declination correction delalfa = (scrval[1] - rcrval[1]) * cos( rcrval[2] /rad_to_deg) deldec = (scrval[2] - rcrval[2]) } else { delalfa = (scrval[1] - rcrval[1]) deldec = (scrval[2] - rcrval[2]) * cos( rcrval[1]/rad_to_deg) } } det = rcd[1,1]*rcd[2,2] - rcd[1,2] * rcd[2,1] adj11 = rcd[2,2]*scd[1,1] - rcd[1,2]*scd[2,1] adj21 = rcd[1,1]*scd[2,1] - rcd[2,1]*scd[1,1] adj12 = rcd[2,2]*scd[1,2] - rcd[1,2]*scd[2,2] adj22 = rcd[1,1]*scd[2,2] - rcd[2,1]*scd[1,2] # Allocate some memory for fit coefficients call salloc (coeff,6,TY_DOUBLE) Memd[coeff+1] = adj11 / det Memd[coeff+2] = adj12 / det Memd[coeff+4] = adj21 / det Memd[coeff+5] = adj22 / det xcentr = scrpix[1] ycentr = scrpix[2] if (det < 0.0) { sign = -1.0 } else { sign = 1.0 } crota = atan2 (sign * rcd[1,2],rcd[2,2]) radcrota = crota / rad_to_deg Memd[coeff] = rcrpix[1] + (delalfa * cos(radcrota) + deldec * sin(radcrota)) / rcd[1,1] - Memd[coeff+1]*scrpix[1] - Memd[coeff+2]*scrpix[2] Memd[coeff+3] = rcrpix[2] + (-delalfa * sin(radcrota) + deldec * cos(radcrota)) / rcd[2,2] - Memd[coeff+4]*scrpix[1] - Memd[coeff+5]*scrpix[2] 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 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