From 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 Mon Sep 17 00:00:00 2001 From: Joe Hunkeler Date: Tue, 11 Aug 2015 16:51:37 -0400 Subject: Repatch (from linux) of OSX IRAF --- pkg/plot/t_surface.x | 501 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 501 insertions(+) create mode 100644 pkg/plot/t_surface.x (limited to 'pkg/plot/t_surface.x') diff --git a/pkg/plot/t_surface.x b/pkg/plot/t_surface.x new file mode 100644 index 00000000..7faa38a2 --- /dev/null +++ b/pkg/plot/t_surface.x @@ -0,0 +1,501 @@ +include +include +include +include +include +include +include + +define DUMMY 6 +define SZ_TLABEL 10 +define CSIZE 24 + + +# SURFACE -- Draw a perspective view of an image section. The altitude +# and azimuth of the viewing angle are variable. Floor and ceiling +# constraints may be applied to the image data before plotting if desired. + +procedure t_surface() + +char imsect[SZ_FNAME] +char device[SZ_FNAME], title[SZ_LINE] +bool label, sub, pre +pointer im, subras, work +int ncols, nlines, mode, wkid, nx, ny, npix +int epa, status, old_onint, tsujmp[LEN_JUMPBUF] +int xres, yres, first +real angh, angv, imcols, imlines +real floor, ceiling, vpx1, vpx2, vpy1, vpy2 + +pointer gp, gopen() +bool clgetb(), streq() +int clgeti(), surf_limits() +real clgetr() +extern tsu_onint() +pointer immap(), plt_getdata() +common /tsucom/ tsujmp +common /noaovp/ vpx1, vpx2, vpy1, vpy2 +common /frstfg/ first + +define exit_ 91 +begin + # First initialize srface common blocks before changing any parameters + first = 1 + call srfabd + + # Get image section string and output device. + call clgstr ("image", imsect, SZ_FNAME) + call clgstr ("device", device, SZ_FNAME) + + # Map image and open graphics device. + im = immap (imsect, READ_ONLY, 0) + + angh = clgetr ("angh") + angv = clgetr ("angv") + floor = clgetr ("floor") + ceiling = clgetr ("ceiling") + label = clgetb ("label") + + call clgstr ("title", title, SZ_LINE) + if (streq (title, "imtitle")) { + call strcpy (imsect, title, SZ_LINE) + call strcat (": ", title, SZ_LINE) + call strcat (IM_TITLE(im), title, SZ_LINE) + } + + # If a label is to be drawn, don't use the full device viewport for + # the surface plot. This doesn't allow room for the axes and labels. + + if (label) { + vpx1 = 0.10 + vpx2 = 0.90 + vpy1 = 0.10 + vpy2 = 0.90 + } + + mode = NEW_FILE + if (clgetb ("append")) + mode = APPEND + + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + imcols = real (ncols) + imlines = real (nlines) + + xres = clgeti ("xres") + yres = clgeti ("yres") + sub = clgetb ("subsample") + pre = clgetb ("preserve") + + # Get data with proper resolution. Procedure plt_getdata returns a + # pointer to the data matrix to be contoured. The resolution is + # decreased by the specified method in this procedure. The image + # header pointer can be unmapped after plt_getdata is called. + + nx = 0 + ny = 0 + subras = plt_getdata (im, sub, pre, xres, yres, nx, ny) + call imunmap (im) + + # Allocate the working storage needed by EZSRFC. + #call malloc (work, (2 * nx * ny) + nx + ny, TY_REAL) + call malloc (work, 2 * ((2 * nx * ny) + nx + ny), TY_REAL) + + # Take floor and ceiling if enabled (nonzero). + npix = nx * ny + if (surf_limits (Memr[subras], npix, floor, ceiling) == ERR) + goto exit_ + + # Open graphics device and make plot. + call gopks (STDERR) + wkid = 1 + gp = gopen (device, mode, STDGRAPH) + call gopwk (wkid, DUMMY, gp) + call gacwk (wkid) + call gtext (gp, 0.5, .96, title, "s=0.8;f=b;h=c") + call set (vpx1, vpx2, vpy1, vpy2, 1.0, 1024., 1.0, 1024., 1) + + # Install interrupt exception handler. + call zlocpr (tsu_onint, epa) + call xwhen (X_INT, epa, old_onint) + + call zsvjmp (tsujmp, status) + if (status == OK) + call ezsrfc (Memr[subras], nx, ny, angh, angv, Memr[work]) + else { + call gcancel (gp) + call fseti (STDOUT, F_CANCEL, OK) + } + + if (label) { + # Establish plotting window in full scale image coordinates. + call gswind (gp, 1.0, imcols, 1.0, imlines) + call gseti (gp, G_CLIP, NO) + call srf_perimeter (gp, Memr[subras], nx, ny, angh, angv) + } + + call gdawk (wkid) + call gclwk (wkid) + call gclks () +exit_ + call mfree (subras, TY_REAL) + call mfree (work, TY_REAL) + +end + + +# TSU_ONINT -- Interrupt handler for the task surface. Branches back to ZSVJMP +# in the main routine to permit shutdown without an error message. + +procedure tsu_onint (vex, next_handler) + +int vex # virtual exception +int next_handler # not used + +int tsujmp[LEN_JUMPBUF] +common /tsucom/ tsujmp + +begin + call xer_reset() + call zdojmp (tsujmp, vex) +end + + +# SURF_LIMITS -- Apply the floor and ceiling constraints to the subraster. +# If both values are exactly zero, they are not applied. + +int procedure surf_limits (ras, npix, floor, ceiling) + +real ras[npix] # Input array of pixels +int npix # npixels in array +real floor, ceiling # user specified parameters +int apply, i +real tfloor, tceiling, delta, dmin, dmax +bool fp_equalr() + +begin + tfloor = floor + tceiling = ceiling + apply = YES + + call alimr (ras, npix, dmin, dmax) + if (fp_equalr (dmin, dmax)) { + call eprintf ("Constant valued array; no plot drawn\n") + return (ERR) + } + + if (fp_equalr (tfloor, INDEF)) + tfloor = dmin + if (fp_equalr (tceiling, INDEF)) + tceiling = dmax + + delta = tceiling - tfloor + if (delta < 0.0) { + # specified ceiling is lower than floor, flip them + floor = tceiling + ceiling = tfloor + } else if (fp_equalr (delta, 0.0)) { + # degenerate values + apply = NO + floor = dmin + ceiling = dmax + call eprintf ( + "Floor and ceiling are degenerate values and will be ignored\n") + } else { + # Non-degenerate, ceiling exceedes floor as expected + floor = tfloor + ceiling = tceiling + } + + if (apply == YES) { + # First verify that floor and ceiling are valid + if (dmax <= floor) { + call eprintf ("Entire image is at or below specified floor\n") + return (ERR) + } + if (dmin >= ceiling) { + call eprintf ("Entire image is at or above specified ceiling\n") + return (ERR) + } + + do i = 1, npix { + # Apply surface limits + ras[i] = max (floor, ras[i]) + ras[i] = min (ceiling, ras[i]) + } + } + return (OK) +end + + +# SRF_PERIMETER -- draw and label axes around the surface plot. + +procedure srf_perimeter (gp, z, ncols, nlines, angh, angv) + +pointer gp # Graphics pointer +int ncols # Number of image columns +int nlines # Number of image lines +real z[ncols, nlines] # Array of intensity values +real angh # Angle of horizontal inclination +real angv # Angle of vertical inclination + +pointer sp, x_val, y_val, kvec +char tlabel[SZ_TLABEL] +real xmin, ymin, delta, fact1, flo, hi, xcen, ycen +real x1_perim, x2_perim, y1_perim, y2_perim, z1, z2 +real wc1, wc2, wl1, wl2, del +int i, j, junk +int itoc() +data fact1 /2.0/ +real vpx1, vpx2, vpy1, vpy2 +common /noaovp/ vpx1, vpx2, vpy1, vpy2 + +begin + call smark (sp) + call salloc (x_val, ncols + 2, TY_REAL) + call salloc (y_val, nlines + 2, TY_REAL) + call salloc (kvec, max (ncols, nlines) + 2, TY_REAL) + + # Get window coordinates set up in calling procedure. + call ggwind (gp, wc1, wc2, wl1, wl2) + + # Set up window, viewport for output. The coordinates returned + # from trn32s are in the range [1-1024]. + call set (vpx1, vpx2, vpy1, vpy2, 1.0, 1024., 1.0, 1024., 1) + + # Find range of z for determining perspective + flo = MAX_REAL + hi = -flo + do j = 1, nlines { + call alimr (z[1,j], ncols, z1, z2) + flo = min (flo, z1) + hi = max (hi, z2) + } + + # Set up linear endpoints and spacing as used in surface. + + delta = (hi-flo) / (max (ncols,nlines) -1.) * fact1 + xmin = -(real (ncols/2) * delta + real (mod (ncols+1, 2)) * delta) + ymin = -(real (nlines/2) * delta + real (mod (nlines+1, 2)) * delta) + del = 2.0 * delta + + # The perimeter is separated from the surface plot by the + # width of delta. + + x1_perim = xmin - delta + y1_perim = ymin - delta + x2_perim = xmin + (real (ncols) * delta) + y2_perim = ymin + (real (nlines) * delta) + # Set up linear arrays over full perimeter range + do i = 1, ncols + 2 + Memr[x_val+i-1] = x1_perim + (i-1) * delta + do i = 1, nlines + 2 + Memr[y_val+i-1] = y1_perim + (i-1) * delta + + # Draw and label axes and tick marks. + # It is important that frame has not been called after calling srface. + # First to draw the perimeter. Which axes get drawn depends on the + # values of angh and angv. Get angles in the range [-180, 180]. + + if (angh > 180.) + angh = angh - 360. + else if (angh < -180.) + angh = angh + 360. + + if (angv > 180.) + angv = angv - 360. + else if (angv < -180.) + angv = angv + 360. + + # Calculate positions for the axis labels + xcen = 0.5 * (x1_perim + x2_perim) + ycen = 0.5 * (y1_perim + y2_perim) + + if (angh >= 0) { + if (angv >= 0) { + # Case 1: xy rotation positive, looking down from above mid Z + + # First draw x axis + call amovkr (y2_perim, Memr[kvec], ncols + 2) + call srf_draw_axis (Memr[x_val+1], Memr[kvec], flo, ncols + 1) + call srf_label_axis (xcen, y2_perim+del, flo, "X-AXIS", -1, -2) + call srf_draw_ticksx (Memr[x_val+1], y2_perim, y2_perim+delta, + flo, ncols) + call srf_label_axis (xmin, y2_perim+del, flo, "1", -1, -2) + junk = itoc (int (wc2), tlabel, SZ_TLABEL) + call srf_label_axis (Memr[x_val+ncols], y2_perim+del, flo, + tlabel, -1, -2) + + # Now draw y axis + call amovkr (x2_perim, Memr[kvec], nlines + 2) + call srf_draw_axis (Memr[kvec], Memr[y_val+1], flo, nlines + 1) + call srf_label_axis (x2_perim+del, ycen, flo, "Y-AXIS", 2, -1) + call srf_draw_ticksy (x2_perim, x2_perim+delta, Memr[y_val+1], + flo, nlines) + call srf_label_axis (x2_perim+del, ymin, flo, "1", 2, -1) + junk = itoc (int (wl2), tlabel, SZ_TLABEL) + call srf_label_axis (x2_perim+del, Memr[y_val+nlines], flo, + tlabel, 2, -1) + } else { + # Case 2: xy rotation positive, looking up from below mid Z + # First draw x axis + call amovkr (y1_perim, Memr[kvec], ncols + 2) + call srf_draw_axis (Memr[x_val], Memr[kvec], flo, ncols + 1) + call srf_label_axis (xcen, y1_perim-del, flo, "X-AXIS", -1, 2) + call srf_draw_ticksx (Memr[x_val+1], y1_perim, y1_perim-delta, + flo, ncols) + call srf_label_axis (xmin, y1_perim-del, flo, "1", -1, 2) + junk = itoc (int (wc2), tlabel, SZ_TLABEL) + call srf_label_axis (Memr[x_val+ncols], y1_perim-del, flo, + tlabel, -1, 2) + + # Now draw y axis + call amovkr (x1_perim, Memr[kvec], nlines + 2) + call srf_draw_axis (Memr[kvec], Memr[y_val], flo, nlines + 1) + call srf_label_axis (x1_perim-del, ycen, flo, "Y-AXIS", 2, 1) + call srf_draw_ticksy (x1_perim, x1_perim-delta, Memr[y_val+1], + flo, nlines) + call srf_label_axis (x1_perim-del, ymin, flo, "1", 2, 1) + junk = itoc (int (wl2), tlabel, SZ_TLABEL) + call srf_label_axis (x1_perim-del, Memr[y_val+nlines], flo, + tlabel, 2, 1) + } + } + + if (angh < 0) { + if (angv > 0) { + # Case 3: xy rotation negative, looking down from above mid Z + # (default). First draw x axis + call amovkr (y1_perim, Memr[kvec], ncols + 2) + call srf_draw_axis (Memr[x_val+1], Memr[kvec], flo, ncols + 1) + call srf_label_axis (xcen, y1_perim-del, flo, "X-AXIS", 1, 2) + call srf_draw_ticksx (Memr[x_val+1], y1_perim, y1_perim-delta, + flo, ncols) + call srf_label_axis (xmin, y1_perim-del, flo, "1", 1, 2) + junk = itoc (int (wc2), tlabel, SZ_TLABEL) + call srf_label_axis (Memr[x_val+ncols], y1_perim-del, flo, + tlabel, 1, 2) + + # Now draw y axis + call amovkr (x2_perim, Memr[kvec], nlines + 2) + call srf_draw_axis (Memr[kvec], Memr[y_val], flo, nlines + 1) + call srf_label_axis (x2_perim+del, ycen, flo, "Y-AXIS", 2, -1) + call srf_draw_ticksy (x2_perim, x2_perim+delta, Memr[y_val+1], + flo, nlines) + call srf_label_axis (x2_perim+del, ymin, flo, "1", 2, -1) + junk = itoc (int (wl2), tlabel, SZ_TLABEL) + call srf_label_axis (x2_perim+del, Memr[y_val+nlines], flo, + tlabel, 2, -1) + } else { + # Case 4: xy rotation negative, looking up from below mid Z + # First draw x axis + call amovkr (y2_perim, Memr[kvec], ncols + 2) + call srf_draw_axis (Memr[x_val], Memr[kvec], flo, ncols + 1) + call srf_label_axis (xcen, y2_perim+del, flo, "X-AXIS", 1, -2) + call srf_draw_ticksx (Memr[x_val+1], y2_perim, y2_perim+delta, + flo, ncols) + call srf_label_axis (xmin, y2_perim+del, flo, "1", 1, -2) + junk = itoc (int (wc2), tlabel, SZ_TLABEL) + call srf_label_axis (Memr[x_val+ncols], y2_perim+del, flo, + tlabel, 1, -2) + + # Now draw y axis + call amovkr (x1_perim, Memr[kvec], nlines + 2) + call srf_draw_axis (Memr[kvec], Memr[y_val+1], flo, nlines + 1) + call srf_label_axis (x1_perim-del, ycen, flo, "Y-AXIS", 2, 1) + call srf_draw_ticksy (x1_perim, x1_perim-delta, Memr[y_val+1], + flo, nlines) + call srf_label_axis (x1_perim-del, ymin, flo, "1", 2, 1) + junk = itoc (int (wl2), tlabel, SZ_TLABEL) + call srf_label_axis (x1_perim-del, Memr[y_val+nlines], flo, + tlabel, 2, 1) + } + } + + # Flush plotit buffer before returning + call plotit (0, 0, 2) + call sfree (sp) +end + + +procedure srf_draw_axis (xvals, yvals, zval, nvals) + +int nvals +real xvals[nvals] +real yvals[nvals] +real zval +pointer sp, xt, yt +int i +real dum + +begin + call smark (sp) + call salloc (xt, nvals, TY_REAL) + call salloc (yt, nvals, TY_REAL) + + do i = 1, nvals + call trn32s (xvals[i], yvals[i], zval, Memr[xt+i-1], Memr[yt+i-1], + dum, 1) + + call gpl (nvals, Memr[xt], Memr[yt]) + call sfree (sp) +end + + +procedure srf_label_axis (xval, yval, zval, sppstr, path, up) + +real xval +real yval +real zval +char sppstr[SZ_LINE] +int path +int up + +int nchars +int strlen() +% character*64 fstr + +begin + nchars = strlen (sppstr) + +% call f77pak (sppstr, fstr, 64) + call pwrzs (xval, yval, zval, fstr, nchars, CSIZE, path, up, 0) +end + + +procedure srf_draw_ticksx (x, y1, y2, zval, nvals) + +int nvals +real x[nvals] +real y1, y2 +real zval + +int i +real tkx[2], tky[2], dum + +begin + do i = 1, nvals { + call trn32s (x[i], y1, zval, tkx[1], tky[1], dum, 1) + call trn32s (x[i], y2, zval, tkx[2], tky[2], dum, 1) + call gpl (2, tkx[1], tky[1]) + } +end + + +procedure srf_draw_ticksy (x1, x2, y, zval, nvals) + +int nvals +real x1, x2 +real y[nvals] +real zval + +int i +real tkx[2], tky[2], dum + +begin + do i = 1, nvals { + call trn32s (x1, y[i], zval, tkx[1], tky[1], dum, 1) + call trn32s (x2, y[i], zval, tkx[2], tky[2], dum, 1) + call gpl (2, tkx[1], tky[1]) + } +end -- cgit