Initial commit

8 files changed, 1402 insertions, 0 deletions

diff --git a/pkg/proto/vol/src/i2sun/cnvimage.x b/pkg/proto/vol/src/i2sun/cnvimage.x
new file mode 100644
index 00000000..59bd4655
--- /dev/null
+++ b/pkg/proto/vol/src/i2sun/cnvimage.x
@@ -0,0 +1,142 @@
+include <imhdr.h>
+include <mach.h>
+include "i2sun.h"
+
+
+# CNV_IMAGE -- Read each line of the input image, apply ztransform, convert
+# to rasterfile form, and write to rasterfile.
+
+procedure cnv_image (im, slice, tr, uptr, rfd)
+pointer	im		# input image
+int	slice		# current slice if n>2 image
+pointer	tr		# spatial & greyscale transform structure
+pointer	uptr		# pointer to user-specified transfer lut
+pointer rfd		# output rasterfile
+
+real	z1, z2, rz1, rz2
+int	ztrans, row, xblk, yblk, ocols, olines
+real	px1, px2, py1, py2	# image coords in fractional image pixels
+pointer	sp, si, bufptr, out, rtemp, packed
+short	z1_s, z2_s, rz1_s, rz2_s
+bool	unitary_greyscale_transformation
+
+bool	fp_equalr()
+pointer	siglns(), siglnr(), sigln_setup()
+errchk	siglns(), siglnr(), sigln_setup()
+
+begin
+	# Set up for scaled image input.
+	px1 = 1
+	px2 = IM_LEN(im,COL)
+	py1 = 1
+	py2 = IM_LEN(im,LINE)
+	ocols = TR_XE(tr) - TR_XS(tr) + 1
+	olines = TR_YE(tr) - TR_YS(tr) + 1
+
+	# Round odd-numbered numbers of columns up due to restrictions on
+	# IRAF binary byte i/o (number of bytes of image data must match
+	# that specified in rasterfile header).
+	if (mod (ocols, 2) == 1)
+	    ocols = ocols + 1
+
+	xblk = INDEFI
+	yblk = INDEFI
+	si = sigln_setup (im, px1,px2,ocols,xblk, py1,py2,olines,yblk, 
+	    TR_ORDER(tr))
+
+	# Type short pixels are treated as a special case to minimize vector
+	# operations for such images (which are common).  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.
+
+	z1 = TR_Z1(tr)
+	z2 = TR_Z2(tr)
+	ztrans = TR_ZTRANS(tr)
+	rz1 = COLORSTART
+	rz2 = COLOREND
+	if (ztrans == Z_UNITARY) {
+	    unitary_greyscale_transformation = true
+	} else if (ztrans == Z_LINEAR) {
+	    unitary_greyscale_transformation =
+		((fp_equalr(z1,rz1) && fp_equalr(z2,rz2)) || fp_equalr(z1,z2))
+	} else
+	    unitary_greyscale_transformation = false
+
+	if (IM_PIXTYPE(im) == TY_SHORT && ztrans != Z_LOG) {
+
+	    call smark (sp)
+	    call salloc (out, ocols, TY_SHORT)
+	    call salloc (packed, ocols, TY_CHAR)
+	    z1_s = z1;  z2_s = z2
+
+	    for (row=olines;  row >= 1;  row=row-1) {
+		bufptr = siglns (si, row, TR_SLICEAXIS(tr), slice)
+
+		if (unitary_greyscale_transformation) {
+		    call amovs (Mems[bufptr], Mems[out], ocols)
+		} else if (ztrans == Z_USER) {
+		    rz1_s = U_Z1; rz2_s = U_Z2
+		    call amaps (Mems[bufptr], Mems[out], ocols, z1_s, z2_s,
+			rz1_s, rz2_s)
+		    call aluts (Mems[out], Mems[out], ocols, Mems[uptr])
+		} else {
+		    rz1_s = rz1; rz2_s = rz2
+		    call amaps (Mems[bufptr], Mems[out], ocols, z1_s, z2_s,
+			rz1_s, rz2_s)
+		}
+
+		# Pack to unsigned byte and write to file.
+		call achtsc (Mems[out], Memc[packed], ocols)
+		call chrpak (Memc[packed], 1, Memc[packed], 1, ocols)
+		call write (rfd, Memc[packed], ocols / SZB_CHAR)
+	    }
+	    call sfree (sp)
+
+	} else if (ztrans == Z_USER) {
+	    call smark (sp)
+	    call salloc (rtemp, ocols, TY_REAL)
+	    call salloc (out, ocols, TY_SHORT)
+	    call salloc (packed, ocols, TY_CHAR)
+
+	    for (row=olines;  row >= 1;  row=row-1) {
+		bufptr = siglnr (si, row, TR_SLICEAXIS(tr), slice)
+		call amapr (Memr[bufptr], Memr[rtemp], ocols, z1, z2, 
+		    real(U_Z1), real(U_Z2))
+		call achtrs (Memr[rtemp], Mems[out], ocols)
+		call aluts (Mems[out], Mems[out], ocols, Mems[uptr])
+		call achtsc (Mems[out], Memc[packed], ocols)
+		call chrpak (Memc[packed], 1, Memc[packed], 1, ocols)
+		call write (rfd, Memc[packed], ocols / SZB_CHAR)
+	    }
+	    call sfree (sp)
+
+	} else {
+	    call smark (sp)
+	    call salloc (rtemp, ocols, TY_REAL)
+	    call salloc (packed, ocols, TY_CHAR)
+
+	    for (row=olines;  row >= 1;  row=row-1) {
+		bufptr = siglnr (si, row, TR_SLICEAXIS(tr), slice)
+
+		if (unitary_greyscale_transformation)  {
+		    call amovr (Memr[bufptr], Memr[rtemp], ocols)
+		} else if (ztrans == Z_LOG) {
+		    call amapr (Memr[bufptr], Memr[rtemp], ocols,
+			z1, z2, 1.0, 10.0 ** MAXLOG)
+		    call alogr (Memr[rtemp], Memr[rtemp], ocols)
+		    call amapr (Memr[rtemp], Memr[rtemp], ocols,
+			1.0, real(MAXLOG), rz1, rz2)
+		} else
+		    call amapr (Memr[bufptr], Memr[rtemp], ocols, z1, z2,
+			rz1, rz2)
+		call achtrc (Memr[rtemp], Memc[packed], ocols)
+		call chrpak (Memc[packed], 1, Memc[packed], 1, ocols)
+		call write (rfd, Memc[packed], ocols / SZB_CHAR)
+	    }
+	    call sfree (sp)
+	}
+
+	call sigln_free (si)
+end
+
diff --git a/pkg/proto/vol/src/i2sun/i2sun.h b/pkg/proto/vol/src/i2sun/i2sun.h
new file mode 100644
index 00000000..73f2ea3f
--- /dev/null
+++ b/pkg/proto/vol/src/i2sun/i2sun.h
@@ -0,0 +1,46 @@
+# I2SUNRAS.H -- Include file for IRAF to Sun rasterfile program i2sunras.
+
+define	COL		1
+define	LINE		2
+define	BAND		3
+define	Z_LINEAR	1		# linear ztransform
+define	Z_LOG		2		# log ztransform
+define	Z_UNITARY	3		# no ztransform
+define	Z_USER		4		# user-specified transform
+define	U_MAXPTS	4096		# max user-specified lut pairs (DISPLAY)
+define	U_Z1		0		# base user-specified transfer val
+define	U_Z2		4095		# upper user-specified transfer val
+define	MAXLOG		3		# if log, map to 1:10**MAXLOG b4 log10
+define	DSP_MIN		0		# minimum display pixel value
+define	DSP_MAX		255		# maximum display pixel value
+define	RAS_HDR_INTS	8		# SunOS4.0 and earlier
+define	RMT_NONE	0		# SunOS4.0 and earlier
+define	RMT_EQUAL_RGB	1		# SunOS4.0 and earlier
+define	RMT_STANDARD	1		# SunOS4.0 and earlier
+define	RAS_MAGIC	1504078485	# SunOS4.0 and earlier
+define	NGREY		256		# SunOS4.0 and earlier, 8bit fb
+define	COLORSTART	1		# IMTOOL
+define	COLOREND	200		# IMTOOL
+define	COLORRANGE	200		# IMTOOL
+define	WHITE		(NGREY-1)	# IMTOOL
+define	BLACK		0		# IMTOOL
+define	NBITS_FB	8
+define	wrapup_		91
+
+# Spatial and greyscale transformation structure.
+define	LEN_TR		20
+define	TR_ZTRANS	Memi[$1]	# Greyscale transformation.
+define	TR_Z1		Memr[P2R($1+1)]	# Minimum data z-value
+define	TR_Z2		Memr[P2R($1+2)]	# Maximum data z-value
+define	TR_XSIZE	Memi[$1+3]	# Output rasterfile size in x
+define	TR_YSIZE	Memi[$1+4]	# Output rasterfile size in y
+define	TR_XMAG		Memr[P2R($1+5)]	# Magnification factor in x
+define	TR_YMAG		Memr[P2R($1+6)]	# Magnification factor in y
+define	TR_ORDER	Memi[$1+7]	# Interpolation order
+define	TR_XS		Memi[$1+8]	# Starting output x pixel index
+define	TR_XE		Memi[$1+9]	# Ending output x pixel index
+define	TR_YS		Memi[$1+10]	# Starting output y pixel index
+define	TR_YE		Memi[$1+11]	# Ending output y pixel index
+define	TR_SLICEAXIS	Memi[$1+12]	# Slice or frame axis when ndim>2 
+define	TR_SWAPBYTES	Memb[$1+13]	# Swap output bytes?
+#					# Reserved space
diff --git a/pkg/proto/vol/src/i2sun/mkpkg b/pkg/proto/vol/src/i2sun/mkpkg
new file mode 100644
index 00000000..b1a8c4f4
--- /dev/null
+++ b/pkg/proto/vol/src/i2sun/mkpkg
@@ -0,0 +1,27 @@
+# Library for the I2SUN task.
+
+$checkout libpkg.a ../
+$update   libpkg.a
+$checkin  libpkg.a ../
+$exit
+
+libpkg.a:
+	t_i2sun.x	<imhdr.h> <ctype.h> i2sun.h
+	trulut.x	<error.h> <ctype.h> i2sun.h
+	trsetup.x	<imhdr.h> i2sun.h
+	cnvimage.x	<imhdr.h> i2sun.h
+	sigln.x		<error.h> <imhdr.h>
+	;
+
+dbx:
+	$set	XFLAGS = "-c -g -F -q"
+	$set	LFLAGS = "-g -q"
+	$omake	x_i2sun.x
+	$omake	t_i2sun.x
+	$omake	trulut.x
+	$omake	trsetup.x
+	$omake	cnvimage.x
+	$omake	sigln.x
+	$link	x_i2sun.o t_i2sun.o trulut.o trsetup.o cnvimage.o sigln.o \
+		-o xx_i2sun.e
+	;
diff --git a/pkg/proto/vol/src/i2sun/sigln.x b/pkg/proto/vol/src/i2sun/sigln.x
new file mode 100644
index 00000000..9d763b3f
--- /dev/null
+++ b/pkg/proto/vol/src/i2sun/sigln.x
@@ -0,0 +1,783 @@
+include	<imhdr.h>
+include	<error.h>
+
+.help sigl2, sigl2_setup
+.nf ___________________________________________________________________________
+SIGLN -- Get a line from a spatially scaled image of any dimensionality.
+This procedure works like the regular IMIO get line procedure, but rescales
+the input image in 1 or two axes upon input (for a resulting 2d output image).
+If the magnification
+ratio required is greater than 0 and less than 2 then linear interpolation is
+used to resample the image.  If the magnification ratio is greater than or
+equal to 2 then the image is block averaged by the smallest factor which
+reduces the magnification to the range 0-2 and then interpolated back up to
+the desired size.  In some cases this will smooth the data slightly, but the
+operation is efficient and avoids aliasing effects.
+
+	si =	sigln_setup (im, x1,x2,nx,xblk, y1,y2,ny,yblk, order)
+		sigln_free (si)
+	ptr =	sigln[sr] (si, linenumber)
+
+SIGLN_SETUP must be called to set up the transformations after mapping the
+image and before performing any scaled i/o to the image.  SIGLN_FREE must be
+called when finished to return buffer space.
+.endhelp ______________________________________________________________________
+
+# Scaled image descriptor for 2-dim images
+
+define	SI_LEN		16
+define	SI_MAXDIM	2		# 2 dimensions of spatial scaling
+define	SI_NBUFS	3		# nbuffers used by SIGLN
+
+define	SI_IM		Memi[$1]	# pointer to input image header
+define	SI_GRID		Memi[$1+1+$2-1]	# pointer to array of X coords
+define	SI_NPIX		Memi[$1+3+$2-1]	# number of X coords
+define	SI_BAVG		Memi[$1+5+$2-1]	# X block averaging factor
+define	SI_INTERP	Memi[$1+7+$2-1]	# interpolate X axis
+define	SI_BUF		Memi[$1+9+$2-1]	# line buffers
+define	SI_ORDER	Memi[$1+12]	# interpolator order, 0 or 1
+define	SI_TYBUF	Memi[$1+13]	# buffer type
+define	SI_XOFF		Memi[$1+14]	# offset in input image to first X
+define	SI_INIT		Memi[$1+15]	# YES until first i/o is done
+
+define	OUTBUF		SI_BUF($1,3)
+
+define	SI_TOL		(1E-5)		# close to a pixel
+define	INTVAL		(abs ($1 - nint($1)) < SI_TOL)
+define	SWAPI		{tempi=$2;$2=$1;$1=tempi}
+define	SWAPP		{tempp=$2;$2=$1;$1=tempp}
+define	NOTSET		(-9999)
+
+# SIGLN_SETUP -- Set up the spatial transformation for SIGLN[SR].  Compute
+# the block averaging factors (1 if no block averaging is required) and
+# the sampling grid points, i.e., pixel coordinates of the output pixels in
+# the input image.
+
+pointer procedure sigln_setup (im, px1,px2,nx,xblk, py1,py2,ny,yblk, order)
+
+pointer	im			# the input image
+real	px1, px2		# range in X to be sampled on an even grid
+int	nx			# number of output pixels in X
+int	xblk			# blocking factor in x
+real	py1, py2		# range in Y to be sampled on an even grid
+int	ny			# number of output pixels in Y
+int	yblk			# blocking factor in y
+int	order			# interpolator order (0=replicate, 1=linear)
+
+int	npix, noldpix, nbavpix, i, j
+int	npts[SI_MAXDIM]		# number of output points for axis
+int	blksize[SI_MAXDIM]	# block averaging factor (npix per block)
+real	tau[SI_MAXDIM]		# tau = p(i+1) - p(i) in fractional pixels
+real	p1[SI_MAXDIM]		# starting pixel coords in each axis
+real	p2[SI_MAXDIM]		# ending pixel coords in each axis
+real	scalar, start
+pointer	si, gp
+
+begin
+	iferr (call calloc (si, SI_LEN, TY_STRUCT))
+	    call erract (EA_FATAL)
+
+	SI_IM(si) = im
+	SI_NPIX(si,1) = nx
+	SI_NPIX(si,2) = ny
+	SI_ORDER(si) = order
+	SI_INIT(si) = YES
+
+	p1[1] = px1			# X = index 1
+	p2[1] = px2
+	npts[1] = nx
+	blksize[1] = xblk
+
+	p1[2] = py1			# Y = index 2
+	p2[2] = py2
+	npts[2] = ny
+	blksize[2] = yblk
+
+	# Compute block averaging factors if not defined.
+	# If there is only one pixel then the block average is the average
+	# between the first and last point.
+
+	do i = 1, SI_MAXDIM {
+	    if ((blksize[i] >= 1) && (blksize[i] != INDEFI)) {
+	        if (npts[i] == 1)
+		    tau[i] = 0.
+	        else
+	            tau[i] = (p2[i] - p1[i]) / (npts[i] - 1)
+	    } else {
+		if (npts[i] == 1) {
+		    tau[i] = 0.
+		    blksize[i] = int (p2[i] - p1[i] + 1)
+	        } else {
+	            tau[i] = (p2[i] - p1[i]) / (npts[i] - 1)
+	    	    if (tau[i] >= 2.0) {
+
+			# If nx or ny is not an integral multiple of the block
+			# averaging factor, noldpix is the next larger number
+			# which is an integral multiple.  When the image is
+			# block averaged pixels will be replicated as necessary
+			# to fill the last block out to this size.  
+
+			blksize[i] = int (tau[i])
+			npix = p2[i] - p1[i] + 1
+			noldpix = (npix+blksize[i]-1) / blksize[i] * blksize[i]
+			nbavpix = noldpix / blksize[i]
+			scalar = real (nbavpix - 1) / real (noldpix - 1)
+			p1[i] = (p1[i] - 1.0) * scalar + 1.0
+			p2[i] = (p2[i] - 1.0) * scalar + 1.0
+			tau[i] = (p2[i] - p1[i]) / (npts[i] - 1)
+		    } else
+			blksize[i] = 1
+		}
+	    }
+	}
+
+	SI_BAVG(si,1) = blksize[1]
+	SI_BAVG(si,2) = blksize[2]
+
+	if (IS_INDEFI (xblk))
+	    xblk = blksize[1]
+	if (IS_INDEFI (yblk))
+	    yblk = blksize[2]
+
+	# Allocate and initialize the grid arrays, specifying the X and Y
+	# coordinates of each pixel in the output image, in units of pixels
+	# in the input (possibly block averaged) image.
+
+	do i = 1, SI_MAXDIM {
+	    # The X coordinate is special.  We do not want to read entire
+	    # input image lines if only a range of input X values are needed.
+	    # Since the X grid vector passed to ALUI (the interpolator) must
+	    # contain explicit offsets into the vector being interpolated,
+	    # we must generate interpolator grid points starting near 1.0.
+	    # The X origin, used to read the block averaged input line, is
+	    # given by XOFF.
+
+	    if (i == 1) {
+		SI_XOFF(si) = int (p1[i])
+		start = p1[1] - int (p1[i]) + 1.0
+	    } else
+	    	start = p1[i]
+
+	    # Do the axes need to be interpolated?
+	    if (INTVAL(start) && INTVAL(tau[i]))
+		SI_INTERP(si,i) = NO
+	    else
+		SI_INTERP(si,i) = YES
+
+	    # Allocate grid buffer and set the grid points.
+	    iferr (call malloc (gp, npts[i], TY_REAL))
+		call erract (EA_FATAL)
+	    SI_GRID(si,i) = gp
+	    if (SI_ORDER(si) <= 0) {
+		do j = 0, npts[i]-1
+		    Memr[gp+j] = int (start + (j * tau[i]) + 0.5)
+	    } else {
+		do j = 0, npts[i]-1
+		    Memr[gp+j] = start + (j * tau[i])
+	    }
+	}
+
+	return (si)
+end
+
+
+# SIGLN_FREE -- Free storage associated with an image opened for scaled
+# input.  This does not close and unmap the image.
+
+procedure sigln_free (si)
+
+pointer	si
+int	i
+
+begin
+	# Free SIGLN buffers.
+	do i = 1, SI_NBUFS
+	    if (SI_BUF(si,i) != NULL)
+		call mfree (SI_BUF(si,i), SI_TYBUF(si))
+
+	# Free GRID buffers.
+	do i = 1, SI_MAXDIM
+	    if (SI_GRID(si,i) != NULL)
+		call mfree (SI_GRID(si,i), TY_REAL)
+
+	call mfree (si, TY_STRUCT)
+end
+
+
+# SIGLNS -- Get a line of type short from a scaled image.  Block averaging is
+# done by a subprocedure; this procedure gets a line from a possibly block
+# averaged image and if necessary interpolates it to the grid points of the
+# output line.
+
+pointer procedure siglns (si, lineno, slice_axis, slice)
+
+pointer	si		# pointer to SI descriptor
+int	lineno
+int	slice_axis	# axis from which to slice section for ndim>2 images
+int	slice		# current slice index
+
+pointer	rawline, tempp, gp
+int	i, buf_y[2], new_y[2], tempi, curbuf, altbuf
+int	npix, nblks_y, ybavg, x1, x2
+real	x, y, weight_1, weight_2
+pointer	si_blkavgs()
+errchk	si_blkavgs
+
+begin
+	npix = SI_NPIX(si,1)
+
+	# Determine the range of X (in pixels on the block averaged input image)
+	# required for the interpolator.
+
+	gp = SI_GRID(si,1)
+	x1 = SI_XOFF(si)
+	x = Memr[gp+npix-1]
+	x2 = x1 + int(x)
+	if (INTVAL(x))
+	    x2 = x2 - 1
+	x2 = max (x1 + 1, x2)
+
+	gp = SI_GRID(si,2)
+	y = Memr[gp+lineno-1]
+
+	# The following is an optimization provided for the case when it is
+	# not necessary to interpolate in either X or Y.  Block averaging is
+	# permitted.
+
+	if (SI_INTERP(si,1) == NO && SI_INTERP(si,2) == NO)
+	    return (si_blkavgs (SI_IM(si), x1, x2, int(y),
+		SI_BAVG(si,1), SI_BAVG(si,2), slice_axis, slice))
+
+	# If we are interpolating in Y two buffers are required, one for each
+	# of the two input image lines required to interpolate in Y.  The lines
+	# stored in these buffers are interpolated in X to the output grid but
+	# not in Y.  Both buffers are not required if we are not interpolating
+	# in Y, but we use them anyhow to simplify the code.
+
+	if (SI_INIT(si) == YES) {
+	    do i = 1, 2 {
+		if (SI_BUF(si,i) != NULL)
+		    call mfree (SI_BUF(si,i), SI_TYBUF(si))
+		call malloc (SI_BUF(si,i), npix, TY_SHORT)
+		SI_TYBUF(si) = TY_SHORT
+		buf_y[i] = NOTSET
+	    }
+	    if (OUTBUF(si) != NULL)
+		call mfree (OUTBUF(si), SI_TYBUF(si))
+	    call malloc (OUTBUF(si), npix, TY_SHORT)
+	    SI_INIT(si) = NO
+	}
+
+	# If the Y value of the new line is not in range of the contents of the
+	# current line buffers, refill one or both buffers.  To refill we must
+	# read a (possibly block averaged) input line and interpolate it onto
+	# the X grid.  The X and Y values herein are in the coordinate system
+	# of the (possibly block averaged) input image.
+
+	new_y[1] = int(y)
+	new_y[2] = int(y) + 1
+
+	# Get the pair of lines whose integral Y values form an interval
+	# containing the fractional Y value of the output line.  Sometimes the
+	# desired line will happen to be in the other buffer already, in which
+	# case we just have to swap buffers.  Often the new line will be the
+	# current line, in which case nothing is done.  This latter case occurs
+	# frequently when the magnification ratio is large.
+
+	curbuf = 1
+	altbuf = 2
+
+	do i = 1, 2 {
+	    if (new_y[i] == buf_y[i]) {
+		;
+	    } else if (new_y[i] == buf_y[altbuf]) {
+		SWAPP (SI_BUF(si,1), SI_BUF(si,2))
+		SWAPI (buf_y[1], buf_y[2])
+
+	    } else {
+		# Get line and interpolate onto output grid.  If interpolation
+		# is not required merely copy data out.  This code is set up
+		# to always use two buffers; in effect, there is one buffer of
+		# look ahead, even when Y[i] is integral.  This means that we
+		# will go out of bounds by one line at the top of the image.
+		# This is handled by copying the last line.
+
+		ybavg = SI_BAVG(si,2)
+		nblks_y = (IM_LEN (SI_IM(si), 2) + ybavg-1) / ybavg
+		if (new_y[i] <= nblks_y)
+		    rawline = si_blkavgs (SI_IM(si), x1, x2, new_y[i],
+			SI_BAVG(si,1), SI_BAVG(si,2), slice_axis, slice)
+
+		if (SI_INTERP(si,1) == NO) {
+		    call amovs (Mems[rawline], Mems[SI_BUF(si,i)], npix)
+		} else if (SI_ORDER(si) <= 0) {
+		    call si_samples (Mems[rawline], Mems[SI_BUF(si,i)],
+			Memr[SI_GRID(si,1)], npix)
+		} else {
+		    call aluis (Mems[rawline], Mems[SI_BUF(si,i)],
+			Memr[SI_GRID(si,1)], npix)
+		}
+
+		buf_y[i] = new_y[i]
+	    }
+
+	    SWAPI (altbuf, curbuf)
+	}
+
+	# We now have two line buffers straddling the output Y value,
+	# interpolated to the X grid of the output line.  To complete the
+	# bilinear interpolation operation we take a weighted sum of the two
+	# lines.  If the range from buf_y[1] to buf_y[2] is repeatedly
+	# interpolated in Y no additional i/o occurs and the linear
+	# interpolation operation (ALUI) does not have to be repeated (only the
+	# weighted sum is required).  If the distance of Y from one of the
+	# buffers is zero then we do not even have to take a weighted sum.
+	# This is not unusual because we may be called with a magnification
+	# of 1.0 in Y.
+
+	weight_1 = 1.0 - (y - buf_y[1])
+	weight_2 = 1.0 - weight_1
+
+	if (weight_1 < SI_TOL)
+	    return (SI_BUF(si,2))
+	else if (weight_2 < SI_TOL || SI_ORDER(si) <= 0) 
+	    return (SI_BUF(si,1))
+	else {
+	    call awsus (Mems[SI_BUF(si,1)], Mems[SI_BUF(si,2)],
+		Mems[OUTBUF(si)], npix, weight_1, weight_2)
+	    return (OUTBUF(si))
+	}
+end
+
+
+# SI_BLKAVGS -- Get a line from a block averaged image of type short.
+# For example, block averaging by a factor of 2 means that pixels 1 and 2
+# are averaged to produce the first output pixel, 3 and 4 are averaged to
+# produce the second output pixel, and so on.  If the length of an axis
+# is not an integral multiple of the block size then the last pixel in the
+# last block will be replicated to fill out the block; the average is still
+# defined even if a block is not full.
+
+pointer procedure si_blkavgs (im, x1, x2, y, xbavg, ybavg, slice_axis, slice)
+
+pointer	im			# input image
+int	x1, x2			# range of x blocks to be read
+int	y			# y block to be read
+int	xbavg, ybavg		# X and Y block averaging factors
+int	slice_axis		# slice dimension if ndim>2 image
+int	slice			# slice if ndim>2 image
+
+short	temp_s
+int	nblks_x, nblks_y, ncols, nlines, xoff, i, j
+int	first_line, nlines_in_sum, npix, nfull_blks, count
+long	vs[IM_MAXDIM], ve[IM_MAXDIM]
+real	sum
+pointer	sp, a, b
+pointer	imgs2s(), imgs3s(), imggss()
+errchk	imgs2s, imgs3s, imggss
+
+begin
+	call smark (sp)
+
+	ncols  = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	xoff   = (x1 - 1) * xbavg + 1
+	npix   = min (ncols, xoff + (x2 - x1 + 1) * xbavg - 1)
+
+	if ((xbavg < 1) || (ybavg < 1))
+	    call error (1, "si_blkavg: illegal block size")
+	else if (x1 < 1 || x2 > ncols)
+	    call error (2, "si_blkavg: column index out of bounds")
+	else if ((xbavg == 1) && (ybavg == 1))
+	    if (IM_NDIM(im) == 2)
+		return (imgs2s (im, xoff, xoff + npix - 1, y, y))
+	    else if (IM_NDIM(im) == 3)
+		return (imgs3s (im, xoff, xoff + npix - 1, y, y, slice, slice))
+	    else {
+		call amovkl (long(1), vs, IM_MAXDIM)
+		call amovkl (long(1), ve, IM_MAXDIM)
+		vs[1] = xoff
+		ve[1] = xoff + npix - 1
+		vs[2] = y
+		ve[2] = y
+		vs[slice_axis] = slice
+		ve[slice_axis] = slice
+		return (imggss (im, vs, ve, 2))
+	    }
+
+	nblks_x = (npix   + xbavg-1) / xbavg
+	nblks_y = (nlines + ybavg-1) / ybavg
+
+	if (y < 1 || y > nblks_y)
+	    call error (2, "si_blkavg: block number out of range")
+
+	call salloc (b, nblks_x, TY_SHORT)
+
+	if (ybavg > 1) {
+	    call aclrs (Mems[b], nblks_x)
+	    nlines_in_sum = 0
+	}
+
+	# Read and accumulate all input lines in the block.
+	first_line = (y - 1) * ybavg + 1
+
+	do i = first_line, min (nlines, first_line + ybavg - 1) {
+	    # Get line from input image.
+	    if (IM_NDIM(im) == 2)
+		a = imgs2s (im, xoff, xoff + npix - 1, i, i)
+	    else if (IM_NDIM(im) == 3)
+		a = imgs3s (im, xoff, xoff + npix - 1, i, i, slice, slice)
+	    else {
+		call amovkl (long(1), vs, IM_MAXDIM)
+		call amovkl (long(1), ve, IM_MAXDIM)
+		vs[1] = xoff
+		ve[1] = xoff + npix - 1
+		vs[2] = i
+		ve[2] = i
+		vs[slice_axis] = slice
+		ve[slice_axis] = slice
+		return (imggss (im, vs, ve, 2))
+	    }
+
+	    # Block average line in X.
+	    if (xbavg > 1) {
+		# First block average only the full blocks.
+		nfull_blks = npix / xbavg
+		call abavs (Mems[a], Mems[a], nfull_blks, xbavg)
+
+		# Now average the final partial block, if any.
+		if (nfull_blks < nblks_x) {
+		    sum = 0.0
+		    count = 0
+		    do j = nfull_blks * xbavg + 1, npix {
+			sum = sum + Mems[a+j-1]
+			count = count + 1
+		    }
+		    Mems[a+nblks_x-1] = sum / count
+		}
+	    }
+
+	    # Add line into block sum.  Keep track of number of lines in sum
+	    # so that we can compute block average later.
+	    if (ybavg > 1) {
+		call aadds (Mems[a], Mems[b], Mems[b], nblks_x)
+		nlines_in_sum = nlines_in_sum + 1
+	    }
+	}
+
+	# Compute the block average in Y from the sum of all lines block
+	# averaged in X.  Overwrite buffer A, the buffer returned by IMIO.
+	# This is kosher because the block averaged line is never longer
+	# than an input line.
+
+	if (ybavg > 1) {
+	    temp_s = nlines_in_sum
+	    call adivks (Mems[b], temp_s, Mems[a], nblks_x)
+	}
+
+	call sfree (sp)
+	return (a)
+end
+
+
+# SI_SAMPLES -- Resample a line via nearest neighbor, rather than linear
+# interpolation (ALUI).  The calling sequence is the same as for ALUIS.
+
+procedure si_samples (a, b, x, npix)
+
+short	a[ARB], b[ARB]		# input, output data arrays
+real	x[ARB]			# sample grid
+int	npix, i
+
+begin
+	do i = 1, npix
+	    b[i] = a[int(x[i])]
+end
+
+
+# SIGLNR -- Get a line of type real from a scaled image.  Block averaging is
+# done by a subprocedure; this procedure gets a line from a possibly block
+# averaged image and if necessary interpolates it to the grid points of the
+# output line.
+
+pointer procedure siglnr (si, lineno, slice_axis, slice)
+
+pointer	si		# pointer to SI descriptor
+int	lineno
+int	slice_axis	# axis from which to slice section if ndim>2
+int	slice		# current slice index
+
+pointer	rawline, tempp, gp
+int	i, buf_y[2], new_y[2], tempi, curbuf, altbuf
+int	npix, nblks_y, ybavg, x1, x2
+real	x, y, weight_1, weight_2
+pointer	si_blkavgr()
+errchk	si_blkavgr
+
+begin
+	npix = SI_NPIX(si,1)
+
+	# Deterine the range of X (in pixels on the block averaged input image)
+	# required for the interpolator.
+
+	gp = SI_GRID(si,1)
+	x1 = SI_XOFF(si)
+	x = Memr[gp+npix-1]
+	x2 = x1 + int(x)
+	if (INTVAL(x))
+	    x2 = x2 - 1
+	x2 = max (x1 + 1, x2)
+
+	gp = SI_GRID(si,2)
+	y = Memr[gp+lineno-1]
+
+	# The following is an optimization provided for the case when it is
+	# not necessary to interpolate in either X or Y.  Block averaging is
+	# permitted.
+
+	if (SI_INTERP(si,1) == NO && SI_INTERP(si,2) == NO)
+	    return (si_blkavgr (SI_IM(si), x1, x2, int(y),
+		SI_BAVG(si,1), SI_BAVG(si,2), slice_axis, slice))
+
+	# If we are interpolating in Y two buffers are required, one for each
+	# of the two input image lines required to interpolate in Y.  The lines
+	# stored in these buffers are interpolated in X to the output grid but
+	# not in Y.  Both buffers are not required if we are not interpolating
+	# in Y, but we use them anyhow to simplify the code.
+
+	if (SI_INIT(si) == YES) {
+	    do i = 1, 2 {
+		if (SI_BUF(si,i) != NULL)
+		    call mfree (SI_BUF(si,i), SI_TYBUF(si))
+		call malloc (SI_BUF(si,i), npix, TY_REAL)
+		SI_TYBUF(si) = TY_REAL
+		buf_y[i] = NOTSET
+	    }
+	    if (OUTBUF(si) != NULL)
+		call mfree (OUTBUF(si), SI_TYBUF(si))
+	    call malloc (OUTBUF(si), npix, TY_REAL)
+	    SI_INIT(si) = NO
+	}
+
+	# If the Y value of the new line is not in range of the contents of the
+	# current line buffers, refill one or both buffers.  To refill we must
+	# read a (possibly block averaged) input line and interpolate it onto
+	# the X grid.  The X and Y values herein are in the coordinate system
+	# of the (possibly block averaged) input image.
+
+	new_y[1] = int(y)
+	new_y[2] = int(y) + 1
+
+	# Get the pair of lines whose integral Y values form an interval
+	# containing the fractional Y value of the output line.  Sometimes the
+	# desired line will happen to be in the other buffer already, in which
+	# case we just have to swap buffers.  Often the new line will be the
+	# current line, in which case nothing is done.  This latter case occurs
+	# frequently when the magnification ratio is large.
+
+	curbuf = 1
+	altbuf = 2
+
+	do i = 1, 2 {
+	    if (new_y[i] == buf_y[i]) {
+		;
+	    } else if (new_y[i] == buf_y[altbuf]) {
+		SWAPP (SI_BUF(si,1), SI_BUF(si,2))
+		SWAPI (buf_y[1], buf_y[2])
+
+	    } else {
+		# Get line and interpolate onto output grid.  If interpolation
+		# is not required merely copy data out.  This code is set up
+		# to always use two buffers; in effect, there is one buffer of
+		# look ahead, even when Y[i] is integral.  This means that we
+		# will go out of bounds by one line at the top of the image.
+		# This is handled by copying the last line.
+
+		ybavg = SI_BAVG(si,2)
+		nblks_y = (IM_LEN (SI_IM(si), 2) + ybavg-1) / ybavg
+		if (new_y[i] <= nblks_y)
+		    rawline = si_blkavgr (SI_IM(si), x1, x2, new_y[i],
+			SI_BAVG(si,1), SI_BAVG(si,2), slice_axis, slice)
+
+		if (SI_INTERP(si,1) == NO) {
+		    call amovr (Memr[rawline], Memr[SI_BUF(si,i)], npix)
+		} else if (SI_ORDER(si) <= 0) {
+		    call si_sampler (Memr[rawline], Memr[SI_BUF(si,i)],
+			Memr[SI_GRID(si,1)], npix)
+		} else {
+		    call aluir (Memr[rawline], Memr[SI_BUF(si,i)],
+			Memr[SI_GRID(si,1)], npix)
+		}
+
+		buf_y[i] = new_y[i]
+	    }
+
+	    SWAPI (altbuf, curbuf)
+	}
+
+	# We now have two line buffers straddling the output Y value,
+	# interpolated to the X grid of the output line.  To complete the
+	# bilinear interpolation operation we take a weighted sum of the two
+	# lines.  If the range from buf_y[1] to buf_y[2] is repeatedly
+	# interpolated in Y no additional i/o occurs and the linear
+	# interpolation operation (ALUI) does not have to be repeated (only the
+	# weighted sum is required).  If the distance of Y from one of the
+	# buffers is zero then we do not even have to take a weighted sum.
+	# This is not unusual because we may be called with a magnification
+	# of 1.0 in Y.
+
+	weight_1 = 1.0 - (y - buf_y[1])
+	weight_2 = 1.0 - weight_1
+
+	if (weight_1 < SI_TOL)
+	    return (SI_BUF(si,2))
+	else if (weight_2 < SI_TOL || SI_ORDER(si) <= 0) 
+	    return (SI_BUF(si,1))
+	else {
+	    call awsur (Memr[SI_BUF(si,1)], Memr[SI_BUF(si,2)],
+		Memr[OUTBUF(si)], npix, weight_1, weight_2)
+	    return (OUTBUF(si))
+	}
+end
+
+
+# SI_BLKAVGR -- Get a line from a block averaged image of type short.
+# For example, block averaging by a factor of 2 means that pixels 1 and 2
+# are averaged to produce the first output pixel, 3 and 4 are averaged to
+# produce the second output pixel, and so on.  If the length of an axis
+# is not an integral multiple of the block size then the last pixel in the
+# last block will be replicated to fill out the block; the average is still
+# defined even if a block is not full.
+
+pointer procedure si_blkavgr (im, x1, x2, y, xbavg, ybavg, slice_axis, slice)
+
+pointer	im			# input image
+int	x1, x2			# range of x blocks to be read
+int	y			# y block to be read
+int	xbavg, ybavg		# X and Y block averaging factors
+int	slice_axis		# axis from which to slice section if ndim>2
+int	slice			# current slice
+
+int	nblks_x, nblks_y, ncols, nlines, xoff, i, j
+int	first_line, nlines_in_sum, npix, nfull_blks, count
+long	vs[IM_MAXDIM], ve[IM_MAXDIM]
+real	sum
+pointer	sp, a, b
+pointer	imgs2r(), imgs3r(), imggsr()
+errchk	imgs2r, imgs3r, imggsr()
+
+begin
+	call smark (sp)
+
+	ncols  = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	xoff   = (x1 - 1) * xbavg + 1
+	npix   = min (ncols, xoff + (x2 - x1 + 1) * xbavg - 1)
+
+	if ((xbavg < 1) || (ybavg < 1))
+	    call error (1, "si_blkavg: illegal block size")
+	else if (x1 < 1 || x2 > ncols)
+	    call error (2, "si_blkavg: column index out of bounds")
+	else if ((xbavg == 1) && (ybavg == 1))
+	    if (IM_NDIM(im) == 2)
+		return (imgs2r (im, xoff, xoff + npix - 1, y, y))
+	    else if (IM_NDIM(im) == 3)
+		return (imgs3r (im, xoff, xoff + npix - 1, y, y, slice, slice))
+	    else {
+		call amovkl (long(1), vs, IM_MAXDIM)
+		call amovkl (long(1), ve, IM_MAXDIM)
+		vs[1] = xoff
+		ve[1] = xoff + npix - 1
+		vs[2] = y
+		ve[2] = y
+		vs[slice_axis] = slice
+		ve[slice_axis] = slice
+		return (imggsr (im, vs, ve, 2))
+	    }
+
+	nblks_x = (npix   + xbavg-1) / xbavg
+	nblks_y = (nlines + ybavg-1) / ybavg
+
+	if (y < 1 || y > nblks_y)
+	    call error (2, "si_blkavg: block number out of range")
+
+	call salloc (b, nblks_x, TY_REAL)
+
+	if (ybavg > 1) {
+	    call aclrr (Memr[b], nblks_x)
+	    nlines_in_sum = 0
+	}
+
+	# Read and accumulate all input lines in the block.
+	first_line = (y - 1) * ybavg + 1
+
+	do i = first_line, min (nlines, first_line + ybavg - 1) {
+	    # Get line from input image.
+	    if (IM_NDIM(im) == 2)
+		a = imgs2r (im, xoff, xoff + npix - 1, i, i)
+	    else if (IM_NDIM(im) == 3)
+		a = imgs3r (im, xoff, xoff + npix - 1, i, i, slice, slice)
+	    else {
+		call amovkl (long(1), vs, IM_MAXDIM)
+		call amovkl (long(1), ve, IM_MAXDIM)
+		vs[1] = xoff
+		ve[1] = xoff + npix - 1
+		vs[2] = i
+		ve[2] = i
+		vs[slice_axis] = slice
+		ve[slice_axis] = slice
+		return (imggsr (im, vs, ve, 2))
+	    }
+
+	    # Block average line in X.
+	    if (xbavg > 1) {
+		# First block average only the full blocks.
+		nfull_blks = npix / xbavg
+		call abavr (Memr[a], Memr[a], nfull_blks, xbavg)
+
+		# Now average the final partial block, if any.
+		if (nfull_blks < nblks_x) {
+		    sum = 0.0
+		    count = 0
+		    do j = nfull_blks * xbavg + 1, npix {
+			sum = sum + Memr[a+j-1]
+			count = count + 1
+		    }
+		    Memr[a+nblks_x-1] = sum / count
+		}
+	    }
+
+	    # Add line into block sum.  Keep track of number of lines in sum
+	    # so that we can compute block average later.
+	    if (ybavg > 1) {
+		call aaddr (Memr[a], Memr[b], Memr[b], nblks_x)
+		nlines_in_sum = nlines_in_sum + 1
+	    }
+	}
+
+	# Compute the block average in Y from the sum of all lines block
+	# averaged in X.  Overwrite buffer A, the buffer returned by IMIO.
+	# This is kosher because the block averaged line is never longer
+	# than an input line.
+
+	if (ybavg > 1)
+	    call adivkr (Memr[b], real(nlines_in_sum), Memr[a], nblks_x)
+
+	call sfree (sp)
+	return (a)
+end
+
+
+# SI_SAMPLER -- Resample a line via nearest neighbor, rather than linear
+# interpolation (ALUI).  The calling sequence is the same as for ALUIR.
+
+procedure si_sampler (a, b, x, npix)
+
+real	a[ARB], b[ARB]		# input, output data arrays
+real	x[ARB]			# sample grid
+int	npix, i
+
+begin
+	do i = 1, npix
+	    b[i] = a[int(x[i])]
+end
diff --git a/pkg/proto/vol/src/i2sun/t_i2sun.x b/pkg/proto/vol/src/i2sun/t_i2sun.x
new file mode 100644
index 00000000..ab8119b7
--- /dev/null
+++ b/pkg/proto/vol/src/i2sun/t_i2sun.x
@@ -0,0 +1,240 @@
+include <imhdr.h>
+include <ctype.h>
+include <mach.h>
+include "i2sun.h"
+
+
+# I2SUN -- IRAF to Sun Rasterfile:  convert either a list of IRAF images 
+# or all slices from a specified axis of a dimension>2 image into a series
+# of Sun rasterfiles.  This format-specific task is primarily used to make
+# movies in the absence of a portable movie/filmloop utility, if a
+# Sun-specific movie task is available.
+# ** The format of the output Sun rasterfiles is hard-coded into this task,
+# ** and thus could diverge from a future Sun format; we do not want to link
+# ** with Sun libraries, as this task should be runnable on other machines.
+
+procedure t_i2sun
+
+pointer	sp, tr, input, im, rfnames, clutfile, transform, cur_rf
+pointer	ulutfile, ulut, colormap, pk_colormap, lut
+int	list, lfd, rfd, nslices, stat, nimages
+int	rheader[RAS_HDR_INTS], ras_maptype, ras_maplength, frame, slice, i, j
+short	lut1, lut2
+bool	use_clut, make_map
+
+pointer immap()
+int	open(), access(), clgeti(), imtopenp(), imtlen(), imtgetim(), read()
+real	clgetr()
+bool	streq(), clgetb()
+
+errchk	open()
+
+begin
+	call smark (sp)
+	call salloc (tr, LEN_TR, TY_STRUCT)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (rfnames, SZ_FNAME, TY_CHAR)
+	call salloc (clutfile, SZ_FNAME, TY_CHAR)
+	call salloc (cur_rf, SZ_FNAME, TY_CHAR)
+	call salloc (transform, SZ_LINE, TY_CHAR)
+	call salloc (pk_colormap, NGREY*3, TY_CHAR)
+	call salloc (colormap, NGREY*3, TY_SHORT)
+	lut = NULL
+	im = NULL
+
+	# Input parameters.
+	list = imtopenp ("input")
+	call clgstr ("output", Memc[rfnames], SZ_FNAME)
+	call clgstr ("clutfile", Memc[clutfile], SZ_FNAME)
+	call clgstr ("ztrans", Memc[transform], SZ_LINE)
+	TR_Z1(tr) = clgetr ("z1")
+	TR_Z2(tr) = clgetr ("z2")
+	TR_ZTRANS(tr) = Z_LINEAR
+	if (streq (Memc[transform], "log"))
+	    TR_ZTRANS(tr) = Z_LOG
+	else if (streq (Memc[transform], "none"))
+	    TR_ZTRANS(tr) = Z_UNITARY
+	else if (streq (Memc[transform], "user")) {
+
+	    # Get user-specified transfer lookup table.
+	    TR_ZTRANS(tr) = Z_USER
+	    call salloc (ulutfile, SZ_FNAME, TY_CHAR)
+	    call clgstr ("ulutfile", Memc[ulutfile], SZ_FNAME)
+
+	    # Borrowed from DISPLAY; mallocs storage for ulut:
+	    call tr_ulut (Memc[ulutfile], TR_Z1(tr), TR_Z2(tr), ulut)
+	}
+	TR_XSIZE(tr) = clgeti ("xsize")
+	TR_YSIZE(tr) = clgeti ("ysize")
+	TR_ORDER(tr) = clgeti ("order")
+	TR_XMAG(tr)  = clgetr ("xmag")
+	TR_YMAG(tr)  = clgetr ("ymag")
+
+	# Get input image axes to map to output frames.  At present we
+	# can only traverse one slice axis.
+	TR_SLICEAXIS(tr) = clgeti ("sliceaxis")
+
+	# Swap bytes in output rasterfile?  (useful when I2SUN run on VAX etc.)
+	TR_SWAPBYTES(tr) = clgetb ("swap")
+
+	# Check if there are no images.
+	nimages = imtlen (list)
+	if (nimages == 0) {
+	    call eprintf (0, "No input images to convert")
+	    goto wrapup_
+	}
+
+	# Open color lookup table file (an existing Sun rasterfile at present)
+	if (access (Memc[clutfile], READ_ONLY, BINARY_FILE) == YES) {
+	    lfd = open (Memc[clutfile], READ_ONLY, BINARY_FILE)
+	    use_clut = true
+	} else
+	    use_clut = false
+
+	# Read color lookup table.
+	make_map = false
+	if (use_clut) {
+	    # Only the color table is used from the rasterfile; ignore all else.
+	    stat = read (lfd, rheader, RAS_HDR_INTS * SZB_CHAR)
+	    if (stat != RAS_HDR_INTS * SZB_CHAR) {
+		call eprintf ("Error reading header from file `%s'\n")
+		    call pargstr (Memc[clutfile])
+		goto wrapup_
+	    }
+	    if (rheader[1] != RAS_MAGIC) {
+		call eprintf ("File `%s' not a valid Sun rasterfile\n")
+		    call pargstr (Memc[clutfile])
+		goto wrapup_
+	    }
+	    ras_maptype = rheader[7]
+	    ras_maplength = rheader[8]
+	    if (ras_maptype != RMT_NONE && ras_maplength > 0) {
+		stat = read (lfd, Memc[colormap], ras_maplength / SZB_CHAR)
+		if (stat != ras_maplength / SZB_CHAR) {
+		    call eprintf ("Error reading colormap from %s\n")
+			call pargstr (Memc[clutfile])
+		    goto wrapup_
+		}
+		# Colormap was already packed on disk.
+		call achtsc (Mems[colormap], Memc[pk_colormap], ras_maplength)
+	    } else {
+		make_map = true
+		call eprintf ("Invalid colormap in %s; using greyscale\n")
+		    call pargstr (Memc[clutfile])
+	    }
+	} else
+	    make_map = true
+
+	if (make_map) {
+	    # Construct a greyscale colormap of same range as IMTOOL.
+	    ras_maptype = RMT_EQUAL_RGB
+	    ras_maplength = NGREY * 3
+	    do i = 1, 3 {
+		Mems[colormap+(i-1)*NGREY] = WHITE
+		do j = COLORSTART+1, COLOREND
+		    Mems[colormap+j-1+(i-1)*NGREY] = j * (WHITE+1) /
+			NGREY
+		Mems[colormap+COLOREND-1+1+(i-1)*NGREY] = WHITE
+		do j = COLOREND+2, NGREY
+		    Mems[colormap+j-1+(i-1)*NGREY] = BLACK
+	    }
+	    call achtsc (Mems[colormap], Memc[pk_colormap], ras_maplength)
+
+	    # Pack to byte stream.
+	    call chrpak (Memc[pk_colormap], 1, Memc[pk_colormap], 1,
+		ras_maplength)
+	}
+
+	# For each IRAF image or band, construct and dispose of a rasterfile.
+	frame = 0
+	while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+	    im = immap (Memc[input], READ_ONLY, 0)
+	    if (IM_NDIM(im) > 2 && TR_SLICEAXIS(tr) > IM_NDIM(im)) {
+		call eprintf ("Specified slice axis invalid for image %s\n")
+		    call pargstr (Memc[input])
+		goto wrapup_
+	    }
+	    nslices = IM_LEN(im, TR_SLICEAXIS(tr))
+	    if (nslices < 1)
+		nslices = 1
+
+	    # Set up spatial transformation (technically, could be different
+	    # for each input image).
+	    call tr_setup (im, tr)
+
+	    # We assume that if any n>2 images are present, the user wants
+	    # all bands dumped out.
+	    do slice = 1, nslices {
+
+		# Construct next rasterfile name and open file; works in
+		# 'append' mode, next higher available frame number.
+		call sprintf (Memc[cur_rf], SZ_FNAME, Memc[rfnames])
+		    call pargi (frame)
+		while (access (Memc[cur_rf], READ_ONLY, BINARY_FILE) == YES) {
+		    frame = frame + 1
+		    call sprintf (Memc[cur_rf], SZ_FNAME, Memc[rfnames])
+			call pargi (frame)
+		}
+		iferr (rfd = open (Memc[cur_rf], NEW_FILE, BINARY_FILE)) {
+		    call eprintf ("Cannot open output rasterfile `%s'\n")
+			call pargstr (Memc[cur_rf])
+		    goto wrapup_
+		}
+		frame = frame + 1
+
+		# Write header to rasterfile:
+		rheader[1] = RAS_MAGIC
+		rheader[2] = TR_XE(tr) - TR_XS(tr) + 1
+		rheader[3] = TR_YE(tr) - TR_YS(tr) + 1
+		rheader[4] = NBITS_FB
+		rheader[5] = rheader[2] * rheader[3]
+		rheader[6] = RMT_STANDARD
+		rheader[7] = ras_maptype
+		rheader[8] = ras_maplength
+		if (TR_SWAPBYTES(tr))
+		    call bswap4 (rheader, 1, rheader, 1, RAS_HDR_INTS*4)
+		call write (rfd, rheader, RAS_HDR_INTS * SZB_CHAR)
+
+		# Write colormap to rasterfile.
+		call write (rfd, Memc[pk_colormap], ras_maplength / SZB_CHAR)
+
+		# Verify user-specified transfer function parameters.
+		if (TR_ZTRANS(tr) == Z_USER) {
+		    call alims (Mems[ulut], U_MAXPTS, lut1, lut2)
+		    if (lut2 < short(DSP_MIN) || lut1 > short(DSP_MAX)) {
+			call eprintf ("User specified greyscales <> range\n")
+			call eprintf ("ulut1=%D, dmin=%D; ulut2=%D, dmax=%D\n")
+			    call pargi (lut1)
+			    call pargi (DSP_MIN)
+			    call pargi (lut2)
+			    call pargi (DSP_MAX)
+		    }
+		    if (!IS_INDEF(TR_Z1(tr)) && !IS_INDEF(TR_Z2(tr)) &&
+			TR_Z2(tr) < IM_MIN(im) || TR_Z1(tr) > IM_MAX(im)) {
+			call eprintf ("User specified intensities <> range\n")
+			call eprintf ("z1=%g, im_min=%g; z2=%g, im_max=%g\n")
+			    call pargr (TR_Z1(tr))
+			    call pargr (IM_MIN(im))
+			    call pargr (TR_Z2(tr))
+			    call pargr (IM_MAX(im))
+			call eprintf ("continuing anyway.\n")
+		    }
+		}
+
+		# Read image pixels and write to rasterfile.
+		call cnv_image (im, slice, tr, ulut, rfd)
+
+		call close (rfd)
+	    }
+	    call imunmap (im)
+	}
+
+wrapup_
+	if (im != NULL)
+	    call imunmap(im)
+	call imtclose (list)
+	call close (rfd)
+	call sfree (sp)
+	if (ulut != NULL)
+	    call mfree (ulut, TY_SHORT)
+end
diff --git a/pkg/proto/vol/src/i2sun/trsetup.x b/pkg/proto/vol/src/i2sun/trsetup.x
new file mode 100644
index 00000000..1b14afb2
--- /dev/null
+++ b/pkg/proto/vol/src/i2sun/trsetup.x
@@ -0,0 +1,32 @@
+include <imhdr.h>
+include "i2sun.h"
+
+
+# TR_SETUP -- Set up spatial transformation parameters.
+
+procedure tr_setup (im, tr)
+
+pointer	im			# An input image descriptor
+pointer	tr			# Transformation structure
+
+int	ncols, nlines
+
+begin
+	ncols  = IM_LEN(im,COL)
+	nlines = IM_LEN(im,LINE)
+
+	# Determine output raster dimensions.  
+	TR_XS(tr) = 1
+	TR_XE(tr) = ncols
+	if (!IS_INDEFI(TR_XSIZE(tr)))
+	    TR_XE(tr) = max (1, TR_XSIZE(tr))
+	else if (TR_XMAG(tr) != 1.0)
+	    TR_XE(tr) = max (1, ncols * int(TR_XMAG(tr)))
+
+	TR_YS(tr) = 1
+	TR_YE(tr) = nlines
+	if (!IS_INDEFI(TR_YSIZE(tr)))
+	    TR_YE(tr) = max (1, TR_YSIZE(tr))
+	else if (TR_YMAG(tr) != 1.0)
+	    TR_YE(tr) = max (1, nlines * int(TR_YMAG(tr)))
+end
diff --git a/pkg/proto/vol/src/i2sun/trulut.x b/pkg/proto/vol/src/i2sun/trulut.x
new file mode 100644
index 00000000..4787b9b3
--- /dev/null
+++ b/pkg/proto/vol/src/i2sun/trulut.x
@@ -0,0 +1,128 @@
+include	<error.h>
+include	<ctype.h>
+include	"i2sun.h"
+
+# TR_ULUT -- Generates a look up table from data supplied by user.  The
+# data is read from a two column text file of intensity, greyscale values.
+# The input data are sorted, then mapped to the x range [0-4095].  A 
+# piecewise linear look up table of 4096 values is then constructed from 
+# the (x,y) pairs given.  A pointer to the look up table, as well as the z1 
+# and z2 intensity endpoints, is returned.
+
+procedure tr_ulut (fname, z1, z2, lut)
+
+char	fname[SZ_FNAME]		# Name of file with intensity, greyscale values
+real	z1			# Intensity mapped to minimum gs value
+real	z2			# Intensity mapped to maximum gs value
+pointer	lut			# Look up table - pointer is returned
+
+pointer	sp, x, y
+int	nvalues, i, j, x1, x2, y1
+real	delta_gs, delta_xv, slope
+errchk	ds_rlut, ds_sort, malloc		
+
+begin
+	call smark (sp)
+	call salloc (x, U_MAXPTS, TY_REAL)	
+	call salloc (y, U_MAXPTS, TY_REAL)
+
+	# Read intensities and greyscales from the user's input file.  The
+	# intensity range is then mapped into a standard range and the 
+	# values sorted.
+
+	call ds_rlut (fname, Memr[x], Memr[y], nvalues)
+	call alimr (Memr[x], nvalues, z1, z2)
+	call amapr (Memr[x], Memr[x], nvalues, z1, z2, real(U_Z1), real(U_Z2))
+	call ds_sort (Memr[x], Memr[y], nvalues)
+
+	# Fill lut in straight line segments - piecewise linear
+	call malloc (lut, U_MAXPTS, TY_SHORT)	
+	do i = 1, nvalues-1 {
+	    delta_gs = Memr[y+i] - Memr[y+i-1]
+	    delta_xv = Memr[x+i] - Memr[x+i-1]
+	    slope = delta_gs / delta_xv
+	    x1 = int (Memr[x+i-1]) 
+	    x2 = int (Memr[x+i])
+	    y1 = int (Memr[y+i-1])
+	    do j = x1, x2
+		Mems[lut+j] = y1 + slope * (j-x1)
+	}
+	Mems[lut+U_MAXPTS-1] = y1 + (slope * U_Z2)
+
+	call sfree (sp)
+end
+
+
+# DS_RLUT -- Read text file of x, y, values.
+
+procedure ds_rlut (utab, x, y, nvalues)
+
+char	utab[SZ_FNAME]		# Name of list file
+real	x[U_MAXPTS]		# Array of x values, filled on return
+real	y[U_MAXPTS]		# Array of y values, filled on return
+int	nvalues			# Number of values in x, y vectors - returned
+
+int	n, fd
+pointer	sp, lbuf, ip
+real	xval, yval
+int	getline(), open()
+errchk	open, sscan, getline, salloc
+
+begin
+	call smark (sp)
+	call salloc (lbuf, SZ_LINE, TY_CHAR)
+
+	iferr (fd = open (utab, READ_ONLY, TEXT_FILE))
+	    call error (1, "Error opening user lookup table")
+
+	n = 0
+	while (getline (fd, Memc[lbuf]) != EOF) {
+	    # Skip comment lines and blank lines.
+	    if (Memc[lbuf] == '#')
+		next
+	    for (ip=lbuf;  IS_WHITE(Memc[ip]);  ip=ip+1)
+		;
+	    if (Memc[ip] == '\n' || Memc[ip] == EOS)
+		next
+
+	    # Decode the points to be plotted.
+	    call sscan (Memc[ip])
+		call gargr (xval)
+		call gargr (yval)
+
+	    n = n + 1
+	    if (n > U_MAXPTS)
+	        call error (2,
+		    "Intensity transformation table cannot exceed 4096 values")
+
+	    x[n] = xval
+	    y[n] = yval
+	}
+
+	nvalues = n
+	call close (fd)
+	call sfree (sp)
+end
+
+
+# DS_SORT -- Bubble sort of paired arrays.  
+
+procedure ds_sort (xvals, yvals, nvals)
+
+real	xvals[nvals]		# Array of x values
+real	yvals[nvals]		# Array of y values
+int	nvals			# Number of values in each array
+
+int	i, j
+real	temp
+define	swap	{temp=$1;$1=$2;$2=temp}
+
+begin
+	for (i=nvals;  i > 1;  i=i-1)
+	    for (j=1;  j < i;  j=j+1) 
+		if (xvals[j] > xvals[j+1]) {
+		    # Out of order; exchange y values
+		    swap (xvals[j], xvals[j+1])
+		    swap (yvals[j], yvals[j+1])
+		}
+end
diff --git a/pkg/proto/vol/src/i2sun/x_i2sun.x b/pkg/proto/vol/src/i2sun/x_i2sun.x
new file mode 100644
index 00000000..20a36169
--- /dev/null
+++ b/pkg/proto/vol/src/i2sun/x_i2sun.x
@@ -0,0 +1,4 @@
+# X_I2SUN -- Task statement for I2SUN, used only for debugging (normally task
+# resides in X_PVOL.E.
+
+task	i2sun = t_i2sun