Initial commit

1 files changed, 335 insertions, 0 deletions

diff --git a/pkg/images/tv/iis/src/load2.x b/pkg/images/tv/iis/src/load2.x
new file mode 100644
index 00000000..5372907f
--- /dev/null
+++ b/pkg/images/tv/iis/src/load2.x
@@ -0,0 +1,335 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+#### load2.x (from load.x) ####
+
+include <mach.h>
+include <imset.h>
+include <imhdr.h>
+include	<error.h>
+include	<gki.h>
+include <fio.h>
+include <fset.h>
+include	"gwindow.h"
+include	"../lib/ids.h"
+include	"cv.h"
+
+# DS_LOAD_DISPLAY -- Map an image into the display window.  In general this
+#   involves independent linear transformations in the X, Y, and Z (greyscale)
+#   dimensions.  If a spatial dimension is larger than the display window then
+#   the image is block averaged.  If a spatial dimension or a block averaged
+#   dimension is smaller than the display window then linear interpolation is
+#   used to expand the image.  Both the input image and the output device appear
+#   to us as images, accessed via IMIO.
+#
+# World coordinate system 0 (WCS 0) defines the position and size of the device
+#   window in NDC coordinates (0-1 in either axis).  WCS 1 assigns a pixel
+#   coordinate system to the same window.  If we convert the NDC coordinates of
+#   the window into device coordinates in pixels, then the ratios of the window
+#   coordinates in pixels to the image coordinates in pixels defines the real
+#   magnification factors for the two spatial axes.  If the pixel coordinates
+#   are out of bounds then the image will be displayed centered in the window
+#   with zero fill at the edges.  If the frame has not been erased then the fill
+#   areas must be explicitly zeroed.
+
+procedure ds_load_display (im, wdes, border_erase)
+
+pointer	im			# input image
+pointer	wdes			# graphics window descriptor
+bool	border_erase
+
+int	wx1, wx2, wy1, wy2	# device window to be filled with image data
+real	px1, px2, py1, py2	# image coords in fractional image pixels
+real	pxsize, pysize		# size of image section in fractional pixels
+real	wxcenter, wycenter	# center of device window in frac device pixels
+real	xmag, ymag		# x,y magnification ratios
+pointer	w0, w1			# world coord systems 0 (NDC) and 1 (pixel)
+
+include	"cv.com"
+
+begin
+	# Compute pointers to WCS 0 and 1.
+	w0 = W_WC(wdes,0)
+	w1 = W_WC(wdes,1)
+
+	# Compute X and Y magnification ratios required to map image into
+	# the device window in device pixel units.
+
+	xmag = (W_XE(w0) - W_XS(w0)) * cv_xres / (W_XE(w1) - W_XS(w1))
+	ymag = (W_YE(w0) - W_YS(w0)) * cv_yres / (W_YE(w1) - W_YS(w1))
+
+	# Compute the coordinates of the image section to be displayed.
+	# This is not necessarily the same as WCS 1 since the WCS coords
+	# need not be inbounds.
+
+	px1 = max (1.0,		 W_XS(w1))
+	px2 = min (real (IM_LEN(im,1)), W_XE(w1))
+	py1 = max (1.0,		 W_YS(w1))
+	py2 = min (real (IM_LEN(im,2)), W_YE(w1))
+
+	# Now compute the coordinates of the image section to be written in
+	# device pixel units.  This section must lie within or on the device
+	# window.
+	# This computation for I2S will give 257, which does differ by one
+	# for the Y center (due to inversion in I2S).  This should not matter,
+	# but if it does, this comment will change!
+
+	pxsize = px2 - px1
+	pysize = py2 - py1
+	wxcenter = (W_XE(w0) + W_XS(w0)) / 2.0 * cv_xres + 1
+	wycenter = (W_YE(w0) + W_YS(w0)) / 2.0 * cv_yres + 1
+
+	wx1 = max (1, int (wxcenter - (pxsize / 2.0 * xmag)))
+	wx2 = max (wx1, min (cv_xres, int (wx1 + (pxsize * xmag))))
+	wy1 = max (1, int (wycenter - (pysize / 2.0 * ymag)))
+	wy2 = max (wy1, min (cv_yres, int (wy1 + (pysize * ymag))))
+
+	# Display the image data, ignoring zero filling at the boundaries.
+
+	call ds_map_image (im, px1,px2,py1,py2, wx1,wx2,wy1,wy2,
+	    W_ZS(w1), W_ZE(w1), W_ZT(w1), W_UPTR(w1))
+
+	# Zero the border of the window if the frame has not been erased,
+	# and if the displayed section does not occupy the full window.
+
+	if (border_erase)
+	    call ds_erase_border (im, wdes, wx1,wx2,wy1,wy2)
+end
+
+
+# DS_MAP_IMAGE -- Map an image section from the input image to a section
+# (window) of the output image (the display device).  All spatial scaling is 
+# handled by the "scaled input" package, i.e., SIGL2[SR].  Our task is to
+# get lines from the scaled input image, transform the greyscale if necessary,
+# and write the lines to the output device.
+
+procedure ds_map_image (im, px1,px2,py1,py2, wx1,wx2,wy1,wy2, z1,z2,zt, uptr)
+
+pointer	im			# input image
+real	px1,px2,py1,py2		# input section
+int	wx1,wx2,wy1,wy2		# output section
+real	z1,z2			# range of input greylevels to be mapped.
+int	zt			# log or linear greylevel transformation
+pointer uptr			# pointer to user transformation table
+
+bool	unitary_greyscale_transformation
+short	lut1, lut2, z1_s, z2_s, dz1_s, dz2_s
+real	dz1, dz2
+int	wy, nx, ny, xblk, yblk
+pointer	in, out, si
+pointer	sigl2s(), sigl2r(), sigl2_setup()
+errchk	sigl2s, sigl2r, sigl2_setup
+real	xs, xe, y
+pointer	sp, outr
+bool    fp_equalr()
+real	if_elogr()
+extern	if_elogr
+
+include "cv.com"
+
+begin
+	call smark (sp)
+
+	# Set up for scaled image input.
+
+	nx = wx2 - wx1 + 1
+	ny = wy2 - wy1 + 1
+	xblk = INDEFI
+	yblk = INDEFI
+	si = sigl2_setup (im, px1,px2,nx,xblk, py1,py2,ny,yblk)
+
+	# Output array, and limiting x values in NDC
+
+	call salloc (out, nx, TY_SHORT)
+	xs = real(wx1 - 1) * cv_xcon / GKI_MAXNDC
+	# Don't subtract 1 from wx2 as we want it to be first one not filled
+	xe = real(wx2) * cv_xcon / GKI_MAXNDC
+	if ( xe > 1.0)
+	    xe = 1.0
+
+	# The device ZMIN and ZMAX parameters define the acceptable range
+	# of greyscale values for the output device (e.g., 0-255 for most 8-bit
+	# display devices).  For the general display, we use 0 and the
+	# device "z" resolution.  Values Z1 and Z2 are mapped linearly or
+	# logarithmically into these.
+
+	dz1 = 0
+	dz2 = cv_zres-1
+
+	# If the user specified the transfer function, see that the
+	# intensity and greyscale values are in range.
+
+	if (zt == W_USER) {
+	    call alims (Mems[uptr], SZ_BUF, lut1, lut2)
+	    dz1_s = short (dz1)
+	    dz2_s = short (dz2)
+	    if (lut2 < dz1_s || lut1 > dz2_s)
+		call eprintf ("User specified greyscales out of range\n")
+	    if (z2 < IM_MIN(im) || z1 > IM_MAX(im))
+		call eprintf ("User specified intensities out of range\n")
+	}
+
+	# Type short pixels are treated as a special case to minimize vector
+	# operations for such images (which are common).  If the image pixels
+	# are either short or real then only the ALTR (greyscale transformation)
+	# vector operation is required.  The ALTR operator linearly maps
+	# greylevels in the range Z1:Z2 to DZ1:DZ2, and does a floor ceiling
+	# of DZ1:DZ2 on all pixels outside the range.  If unity mapping is
+	# employed the data is simply copied, i.e., floor ceiling constraints
+	# are not applied.  This is very fast and will produce a contoured
+	# image on the display which will be adequate for some applications.
+
+	if (zt == W_UNITARY)
+	    unitary_greyscale_transformation = true
+	else
+	    unitary_greyscale_transformation =
+		(fp_equalr (dz1,z1) && fp_equalr (dz2,z2)) || fp_equalr (z1,z2)
+
+	if (IM_PIXTYPE(im) == TY_SHORT && zt != W_LOG) {
+
+	    # Set dz1_s and dz2_s depending on transformation
+	    if (zt != W_USER) {
+	        dz1_s = short (dz1)
+	        dz2_s = short (dz2)
+	    } else {
+	        dz1_s = short (STARTPT)
+	        dz2_s = short (ENDPT)
+	    }
+	    z1_s = short (z1)
+	    z2_s = short (z2)
+
+	    for (wy=wy1;  wy <= wy2;  wy=wy+1) {
+		in  = sigl2s (si, wy - wy1 + 1)
+		y = real(wy-1) * cv_ycon / GKI_MAXNDC
+		if (unitary_greyscale_transformation)
+		    call gpcell (cv_gp, Mems[in], nx, 1, xs, y, xe, y) 
+		else if (zt == W_USER) {
+		    call amaps (Mems[in], Mems[out], nx, z1_s,z2_s, dz1_s,dz2_s)
+		    call aluts (Mems[out], Mems[out], nx, Mems[uptr])
+		    call gpcell (cv_gp, Mems[out], nx, 1, xs, y, xe, y) 
+		} else {
+		    call amaps (Mems[in], Mems[out], nx, z1_s,z2_s, dz1_s,dz2_s)
+		    call gpcell (cv_gp, Mems[out], nx, 1, xs, y, xe, y) 
+		}
+	    }
+	} else {
+	    call salloc (outr, nx, TY_REAL)
+	    for (wy=wy1;  wy <= wy2;  wy=wy+1) {
+		in  = sigl2r (si, wy - wy1 + 1)
+		y = real(wy - 1) * cv_ycon / GKI_MAXNDC
+
+		if (zt == W_LOG) {
+		    call amapr (Memr[in], Memr[outr], nx,
+			z1, z2, 1.0, 10.0 ** MAXLOG)
+		    call alogr (Memr[outr], Memr[outr], nx, if_elogr)
+		    call amapr (Memr[outr], Memr[outr], nx,
+			1.0, real(MAXLOG), dz1, dz2)
+		    call achtrs (Memr[outr], Mems[out], nx)
+		} else if (unitary_greyscale_transformation) {
+		    call achtrs (Memr[in], Mems[out], nx)
+		} else if (zt == W_USER) {
+		    call amapr (Memr[in], Memr[outr], nx, z1,z2, STARTPT,ENDPT) 
+		    call achtrs (Memr[outr], Mems[out], nx)
+		    call aluts (Mems[out], Mems[out], nx, Mems[uptr])
+		} else {
+		    call amapr (Memr[in], Memr[outr], nx, z1, z2, dz1, dz2)
+		    call achtrs (Memr[outr], Mems[out], nx)
+		}
+		call gpcell (cv_gp, Mems[out], nx, 1, xs, y, xe, y)
+	    }
+	}
+
+	call sfree (sp)
+	call sigl2_free (si)
+end
+
+
+# DS_ERASE_BORDER -- Zero the border of the window if the frame has not been
+# erased, and if the displayed section does not occupy the full window.
+# It would be more efficient to do this while writing the greyscale data to
+# the output image, but that would complicate the display procedures and frames
+# are commonly erased before displaying an image.
+
+procedure ds_erase_border (im, wdes, wx1,wx2,wy1,wy2)
+
+pointer	im			# input image
+pointer	wdes			# window descriptor
+int	wx1,wx2,wy1,wy2		# section of display window filled by image data
+
+int	dx1,dx2,dy1,dy2		# coords of full display window in device pixels
+int	j, n, n1
+pointer	w0
+pointer	sp, zero
+real	xls, xle, xrs, xre, y
+
+include	"cv.com"
+
+begin
+	call smark (sp)
+	call salloc (zero, cv_xres, TY_SHORT)
+	call aclrs (Mems[zero], cv_xres)
+
+	# Compute device pixel coordinates of the full display window.
+	w0 = W_WC(wdes,0)
+	dx1 = W_XS(w0) * (cv_xres - 1) + 1
+	dx2 = W_XE(w0) * (cv_xres - 1) + 1
+	dy1 = W_YS(w0) * (cv_yres - 1) + 1
+	dy2 = W_YE(w0) * (cv_yres - 1) + 1
+
+	# Determine left and right (exclusive), start and end, x values in NDC
+	# for pixels not already filled.
+	# If, say, dx1 < wx1, we want to clear dx1 through wx1-1, which means
+	# that for gpcell, we want the (right) end points to be the first
+	# pixel not cleared.
+	xls = real(dx1 - 1) * cv_xcon / GKI_MAXNDC
+	xle = real(wx1) * cv_xcon / GKI_MAXNDC
+	if (xle > 1.0)
+	    xle = 1.0
+	xre = real(dx2 - 1) * cv_xcon / GKI_MAXNDC
+	xrs = real(wx2) * cv_xcon / GKI_MAXNDC
+	if (xre > 1.0)
+	    xre = 1.0
+
+	# Erase lower margin.
+	n = dx2 - dx1 + 1
+	for (j=dy1;  j < wy1;  j=j+1) {
+	    y = real(j-1) * cv_ycon / GKI_MAXNDC
+	    call gpcell (cv_gp, Mems[zero], n, 1, xls, y, xre, y)
+	}
+
+	# Erase left and right margins.  By doing the right margin of a line
+	# immediately after the left margin we have a high liklihood that the
+	# display line will still be in the FIO buffer.
+
+	n = wx1 - dx1
+	n1 = dx2 - wx2
+	for (j=wy1;  j <= wy2;  j=j+1) {
+	    y = real(j-1) * cv_ycon / GKI_MAXNDC
+	    if (dx1 < wx1)
+		call gpcell (cv_gp, Mems[zero], n, 1, xls, y, xle, y)
+	    if (wx2 < dx2)
+		call gpcell (cv_gp, Mems[zero], n1, 1, xrs, y, xre, y)
+	}
+
+	# Erase upper margin.
+	n = dx2 - dx1 + 1
+	for (j=wy2+1;  j <= dy2;  j=j+1) {
+	    y = real(j-1) * cv_ycon / GKI_MAXNDC
+	    call gpcell (cv_gp, Mems[zero], n, 1, xls, y, xre, y)
+	}
+
+	call sfree (sp)
+end
+
+
+# IF_ELOG -- The error function for log10. Note that MAX_EXPONENT is
+# currently an integer so it is converted to the appropriate data type
+# before being returned.
+
+real procedure if_elogr (x)
+
+real   x                               # the input pixel value
+
+begin
+        return (real(-MAX_EXPONENT))
+end
+