aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/tv/wcslab/wlutil.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/tv/wcslab/wlutil.x')
-rw-r--r--pkg/images/tv/wcslab/wlutil.x390
1 files changed, 390 insertions, 0 deletions
diff --git a/pkg/images/tv/wcslab/wlutil.x b/pkg/images/tv/wcslab/wlutil.x
new file mode 100644
index 00000000..c79b8f5e
--- /dev/null
+++ b/pkg/images/tv/wcslab/wlutil.x
@@ -0,0 +1,390 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imio.h>
+include <imhdr.h>
+include <gset.h>
+include <math.h>
+
+# WL_IMD_VIEWPORT -- Map the viewport and window of the image display.
+
+procedure wl_imd_viewport (frame, im, c1, c2, l1, l2, vl, vr, vb, vt)
+
+int frame # I: display frame to be overlayed
+pointer im # I: pointer to the input image
+real c1, c2, l1, l2 # I/O: input/output window
+real vl, vr, vb, vt # I/O: input/output viewport
+
+int wcs_status, dim1, dim2, step1, step2
+pointer sp, frimage, frim, iw
+real x1, x2, y1, y2, fx1, fx2, fy1, fy2, junkx, junky
+real vx1, vx2, vy1, vy2, nx1, nx2, ny1, ny2
+pointer imd_mapframe(), iw_open()
+
+
+begin
+ # If all of the viewport parameters were defined by the user
+ # use the default viewport and window.
+ if (! IS_INDEFR(vl) && ! IS_INDEFR(vr) && ! IS_INDEFR(vb) &&
+ ! IS_INDEFR(vt))
+ return
+
+ # Allocate some memory.
+ call smark (sp)
+ call salloc (frimage, SZ_FNAME, TY_CHAR)
+
+ # Open the requested display frame and get the loaded image name.
+ # If this name is blank, use the default viewport and window.
+
+ frim = imd_mapframe (frame, READ_ONLY, YES)
+ iw = iw_open (frim, frame, Memc[frimage], SZ_FNAME, wcs_status)
+ if (Memc[frimage] == EOS || wcs_status == ERR) {
+ call iw_close (iw)
+ call imunmap (frim)
+ call sfree (sp)
+ return
+ }
+
+ # Find the beginning and end points of the requested image section.
+ # We already know at this point that the input logical image is
+ # 2-dimensional. However this 2-dimensional section may be part of
+ # n-dimensional image.
+
+ # X dimension.
+ dim1 = IM_VMAP(im,1)
+ step1 = IM_VSTEP(im,1)
+ if (step1 >= 0) {
+ x1 = IM_VOFF(im,dim1) + 1
+ x2 = x1 + IM_LEN(im,1) - 1
+ } else {
+ x1 = IM_VOFF(im,dim1) - 1
+ x2 = x1 - IM_LEN(im,1) + 1
+ }
+
+ # Y dimension.
+ dim2 = IM_VMAP(im,2)
+ step2 = IM_VSTEP(im,2)
+ if (step2 >= 0) {
+ y1 = IM_VOFF(im,dim2) + 1
+ y2 = y1 + IM_LEN(im,2) - 1
+ } else {
+ y1 = IM_VOFF(im,dim2) - 1
+ y2 = y1 - IM_LEN(im,2) + 1
+ }
+
+ # Get the frame buffer coordinates corresponding to the lower left
+ # and upper right corners of the image section.
+
+ call iw_im2fb (iw, x1, y1, fx1, fy1)
+ call iw_im2fb (iw, x2, y2, fx2, fy2)
+ if (fx1 > fx2) {
+ junkx = fx1
+ fx1 = fx2
+ fx2 = junkx
+ }
+ if (fy1 > fy2) {
+ junky = fy1
+ fy1 = fy2
+ fy2 = junky
+ }
+
+ # Check that some portion of the input image is in the display.
+ # If not select the default viewport and window coordinates.
+ if (fx1 > IM_LEN(frim,1) || fx2 < 1.0 || fy1 > IM_LEN(frim,2) ||
+ fy2 < 1.0) {
+ call iw_close (iw)
+ call imunmap (frim)
+ call sfree (sp)
+ return
+ }
+
+ # Compute a new viewport and window for X.
+ if (fx1 >= 1.0) {
+ vx1 = max (0.0, min (1.0, (fx1 - 0.5) / IM_LEN(frim,1)))
+ nx1 = 1.0
+ } else {
+ vx1 = 0.0
+ call iw_fb2im (iw, 1.0, 1.0, junkx, junky)
+ if (step1 >= 0)
+ nx1 = max (1.0, junkx - x1 + 1.0)
+ else
+ nx2 = max (1.0, junkx - x2 + 1.0)
+ }
+ if (fx2 <= IM_LEN(frim,1)) {
+ vx2 = max (0.0, min (1.0, (fx2 + 0.5) / IM_LEN(frim,1)))
+ nx2 = IM_LEN(im,1)
+ } else {
+ vx2 = 1.0
+ call iw_fb2im (iw, real(IM_LEN(frim,1)), real (IM_LEN(frim,2)),
+ junkx, junky)
+ if (step1 >= 0)
+ nx2 = min (real (IM_LEN(im,1)), junkx - x1 + 1.0)
+ else
+ nx1 = min (real (IM_LEN(im,1)), junkx - x2 + 1.0)
+ }
+
+ # Compute a new viewport and window for Y.
+ if (fy1 >= 1.0) {
+ vy1 = max (0.0, min (1.0, (fy1 - 0.5) / IM_LEN(frim,2)))
+ ny1 = 1.0
+ } else {
+ vy1 = 0.0
+ call iw_fb2im (iw, 1.0, 1.0, junkx, junky)
+ if (step2 >= 0)
+ ny1 = max (1.0, junky - y1 + 1)
+ else
+ ny2 = max (1.0, junky - y2 + 1)
+ }
+ if (fy2 <= IM_LEN(frim,2)) {
+ vy2 = max (0.0, min (1.0, (fy2 + 0.5) / IM_LEN(frim,2)))
+ ny2 = IM_LEN(im,2)
+ } else {
+ vy2 = 1.0
+ call iw_fb2im (iw, real (IM_LEN(frim,1)), real (IM_LEN(frim,2)),
+ junkx, junky)
+ if (step2 >= 0)
+ ny2 = min (real (IM_LEN(im,2)), junky - y1 + 1.0)
+ else
+ ny1 = min (real (IM_LEN(im,2)), junky - y2 + 1.0)
+ }
+
+ # Define a the new viewport and window.
+ if (IS_INDEFR(vl)) {
+ vl = vx1
+ c1 = nx1
+ }
+ if (IS_INDEFR(vr)) {
+ vr = vx2
+ c2 = nx2
+ }
+ if (IS_INDEFR(vb)) {
+ vb = vy1
+ l1 = ny1
+ }
+ if (IS_INDEFR(vt)) {
+ vt = vy2
+ l2 = ny2
+ }
+
+ # Clean up.
+ call iw_close (iw)
+ call imunmap (frim)
+ call sfree (sp)
+end
+
+
+define EDGE1 0.1
+define EDGE2 0.9
+define EDGE3 0.12
+define EDGE4 0.92
+
+# WL_MAP_VIEWPORT -- Set device viewport wcslab plots. If not specified by
+# user, a default viewport centered on the device is used.
+
+procedure wl_map_viewport (gp, c1, c2, l1, l2, ux1, ux2, uy1, uy2, fill)
+
+pointer gp # I: pointer to graphics descriptor
+real c1, c2, l1, l2 # I: the column and line limits
+real ux1, ux2, uy1, uy2 # I/O: NDC coordinates of requested viewort
+bool fill # I: fill viewport (vs preserve aspect ratio)
+
+int ncols, nlines
+real xcen, ycen, ncolsr, nlinesr, ratio, aspect_ratio
+real x1, x2, y1, y2, ext, xdis, ydis
+bool fp_equalr()
+real ggetr()
+data ext /0.0625/
+
+begin
+ ncols = nint (c2 - c1) + 1
+ ncolsr = real (ncols)
+ nlines = nint (l2 - l1) + 1
+ nlinesr = real (nlines)
+
+ # Determine the standard window sizes.
+ if (fill) {
+ x1 = 0.0; x2 = 1.0
+ y1 = 0.0; y2 = 1.0
+ } else {
+ x1 = EDGE1; x2 = EDGE2
+ y1 = EDGE3; y2 = EDGE4
+ }
+
+ # If any values were specified, then replace them here.
+ if (! IS_INDEFR(ux1))
+ x1 = ux1
+ if (! IS_INDEFR(ux2))
+ x2 = ux2
+ if (! IS_INDEFR(uy1))
+ y1 = uy1
+ if (! IS_INDEFR(uy2))
+ y2 = uy2
+
+ # Calculate optimum viewport, as in NCAR's conrec, hafton.
+ if (! fill) {
+ ratio = min (ncolsr, nlinesr) / max (ncolsr, nlinesr)
+ if (ratio >= ext) {
+ if (ncols > nlines)
+ y2 = (y2 - y1) * nlinesr / ncolsr + y1
+ else
+ x2 = (x2 - x1) * ncolsr / nlinesr + x1
+ }
+ }
+
+ xdis = x2 - x1
+ ydis = y2 - y1
+ xcen = (x2 + x1) / 2.
+ ycen = (y2 + y1) / 2.
+
+ # So far, the viewport has been calculated so that equal numbers of
+ # image pixels map to equal distances in NDC space, regardless of
+ # the aspect ratio of the device. If the parameter "fill" has been
+ # set to no, the user wants to compensate for a non-unity aspect
+ # ratio and make equal numbers of image pixels map to into the same
+ # physical distance on the device, not the same NDC distance.
+
+ if (! fill) {
+ aspect_ratio = ggetr (gp, "ar")
+ if (fp_equalr (aspect_ratio, 0.0))
+ aspect_ratio = 1.0
+
+ if (aspect_ratio < 1.0)
+ # Landscape
+ xdis = xdis * aspect_ratio
+ else if (aspect_ratio > 1.0)
+ # Portrait
+ ydis = ydis / aspect_ratio
+ }
+
+ ux1 = xcen - (xdis / 2.0)
+ ux2 = xcen + (xdis / 2.0)
+ uy1 = ycen - (ydis / 2.0)
+ uy2 = ycen + (ydis / 2.0)
+
+ call gsview (gp, ux1, ux2, uy1, uy2)
+ call gswind (gp, c1, c2, l1, l2)
+end
+
+
+# WL_W2LD -- Transform world coordinates to logical coordinates.
+
+procedure wl_w2ld (wlct, flip, wx, wy, lx, ly, npts)
+
+pointer wlct # I: the MWCS coordinate transformation descriptor
+int flip # I: true if the axes are transposed
+double wx[npts], wy[npts] # I: the world coordinates
+double lx[npts], ly[npts] # O: the logical coordinates
+int npts # I: the number of points to translate
+
+begin
+ if (flip == YES)
+ call mw_v2trand (wlct, wx, wy, ly, lx, npts)
+ else
+ call mw_v2trand (wlct, wx, wy, lx, ly, npts)
+end
+
+
+# WL_L2WD -- Transform logical coordinates to world coordinates.
+
+procedure wl_l2wd (lwct, flip, lx, ly, wx, wy, npts)
+
+pointer lwct # I: the MWCS coordinate transformation descriptor
+int flip # I: true if the axes are transposed
+double lx[npts], ly[npts] # I: the logical coordinates
+double wx[npts], wy[npts] # O: the world coordinates
+int npts # I: the number of points to translate
+
+begin
+ if (flip == YES)
+ call mw_v2trand (lwct, ly, lx, wx, wy, npts)
+ else
+ call mw_v2trand (lwct, lx, ly, wx, wy, npts)
+end
+
+
+# WL_MAX_ELEMENT_ARRAY -- Return the index of the maximum array element.
+#
+# Description
+# This function returns the index of the maximum value of the input array.
+
+int procedure wl_max_element_array (array, npts)
+
+double array[ARB] # I: the array to look through for the maximum
+int npts # I: the number of points in the array
+
+int i, maximum
+
+begin
+ maximum = 1
+ for (i = 2; i <= npts; i = i + 1)
+ if (array[i] > array[maximum])
+ maximum = i
+
+ return (maximum)
+end
+
+
+# WL_DISTANCED - Determine the distance between two points.
+
+double procedure wl_distanced (x1, y1, x2, y2)
+
+double x1, y1 # I: coordinates of point 1
+double x2, y2 # I: coordinates of point 2
+
+double a, b
+
+begin
+ a = x1 - x2
+ b = y1 - y2
+ return (sqrt ((a * a) + (b * b)))
+end
+
+
+# WL_DISTANCER -- Determine the distance between two points.
+
+real procedure wl_distancer (x1, y1, x2, y2)
+
+real x1, y1 # I: coordinates of point 1
+real x2, y2 # I: coordinates of point 2
+
+real a, b
+
+begin
+ a = x1 - x2
+ b = y1 - y2
+ return (sqrt ((a * a) + (b * b)))
+end
+
+
+# The dimensionality.
+define N_DIM 2
+
+# Define some memory management.
+define ONER Memr[$1+$2-1]
+
+# WL_ROTATE -- Rotate a vector.
+
+procedure wl_rotate (x, y, npts, angle, nx, ny)
+
+real x[npts], y[npts] # I: the vectors to rotate
+int npts # I: the number of points in the vectors
+real angle # I: the angle to rotate (radians)
+real nx[npts], ny[npts] # O: the transformed vectors
+
+pointer sp, center, mw
+pointer mw_open(), mw_sctran()
+
+begin
+ # Get some memory.
+ call smark (sp)
+ call salloc (center, N_DIM, TY_REAL)
+
+ mw = mw_open (NULL, N_DIM)
+ ONER(center,1) = 0.
+ ONER(center,2) = 0.
+ call mw_rotate (mw, -DEGTORAD( angle ), ONER(center,1), 3b)
+ call mw_v2tranr (mw_sctran (mw, "physical", "logical", 3b),
+ x, y, nx, ny, npts)
+
+ call mw_close (mw)
+ call sfree (sp)
+end
generated by cgit (git 2.34.1) at 2025-09-20 09:20:59 -0400