include include # t_reflresamp -- resamples image data for a given amount of translation, #rotation (images only), and magnification of the science data based on #parameters determined by the register package. Also can resample using #user defined resampling parameters and can reflect data horizontally, #verically, or both ways. # # R.L. Williamson II, 14-Aug-92 Task derived from RESAMP_SX in old SDAS # # This routine resamples time series, spectral, or image data for a given # amount of translation, rotation, and/or magnification of the data in the # data window. This includes image rotation, scale changes, reduction or # magnification, and reflection. For rotation, the axis can be # cursor-designated and different algorithms will be used for large and # for small angle rotations. Interpolation techniques include nearest # neighbor, n-th order polynomial, cubic spline, and Aiken- Lagrange. # Scale changes require user to input magnification (or demagnification) # factor. In the case of image data, independent scale # changes can be done for each axis. # # H.A. Bushouse, 25-Jun-98 Changed "INDEF" to "INDEFD" for v2.11 # procedure t_reflresamp() bool nullflag[6] # null pointer bool mask # mask file logical bool regrsp # registration table logical char infile[SZ_FNAME] # input image file name char outfile[SZ_FNAME] # output image file name char reffile[SZ_FNAME] # reference image file name char maskfile[SZ_FNAME] # input mask file name char regtable[SZ_FNAME] # input registration table char tem[SZ_FNAME] # place holder double coeff[6] # scratch for fit coefficients double determ # determinant of coefficient matrix double incoeff[6] # scratch for inverse transfm coeff double xcentr, ycentr # x,y of center of secondary image double rad_to_deg # convert from radians to degrees double crpix[2] # image reference pixel double crval[2] # world coordinates at reference pixel char ctype[SZ_CTYPE,2] # Ctype matrix double cd[2,2] # CD Matrix int arrsz # Size of array of column names int option # Resample or reflection? int rflopt # reflection option int intp # which interpolation to use? int rnaxis # number of axes in reference image int naxis # number of axes in secondary image int numcoef # Number of coefficients (2 or 6) int dim[2] # Dimensions of input image int outdim[2] # Dimensions of output image int i # Counter int degree # # degrees of freedom in fit real deltax, deltay, xmag, ymag, theta # output fit parameters real theta_rad,costh,sinth real badpix # Window mask value pointer sp pointer rp # reference image descr]tor pointer ip # input image descriptor pointer op # output image descriptor pointer mp # mask image descriptor pointer bpt # input data array pointer pointer mpt # mask data array pointer pointer opt # output data array pointer pointer rtp # input table pointer pointer cp # pointer for column name strings pointer rcolptr[6] # pointers to columns of input table pointer msg pointer str string dimmism "Dimensional mismatch between input and reference file" string cmsg " C(%d) = %15.8e \n" string lmsg1 " The output file will be %d points long\n " string lmsg2 " The output file will be a %d X %d image\n" int clgeti() #Functions int clgwrd() real clgetr() double clgetd() pointer tbtopn() pointer immap() pointer mwr pointer mw pointer mw_open() pointer mw_openim() pointer imgs1r(),imgs2r(),imps1r(), imps2r() bool streq() bool isblank errchk mw_gwattrs() data rad_to_deg/57.2957795130823/ begin call smark (sp) regrsp = false mask = false # Get the name of the input files from the parameter file call clgstr ("infile", infile, SZ_FNAME) call clgstr ("outfile", outfile, SZ_FNAME) call clgstr ("reffile", reffile, SZ_FNAME) call clgstr ("maskfile", maskfile, SZ_FNAME) call clgstr ("regtable", regtable, SZ_FNAME) if (streq(infile,outfile)) call error (1, "Output file cannot have same name as input!") ip = immap (infile, READ_ONLY, 0) if (!isblank(maskfile)) mask = true if (!isblank(regtable)) regrsp = true # Get number of dimensions naxis = IM_NDIM(ip) if (naxis != 1 && naxis != 2) call error(1, "input file not vector or image!") do i = 1, naxis { dim[i] = IM_LEN(ip, i) } if (naxis == 1) dim[2] = 1 if (mask) { mp = immap(maskfile, READ_ONLY, 0) } # Get relection or resample option option = clgwrd ("option", tem, SZ_LINE, "|resample|reflect") # Task options # 1 = resample # 2 = reflect if (option == 0) { call error (1, "Illegal task option") } # If option is resample, set reflection option to none if (option == 1) rflopt = 1 # If relection is chosen, get the reflection option if (option == 2) { rflopt = clgwrd ("rflopt", tem, SZ_LINE, "|none|horizontal|vertical|both") # Task options # 1 = none # 2 = horizontal # 3 = vertical [2-D only] # 4 = both [2-D only] if (rflopt == 0) { call error (1, "Illegal reflection option") } if (rflopt == 3 && naxis != 2) { call error (1, "Cannot do vertical reflection on a vector!") } if (rflopt == 4 && naxis != 2) { call error (1, "Both option is only for 2-D images!") } call aclrd (incoeff,6) #Initialize inverse coefficients do i = 1, naxis { outdim[i] = dim[i] } } # Get the interpolation option (resample option only) if (option == 1) { intp = clgwrd ("inttyp", tem, SZ_LINE, "|nearneig|polynom|4point|") # Fit options # 1 = nearest neighbor # 2 = polynomial fit # 3 = 4 point interpolation if (intp == 0) { call error (1, "Illegal interpolation option") } if (intp == 3 && naxis != 2) { call error (1, "Cannot do 4-point interpolation on a vector!") } if (intp == 2) { degree = clgeti ("degree") } # If using a registration table, need to open reference image, # open registration table and get the fit # coefficients. call aclrd (coeff,6) #Initialize coefficients if (regrsp) { # If there is a registration table #Open reference image rp = immap (reffile, READ_ONLY, 0) # Get number of dimensions rnaxis = IM_NDIM(rp) # Test to see if dimensionality of images the same if (naxis != rnaxis) { call error (1, dimmism) } do i = 1, naxis { outdim[i] = IM_LEN(rp, i) } if (naxis == 1) outdim[2] = 1 numcoef = 6 if (naxis == 1) numcoef = 2 # Open registration table rtp = tbtopn (regtable, READ_ONLY, 0) arrsz = (SZ_COLNAME + 1) * 6 call malloc (cp, arrsz, TY_CHAR) call strcpy("coeff1", Memc[cp], SZ_COLNAME) call strcpy("coeff2", Memc[cp+SZ_COLNAME+1], SZ_COLNAME) call strcpy("coeff3", Memc[cp+2*(SZ_COLNAME+1)],SZ_COLNAME) call strcpy("coeff4", Memc[cp+3*(SZ_COLNAME+1)],SZ_COLNAME) call strcpy("coeff5", Memc[cp+4*(SZ_COLNAME+1)],SZ_COLNAME) call strcpy("coeff6", Memc[cp+5*(SZ_COLNAME+1)],SZ_COLNAME) call tbcfnd (rtp, Memc[cp], rcolptr, numcoef) do i=1,numcoef { if (rcolptr[i] == NULL) call error (1, "All coefficient columns not present!") } call tbrgtd (rtp, rcolptr, coeff, nullflag, numcoef, 1) #Can't do fit if any of the coefficients is missing, check do i=1,numcoef { if (nullflag[i]) { call error(1, "All coefficients not present in registration file!") } # Close the table call tbtclo(rtp) } } else { # User defined parameters used numcoef = 2 deltax = clgetr ("deltax") xmag = clgetr ("xmag") coeff[1] = deltax coeff[2] = xmag # 2-D case if (naxis ==2) { numcoef = 6 deltay = clgetr ("deltay") theta = clgetr ("theta") ymag = clgetr ("ymag") do i = 1, naxis { outdim[i] = dim[i] } # Get the central pixel coords xcentr = clgetd ("xcentr") ycentr = clgetd ("ycentr") # If these values are INDEF, then set these to the # center pixel of the input image. if (xcentr == INDEFD) { xcentr = (dim[1] + 1) /2 } if (ycentr == INDEFD) { ycentr = (dim[2] + 1) /2 } # Now get coordinate transform coefficients theta_rad = theta/rad_to_deg costh = cos (theta_rad) sinth = sin (theta_rad) coeff[2] = xmag*costh coeff[3] = -ymag*sinth coeff[1] = xcentr + deltax - coeff[2]*xcentr - coeff[3]*ycentr coeff[5] = xmag*sinth coeff[6] = ymag*costh coeff[4] = ycentr + deltay - coeff[5]*xcentr - coeff[6]*ycentr } } # Set up the inverse transform if (naxis == 1) { incoeff[1] = -coeff[1] / coeff[2] incoeff[2] = 1.0d0/coeff[2] } else { determ = coeff[2]*coeff[6] - coeff[3]*coeff[5] if (abs(determ) < 1.E-6) { call error(1,"Fatal:Coefficients transform is singular") } incoeff[2] = coeff[6] / determ incoeff[3] = -coeff[3] / determ # See E. Kreyszig, Advanced incoeff[5] = -coeff[5] / determ # Engineering Mathematics incoeff[6] = coeff[2] / determ incoeff[1] = -(incoeff[2]*coeff[1] + incoeff[3]*coeff[4]) incoeff[4] = -(incoeff[5]*coeff[1] + incoeff[6]*coeff[4]) } # Finally tell the user what the coordinate transform # coefficients are and how big the resulting output will be # If he wants now is the time to abort the mess! call eprintf(" \n") call eprintf(" \n") call eprintf(" \n") call eprintf(" Through your selections the following\n") call eprintf(" coordinate transform coefficients\n") call eprintf(" will be used to resample your input\n") if (naxis == 1) { call eprintf(" Xout = C(1) + C(2) * Xin\n") } else { call eprintf(" Xout = C(1) + C(2) * Xin + C(3) * Yin\n") call eprintf(" Yout = C(4) + C(5) * Xin + C(6) * Yin\n") } # Allocate message array call salloc(msg, SZ_LINE, TY_CHAR) do i=1,numcoef { call sprintf(Memc[msg], SZ_LINE, cmsg) call pargi(i) call pargd(coeff[i]) call eprintf(Memc[msg]) } if (naxis == 1){ call sprintf(Memc[msg], SZ_LINE, lmsg1) call pargi(outdim[1]) call eprintf(Memc[msg]) } else { call sprintf(Memc[msg], SZ_LINE, lmsg2) call pargi(outdim[1]) call pargi(outdim[2]) call eprintf(Memc[msg]) } } # Open the output file op = immap (outfile, NEW_COPY, ip) IM_LEN (op,1) = outdim[1] IM_LEN (op,2) = outdim[2] # If regrsp option is chosen, want to copy WCS of reference image, # so get them, and add them to the output file. if (regrsp){ call stx_getcoord (rp, crpix, crval, cd, 2, ctype, SZ_CTYPE) # Set the MWCS info to an identity with dimension equal to # that of the image. mw = mw_open (NULL, naxis) #Set axis type and attribute to that of reference image mwr = mw_openim(rp) #Allocate String pointer call salloc (str, SZ_LINE, TY_CHAR) iferr (call mw_gwattrs (mwr, 0, "system", Memc[str], SZ_LINE)) Memc[str] = EOS call mw_swattrs (mw, 0, "system", Memc[str]) do i = 1, naxis { iferr (call mw_gwattrs (mwr, i, "wtype", Memc[str], SZ_LINE)) Memc[str] = EOS call mw_swattrs (mw, i, "wtype", Memc[str]) iferr (call mw_gwattrs (mwr, i, "axtype", Memc[str], SZ_LINE)) Memc[str] = EOS call mw_swattrs (mw, i, "axtype", Memc[str]) iferr (call mw_gwattrs (mwr, i, "units", Memc[str], SZ_LINE)) Memc[str] = EOS call mw_swattrs (mw, i, "units", Memc[str]) iferr (call mw_gwattrs (mwr, i, "label", Memc[str], SZ_LINE)) Memc[str] = EOS call mw_swattrs (mw, i, "label", Memc[str]) iferr (call mw_gwattrs (mwr, i, "format", Memc[str], SZ_LINE)) Memc[str] = EOS call mw_swattrs (mw, i, "format", Memc[str]) } call mw_saveim (mw, op) call mw_close (mw) call mw_close (mwr) #Put in new values for WCS stuff call imaddd (op, "cd1_1", cd[1,1]) call imaddd (op, "cd1_2", cd[1,2]) call imaddd (op, "cd2_1", cd[2,1]) call imaddd (op, "cd2_2", cd[2,2]) call imaddd (op, "crpix1", crpix[1]) call imaddd (op, "crpix2", crpix[2]) call imaddd (op, "crval1", crval[1]) call imaddd (op, "crval2", crval[2]) call imastr (op, "ctype1", ctype[1,1]) call imastr (op, "ctype2", ctype[1,2]) } # Convert output file to REAL IM_PIXTYPE(op) = TY_REAL # Grab the fill value from the par file (badpix) badpix = clgetr("badpix") # Now call either the vector or image resampling subroutines if(mask) { # Is there a data mask? if (naxis == 1) { # Vector? bpt = imgs1r (ip, 1, dim[1]) mpt = imgs1r (mp, 1, dim[1]) opt = imps1r (op, 1, outdim[1]) call vcresw (Memr[bpt], dim[1],intp, rflopt, incoeff, degree, Memr[mpt], badpix, Memr[opt], outdim[1]) } else { # Image bpt = imgs2r (ip, 1, dim[1],1,dim[2]) mpt = imgs2r (mp, 1, dim[1],1,dim[2]) opt = imps2r (op, 1, outdim[1],1,outdim[2]) call imresw (Memr[bpt], dim[1],dim[2],intp, rflopt, incoeff,degree, Memr[mpt], badpix, Memr[opt],outdim[1],outdim[2]) } } else { # No mask if (naxis == 1) { # Vector? bpt = imgs1r (ip, 1, dim[1]) opt = imps1r (op, 1, outdim[1]) call vecres (Memr[bpt], dim[1],intp, rflopt, incoeff, degree, badpix,Memr[opt], outdim[1]) } else { # Image bpt = imgs2r (ip, 1, dim[1],1,dim[2]) opt = imps2r (op, 1, outdim[1],1,outdim[2]) call imgres (Memr[bpt], dim[1],dim[2],intp, rflopt, incoeff,degree, badpix, Memr[opt],outdim[1],outdim[2]) } } # if the reflection option is chosen, need to update the MWCS to # "reflect" the changes if (option == 2) { call stx_getcoord (op, crpix, crval, cd, 2, ctype, SZ_CTYPE) # 1 = none # 2 = horizontal # 3 = vertical [2-D only] # 4 = both [2-D only] if (rflopt == 2 || rflopt == 4) { call imaddd (op, "cd1_1", -cd[1,1]) if (naxis == 2) { call imaddd (op, "cd2_1", -cd[2,1]) } } if (rflopt == 3 || rflopt == 4) { call imaddd (op, "cd1_2", -cd[1,2]) call imaddd (op, "cd2_2", -cd[2,2]) } } # close the images call imunmap(ip) if (regrsp) call imunmap(rp) call imunmap(op) call sfree(sp) end