diff options
Diffstat (limited to 'noao/twodspec/longslit/transform')
27 files changed, 3755 insertions, 0 deletions
diff --git a/noao/twodspec/longslit/transform/Notes b/noao/twodspec/longslit/transform/Notes new file mode 100644 index 00000000..16f5a7a3 --- /dev/null +++ b/noao/twodspec/longslit/transform/Notes @@ -0,0 +1,6 @@ +May 29, 1987 + +If a user accidentally leaves the user coordinate as INDEF in tracing +the spatial distortion then FITCOORDS uses the fitted coordinate +which is the same as the pixel coordinate. This causes incorrect +results. Some thought should be given to this situation. diff --git a/noao/twodspec/longslit/transform/fcdbio.x b/noao/twodspec/longslit/transform/fcdbio.x new file mode 100644 index 00000000..caf4ac5d --- /dev/null +++ b/noao/twodspec/longslit/transform/fcdbio.x @@ -0,0 +1,99 @@ +include <error.h> +include <math/gsurfit.h> +include <pkg/dttext.h> +include <units.h> + +# FC_DBWRITE -- Write an fitcoords database entry. + +procedure fc_dbwrite (database, fitname, axis, un, sf) + +char database[ARB] # Database +char fitname[ARB] # Database fit name +int axis # Axis for surface +pointer un # Units pointer +pointer sf # Surface pointer + +int i, nsave +pointer dt, coeffs, sp, dbfile + +int xgsgeti() +pointer dtmap1() + +begin + if (sf == NULL) + return + + call smark (sp) + call salloc (dbfile, SZ_FNAME, TY_CHAR) + call strcpy ("fc", Memc[dbfile], SZ_FNAME) + call imgcluster (fitname, Memc[dbfile+2], SZ_FNAME-2) + dt = dtmap1 (database, Memc[dbfile], APPEND) + + call dtptime (dt) + call dtput (dt, "begin\t%s\n") + call pargstr (fitname) + call dtput (dt, "\ttask\tfitcoords\n") + call dtput (dt, "\taxis\t%d\n") + call pargi (axis) + if (un != NULL) { + call dtput (dt, "\tunits\t%s\n") + call pargstr (UN_UNITS(un)) + } + + nsave = xgsgeti (sf, GSNSAVE) + call salloc (coeffs, nsave, TY_DOUBLE) + call xgssave (sf, Memd[coeffs]) + call dtput (dt, "\tsurface\t%d\n") + call pargi (nsave) + do i = 1, nsave { + call dtput (dt, "\t\t%g\n") + call pargd (Memd[coeffs+i-1]) + } + + call sfree (sp) + call dtunmap (dt) +end + + +# LM_DBREAD -- Read an lsmap database entry. + +procedure lm_dbread (database, fitname, axis, un, sf) + +char database[ARB] # Database +char fitname[ARB] # Fit name +int axis # Axis for surface +pointer un # Units pointer +pointer sf # Surface pointer + +int rec, ncoeffs +pointer dt, coeffs, sp, dbfile, units + +int dtlocate(), dtgeti() +pointer dtmap1(), un_open() + +errchk dtlocate(), dtgeti(), dtgad(), un_open() + +begin + un = NULL + sf = NULL + coeffs = NULL + + call smark (sp) + call salloc (dbfile, SZ_FNAME, TY_CHAR) + call salloc (units, SZ_FNAME, TY_CHAR) + call strcpy ("fc", Memc[dbfile], SZ_FNAME) + call imgcluster (fitname, Memc[dbfile+2], SZ_FNAME-2) + dt = dtmap1 (database, Memc[dbfile], READ_ONLY) + + rec = dtlocate (dt, fitname) + axis = dtgeti (dt, rec, "axis") + ifnoerr (call dtgstr (dt, rec, "units", Memc[units], SZ_FNAME)) + un = un_open (Memc[units]) + ncoeffs = dtgeti (dt, rec, "surface") + call salloc (coeffs, ncoeffs, TY_DOUBLE) + call dtgad (dt, rec, "surface", Memd[coeffs], ncoeffs, ncoeffs) + call xgsrestore (sf, Memd[coeffs]) + + call sfree (sp) + call dtunmap (dt) +end diff --git a/noao/twodspec/longslit/transform/fcdlist.x b/noao/twodspec/longslit/transform/fcdlist.x new file mode 100644 index 00000000..7b9816a7 --- /dev/null +++ b/noao/twodspec/longslit/transform/fcdlist.x @@ -0,0 +1,91 @@ +include <mach.h> +include <error.h> + +# FC_DLIST -- Fit Coordinates Deletion List Procedures. + +# FC_DLREAD -- Fit Coordinates Deletion List Read. +# Read the deletion list file and match points in the list with the data +# and delete them. + +procedure fc_dlread (x, y, w, npts) + +real x[npts] # First coordinate to match +real y[npts] # Second coordinate to match +real w[npts] # Weight of coordinate +int npts # Number of coordinates + +int i, fd +real r +char file[SZ_FNAME] +real xdel, ydel + +int access(), open(), fscan(), nscan() + +begin + call clgstr ("deletions", file, SZ_FNAME) + + if (access (file, READ_ONLY, TEXT_FILE) == NO) + return + + fd = open (file, READ_ONLY, TEXT_FILE) + + while (fscan (fd) != EOF) { + call gargr (xdel) + call gargr (ydel) + + if (nscan() != 2) + next + + do i = 1, npts { + r = sqrt ((x[i]-xdel)**2 + (y[i]-ydel)**2) + if (r < 10*EPSILONR) + w[i] = 0. +# if (x[i] != xdel) +# next +# if (y[i] != ydel) +# next +# w[i] = 0. + } + } + + call close (fd) +end + + +# FC_DLWRITE -- Fit Coordinates Deletion List Write. + +procedure fc_dlwrite (x, y, w, npts) + +real x[npts] # First coordinate to match +real y[npts] # Second coordinate to match +real w[npts] # Weight of coordinate +int npts # Number of coordinates + +int i, fd +char file[SZ_FNAME] + +int open() + +begin + call clgstr ("deletions", file, SZ_FNAME) + + if (file[1] == EOS) + return + + iferr (call delete (file)) + ; + iferr (fd = open (file, NEW_FILE, TEXT_FILE)) { + call erract (EA_WARN) + return + } + + do i = 1, npts { + if (w[i] == 0.) { + call fprintf (fd, "%g %g\n") + call pargr (x[i]) + call pargr (y[i]) + } + } + + call close (fd) +end diff --git a/noao/twodspec/longslit/transform/fcfitcoords.x b/noao/twodspec/longslit/transform/fcfitcoords.x new file mode 100644 index 00000000..13943302 --- /dev/null +++ b/noao/twodspec/longslit/transform/fcfitcoords.x @@ -0,0 +1,211 @@ +include <pkg/gtools.h> +include <pkg/igsfit.h> +include <pkg/xtanswer.h> + +# FC_FITCOORDS -- Fit a surface to the user coordinates. + +procedure fc_fitcoords (fitname, database, list, logfiles, interactive) + +char fitname[SZ_FNAME] # Fitname +char database[SZ_FNAME] # Database +int list # List of images +int logfiles # List of log files +int interactive # Interactive? + +int axis # Axis of surface fit +pointer sf # Surface pointer +char logfile[SZ_FNAME], labels[SZ_LINE, IGSPARAMS] +bool answer +int ncoords, logfd, axes[2] +real xmin, xmax, ymin, ymax +pointer gp, gplog, gt, coords, title, un + +int imtgetim(), fntgfntb(), open(), igs_geti(), scan() +real xgseval() +pointer gopen(), gt_init() + +errchk fc_getcoords + +begin + # Print a header to the log files giving the inputs. This is + # done first so that if one of the logfiles is STDOUT the user + # will see that something is happening. + + axis = 0 + while (fntgfntb (logfiles, logfile, SZ_FNAME) != EOF) { + logfd = open (logfile, APPEND, TEXT_FILE) + call sysid (logfile, SZ_FNAME) + call fprintf (logfd, "\n%s\n") + call pargstr (logfile) + call fprintf (logfd, " Longslit coordinate fit name is %s.\n") + call pargstr (fitname) + call fprintf (logfd, " Longslit database is %s.\n") + call pargstr (database) + call fprintf (logfd, " Features from images:\n") + while (imtgetim (list, logfile, SZ_FNAME) != EOF) { + call fprintf (logfd, " %s\n") + call pargstr (logfile) + } + call imtrew (list) + call close (logfd) + } + call fntrewb (logfiles) + + # Get the coordinates for the specified images and axis. The + # coordinates are returned in a pointer which must be explicitly + # freed. + + call fc_getcoords (database, list, axis, xmin, xmax, ymin, ymax, + coords, ncoords, labels, un) + + # Read points from the deletion list. + + switch (axis) { + case 1: + call fc_dlread (Memr[coords+(Z-1)*ncoords], + Memr[coords+(Y-1)*ncoords], Memr[coords+(W-1)*ncoords], ncoords) + case 2: + call fc_dlread (Memr[coords+(Z-1)*ncoords], + Memr[coords+(X-1)*ncoords], Memr[coords+(W-1)*ncoords], ncoords) + } + + # Initialize the graphics. + + if ((interactive == YES) || (interactive == ALWAYSYES)) { + call clgstr ("graphics", logfile, SZ_FNAME) + gp = gopen (logfile, NEW_FILE, STDGRAPH) + } + + # Set plot log. + + gplog = NULL + call clgstr ("plotfile", logfile, SZ_FNAME) + if (logfile[1] != EOS) { + logfd = open (logfile, APPEND, BINARY_FILE) + gplog = gopen ("stdplot", APPEND, logfd) + } else + gplog = NULL + + gt = gt_init () + call malloc (title, SZ_LINE, TY_CHAR) + call sprintf (Memc[title], SZ_LINE, + "Fit User Coordinates to Image Coordinates for %s") + call pargstr (fitname) + call gt_sets (gt, GTTITLE, Memc[title]) + call mfree (title, TY_CHAR) + + # Fit the surface. The surface is defined over the full range of + # image coordinates. + + call igs_setr (IGS_XMIN, xmin) + call igs_setr (IGS_XMAX, xmax) + call igs_setr (IGS_YMIN, ymin) + call igs_setr (IGS_YMAX, ymax) + + switch (axis) { + case 1: + if (Memr[coords+ncoords-1] == 1) { + axes[1] = Y + axes[2] = R + call igs_fit2 (sf, gp, gplog, gt, axes, Memr[coords], ncoords, + labels, interactive) + } else { + axes[1] = X + axes[2] = R + call igs_fit1 (sf, gp, gplog, gt, axes, Memr[coords], ncoords, + labels, interactive) + } + case 2: + if (Memr[coords+ncoords-1] == 1) { + axes[1] = X + axes[2] = R + call igs_fit3 (sf, gp, gplog, gt, axes, Memr[coords], ncoords, + labels, interactive) + } else { + axes[1] = Y + axes[2] = R + call igs_fit1 (sf, gp, gplog, gt, axes, Memr[coords], ncoords, + labels, interactive) + } + } + + # Close graphics. + + if (gp != NULL) + call gclose (gp) + if (gplog != NULL) { + call gclose (gplog) + call close (logfd) + } + call gt_free (gt) + + # Print logs. + + while (fntgfntb (logfiles, logfile, SZ_FNAME) != EOF) { + logfd = open (logfile, APPEND, TEXT_FILE) + call fprintf (logfd, + " Map %s coordinates for axis %d using image features:\n") + call pargstr (labels[1, Z]) + call pargi (axis) + call fprintf (logfd, " Number of feature coordnates = %d\n") + call pargi (ncoords) + call igs_gets (IGS_FUNCTION, logfile, SZ_FNAME) + call fprintf (logfd, " Mapping function = %s\n") + call pargstr (logfile) + call fprintf (logfd, " X order = %d\n Y order = %d\n") + call pargi (igs_geti (IGS_XORDER)) + call pargi (igs_geti (IGS_YORDER)) + call fprintf (logfd, + " Fitted coordinates at the corners of the images:\n") + call fprintf (logfd, " (%d, %d) = %g (%d, %d) = %g\n") + call pargr (xmin) + call pargr (ymin) + call pargr (xgseval (sf, xmin, ymin)) + call pargr (xmax) + call pargr (ymin) + call pargr (xgseval (sf, xmax, xmin)) + call fprintf (logfd, " (%d, %d) = %g (%d, %d) = %g\n") + call pargr (xmin) + call pargr (ymax) + call pargr (xgseval (sf, xmin, ymax)) + call pargr (xmax) + call pargr (ymax) + call pargr (xgseval (sf, xmax, ymax)) + call close (logfd) + } + call fntrewb (logfiles) + + # Write the fit to the database. + + answer = true + if ((interactive == YES) || (interactive == ALWAYSYES)) { + call printf ("Write coordinate map to the database (yes)? ") + call flush (STDOUT) + if (scan() != EOF) + call gargb (answer) + } + if (answer) + call fc_dbwrite (database, fitname, axis, un, sf) + + # Write list of deleted points. + + if ((interactive == YES) || (interactive == ALWAYSYES)) { + switch (axis) { + case 1: + call fc_dlwrite (Memr[coords+(Z-1)*ncoords], + Memr[coords+(Y-1)*ncoords], + Memr[coords+(W-1)*ncoords], ncoords) + case 2: + call fc_dlwrite (Memr[coords+(Z-1)*ncoords], + Memr[coords+(X-1)*ncoords], + Memr[coords+(W-1)*ncoords], ncoords) + } + } + + # Free memory. + + call mfree (coords, TY_REAL) + if (un != NULL) + call un_close (un) + call xgsfree (sf) +end diff --git a/noao/twodspec/longslit/transform/fcgetcoords.x b/noao/twodspec/longslit/transform/fcgetcoords.x new file mode 100644 index 00000000..dda1c0f0 --- /dev/null +++ b/noao/twodspec/longslit/transform/fcgetcoords.x @@ -0,0 +1,212 @@ +include <imio.h> +include <mach.h> +include <mwset.h> +include <pkg/dttext.h> +include <pkg/igsfit.h> + +# FC_GETCOORDS -- Get feature coordinates for the specified axis and list +# of images. Determine the image dimensions. + +procedure fc_getcoords (database, list, axis, xmin, xmax, ymin, ymax, + coords, ncoords, labels, un) + +char database[ARB] # Database +int list # List of images +int axis # Image axis +real xmin, xmax # Image X limits +real ymin, ymax # Image Y limits +pointer coords # Coordinate data pointer +pointer ncoords # Number of coordinate points +char labels[SZ_LINE,IGSPARAMS] # Axis labels +pointer un # Units pointer + +char image1[SZ_FNAME], image2[SZ_FNAME], root[SZ_FNAME], units[SZ_FNAME] +int i, j, rec, index, imin, imax, nfeatures, ntotal +real value, wt, ltm[2,2], ltv[2] +pointer dt, im, mw, ct, x, y, user + +int fc_getim(), dtgeti(), dtscan(), mw_stati() +real mw_c1tranr() +bool strne() +pointer dtmap1(), immap(), mw_openim(), mw_sctran(), un_open() + +errchk dtmap1, dtgstr, immap + +begin + x = NULL + ncoords = 0 + ntotal = 0 + axis = 0 + imin = MAX_INT + imax = -MAX_INT + un = NULL + + while (fc_getim (list, image1, SZ_FNAME) != EOF) { + call strcpy ("id", root, SZ_FNAME) + call imgcluster (image1, root[3], SZ_FNAME-2) + dt = dtmap1 (database, root, READ_ONLY) + do rec = 1, DT_NRECS(dt) { + + iferr (call dtgstr (dt, rec, "task", image2, SZ_FNAME)) + next + if (strne ("identify", image2)) + next + + call dtgstr (dt, rec, "image", image2, SZ_FNAME) + call get_root (image2, root, SZ_FNAME) + if (strne (image1, root)) + next + + # Map the 1D image section and determine the axis, the + # line or column in the 2D image, and the 2D image size. + + im = immap (image2, READ_ONLY, 0) + j = IM_VMAP(im, 1) + switch (j) { + case 1: + index = IM_VOFF (im, 2) + 1 + case 2: + index = IM_VOFF (im, 1) + 1 + } + imin = min (imin, index) + imax = max (imax, index) + + xmin = 1. + xmax = IM_SVLEN (im, 1) + ymin = 1. + ymax = IM_SVLEN (im, 2) + + if (axis == 0) + axis = j + + if (j != axis) { + call imunmap (im) + call eprintf ( + "Warning: Fit axes don't agree for combine option. Ignoring %s.\n") + call pargstr (image1) + break + } + + # Set the WCS to convert the feature positions from + # IDENTIFY/REIDENTIFY which are in "physical" coordinates + # to "logical" coordinates currently used by TRANSFORM. + + mw = mw_openim (im) + call mw_seti (mw, MW_USEAXMAP, NO) + i = mw_stati (mw, MW_NPHYSDIM) + call mw_gltermr (mw, ltm, ltv, i) + if (ltm[1,1] == 0. && ltm[2,2] == 0.) { + ltm[1,1] = ltm[2,1] + ltm[2,1] = 0. + ltm[2,2] = ltm[1,2] + ltm[1,2] = 0. + call mw_sltermr (mw, ltm, ltv, i) + } else if (ltm[1,2] != 0. || ltm[2,1] != 0.) { + ltv[1] = 0. + ltv[2] = 0. + ltm[1,1] = 1. + ltm[2,1] = 0. + ltm[2,2] = 1. + ltm[1,2] = 0. + call mw_sltermr (mw, ltm, ltv, i) + } + call mw_seti (mw, MW_USEAXMAP, YES) + ct = mw_sctran (mw, "physical", "logical", 1) + + # Allocate memory for the feature information and read + # the database. + + ifnoerr (call dtgstr (dt, rec, "units", units, SZ_FNAME)) + un = un_open (units) + nfeatures = dtgeti (dt, rec, "features") + if (x == NULL) { + call malloc (x, nfeatures, TY_REAL) + call malloc (y, nfeatures, TY_REAL) + call malloc (user, nfeatures, TY_REAL) + } else { + call realloc (x, ncoords+nfeatures, TY_REAL) + call realloc (y, ncoords+nfeatures, TY_REAL) + call realloc (user, ncoords+nfeatures, TY_REAL) + } + + do i = 1, nfeatures { + j = dtscan (dt) + call gargr (value) + switch (axis) { + case 1: + Memr[x+ncoords] = mw_c1tranr (ct, value) + Memr[y+ncoords] = index + case 2: + Memr[x+ncoords] = index + Memr[y+ncoords] = mw_c1tranr (ct, value) + } + call gargr (value) + call gargr (value) + call gargr (wt) + call gargr (wt) + call gargr (wt) + if (!IS_INDEF (value) && wt > 0.) { + Memr[user+ncoords] = value + ncoords = ncoords + 1 + } + ntotal = ntotal + 1 + } + call mw_close (mw) + call imunmap (im) + } + + # Finish up + call dtunmap (dt) + } + + # Set coordinates. Take error action if no features are found. + + if (ncoords > 0) { + call xt_sort3 (Memr[user], Memr[x], Memr[y], ncoords) + call malloc (coords, ncoords*IGSPARAMS, TY_REAL) + call amovr (Memr[x], Memr[coords+(X-1)*ncoords], ncoords) + call amovr (Memr[y], Memr[coords+(Y-1)*ncoords], ncoords) + call amovr (Memr[user], Memr[coords+(Z-1)*ncoords], ncoords) + call amovkr (1., Memr[coords+(W-1)*ncoords], ncoords) + + call fc_setfeatures (Memr[coords], Memr[coords+(Z-1)*ncoords], + ncoords) + + call strcpy ("X (pixels)", labels[1,X], SZ_LINE) + call strcpy ("Y (pixels)", labels[1,Y], SZ_LINE) + call strcpy ("User", labels[1,Z], SZ_LINE) + call strcpy ("Surface", labels[1,S], SZ_LINE) + call strcpy ("Residuals", labels[1,R], SZ_LINE) + } + + call mfree (x, TY_REAL) + call mfree (y, TY_REAL) + call mfree (user, TY_REAL) + + if (ncoords == 0) { + if (ntotal == 0) + call error (1, "No coordinates found in database") + else + call error (1, "Only INDEF coordinates found in database") + } +end + + +# FC_SETFEATURES -- Set the feature numbers. + +procedure fc_setfeatures (features, user, npts) + +real features[npts] # Feature numbers +real user[npts] # User coordinates +int npts # Number of points + +int i + +begin + features[1] = 1 + do i = 2, npts { + features[i] = features[i-1] + if (user[i] != user[i-1]) + features[i] = features[i] + 1 + } +end diff --git a/noao/twodspec/longslit/transform/fcgetim.x b/noao/twodspec/longslit/transform/fcgetim.x new file mode 100644 index 00000000..e76ba25a --- /dev/null +++ b/noao/twodspec/longslit/transform/fcgetim.x @@ -0,0 +1,32 @@ +# FC_GETIM -- Get next image name with standard image extensions removed. +# This is necessary to avoid having two legal image names refering to the +# same image. + +int procedure fc_getim (list, image, maxchar) + +int list # Image list +char image[maxchar] # Image name +int maxchar # Maximum number of chars in image name + +int i, stat, imtgetim(), strmatch() + +begin + stat = imtgetim (list, image, maxchar) + + if (stat == EOF) + return (stat) + + i = strmatch (image, ".imh") + if (i > 0) { + call strcpy (image[i], image[i-4], maxchar) + return (stat) + } + + i = strmatch (image, ".hhh") + if (i > 0) { + call strcpy (image[i], image[i-4], maxchar) + return (stat) + } + + return (stat) +end diff --git a/noao/twodspec/longslit/transform/fitcoords.x b/noao/twodspec/longslit/transform/fitcoords.x new file mode 100644 index 00000000..e849caf2 --- /dev/null +++ b/noao/twodspec/longslit/transform/fitcoords.x @@ -0,0 +1,83 @@ +include <error.h> +include <pkg/igsfit.h> +include <pkg/xtanswer.h> + +# T_FITCOORDS -- Fit a surface to the coordinates of longslit images. +# +# This is the CL entry for this task. All the real work is done by +# fc_fitcoords. + +procedure t_fitcoords () + +int list1 # Image list +char fitname[SZ_FNAME] # Database name for coordinate fit +char database[SZ_FNAME] # Database +int logfiles # List of log files +bool combine # Combine input data? +int interactive # Interactive? + +char image[SZ_FNAME], prompt[SZ_LINE] +int list2 + +int clgeti(), clpopnu(), imtopen(), fc_getim() +bool clgetb() + +begin + # Get the task parameters. + + call clgstr ("fitname", fitname, SZ_FNAME) + call xt_stripwhite (fitname) + combine = clgetb ("combine") + + if (combine && (fitname[1] == EOS)) + call error (0, "Fit name not specified") + + call clgstr ("images", prompt, SZ_LINE) + list1 = imtopen (prompt) + call clgstr ("database", database, SZ_FNAME) + logfiles = clpopnu ("logfiles") + if (clgetb ("interactive")) + interactive = YES + else + interactive = ALWAYSNO + + # Set the initial surface in the igsfit package. + + call clgstr ("function", prompt, SZ_LINE) + call igs_sets (IGS_FUNCTION, prompt) + call igs_seti (IGS_XORDER, clgeti ("xorder")) + call igs_seti (IGS_YORDER, clgeti ("yorder")) + + # For each fit ask the user whether to do the fit interactively. + # If combining the coordinates from all the images in the + # input list then pass the list directly to fc_fitcoords. + # Otherwise for each image in the list create a second list + # containing just that image. A second list is needed because + # fc_fitcoords expects a list. + + if (combine) { + call sprintf (prompt, SZ_LINE, "Fit interactively") + call xt_answer (prompt, interactive) + call fc_fitcoords (fitname, database, list1, logfiles, interactive) + + } else { + while (fc_getim (list1, image, SZ_FNAME) != EOF) { + list2 = imtopen (image) + call sprintf (prompt, SZ_LINE, "Fit %s interactively") + call pargstr (image) + call xt_answer (prompt, interactive) + call sprintf (prompt, SZ_LINE, "%s%s") + call pargstr (fitname) + call pargstr (image) + iferr (call fc_fitcoords (prompt, database, list2, logfiles, + interactive)) + call erract (EA_WARN) + call imtclose (list2) + } + } + + # Finish up. + + call clpcls (logfiles) + call imtclose (list1) +end diff --git a/noao/twodspec/longslit/transform/igsfit/Revisions b/noao/twodspec/longslit/transform/igsfit/Revisions new file mode 100644 index 00000000..92b36cca --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/Revisions @@ -0,0 +1,42 @@ +.help revisions Jun88 noao.twodspec.longslit.transform.igsfit +.nf + igsfit.x + igsnearest.x + GSCUR was being called with DOUBLE precision values. (12/22/87) + + igsfit.x + igscolon.x + igsget.x + Added colon options to print fit at corners of surface. (8/10/87 Valdes) + + ==== + V2.5 + ==== + +noao$twodspec/longslit/transform/igsfit/*.x + Valdes, February 17, 1987 + 1. GIO changes. + +noao$twodspec/longslit/transform/igsfit/igsfit.x +noao$twodspec/longslit/transform/igsfit/igscolon.x + Valdes, January 16, 1987 + 1. '?' now uses system page facility. + 2. Colon command dictionary and switch modified to use macro definitions. + +noao$twodspec/longslit/transform/igsfit/igsdelete.x +noao$twodspec/longslit/transform/igsfit/igsundelete.x + Valdes, October 16, 1986 + 1. Real line type specified in gseti call changed to integer. + This caused a crash on AOS/IRAF. + +======================================================== + +From Valdes on Feb 7, 1986: + +1. Bug fixed in deleting and undeleting points. +------ +From Valdes on Jan 3, 1986: + +1. Modified IGSFIT to allow zooming on constant x, constant y, constant z, +and constant feature. +.endhelp diff --git a/noao/twodspec/longslit/transform/igsfit/igscolon.x b/noao/twodspec/longslit/transform/igsfit/igscolon.x new file mode 100644 index 00000000..6847974a --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/igscolon.x @@ -0,0 +1,115 @@ +include <gset.h> + +# List of colon commands +define CMDS "|show|function|xorder|yorder|corners|" + +define SHOW 1 # Show parameters +define FUNCTION 2 # Set or show function type +define XORDER 3 # Set or show x order of function +define YORDER 4 # Set or show y order of function +define CORNERS 5 # Show corners + +# IGS_COLON -- Processes colon commands. + +procedure igs_colon (cmdstr, gp, sf) + +char cmdstr[ARB] # Command string +pointer gp # GIO pointer +pointer sf # Surface pointer + +char cmd[SZ_LINE] +int ncmd, ival + +int nscan(), strdic() +real xgseval() + +string funcs "|chebyshev|legendre|" + +include "igsfit.com" + +begin + # Use formated scan to parse the command string. + # The first word is the command and it may be minimum match + # abbreviated with the list of commands. + + call sscan (cmdstr) + call gargwrd (cmd, SZ_LINE) + ncmd = strdic (cmd, cmd, SZ_LINE, CMDS) + + switch (ncmd) { + case SHOW: # :show - Show the values of the fitting parameters. + call gdeactivate (gp, AW_CLEAR) + call printf ("function %s\n") + call pargstr (function) + call printf ("xorder %d\n") + call pargi (xorder) + call printf ("yorder %d\n") + call pargi (yorder) + call printf ("Fitted coordinates at the corners of the images:\n") + call printf (" (%d, %d) = %g (%d, %d) = %g\n") + call pargr (xmin) + call pargr (ymin) + call pargr (xgseval (sf, xmin, ymin)) + call pargr (xmax) + call pargr (ymin) + call pargr (xgseval (sf, xmax, xmin)) + call printf (" (%d, %d) = %g (%d, %d) = %g\n") + call pargr (xmin) + call pargr (ymax) + call pargr (xgseval (sf, xmin, ymax)) + call pargr (xmax) + call pargr (ymax) + call pargr (xgseval (sf, xmax, ymax)) + call printf ("rms %g\n") + call pargr (rms) + call greactivate (gp, AW_PAUSE) + + case FUNCTION: # :function - List or set the fitting function. + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + call printf ("function = %s\n") + call pargstr (function) + } else { + if (strdic (cmd, cmd, SZ_LINE, funcs) > 0) + call strcpy (cmd, function, SZ_LINE) + else + call printf ("Unknown or ambiguous function\n") + } + + case XORDER: # xorder: List or set the function order. + call gargi (ival) + if (nscan() == 1) { + call printf ("xorder %d\n") + call pargi (xorder) + } else if (ival < 2) + call printf ("xorder must be at least 2\n") + else + xorder = ival + + case YORDER: # yorder: List or set the function order. + call gargi (ival) + if (nscan() == 1) { + call printf ("yorder %d\n") + call pargi (yorder) + } else if (ival < 2) + call printf ("yorder must be at least 2\n") + else + yorder = ival + case CORNERS: # corners: List coordinates at corners. + call printf ("(%d,%d)=%g (%d,%d)=%g (%d,%d)=%g (%d,%d)=%g\n") + call pargr (xmin) + call pargr (ymin) + call pargr (xgseval (sf, xmin, ymin)) + call pargr (xmax) + call pargr (ymin) + call pargr (xgseval (sf, xmax, xmin)) + call pargr (xmin) + call pargr (ymax) + call pargr (xgseval (sf, xmin, ymax)) + call pargr (xmax) + call pargr (ymax) + call pargr (xgseval (sf, xmax, ymax)) + default: + call printf ("Unrecognized or ambiguous command\007") + } +end diff --git a/noao/twodspec/longslit/transform/igsfit/igsdelete.x b/noao/twodspec/longslit/transform/igsfit/igsdelete.x new file mode 100644 index 00000000..3de2fb25 --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/igsdelete.x @@ -0,0 +1,103 @@ +include <mach.h> +include <gset.h> +include <pkg/gtools.h> +include <pkg/igsfit.h> + +# IGS_NEARESTD -- Nearest point to delete. + +int procedure igs_nearestd (gp, ztype, refpt, axis, pts, npts, wx, wy, wcs) + +pointer gp # GIO pointer +int ztype # Zoom type +int refpt # Reference point +int axis[2] # Axes +real pts[npts, ARB] # Data points +int npts # Number of data points +real wx, wy # Cursor coordinates +int wcs # WCS + +int i, j, x, y +real r2, r2min, x0, y0 + +begin + x = axis[1] + y = axis[2] + + call gctran (gp, wx, wy, wx, wy, wcs, 0) + r2min = MAX_REAL + j = 0 + + if (IS_INDEFI (ztype)) { + do i = 1, npts { + if (pts[i,W] == 0.) + next + call gctran (gp, pts[i, x], pts[i, y], x0, y0, wcs, 0) + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + } else { + do i = 1, npts { + if ((pts[i,ztype] != pts[refpt,ztype]) || (pts[i,W] == 0.)) + next + call gctran (gp, pts[i, x], pts[i, y], x0, y0, wcs, 0) + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + } + + return (j) +end + +# IGS_DELETE -- Delete points or subsets. + +procedure igs_delete (gp, gt, ztype, refpt, axis, pts, npts, dtype) + +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +int ztype # Zoom type +int refpt # Reference point for deletion +int axis[2] # Axes +real pts[npts, ARB] # Data points +int npts # Number of data points +int dtype # Deletion type + +int i, x, y +real xsize, ysize + +real gt_getr() + +begin + x = axis[1] + y = axis[2] + + xsize = gt_getr (gt, GTXSIZE) + ysize = gt_getr (gt, GTYSIZE) + + switch (dtype) { + case X, Y, Z: + do i = 1, npts { + if (!IS_INDEFI (ztype)) + if (pts[i,ztype] != pts[refpt,ztype]) + next + if (pts[i,dtype] != pts[refpt,dtype]) + next + call gseti (gp, G_PMLTYPE, 0) + call gmark (gp, pts[i,x], pts[i,y], GM_PLUS, xsize, ysize) + call gseti (gp, G_PMLTYPE, 1) + call gmark (gp, pts[i,x], pts[i,y], GM_CROSS, xsize, ysize) + pts[i,W] = 0. + } + default: + call gseti (gp, G_PMLTYPE, 0) + call gmark (gp, pts[refpt,x], pts[refpt,y], GM_PLUS, xsize, ysize) + call gseti (gp, G_PMLTYPE, 1) + call gmark (gp, pts[refpt,x], pts[refpt,y], GM_CROSS, xsize, ysize) + pts[refpt,W] = 0. + } +end diff --git a/noao/twodspec/longslit/transform/igsfit/igsfit.com b/noao/twodspec/longslit/transform/igsfit/igsfit.com new file mode 100644 index 00000000..90bf90aa --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/igsfit.com @@ -0,0 +1,10 @@ +# Common parameters. + +char function[SZ_LINE] # Surface function +int xorder # X order of surface function +int yorder # Y order of surface function +real xmin, xmax # X range +real ymin, ymax # Y range +real mean, rms # Mean and RMS of fit + +common /igscom/ xmin, xmax, ymin, ymax, xorder, yorder, function, mean, rms diff --git a/noao/twodspec/longslit/transform/igsfit/igsfit.x b/noao/twodspec/longslit/transform/igsfit/igsfit.x new file mode 100644 index 00000000..14e8e51e --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/igsfit.x @@ -0,0 +1,373 @@ +include <mach.h> +include <pkg/gtools.h> +include <pkg/igsfit.h> + +define HELP "noao$lib/scr/igsfit.key" +define PROMPT "fitcoords surface fitting options" + + +# IGS_FIT1 -- Fit z = f(x, y) + +procedure igs_fit1 (sf, gp, gplog, gt, axis, pts, npts, labels, interactive) + +pointer sf # GSURFIT pointer +pointer gp # GIO pointer +pointer gplog # GIO pointer for plot log +pointer gt # GTOOLS pointer +int axis[2] # Axis definitions +real pts[npts, ARB] # Data +int npts # Number of pts points +char labels[SZ_LINE, ARB] # Identification labels +int interactive # Interactive? + +extern igs_solve1() + +begin + call igs_fit (sf, gp, gplog, gt, axis, pts, npts, labels, interactive, + igs_solve1) +end + + +# IGS_FIT2 -- Fit z = x + f(y) + +procedure igs_fit2 (sf, gp, gplog, gt, axis, pts, npts, labels, interactive) + +pointer sf # GSURFIT pointer +pointer gp # GIO pointer +pointer gplog # GIO pointer for plot log +pointer gt # GTOOLS pointer +int axis[2] # Axis definitions +real pts[npts, ARB] # Data +int npts # Number of pts points +char labels[SZ_LINE, ARB] # Identification labels +int interactive # Interactive? + +extern igs_solve2() + +begin + call igs_fit (sf, gp, gplog, gt, axis, pts, npts, labels, interactive, + igs_solve2) +end + + +# IGS_FIT3 -- Fit z = y + f(x) + +procedure igs_fit3 (sf, gp, gplog, gt, axis, pts, npts, labels, interactive) + +pointer sf # GSURFIT pointer +pointer gp # GIO pointer +pointer gplog # GIO pointer for plot log +pointer gt # GTOOLS pointer +int axis[2] # Axis definitions +real pts[npts, ARB] # Data +int npts # Number of pts points +char labels[SZ_LINE, ARB] # Identification labels +int interactive # Interactive? + +extern igs_solve3() + +begin + call igs_fit (sf, gp, gplog, gt, axis, pts, npts, labels, interactive, + igs_solve3) +end + + +# IGS_FIT -- Interactive surface fitting. + +procedure igs_fit (sf, gp, gplog, gt, axis, pts, npts, labels, interactive, igs_solve) + +pointer sf # GSURFIT pointer +pointer gp # GIO pointer +pointer gplog # GIO pointer for plot log +pointer gt # GTOOLS pointer +int axis[2] # Axis definitions +real pts[npts, ARB] # Data +int npts # Number of pts points +char labels[SZ_LINE, ARB] # Identification labels +int interactive # Interactive? +extern igs_solve() # Surface solution routine + +int i, newgraph, ztype, dtype, refpt, refpt1 +real zval, zval1 +pointer wts + +real wx, wy +int wcs, key +char cmd[SZ_LINE] + +int clgcur(), gt_gcur(), igs_nearest(), igs_nearestd(), igs_nearestu() +errchk igs_solve + +include "igsfit.com" + +begin + # Compute a solution and set the residuals. + + call igs_solve (sf, pts[1,X], pts[1,Y], pts[1,Z], pts[1,W], npts) + call xgsvector (sf, pts[1,X], pts[1,Y], pts[1,S], npts) + call asubr (pts[1,Z], pts[1,S], pts[1,R], npts) + call aavgr (pts[1,R], npts, mean, rms) + call igs_params (gt) + + # Return if not interactive. + + ztype = INDEFI + if ((gp == NULL) || (interactive == NO)) + goto 30 + + call malloc (wts, npts, TY_REAL) + call amovr (pts[1,W], Memr[wts], npts) + + call igs_graph (gp, gt, ztype, refpt, axis, pts, npts, labels) + newgraph = NO + + # Read cursor commands. + +10 while (gt_gcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) != EOF) { + switch (key) { + case '?': + # Print help text. + + call gpagefile (gp, HELP, PROMPT) + + case ':': + # List or set parameters + + if (cmd[1] == '/') + call gt_colon (cmd, gp, gt, newgraph) + else + call igs_colon (cmd, gp, sf) + + # Set abscissa + + case 'x': + call printf ("Select abscissa (x, y, z, s, r): ") + if (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF) + goto 10 + call printf ("\n") + + switch (key) { + case 'x': + i = X + case 'y': + i = Y + case 'z': + i = Z + case 's': + i = S + case 'r': + i = R + default: + call printf ("\07\n") + goto 10 + } + + if (axis[1] != i) { + axis[1] = i + call gt_setr (gt, GTXMIN, INDEF) + call gt_setr (gt, GTXMAX, INDEF) + } + + # Set ordinate + + case 'y': + call printf ("Select ordinate (x, y, z, s, r): ") + if(clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF) + goto 10 + call printf ("\n") + + switch (key) { + case 'x': + i = X + case 'y': + i = Y + case 'z': + i = Z + case 's': + i = S + case 'r': + i = R + default: + call printf ("\07\n") + goto 10 + } + + if (axis[2] != i) { + axis[2] = i + call gt_setr (gt, GTYMIN, INDEF) + call gt_setr (gt, GTYMAX, INDEF) + } + + case 'r': + newgraph = YES + + case 'z': + if (IS_INDEFI (ztype)) { + refpt = igs_nearest (gp, ztype, refpt, axis, pts, npts, wx, + wy, wcs) + + call printf ("Zoom type (x, y, z): ") + if (clgcur ("cursor",wx,wy,wcs,key,cmd,SZ_LINE) == EOF) + goto 10 + call printf ("\n") + + switch (key) { + case 'x': + ztype = X + case 'y': + ztype = Y + case 'z': + ztype = Z + default: + call printf ("\07\n") + goto 10 + } + + newgraph = YES + } + + case 'p': + if (!IS_INDEFI (ztype)) { + ztype = INDEFI + newgraph = YES + } + + case 'l': + if (!IS_INDEFI (ztype)) { + refpt1 = 0 + zval = pts[refpt, ztype] + zval1 = -MAX_REAL + do i = 1, npts { + if ((pts[i,ztype] < zval) && (pts[i,ztype] > zval1)) { + refpt1 = i + zval1 = pts[refpt1,ztype] + } + } + + if (refpt1 != 0) { + refpt = refpt1 + newgraph = YES + } + } + + case 'n': + if (!IS_INDEFI (ztype)) { + refpt1 = 0 + zval = pts[refpt, ztype] + zval1 = MAX_REAL + do i = 1, npts { + if ((pts[i,ztype] > zval) && (pts[i,ztype] < zval1)) { + refpt1 = i + zval1 = pts[refpt1,ztype] + } + } + + if (refpt1 != 0) { + refpt = refpt1 + newgraph = YES + } + } + + case 'c': + # cursor read + i = igs_nearest (gp, ztype, refpt, axis, pts, npts, wx, wy, wcs) + call printf ("%g %g %g %g %g %g\n") + call pargr (pts[i, X]) + call pargr (pts[i, Y]) + call pargr (pts[i, Z]) + call pargr (pts[i, W]) + call pargr (pts[i, S]) + call pargr (pts[i, R]) + + case 'd': + i = igs_nearestd (gp, ztype, refpt, axis, pts, npts, wx, wy, + wcs) + if (i == 0) + goto 10 + + call gscur (gp, real (pts[i,axis[1]]), real (pts[i,axis[2]])) + + call printf ( "Delete 'p'oint or constant 'x', 'y', or 'z': ") + if (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF) + goto 10 + call printf ("\n") + + switch (key) { + case 'p': + dtype = 0 + case 'x': + dtype = X + case 'y': + dtype = Y + case 'z': + dtype = Z + default: + call printf ("\07\n") + goto 10 + } + + call igs_delete (gp, gt, ztype, i, axis, pts, npts, dtype) + + case 'u': + i = igs_nearestu (gp, ztype, refpt, axis, pts, npts, wx, wy, + wcs) + if (i == 0) + goto 10 + + call gscur (gp, real (pts[i,axis[1]]), real (pts[i,axis[2]])) + + call printf ( "Undelete 'p'oint or constant 'x', 'y', or 'z': ") + if (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF) + goto 10 + call printf ("\n") + + switch (key) { + case 'p': + dtype = 0 + case 'x': + dtype = X + case 'y': + dtype = Y + case 'z': + dtype = Z + default: + call printf ("\07\n") + goto 10 + } + + call igs_undelete (gp, gt, ztype, i, axis, pts, Memr[wts], + npts, dtype) + + case 'f': + #call printf ("Fitting ...") + #call flush (STDOUT) + call igs_solve (sf,pts[1,X],pts[1,Y],pts[1,Z],pts[1,W],npts) + call xgsvector (sf, pts[1,X], pts[1,Y], pts[1,S], npts) + call asubr (pts[1,Z], pts[1,S], pts[1,R], npts) + call aavgr (pts[1,R], npts, mean, rms) + call igs_params (gt) + newgraph = YES + + case 'w': + call gt_window (gt, gp, "cursor", newgraph) + + case 'I': + call fatal (0, "Interrupt") + + default: + # Ring the bell. + + call printf ("\07\n") + } + + if (newgraph == YES) { + call igs_graph (gp, gt, ztype, refpt, axis, pts, npts, labels) + newgraph = NO + } + } + + call mfree (wts, TY_REAL) + +30 call igs_graph (gplog, gt, ztype, refpt, axis, pts, npts, labels) + +end diff --git a/noao/twodspec/longslit/transform/igsfit/igsget.x b/noao/twodspec/longslit/transform/igsfit/igsget.x new file mode 100644 index 00000000..ccd1fb6c --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/igsget.x @@ -0,0 +1,62 @@ +include <pkg/igsfit.h> + +# IGS_GETI -- Get the value of an integer parameter. + +int procedure igs_geti (param) + +int param # IGS parameter + +include "igsfit.com" + +begin + switch (param) { + case IGS_XORDER: + return (xorder) + case IGS_YORDER: + return (yorder) + default: + call error (0, "igs_geti: Unknown parameter") + } +end + + +# IGS_GETS -- Get the value of a string parameter. + +procedure igs_gets (param, str, maxchar) + +int param # IGS parameter +char str[maxchar] # String +int maxchar # Maximum number of characters + +include "igsfit.com" + +begin + switch (param) { + case IGS_FUNCTION: + call strcpy (function, str, maxchar) + default: + call error (0, "igs_gets: Unknown parameter") + } +end + + +# IGS_GETR -- Get the values of real valued fitting parameters. + +real procedure igs_getr (param) + +int param # Parameter to be get + +include "igsfit.com" + +begin + switch (param) { + case IGS_XMIN: + return (xmin) + case IGS_XMAX: + return (xmax) + case IGS_YMIN: + return (ymin) + case IGS_YMAX: + return (ymax) + } +end diff --git a/noao/twodspec/longslit/transform/igsfit/igsgraph.x b/noao/twodspec/longslit/transform/igsfit/igsgraph.x new file mode 100644 index 00000000..83eba7e1 --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/igsgraph.x @@ -0,0 +1,73 @@ +include <mach.h> +include <gset.h> +include <pkg/gtools.h> +include <pkg/igsfit.h> + +procedure igs_graph (gp, gt, ztype, refpt, axis, pts, npts, labels) + +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +int ztype # Zoom type +int refpt # Reference point +int axis[2] # Axis definitions +real pts[npts, ARB] # Data +int npts # Number of pts points +char labels[SZ_LINE, ARB] # Data labels + +int i, x, y +real xmin, xmax, ymin, ymax, xsize, ysize, gt_getr() + +begin + if (gp == NULL) + return + + x = axis[1] + y = axis[2] + + call gt_sets (gt, GTXLABEL, labels[1, x]) + call gt_sets (gt, GTYLABEL, labels[1, y]) + xsize = gt_getr (gt, GTXSIZE) + ysize = gt_getr (gt, GTYSIZE) + + call gclear (gp) + + if (IS_INDEFI (ztype)) { + call gascale (gp, pts[1, x], npts, 1) + call gascale (gp, pts[1, y], npts, 2) + } else { + xmin = MAX_REAL + xmax = -MAX_REAL + ymin = MAX_REAL + ymax = -MAX_REAL + do i = 1, npts { + if (pts[i,ztype] != pts[refpt,ztype]) + next + xmin = min (xmin, pts[i,x]) + xmax = max (xmax, pts[i,x]) + ymin = min (ymin, pts[i,y]) + ymax = max (ymax, pts[i,y]) + } + call gswind (gp, xmin, xmax, ymin, ymax) + } + + call gt_swind (gp, gt) + call gt_labax (gp, gt) + + if (IS_INDEFI (ztype)) { + do i = 1, npts { + if (pts[i,W] == 0.) + call gmark (gp, pts[i,x], pts[i,y], GM_CROSS, xsize, ysize) + else + call gmark (gp, pts[i,x], pts[i,y], GM_PLUS, xsize, ysize) + } + } else { + do i = 1, npts { + if (pts[i,ztype] != pts[refpt,ztype]) + next + if (pts[i,W] == 0.) + call gmark (gp, pts[i,x], pts[i,y], GM_CROSS, xsize, ysize) + else + call gmark (gp, pts[i,x], pts[i,y], GM_PLUS, xsize, ysize) + } + } +end diff --git a/noao/twodspec/longslit/transform/igsfit/igsinit.x b/noao/twodspec/longslit/transform/igsfit/igsinit.x new file mode 100644 index 00000000..f084e7ff --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/igsinit.x @@ -0,0 +1,21 @@ +include <pkg/igsfit.h> + +# IGS_INIT -- Initialize the surface fitting parameters. + +procedure igs_init (function, xorder, yorder, xmin, xmax, ymin, ymax) + +char function[ARB] # Function +int xorder # X order +int yorder # Y order +real xmin, xmax # X range +real ymin, ymax # Y range + +begin + call igs_sets (IGS_FUNCTION, function) + call igs_seti (IGS_XORDER, xorder) + call igs_seti (IGS_YORDER, yorder) + call igs_setr (IGS_XMIN, xmin) + call igs_setr (IGS_XMAX, xmax) + call igs_setr (IGS_YMIN, ymin) + call igs_setr (IGS_YMAX, ymax) +end diff --git a/noao/twodspec/longslit/transform/igsfit/igsnearest.x b/noao/twodspec/longslit/transform/igsfit/igsnearest.x new file mode 100644 index 00000000..69888509 --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/igsnearest.x @@ -0,0 +1,51 @@ +include <mach.h> +include <gset.h> +include <pkg/igsfit.h> + +int procedure igs_nearest (gp, ztype, refpt, axis, pts, npts, wx, wy, wcs) + +pointer gp # GIO pointer +int ztype # Zoom type +int refpt # Reference point +int axis[2] # Axes +real pts[npts, ARB] # Data points +int npts # Number of data points +real wx, wy # Cursor coordinates +int wcs # WCS + +int i, j, x, y +real r2, r2min, x0, y0 + +begin + x = axis[1] + y = axis[2] + + call gctran (gp, wx, wy, wx, wy, wcs, 0) + r2min = MAX_REAL + j = 0 + + if (IS_INDEFI (ztype)) { + do i = 1, npts { + call gctran (gp, pts[i,x], pts[i,y], x0, y0, wcs, 0) + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + } else { + do i = 1, npts { + if (pts[i,ztype] != pts[refpt,ztype]) + next + call gctran (gp, pts[i,x], pts[i,y], x0, y0, wcs, 0) + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + } + + call gscur (gp, real (pts[j,x]), real (pts[j,y])) + return (j) +end diff --git a/noao/twodspec/longslit/transform/igsfit/igsparams.x b/noao/twodspec/longslit/transform/igsfit/igsparams.x new file mode 100644 index 00000000..9ecdd422 --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/igsparams.x @@ -0,0 +1,23 @@ +include <pkg/gtools.h> + +# IGS_PARAMS -- Set the GTOOLS parameter string. + +procedure igs_params (gt) + +pointer gt # GTOOLS pointer + +pointer params + +include "igsfit.com" + +begin + call malloc (params, SZ_LINE, TY_CHAR) + call sprintf (Memc[params], SZ_LINE, + "Function = %s, xorder = %d, yorder = %d, rms = %.4g") + call pargstr (function) + call pargi (xorder) + call pargi (yorder) + call pargr (rms) + call gt_sets (gt, GTPARAMS, Memc[params]) + call mfree (params, TY_CHAR) +end diff --git a/noao/twodspec/longslit/transform/igsfit/igsset.x b/noao/twodspec/longslit/transform/igsfit/igsset.x new file mode 100644 index 00000000..ea74e8c9 --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/igsset.x @@ -0,0 +1,59 @@ +include <pkg/igsfit.h> + +# IGS_SETS -- Set the values of string valued fitting parameters. + +procedure igs_sets (param, str) + +int param # Parameter to be set +char str[ARB] # String value + +include "igsfit.com" + +begin + switch (param) { + case IGS_FUNCTION: + call strcpy (str, function, SZ_LINE) + } +end + + +# IGS_SETI -- Set the values of integer valued fitting parameters. + +procedure igs_seti (param, ival) + +int param # Parameter to be set +int ival # Integer value + +include "igsfit.com" + +begin + switch (param) { + case IGS_XORDER: + xorder = ival + case IGS_YORDER: + yorder = ival + } +end + + +# IGS_SETR -- Set the values of real valued fitting parameters. + +procedure igs_setr (param, rval) + +int param # Parameter to be set +real rval # Real value + +include "igsfit.com" + +begin + switch (param) { + case IGS_XMIN: + xmin = rval + case IGS_XMAX: + xmax = rval + case IGS_YMIN: + ymin = rval + case IGS_YMAX: + ymax = rval + } +end diff --git a/noao/twodspec/longslit/transform/igsfit/igssolve.x b/noao/twodspec/longslit/transform/igsfit/igssolve.x new file mode 100644 index 00000000..a7e39354 --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/igssolve.x @@ -0,0 +1,173 @@ +include <math/gsurfit.h> + + +# IGS_SOLVE1 -- Fit z = f(x, y). + +define SFTYPES "|chebyshev|legendre|" # Surface types + +procedure igs_solve1 (sf, x, y, z, w, npts) + +pointer sf # GSURFIT pointer +real x[npts] # X points +real y[npts] # Y points +real z[npts] # Z points +real w[npts] # Weights +int npts # Number of points + +int i, nfunc, ix, iy +pointer sf1, sf2, resids + +int strdic() + +include "igsfit.com" + +begin + # Determine the function type. + + nfunc = strdic (function, function, SZ_LINE, SFTYPES) + + # Fit the first surface. + + ix = min (2, xorder) + iy = min (2, yorder) + call xgsinit (sf1, nfunc, ix, iy, NO, xmin, xmax, ymin, ymax) + call xgsfit (sf1, x, y, z, w, npts, WTS_USER, i) + + switch (i) { + case SINGULAR: + call eprintf ("Singular solution\n") + case NO_DEG_FREEDOM: + call error (0, "No degrees of freedom") + } + + # Evaluate the first surface and fit the residuals. + + call malloc (resids, npts, TY_REAL) + call xgsvector (sf1, x, y, Memr[resids], npts) + call asubr (z, Memr[resids], Memr[resids], npts) + + call xgsinit (sf2, nfunc, xorder, yorder, YES, xmin,xmax,ymin,ymax) + call xgsfit (sf2, x, y, Memr[resids], w, npts, WTS_USER, i) + + switch (i) { + case SINGULAR: + call eprintf ("Singular solution\n") + case NO_DEG_FREEDOM: + call error (0, "No degrees of freedom") + } + + # Add the two surfaces and free memory. + + call xgsadd (sf1, sf2, sf) + call xgsfree (sf1) + call xgsfree (sf2) + call mfree (resids, TY_REAL) +end + + +# IGS_SOLVE2 -- Fit z = x + f(y). + + +procedure igs_solve2 (sf, x, y, z, w, npts) + +pointer sf # GSURFIT pointer +real x[npts] # X points +real y[npts] # Y points +real z[npts] # Z points +real w[npts] # Weights +int npts # Number of points + +int i, nfunc +real a +pointer sf1 + +int strdic() +real xgsgcoeff() + +include "igsfit.com" + +begin + nfunc = strdic (function, function, SZ_LINE, SFTYPES) + call xgsinit (sf1, nfunc, 1, yorder, NO, xmin, xmax, ymin, ymax) + + call asubr (z, x, z, npts) + call xgsfit (sf1, x, y, z, w, npts, WTS_USER, i) + call aaddr (z, x, z, npts) + + switch (i) { + case SINGULAR: + call eprintf ("Singular solution\n") + case NO_DEG_FREEDOM: + call error (0, "No degrees of freedom") + } + + call xgsfree (sf) + call xgsinit (sf, nfunc, 2, yorder, NO, xmin, xmax, ymin, ymax) + a = xgsgcoeff (sf1, 1, 1) + + a = a + (xmin + xmax) / 2 + call xgsscoeff (sf, 1, 1, a) + + a = (xmax - xmin) / 2 + call xgsscoeff (sf, 2, 1, a) + + do i = 2, yorder { + a = xgsgcoeff (sf1, 1, i) + call xgsscoeff (sf, 1, i, a) + } + + call xgsfree (sf1) +end + +# IGS_SOLVE3 -- Fit z = y + f(x). + +procedure igs_solve3 (sf, x, y, z, w, npts) + +pointer sf # GSURFIT pointer +real x[npts] # X points +real y[npts] # Y points +real z[npts] # Z points +real w[npts] # Weights +int npts # Number of points + +int i, nfunc +real a +pointer sf1 + +int strdic() +real xgsgcoeff() + +include "igsfit.com" + +begin + nfunc = strdic (function, function, SZ_LINE, SFTYPES) + call xgsinit (sf1, nfunc, xorder, 1, NO, xmin, xmax, ymin, ymax) + + call asubr (z, y, z, npts) + call xgsfit (sf1, x, y, z, w, npts, WTS_USER, i) + call aaddr (z, y, z, npts) + + switch (i) { + case SINGULAR: + call eprintf ("Singular solution\n") + case NO_DEG_FREEDOM: + call error (0, "No degrees of freedom") + } + + call xgsfree (sf) + call xgsinit (sf, nfunc, xorder, 2, NO, xmin, xmax, ymin, ymax) + a = xgsgcoeff (sf1, 1, 1) + + a = a + (ymin + ymax) / 2 + call xgsscoeff (sf, 1, 1, a) + + a = (ymax - ymin) / 2 + call xgsscoeff (sf, 1, 2, a) + + do i = 2, xorder { + a = xgsgcoeff (sf1, i, 1) + call xgsscoeff (sf, i, 1, a) + } + + call xgsfree (sf1) +end diff --git a/noao/twodspec/longslit/transform/igsfit/igsundelete.x b/noao/twodspec/longslit/transform/igsfit/igsundelete.x new file mode 100644 index 00000000..dc7b802e --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/igsundelete.x @@ -0,0 +1,107 @@ +include <mach.h> +include <gset.h> +include <pkg/gtools.h> +include <pkg/igsfit.h> + +int procedure igs_nearestu (gp, ztype, refpt, axis, pts, npts, wx, wy, wcs) + +pointer gp # GIO pointer +int ztype # Zoom type +int refpt # Reference point +int axis[2] # Axes +real pts[npts, ARB] # Data points +int npts # Number of data points +real wx, wy # Cursor coordinates +int wcs # WCS + +int i, j, x, y +real r2, r2min, x0, y0 + +begin + x = axis[1] + y = axis[2] + + call gctran (gp, wx, wy, wx, wy, wcs, 0) + r2min = MAX_REAL + j = 0 + + if (IS_INDEFI (ztype)) { + do i = 1, npts { + if (pts[i,W] != 0.) + next + call gctran (gp, pts[i, x], pts[i, y], x0, y0, wcs, 0) + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + } else { + do i = 1, npts { + if ((pts[i,ztype] != pts[refpt,ztype]) || (pts[i,W] != 0.)) + next + call gctran (gp, pts[i, x], pts[i, y], x0, y0, wcs, 0) + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + } + + return (j) +end + + +# IGS_UNDELETE - Undelete point or subset. + +procedure igs_undelete (gp, gt, ztype, refpt, axis, pts, wts, npts, dtype) + +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +int ztype # Zoom type +int refpt # Reference point for undeletion +int axis[2] # Axes +real pts[npts, ARB] # Data points +real wts[npts] # Original weights +int npts # Number of data points +int dtype # Undeletion type + +int i, x, y +real xsize, ysize + +real gt_getr() + +begin + x = axis[1] + y = axis[2] + + xsize = gt_getr (gt, GTXSIZE) + ysize = gt_getr (gt, GTYSIZE) + + switch (dtype) { + case X, Y, Z: + do i = 1, npts { + if (!IS_INDEFI (ztype)) + if (pts[refpt,ztype] != pts[i,ztype]) + next + if (pts[refpt,dtype] != pts[i,dtype]) + next + call gseti (gp, G_PMLTYPE, 0) + call gmark (gp, pts[i,x], pts[i,y], GM_CROSS, xsize, ysize) + call gseti (gp, G_PMLTYPE, 1) + call gmark (gp, pts[i,x], pts[i,y], GM_PLUS, xsize, ysize) + if (wts[i] == 0) + wts[i] = 1 + pts[i,W] = wts[i] + } + default: + call gseti (gp, G_PMLTYPE, 0) + call gmark (gp, pts[refpt,x], pts[refpt,y], GM_CROSS, xsize, ysize) + call gseti (gp, G_PMLTYPE, 1) + call gmark (gp, pts[refpt,x], pts[refpt,y], GM_PLUS, xsize, ysize) + if (wts[refpt] == 0) + wts[refpt] = 1 + pts[refpt,W] = wts[refpt] + } +end diff --git a/noao/twodspec/longslit/transform/igsfit/mkpkg b/noao/twodspec/longslit/transform/igsfit/mkpkg new file mode 100644 index 00000000..ac5a6ca9 --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/mkpkg @@ -0,0 +1,21 @@ +# Interactive General Surface Fitting Package + +$checkout libpkg.a ../../ +$update libpkg.a +$checkin libpkg.a ../../ +$exit + +libpkg.a: + igscolon.x igsfit.com <gset.h> + igsdelete.x <gset.h> <mach.h> <pkg/gtools.h> <pkg/igsfit.h> + igsfit.x igsfit.com <mach.h> <pkg/gtools.h> <pkg/igsfit.h> + igsget.x igsfit.com <pkg/igsfit.h> + igsgraph.x <gset.h> <mach.h> <pkg/gtools.h> <pkg/igsfit.h> + igsinit.x <pkg/igsfit.h> + igsnearest.x <gset.h> <mach.h> <pkg/igsfit.h> + igsparams.x igsfit.com <pkg/gtools.h> + igsset.x igsfit.com <pkg/igsfit.h> + igssolve.x igsfit.com <math/gsurfit.h> + igsundelete.x <gset.h> <mach.h> <pkg/gtools.h> <pkg/igsfit.h> + xgs.x <math/gsurfit.h> + ; diff --git a/noao/twodspec/longslit/transform/igsfit/xgs.x b/noao/twodspec/longslit/transform/igsfit/xgs.x new file mode 100644 index 00000000..7d2ea331 --- /dev/null +++ b/noao/twodspec/longslit/transform/igsfit/xgs.x @@ -0,0 +1,243 @@ +include <math/gsurfit.h> + +# XGS -- These routines provide an interface between real input data and +# the double precision surface fitting. Rather than make the input data +# be double precision we only want the internal surface fitting arithmetic +# to be double. But the surface fitting package only provides real +# arithmetic for real input and double precision arithmetic for double +# precision input. Hence these interfaces. Note that the save and restore +# functions use double precision. + +# XGSINIT -- Procedure to initialize the surface descriptor. + +procedure xgsinit (sf, surface_type, xorder, yorder, xterms, xmin, xmax, + ymin, ymax) + +pointer sf # surface descriptor +int surface_type # type of surface to be fitted +int xorder # x order of surface to be fit +int yorder # y order of surface to be fit +int xterms # presence of cross terms +real xmin # minimum value of x +real xmax # maximum value of x +real ymin # minimum value of y +real ymax # maximum value of y + +begin + call dgsinit (sf, surface_type, xorder, yorder, xterms, double (xmin), + double (xmax), double (ymin), double (ymax)) +end + + +# XGSFIT -- Procedure to solve the normal equations for a surface. + +procedure xgsfit (sf, x, y, z, w, npts, wtflag, ier) + +pointer sf # surface descriptor +real x[npts] # array of x values +real y[npts] # array of y values +real z[npts] # data array +real w[npts] # array of weights +int npts # number of data points +int wtflag # type of weighting +int ier # ier = OK, everything OK + # ier = SINGULAR, matrix is singular, 1 or more + # coefficients are 0. + # ier = NO_DEG_FREEDOM, too few points to solve matrix + +pointer sp, xd, yd, zd, wd +errchk salloc + +begin + call smark (sp) + call salloc (xd, npts, TY_DOUBLE) + call salloc (yd, npts, TY_DOUBLE) + call salloc (zd, npts, TY_DOUBLE) + call salloc (wd, npts, TY_DOUBLE) + call achtrd (x, Memd[xd], npts) + call achtrd (y, Memd[yd], npts) + call achtrd (z, Memd[zd], npts) + call achtrd (w, Memd[wd], npts) + call dgsfit (sf, Memd[xd], Memd[yd], Memd[zd], Memd[wd], npts, + wtflag, ier) + call sfree (sp) +end + + +# XGSVECTOR -- Procedure to evaluate the fitted surface at an array of points. + +procedure xgsvector (sf, x, y, zfit, npts) + +pointer sf # pointer to surface descriptor structure +real x[ARB] # x value +real y[ARB] # y value +real zfit[ARB] # fits surface values +int npts # number of data points + +pointer sp, xd, yd, zd +errchk salloc + +begin + call smark (sp) + call salloc (xd, npts, TY_DOUBLE) + call salloc (yd, npts, TY_DOUBLE) + call salloc (zd, npts, TY_DOUBLE) + call achtrd (x, Memd[xd], npts) + call achtrd (y, Memd[yd], npts) + call dgsvector (sf, Memd[xd], Memd[yd], Memd[zd], npts) + call achtdr (Memd[zd], zfit, npts) + call sfree (sp) +end + + +# XGSEVAL -- Procedure to evaluate the fitted surface at a single point. + +real procedure xgseval (sf, x, y) + +pointer sf # pointer to surface descriptor structure +real x # x value +real y # y value + +double dgseval() + +begin + return (real (dgseval (sf, double (x), double (y)))) +end + + +# XGSADD -- Procedure to add the fits from two surfaces together. + +procedure xgsadd (sf1, sf2, sf3) + +pointer sf1 # pointer to the first surface +pointer sf2 # pointer to the second surface +pointer sf3 # pointer to the output surface + +begin + call dgsadd (sf1, sf2, sf3) +end + + +# XGSFREE -- Procedure to free the surface descriptor + +procedure xgsfree (sf) + +pointer sf # the surface descriptor + +begin + call dgsfree (sf) +end + + +# XGSGCOEFF -- Procedure to fetch a particular coefficient. + +real procedure xgsgcoeff (sf, xorder, yorder) + +pointer sf # pointer to the surface fitting descriptor +int xorder # X order of desired coefficent +int yorder # Y order of desired coefficent + +double dgsgcoeff() + +begin + return (real (dgsgcoeff (sf, xorder, yorder))) +end + + +# XGSSCOEFF -- Procedure to set a particular coefficient. + +procedure xgsscoeff (sf, xorder, yorder, coeff) + +pointer sf # pointer to the surface fitting descriptor +int xorder # X order of desired coefficent +int yorder # Y order of desired coefficent +real coeff # Coefficient value + +begin + call dgsscoeff (sf, xorder, yorder, double (coeff)) +end + + +# XGSGETR -- Procedure to fetch a real gsurfit parameter + +real procedure xgsgetr (sf, parameter) + +pointer sf # pointer to the surface fit +int parameter # parameter to be fetched + +double dgsgetd() + +begin + return (real (dgsgetd (sf, parameter))) +end + + +# XGSGETI -- Procedure to fetch an integer parameter + +int procedure xgsgeti (sf, parameter) + +pointer sf # pointer to the surface fit +int parameter # integer parameter + +int dgsgeti() + +begin + return (dgsgeti (sf, parameter)) +end + + +# XGSSAVE -- Procedure to save the surface fit for later use by the +# evaluate routines. +# +# NOTE THAT THIS USES DOUBLE PRECISION FOR THE COEFFICIENTS. + +procedure xgssave (sf, fit) + +pointer sf # pointer to the surface descriptor +double fit[ARB] # array for storing fit + +begin + call dgssave (sf, fit) +end + + +# XGSRESTORE -- Procedure to restore the surface fit stored by GSSAVE +# to the surface descriptor for use by the evaluating routines. +# +# NOTE THAT THIS USES DOUBLE PRECISION FOR THE COEFFICIENTS. + +procedure xgsrestore (sf, fit) + +pointer sf # surface descriptor +double fit[ARB] # array containing the surface parameters and + +begin + call dgsrestore (sf, fit) +end + + +# XGSDER -- Procedure to calculate a new surface which is a derivative of +# the previous surface + +procedure xgsder (sf1, x, y, zfit, npts, nxd, nyd) + +pointer sf1 # pointer to the previous surface +real x[npts] # x values +real y[npts] # y values +real zfit[npts] # fitted values +int npts # number of points +int nxd, nyd # order of the derivatives in x and y + +pointer sp, xd, yd, zd + +begin + call smark (sp) + call salloc (xd, npts, TY_DOUBLE) + call salloc (yd, npts, TY_DOUBLE) + call salloc (zd, npts, TY_DOUBLE) + call achtrd (x, Memd[xd], npts) + call achtrd (y, Memd[yd], npts) + call dgsder (sf1, Memd[xd], Memd[yd], Memd[zd], npts, nxd, nyd) + call achtdr (Memd[zd], zfit, npts) + call sfree (sp) +end diff --git a/noao/twodspec/longslit/transform/mkpkg b/noao/twodspec/longslit/transform/mkpkg new file mode 100644 index 00000000..8ea1b584 --- /dev/null +++ b/noao/twodspec/longslit/transform/mkpkg @@ -0,0 +1,20 @@ +# Coordinate Transformation Tasks + +$checkout libpkg.a ../ +$update libpkg.a +$checkin libpkg.a ../ +$exit + +libpkg.a: + @igsfit + + fcdbio.x <error.h> <math/gsurfit.h> <pkg/dttext.h> <units.h> + fcdlist.x <error.h> <mach.h> + fcfitcoords.x <pkg/gtools.h> <pkg/igsfit.h> <pkg/xtanswer.h> + fcgetcoords.x <imio.h> <mach.h> <pkg/dttext.h> <pkg/igsfit.h> + fcgetim.x + fitcoords.x <error.h> <pkg/igsfit.h> <pkg/xtanswer.h> + trsetup.x <math.h> <math/gsurfit.h> <math/iminterp.h> + t_fceval.x + t_transform.x transform.com <imhdr.h> <math/iminterp.h> <units.h> + ; diff --git a/noao/twodspec/longslit/transform/t_fceval.x b/noao/twodspec/longslit/transform/t_fceval.x new file mode 100644 index 00000000..a9c5cc75 --- /dev/null +++ b/noao/twodspec/longslit/transform/t_fceval.x @@ -0,0 +1,107 @@ +# T_FCEVAL -- Evaluate FITCOORDS solutions. +# Input consists of a text file of pixel coordinates to be evaluated and the +# user coordinate surfaces from FITCOORDS. The output is a text file of the +# input coordinates followed by the output coordinates. When there is no fit +# for an axis the unit transformation is used and when there is more than one +# fit for an axis the average is used. + +procedure t_fceval () + +pointer input # File of input coordinates +pointer output # File of output coordinates +int fitnames # List of user coordinate fits +pointer database # Database + +int i, j, in, out, nsf[2] +double x[2], y[2] +pointer sp, fitname, sf[2], un[2], sf1, un1 + +bool un_compare() +int open(), fscan(), nscan() +int clpopnu(), clplen(), clgfil() +double dgseval() +errchk open, lm_dbread + +begin + call smark (sp) + call salloc (input, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (database, SZ_FNAME, TY_CHAR) + call salloc (fitname, SZ_FNAME, TY_CHAR) + + # Get parameters. + call clgstr ("input", Memc[input], SZ_FNAME) + call clgstr ("output", Memc[output], SZ_FNAME) + fitnames = clpopnu ("fitnames") + call clgstr ("database", Memc[database], SZ_FNAME) + + # Open the input and output files. + in = open (Memc[input], READ_ONLY, TEXT_FILE) + out = open (Memc[output], NEW_FILE, TEXT_FILE) + + # Read the solutions. + i = max (1, clplen (fitnames)) + call salloc (sf[1], i, TY_INT) + call salloc (sf[2], i, TY_INT) + + nsf[1] = 0; nsf[2] = 0; un[1] = NULL; un[2] = NULL + while (clgfil (fitnames, Memc[fitname], SZ_FNAME) != EOF) { + call lm_dbread (Memc[database], Memc[fitname], j, un1, sf1) + if (un1 != NULL) { + if (un[j] == NULL) + un[j] = un1 + else if (un_compare (un1, un[j])) + call un_close (un1) + else + call error (1, "Input units disagree") + } + + if (sf1 != NULL) { + Memi[sf[j]+nsf[j]] = sf1 + nsf[j] = nsf[j] + 1 + } + } + + if (nsf[1] + nsf[2] == 0) + call error (0, "No user coordinates") + + # Evaluate the fits at each input coordinate. + while (fscan (in) != EOF) { + call gargd (x[1]) + call gargd (x[2]) + if (nscan() != 2) + next + + do j = 1, 2 { + if (nsf[j] == 0) + y[j] = x[j] + else { + y[j] = dgseval (Memi[sf[j]], x[1], x[2]) + do i = 2, nsf[1] + y[j] = y[j] + dgseval (Memi[sf[j]+i-1], x[1], y[2]) + y[j] = y[j] / nsf[j] + } + } + + call fprintf (out, "%g %g %g %g\n") + call pargd (x[1]) + call pargd (x[2]) + call pargd (y[1]) + call pargd (y[2]) + call flush (out) + } + + # Free the surfaces and units structures. + do j = 1, 2 { + for (i=1; i<=nsf[j]; i=i+1) + call dgsfree (Memi[sf[j]+i-1]) + if (un[j] != NULL) + call un_close (un[j]) + } + + # Finish up. + call clpcls (fitnames) + call close (out) + call close (in) + call sfree (sp) +end diff --git a/noao/twodspec/longslit/transform/t_transform.x b/noao/twodspec/longslit/transform/t_transform.x new file mode 100644 index 00000000..5610858e --- /dev/null +++ b/noao/twodspec/longslit/transform/t_transform.x @@ -0,0 +1,741 @@ +include <imhdr.h> +include <math/iminterp.h> +include <units.h> + +define ITYPES "|nearest|linear|poly3|poly5|spline3|" + +# T_TRANSFORM -- Transform longslit images. +# Input consists of images to be transformed, the user coordinate surfaces +# describing the output coordinates in terms of the input coordinates, +# and the desired coordinates for the output images. The type of image +# interpolation is also input. There is a log output as well as the +# transformed images. The output image may replace the input image. + +procedure t_transform () + +int input # List of input images +int output # List of output images +int minput # List of input masks +int moutput # List of output masks +int fitnames # List of user coordinate fits +pointer database # Database +char interp[10] # Interpolation type +int logfiles # List of log files + +int itypes[II_NTYPES2D], logfd, nusf, nvsf +pointer in, out, pmin, pmout +pointer un[2], mw, ct, usf, vsf, xmsi, ymsi, jmsi, xout, yout, dxout, dyout +pointer sp, image1, image2, image3, minname, moutname, mname, str + +int clpopnu(), clgfil(), clplen(), clgeti(), clgwrd(), open() +int imtopenp(), imtlen(), imtgetim() +bool clgetb() +real clgetr() +pointer immap(), mw_openim(), yt_mappm() +errchk tr_gsf, tr_setup, open, mw_openim, yt_mappm + +data itypes /II_BINEAREST, II_BILINEAR, II_BIPOLY3, II_BIPOLY5, + II_BISPLINE3, II_SINC, II_LSINC, II_DRIZZLE/ + +include "transform.com" + + +begin + call smark (sp) + call salloc (database, SZ_FNAME, TY_CHAR) + call salloc (image1, SZ_FNAME, TY_CHAR) + call salloc (image2, SZ_FNAME, TY_CHAR) + call salloc (image3, SZ_FNAME, TY_CHAR) + call salloc (minname, SZ_FNAME, TY_CHAR) + call salloc (moutname, SZ_FNAME, TY_CHAR) + call salloc (mname, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + + # Get and error check the input and output image lists and the other + # task parameters. + + input = imtopenp ("input") + output = imtopenp ("output") + if (imtlen (input) != imtlen (output)) { + call imtclose (input) + call imtclose (output) + call error (1, "Number of input and output images differ") + } + minput = imtopenp ("minput") + moutput = imtopenp ("moutput") + if (imtlen (minput) > 1 && imtlen (minput) != imtlen (input)) { + call imtclose (input) + call imtclose (output) + call imtclose (minput) + call imtclose (moutput) + call error (1, "Can't associate input masks with input images") + } + if (imtlen (moutput) > 0 && imtlen (input) != imtlen (moutput)) { + call imtclose (input) + call imtclose (output) + call imtclose (minput) + call imtclose (moutput) + call error (1, "Number output masks differ from input") + } + + fitnames = clpopnu ("fitnames") + call clgstr ("database", Memc[database], SZ_FNAME) + itype = itypes[clgwrd ("interptype", interp, 10, II_FUNCTIONS)] + logfiles = clpopnu ("logfiles") + + u1 = clgetr ("x1") + u2 = clgetr ("x2") + du = clgetr ("dx") + nu = clgeti ("nx") + v1 = clgetr ("y1") + v2 = clgetr ("y2") + dv = clgetr ("dy") + nv = clgeti ("ny") + + ulog = clgetb ("xlog") + vlog = clgetb ("ylog") + flux = clgetb ("flux") + blank = clgetr ("blank") + + usewcs = (clplen (fitnames) == 0) + + # Transform each input image to the output image. + Memc[minname] = EOS + Memc[moutname] = EOS + Memc[mname] = EOS + xmsi = NULL + while ((imtgetim (input, Memc[image1], SZ_FNAME) != EOF) && + (imtgetim (output, Memc[image2], SZ_FNAME) != EOF)) { + + # Get mask names. + if (imtgetim (minput, Memc[image3], SZ_FNAME) != EOF) + call strcpy (Memc[image3], Memc[minname], SZ_FNAME) + if (imtgetim (moutput, Memc[image3], SZ_FNAME) != EOF) + call strcpy (Memc[image3], Memc[moutname], SZ_FNAME) + + # Map the input and output images. + call xt_mkimtemp (Memc[image1], Memc[image2], Memc[image3],SZ_FNAME) + in = immap (Memc[image1], READ_ONLY, 0) + out = immap (Memc[image2], NEW_COPY, in) + + # Map masks. + pmin = NULL; pmout = NULL + if (Memc[minname] != EOS) + pmin = yt_mappm (Memc[minname], in, "logical", Memc[mname], + SZ_FNAME) + if (Memc[moutname] != EOS) { + call xt_maskname (Memc[moutname], "", NEW_IMAGE, + Memc[moutname], SZ_FNAME) + pmout = immap (Memc[moutname], NEW_COPY, in) + call imastr (out, "BPM", Memc[moutname]) + } + + # Get the coordinate transformation surfaces from the database + # and setup the transformations. + # Do this only on the first pass. + + if (xmsi == NULL) { + if (usewcs) { + mw = mw_openim (in) + call tr_gwcs (mw, un, IM_LEN(in,1), IM_LEN(in,2), ct, + usf, nusf, vsf, nvsf) + } else { + mw = NULL + ct = NULL + call tr_gsf (Memc[database], fitnames, un, usf, nusf, + vsf, nvsf) + } + call tr_setup (ct, usf, nusf, vsf, nvsf, un, xmsi, ymsi, jmsi, + xout, yout, dxout, dyout) + if (mw != NULL) + call mw_close (mw) + } + + # Write log information. + while (clgfil (logfiles, Memc[str], SZ_LINE) != EOF) { + logfd = open (Memc[str], APPEND, TEXT_FILE) + call sysid (Memc[str], SZ_LINE) + call fprintf (logfd, "\n%s\n") + call pargstr (Memc[str]) + call fprintf (logfd, " Transform %s to %s.\n") + call pargstr (Memc[image1]) + call pargstr (Memc[image3]) + if (pmout != EOS) { + if (pmin != EOS) { + call fprintf (logfd, " Transform mask %s to %s.\n") + call pargstr (Memc[mname]) + call pargstr (Memc[moutname]) + } else { + call fprintf (logfd, " Output mask is %s.\n") + call pargstr (Memc[moutname]) + } + } + if (flux) + call fprintf (logfd, " Conserve flux per pixel.\n") + if (usewcs) + call fprintf (logfd, " Transforming using image WCS.\n") + else { + call fprintf (logfd, " User coordinate transformations:\n") + while (clgfil (fitnames, Memc[str], SZ_LINE) != EOF) { + call fprintf (logfd, " %s\n") + call pargstr (Memc[str]) + } + } + call fprintf (logfd, " Interpolation is %s.\n") + call pargstr (interp) + if (!IS_INDEFR(blank)) { + call fprintf (logfd, " Out of bounds pixel value is %g.\n") + call pargr (blank) + } else + call fprintf (logfd, + " Using edge extension for out of bounds pixel values.\n") + call fprintf (logfd, " Output coordinate parameters are:\n") + call fprintf (logfd, + " x1 = %10.4g, x2 = %10.4g, dx = %10.4g, nx = %4d, xlog = %b\n") + call pargr (u1) + call pargr (u2) + call pargr (du) + call pargi (nu) + call pargb (ulog) + call fprintf (logfd, + " y1 = %10.4g, y2 = %10.4g, dy = %10.4g, ny = %4d, ylog = %b\n") + call pargr (v1) + call pargr (v2) + call pargr (dv) + call pargi (nv) + call pargb (vlog) + call close (logfd) + } + call clprew (logfiles) + + call tr_transform (in, out, pmin, pmout, un, xmsi, ymsi, jmsi, + Memr[xout], Memr[yout], Memr[dxout], Memr[dyout]) + + if (pmout != NULL) + call imunmap (pmout) + if (pmin != NULL) + call xt_pmunmap (pmin) + call imunmap (in) + call imunmap (out) + call xt_delimtemp (Memc[image2], Memc[image3]) + + if (usewcs) { + call mfree (xout, TY_REAL) + call mfree (yout, TY_REAL) + call mfree (dxout, TY_REAL) + call mfree (dyout, TY_REAL) + if (xmsi != NULL) + call msifree (xmsi) + if (ymsi != NULL) + call msifree (ymsi) + if (jmsi != NULL) + call msifree (jmsi) + if (un[1] != NULL) + call un_close (un[1]) + if (un[2] != NULL) + call un_close (un[2]) + xmsi = NULL + } + + } + + call mfree (xout, TY_REAL) + call mfree (yout, TY_REAL) + call mfree (dxout, TY_REAL) + call mfree (dyout, TY_REAL) + if (xmsi != NULL) + call msifree (xmsi) + if (ymsi != NULL) + call msifree (ymsi) + if (jmsi != NULL) + call msifree (jmsi) + if (un[1] != NULL) + call un_close (un[1]) + if (un[2] != NULL) + call un_close (un[2]) + call imtclose (minput) + call imtclose (moutput) + call imtclose (input) + call imtclose (output) + call clpcls (fitnames) + call clpcls (logfiles) + call sfree (sp) +end + + +# TR_SETOUTPUT -- Set the output coordinates in the common block. +# This procedure allows the user to specifying a part of the output +# coordinates and let the rest default based on the full limits of +# the user coordinate surfaces. + +procedure tr_setoutput (xmin, xmax, ymin, ymax, umin, umax, vmin, vmax) + +real xmin, xmax, ymin, ymax +real umin, umax, vmin, vmax + +int nua, nva +real u1a, u2a, dua, v1a, v2a, dva + +include "transform.com" + +begin + # Save the original values of the user parameters. + u1a = u1 + u2a = u2 + dua = du + nua = nu + v1a = v1 + v2a = v2 + dva = dv + nva = nv + + # If the output coordinate limits are not defined then use the + # transformation surface limits. + + if (IS_INDEF (u1)) + u1 = umin + if (IS_INDEF (u2)) + u2 = umax + if (IS_INDEF (v1)) + v1 = vmin + if (IS_INDEF (v2)) + v2 = vmax + + # If the number of output pixels are not defined then use the number + # of pixels in the input image. + + if (IS_INDEFI (nu)) + nu = xmax - xmin + 1 + if (IS_INDEFI (nv)) + nv = ymax - ymin + 1 + + # If the coordinate interval is not defined determine it from the + # number of pixels and the coordinate limits. If the interval is + # defined then override the number of pixels. + + if (ulog) { + if (IS_INDEF (du)) + du = (log10 (u2) - log10 (u1)) / (nu - 1) + else if (IS_INDEFI (nua)) + nu = nint ((log10 (u2) - log10 (u1)) / du + 1) + else if (IS_INDEF (u1a)) + u1 = 10.0 ** (log10 (u2) - du * (nu - 1)) + else + u2 = 10.0 ** (log10 (u1) + du * (nu - 1)) + } else { + if (IS_INDEF (du)) + du = (u2 - u1) / (nu - 1) + else if (IS_INDEFI (nua)) + nu = nint ((u2 - u1) / du + 1) + else if (IS_INDEF (u1a)) + u1 = u2 - du * (nu - 1) + else + u2 = u1 + du * (nu - 1) + } + + if (vlog) { + if (IS_INDEF (dv)) + dv = (log10 (v2) - log10 (v1)) / (nv - 1) + else if (IS_INDEFI (nva)) + nv = nint ((log10 (v2) - log10 (v1)) / dv + 1) + else if (IS_INDEF (v1a)) + v1 = 10.0 ** (log10 (v2) - dv * (nv - 1)) + else + v2 = 10.0 ** (log10 (v1) + dv * (nv - 1)) + } else { + if (IS_INDEF (dv)) + dv = (v2 - v1) / (nv - 1) + else if (IS_INDEFI (nva)) + nv = nint ((v2 - v1) / dv + 1) + else if (IS_INDEF (v1a)) + v1 = v2 - dv * (nv - 1) + else + v2 = v1 + dv * (nv - 1) + } +end + + +define NBUF 16 # Additional buffer for interpolation +define NEDGE 2 # Number of edge lines to add for interpolation +define MINTERP 100 # Mask value for input mask interpolation +define MTHRESH 10 # Interpolated mask value for bad pixels +define MBAD 1 # Mask value for output bad pixels +define MBLANK 1 # Mask value for out of bounds pixels + +# TR_TRANSFORM -- Perform the image transformation using a user specified +# image interpolator. If an input and output mask are included the input +# mask values are set to MINTERP, interpolated in the same way, and any values +# greater than MTHRESH are set to MBAD. Note that currently the input mask +# values are not used in computing the input data interpolation value. +# The masks MUST be the same size as the input data and are assumed to +# be registered in logical pixel coordinates. + +procedure tr_transform (in, out, pmin, pmout, un, xmsi, ymsi, jmsi, xout, yout, + dxout, dyout) + +pointer in, out #I IMIO data pointers +pointer pmin, pmout #I IMIO mask pointers (NULL if not used) +pointer un[2] #I Units +pointer xmsi, ymsi #I Coordinate interpolation pointers +pointer jmsi #I Jacobian interpolation pointer +real xout[ARB], yout[ARB] #I Output grid relative to interpolation surface +real dxout[ARB], dyout[ARB] #I Output coordinate intervals + +int i, j, nxin, nyin, line1, line2, line3, line4, nlines, laxis, paxis +bool xofb, yofb +real a, b, c, r[2], w[2], cd[2,2] +pointer zmsi, mzmsi, buf, mbuf, bufout +pointer sp, xin, yin, jbuf, xin1, yin1, y, mw + +pointer mw_open(), impl2r() +errchk get_daxis + +include "transform.com" + +begin + # Initialize the output image header. + + IM_LEN(out, 1) = nu + IM_LEN(out, 2) = nv + if (pmout != NULL) { + IM_LEN(pmout, 1) = nu + IM_LEN(pmout, 2) = nv + } + + mw = mw_open (NULL, 2) + call mw_newsystem (mw, "world", 2) + do i = 1, 2 { + call mw_swtype (mw, i, 1, "linear", "") + if (un[i] != NULL) { + call mw_swattrs (mw, i, "label", UN_LABEL(un[i])) + call mw_swattrs (mw, i, "units", UN_UNITS(un[i])) + } + } + + r[1] = 1. + if (ulog) + w[1] = log10 (u1) + else + w[1] = u1 + cd[1,1] = du + cd[1,2] = 0. + r[2] = 1. + if (vlog) + w[2] = log10 (v1) + else + w[2] = v1 + cd[2,2] = dv + cd[2,1] = 0. + call mw_swtermr (mw, r, w, cd, 2) + + # The following image parameters are for compatibility with the + # ONEDSPEC package if using database solutions. + + if (!usewcs) { + call imastr (out, "DCLOG1", "Transform") + iferr (call imdelf (out, "REFSPEC1")) + ; + iferr (call imdelf (out, "REFSPEC2")) + ; + call get_daxis (in, laxis, paxis) + call imaddi (out, "dispaxis", laxis) + switch (laxis) { + case 1: + if (ulog) + call imaddi (out, "dc-flag", 1) + else + call imaddi (out, "dc-flag", 0) + if (un[laxis] == NULL) { + call mw_swattrs (mw, laxis, "label", "Wavelength") + call mw_swattrs (mw, laxis, "units", "Angstroms") + } + case 2: + if (vlog) + call imaddi (out, "dc-flag", 1) + else + call imaddi (out, "dc-flag", 0) + if (un[laxis] == NULL) { + call mw_swattrs (mw, laxis, "label", "Wavelength") + call mw_swattrs (mw, laxis, "units", "Angstroms") + } + } + } + call mw_saveim (mw, out) + if (pmout != NULL) + call mw_saveim (mw, pmout) + call mw_close (mw) + + # Allocate memory for the input coordinates and a vector for the + # output y coordinates. Also initialize the image data buffer. + + call smark (sp) + call salloc (xin, nu, TY_REAL) + call salloc (yin, nu, TY_REAL) + call salloc (y, nu, TY_REAL) + if (flux) + call salloc (jbuf, nu, TY_REAL) + if (!IS_INDEFR(blank) || pmout != NULL) { + call salloc (xin1, nu, TY_REAL) + call salloc (yin1, nu, TY_REAL) + } + + buf = NULL + mbuf = NULL + nlines = 0 + + # Initialize the interpolator. + + call msiinit (zmsi, itype) + if (pmin != NULL) + call msiinit (mzmsi, itype) + + # Do each line of the output image. + + nxin = IM_LEN(in, 1) + nyin = IM_LEN(in, 2) + + do i = 1, nv { + + # Evaluate the input coordinates at the output grid for a line + # of the output image using the interpolation surfaces. + + call amovkr (yout[i], Memr[y], nu) + if (!IS_INDEFR(blank) || pmout != NULL) { + call msivector (xmsi, xout, Memr[y], Memr[xin1], nu) + call msivector (ymsi, xout, Memr[y], Memr[yin1], nu) + call amovr (Memr[xin1], Memr[xin], nu) + call amovr (Memr[yin1], Memr[yin], nu) + } else { + call msivector (xmsi, xout, Memr[y], Memr[xin], nu) + call msivector (ymsi, xout, Memr[y], Memr[yin], nu) + } + + # Determine the coordinate ranges and check for out of bounds. + + call alimr (Memr[xin], nu, a, b) + xofb = (a < 1 || b > nxin) + if (xofb) { + if (a < 1) + call arltr (Memr[xin], nu, 1., 1.) + if (b > nxin) + call argtr (Memr[xin], nu, real (nxin), real (nxin)) + } + + call alimr (Memr[yin], nu, a, b) + yofb = (a < 1 || b > nyin) + if (yofb) { + if (a < 1) { + call arltr (Memr[yin], nu, 1., 1.) + a = 1. + b = max (a, b) + } + if (b > nyin) { + call argtr (Memr[yin], nu, real (nyin), real (nyin)) + b = nyin + a = min (a, b) + } + } + + # Get the input image data and fit an interpolator to the data. + + if ((buf == NULL) || (b > line2) || (a < line1)) { + nlines = max (nlines, int (b - a + 2 + NBUF)) + if (buf == NULL) { + if (a < nyin / 2) { + line1 = max (1, int (a)) + line2 = min (nyin, line1 + nlines - 1) + } else { + line2 = min (nyin, int (b+1.)) + line1 = max (1, line2 - nlines + 1) + } + } else if (b > line2) { + line1 = max (1, int (a)) + line2 = min (nyin, line1 + nlines - 1) + line1 = max (1, line2 - nlines + 1) + } else { + line2 = min (nyin, int (b+1.)) + line1 = max (1, line2 - nlines + 1) + line2 = min (nyin, line1 + nlines - 1) + } + line3 = max (1, line1 - NEDGE) + line4 = min (nyin, line2 + NEDGE) + call tr_bufl2r (in, pmin, line3, line4, buf, mbuf) + call msifit (zmsi, Memr[buf], nxin, line4 - line3 + 1, nxin) + if (pmin != NULL) + call msifit (mzmsi, Memr[mbuf], nxin, line4 - line3 + 1, + nxin) + } + + # The input coordinates must be offset to interpolation data grid. + call asubkr (Memr[yin], real (line3 - 1), Memr[yin], nu) + + # Evaluate output image pixels, conserve flux (if requested) using + # the Jacobian, and set the out of bounds values. + + bufout = impl2r (out, i) + call msivector (zmsi, Memr[xin], Memr[yin], Memr[bufout], nu) + if (flux) { + call msivector (jmsi, xout, Memr[y], Memr[jbuf], nu) + call amulr (dxout, Memr[jbuf], Memr[jbuf], nu) + call amulkr (Memr[jbuf], dyout[i], Memr[jbuf], nu) + call amulr (Memr[bufout], Memr[jbuf], Memr[bufout], nu) + } + if (!IS_INDEFR(blank)) { + if (xofb) { + do j = 0, nu-1 { + if (Memr[xin1+j] < 1 || Memr[xin1+j] > nxin) + Memr[bufout+j] = blank + } + } + if (yofb) { + do j = 0, nu-1 { + if (Memr[yin1+j] < 1 || Memr[yin1+j] > nyin) + Memr[bufout+j] = blank + } + } + } + + # Evaluate output mask pixels and set output bad values. + + if (pmout != NULL) { + bufout = impl2r (pmout, i) + if (pmin != NULL) { + call msivector (mzmsi, Memr[xin], Memr[yin], Memr[bufout], + nu) + do j = 0, nu-1 { + c = Memr[bufout+j] + if (Memr[xin1+j] < 1 || Memr[xin1+j] > nxin || + Memr[yin1+j] < 1 || Memr[yin1+j] > nyin) + Memr[bufout+j] = MBLANK + else if (c > 0.) { + if (c > MTHRESH) + Memr[bufout+j] = MBAD + else + Memr[bufout+j] = 0 + } + } + } else { + call aclrr (Memr[bufout], nu) + if (xofb) { + do j = 0, nu-1 { + if (Memr[xin1+j] < 1 || Memr[xin1+j] > nxin) + Memr[bufout+j] = MBLANK + } + } + if (yofb) { + do j = 0, nu-1 { + if (Memr[yin1+j] < 1 || Memr[yin1+j] > nyin) + Memr[bufout+j] = MBLANK + } + } + } + } + } + + # Free memory. + + call mfree (buf, TY_REAL) + call mfree (mbuf, TY_REAL) + call msifree (zmsi) + if (pmin != NULL) + call msifree (mzmsi) + call sfree (sp) +end + + +# TR_BUFL2R -- Maintain buffer of image lines. A new buffer is created when +# the buffer pointer is null or if the number of lines requested is changed. +# The minimum number of image reads is used. + +procedure tr_bufl2r (im, pmin, line1, line2, buf, mbuf) + +pointer im #I Image pointer +pointer pmin #I Mask pointer +int line1 #I First image line of buffer +int line2 #I Last image line of buffer +pointer buf #U Output data buffer +pointer mbuf #U Output mask buffer + +int i, nlines, nx, last1, last2, nlast +pointer buf1, buf2 + +pointer imgl2r() + +begin + nlines = line2 - line1 + 1 + + # If the buffer pointer is undefined then allocate memory for the + # buffer. If the number of lines requested changes reallocate + # the buffer. Initialize the last line values to force a full + # buffer image read. + + if (buf == NULL) { + nx = IM_LEN(im, 1) + call malloc (buf, nx * nlines, TY_REAL) + if (pmin != NULL) + call malloc (mbuf, nx * nlines, TY_REAL) + last1 = line1 - nlines + last2 = line2 - nlines + } else if (nlines != nlast) { + call realloc (buf, nx * nlines, TY_REAL) + if (pmin != NULL) + call realloc (mbuf, nx * nlines, TY_REAL) + last1 = line1 - nlines + last2 = line2 - nlines + } + + # Read only the image lines with are different from the last buffer. + + if (line1 < last1) { + do i = line2, line1, -1 { + if (i > last1) + buf1 = buf + (i - last1) * nx + else + buf1 = imgl2r (im, i) + + buf2 = buf + (i - line1) * nx + call amovr (Memr[buf1], Memr[buf2], nx) + } + } else if (line2 > last2) { + do i = line1, line2 { + if (i < last2) + buf1 = buf + (i - last1) * nx + else + buf1 = imgl2r (im, i) + + buf2 = buf + (i - line1) * nx + call amovr (Memr[buf1], Memr[buf2], nx) + } + } + if (pmin != NULL) { + if (line1 < last1) { + do i = line2, line1, -1 { + if (i > last1) + buf1 = mbuf + (i - last1) * nx + else + buf1 = imgl2r (pmin, i) + + buf2 = mbuf + (i - line1) * nx + call amovr (Memr[buf1], Memr[buf2], nx) + call argtr (Memr[buf2], nx, 0.1, real(MINTERP)) + } + } else if (line2 > last2) { + do i = line1, line2 { + if (i < last2) + buf1 = mbuf + (i - last1) * nx + else + buf1 = imgl2r (pmin, i) + + buf2 = mbuf + (i - line1) * nx + call amovr (Memr[buf1], Memr[buf2], nx) + call argtr (Memr[buf2], nx, 0.1, real(MINTERP)) + } + } + } + + # Save the buffer parameters. + + last1 = line1 + last2 = line2 + nlast = nlines +end diff --git a/noao/twodspec/longslit/transform/transform.com b/noao/twodspec/longslit/transform/transform.com new file mode 100644 index 00000000..baaae3ab --- /dev/null +++ b/noao/twodspec/longslit/transform/transform.com @@ -0,0 +1,14 @@ +# TRANSFORM -- Common task parameters. + +int itype # Interpolation type +real u1, v1 # Starting coordinates +real u2, v2 # Ending coordinates +real du, dv # Coordinate intervals +int nu, nv # Number of pixels +bool ulog, vlog # Logrithmic coordinates? +bool flux # Conserve flux per pixel? +bool usewcs # Use WCS? +real blank # Blank value + +common /trcom/ u1, v1, u2, v2, du, dv, nu, nv, itype, ulog, vlog, + flux, usewcs, blank diff --git a/noao/twodspec/longslit/transform/trsetup.x b/noao/twodspec/longslit/transform/trsetup.x new file mode 100644 index 00000000..72db570d --- /dev/null +++ b/noao/twodspec/longslit/transform/trsetup.x @@ -0,0 +1,663 @@ +include <math.h> +include <math/gsurfit.h> +include <math/iminterp.h> + +# Wrapper for MWCS CT pointer to include the image pixel range. + +define CT_LW Memi[$1] # MWCS CT (logical -> world) +define CT_WL Memi[$1+1] # MWCS CT (world -> logical) +define CT_NX Memi[$1+2] # Number of pixels in X +define CT_NY Memi[$1+3] # Number of pixels Y + + +# TR_GSF -- Get coordinate surface fits from the database. + +procedure tr_gsf (database, sflist, un, usf, nusf, vsf, nvsf) + +char database #I Database containing coordinate surfaces +int sflist #I List of user coordinate surfaces +pointer un[2] #O Units pointers +pointer usf #O Pointer to array of U surface fits +int nusf #O Number of U surface fits +pointer vsf #O Pointer to array of V surface fits +int nvsf #O Number of U surface fits + +int i, nsf +pointer sp, sfname, un1, sf + +bool un_compare() +int clgfil(), clplen() + +begin + # Get the user coordinate surfaces and separate them into U and V. + # Check that all surfaces have the same range of X and Y and determine + # the range of U and V. + + call smark (sp) + call salloc (sfname, SZ_FNAME, TY_CHAR) + + nsf = max (1, clplen (sflist)) + call malloc (usf, nsf, TY_INT) + call malloc (vsf, nsf, TY_INT) + + un[1] = NULL + un[2] = NULL + Memi[usf] = NULL + Memi[vsf] = NULL + nusf = 0 + nvsf = 0 + while (clgfil (sflist, Memc[sfname], SZ_FNAME) != EOF) { + call lm_dbread (database, Memc[sfname], i, un1, sf) + if (un1 != NULL) { + if (un[i] == NULL) + un[i] = un1 + else if (un_compare (un1, un[i])) + call un_close (un1) + else { + call un_close (un1) + call un_close (un[i]) + call sfree (sp) + call error (1, "Input units disagree") + } + } + + if (sf != NULL) { + if (i == 1) { + nusf = nusf+1 + Memi[usf+nusf-1] = sf + } else if (i == 2) { + nvsf = nvsf+1 + Memi[vsf+nvsf-1] = sf + } + } + } + call clprew (sflist) + + if (nusf + nvsf == 0) + call error (0, "No user coordinates") + + call sfree (sp) +end + + +# TR_GWCS -- Get WCS. + +procedure tr_gwcs (mw, un, nx, ny, ct, usf, nusf, vsf, nvsf) + +pointer mw #I MWCS pointer +pointer un[2] #O Units pointers +int nx, ny #I Image size + +pointer ct #O CT pointer +pointer usf #O Pointer to array of U surface fits +int nusf #O Number of U surface fits +pointer vsf #O Pointer to array of V surface fits +int nvsf #O Number of U surface fits + +int i +pointer sp, units, un_open(), mw_sctran() +errchk un_open + +begin + call smark (sp) + call salloc (units, SZ_FNAME, TY_CHAR) + + call malloc (ct, 4, TY_STRUCT) + nusf = 1 + call calloc (usf, nusf, TY_INT) + nvsf = 1 + call calloc (vsf, nvsf, TY_INT) + + CT_LW(ct) = mw_sctran (mw, "logical", "world", 3) + CT_WL(ct) = mw_sctran (mw, "world", "logical", 3) + CT_NX(ct) = nx + CT_NY(ct) = ny + + do i = 1, 2 { + ifnoerr (call mw_gwattrs (mw, i, "units", Memc[units], SZ_FNAME)) + un[i] = un_open (Memc[units]) + else + un[i] = NULL + } +end + + +# TR_SETUP -- Setup the transformation interpolation. +# +# At each point (U,V) in the output image we need to know the coordinate +# (X,Y) of the input images to be interpolated. This means we need +# to determine X(U,V) and Y(U,V). The input user coordinate surfaces, +# however, are U(X,Y) and V(X,Y) (a missing surface implies a one to one +# mapping of U=X or V=Y). This requires simultaneously inverting the user +# coordinate surfaces. This is a slow process using a gradient following +# iterative technique. +# +# Note that when an WCS is used, the MWCS routines already provide the +# inverse mapping. But even in this case it may be slow and so we use the +# same sampling and surface fitting technique for setting up the inversion +# mapping. +# +# The inverted coordinates are determined on a evenly subsampled grid of +# linear output coordinates. A linear interpolation surface can then be fit +# to this grid which is much faster to evaluate at each output coordinate. +# These interpolation surfaces are returned. If flux is to be conserved a +# similar interpolation surface for the Jacobian, J(U,V) is also returned. +# There may also be a mapping of the output image into logrithmic intervals +# which maps to the linearly sampled interpolation surfaces. The mappings +# of the output U and V intervals to the subsampled interpolation coordinates +# are also returned. +# +# 1. Set the output coordinate system based on the ranges of X, Y, U, and V. +# 2. Determine X(U,V), Y(U,V), and J(U,V) on a evenly subsampled grid of +# U and V. +# 3. Fit linear interpolation surfaces to these data. +# 4. Compute the mapping between output coordinates along each axis, which +# may be logrithmic, into the subsampling interpolation coordinates. + +procedure tr_setup (ct, usf, nusf, vsf, nvsf, un, xmsi, ymsi, jmsi, + uout, vout, duout, dvout) + +pointer ct #I CT pointer +pointer usf #U Pointers to U surface fits: freed upon return +int nusf #I Number of U surface fits +pointer vsf #U Pointers to V surface fits: freed upon return +int nvsf #I Number of V surface fits +pointer un[2] #O Units pointers +pointer xmsi, ymsi, jmsi #O Surface interpolators for X, Y and Jacobian +pointer uout, vout #O Output coordinates relative to interpolator +pointer duout, dvout #O Output coordinate intervals + +int i, j, step, nu1, nv1 +real xmin, xmax, ymin, ymax, umin, umax, vmin, vmax +real u, v, x, y, du1, dv1, der[8] +double dval +pointer xgrid, ygrid, zgrid, ptr1, ptr2, ptr3 + +real tr_getr(), tr_eval() + +include "transform.com" + +begin + #step = clgeti ("step") + step = 10 + + xmin = INDEF + xmax = INDEF + ymin = INDEF + ymax = INDEF + umin = INDEF + umax = INDEF + vmin = INDEF + vmax = INDEF + do i = 1, nusf { + if (IS_INDEF (xmin)) { + xmin = tr_getr (ct, Memi[usf+i-1], GSXMIN) + xmax = tr_getr (ct, Memi[usf+i-1], GSXMAX) + ymin = tr_getr (ct, Memi[usf+i-1], GSYMIN) + ymax = tr_getr (ct, Memi[usf+i-1], GSYMAX) + } else { + if ((xmin != tr_getr (ct, Memi[usf+i-1], GSXMIN)) || + (xmax != tr_getr (ct, Memi[usf+i-1], GSXMAX)) || + (ymin != tr_getr (ct, Memi[usf+i-1], GSYMIN)) || + (ymax != tr_getr (ct, Memi[usf+i-1], GSYMAX))) + call error (0, "tr_setup: Inconsistent coordinate fits") + } + + if (IS_INDEF (umin)) { + umin = tr_eval (ct, Memi[usf+i-1], 1, xmin, ymin) + umax = umin + } + u = tr_eval (ct, Memi[usf+i-1], 1, xmin, ymin) + umin = min (u, umin) + umax = max (u, umax) + u = tr_eval (ct, Memi[usf+i-1], 1, xmax, ymin) + umin = min (u, umin) + umax = max (u, umax) + u = tr_eval (ct, Memi[usf+i-1], 1, xmin, ymax) + umin = min (u, umin) + umax = max (u, umax) + u = tr_eval (ct, Memi[usf+i-1], 1, xmax, ymax) + umin = min (u, umin) + umax = max (u, umax) + } + do i = 1, nvsf { + if (IS_INDEF (xmin)) { + xmin = tr_getr (ct, Memi[vsf+i-1], GSXMIN) + xmax = tr_getr (ct, Memi[vsf+i-1], GSXMAX) + ymin = tr_getr (ct, Memi[vsf+i-1], GSYMIN) + ymax = tr_getr (ct, Memi[vsf+i-1], GSYMAX) + } else { + if ((xmin != tr_getr (ct, Memi[vsf+i-1], GSXMIN)) || + (xmax != tr_getr (ct, Memi[vsf+i-1], GSXMAX)) || + (ymin != tr_getr (ct, Memi[vsf+i-1], GSYMIN)) || + (ymax != tr_getr (ct, Memi[vsf+i-1], GSYMAX))) + call error (0, "tr_setup: Inconsistent coordinate fits") + } + + if (IS_INDEF (vmin)) { + vmin = tr_eval (ct, Memi[vsf+i-1], 2, xmin, ymin) + vmax = vmin + } + v = tr_eval (ct, Memi[vsf+i-1], 2, xmin, ymin) + vmin = min (v, vmin) + vmax = max (v, vmax) + v = tr_eval (ct, Memi[vsf+i-1], 2, xmax, ymin) + vmin = min (v, vmin) + vmax = max (v, vmax) + v = tr_eval (ct, Memi[vsf+i-1], 2, xmin, ymax) + vmin = min (v, vmin) + vmax = max (v, vmax) + v = tr_eval (ct, Memi[vsf+i-1], 2, xmax, ymax) + vmin = min (v, vmin) + vmax = max (v, vmax) + } + if (IS_INDEF (umin)) { + umin = xmin + umax = xmax + } + if (IS_INDEF (vmin)) { + vmin = ymin + vmax = ymax + } + + # Set the output coordinate system which is in a common block. + call tr_setoutput (xmin, xmax, ymin, ymax, umin, umax, vmin, vmax) + + # Subsample the inverted coordinates and fit an interpolation + # surface. The grid is evaluated in a back and forth pattern to + # use the last point evaluated and the starting point for the next + # point. This allows the interative inversion routine to work most + # efficiently with typically only two evaluations per step. + + nu1 = max (2, nu / step) + nv1 = max (2, nv / step) + du1 = (u2 - u1) / (nu1 - 1) + dv1 = (v2 - v1) / (nv1 - 1) + + call malloc (xgrid, nu1 * nv1, TY_REAL) + call malloc (ygrid, nu1 * nv1, TY_REAL) + call malloc (zgrid, nu1 * nv1, TY_REAL) + + call tr_init (ct, Memi[usf], nusf, Memi[vsf], nvsf, xmin, ymin, der) + do i = 1, nv1, 2 { + # Do this line from left to right. + ptr1 = xgrid + (i - 1) * nu1 - 1 + ptr2 = ygrid + (i - 1) * nu1 - 1 + ptr3 = zgrid + (i - 1) * nu1 - 1 + v = v1 + (i - 1) * dv1 + do j = 1, nu1 { + u = u1 + (j - 1) * du1 + call tr_invert (ct, Memi[usf], nusf, Memi[vsf], nvsf, u, v, + x, y, der, xmin, xmax, ymin, ymax) + # V2.10.2 + #Memr[ptr1+j] = der[1] + #Memr[ptr2+j] = der[2] + # After V2.10.3 + Memr[ptr1+j] = x + Memr[ptr2+j] = y + + Memr[ptr3+j] = 1. / abs (der[4] * der[8] - der[5] * der[7]) + } + if (i == nv1) + break + + # Do the next line from right to left. + ptr1 = xgrid + i * nu1 - 1 + ptr2 = ygrid + i * nu1 - 1 + ptr3 = zgrid + i * nu1 - 1 + v = v1 + i * dv1 + do j = nu1, 1, -1 { + u = u1 + (j - 1) * du1 + call tr_invert (ct, Memi[usf], nusf, Memi[vsf], nvsf, u, v, + x, y, der, xmin, xmax, ymin, ymax) + # V2.10.2 + #Memr[ptr1+j] = der[1] + #Memr[ptr2+j] = der[2] + # V2.10.3 + Memr[ptr1+j] = x + Memr[ptr2+j] = y + Memr[ptr3+j] = 1. / abs (der[4] * der[8] - der[5] * der[7]) + } + } + + # Free the surfaces since we are now done with them. + if (ct != NULL) + call mfree (ct, TY_STRUCT) + for (i=1; i<=nusf; i=i+1) + if (Memi[usf+i-1] != NULL) + call xgsfree (Memi[usf+i-1]) + call mfree (usf, TY_POINTER) + for (i=1; i<=nvsf; i=i+1) + if (Memi[vsf+i-1] != NULL) + call xgsfree (Memi[vsf+i-1]) + call mfree (vsf, TY_POINTER) + + # Fit a linear interpolator to the subsampled grids of X(U,V), Y(U,V), + # and J(U,V) to avoid having to evaluate the inverse at each point in + # the output image. The inversion is slow because of the many + # evaluations of the surfaces coordinates. Also compute an return + # arrays mapping the output coordinates to the subsampled coordinates. + # This may include a transformation to logrithmic intervals. + + call msiinit (xmsi, II_BILINEAR) + call msifit (xmsi, Memr[xgrid], nu1, nv1, nu1) + call mfree (xgrid, TY_REAL) + + call msiinit (ymsi, II_BILINEAR) + call msifit (ymsi, Memr[ygrid], nu1, nv1, nu1) + call mfree (ygrid, TY_REAL) + + if (flux) { + call msiinit (jmsi, II_BILINEAR) + call msifit (jmsi, Memr[zgrid], nu1, nv1, nu1) + } + call mfree (zgrid, TY_REAL) + + # Compute the mapping between output coordinates and the subsampled + # interpolation surface. Also compute the intervals used to define + # the pixel areas for conserving flux. + + call malloc (uout, nu, TY_REAL) + call malloc (duout, nu, TY_REAL) + if (ulog) { + dval = log10 (double(u1)) + do i = 0, nu - 1 + Memr[uout+i] = 10.**(dval+i*du) + call amulkr (Memr[uout], du * LN_10, Memr[duout], nu) + } else { + do i = 0, nu - 1 + Memr[uout+i] = u1 + i * du + call amovkr (du, Memr[duout], nu) + } + u2 = Memr[uout+nu-1] + + call malloc (vout, nv, TY_REAL) + call malloc (dvout, nv, TY_REAL) + if (vlog) { + dval = log10 (double(v1)) + do i = 0, nv - 1 + Memr[vout+i] = 10.**(dval+i*dv) + call amulkr (Memr[vout], dv * LN_10, Memr[dvout], nv) + } else { + do i = 0, nv - 1 + Memr[vout+i] = v1 + i * dv + call amovkr (dv, Memr[dvout], nv) + } + v2 = Memr[vout+nv-1] + + # Convert to interpolation coordinates. + umin = 1.; umax = nu + do i = 0, nu - 1 + Memr[uout+i] = max (umin, min (umax, (Memr[uout+i]-u1)/du1+1)) + vmin = 1.; vmax = nv + do i = 0, nv - 1 + Memr[vout+i] = max (vmin, min (vmax, (Memr[vout+i]-v1)/dv1+1)) +end + + +define MAX_ITERATE 10 +define ERROR 0.05 +define FUDGE 0.5 + +# TR_INVERT -- Given user coordinate surfaces U(X,Y) and V(X,Y) +# (if none use one-to-one mapping and if more than one average) +# corresponding to a given U and V and also the various partial +# derivatives. This is done using a gradient following interative +# method based on evaluating the partial derivative at each point +# and solving the linear Taylor expansions simultaneously. The last +# point sampled is used as the starting point. Thus, if the +# input U and V progress smoothly then the number of iterations +# can be small. The output is returned in x and y and in the derivative array +# DER. A point outside of the surfaces is returned as the nearest +# point at the edge of the surfaces in the DER array. +# +# If a WCS is used then we let MWCS do the inversion and compute the +# derivatives numerically. + +procedure tr_invert (ct, usf, nusf, vsf, nvsf, u, v, x, y, der, + xmin, xmax, ymin, ymax) + +pointer ct #I CT pointer +pointer usf[ARB], vsf[ARB] #I User coordinate surfaces U(X,Y) and V(X,Y) +int nusf, nvsf #I Number of surfaces for each coordinate +real u, v #I Input U and V to determine X and Y +real x, y #O Output X and Y +real der[8] #U Last result as input, new result as output + # 1=X, 2=Y, 3=U, 4=DUDX, 5=DUDY, 6=V, + # 7=DVDX, 8=DVDY +real xmin, xmax, ymin, ymax #I Limits of coordinate surfaces. + +int i, j, nedge +real fudge, du, dv, dx, dy, a, b, tmp[4] + +begin + # If using a WCS we let MWCS do the inversion. + if (ct != NULL) { + call mw_c2tranr (CT_WL(ct), u, v, x, y) + call mw_c2tranr (CT_LW(ct), x-0.5, y, tmp[1], tmp[3]) + call mw_c2tranr (CT_LW(ct), x+0.5, y, tmp[2], tmp[4]) + der[4] = tmp[2] - tmp[1] + der[7] = tmp[4] - tmp[3] + call mw_c2tranr (CT_LW(ct), x, y-0.5, tmp[1], tmp[3]) + call mw_c2tranr (CT_LW(ct), x, y+0.5, tmp[2], tmp[4]) + der[5] = tmp[2] - tmp[1] + der[8] = tmp[4] - tmp[3] + return + } + + # Use the last result as the starting point for the next position. + # If this is near the desired value then the interation will converge + # quickly. Allow a iteration to go off the surface twice. + # Quit when DX and DY are within ERROR. + + nedge = 0 + do i = 1, MAX_ITERATE { + du = u - der[3] + dv = v - der[6] + a = der[8] * du - der[5] * dv + b = der[8] * der[4] - der[5] * der[7] + if (b == 0.) { + if (a < 0.) + dx = -2. + else + dx = 2. + } else + dx = a / b + a = dv - der[7] * dx + b = der[8] + if (b == 0.) { + if (a < 0.) + dy = -2. + else + dy = 2. + } else + dy = a / b + fudge = 1 - FUDGE / i + x = der[1] + fudge * dx + y = der[2] + fudge * dy + der[1] = max (xmin, min (xmax, x)) + der[2] = max (ymin, min (ymax, y)) +# if (x < xmin || x > xmax) +# nedge = nedge + 1 +# if (y < ymin || y > ymax) +# nedge = nedge + 1 +# if (nedge > 2) +# break + if ((abs (dx) < ERROR) && (abs (dy) < ERROR)) + break + + if (nusf == 0) + der[3] = der[1] + else if (nusf == 1) { + call xgsder (usf[1], der[1], der[2], der[3], 1, 0, 0) + call xgsder (usf[1], der[1], der[2], der[4], 1, 1, 0) + call xgsder (usf[1], der[1], der[2], der[5], 1, 0, 1) + } else { + call xgsder (usf[1], der[1], der[2], der[3], 1, 0, 0) + call xgsder (usf[1], der[1], der[2], der[4], 1, 1, 0) + call xgsder (usf[1], der[1], der[2], der[5], 1, 0, 1) + do j = 2, nusf { + call xgsder (usf[j], der[1], der[2], tmp[1], 1, 0, 0) + call xgsder (usf[j], der[1], der[2], tmp[2], 1, 1, 0) + call xgsder (usf[j], der[1], der[2], tmp[3], 1, 0, 1) + der[3] = der[3] + tmp[1] + der[4] = der[4] + tmp[2] + der[5] = der[5] + tmp[3] + } + der[3] = der[3] / nusf + der[4] = der[4] / nusf + der[5] = der[5] / nusf + } + + if (nvsf == 0) + der[6] = der[2] + else if (nvsf == 1) { + call xgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0) + call xgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0) + call xgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1) + } else { + call xgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0) + call xgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0) + call xgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1) + do j = 2, nvsf { + call xgsder (vsf[j], der[1], der[2], tmp[1], 1, 0, 0) + call xgsder (vsf[j], der[1], der[2], tmp[2], 1, 1, 0) + call xgsder (vsf[j], der[1], der[2], tmp[3], 1, 0, 1) + der[6] = der[6] + tmp[1] + der[7] = der[7] + tmp[2] + der[8] = der[8] + tmp[3] + } + der[6] = der[6] / nvsf + der[7] = der[7] / nvsf + der[8] = der[8] / nvsf + } + } +end + + +# TR_INIT -- Since the inversion iteration always begins from the last +# point we need to initialize before the first call to TR_INVERT. +# When using a WCS this simply returns. + +procedure tr_init (ct, usf, nusf, vsf, nvsf, x, y, der) + +pointer ct #I CT pointer +pointer usf[ARB], vsf[ARB] #I User coordinate surfaces +int nusf, nvsf #I Number of surfaces for each coordinate +real x, y #I Starting X and Y +real der[8] #O Inversion data + +int j +real tmp[3] + +begin + if (ct != NULL) + return + + der[1] = x + der[2] = y + if (nusf == 0) { + der[3] = der[1] + der[4] = 1. + der[5] = 0. + } else if (nusf == 1) { + call xgsder (usf[1], der[1], der[2], der[3], 1, 0, 0) + call xgsder (usf[1], der[1], der[2], der[4], 1, 1, 0) + call xgsder (usf[1], der[1], der[2], der[5], 1, 0, 1) + } else { + call xgsder (usf[1], der[1], der[2], der[3], 1, 0, 0) + call xgsder (usf[1], der[1], der[2], der[4], 1, 1, 0) + call xgsder (usf[1], der[1], der[2], der[5], 1, 0, 1) + do j = 2, nusf { + call xgsder (usf[j], der[1], der[2], tmp[1], 1, 0, 0) + call xgsder (usf[j], der[1], der[2], tmp[2], 1, 1, 0) + call xgsder (usf[j], der[1], der[2], tmp[3], 1, 0, 1) + der[3] = der[3] + tmp[1] + der[4] = der[4] + tmp[2] + der[5] = der[5] + tmp[3] + } + der[3] = der[3] / nusf + der[4] = der[4] / nusf + der[5] = der[5] / nusf + } + + if (nvsf == 0) { + der[6] = der[2] + der[7] = 0. + der[8] = 1. + } else if (nvsf == 1) { + call xgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0) + call xgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0) + call xgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1) + } else { + call xgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0) + call xgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0) + call xgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1) + do j = 2, nvsf { + call xgsder (vsf[j], der[1], der[2], tmp[1], 1, 0, 0) + call xgsder (vsf[j], der[1], der[2], tmp[2], 1, 1, 0) + call xgsder (vsf[j], der[1], der[2], tmp[3], 1, 0, 1) + der[6] = der[6] + tmp[1] + der[7] = der[7] + tmp[2] + der[8] = der[8] + tmp[3] + } + der[6] = der[6] / nvsf + der[7] = der[7] / nvsf + der[8] = der[8] / nvsf + } +end + + +# TR_EVAL -- Evalute coordinate function. +# +# This is an interface routine to allow using either an MWCS CT (coordinate +# transform) pointer or a GSURFIT SF (2D surface function) pointer. The +# surface method is used with a FITCOORDS database. The MWCS method is +# used to retransform an image with a WCS. + +real procedure tr_eval (ct, sf, axis, x, y) + +pointer ct #I CT pointer +pointer sf #I SF pointer +int axis #I World coordinate axis to return +real x, y #I Pixel coordinate to transform + +real w[2], xgseval() + +begin + if (sf != NULL) + return (xgseval (sf, x, y)) + + call mw_c2tranr (CT_LW(ct), x, y, w[1], w[2]) + return (w[axis]) +end + + +# TR_GETR -- Get real valued parameter. +# +# This is an interface routine to allow using either an MWCS CT (coordinate +# transform) pointer or a GSURFIT SF (2D surface function) pointer. The +# surface method is used with a FITCOORDS database. The MWCS method is +# used to retransform an image with a WCS. + +real procedure tr_getr (ct, sf, param) + +pointer ct #I CT pointer +pointer sf #I SF pointer +int param #I Parameter code + +real xgsgetr() + +begin + if (sf != NULL) + return (xgsgetr (sf, param)) + + switch (param) { + case GSXMIN, GSYMIN: + return (real (1)) + case GSXMAX: + return (real (CT_NX(ct))) + case GSYMAX: + return (real (CT_NY(ct))) + } +end |