1 files changed, 117 insertions, 0 deletions
diff --git a/pkg/plot/crtpict/xformimage.x b/pkg/plot/crtpict/xformimage.x
new file mode 100644
index 00000000..0400cc91
--- /dev/null
+++ b/pkg/plot/crtpict/xformimage.x
@@ -0,0 +1,117 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<mach.h>
+include	<gset.h>
+include	"wdes.h"
+include	"crtpict.h"
+
+# CRT_TRANSFORM_IMAGE -- Map an image into the output device.  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.  The input image is accessed via IMIO; the 
+# transformed image is output to a greyscale device via GIO.    
+
+procedure crt_transform_image (gp, im, wdes, cl)
+
+pointer	gp			# graphics descriptor for output
+pointer	im			# input image
+pointer	wdes			# graphics window descriptor
+pointer	cl
+
+real	ndc_xs,ndc_xe,ndc_ys,ndc_ye	# NDC of output device window
+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)
+int	ndev_cols, ndev_rows, nx_output, ny_output
+int	ggeti()
+errchk	ggeti, gseti, gsview, gswind, draw_perimeter
+errchk	crt_map_image, crt_replicate_image
+
+begin
+	# Compute pointers to WCS 0 and 1.
+	w0 = W_WC(wdes,0)
+	w1 = W_WC(wdes,1)
+
+	call printf ("%s: %s\n")
+	    call pargstr (W_IMSECT(wdes))
+	    call pargstr (IM_TITLE(im))
+
+	ndev_cols = ggeti (gp, "xr")
+	ndev_rows = ggeti (gp, "yr")
+
+	# Compute X and Y magnification ratios required to map image into
+	# the device window.
+
+	xmag = ((W_XE(w0) - W_XS(w0)) * ndev_cols) / (W_XE(w1) - W_XS(w1))
+	ymag = ((W_YE(w0) - W_YS(w0)) * ndev_rows) / (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, but they are integral pixels.
+
+	px1 = max (1.0,		        real (int (W_XS(w1) + 0.5)))
+	px2 = min (real (IM_LEN(im,1)), real (int (W_XE(w1) + 0.5)))
+	py1 = max (1.0,		        real (int (W_YS(w1) + 0.5)))
+	py2 = min (real (IM_LEN(im,2)), real (int (W_YE(w1) + 0.5)))
+
+	# The extent of the output image in NDC coords
+	pxsize = ((px2 - px1) + 1.0) / ndev_cols
+	pysize = ((py2 - py1) + 1.0) / ndev_rows
+
+	# The NDC coordinates that map to the central image pixel output
+	wxcenter = (W_XE(w0) + W_XS(w0)) / 2.0
+	wycenter = (W_YE(w0) + W_YS(w0)) / 2.0
+
+	# Now compute the NDC coordinates of the output device that the
+	# image will occupy.  
+	ndc_xs = max (W_XS(w0), wxcenter - (pxsize / 2.0 * xmag))
+	ndc_xe = max (ndc_xs, min (W_XE(w0), ndc_xs + (pxsize * xmag)))
+	ndc_ys = max (W_YS(w0), wycenter - (pysize / 2.0 * ymag))
+	ndc_ye = max (ndc_ys, min (W_YE(w0), ndc_ys + (pysize * ymag)))
+
+	# To avoid possible truncation errors down the line, make sure
+	# the ndc coordinates passed to the output procedures represent
+	# integer pixels.
+	ndc_xs = real (int (ndc_xs * ndev_cols)) / ndev_cols
+	ndc_xe = real (int (ndc_xe * ndev_cols)) / ndev_cols
+	ndc_ys = real (int (ndc_ys * ndev_rows)) / ndev_rows
+	ndc_ye = real (int (ndc_ye * ndev_rows)) / ndev_rows
+
+	# Output the image data in WCS 0.  The number of image pixels that
+	# will be put out across the device is calculated first.  
+
+	call gseti (gp, G_WCS, 0)
+	if (REPLICATE(cl) == YES) {
+	    nx_output = (int(px2) - int(px1)) + 1
+	    ny_output = (int(py2) - int(py1)) + 1
+	} else {
+	    # Image pixels will be scaled to number of device pixels 
+	    nx_output = ((ndc_xe * ndev_cols) - (ndc_xs * ndev_cols)) + 1
+	    ny_output = ((ndc_ye * ndev_rows) - (ndc_ys * ndev_rows)) + 1
+	    nx_output = nx_output / X_BA(cl)
+	    ny_output = ny_output / Y_BA(cl)
+	}
+
+	# Tweak the ndc coordinates to insure integer replication factors.
+	# This may change the ndc coordinates.
+	call crt_tweak_ndc (nx_output, ndc_xs, ndc_xe, ndev_cols)
+	call crt_tweak_ndc (ny_output, ndc_ys, ndc_ye, ndev_rows)
+
+	# Call routine that actually puts out the pixels row by row
+	call crt_map_image (im, gp, px1,px2,py1,py2, ndc_xs,ndc_xe,
+	       ndc_ys,ndc_ye, nx_output, ny_output, W_ZS(w1),W_ZE(w1),W_ZT(w1),
+	           cl)
+
+	# Change to pixel coordinates and draw perimeter axes if requested
+	call gseti  (gp, G_WCS, 2)
+	call gsview (gp, ndc_xs, ndc_xe, ndc_ys, ndc_ye)
+	call gswind (gp, px1, px2, py1, py2)
+
+	if (PERIM(cl) == YES)
+	    call draw_perimeter (gp)
+end