Repatch (from linux) of OSX IRAF

48 files changed, 6293 insertions, 0 deletions

diff --git a/pkg/proto/vol/src/doc/concept.hlp b/pkg/proto/vol/src/doc/concept.hlp
new file mode 100644
index 00000000..26a814b9
--- /dev/null
+++ b/pkg/proto/vol/src/doc/concept.hlp
@@ -0,0 +1,177 @@
+.help volumes
+
+[OUT OF DATE (Jan 89); this is an original pre-design document]
+
+.ce
+Conceptual Model for 2D Projections of Rotating 3D Images
+
+Consider the problem of visualizing a volume containing emissive and
+absorptive material.  If we had genuine 3D display tools, we could imagine
+something like a translucent 3D image display that we could hold in our
+hand, peer into, and rotate around at will to study the spatial distribution
+of materials inside the volume.
+
+Lacking a 3D display device, we can resort to 2D projections of the interior
+of the volume.  In order to render absorptive material, we need a light
+source behind the volume; this light gets attenuated by absorption as it
+passes through the volume toward the projection plane.  In the general case
+light emitted within the volume contributes positively to the projected
+light intensity, but it also gets attenuated by absorbing material between
+it and the projection plane.  At this point the projection plane has a range
+of intensities representing the combined absorption and emission of material
+along columns through the volume.  But looking at a single 2D projection,
+there would be no way to determine how deep in the original volume a particular
+emitting or absorbing region lay.  One way around this is to cause the volume
+to rotate, making a series of 2D projections.  Playing the projections back
+as a movie gives the appearance of seeing inside a translucent rotating
+volume.
+
+Modelling the full physics of light transmission, absorption, refraction,
+etc. with arbitrary projective geometries would be quite computationally
+intensive and could rival many supercomputer simulations.  However, it is
+possible to constrain the model such that an effective display can be
+generated allowing the viewer to grasp the essential nature of the spatial
+relationships among the volume data in reasonable computational time.  This
+is called volume visualization, which can include a range of display
+techniques approximating the actual physics to varying extents.  There is
+some debate whether visualization problems can best be attacked by simplified
+direct physical models or by different models, such as ones that might better
+enhance the \fBperception\fR of depth.  We will stick with direct physical
+models here, though simplified for computer performance reasons.
+
+For computational purposes we will constrain the projection to be orthogonal,
+i.e. the light source is at infinity, so the projection rays are all parallel.
+With the light source at infinity behind the volume (a datacube), we need not
+model reflection at all.  We will also ignore refraction (and certainly
+diffraction effects).
+
+We can now determine a pixel intensity on the output plane by starting
+at the rear of the column of voxels (volume elements) that project from
+the datacube onto that pixel.  At each successive voxel along that column
+we will attenuate the light we started with by absorption, and add to it
+any light added by emission.  If we consider emission (voxel intensity)
+alone, the projection would just be the sum of the contributing intensities.
+Absorption alone would simply decrease the remaining transmitted light
+proportionally to the opacity of each of the voxels along the column.
+Since we are combining the effects of absorption and emission, the ratio
+of the intensity of the original incident light to that of the interior
+voxels is important, so we will need a beginning intensity.
+
+The opacities have a physical meaning in the model.  However, we are more
+interested here in visualizing the volume interior than in treating it as
+a pure physical model, so we add an opacity transformation function.  This
+allows us to generate different views of the volume interior without having
+to modify all the raw opacity values in the datacube.  For maximum flexibility
+we would like to be able to modify the opacity function interactively, e.g.
+with a mouse, but this would amount to computing the projections in real
+time and is not likely at present.
+
+.nf
+
+Let:	i     = projected intensity before considering the current
+	        voxel
+	i'    = intensity of light after passing through the current
+	        voxel
+	I0    = initial light intensity (background iillumination
+		before encountering the volume)
+	Vo    = opacity at the current voxel, range 0:1 with
+	        0=transparent, 1=opaque
+	Vi    = intensity at the current voxel
+	f(Vo) = function of the opacity at the current voxel,
+		normalized to the range 0:1
+	g(Vi) = function of the voxel's intensity, normalized
+		to the range 0:1
+
+Then:	i' = i * (1 - f(Vo)) + g(Vi)
+		[initial i = Imax, then iterate over all voxels in path]
+
+.fi
+
+We want to choose the opacity and intensity transformation functions in such
+a way that we can easily control the appearance of the final projection.
+In particular, we want to be able to adjust both the opacity and intensity
+functions to best reveal the interior details of the volume during a
+rotation sequence.  For example, we might want to eliminate absorption
+"noise" so that we can see through it to details of more interest, so we
+need a lower opacity cutoff.  Likewise, we would want an upper opacity
+cutoff above which all voxels would appear opaque.  We will need the same
+control over intensity.
+
+.nf
+Let:	o1   = lower voxel opacity cutoff
+	o2   = upper voxel opacity cutoff
+	t1   = lower transmitted intensity cutoff
+	t2   = upper transmitted intensity cutoff
+	i1   = lower voxel intensity cutoff
+	i2   = upper voxel intensity cutoff
+	Imax = output intensity maximum for int transform function
+.fi
+
+Now all we need is the form of the opacity and intensity functions between
+their input cutoffs.  A linear model would seem to be useful, perhaps with
+logarithmic and exponential options later to enhance the lower or upper
+end of the range.  f(Vo) is constrained to run between 0 and 1, because
+after being subtracted from 1.0 it is the intensity attenuation factor
+for the current voxel.
+
+.nf
+
+		Opacity Transformation Function f(Vo):
+
+	        { Vo < o1       : 	    0.0		    }
+	        {					    }
+	        { o1 <= Vo < o2 : (t2 - (Vo - o1)(t2 - t1)) }
+	        {		  (	---------	  ) }
+	f(Vo) = {		  (	(o2 - o1)	  ) }
+	        {		  ------------------------- }
+	        {			    I0	    }
+	        {					    }
+	        { o2 <= Vo      : 	    1.0		    }
+
+
+backg. int.  I0-|
+		|
+	     t2-|------			Transmitted Intensity
+		|        `		as function of opacity
+	        |           `		(ignoring independent
+i * (1 - f(Vo))-|..............`	voxel intensity contri-
+                |               . `	bution)
+    		|               .    `
+		|               .       `
+	     t1-|		.	  |
+		|		.	  |
+		+____________________________________________
+		       |        |         |
+		       o1	Vo 	  o2
+			   Voxel Opacity
+
+       ------------------------------------------------------------
+
+		Intensity Transformation Function g(Vi):
+
+		{ Vi < i1       :    0.0
+		{
+		{ i1 <= Vi < i2 : (Vi - i1) * Imax
+		{		  ---------
+	g(Vi) = {		  (i2 - i1)
+		{
+		{ i2 <= Vi	:    Imax
+		{
+		{
+		{
+
+
+	        |
+		|
+	   Imax-|                       ---------------------
+		|                    /
+	  g(Vi)-|................./
+      		|              /  .
+		|           /     .
+		|        /        .
+	    0.0	+___________________________________________
+		     |            |     |
+		     i1	          Vi    i2
+			Voxel Intensity
+.fi
+
diff --git a/pkg/proto/vol/src/doc/i2sun.hlp b/pkg/proto/vol/src/doc/i2sun.hlp
new file mode 100644
index 00000000..70448d64
--- /dev/null
+++ b/pkg/proto/vol/src/doc/i2sun.hlp
@@ -0,0 +1,152 @@
+.help i2sun Oct88 local
+.ih
+NAME
+i2sun -- convert IRAF images to Sun rasterfiles
+.ih
+USAGE
+i2sun input output z1 z2
+.ih
+PARAMETERS
+.ls input
+Input image template, @file, n-dimensional image, or combination.
+.le
+.ls output
+Root template for output images, e.g. "home$ras/frame.%d".
+.le
+.ls clutfile
+Previously saved Sun rasterfile (e.g. output from IMTOOL), containing the
+color/greyscale lookup table information to be passed along to each output
+frame.  Standard ones can be saved and used with any number of images (e.g.
+"pseudo.ras").
+.le
+.ls z1 = INDEF, z2 = INDEF
+Minimum and maximum pixel/voxel intensities to scale to full output
+color/greyscale range.  Both are required parameters, and will apply to all
+images in the sequence.
+.le
+.ls ztrans = "linear"
+Intensity transformation on input data (linear|log|none|user).
+If "user", you must also specify \fIulutfile\fR.
+.le
+.ls ulutfile
+Name of text file containing the look up table when \fIztrans\fR = user.
+The table should contain two columns per line; column 1 contains the
+intensity, column 2 the desired greyscale output.
+.le
+.ls xsize = INDEF, ysize = INDEF
+If specified, these will be the dimensions of all output Sun rasterfiles
+in pixels.  The default will be the same size as the input images (which
+could vary, though this would create a jittery movie).
+.le
+.ls xmag = 1.0, ymag = 1.0
+Another way to specify output rasterfile dimensions.  These are the 
+magnification factors to apply to the input image dimensions.
+.le
+.ls order = 1
+Order of the interpolator to be used for spatially interpolating the image.
+The current choices are 0 for pixel replication, and 1 for bilinear
+interpolation.
+.le
+.ls sliceaxis = 3
+Image axis from which to cut multiple slices when input image dimension is
+greater than 2.  Only x-y sections are allowed, so \fIsliceaxis\fR must
+be 3 or greater.
+.le
+.ls swap = no
+Swap rasterfile bytes on output?  Used when rasterfiles are being written
+to a computer with opposite byte-swapping from that of the home computer
+(e.g. between VAX and Sun).
+.le
+
+
+.ih
+DESCRIPTION
+
+Given a series of IRAF images, an intensity transformation, and a file
+containing color/greyscale lookup table information, produces one 2d image
+in Sun rasterfile format for each 2D IRAF image.  This is a temporary task
+usually used as a step in creating filmloops for playback by a Sun Movie
+program.
+
+The input images may be specified as an image template ("zoom*.imh"),
+an "@" file ("@movie.list"), or as an n-dimensional image from which to
+create multiple 2d rasterfiles.  If any images in a list are nD images,
+all 2d sections from the specified \fIsliceaxis\fR will be written out
+(default = band or z axis).  At present, only x-y sections may be made,
+i.e. the slice axis must be axis 3 or higher.
+
+The minimum and maximum pixel/voxel intensities, z1 and z2, must be specified
+as it would be not only inefficient to calculate the full zrange of
+each image in a sequence, but would also make very jumpy movies.
+Between input intensities z1 and z2, the pixel intensities may be transformed
+according to the \fIztrans\fR parameter: "linear", "log10", "none",
+or "user".
+
+When \fIztrans\fR = "user", a look up table of intensity values and their
+corresponding greyscale levels is read from the file specified by the
+\fIulutfile\fR parameter.  From this information, a piecewise linear
+look up table containing 4096 discrete values is composed.  The text
+format table contains two columns per line; column 1 contains the
+intensity, column 2 the desired greyscale output.  The greyscale values
+specified by the user must match those available on the output device.
+Task \fIshowcap\fR can be used to determine the range of acceptable
+greyscale levels.  
+
+A color table file (\fIclutfile\fR) may be produced on a Sun workstation from
+IMTOOL (see IMTOOL manual page, R_RASTERFILE parameter and Imcopy function).
+This file may be specified to I2SUN as the \fIclutfile\fR parameter.
+Likewise, any rasterfiles previously created with
+I2SUN may be used as input clutfiles.
+
+The output rasterfile dimensions may be larger or smaller than the input 
+images (see parameters \fIxsize\fR and \fIysize\fR, or \fIxmag\fR and
+\fIymag\fR).  The parameter \fIorder\fR controls the mode of interpolation;
+0=pixel replication, 1=bilinear.
+
+If the output rasterfiles are being sent to a computer with opposite
+byte-swapping characteristics, set \fIswap\fR = yes (e.g., when running
+I2SUN on a VAX, with output to a Sun).
+
+
+.ih
+EXAMPLES
+
+.nf
+1.  Produce a series of Sun rasterfiles in tmp$mydir/movie/,
+    using a pseudocolor color table file saved earlier, with
+    input greylevels scaled between 10 and 100.
+
+    cl> i2sun nzoom*.imh tmp$mydir/movie/frame.%d \
+	home$colors/pseudo.ras 10 100
+
+2.  Make a movie through the z, or band, axis of a datacube.
+
+    cl> i2sun cube tmp$cubemovie/frame.%d 1 256 
+
+3.  Make a movie through the 4th, or hyper-axis of a datacube,
+    holding image band 10 constant.
+
+    cl> i2sun hypercube[*,*,10,*] tmp$movie/frame.%d 1 256 \
+	sliceaxis=4
+
+4.  Run I2SUN on a VAX, with output to a Sun.
+
+    cl> i2sun @imlist sunnode!home$ras/frame.%d 1 256 swap+
+
+.fi
+
+.ih
+TIMINGS
+49 seconds (1 sec/frame) to produce 50 100*100 rasterfiles from a
+100*100*50 datacube with no magnification, on a diskless Sun-3/110
+using NFS to Eagle disks on a lightly loaded Sun-3/160 fileserver
+(load factor < 1.5).  
+5 minutes for the same with a magnification factor of 2 in both x and y,
+bilinear interpolation.
+20 minutes for the same with a magnification factor of 5 in both x and y.
+.ih
+BUGS
+.ih
+SEE ALSO
+display, imtool, volumes.pvol
+.endhelp
diff --git a/pkg/proto/vol/src/doc/im3dtran.hlp b/pkg/proto/vol/src/doc/im3dtran.hlp
new file mode 100644
index 00000000..75fd85fe
--- /dev/null
+++ b/pkg/proto/vol/src/doc/im3dtran.hlp
@@ -0,0 +1,85 @@
+.help im3dtran Jan89 volumes
+.ih
+NAME
+im3dtran -- 3d image transpose, any axis to any other axis
+.ih
+USAGE
+im3dtran input output 
+.ih
+PARAMETERS
+.ls input
+Input 3d image (datacube).
+.le
+.ls output
+Transposed datacube.
+.le
+.ls len_blk = 128
+Size in pixels of linear internal subraster.  IM3DTRAN will try to transpose
+a subraster up to len_blk cubed at one time.  Runtime is much faster with
+larger \fBlen_blk\fR, but the task may run out of memory.
+.le
+.ls new_x = 3
+New x axis = old axis (1=x, 2=y, 3=z).  Default (3) replaces new x with old z.
+.le
+.ls new_y = 2
+New y axis = old axis.  Default (2) is identity.
+.le
+.ls new_z = 1
+New z axis = old axis.  Default (1) replaces new z with old x.
+.le
+
+.ih
+DESCRIPTION
+
+IM3DTRAN is very similar to IMAGES.IMTRANSPOSE, except that it can accomplish
+3d image transposes.  In 3 dimensions, it is necessary to specify which old
+axes map to the new axes.  In all cases, IM3DTRAN maps old axis element 1 to
+new axis element 1, i.e. it does not flip axes, just transposes them.
+
+If one wants to use IM3DTRAN to rotate a datacube 90 degrees in any direction,
+it is necessary to use a combination of flip and transpose (just like in the
+2d case).  For example, let the original datacube be visualized with its
+origin at the lower left front when seen by the viewer, with the abscissa
+being the x axis (dim1), ordinate the y axis (dim2), and depth being the
+z axis (dim3), z increasing away from the viewer or into the datacube [this
+is a left-handed coordinate system].  One then wants to rotate the datacube
+by 90 degrees clockwise about the y axis when viewed from +y (the "top");
+this means the old z axis becomes the new x axis, and the old x axis becomes
+the new z axis, while new y remains old y.  In this case the axis that must
+be flipped prior to transposition is the \fBx axis\fR; see Example 1.
+
+The parameter \fBlen_blk\fR controls how much memory is used during the
+transpose operation.  \fBlen_blk\fR elements are used in each axis at a
+time, or a cube len_blk elements on a side.  If \fBlen_blk\fR is too large,
+the task will abort with an "out of memory" error.  If it is too small,
+the task can take a very long time to run.  The maximum size of len_blk
+depends on how much memory is available at the time IM3DTRAN is run,
+and the size and datatype of the image to be transposed.
+
+.ih
+EXAMPLES
+
+.nf
+1.  For an input datacube with columns = x = abscissa, lines = y = ordinate,
+    and bands = z = depth increasing away from viewer, and with the image
+    origin at the lower left front, rotate datacube 90 degrees clockwise
+    around the y axis when viewed from +y (top):
+
+    cl> im3dtran input[-*,*,*] output 3 2 1
+
+.fi
+
+.ih
+TIMINGS
+
+[Not available yet]
+
+.ih
+BUGS
+
+[Not available yet]
+
+.ih
+SEE ALSO
+pvol i2sun
+.endhelp
diff --git a/pkg/proto/vol/src/doc/imjoin.hlp b/pkg/proto/vol/src/doc/imjoin.hlp
new file mode 100644
index 00000000..6d7a59a1
--- /dev/null
+++ b/pkg/proto/vol/src/doc/imjoin.hlp
@@ -0,0 +1,76 @@
+.help imjoin Jan89 images
+.ih
+NAME
+imjoin -- join input images into output image along specified axis
+.ih
+USAGE
+imjoin input output 
+.ih
+PARAMETERS
+.ls input
+Input images or @file
+.le
+.ls output
+Output joined image
+.le
+.ls joindim = 1
+Image dimension along which the input images will be joined.
+.le
+.ls outtype = ""
+Output image datatype.  If not specified, defaults to highest precedence
+input image datatype.
+.le
+
+.ih
+DESCRIPTION
+
+IMJOIN concatenates a set of input images into a single output image,
+in a specified dimension only.  For example, it can join a set of one
+dimensional images into a single, long one dimensional image, or a
+set of one dimensional images into a single two dimensional image.
+IMJOIN may be used to piece together datacubes into larger
+datacubes, either in x, y, or z; likewise with higher dimensional images.
+
+For joining a set of 1 or 2 dimensional images in both x and y at the same
+time, see IMMOSAIC.  For stacking images of any dimension into an image
+of the next higher dimension, see IMSTACK.  Although IMJOIN can also
+stack a set of images into a single higher dimensional image, IMSTACK
+is more efficient for that operation.  In most cases, IMJOIN must keep
+all input images open at the same time, while IMSTACK does not (there may
+be limitations on the number of files that can be kept open at one time).
+Use IMJOIN primarily when joining a set of images along any dimension that
+is not the next higher one from that of the input images.
+
+.ih
+EXAMPLES
+
+.nf
+1.  Join a list of one dimensional spectra into a single long image.
+
+    cl> imjoin @inlist output 1
+
+2.  Join three datacubes along the z direction.
+
+    cl> imjoin c1,c2,c3 fullxcube 3
+
+.fi
+
+.ih
+TIMINGS
+
+Join 10 5000 column type short spectra into one 50000 column image:
+6 seconds on a diskless Sun-3.  
+
+Join 2 512*512 images:  28 seconds on diskless Sun-3.  Join 2 50*50*50
+datacubes in x, y, or z:  15 seconds.
+
+.ih
+BUGS
+
+There may be limitations on the number of input images that can be handled
+in one execution on some systems.
+
+.ih
+SEE ALSO
+immosaic, imstack, imslice
+.endhelp
diff --git a/pkg/proto/vol/src/doc/proj.hlp b/pkg/proto/vol/src/doc/proj.hlp
new file mode 100644
index 00000000..f0ed8a3e
--- /dev/null
+++ b/pkg/proto/vol/src/doc/proj.hlp
@@ -0,0 +1,139 @@
+.help volumes Jan89 "Volume Rotation-Projection Algorithm"
+
+.ce
+Volume Rotation-Projection Algorithm
+.ce
+January 1989
+
+.sh
+Introduction
+
+See help for VOLUMES and PVOL for general information.  Here we describe
+the volume projection algorithm used in PVOL.
+
+.sh
+Algorithms for Collecting Object Voxels that Project onto Image Plane
+
+PVOL is a task for making successive projections through a 3d image onto
+2d images placed along a great circle around an input datacube, with varying
+degrees of translucency.  The technique of viewing successive projections
+around the input datacube causes interior features to appear to "orbit"
+the axis of datacube rotation; the apparent orbital radii generate the
+illusion of seeing in three dimensions.  We limit ourselves to parallel rather
+than perspective projections as the computations are simpler and the resulting
+images preserve distance ratios.
+
+When we are considering orthogonal projections only, the 3D problem becomes
+a 2D problem geometrically, collapsed into a plane at right angles to the
+datacube rotation axis.  Otherwise a full 3D solution would be needed.
+To keep things straight, I will use "object voxel"
+to represent voxels from the input volume image and "image pixel" to represent
+output pixels in the projection plane.
+
+In addition to the projections being parallel, we also want them centered
+and the projection plane perpendicular to the projection rays (we always want
+to be looking toward the center of the volume regardless of the rotation angle).
+Thus we will always orient the center of the projection plane perpendicular
+to the ray passing through the center of the volume for the given rotation
+angle.  
+
+Methods in the literature include back-to-front (BTF) and front-to-back (FTB)
+traversals, digital differential analyzer (DDA) techniques, and octree
+encoding.  Because of the nature of our light-transmission algorithm, we
+must choose a BTF approach.  For standard ray-tracing applications, involving
+discrete objects within the volume image space, octree techniques can be
+the most efficient, depending on the ratio of filled to un-filled space and
+number of objects.  However, for arbitrary voxel images (no explicit geometric
+surfaces included, so every voxel must be examined) simpler techniques are
+considered more efficient.  There are basically two approaches:
+[1] image-plane order:  build up the output image one line at a time by
+computing all contributing voxels, and
+[2] volume-image order:  traverse the voxels one line at a time, building
+up the output image in successive overlapping sheets.
+
+The image-plane order approach is similar to rasterizing a line segment, namely
+the projection ray through the lattice of voxels.  Examples are the incremental
+algorithm discussed in Foley and Van Dam (p. 432), implemented with
+modifications in the IRAF SGI kernel, and Bresenham's algorithm, outlined in
+the same place.  Both methods can be extended to include information from
+extra surrounding voxels, similar to anti-aliasing problems, and this may
+be necessary for effective volume projections, especially of small spatial
+resolution volumes.  This approach may not necessarily be the most efficient
+if the volume image cannot be held in memory and must be accessed randomly
+from disk.  Initially, we will code this algorithm only for the case where the
+rotation is around the X axis of the volume and the viewing direction is
+perpendicular to that axis.
+
+[Discussion of various algorithms for determining which set of voxels gets
+included along a given projection ray follows.  After this was coded, it
+became apparent that runtime was largely dominated by the voxel memory
+accesses after the voxel lists have been prepared.  Consequently, the 
+incremental algorithm is all that is now used.]
+
+The straightforward incremental algorithm would be the simplest to implement,
+though not the most efficient.  Bresenham's algorithm, extended to include
+information from fractionally pierced neighboring voxels, would be more
+efficient as it need not use any real variables, and therefore does not
+require rounding.  Both these methods choose a single ray at a time hitting
+the projection plane, and proceed along that ray, determining which voxels
+contribute, and their weights, which are proportional to the path length
+of the ray through the object voxels.  By proceeding from back to front, we are
+guaranteed that each contributing voxel from the volume succeeds any previous
+one arriving at the current output pixel.  Thus, we can use the output
+pixel to store the results of any previous light transmission and absorption
+operation, and feed that value back in to combine with the properties of
+the next contributing volume voxel.  This method fills up the image plane
+in line-sequential order.  Of course, we determine the list of object voxels
+contributing to a given line of image pixels only once per rotation.
+
+In the volume-image order approach the input voxels are traversed line by line
+in any correct BTF order; they can always be accessed band by band if that is
+the disk storage order.  This method fills up the image plane in successive
+sheets, continually updating the image pixels previously written as it goes.
+Determining which image pixel should be hit by the current object voxel
+requires a transformation matrix.  However, the information in the matrix can
+be pre-multiplied with all possible values of voxel coordinates and stored in
+a lookup table, resulting in much more efficient code than a straightforward
+matrix multiplication for each object voxel (Frieder, Gordon, and Reynolds,
+IEEE CG&A, Jan 1985, p. 52-60).  Due to the significantly increased 
+computation time, this approach should only be used when datacube projections
+are desired along any arbitrary 3D orientation.
+
+In the current implementation only rotations by PVOL around the image X
+axis are allowed.  If rotation is desired about either Y or Z, it is easy
+to first rotate the input image, then run PVOL around the new X axis.
+See D3TRANSPOSE [IMTRANS3D?] for help in rotating datacubes.
+
+.sh
+Memory Management
+
+Now we know how to construct a list of indices of input voxels in
+BTF order that impinge upon a given pixel in the projection plane.
+The original PVOL prototype used line-oriented image i/o to access
+the datacube.  Profiles showed 90% of task execution time spent in
+OS-level reads.  Various other approaches were investigated, which
+determined that actual voxel-value i/o was the most important factor
+in performance.  Since "in-core" i/o is the fastest, the problem became
+one of getting as much of the input datacube into main memory as possible.  
+
+A task maximum working set size parameter was added, and code for attempting
+to grab this much memory, then cascading down to a reasonable amount if
+the requested amount was too much (had adverse effects on PVOL or other
+processes).  Given a fixed amount of available memory smaller than that 
+required to hold the entire datacube in memory, the fastest way is to
+volume-project through successive groups of YZ slices.  A single YZ slice
+of the datacube is sufficient for projecting any and all great-circle
+orientations (360 degrees around the X axis).  The more YZ slices that
+can be held in memory, the better.  If there is room for N YZ slices at
+a time, and there are COLUMNS voxels in the X direction, then all volume
+rotations must be made in each of (COLUMNS/N) passes.  
+
+This approach sped things up by about a factor of 20 over random 
+line-oriented i/o.  For very large datacubes (order of 500 voxels on
+a side) there are on the order of 10 passes required when the task
+working set is in the 10Mb range.  Clearly available memory and/or super
+fast disk i/o, dominates volume rotations.  A general purpose workstation
+with enough main memory can apparently approach the speed of the specialized
+processors usually used in volume rendering.
+
+
diff --git a/pkg/proto/vol/src/doc/pvol.hlp b/pkg/proto/vol/src/doc/pvol.hlp
new file mode 100644
index 00000000..30ae4f38
--- /dev/null
+++ b/pkg/proto/vol/src/doc/pvol.hlp
@@ -0,0 +1,398 @@
+.help pvol Jan89 volumes
+.ih
+NAME
+pvol -- project rotations of a volume datacube onto series of 2d images
+.ih
+USAGE
+pvol input output 
+.ih
+PARAMETERS
+.ls input
+Input 3d or 4d image (datacube).
+.le
+.ls output
+Output datacube, one image band per rotation (type real only).
+.le
+.ls nframes = (360 / \fBdegrees\fR)
+Number of frames to generate, 1 per rotation.
+.le
+.ls degrees = 10
+Number of degrees to rotate datacube for each successive projection.
+.le
+.ls theta0 = 0.0
+Initial projection angle for rotation sequence by \fBdegrees\fR increments.
+Measured counterclockwise from +x axis when looking back toward the image
+origin.
+.le
+.ls ptype = 2
+Projection type;
+1 = opacity:  attenuation along projection column by voxel opacity value.
+2 = average voxel intensity along projection column.
+3 = sum of voxel intensities.
+4 = proportional distance weighting: voxel intensity
+along projection column weighted by (curvoxel / voxels_in_column)
+**\fBdispower\fR.
+5 = mod(n):  same as proportional distance weighting, but use only voxel values
+which match mod(normalized_voxel * 100) = \fBmodn\fR.
+6 = use last voxel value within cutoffs only.
+.le
+.ls imin, imax = INDEF
+Input voxel intensity ranges within which to apply intensity transformation.
+Defaults to input image min and max if not specified (see comments below).
+.le
+.ls omin, omax = INDEF
+Input voxel opacity ranges within which to apply opacity transformation.
+Defaults to input image min and max if not specified (see comments below).
+.le
+.ls amin, amax = 0.0, 1.0
+Attenuation factor minimum and maximum for ptype=1 (opacity).  Voxel values
+<= omin map to attenuation factor amin, >= omax map to attenuation amax.
+.le
+.ls izero = 1.0
+Initial background iillumination intensity when \fBptype\fR = 1 (opacity).
+This intensity will be attenuated consecutively by (transformed voxel_value *
+\fBoscale\fR)
+along the projection column toward the projection plane.
+.le
+.ls oscale = 1.0
+Voxel opacity scale factor.  Multiplied by voxel value before attenuating
+remaining light along projection column for \fBptype\fR = 1.
+.le
+.ls opacelem = 1
+Opacity element in 4th dimension of input image.  When input image is 4d,
+and there are two elements in the 4th dimension, the \fBopacelem\fR element
+will be treated as opacity and the other will be considered intensity.
+.le
+.ls dispower = 2.0
+Inverse distance weighting power for \fBptype\fR = 4,5.  Voxel intensities will
+be multiplied by (voxel position in column / voxels in column) **
+\fBdispower\fR before being summed into the output projection pixel.
+.le
+.ls discutoff = no
+When distance weighting, measure the distance within that set of projecting
+voxels that lies between the intensity cutoffs rather than from
+the edges of the datacube.  Usually results in faster run times and is
+appropriate when the interior of a well-defined object is of interest
+rather than its placement inside the datacube.
+.le
+.ls modn = 10
+For ptype=5, only voxel values satisfying mod (int (voxval * 100.0)) =
+\fBmodn\fR will be proportional distance-weighted and summed into
+projection pixel.  Useful for viewing volume interiors with high contrast
+voxel values (like solid objects in an otherwise empty datacube).
+.le
+.ls vecx = 1.0
+Rotation axis X vector.  Part of the specification of a three-dimensional
+orientation vector around which the datacube will appear to rotate when
+viewed from the front.  PROTOTYPE only supports rotations around the x axis.
+.le
+.ls vecy, vecz = 0.0
+Rotation axis Y and Z vectors.  In prototype, must be zero.
+.le
+.ls title = ""
+Output datacube title for rotation sequence.
+.le
+.ls maxws = 2000000
+Maximum workingset size in chars (usually 2 bytes).  Decrease if machine
+performance degrades noticeably during a run.  Increase if the machine has
+lots of memory and PVOL does not affect other processes.
+.le
+.ls abs = no
+If yes, take absolute value of voxel before applying any transformation.
+.le
+.ls verbose = yes
+Report memory usage, progress around the rotation, and more detail on
+errors if yes.
+.le
+
+
+.ih
+DESCRIPTION
+
+PVOL is used for visualizing the interiors of three-dimensional images.
+Opacity and intensity information is used to construct projected 2d images
+approximating an "xray" view through the original "solid", with varying
+amounts of apparent translucency.  Playing the resulting 2d images back
+rapidly as a filmloop generates the impression of a rotating translucent
+datacube inside of which you can view much of the original information with
+the illusion of seeing it in 3 dimensions.
+
+Given an input datacube plus rotation and projection parameters, PVOL
+produces a series of projected 2d images written out as another datacube.
+Rotation parameters control the number of frames to project, their
+angular separation, and the 3 vectors comprising the axis of rotation.
+In the prototype, only one rotation axis is allowed, counterclockwise
+about the X-axis when viewed facing the origin from +X (however, the user
+is viewing the datacube from -Z, and so sees the datacube rotating toward
+him/her).  When off-axis rotations are added, the view angle will still be
+from the front of the datacube.
+Non-orthogonal rotations in the prototype will have to be accomplished by
+first rotating the input datacube appropriately with other tools.
+
+Projection parameters
+provide control over the appearance of the projected images.  They may be
+tuned to visually enhance the apparent placement of interior regions in three
+dimensions during the rotation sequence.  Frames from the output datacube
+may be viewed individually on standard image display devices, may be
+played back rapidly with filmloop tools, or may be recorded to video as
+smooth, rotating volumes.  [At present the only filmloop tool available to us
+is MOVIE on Sun workstations, which requires preprocessing the datacube
+output from this task with another task called I2SUN].
+
+Sequences where the volume's rotation axis is the same as the viewing or
+projection axis are little more useful than a block average of the datacube,
+as hidden regions never rotate into view.  Volume rotations about the cube's
+X-axis (viewed from the front, or -Z) are the fastest and the only type
+implemented in the prototype.
+
+The \fBptype\fR parameter provides control over the type of projection.
+There are three main types of projection:  opacity, intensity, and both
+together.  If the
+input datacube is 4-dimensional, with two elements in the 4th dimension,
+both opacity and intensity information will be used -- first the remaining
+light along the projection will be attenuated by the opacity function, then
+the new voxel's intensity contribution added, according to \fBptype\fR.  Before
+the projection function is applied, the raw voxel intensity or opacity is
+clipped and scaled by transformation functions under control of task 
+parameters.
+.PP
+The image MIN and MAX must be present in the input image header, or they
+will default to 0.0 and 1.0 and a warning will be issued (run IMAGES.MINMAX
+with \fBupdate\fR=yes to set them if not already present).
+If intensity information is being used, \fBimin\fR and \fBimax\fR
+must be specified, or they will default to the image min and max.
+First we consider the intensity/opacity transformation functions, then we
+discuss how the transformed value contributes to the final projected image.
+
+.nf
+	Intensity transformation:
+
+	if (voxval < imin)
+	    newval = imin
+	else if (imin <= voxval && voxval < imax)
+	    newval = im_min + (im_max-im_min) * (voxval-imin)/(imax-imin)
+	else
+	    newval = imax
+	
+	Opacity transformation (0.0 <= attenuation <= 1.0):
+	if (voxval < omin)	# let maximum amount of light through
+	    attenuation = amax
+	else if (omin <= voxval && voxval < omax)
+	    attenuation = amin + (amax-amin) * (voxval*oscale - omin) /
+		(omax-omin)
+	else			# let minimum amount of light through
+	    attenuation = amin
+
+.fi
+
+The intensity class of projections includes \fBptype\fR = 2, 3, 4, 5, and 6.
+The default, \fBptype\fR 2, results in the AVERAGE transformed intensity along
+the projection column, while type 3 yields the SUM of transformed intensities.
+
+Type 4, PROPORTIONAL DISTANCE WEIGHTING, is used in conjunction with the 
+\fBdispower\fR parameter to weight the transformed voxel intensities by
+their inverse proportional depth along the projection column.
+If \fBdiscutoff\fR is no, the default, the distance will be that portion of
+the datacube intersected by the projection ray, measured starting at the
+rear (far side from the projection plane).  If \fBdiscutoff\fR is yes,
+the distance will be measured between the first and last voxels that fell
+between the cutoffs \fBimin\fR and \fBimax\fR.
+This projection generates a kind
+of depth cueing often useful in determining visually during filmloop playback
+which portions of the rotating image are in the foreground and which in the
+background (and how far).  The distance weighting is accomplished as follows,
+where voxposition and totvoxels are determined according to \fBdiscutoff\fR:
+
+.nf
+	\fBptype\fR = 4 (distance weighting):
+	newval = newval * (voxposition / voxelsincolumn) ** \fBdispower\fR
+.fi
+
+\fBptype\fR = 5, MODULAR PROPORTIONAL DISTANCE WEIGHTING, is useful for better
+seeing into the interiors of high-contrast datacubes.  Rather than using each
+voxel value along the projection column, only certain voxel values contribute,
+based on the \fBmodn\fR parameter (sometimes it is necessary to artificially
+"thin out" the data to see far enough into or through it).
+
+.nf
+	\fBptype\fR = 5 (modular distance weighting):
+	if (mod (int (newval/val_range * 100)) = \fBmodn\fR)
+	    use newval as in normal distance weighting
+	else
+	    ignore newval
+.fi
+
+\fBptype\fR = 6 results in only the LAST transformed voxel intensity that
+is between the \fBimin\fR and \fBimax\fR cutoffs being used.  This corresponds
+to seeing only the outer surface of datacube interior regions between the
+cutoffs (though since not every projection ray will pass through voxels
+right on the cutoff boundary, this will not necessarily result in a three
+dimensional intensity contour of an interior object; i.e. the intensities
+of those outer voxels can vary).
+
+OPACITY information can be used in viewing the interiors of 3d images, unlike
+in 2d images.  For \fBptype=1\fR parallel rays of light may be pictured
+shining through the datacube toward the projection plane, along the normal
+to that plane.  The voxel values in this
+case are considered to represent a degree of opacity, and a column of light
+will be attenuated by each voxel according to a function of its opacity value
+as the ray proceeds through the volume.  The \fBizero\fR parameter provides
+the initial incident "light" intensity before any attenuation.  The
+amount of remaining light after projection through the datacube is very
+sensitive to the voxel opacities and the number of voxels in each projection
+column.  Consequently, the \fBoscale\fR parameter is supplied to enable
+adjusting the relative attenuation in a single step while scouting for
+the right opacity transformation function to generate the desired effect
+during playback rotation.  Given the amount of attenuation
+as determined in the opacity transformation function above, for each 
+contributing voxel along the projection column:
+
+.nf
+	projection pixel = projection pixel * attenuation
+.fi
+
+If the input image is 4-dimensional, with 2 elements in the 4th dimension,
+voxel intensities will be added after attenuation 
+to contribute to the total projected pixel value (like a cloud
+with both absorption and emission).  For
+purposes of visualization only, it is not necessary that the voxel value
+represent a physically real opacity; any data value may be treated as
+attenuating an imaginary xray passing through the solid in order to help
+image the volume in three apparent dimensions.
+
+For all of the projection types, once the modified intensity
+has been determined, it contributes to the output pixel onto which the
+current, arbitrarily-oriented column of voxels projects.  To summarize:
+
+.nf
+	1 OPACITY:
+	    proj_pix = proj_pix * attenuation
+	2 AVERAGE:
+	    proj_pix = proj_pix + newval / nvox
+	3 SUM:
+	    proj_pix = proj_pix + newval
+	4 INVDISPOW:
+	    proj_pix = proj_pix + newval * (vox/voxincol)**dispow
+	5 MOD:
+	    if mod (int (newval/val_range * 100.0)) = \fBmodn\fR
+		proj_pix = proj_pix + newval * (vox/voxincol)**dispow
+	6 LASTONLY:
+	    if (\fBimin\fR < newval && newval <= \fBimax\fR)
+		proj_pix = newval
+
+.fi
+
+.ih
+PERFORMANCE AND SIZE CONSTRAINTS
+
+Projections through 3d images inherently require large amounts of memory,
+or else the tasks will spend all their time thrashing with I/O.  In volume
+rotations about the X-axis, each output pixel is derived by projecting at
+an arbitrary angle through a YZ slice of the input image.  Because of otherwise
+excessive thrashing, PVOL requires sufficient memory for at least one YZ
+slice.  The more YZ slices that will fit into memory at one time, the better,
+because I/O is more efficient the larger the chunk of the image that can
+be read at one time.  It is best if the entire image will fit into memory,
+as the output image (all rotations) will not have to be reread for each
+successive chunk of YZ slices.  Available memory is that actually allocable
+by PVOL for the slices plus one line of the output image.  On a workstation
+there will usually be considerably less memory available for PVOL than
+the amount physically in the machine if running in a window environment.
+Examples of the number of YZ slices that will fit based on image size and
+available memory follow; image datatype is assumed to be REAL -- multiply
+number of YZ slices by 2 for SHORT images.
+
+.nf
+	Usable Memory	Image Size	Approx YZ Slices
+	------------------------------------------------
+	1 Mb		64*64*64	64 (whole image)
+	1 Mb		512*512*512	1
+	4 Mb		101*101*101	101 (whole image)
+	4 Mb		1024*1024*1024	1
+	8 Mb		128*128*128	128 (whole image)
+	8 Mb		1448*1448*1448	1
+	16 Mb		161*161*161	161 (whole image)
+	16 Mb		2048*2048*2048	1
+	32 Mb		203*203*203	203 (whole image)
+	32 Mb		2896*2896*2896	1
+	64 Mb		256*256*256	256 (whole image)
+	128 Mb		322*322*322	322 (whole image)
+	512 Mb		512*512*512	512 (whole image)
+.fi
+
+PVOL checks to see how much memory it can grab, then actually allocates
+somewhat less than this (otherwise you wouldn't be able to do anything 
+except run IRAF tasks already loaded in the process cache until PVOL
+finishes).  With \fBverbose\fR on, the task reports memory usage figures.  
+On some machines the system will continue to allocate more memory for a
+task even above that reported by PVOL.  This can be a problem if you fire
+up PVOL from a workstation (even with lots of windows already open);
+after you log out, the system may grab that extra memory you were using,
+and not even let you back in later.  This is why the \fBmaxws\fR
+parameter is supplied -- lower it if this type of behavior is experienced.
+
+.ih
+EXAMPLES
+
+.nf
+1.  Produce 36 rotation projections (one every 10 degrees) around the
+    x-axis of a datacube, viewed from the front (negative z
+    direction).  Assume that the single-valued input voxel values
+    are intensities, and that the image header contains MIN and MAX.
+
+    cl> pvol input output
+
+2.  Generate 180 frames, one every two degrees.
+
+    cl> pvol input output nframes=180 degrees=2
+
+3.  Use inverse proportional distance cubed weighting in two
+    subsampled projections for a quick look.  Distance-weight
+    only between projection voxels falling within the specified
+    cutoffs (0.1 to 1.0).
+
+    cl> pvol input[*:4,*:4,*:4] output nfr=2 deg=90 ptype=4 \
+	dispower=3 discutoff+ imin=.1 imax=1.0
+
+4.  Project through a 4d image containing opacity information in
+    element 2 of the 4th axis and intensity in element 1.  Scale
+    the voxel opacities by 0.1 to allow more light through.  Use
+    the SUM of the voxel intensity values (which will be attenuated
+    by subsequent opacities), with no distance weighting.
+
+    cl> pvol input output ptype=3 opacelem=2
+
+.fi
+
+.ih
+TIMINGS
+
+1min 12sec cpu on an unloaded Sun-4 to produce
+36 rotation increments around a 50*50*50 datacube with \fBptype\fR=2
+(uses less than 1 Mb of memory for image data); 46sec for \fBptype\fR=1;
+2min 19sec for \fBptype\fR=4.
+
+4min 32sec cpu on an unloaded Sun-3 with 8 Mb memory to do 36 steps around a
+50*50*50 datacube with \fBptype\fR=2 (also uses less than 1 Mb);
+3min 20sec for \fBptype\fR=1; 10min 51sec for \fBptype\fR=4.
+
+17hr 20 min cpu on a Sun-4 to do 36 rotation steps around a 450*450*450
+datacube with \fBptype\fR=4.
+
+.ih
+BUGS
+
+Maximizing memory usage without adversely impacting other functions can be
+tricky.  Adverse effects may result from using too high a \fBmaxws\fR.
+
+Cannot rotate around arbitrary axis yet.
+
+Lacks shading algorithm.
+
+Needs easier user interface to adjust translucency parameters (e.g. with
+mouse when workstations become fast enough to do this in real time).
+
+.ih
+SEE ALSO
+i2sun, im3dtran, im3dstack
+.endhelp
diff --git a/pkg/proto/vol/src/doc/volumes.hlp b/pkg/proto/vol/src/doc/volumes.hlp
new file mode 100644
index 00000000..4ebb4aeb
--- /dev/null
+++ b/pkg/proto/vol/src/doc/volumes.hlp
@@ -0,0 +1,56 @@
+.help volumes Jan89 "Volumes Package"
+
+***** NOTE:  This is just a suggested package organization and will
+*****        definitely NOT be the final one chosen.
+
+.ce
+Volume or 3d Image Applications in IRAF
+.ce
+January 1989
+
+.sh
+Introduction
+
+The Volumes package collects tasks related to manipulating and displaying
+volume images (3d images, or datacubes).  Although all IRAF images can be
+multidimensional (currently up to 7 dimensions), not all applications tasks
+are equipped to handle images of dimension greater than 2.  Examples of
+tasks that are so equipped are IMARITH, IMSTATISTICS, BLKAVG, and DISPLAY
+for looking at arbitrary 2d sections of higher dimensional images.
+
+Volumes applications include tasks for manipulating the orientation of
+a 3d image, joining 3d images, projections of datacube contents
+onto 2d images, and tasks related to viewing a datacube or its projections
+as a movie.
+
+.ih
+Datacube Manipulation Tasks
+
+D3TRANSPOSE	3d transpose, any axis to any other axis
+IMJOIN		join 2 or more 3d images together along specified axis
+IMCOPY
+BLKAVG
+IMSLICE
+
+.ih
+Datacube Generation Tasks
+
+BINTOIM		[not in VOLUMES; probably still PROTO after upgrade to 3d?]
+POINTOIM	convert n-dimensional point data into volumes in datacube
+MANDEL4		4d Mandelbrot set generator
+
+.ih
+Volume Projection Tasks
+
+PVOL		project volume contents onto series of 2d images
+SLICEVOL*	"cubetool" -- slice off faces of datacube rendered from
+		    arbitrary angle w/translucency
+
+.ih
+Movie-Related Tasks
+
+IMTOSUN		convert datacube or list of 2d images into Sun rasterfiles
+IMTOVID		(script) record set of 2d images onto panasonic video recorder
+CUBETOVID	(script) record sliced from databube onto video recorder
+
+* = [not yet implemented] 
diff --git a/pkg/proto/vol/src/i2sun.par b/pkg/proto/vol/src/i2sun.par
new file mode 100644
index 00000000..d28d887c
--- /dev/null
+++ b/pkg/proto/vol/src/i2sun.par
@@ -0,0 +1,14 @@
+input,s,a,,,,Input image template or 3d image
+output,s,a,,,,Output rasterfile template
+z1,r,a,,,,Minimum greylevel to be displayed
+z2,r,a,,,,Maximum greylevel to be displayed
+clutfile,f,h,"",,,Input rasterfile containing color lookup table
+ztrans,s,h,linear,,,Greylevel transformation (linear|log|none|user)
+ulutfile,f,h,"",,,File containing user defined look up table
+xsize,i,h,INDEF,1,,Output rasterfile horizontal size
+ysize,i,h,INDEF,1,,Output rasterfile vertical size
+xmag,r,h,1.,,,Output rasterfile horizontal magnification
+ymag,r,h,1.,,,Output rasterfile vertical magnification
+order,i,h,1,0,1,"Spatial interpolator order; 0=replic., 1=linear"
+sliceaxis,i,h,3,,,"Slice a 3d or higher image through this axis"
+swap,b,h,no,,,"Swap bytes in output rasterfiles?"
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
diff --git a/pkg/proto/vol/src/im3dtran.par b/pkg/proto/vol/src/im3dtran.par
new file mode 100644
index 00000000..3c24953b
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran.par
@@ -0,0 +1,6 @@
+input,s,a,,,,Input 3d image (datacube)
+output,s,a,,,,Output 3d image
+new_x,i,a,3,1,3,"New x axis = old axis (1=x, 2=y, 3=z)"
+new_y,i,a,2,1,3,"New y axis = old axis (1=x, 2=y, 3=z)"
+new_z,i,a,1,1,3,"New z axis = old axis (1=x, 2=y, 3=z)"
+len_blk,i,h,128,,,Size in pixels of internal subraster
diff --git a/pkg/proto/vol/src/im3dtran/mkpkg b/pkg/proto/vol/src/im3dtran/mkpkg
new file mode 100644
index 00000000..19b0b52d
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/mkpkg
@@ -0,0 +1,52 @@
+# Library for the 3DTRANSPOSE task.
+
+$checkout libpkg.a /u2/rooke/vol/
+$update   libpkg.a
+$checkin  libpkg.a /u2/rooke/vol/
+$exit
+
+tfiles:
+	$ifolder (txyz3.x, txyz3.gx)
+	    $generic -k txyz3.gx -o txyz3.x		$endif
+	$ifolder (txzy3.x, txzy3.gx)
+	    $generic -k txzy3.gx -o txzy3.x		$endif
+	$ifolder (tyxz3.x, tyxz3.gx)
+	    $generic -k tyxz3.gx -o tyxz3.x		$endif
+	$ifolder (tyzx3.x, tyzx3.gx)
+	    $generic -k tyzx3.gx -o tyzx3.x		$endif
+	$ifolder (tzxy3.x, tzxy3.gx)
+	    $generic -k tzxy3.gx -o tzxy3.x		$endif
+	$ifolder (tzyx3.x, tzyx3.gx)
+	    $generic -k tzyx3.gx -o tzyx3.x		$endif
+	;
+
+libpkg.a:
+	$ifeq (USE_GENERIC, yes) $call tfiles $endif
+
+	t_im3dtran.x	<imhdr.h>
+	txyz3.x
+	txzy3.x
+	tyxz3.x
+	tyzx3.x
+	tzxy3.x
+	tzyx3.x
+	;
+
+dbx:
+	$set	XFLAGS = "-c -g -F -q"
+	$set	LFLAGS = "-g -q"
+	$set	LIBS = "-lxtools"
+
+	$ifeq (USE_GENERIC, yes) $call tfiles $endif
+
+	$omake	x_im3dtran.x
+	$omake	t_im3dtran.x
+	$omake	txyz3.x
+	$omake	txzy3.x
+	$omake	tyxz3.x
+	$omake	tyzx3.x
+	$omake	tzxy3.x
+	$omake	tzyx3.x
+	$link	x_im3dtran.o t_im3dtran.o txyz3.o txzy3.o tyxz3.o tyzx3.o \
+		tzxy3.o tzyx3.o $(LIBS) -o xx_im3dtran.e
+	;
diff --git a/pkg/proto/vol/src/im3dtran/t_im3dtran.x b/pkg/proto/vol/src/im3dtran/t_im3dtran.x
new file mode 100644
index 00000000..a77c0703
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/t_im3dtran.x
@@ -0,0 +1,307 @@
+include	<imhdr.h>
+include	<error.h>
+
+define	XYZ	1	# xyz -> xyz (identity)
+define	XZY	2	# xyz -> xzy
+define	YXZ	3	# xyz -> yxz
+define	YZX	4	# xyz -> yzx
+define	ZXY	5	# xyz -> zxy
+define	ZYX	6	# xyz -> zyx
+
+
+# T_IM3DTRAN -- Transpose 3d images.
+#
+# The input and output images are given by image template lists.  The
+# number of output images must match the number of input images.  Image
+# sections are allowed in the input images and are ignored in the output
+# images.  If the input and output image names are the same then the transpose
+# is performed to a temporary file which then replaces the input image.
+
+procedure t_im3dtran ()
+
+char	imtlist1[SZ_LINE]			# Input image list
+char	imtlist2[SZ_LINE]			# Output image list
+int	len_blk					# 1D length of transpose block
+
+char	image1[SZ_FNAME]			# Input image name
+char	image2[SZ_FNAME]			# Output image name
+char	imtemp[SZ_FNAME]			# Temporary file
+
+int	list1, list2, new_ax[3], which3d
+pointer	im1, im2
+
+int	clgeti(), imtopen(), imtgetim(), imtlen(), whichtran()
+pointer	immap()
+
+begin
+	# Get input and output image template lists, the size of the transpose
+	# block, and the transpose mapping.
+
+	call clgstr ("input", imtlist1, SZ_LINE)
+	call clgstr ("output", imtlist2, SZ_LINE)
+	len_blk = clgeti ("len_blk")
+	new_ax[1] = clgeti ("new_x")
+	new_ax[2] = clgeti ("new_y")
+	new_ax[3] = clgeti ("new_z")
+
+	# Determine the type of 3d transpose.
+	which3d = whichtran (new_ax)
+	if (which3d <= 0)
+	    call error (0, "Invalid mapping of new_x, new_y, new_z")
+
+	# Expand the input and output image lists.
+
+	list1 = imtopen (imtlist1)
+	list2 = imtopen (imtlist2)
+	if (imtlen (list1) != imtlen (list2)) {
+	    call imtclose (list1)
+	    call imtclose (list2)
+	    call error (1, "Number of input and output images not the same")
+	}
+
+	# Do each set of input/output images.
+
+	while ((imtgetim (list1, image1, SZ_FNAME) != EOF) &&
+	    (imtgetim (list2, image2, SZ_FNAME) != EOF)) {
+
+	    call xt_mkimtemp (image1, image2, imtemp, SZ_FNAME)
+
+	    im1 = immap (image1, READ_ONLY, 0)
+	    im2 = immap (image2, NEW_COPY, im1)
+
+	    # Do the transpose.
+	    call im3dtranspose (im1, im2, len_blk, which3d, new_ax)
+
+	    # Unmap the input and output images.
+	    call imunmap (im1)
+	    call imunmap (im2)
+
+	    call xt_delimtemp (image2, imtemp)
+	}
+
+	call imtclose (list1)
+	call imtclose (list2)
+end
+
+
+# IM3DTRANSPOSE -- Transpose an image. 
+#
+# Divide the image into square blocks of size len_blk by len_blk.
+# Transpose each block with a generic array transpose operator.
+
+procedure im3dtranspose (im_in, im_out, len_blk, which3d, new_ax)
+
+pointer	im_in				# Input image descriptor
+pointer	im_out				# Output image descriptor
+int	len_blk				# 1D length of transpose block
+int	which3d				# Parameterized transpose order
+int	new_ax[3]			# Map old axis[index] to new value
+
+int	x1, x2, nx
+int	y1, y2, ny
+int	z1, z2, nz
+pointer	buf_in, buf_out
+
+pointer	imgs3s(), imps3s(), imgs3i(), imps3i(), imgs3l(), imps3l()
+pointer	imgs3r(), imps3r(), imgs3d(), imps3d(), imgs3x(), imps3x()
+
+begin
+	# Output image is a copy of input image with dims transposed.
+
+	IM_LEN (im_out, 1) = IM_LEN (im_in, new_ax[1])
+	IM_LEN (im_out, 2) = IM_LEN (im_in, new_ax[2])
+	IM_LEN (im_out, 3) = IM_LEN (im_in, new_ax[3])
+
+	# Break the input image into blocks of at most (len_blk)**3 .
+
+	do x1 = 1, IM_LEN (im_in, 1), len_blk {
+	    x2 = x1 + len_blk - 1
+	    if (x2 > IM_LEN(im_in, 1))
+	       x2 = IM_LEN(im_in, 1)
+	    nx = x2 - x1 + 1
+
+	    do y1 = 1, IM_LEN (im_in, 2), len_blk {
+	        y2 = y1 + len_blk - 1
+	        if (y2 > IM_LEN(im_in, 2))
+		   y2 = IM_LEN(im_in, 2)
+		ny = y2 - y1 + 1
+
+		do z1 = 1, IM_LEN (im_in, 3), len_blk {
+		    z2 = z1 + len_blk - 1
+		    if (z2 > IM_LEN(im_in, 3))
+			z2 = IM_LEN(im_in, 3)
+		    nz = z2 - z1 + 1
+
+		    # Switch on the pixel type to optimize IMIO.
+
+		    switch (IM_PIXTYPE (im_in)) {
+		    case TY_SHORT:
+			buf_in = imgs3s (im_in, x1, x2, y1, y2, z1, z2)
+			switch (which3d) {
+			case XYZ:
+			    buf_out = imps3s (im_out, x1, x2, y1, y2, z1, z2)
+			    call txyz3s (Mems[buf_in], Mems[buf_out], nx,ny,nz)
+			case XZY:
+			    buf_out = imps3s (im_out, x1, x2, z1, z2, y1, y2)
+			    call txzy3s (Mems[buf_in], Mems[buf_out], nx,ny,nz)
+			case YXZ:
+			    buf_out = imps3s (im_out, y1, y2, x1, x2, z1, z2)
+			    call tyxz3s (Mems[buf_in], Mems[buf_out], nx,ny,nz)
+			case YZX:
+			    buf_out = imps3s (im_out, y1, y2, z1, z2, x1, x2)
+			    call tyzx3s (Mems[buf_in], Mems[buf_out], nx,ny,nz)
+			case ZXY:
+			    buf_out = imps3s (im_out, z1, z2, x1, x2, y1, y2)
+			    call tzxy3s (Mems[buf_in], Mems[buf_out], nx,ny,nz)
+			case ZYX:
+			    buf_out = imps3s (im_out, z1, z2, y1, y2, x1, x2)
+			    call tzyx3s (Mems[buf_in], Mems[buf_out], nx,ny,nz)
+			}
+		    case TY_INT:
+			buf_in = imgs3i (im_in, x1, x2, y1, y2, z1, z2)
+			switch (which3d) {
+			case XYZ:
+			    buf_out = imps3i (im_out, x1, x2, y1, y2, z1, z2)
+			    call txyz3i (Memi[buf_in], Memi[buf_out], nx,ny,nz)
+			case XZY:
+			    buf_out = imps3i (im_out, x1, x2, z1, z2, y1, y2)
+			    call txzy3i (Memi[buf_in], Memi[buf_out], nx,ny,nz)
+			case YXZ:
+			    buf_out = imps3i (im_out, y1, y2, x1, x2, z1, z2)
+			    call tyxz3i (Memi[buf_in], Memi[buf_out], nx,ny,nz)
+			case YZX:
+			    buf_out = imps3i (im_out, y1, y2, z1, z2, x1, x2)
+			    call tyzx3i (Memi[buf_in], Memi[buf_out], nx,ny,nz)
+			case ZXY:
+			    buf_out = imps3i (im_out, z1, z2, x1, x2, y1, y2)
+			    call tzxy3i (Memi[buf_in], Memi[buf_out], nx,ny,nz)
+			case ZYX:
+			    buf_out = imps3i (im_out, z1, z2, y1, y2, x1, x2)
+			    call tzyx3i (Memi[buf_in], Memi[buf_out], nx,ny,nz)
+			}
+		    case TY_LONG:
+			buf_in = imgs3l (im_in, x1, x2, y1, y2, z1, z2)
+			switch (which3d) {
+			case XYZ:
+			    buf_out = imps3l (im_out, x1, x2, y1, y2, z1, z2)
+			    call txyz3l (Meml[buf_in], Meml[buf_out], nx,ny,nz)
+			case XZY:
+			    buf_out = imps3l (im_out, x1, x2, z1, z2, y1, y2)
+			    call txzy3l (Meml[buf_in], Meml[buf_out], nx,ny,nz)
+			case YXZ:
+			    buf_out = imps3l (im_out, y1, y2, x1, x2, z1, z2)
+			    call tyxz3l (Meml[buf_in], Meml[buf_out], nx,ny,nz)
+			case YZX:
+			    buf_out = imps3l (im_out, y1, y2, z1, z2, x1, x2)
+			    call tyzx3l (Meml[buf_in], Meml[buf_out], nx,ny,nz)
+			case ZXY:
+			    buf_out = imps3l (im_out, z1, z2, x1, x2, y1, y2)
+			    call tzxy3l (Meml[buf_in], Meml[buf_out], nx,ny,nz)
+			case ZYX:
+			    buf_out = imps3l (im_out, z1, z2, y1, y2, x1, x2)
+			    call tzyx3l (Meml[buf_in], Meml[buf_out], nx,ny,nz)
+			}
+		    case TY_REAL:
+			buf_in = imgs3r (im_in, x1, x2, y1, y2, z1, z2)
+			switch (which3d) {
+			case XYZ:
+			    buf_out = imps3r (im_out, x1, x2, y1, y2, z1, z2)
+			    call txyz3r (Memr[buf_in], Memr[buf_out], nx,ny,nz)
+			case XZY:
+			    buf_out = imps3r (im_out, x1, x2, z1, z2, y1, y2)
+			    call txzy3r (Memr[buf_in], Memr[buf_out], nx,ny,nz)
+			case YXZ:
+			    buf_out = imps3r (im_out, y1, y2, x1, x2, z1, z2)
+			    call tyxz3r (Memr[buf_in], Memr[buf_out], nx,ny,nz)
+			case YZX:
+			    buf_out = imps3r (im_out, y1, y2, z1, z2, x1, x2)
+			    call tyzx3r (Memr[buf_in], Memr[buf_out], nx,ny,nz)
+			case ZXY:
+			    buf_out = imps3r (im_out, z1, z2, x1, x2, y1, y2)
+			    call tzxy3r (Memr[buf_in], Memr[buf_out], nx,ny,nz)
+			case ZYX:
+			    buf_out = imps3r (im_out, z1, z2, y1, y2, x1, x2)
+			    call tzyx3r (Memr[buf_in], Memr[buf_out], nx,ny,nz)
+			}
+		    case TY_DOUBLE:
+			buf_in = imgs3d (im_in, x1, x2, y1, y2, z1, z2)
+			switch (which3d) {
+			case XYZ:
+			    buf_out = imps3d (im_out, x1, x2, y1, y2, z1, z2)
+			    call txyz3d (Memd[buf_in], Memd[buf_out], nx,ny,nz)
+			case XZY:
+			    buf_out = imps3d (im_out, x1, x2, z1, z2, y1, y2)
+			    call txzy3d (Memd[buf_in], Memd[buf_out], nx,ny,nz)
+			case YXZ:
+			    buf_out = imps3d (im_out, y1, y2, x1, x2, z1, z2)
+			    call tyxz3d (Memd[buf_in], Memd[buf_out], nx,ny,nz)
+			case YZX:
+			    buf_out = imps3d (im_out, y1, y2, z1, z2, x1, x2)
+			    call tyzx3d (Memd[buf_in], Memd[buf_out], nx,ny,nz)
+			case ZXY:
+			    buf_out = imps3d (im_out, z1, z2, x1, x2, y1, y2)
+			    call tzxy3d (Memd[buf_in], Memd[buf_out], nx,ny,nz)
+			case ZYX:
+			    buf_out = imps3d (im_out, z1, z2, y1, y2, x1, x2)
+			    call tzyx3d (Memd[buf_in], Memd[buf_out], nx,ny,nz)
+			}
+		    case TY_COMPLEX:
+			buf_in = imgs3x (im_in, x1, x2, y1, y2, z1, z2)
+			switch (which3d) {
+			case XYZ:
+			    buf_out = imps3x (im_out, x1, x2, y1, y2, z1, z2)
+			    call txyz3x (Memx[buf_in], Memx[buf_out], nx,ny,nz)
+			case XZY:
+			    buf_out = imps3x (im_out, x1, x2, z1, z2, y1, y2)
+			    call txzy3x (Memx[buf_in], Memx[buf_out], nx,ny,nz)
+			case YXZ:
+			    buf_out = imps3x (im_out, y1, y2, x1, x2, z1, z2)
+			    call tyxz3x (Memx[buf_in], Memx[buf_out], nx,ny,nz)
+			case YZX:
+			    buf_out = imps3x (im_out, y1, y2, z1, z2, x1, x2)
+			    call tyzx3x (Memx[buf_in], Memx[buf_out], nx,ny,nz)
+			case ZXY:
+			    buf_out = imps3x (im_out, z1, z2, x1, x2, y1, y2)
+			    call tzxy3x (Memx[buf_in], Memx[buf_out], nx,ny,nz)
+			case ZYX:
+			    buf_out = imps3x (im_out, z1, z2, y1, y2, x1, x2)
+			    call tzyx3x (Memx[buf_in], Memx[buf_out], nx,ny,nz)
+			}
+		    default:
+			call error (3, "unknown pixel type")
+		    }
+		}
+	    }
+	}
+end
+
+
+# WHICHTRAN -- Return transpose type.
+
+int procedure whichtran (new_ax)
+int	new_ax[3]
+
+int	which
+
+begin
+	which = 0
+
+	if (new_ax[1] == 1) {
+	    if (new_ax[2] == 2)
+		which = XYZ
+	    else if (new_ax[2] == 3)
+		which = XZY
+	} else if (new_ax[1] == 2) {
+	    if (new_ax[2] == 1)
+		which = YXZ
+	    else if (new_ax[2] == 3)
+		which = YZX
+	} else if (new_ax[1] == 3) {
+	    if (new_ax[2] == 1)
+		which = ZXY
+	    else if (new_ax[2] == 2)
+		which = ZYX
+	}
+	
+	return (which)
+end
diff --git a/pkg/proto/vol/src/im3dtran/txyz3.gx b/pkg/proto/vol/src/im3dtran/txyz3.gx
new file mode 100644
index 00000000..619734a1
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/txyz3.gx
@@ -0,0 +1,18 @@
+$for (silrdx)
+
+# TXYZ3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txyz3$t (a, b, nx, ny, nz)
+
+PIXEL	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, z, y]
+end
+
+$endfor
diff --git a/pkg/proto/vol/src/im3dtran/txyz3.x b/pkg/proto/vol/src/im3dtran/txyz3.x
new file mode 100644
index 00000000..1cc8ca92
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/txyz3.x
@@ -0,0 +1,103 @@
+
+
+# TXYZ3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txyz3s (a, b, nx, ny, nz)
+
+short	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, z, y]
+end
+
+
+
+# TXYZ3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txyz3i (a, b, nx, ny, nz)
+
+int	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, z, y]
+end
+
+
+
+# TXYZ3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txyz3l (a, b, nx, ny, nz)
+
+long	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, z, y]
+end
+
+
+
+# TXYZ3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txyz3r (a, b, nx, ny, nz)
+
+real	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, z, y]
+end
+
+
+
+# TXYZ3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txyz3d (a, b, nx, ny, nz)
+
+double	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, z, y]
+end
+
+
+
+# TXYZ3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txyz3x (a, b, nx, ny, nz)
+
+complex	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, z, y]
+end
+
+
diff --git a/pkg/proto/vol/src/im3dtran/txzy3.gx b/pkg/proto/vol/src/im3dtran/txzy3.gx
new file mode 100644
index 00000000..a6d18e4a
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/txzy3.gx
@@ -0,0 +1,18 @@
+$for (silrdx)
+
+# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txzy3$t (a, b, nx, ny, nz)
+
+PIXEL	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, y, z]
+end
+
+$endfor
diff --git a/pkg/proto/vol/src/im3dtran/txzy3.x b/pkg/proto/vol/src/im3dtran/txzy3.x
new file mode 100644
index 00000000..ad6096bf
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/txzy3.x
@@ -0,0 +1,103 @@
+
+
+# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txzy3s (a, b, nx, ny, nz)
+
+short	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, y, z]
+end
+
+
+
+# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txzy3i (a, b, nx, ny, nz)
+
+int	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, y, z]
+end
+
+
+
+# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txzy3l (a, b, nx, ny, nz)
+
+long	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, y, z]
+end
+
+
+
+# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txzy3r (a, b, nx, ny, nz)
+
+real	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, y, z]
+end
+
+
+
+# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txzy3d (a, b, nx, ny, nz)
+
+double	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, y, z]
+end
+
+
+
+# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txzy3x (a, b, nx, ny, nz)
+
+complex	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, y, z]
+end
+
+
diff --git a/pkg/proto/vol/src/im3dtran/tyxz3.gx b/pkg/proto/vol/src/im3dtran/tyxz3.gx
new file mode 100644
index 00000000..75c2244f
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/tyxz3.gx
@@ -0,0 +1,18 @@
+$for (silrdx)
+
+# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z.  The arrays need not be
+# identical.
+
+procedure tyxz3$t (a, b, nx, ny, nz)
+
+PIXEL	a[nx, ny, nz], b[ny, nx, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, x, z] = a[x, y, z]
+end
+
+$endfor
diff --git a/pkg/proto/vol/src/im3dtran/tyxz3.x b/pkg/proto/vol/src/im3dtran/tyxz3.x
new file mode 100644
index 00000000..166ae8de
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/tyxz3.x
@@ -0,0 +1,103 @@
+
+
+# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z.  The arrays need not be
+# identical.
+
+procedure tyxz3s (a, b, nx, ny, nz)
+
+short	a[nx, ny, nz], b[ny, nx, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, x, z] = a[x, y, z]
+end
+
+
+
+# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z.  The arrays need not be
+# identical.
+
+procedure tyxz3i (a, b, nx, ny, nz)
+
+int	a[nx, ny, nz], b[ny, nx, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, x, z] = a[x, y, z]
+end
+
+
+
+# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z.  The arrays need not be
+# identical.
+
+procedure tyxz3l (a, b, nx, ny, nz)
+
+long	a[nx, ny, nz], b[ny, nx, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, x, z] = a[x, y, z]
+end
+
+
+
+# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z.  The arrays need not be
+# identical.
+
+procedure tyxz3r (a, b, nx, ny, nz)
+
+real	a[nx, ny, nz], b[ny, nx, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, x, z] = a[x, y, z]
+end
+
+
+
+# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z.  The arrays need not be
+# identical.
+
+procedure tyxz3d (a, b, nx, ny, nz)
+
+double	a[nx, ny, nz], b[ny, nx, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, x, z] = a[x, y, z]
+end
+
+
+
+# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z.  The arrays need not be
+# identical.
+
+procedure tyxz3x (a, b, nx, ny, nz)
+
+complex	a[nx, ny, nz], b[ny, nx, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, x, z] = a[x, y, z]
+end
+
+
diff --git a/pkg/proto/vol/src/im3dtran/tyzx3.gx b/pkg/proto/vol/src/im3dtran/tyzx3.gx
new file mode 100644
index 00000000..108910aa
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/tyzx3.gx
@@ -0,0 +1,18 @@
+$for (silrdx)
+
+# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x.  The arrays need not be
+# identical.
+
+procedure tyzx3$t (a, b, nx, ny, nz)
+
+PIXEL	a[nx, ny, nz], b[ny, nz, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, z, x] = a[x, y, z]
+end
+
+$endfor
diff --git a/pkg/proto/vol/src/im3dtran/tyzx3.x b/pkg/proto/vol/src/im3dtran/tyzx3.x
new file mode 100644
index 00000000..6b05e748
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/tyzx3.x
@@ -0,0 +1,103 @@
+
+
+# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x.  The arrays need not be
+# identical.
+
+procedure tyzx3s (a, b, nx, ny, nz)
+
+short	a[nx, ny, nz], b[ny, nz, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, z, x] = a[x, y, z]
+end
+
+
+
+# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x.  The arrays need not be
+# identical.
+
+procedure tyzx3i (a, b, nx, ny, nz)
+
+int	a[nx, ny, nz], b[ny, nz, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, z, x] = a[x, y, z]
+end
+
+
+
+# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x.  The arrays need not be
+# identical.
+
+procedure tyzx3l (a, b, nx, ny, nz)
+
+long	a[nx, ny, nz], b[ny, nz, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, z, x] = a[x, y, z]
+end
+
+
+
+# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x.  The arrays need not be
+# identical.
+
+procedure tyzx3r (a, b, nx, ny, nz)
+
+real	a[nx, ny, nz], b[ny, nz, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, z, x] = a[x, y, z]
+end
+
+
+
+# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x.  The arrays need not be
+# identical.
+
+procedure tyzx3d (a, b, nx, ny, nz)
+
+double	a[nx, ny, nz], b[ny, nz, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, z, x] = a[x, y, z]
+end
+
+
+
+# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x.  The arrays need not be
+# identical.
+
+procedure tyzx3x (a, b, nx, ny, nz)
+
+complex	a[nx, ny, nz], b[ny, nz, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, z, x] = a[x, y, z]
+end
+
+
diff --git a/pkg/proto/vol/src/im3dtran/tzxy3.gx b/pkg/proto/vol/src/im3dtran/tzxy3.gx
new file mode 100644
index 00000000..3fbed0b5
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/tzxy3.gx
@@ -0,0 +1,18 @@
+$for (silrdx)
+
+# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y.  The arrays need not be
+# identical.
+
+procedure tzxy3$t (a, b, nx, ny, nz)
+
+PIXEL	a[nx, ny, nz], b[nz, nx, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, x, y] = a[x, y, z]
+end
+
+$endfor
diff --git a/pkg/proto/vol/src/im3dtran/tzxy3.x b/pkg/proto/vol/src/im3dtran/tzxy3.x
new file mode 100644
index 00000000..d3b30f31
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/tzxy3.x
@@ -0,0 +1,103 @@
+
+
+# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y.  The arrays need not be
+# identical.
+
+procedure tzxy3s (a, b, nx, ny, nz)
+
+short	a[nx, ny, nz], b[nz, nx, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, x, y] = a[x, y, z]
+end
+
+
+
+# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y.  The arrays need not be
+# identical.
+
+procedure tzxy3i (a, b, nx, ny, nz)
+
+int	a[nx, ny, nz], b[nz, nx, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, x, y] = a[x, y, z]
+end
+
+
+
+# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y.  The arrays need not be
+# identical.
+
+procedure tzxy3l (a, b, nx, ny, nz)
+
+long	a[nx, ny, nz], b[nz, nx, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, x, y] = a[x, y, z]
+end
+
+
+
+# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y.  The arrays need not be
+# identical.
+
+procedure tzxy3r (a, b, nx, ny, nz)
+
+real	a[nx, ny, nz], b[nz, nx, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, x, y] = a[x, y, z]
+end
+
+
+
+# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y.  The arrays need not be
+# identical.
+
+procedure tzxy3d (a, b, nx, ny, nz)
+
+double	a[nx, ny, nz], b[nz, nx, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, x, y] = a[x, y, z]
+end
+
+
+
+# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y.  The arrays need not be
+# identical.
+
+procedure tzxy3x (a, b, nx, ny, nz)
+
+complex	a[nx, ny, nz], b[nz, nx, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, x, y] = a[x, y, z]
+end
+
+
diff --git a/pkg/proto/vol/src/im3dtran/tzyx3.gx b/pkg/proto/vol/src/im3dtran/tzyx3.gx
new file mode 100644
index 00000000..61d32e6d
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/tzyx3.gx
@@ -0,0 +1,18 @@
+$for (silrdx)
+
+# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x.  The arrays need not be
+# identical.
+
+procedure tzyx3$t (a, b, nx, ny, nz)
+
+PIXEL	a[nx, ny, nz], b[nz, ny, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, y, x] = a[x, y, z]
+end
+
+$endfor
diff --git a/pkg/proto/vol/src/im3dtran/tzyx3.x b/pkg/proto/vol/src/im3dtran/tzyx3.x
new file mode 100644
index 00000000..8cc4c877
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/tzyx3.x
@@ -0,0 +1,103 @@
+
+
+# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x.  The arrays need not be
+# identical.
+
+procedure tzyx3s (a, b, nx, ny, nz)
+
+short	a[nx, ny, nz], b[nz, ny, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, y, x] = a[x, y, z]
+end
+
+
+
+# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x.  The arrays need not be
+# identical.
+
+procedure tzyx3i (a, b, nx, ny, nz)
+
+int	a[nx, ny, nz], b[nz, ny, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, y, x] = a[x, y, z]
+end
+
+
+
+# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x.  The arrays need not be
+# identical.
+
+procedure tzyx3l (a, b, nx, ny, nz)
+
+long	a[nx, ny, nz], b[nz, ny, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, y, x] = a[x, y, z]
+end
+
+
+
+# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x.  The arrays need not be
+# identical.
+
+procedure tzyx3r (a, b, nx, ny, nz)
+
+real	a[nx, ny, nz], b[nz, ny, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, y, x] = a[x, y, z]
+end
+
+
+
+# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x.  The arrays need not be
+# identical.
+
+procedure tzyx3d (a, b, nx, ny, nz)
+
+double	a[nx, ny, nz], b[nz, ny, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, y, x] = a[x, y, z]
+end
+
+
+
+# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x.  The arrays need not be
+# identical.
+
+procedure tzyx3x (a, b, nx, ny, nz)
+
+complex	a[nx, ny, nz], b[nz, ny, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, y, x] = a[x, y, z]
+end
+
+
diff --git a/pkg/proto/vol/src/im3dtran/x_im3dtran.x b/pkg/proto/vol/src/im3dtran/x_im3dtran.x
new file mode 100644
index 00000000..b1610b21
--- /dev/null
+++ b/pkg/proto/vol/src/im3dtran/x_im3dtran.x
@@ -0,0 +1,4 @@
+# X_IM3DTRANS.X -- Task statement for IM3DTRANSPOSE.  Used only for debugging
+# (see entry 'dbx:' in mkpkg).
+
+task	im3dtrans = t_im3dtrans
diff --git a/pkg/proto/vol/src/imjoin.gx b/pkg/proto/vol/src/imjoin.gx
new file mode 100644
index 00000000..04d2cc93
--- /dev/null
+++ b/pkg/proto/vol/src/imjoin.gx
@@ -0,0 +1,86 @@
+include <imhdr.h>
+
+define	VPTR		Memi[$1+$2-1]	# Array of axis vector pointers
+
+
+# IMJOIN -- Join the set of input images into an output image along the
+# specified axis, any dimension up to 7 (NOT necessarily IM_MAXDIM!).
+
+$for (silrdx)
+procedure imjoin$t (inptr, nimages, out, joindim, outtype)
+pointer	inptr[nimages]		# Input IMIO pointers
+int	nimages			# Number of input images
+pointer	out			# Output IMIO pointer
+int	joindim			# Dimension along which to join images
+int	outtype			# Output datatype
+
+pointer	in, inbuf, outbuf, sp, vin
+int	stat, image, cum_len, line, band, dim4, dim5, dim6, dim7
+long	vout[IM_MAXDIM]
+
+pointer	imgnl$t()
+pointer	impnl$t()
+
+begin
+	call smark (sp)
+	call salloc (vin, nimages, TY_INT)
+
+	call amovkl (long(1), vout, IM_MAXDIM)
+	do image = 1, nimages {
+	    call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG)
+	    call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM)
+	}
+
+	# Join input images along specified dimension.  Joins along columns
+	# and lines require processing in special order, all others in the
+	# same order.  In the first two cases we process all input images
+	# in inner loops, so we have to keep all those imdescriptors open.
+
+	switch (joindim) {
+	case 1:						# join columns
+	    do band = 1, IM_LEN(out,3)
+		do line = 1, IM_LEN(out,2) {
+		    stat = impnl$t (out, outbuf, vout)
+		    cum_len = 0
+		    do image = 1, nimages {
+			in = inptr[image]
+			stat = imgnl$t (in, inbuf, Meml[VPTR(vin,image)])
+			call amov$t (Mem$t[inbuf], Mem$t[outbuf+cum_len],
+			    IM_LEN(in,1))
+			cum_len = cum_len + IM_LEN(in,1)
+		    }
+		}
+	case 2:						# join lines
+	    do band = 1, IM_LEN(out,3)
+		do image = 1, nimages {
+		    in = inptr[image]
+		    do line = 1, IM_LEN(in,2) {
+			stat = impnl$t (out, outbuf, vout)
+			stat = imgnl$t (in, inbuf, Meml[VPTR(vin,image)])
+			call amov$t (Mem$t[inbuf], Mem$t[outbuf], IM_LEN(in,1))
+		    }
+		}
+	case 3,4,5,6,7:					# join bands or higher
+	    do image = 1, nimages {
+		in = inptr[image]
+		do dim7 = 1, IM_LEN(in,7)
+		    do dim6 = 1, IM_LEN(in,6)
+			do dim5 = 1, IM_LEN(in,5)
+			    do dim4 = 1, IM_LEN(in,4)
+				do band = 1, IM_LEN(in,3)
+				    do line = 1, IM_LEN(in,2) {
+					stat = impnl$t (out, outbuf, vout)
+					stat = imgnl$t (in, inbuf,
+					    Meml[VPTR(vin,image)])
+					call amov$t (Mem$t[inbuf],
+					    Mem$t[outbuf], IM_LEN(in,1))
+				    }
+		# Unmap last image to free resources.
+		call imunmap (in)
+	    }
+	}
+
+	call sfree (sp)
+end
+
+$endfor
diff --git a/pkg/proto/vol/src/imjoin.par b/pkg/proto/vol/src/imjoin.par
new file mode 100644
index 00000000..3acb6e02
--- /dev/null
+++ b/pkg/proto/vol/src/imjoin.par
@@ -0,0 +1,4 @@
+input,s,a,,,,"Input images or @file"
+output,s,a,,,,"Output joined image"
+joindim,i,h,1,1,7,"Splice dimension (1=x, 2=y, 3=z, ...)"
+outtype,s,h,"",,,"Output datatype (defaults to highest intype)"
diff --git a/pkg/proto/vol/src/imjoin.x b/pkg/proto/vol/src/imjoin.x
new file mode 100644
index 00000000..77ce47f3
--- /dev/null
+++ b/pkg/proto/vol/src/imjoin.x
@@ -0,0 +1,471 @@
+include <imhdr.h>
+
+define	VPTR		Memi[$1+$2-1]	# Array of axis vector pointers
+
+
+# IMJOIN -- Join the set of input images into an output image along the
+# specified axis, any dimension up to 7 (NOT necessarily IM_MAXDIM!).
+
+
+procedure imjoins (inptr, nimages, out, joindim, outtype)
+pointer	inptr[nimages]		# Input IMIO pointers
+int	nimages			# Number of input images
+pointer	out			# Output IMIO pointer
+int	joindim			# Dimension along which to join images
+int	outtype			# Output datatype
+
+pointer	in, inbuf, outbuf, sp, vin
+int	stat, image, cum_len, line, band, dim4, dim5, dim6, dim7
+long	vout[IM_MAXDIM]
+
+pointer	imgnls()
+pointer	impnls()
+
+begin
+	call smark (sp)
+	call salloc (vin, nimages, TY_INT)
+
+	call amovkl (long(1), vout, IM_MAXDIM)
+	do image = 1, nimages {
+	    call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG)
+	    call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM)
+	}
+
+	# Join input images along specified dimension.  Joins along columns
+	# and lines require processing in special order, all others in the
+	# same order.  In the first two cases we process all input images
+	# in inner loops, so we have to keep all those imdescriptors open.
+
+	switch (joindim) {
+	case 1:						# join columns
+	    do band = 1, IM_LEN(out,3)
+		do line = 1, IM_LEN(out,2) {
+		    stat = impnls (out, outbuf, vout)
+		    cum_len = 0
+		    do image = 1, nimages {
+			in = inptr[image]
+			stat = imgnls (in, inbuf, Meml[VPTR(vin,image)])
+			call amovs (Mems[inbuf], Mems[outbuf+cum_len],
+			    IM_LEN(in,1))
+			cum_len = cum_len + IM_LEN(in,1)
+		    }
+		}
+	case 2:						# join lines
+	    do band = 1, IM_LEN(out,3)
+		do image = 1, nimages {
+		    in = inptr[image]
+		    do line = 1, IM_LEN(in,2) {
+			stat = impnls (out, outbuf, vout)
+			stat = imgnls (in, inbuf, Meml[VPTR(vin,image)])
+			call amovs (Mems[inbuf], Mems[outbuf], IM_LEN(in,1))
+		    }
+		}
+	case 3,4,5,6,7:					# join bands or higher
+	    do image = 1, nimages {
+		in = inptr[image]
+		do dim7 = 1, IM_LEN(in,7)
+		    do dim6 = 1, IM_LEN(in,6)
+			do dim5 = 1, IM_LEN(in,5)
+			    do dim4 = 1, IM_LEN(in,4)
+				do band = 1, IM_LEN(in,3)
+				    do line = 1, IM_LEN(in,2) {
+					stat = impnls (out, outbuf, vout)
+					stat = imgnls (in, inbuf,
+					    Meml[VPTR(vin,image)])
+					call amovs (Mems[inbuf],
+					    Mems[outbuf], IM_LEN(in,1))
+				    }
+		# Unmap last image to free resources.
+		call imunmap (in)
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+procedure imjoini (inptr, nimages, out, joindim, outtype)
+pointer	inptr[nimages]		# Input IMIO pointers
+int	nimages			# Number of input images
+pointer	out			# Output IMIO pointer
+int	joindim			# Dimension along which to join images
+int	outtype			# Output datatype
+
+pointer	in, inbuf, outbuf, sp, vin
+int	stat, image, cum_len, line, band, dim4, dim5, dim6, dim7
+long	vout[IM_MAXDIM]
+
+pointer	imgnli()
+pointer	impnli()
+
+begin
+	call smark (sp)
+	call salloc (vin, nimages, TY_INT)
+
+	call amovkl (long(1), vout, IM_MAXDIM)
+	do image = 1, nimages {
+	    call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG)
+	    call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM)
+	}
+
+	# Join input images along specified dimension.  Joins along columns
+	# and lines require processing in special order, all others in the
+	# same order.  In the first two cases we process all input images
+	# in inner loops, so we have to keep all those imdescriptors open.
+
+	switch (joindim) {
+	case 1:						# join columns
+	    do band = 1, IM_LEN(out,3)
+		do line = 1, IM_LEN(out,2) {
+		    stat = impnli (out, outbuf, vout)
+		    cum_len = 0
+		    do image = 1, nimages {
+			in = inptr[image]
+			stat = imgnli (in, inbuf, Meml[VPTR(vin,image)])
+			call amovi (Memi[inbuf], Memi[outbuf+cum_len],
+			    IM_LEN(in,1))
+			cum_len = cum_len + IM_LEN(in,1)
+		    }
+		}
+	case 2:						# join lines
+	    do band = 1, IM_LEN(out,3)
+		do image = 1, nimages {
+		    in = inptr[image]
+		    do line = 1, IM_LEN(in,2) {
+			stat = impnli (out, outbuf, vout)
+			stat = imgnli (in, inbuf, Meml[VPTR(vin,image)])
+			call amovi (Memi[inbuf], Memi[outbuf], IM_LEN(in,1))
+		    }
+		}
+	case 3,4,5,6,7:					# join bands or higher
+	    do image = 1, nimages {
+		in = inptr[image]
+		do dim7 = 1, IM_LEN(in,7)
+		    do dim6 = 1, IM_LEN(in,6)
+			do dim5 = 1, IM_LEN(in,5)
+			    do dim4 = 1, IM_LEN(in,4)
+				do band = 1, IM_LEN(in,3)
+				    do line = 1, IM_LEN(in,2) {
+					stat = impnli (out, outbuf, vout)
+					stat = imgnli (in, inbuf,
+					    Meml[VPTR(vin,image)])
+					call amovi (Memi[inbuf],
+					    Memi[outbuf], IM_LEN(in,1))
+				    }
+		# Unmap last image to free resources.
+		call imunmap (in)
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+procedure imjoinl (inptr, nimages, out, joindim, outtype)
+pointer	inptr[nimages]		# Input IMIO pointers
+int	nimages			# Number of input images
+pointer	out			# Output IMIO pointer
+int	joindim			# Dimension along which to join images
+int	outtype			# Output datatype
+
+pointer	in, inbuf, outbuf, sp, vin
+int	stat, image, cum_len, line, band, dim4, dim5, dim6, dim7
+long	vout[IM_MAXDIM]
+
+pointer	imgnll()
+pointer	impnll()
+
+begin
+	call smark (sp)
+	call salloc (vin, nimages, TY_INT)
+
+	call amovkl (long(1), vout, IM_MAXDIM)
+	do image = 1, nimages {
+	    call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG)
+	    call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM)
+	}
+
+	# Join input images along specified dimension.  Joins along columns
+	# and lines require processing in special order, all others in the
+	# same order.  In the first two cases we process all input images
+	# in inner loops, so we have to keep all those imdescriptors open.
+
+	switch (joindim) {
+	case 1:						# join columns
+	    do band = 1, IM_LEN(out,3)
+		do line = 1, IM_LEN(out,2) {
+		    stat = impnll (out, outbuf, vout)
+		    cum_len = 0
+		    do image = 1, nimages {
+			in = inptr[image]
+			stat = imgnll (in, inbuf, Meml[VPTR(vin,image)])
+			call amovl (Meml[inbuf], Meml[outbuf+cum_len],
+			    IM_LEN(in,1))
+			cum_len = cum_len + IM_LEN(in,1)
+		    }
+		}
+	case 2:						# join lines
+	    do band = 1, IM_LEN(out,3)
+		do image = 1, nimages {
+		    in = inptr[image]
+		    do line = 1, IM_LEN(in,2) {
+			stat = impnll (out, outbuf, vout)
+			stat = imgnll (in, inbuf, Meml[VPTR(vin,image)])
+			call amovl (Meml[inbuf], Meml[outbuf], IM_LEN(in,1))
+		    }
+		}
+	case 3,4,5,6,7:					# join bands or higher
+	    do image = 1, nimages {
+		in = inptr[image]
+		do dim7 = 1, IM_LEN(in,7)
+		    do dim6 = 1, IM_LEN(in,6)
+			do dim5 = 1, IM_LEN(in,5)
+			    do dim4 = 1, IM_LEN(in,4)
+				do band = 1, IM_LEN(in,3)
+				    do line = 1, IM_LEN(in,2) {
+					stat = impnll (out, outbuf, vout)
+					stat = imgnll (in, inbuf,
+					    Meml[VPTR(vin,image)])
+					call amovl (Meml[inbuf],
+					    Meml[outbuf], IM_LEN(in,1))
+				    }
+		# Unmap last image to free resources.
+		call imunmap (in)
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+procedure imjoinr (inptr, nimages, out, joindim, outtype)
+pointer	inptr[nimages]		# Input IMIO pointers
+int	nimages			# Number of input images
+pointer	out			# Output IMIO pointer
+int	joindim			# Dimension along which to join images
+int	outtype			# Output datatype
+
+pointer	in, inbuf, outbuf, sp, vin
+int	stat, image, cum_len, line, band, dim4, dim5, dim6, dim7
+long	vout[IM_MAXDIM]
+
+pointer	imgnlr()
+pointer	impnlr()
+
+begin
+	call smark (sp)
+	call salloc (vin, nimages, TY_INT)
+
+	call amovkl (long(1), vout, IM_MAXDIM)
+	do image = 1, nimages {
+	    call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG)
+	    call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM)
+	}
+
+	# Join input images along specified dimension.  Joins along columns
+	# and lines require processing in special order, all others in the
+	# same order.  In the first two cases we process all input images
+	# in inner loops, so we have to keep all those imdescriptors open.
+
+	switch (joindim) {
+	case 1:						# join columns
+	    do band = 1, IM_LEN(out,3)
+		do line = 1, IM_LEN(out,2) {
+		    stat = impnlr (out, outbuf, vout)
+		    cum_len = 0
+		    do image = 1, nimages {
+			in = inptr[image]
+			stat = imgnlr (in, inbuf, Meml[VPTR(vin,image)])
+			call amovr (Memr[inbuf], Memr[outbuf+cum_len],
+			    IM_LEN(in,1))
+			cum_len = cum_len + IM_LEN(in,1)
+		    }
+		}
+	case 2:						# join lines
+	    do band = 1, IM_LEN(out,3)
+		do image = 1, nimages {
+		    in = inptr[image]
+		    do line = 1, IM_LEN(in,2) {
+			stat = impnlr (out, outbuf, vout)
+			stat = imgnlr (in, inbuf, Meml[VPTR(vin,image)])
+			call amovr (Memr[inbuf], Memr[outbuf], IM_LEN(in,1))
+		    }
+		}
+	case 3,4,5,6,7:					# join bands or higher
+	    do image = 1, nimages {
+		in = inptr[image]
+		do dim7 = 1, IM_LEN(in,7)
+		    do dim6 = 1, IM_LEN(in,6)
+			do dim5 = 1, IM_LEN(in,5)
+			    do dim4 = 1, IM_LEN(in,4)
+				do band = 1, IM_LEN(in,3)
+				    do line = 1, IM_LEN(in,2) {
+					stat = impnlr (out, outbuf, vout)
+					stat = imgnlr (in, inbuf,
+					    Meml[VPTR(vin,image)])
+					call amovr (Memr[inbuf],
+					    Memr[outbuf], IM_LEN(in,1))
+				    }
+		# Unmap last image to free resources.
+		call imunmap (in)
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+procedure imjoind (inptr, nimages, out, joindim, outtype)
+pointer	inptr[nimages]		# Input IMIO pointers
+int	nimages			# Number of input images
+pointer	out			# Output IMIO pointer
+int	joindim			# Dimension along which to join images
+int	outtype			# Output datatype
+
+pointer	in, inbuf, outbuf, sp, vin
+int	stat, image, cum_len, line, band, dim4, dim5, dim6, dim7
+long	vout[IM_MAXDIM]
+
+pointer	imgnld()
+pointer	impnld()
+
+begin
+	call smark (sp)
+	call salloc (vin, nimages, TY_INT)
+
+	call amovkl (long(1), vout, IM_MAXDIM)
+	do image = 1, nimages {
+	    call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG)
+	    call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM)
+	}
+
+	# Join input images along specified dimension.  Joins along columns
+	# and lines require processing in special order, all others in the
+	# same order.  In the first two cases we process all input images
+	# in inner loops, so we have to keep all those imdescriptors open.
+
+	switch (joindim) {
+	case 1:						# join columns
+	    do band = 1, IM_LEN(out,3)
+		do line = 1, IM_LEN(out,2) {
+		    stat = impnld (out, outbuf, vout)
+		    cum_len = 0
+		    do image = 1, nimages {
+			in = inptr[image]
+			stat = imgnld (in, inbuf, Meml[VPTR(vin,image)])
+			call amovd (Memd[inbuf], Memd[outbuf+cum_len],
+			    IM_LEN(in,1))
+			cum_len = cum_len + IM_LEN(in,1)
+		    }
+		}
+	case 2:						# join lines
+	    do band = 1, IM_LEN(out,3)
+		do image = 1, nimages {
+		    in = inptr[image]
+		    do line = 1, IM_LEN(in,2) {
+			stat = impnld (out, outbuf, vout)
+			stat = imgnld (in, inbuf, Meml[VPTR(vin,image)])
+			call amovd (Memd[inbuf], Memd[outbuf], IM_LEN(in,1))
+		    }
+		}
+	case 3,4,5,6,7:					# join bands or higher
+	    do image = 1, nimages {
+		in = inptr[image]
+		do dim7 = 1, IM_LEN(in,7)
+		    do dim6 = 1, IM_LEN(in,6)
+			do dim5 = 1, IM_LEN(in,5)
+			    do dim4 = 1, IM_LEN(in,4)
+				do band = 1, IM_LEN(in,3)
+				    do line = 1, IM_LEN(in,2) {
+					stat = impnld (out, outbuf, vout)
+					stat = imgnld (in, inbuf,
+					    Meml[VPTR(vin,image)])
+					call amovd (Memd[inbuf],
+					    Memd[outbuf], IM_LEN(in,1))
+				    }
+		# Unmap last image to free resources.
+		call imunmap (in)
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+procedure imjoinx (inptr, nimages, out, joindim, outtype)
+pointer	inptr[nimages]		# Input IMIO pointers
+int	nimages			# Number of input images
+pointer	out			# Output IMIO pointer
+int	joindim			# Dimension along which to join images
+int	outtype			# Output datatype
+
+pointer	in, inbuf, outbuf, sp, vin
+int	stat, image, cum_len, line, band, dim4, dim5, dim6, dim7
+long	vout[IM_MAXDIM]
+
+pointer	imgnlx()
+pointer	impnlx()
+
+begin
+	call smark (sp)
+	call salloc (vin, nimages, TY_INT)
+
+	call amovkl (long(1), vout, IM_MAXDIM)
+	do image = 1, nimages {
+	    call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG)
+	    call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM)
+	}
+
+	# Join input images along specified dimension.  Joins along columns
+	# and lines require processing in special order, all others in the
+	# same order.  In the first two cases we process all input images
+	# in inner loops, so we have to keep all those imdescriptors open.
+
+	switch (joindim) {
+	case 1:						# join columns
+	    do band = 1, IM_LEN(out,3)
+		do line = 1, IM_LEN(out,2) {
+		    stat = impnlx (out, outbuf, vout)
+		    cum_len = 0
+		    do image = 1, nimages {
+			in = inptr[image]
+			stat = imgnlx (in, inbuf, Meml[VPTR(vin,image)])
+			call amovx (Memx[inbuf], Memx[outbuf+cum_len],
+			    IM_LEN(in,1))
+			cum_len = cum_len + IM_LEN(in,1)
+		    }
+		}
+	case 2:						# join lines
+	    do band = 1, IM_LEN(out,3)
+		do image = 1, nimages {
+		    in = inptr[image]
+		    do line = 1, IM_LEN(in,2) {
+			stat = impnlx (out, outbuf, vout)
+			stat = imgnlx (in, inbuf, Meml[VPTR(vin,image)])
+			call amovx (Memx[inbuf], Memx[outbuf], IM_LEN(in,1))
+		    }
+		}
+	case 3,4,5,6,7:					# join bands or higher
+	    do image = 1, nimages {
+		in = inptr[image]
+		do dim7 = 1, IM_LEN(in,7)
+		    do dim6 = 1, IM_LEN(in,6)
+			do dim5 = 1, IM_LEN(in,5)
+			    do dim4 = 1, IM_LEN(in,4)
+				do band = 1, IM_LEN(in,3)
+				    do line = 1, IM_LEN(in,2) {
+					stat = impnlx (out, outbuf, vout)
+					stat = imgnlx (in, inbuf,
+					    Meml[VPTR(vin,image)])
+					call amovx (Memx[inbuf],
+					    Memx[outbuf], IM_LEN(in,1))
+				    }
+		# Unmap last image to free resources.
+		call imunmap (in)
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
diff --git a/pkg/proto/vol/src/imminmax.x b/pkg/proto/vol/src/imminmax.x
new file mode 100644
index 00000000..dea537c6
--- /dev/null
+++ b/pkg/proto/vol/src/imminmax.x
@@ -0,0 +1,73 @@
+include	<imhdr.h>
+
+# IM_MINMAX -- Compute the minimum and maximum pixel values of an image.
+# Works for images of any dimensionality, size, or datatype, although
+# the min and max values can currently only be stored in the image header
+# as real values.
+
+procedure im_minmax (im, min_value, max_value)
+
+pointer	im				# image descriptor
+real	min_value			# minimum pixel value in image (out)
+real	max_value			# maximum pixel value in image (out)
+
+pointer	buf
+bool	first_line
+long	v[IM_MAXDIM]
+short	minval_s, maxval_s
+long	minval_l, maxval_l
+real	minval_r, maxval_r
+int	imgnls(), imgnll(), imgnlr()
+
+begin
+	call amovkl (long(1), v, IM_MAXDIM)		# start vector
+	first_line = true
+	min_value = INDEF
+	max_value = INDEF
+
+	switch (IM_PIXTYPE(im)) {
+	case TY_SHORT:
+	    while (imgnls (im, buf, v) != EOF) {
+		call alims (Mems[buf], IM_LEN(im,1), minval_s, maxval_s)
+		if (first_line) {
+		    min_value = minval_s
+		    max_value = maxval_s
+		    first_line = false
+		} else {
+		    if (minval_s < min_value)
+			min_value = minval_s
+		    if (maxval_s > max_value)
+			max_value = maxval_s
+		}
+	    }
+	case TY_USHORT, TY_INT, TY_LONG:
+	    while (imgnll (im, buf, v) != EOF) {
+		call aliml (Meml[buf], IM_LEN(im,1), minval_l, maxval_l)
+		if (first_line) {
+		    min_value = minval_l
+		    max_value = maxval_l
+		    first_line = false
+		} else {
+		    if (minval_l < min_value)
+			min_value = minval_l
+		    if (maxval_l > max_value)
+			max_value = maxval_l
+		}
+	    }
+	default:
+	    while (imgnlr (im, buf, v) != EOF) {
+		call alimr (Memr[buf], IM_LEN(im,1), minval_r, maxval_r)
+		if (first_line) {
+		    min_value = minval_r
+		    max_value = maxval_r
+		    first_line = false
+		} else {
+		    if (minval_r < min_value)
+			min_value = minval_r
+		    if (maxval_r > max_value)
+			max_value = maxval_r
+		}
+	    }
+	}
+end
+
diff --git a/pkg/proto/vol/src/mkpkg b/pkg/proto/vol/src/mkpkg
new file mode 100644
index 00000000..c0930db7
--- /dev/null
+++ b/pkg/proto/vol/src/mkpkg
@@ -0,0 +1,44 @@
+# Make the VOLUMES tasks.
+
+$call	relink
+$exit
+
+update:
+	$call relink
+	$call install
+	;
+
+relink:
+	$set	LIBS = "-lxtools"
+	$update	libpkg.a
+	$omake	x_vol.x
+	$link	x_vol.o libpkg.a $(LIBS) -o xx_vol.e
+	;
+
+install:
+	$move	xx_vol.e bin$x_vol.e
+	;
+
+tfiles:
+	$ifolder (vtransmit.x, vtransmit.gx)
+	    $generic -k vtransmit.gx -o vtransmit.x		$endif
+	$ifolder (imjoin.x, imjoin.gx)
+	    $generic -k imjoin.gx -o imjoin.x			$endif
+	;
+
+libpkg.a:
+	$ifeq (USE_GENERIC, yes) $call tfiles $endif
+
+	t_pvol.x	<ctype.h> <imhdr.h> pvol.h
+	vproject.x	<math.h> <imhdr.h> pvol.h
+	vmatrix.x	<imhdr.h> pvol.h
+	vtransmit.x	<imhdr.h> pvol.h
+	vgetincr.x	pvol.h
+	pv_gmem.x	<imhdr.h>
+	t_imjoin.x	<imhdr.h> <error.h> <syserr.h>
+	imjoin.x	<imhdr.h>
+	imminmax.x	<imhdr.h>
+	@im3dtran
+	@i2sun
+	;
+
diff --git a/pkg/proto/vol/src/pv_gmem.x b/pkg/proto/vol/src/pv_gmem.x
new file mode 100644
index 00000000..646baf47
--- /dev/null
+++ b/pkg/proto/vol/src/pv_gmem.x
@@ -0,0 +1,109 @@
+include <imhdr.h>
+
+# PV_GMEM -- Determine how much memory we can get and actually use for the
+# volume rotation sequence.  We only allocate actual memory temporarily in
+# order to see how much is really available; IMIO will later take care of
+# the actual io buffer allocation.
+
+define	DECR_MEM		0.8	# decrement mem by magic factor
+
+
+procedure pv_gmem (im1, im2, use_both, verbose, max_ws, len_x, oldsize)
+pointer	im1			# Input 3d image.
+pointer	im2			# Output projected image(s)
+bool	use_both		# Use both opacity and intensity voxels
+int	verbose			# Report memory usage?
+int	max_ws			# Maximum working set to allocate
+int	len_x			# (output) safe amount of memory to use
+int	oldsize			# (output) old memory to be reset at termination
+
+int	intype, outtype, reqmem, gotmem, needmem, maxsize
+int	yzslice_pix, yzreq
+pointer	buf_in
+bool	topped_out
+
+int	begmem(), sizeof()
+errchk	begmem(), malloc()
+
+begin
+	# See how much memory we can get; if we cannot get whole input image
+	# into memory, do it in chunks of yz slices.
+	intype = IM_PIXTYPE(im1)
+	outtype = IM_PIXTYPE(im2)
+	reqmem = IM_LEN(im1,1) * IM_LEN(im1,2) * IM_LEN(im1,3)
+	reqmem = reqmem * sizeof (intype)
+	if (use_both)
+	    reqmem = 2 * reqmem
+
+	# Add output buffer.
+	reqmem = reqmem + IM_LEN(im2,1) * sizeof (outtype)
+
+	# Decrease to max_ws (a task parameter in CHAR units).
+	reqmem = min (reqmem, max_ws)
+
+	repeat {
+	    iferr (gotmem = begmem (reqmem, oldsize, maxsize)) {
+		reqmem = reqmem * DECR_MEM
+		if (verbose == YES) {
+		    call eprintf ("ERR gotmem=begmem(); retrying at size %d\n")
+			call pargi (reqmem)
+		}
+	    } else {
+		if (verbose == YES) {
+		    call eprintf ("gotmem=%d, oldsize=%d, maxsize=%d\n")
+		    call pargi (gotmem)
+		    call pargi (oldsize)
+		    call pargi (maxsize)
+		}
+		break
+	    }
+	}
+	
+	# Make sure it is really available, and if not, decrement to largest
+	# number of yz slices possible.
+	needmem = gotmem
+	yzslice_pix = IM_LEN(im1,2) * IM_LEN(im1,3)
+	yzreq = yzslice_pix * sizeof(intype)
+	if (yzreq - IM_LEN(im1,1) * sizeof(TY_REAL) > needmem) {
+	    call eprintf ("Not enough memory for 1 yz slice of input image\n")
+	    call error (0, "Out of memory")
+	}
+	topped_out = false
+	repeat {
+	    iferr (call malloc (buf_in, needmem, intype)) {
+		needmem = needmem - yzreq
+		if (needmem < yzreq) {
+		    call eprintf ("Had to decrease memory too much")
+		    call error (0, "Memory allocation error (yzslice_pix)")
+		}
+		topped_out = true
+	    } else {
+		call mfree (buf_in, intype)
+		break
+	    }
+	}
+
+	# Experiments show that horrible things happen if we actually use
+	# this much memory.  Decrease by magic factor.
+	if (topped_out) {
+	    call fixmem (max (needmem, oldsize))
+	    if (verbose == YES) {
+		call eprintf ("Had to decrease memory for malloc().")
+		call eprintf ("  Working set now %d\n")
+		    call pargi (needmem)
+	    }
+	    needmem = needmem * DECR_MEM
+	    if (verbose == YES) {
+		call eprintf ("Remaining memory for image buffers = %d\n")
+		    call pargi (needmem)
+	    }
+	}
+	if (needmem < yzreq) {
+	    call eprintf ("Not enough memory for 1 yz slice of input image\n")
+	    call error (0, "Out of memory")
+	}
+
+	# We return the number of columns to gulp from the input image at one
+	# time and oldmem so the task can release its memory on completion.
+	len_x = needmem / yzreq
+end
diff --git a/pkg/proto/vol/src/pvol.h b/pkg/proto/vol/src/pvol.h
new file mode 100644
index 00000000..e232e8b1
--- /dev/null
+++ b/pkg/proto/vol/src/pvol.h
@@ -0,0 +1,58 @@
+# PVOL.H -- PVOL definitions.
+
+define	COL		1			# image column index
+define	LINE		2			# image line index
+define	BAND		3			# image band index
+define	DIS		(sqrt (($1)*($1) + ($2)*($2)))
+define	DRADIAN		(57.295779513082320877D0)
+define	DDEGTORAD	(double($1)/DRADIAN)
+define	DRADTODEG	(double($1)*DRADIAN)
+define	DPI		3.1415926535897932385D0
+define	DTWOPI		6.2831853071795864769D0
+define	DHALFPI		1.5707963267948966192D0
+define	DTHREEPIOVER2	(1.5D0 * DPI)
+define	DEF_IMIN	(0.0)
+define	DEF_IMAX	(1.0)
+
+define	ALG_INCREM	1			# incremental dda proj. algor.
+define	ALG_BRESEN	2			# bresenham dda proj. algor.
+define	ALG_MATRIX	3			# rotation matrix prol. algor.
+define	P_ATTENUATE	1			# attenuate by voxval (opacity)
+define	P_AVERAGE	2			# average projected intensities
+define	P_SUM		3			# sum voxel intensities
+define	P_INVDISPOW	4			# wt int. by inverse dis power
+define	P_MODN		5			# use only f(ndecades) voxels
+define	P_LASTONLY	6			# use only last voxval > cutoff
+
+# Volume rotation descriptor
+define	LEN_VP		30
+
+# Projection geometry elements:
+define	P_ALGORITHM	Memi[$1]		# Projection algorithm
+define	DEGREES		Memr[P2R($1+1)]		# Degrees per rotation increment
+define	NFRAMES		Memi[$1+2]		# Number of rotation increments
+define	VECX		Memr[P2R($1+3)]		# Rotation axis X vector
+define	VECY		Memr[P2R($1+4)]		# Rotation axis Y vector
+define	VECZ		Memr[P2R($1+5)]		# Rotation axis Z vector
+define	INIT_THETA	Memr[P2R($1+6)]		# Initial rotation angle
+define	MAX_WS		Memi[$1+7]		# Maximum working set size
+#			reserved space
+
+# Light transmission elements:
+define	OPACELEM	Memi[$1+10]		# Opacity element in 4th dimen
+define	PTYPE		Memi[$1+11]		# Projection type, voxel val 1
+define	OMIN		Memr[P2R($1+12)]	# Voxel opacity minimum
+define	OMAX		Memr[P2R($1+13)]	# Voxel opacity maximum
+define	OSCALE		Memr[P2R($1+14)]	# Opacity scale factor
+define	AMIN		Memr[P2R($1+15)]	# Attenuation factor minimum
+define	AMAX		Memr[P2R($1+16)]	# Attenuation factor maximum
+define	VIMIN		Memr[P2R($1+17)]	# Voxel intensity minimum
+define	VIMAX		Memr[P2R($1+18)]	# Voxel intensity maximum
+define	IZERO		Memr[P2R($1+19)]	# Background illumination
+define	DISPOWER	Memr[P2R($1+20)]	# Distance weighting power
+define	MODN		Memi[$1+21]		# Use vox w/ (mod(val*100,modn))
+define	IIMIN		Memr[P2R($1+22)]	# Input intensity minimum
+define	IIMAX		Memr[P2R($1+23)]	# Input intensity maximum
+define	VERBOSE		Memi[$1+24]		# Write verbose output?
+define	DISCUTOFF	Memi[$1+25]		# Measure distance w/in cutoffs
+#			reserved space
diff --git a/pkg/proto/vol/src/pvol.par b/pkg/proto/vol/src/pvol.par
new file mode 100644
index 00000000..3fbc38ba
--- /dev/null
+++ b/pkg/proto/vol/src/pvol.par
@@ -0,0 +1,25 @@
+input,s,a,,,,"Input 3d or 4d image"
+output,s,a,,,,"Output datacube"
+nframes,i,h,INDEF,1,,"Number of rotation frames to compute"
+degrees,r,h,10.0,,,"Degrees per rotation increment"
+theta0,r,h,0.0,,,"Initial rotation angle (ccw from +X)"
+ptype,i,h,2,1,6,"Projection (1=opc 2=av 3=sum 4=invd 5=mod 6=lst)"
+imin,r,h,INDEF,,,"Voxel intensity minimum cutoff"
+imax,r,h,INDEF,,,"Voxel intensity maximum cutoff"
+omin,r,h,INDEF,,,"Voxel opacity minimum cutoff (ptype=1)"
+omax,r,h,INDEF,,,"Voxel opacity maximum cutoff (ptype=1)"
+amin,r,h,0.0,0.0,1.0,"Minimum attenuation factor (ptype=1)"
+amax,r,h,1.0,0.0,1.0,"Maximum attenuation factor (ptype=1)"
+izero,r,h,1.0,,,"Initial intensity (background illumination, ptype=1)"
+oscale,r,h,1.0,,,"Voxel opacity scale factor (ptype=1)"
+opacelem,i,h,1,1,2,"4th dim. opacity element (other=intensity)"
+dispower,r,h,2.0,,,"Inverse distance weighting power (ptype=4,5)"
+discutoff,b,h,n,,,"Measure distance from 1st vox inside cutoff"
+modn,i,h,10,1,100,"Mod(n) for ptype=6; used for high-contrast input"
+vecx,r,h,1.0,-1.0,1.0,"Rotation axis X vector (right hand rule)"
+vecy,r,h,0.0,-1.0,1.0,"Rotation axis Y vector"
+vecz,r,h,0.0,-1.0,1.0,"Rotation axis Z vector"
+title,s,h,"",,,"Title for rotation sequence"
+maxws,i,h,2000000,256000,,"Max workingset size in CHARS (2 bytes usually)"
+abs,b,h,no,,,"Take absolute value of pixels?"
+verbose,b,h,yes,,,"Verbose?  (report progress, memory usage)"
diff --git a/pkg/proto/vol/src/t_imjoin.x b/pkg/proto/vol/src/t_imjoin.x
new file mode 100644
index 00000000..70708571
--- /dev/null
+++ b/pkg/proto/vol/src/t_imjoin.x
@@ -0,0 +1,190 @@
+include	<imhdr.h>
+include	<error.h>
+include	<syserr.h>
+
+define	DEFBUFSIZE		65536		# default IMIO buffer size
+define	FUDGE			0.8		# fudge factor
+
+
+# T_IMJOIN -- Task driver for imjoin:  up to IM_MAXDIM image join, along
+# any one specified axis, from multiple input images.  The set of input
+# images need have the same number of dimensions and elements per dimension
+# ONLY along the axes not being joined.  Datatype will be converted to 
+# highest precedence type if not all the same.
+
+procedure t_imjoin()
+
+int	list			# List of input images
+pointer	output			# Output image
+char	outtype			# Output datatype
+
+int	i, j, nimages, intype, ndim, joindim, outdtype, nelems[IM_MAXDIM]
+int	bufsize, maxsize, memory, oldsize
+pointer	sp, in, out, im, im1, input
+
+int	imtopenp(), imtlen(), imtgetim(), clgeti()
+int	ty_max(), sizeof(), begmem(), errcode()
+char	clgetc()
+pointer	immap()
+errchk	immap
+define	retry_	99
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+
+	# Get the parameters.  Some parameters are obtained later.
+	list = imtopenp ("input")
+	call clgstr ("output", Memc[output], SZ_FNAME)
+	joindim = clgeti ("joindim")
+	outtype = clgetc ("outtype")
+
+	# Check if there are no images.
+	nimages = imtlen (list)
+	if (nimages == 0) {
+	    call imtclose (list)
+	    call sfree (sp)
+	    call error (0, "No input images to join")
+	}
+	call salloc (in, nimages, TY_POINTER)
+
+	# Map the input images.
+	bufsize = 0
+retry_
+	nimages = 0
+	while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+	    nimages = nimages + 1
+	    Memi[in+nimages-1] = immap (Memc[input], READ_ONLY, 0)
+	}
+
+	# Determine the dimensionality, size, and datatype of the output image.
+	im = Memi[in]
+	intype = IM_PIXTYPE(im)
+	ndim = max (IM_NDIM(im), joindim)
+	do j = 1, ndim
+	    nelems[j] = IM_LEN(im,j)
+
+	do i = 2, nimages {
+	    im1 = Memi[in+i-1]
+	    ndim = max (ndim, IM_NDIM(im1))
+	    do j = 1, ndim {
+		if (j == joindim)
+		    nelems[j] = nelems[j] + IM_LEN(im1,j)
+		else {
+		    if (IM_LEN(im1,j) != nelems[j]) {
+			call eprintf ("Image %d different size in dimen %d\n")
+			    call pargi (i)
+			    call pargi (IM_LEN(im1,j))
+			call error (1, "Non-joindim image sizes must match")
+		    }
+		}
+	    }
+	    intype = ty_max (intype, IM_PIXTYPE(im1))
+	}
+
+	# Open the output image and set its pixel datatype.
+	# If outtype was not specified (the default), set it to intype.
+
+	out = immap (Memc[output], NEW_COPY, Memi[in])
+	switch (outtype) {
+	case 's':
+	    outdtype = TY_SHORT
+	case 'i':
+	    outdtype = TY_INT
+	case 'l':
+	    outdtype = TY_LONG
+	case 'r':
+	    outdtype = TY_REAL
+	case 'd':
+	    outdtype = TY_DOUBLE
+	case 'x':
+	    outdtype = TY_COMPLEX
+	default:
+	    outdtype = intype
+	}
+	IM_PIXTYPE(out) = outdtype
+
+	# Set output image dimensionality and axis lengths.
+	IM_NDIM(out) = ndim
+	do j = 1, ndim
+	    IM_LEN(out,j) = nelems[j]
+
+	if (bufsize == 0) {
+	    # Set initial IMIO buffer size based on the number of images
+	    # and maximum amount of working memory available.  The buffer
+	    # size may be adjusted later if the task runs out of memory.
+	    # The FUDGE factor is used to allow for the size of the
+	    # program, memory allocator inefficiencies, and any other
+	    # memory requirements besides IMIO.
+
+	    bufsize = 1
+	    do i = 1, IM_NDIM(out)
+		bufsize = bufsize * IM_LEN(out,i)
+	    bufsize = bufsize * sizeof (intype)
+	    bufsize = min (bufsize, DEFBUFSIZE)
+	    memory = begmem ((nimages + 1) * bufsize, oldsize, maxsize)
+	    memory = min (memory, int (FUDGE * maxsize))
+	    bufsize = memory / (nimages + 1)
+	}
+
+	# Join the images along joindim.  If an out of memory error occurs
+	# close all images and files, divide the IMIO buffer size in half
+	# and try again.
+	iferr {
+	    switch (intype) {
+	    case TY_SHORT:
+		call imjoins (Memi[in], nimages, out, joindim, outdtype)
+	    case TY_INT:
+		call imjoini (Memi[in], nimages, out, joindim, outdtype)
+	    case TY_LONG:
+		call imjoinl (Memi[in], nimages, out, joindim, outdtype)
+	    case TY_REAL:
+		call imjoinr (Memi[in], nimages, out, joindim, outdtype)
+	    case TY_DOUBLE:
+		call imjoind (Memi[in], nimages, out, joindim, outdtype)
+	    case TY_COMPLEX:
+		call imjoinx (Memi[in], nimages, out, joindim, outdtype)
+	    }
+	} then {
+	    switch (errcode()) {
+	    case SYS_MFULL:
+		do j = 1, nimages
+		    call imunmap (Memi[in+j-1])
+		call imunmap (out)
+		call imdelete (Memc[output])
+		call imtrew (list)
+		bufsize = bufsize / 2
+		goto retry_
+	    default:
+		call erract (EA_ERROR)
+	    }
+	}
+
+	# Unmap all the images and restore memory.
+	call imunmap (out)
+	do i = 1, nimages
+	    if (joindim < 3)
+		call imunmap (Memi[in+i-1])
+
+	call sfree (sp)
+	call fixmem (oldsize)
+end
+
+
+# TY_MAX -- Return the datatype of highest precedence.
+
+int procedure ty_max (type1, type2)
+
+int	type1, type2		# Datatypes
+
+int	i, j, order[7]
+data	order/TY_SHORT,TY_INT,TY_LONG,TY_REAL,TY_DOUBLE,TY_COMPLEX,TY_REAL/
+
+begin
+	for (i=1; (i<=6) && (type1!=order[i]); i=i+1)
+	    ;
+	for (j=1; (j<=6) && (type2!=order[j]); j=j+1)
+	    ;
+	return (order[max(i,j)])
+end
diff --git a/pkg/proto/vol/src/t_pvol.x b/pkg/proto/vol/src/t_pvol.x
new file mode 100644
index 00000000..251012c5
--- /dev/null
+++ b/pkg/proto/vol/src/t_pvol.x
@@ -0,0 +1,284 @@
+include <ctype.h>
+include <imhdr.h>
+include <time.h>
+include "pvol.h"
+
+define	iwrapup_	91
+define	mwrapup_	92
+
+
+# PVOL -- Project Volume.  Given an input datacube, produce a series of
+# frames representing projections at stepped rotations around the cube,
+# using voxel intensity and/or opacity information.  This is a form of
+# volume rendering.
+
+procedure t_pvol
+
+pointer	input, output, sp, tmpstr, vp, timestr, im1, im2
+long	clock1, clock2, elapclock, cpu1, cpu2, elapcpu
+bool	need_lims, use_both
+real	tmpmin, tmpmax
+
+pointer immap()
+int	clgeti(), clktime(), cputime()
+bool	clgetb()
+real	clgetr()
+
+begin
+	call smark (sp)
+	call salloc (tmpstr, SZ_LINE, TY_CHAR)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (timestr, SZ_FNAME, TY_CHAR)
+
+	# Allocate storage for volume projection descriptor.
+	call malloc (vp, LEN_VP, TY_STRUCT)
+
+	# Input parameters.
+	if (clgetb ("verbose"))
+	    VERBOSE(vp) = YES
+	else
+	    VERBOSE(vp) = NO
+	call clgstr ("input", Memc[input], SZ_FNAME)
+	call clgstr ("output", Memc[output], SZ_FNAME)
+
+	# Geometric projection parameters:
+	DEGREES(vp) = clgetr ("degrees")
+	INIT_THETA(vp) = clgetr ("theta0")
+	NFRAMES(vp) = clgeti ("nframes")
+	if (IS_INDEFI(NFRAMES(vp)) && !IS_INDEFR(DEGREES(vp)))
+	    NFRAMES(vp) = int (360.0 / DEGREES(vp))
+	else if (IS_INDEFR(DEGREES(vp)) && !IS_INDEFI(NFRAMES(vp)))
+	    DEGREES(vp) = 360.0 / NFRAMES(vp)
+	else if (IS_INDEFR(DEGREES(vp)) && IS_INDEFI(NFRAMES(vp))) {
+	    NFRAMES(vp) = 36
+	    DEGREES(vp) = 10.0
+	}
+	PTYPE(vp) = clgeti ("ptype")
+	VIMIN(vp) = clgetr ("imin")
+	VIMAX(vp) = clgetr ("imax")
+	IZERO(vp) = clgetr ("izero")
+	OSCALE(vp) = clgetr ("oscale")
+	OMIN(vp) = clgetr ("omin")
+	OMAX(vp) = clgetr ("omax")
+	AMIN(vp) = clgetr ("amin")
+	AMAX(vp) = clgetr ("amax")
+	DISCUTOFF(vp) = NO
+	if (PTYPE(vp) == P_INVDISPOW || PTYPE(vp) == P_MODN) {
+	    DISPOWER(vp) = clgetr ("dispower")
+	    if (clgetb ("discutoff"))
+		DISCUTOFF(vp) = YES
+	}
+	if (PTYPE(vp) == P_MODN)
+	    MODN(vp) = clgeti ("modn")
+	VECX(vp) = clgetr ("vecx")
+	VECY(vp) = clgetr ("vecy")
+	VECZ(vp) = clgetr ("vecz")
+
+	MAX_WS(vp) = clgeti ("maxws")
+
+	# In prototype, the incremental algorithm is only implemented for
+	# rotations about the X axis, counterclockwise when viewed from +X
+	# looking back toward the origin.
+
+	if (!(VECX(vp) == +1.0 && VECY(vp) == 0.0 && VECZ(vp) == 0.0)) {
+	    call eprintf ("ERROR:  Only +X axis rotations supported with")
+	    call eprintf (" incremental alg. at present\n")
+	    call error (0, "Unsupported feature")
+	}
+
+	# Open images.
+	im1 = immap (Memc[input], READ_ONLY, 0)
+	im2 = immap (Memc[output], NEW_IMAGE, 0)
+	call clgstr ("title", IM_TITLE(im1), SZ_IMTITLE)
+
+	# If input image is 4d, with 2 elements in 4th dimension, one of them
+	# must be opacity and the other intensity.  If someone wants to merge
+	# two or more sets of intensity data, they can make independent runs
+	# of PVOL and merge the outputs using RGB displays.
+
+	use_both = false
+	OPACELEM(vp) = INDEFI
+	if (IM_NDIM(im1) == 4 && PTYPE(vp) != P_ATTENUATE) {
+	    if (IM_LEN(im1,4) > 2)
+		call error (0, "Don't know how to handle 4d image w/ >2 elems")
+	    else if (IM_LEN(im1,4) == 2) {
+		OPACELEM(vp) = clgeti ("opacelem")
+		if (PTYPE(vp) == P_LASTONLY) {
+		    call eprintf ("Warning: cannot use ptype LASTONLY with ")
+		    call eprintf ("combined opacity/intensity data.\n")
+		    PTYPE(vp) = P_SUM
+		    call eprintf ("         resetting ptype = %d (SUM)\n")
+			call pargi (PTYPE(vp))
+		}
+		use_both = true
+		if (VERBOSE(vp) == YES)
+		    call eprintf ("4D image, using both opacity & intensity.\n")
+	    } else
+		OPACELEM(vp) = INDEFI
+	} else if (IM_NDIM(im1) > 4)
+	    call error (0, "Don't know how to handle > 4d image")
+
+	# Determine voxel intensity minimum & maximum for all intensity
+	# transformations.  Both a specified intensity min & max and an
+	# image min & max are required in the intensity transformation step
+	# function:  if image min & max are up to date in the image header,
+	# they will be used for image min & max; and if task parameters
+	# imin, imax are NOT supplied, they will be set equal to image min
+	# & max.  Likewise, if image min & max are not present, but
+	# task params imin,imax are, the image min & max will be set to
+	# imin,imax for duration of PVOL execution.  If neither are supplied,
+	# the image min & max will be calculated but not updated (might not
+	# have write access, user might not want them updated); however, if
+	# verbose is on, the user will be warned to run MINMAX on the image
+	# in the future to save time.
+	
+	if (PTYPE(vp) == P_ATTENUATE || use_both) {
+	    # Get opacity transformation function parameters.
+	    if (IS_INDEFR(OMIN(vp)))
+		OMIN(vp) = IIMIN(vp)
+	    if (IS_INDEFR(OMAX(vp)))
+		OMAX(vp) = IIMAX(vp)
+	    if (OMAX(vp) - OMIN(vp) <= 0.0) {
+		call eprintf ("Error:  Invalid omin / omax (%g : %g)\n")
+		    call pargr (OMIN(vp))
+		    call pargr (OMAX(vp))
+		goto iwrapup_
+	    }
+	}
+
+	if (PTYPE(vp) != P_ATTENUATE || use_both) {
+	    # Get intensity transformation function parameters & image minmax.
+	    need_lims = false
+	    if (IM_LIMTIME(im1) < IM_MTIME(im1))
+		need_lims = true
+	    else {
+		tmpmin = IM_MIN(im1)
+		tmpmax = IM_MAX(im1)
+	    }
+	    if (IS_INDEFR(VIMIN(vp))) {
+		if (need_lims) {
+		    call imminmax (im1, tmpmin, tmpmax)
+		    need_lims = false
+		    if (VERBOSE(vp) == YES) {
+			call eprintf ("Must compute input image min & max...\n")
+			call eprintf ("NOTE:  run MINMAX with force+ & update+")
+			call eprintf (" on input image in the future.\n")
+		    }
+		}
+		IIMIN(vp) = tmpmin
+		VIMIN(vp) = IIMIN(vp)
+	    } else {
+		if (need_lims) {
+		    IIMIN(vp) = VIMIN(vp)
+		    if (VERBOSE(vp) == YES)
+			call eprintf ("Image MIN not present; using IMIN\n")
+		} else 
+		    IIMIN(vp) = tmpmin
+	    }
+
+	    if (IS_INDEFR(VIMAX(vp))) {
+		if (need_lims) {
+		    call imminmax (im1, tmpmin, tmpmax)
+		    if (VERBOSE(vp) == YES) {
+			call eprintf ("Must compute input image min & max...\n")
+			call eprintf ("NOTE:  run MINMAX with force+ & update+")
+			call eprintf (" on input image in the future.\n")
+		    }
+		}
+		IIMAX(vp) = tmpmax
+		VIMAX(vp) = IIMAX(vp)
+	    } else {
+		if (need_lims) {
+		    IIMAX(vp) = VIMAX(vp)
+		    if (VERBOSE(vp) == YES)
+			call eprintf ("Image MAX not present; using IMAX\n")
+		} else
+		    IIMAX(vp) = tmpmax
+	    }
+
+	    if (VIMAX(vp) - VIMIN(vp) <= 0.0 && PTYPE(vp) != P_ATTENUATE) {
+		call eprintf ("Error:  Invalid imin / imax (%g : %g)\n")
+		    call pargr (VIMIN(vp))
+		    call pargr (VIMAX(vp))
+		goto iwrapup_
+	    }
+
+	}
+
+	# Load the relevant output header parameters.
+	IM_PIXTYPE(im2) = TY_REAL
+	IM_NDIM(im2) = 3
+	IM_LEN(im2,COL) = IM_LEN(im1,COL)
+
+	# Store run parameters in output image header.
+	call imastr (im2, "V_OIMAGE", Memc[input])
+	call imastr (im2, "V_OTITLE", IM_TITLE(im1))
+	call imaddr (im2, "V_OLDMIN", IIMIN(vp))
+	call imaddr (im2, "V_OLDMAX", IIMAX(vp))
+	call imaddr (im2, "V_DEGREES", DEGREES(vp))
+	call imaddr (im2, "V_THETA0", INIT_THETA(vp))
+	call sprintf (Memc[tmpstr], SZ_LINE, "x=%5.2f, y=%5.2f, z=%5.2f")
+	    call pargr (VECX(vp)); call pargr (VECY(vp)); call pargr (VECZ(vp)) 
+	call imastr (im2, "V_ROTVECT", Memc[tmpstr])
+	call imaddi (im2, "V_PTYPE", PTYPE(vp))
+	call sprintf (Memc[tmpstr], SZ_LINE, "%g : %g")
+	    call pargr (VIMIN(vp)); call pargr (VIMAX(vp))
+	call imastr (im2, "V_IMINMX", Memc[tmpstr])
+	call imaddr (im2, "V_IZERO", IZERO(vp))
+	call sprintf (Memc[tmpstr], SZ_LINE, "%g : %g")
+	    call pargr (OMIN(vp)); call pargr (OMAX(vp))
+	call imastr (im2, "V_OMINMX", Memc[tmpstr])
+	call imaddr (im2, "V_OSCALE", OSCALE(vp))
+	call sprintf (Memc[tmpstr], SZ_LINE, "%g : %g")
+	    call pargr (AMIN(vp)); call pargr (AMAX(vp))
+	call imastr (im2, "V_ATTEN", Memc[tmpstr])
+	if (PTYPE(vp) == P_INVDISPOW || PTYPE(vp) == P_MODN)
+	    call imaddr (im2, "V_DISPOW", DISPOWER(vp))
+	call imaddb (im2, "V_DISCUT", (DISCUTOFF(vp) == YES))
+	if (PTYPE(vp) == P_MODN)
+	    call imaddi (im2, "V_MODN", MODN(vp))
+	if (use_both) {
+	    call imastr (im2, "V_BOTH", "4D: Both opacity and intensity used")
+	    call imaddi (im2, "V_OPELEM", OPACELEM(vp))
+	}
+
+	# Initialize timers.
+	clock1 = clktime (long (0))
+	call cnvtime (clock1, Memc[timestr], SZ_TIME)
+	cpu1 = cputime (long (0))
+
+	# Do all the work.
+	call vproject (im1, im2, vp, use_both)
+
+	call sysid (Memc[tmpstr], SZ_LINE)
+	call imastr (im2, "P_SYSTEM", Memc[tmpstr])
+
+	clock2 = clktime (long (0))
+	elapclock = (clock2 - clock1)
+	cpu2 = cputime (long (0))
+	elapcpu = (cpu2 - cpu1) / 1000
+
+	call imastr (im2, "P_STIME", Memc[timestr])
+	clock1 = clktime (long (0))
+	call cnvtime (clock1, Memc[timestr], SZ_TIME)
+	call imastr (im2, "P_ETIME", Memc[timestr])
+	call sprintf (Memc[tmpstr], SZ_LINE,
+	    "Elapsed cpu = %02d %02s:%02s:%02s, clock = %02d %02s:%02s:%02s")
+	    call pargi (elapcpu/86400)
+	    call pargi (mod (elapcpu, 86400) / 3600)
+	    call pargi (mod (elapcpu, 3600) / 60)
+	    call pargi (mod (elapcpu, 60))
+	    call pargi (elapclock/86400)
+	    call pargi (mod (elapclock, 86400) / 3600)
+	    call pargi (mod (elapclock, 3600) / 60)
+	    call pargi (mod (elapclock, 60))
+	call imastr (im2, "P_ELAPSED", Memc[tmpstr])
+
+iwrapup_
+	call imunmap (im1)
+	call imunmap (im2)
+mwrapup_	
+	call mfree (vp, TY_STRUCT)
+	call sfree (sp)
+end
diff --git a/pkg/proto/vol/src/vgetincr.x b/pkg/proto/vol/src/vgetincr.x
new file mode 100644
index 00000000..c96ff2bb
--- /dev/null
+++ b/pkg/proto/vol/src/vgetincr.x
@@ -0,0 +1,92 @@
+include "pvol.h"
+
+
+# VGETINCREM -- Get list of input voxel band & line indices that contribute to
+# current ray, using simple incremental digital differential analyzer.
+
+procedure vgetincrem (tx1,ty1, tx2,ty2, nx,ny, maxvox, nvox, xind, yind)
+double	tx1,ty1		# (in)  starting coordinate of ray
+double	tx2,ty2		# (in)  ending coordinate of ray
+int	nx,ny		# (in)  dimensions of working plane (1:nx, 1:ny)
+int	maxvox		# (in)  max dimension of output index arrays
+int	nvox		# (out) count of indices for current ray
+int	xind[ARB]	# (out) array of input voxel band indices
+int	yind[ARB]	# (out) array of input voxel line indices
+
+real	x1,y1, x2,y2, dy,dx, adx,ady, x,y, length
+int	i, tvox, xi, yi
+
+int	vsign()
+
+begin
+	# Going between integer and floating point grid representations
+	# is tricky, especially for symmetrical rotation angles aligned with
+	# the grid nodes.  Rounding from double to single precision here
+	# is the only way I could get things to work for all possible angles
+	# and grid dimensions.
+
+	x1 = tx1
+	y1 = ty1
+	x2 = tx2
+	y2 = ty2
+	dx = x2 - x1
+	dy = y2 - y1
+	adx = abs (dx)
+	ady = abs (dy)
+
+	# Approximate the line length.
+	if (adx >= ady)
+	    length = adx
+	else
+	    length = ady
+	tvox = int (length) + 1
+	if (tvox > maxvox)
+	    call error (0, "VGETINCREM:  nvox > maxvox")
+	
+	# Select the larger of dx or dy to be one raster unit.
+	dx = dx / length
+	dy = dy / length
+
+	# Round values; using vsign function makes this work in all quadrants.
+	x = x1 + 0.5 * vsign (dx)
+	y = y1 + 0.5 * vsign (dy)
+
+	# Boundary-extend if coming FROM +x or +y and if rounding would not
+	# take us out of range.
+	if (dx == -1.0 || dy == -1.0) {
+	    if (!((int(x-dx) <= 0 || int(x-dx) > nx) ||
+		(int(y-dy) <= 0 || int(y-dy) > ny))) {
+		x = x - dx
+		y = y - dy
+	    }
+	}
+
+	# Fill in the integer grid coordinates.
+	nvox = 0
+	do i = 1, tvox {
+	    xi = nx - int(x) + 1	# yes, we want to truncate here
+	    yi = int (y)
+	    if (1 <= xi && xi <= nx && 1 <= yi && yi <= ny) {
+		nvox = nvox + 1
+		xind[nvox] = xi
+		yind[nvox] = yi
+	    }
+	    x = x + dx
+	    y = y + dy
+	}
+end
+
+
+# VSIGN -- Return -1, 0, +1 if val is <0, =0, >0.
+
+int procedure vsign (val)
+real	val
+
+begin
+	if (val < 0.0)
+	    return (-1)
+	else if (val == 0.0)
+	    return (0)
+	else
+	    return (1)
+end
diff --git a/pkg/proto/vol/src/vmatrix.x b/pkg/proto/vol/src/vmatrix.x
new file mode 100644
index 00000000..bfc01d63
--- /dev/null
+++ b/pkg/proto/vol/src/vmatrix.x
@@ -0,0 +1,31 @@
+include <imhdr.h>
+include "pvol.h"
+
+
+# VMATRIX -- Volume rotation, rotation matrix projection algorithm.
+# Proceeds from origin at back of volume image toward front, writing
+# output image lines in successive overlapping sheets.  See "Back to 
+# Front Display of Voxel-Based Objects", G.Frieder, D.Gordon, R.Reynolds,
+# IEEE Computer Graphics & Applications Jan. 85, p 52-60.
+
+procedure vmatrix (im1, im2, vp)
+pointer im1		# Input volume image
+pointer	im2		# Output projection image
+pointer	vp		# Volume projection descriptor
+
+real	v, vx, vy, vz
+real	dcosa, dcosb, dcosc
+#real	t11,t21,t31, t12,t22,t32, t13,t23,t33
+
+begin
+	vx = VECX(vp)
+	vy = VECY(vp)
+	vz = VECZ(vp)
+	v = sqrt (vx*vx + vy*vy + vz*vz)
+	dcosa = vx / v
+	dcosb = vy / v
+	dcosc = vz / v
+
+	# ???????
+end
+
diff --git a/pkg/proto/vol/src/vproject.x b/pkg/proto/vol/src/vproject.x
new file mode 100644
index 00000000..009f564d
--- /dev/null
+++ b/pkg/proto/vol/src/vproject.x
@@ -0,0 +1,224 @@
+include <math.h>
+include <imhdr.h>
+include "pvol.h"
+
+define	incr_	91
+
+
+# VPROJECT -- Volume rotation, incremental projection algorithm.
+# Routine attempts to hold as much of the input image in memory as possible.
+# Constructs output image one complete line at a time, determining the
+# contributing voxels for each ray by an incremental rasterizer-like algorithm.
+
+procedure vproject (im1, im2, vp, use_both)
+pointer im1		# Input volume image
+pointer	im2		# Output projection image
+pointer	vp		# Volume projection descriptor
+bool	use_both	# Use both opacity and intensity from 4D image
+
+int	plines, pcols, oline, oband, rot, nvox, oldsize
+int	pnx,pny, len_x, px1,px2, ix,iy,iz,ih, ndims
+real	phalf
+double	rx1,ry1,rx2,ry2, orx1,ory1,orx2,ory2, xc,yc, pdx,pdy, pstep_dx,pstep_dy
+double	astep, theta, theta0, uc_theta
+pointer	sp, vox_x,vox_y, optr, ioptr, buf_in
+bool	first_pass
+long	vs[3], ve[3], ivs[4], ive[4]
+
+pointer	imggss(), imggsi(), imggsl(), imggsr(), imggsd(), imggsx()
+pointer	impgsr()
+
+begin
+	ix = IM_LEN(im1,1)
+	iy = IM_LEN(im1,2)
+	iz = IM_LEN(im1,3)
+	if (use_both) {
+	    ih = 2
+	    ndims = 4
+	} else {
+	    ih = 1
+	    ndims = 3
+	}
+
+	# Set up coordinates for rotation by aligning the center of the working
+	# projection plane ("p-plane") with the center of the volume image.
+
+	pnx = iz			# volume image bands become p-plane X
+	pny = iy			# volume image lines become p-plane Y
+	plines = int (DIS(double(pnx),double(pny)))
+	if (mod (plines, 2) == 0)
+	    plines = plines + 1
+	pcols = IM_LEN(im2,1)
+	phalf = (plines - 1) / 2	# distance from center to bottom pline
+	IM_LEN(im2,2) = plines
+	IM_LEN(im2,3) = NFRAMES(vp)
+	xc = 0.5 * (pnx + 1)
+	yc = 0.5 * (pny + 1)
+
+	# Allocate index arrays for contributing voxels.
+	call smark (sp)
+	call salloc (vox_x, plines, TY_INT)
+	call salloc (vox_y, plines, TY_INT)
+
+	astep = DDEGTORAD (DEGREES(vp))	# angular increment in radians
+	
+	# Determine how much memory we can use, and adjust working set.
+	call pv_gmem (im1, im2, use_both, VERBOSE(vp), MAX_WS(vp), len_x,
+	    oldsize)
+
+	# Read as much of the input image as possible into memory, in column
+	# blocks so we can project through all lines and bands in memory; we
+	# only want to read each voxel of the input image once.
+
+	ivs[2] = 1
+	ive[2] = iy
+	ivs[3] = 1
+	ive[3] = iz
+	ivs[4] = 1
+	ive[4] = 2
+	first_pass = true
+	do px1 = 1, ix, len_x {
+	    px2 = px1 + len_x - 1
+	    if (px2 > ix)
+		px2 = ix
+	    if (VERBOSE(vp) == YES) {
+		call eprintf ("px1=%d, px2=%d, len_x=%d\n")
+		call pargi (px1); call pargi (px2); call pargi (px2-px1+1)
+	    }
+	    ivs[1] = px1
+	    ive[1] = px2
+	    switch (IM_PIXTYPE (im1)) {
+	    case TY_SHORT:
+		buf_in = imggss (im1, ivs, ive, ndims)
+	    case TY_INT:
+		buf_in = imggsi (im1, ivs, ive, ndims)
+	    case TY_LONG:
+		buf_in = imggsl (im1, ivs, ive, ndims)
+	    case TY_REAL:
+		buf_in = imggsr (im1, ivs, ive, ndims)
+	    case TY_DOUBLE:
+		buf_in = imggsd (im1, ivs, ive, ndims)
+	    case TY_COMPLEX:
+		buf_in = imggsx (im1, ivs, ive, ndims)
+	    default:
+		call error (3, "unknown pixel type")
+	    }
+
+	    # Invariant part of output image section:
+	    vs[1] = 1
+	    ve[1] = pcols
+
+	    # Produce one output image band per rotation step around vol image.
+	    theta0 = DDEGTORAD (INIT_THETA(vp))
+	    oband = 1
+	    do rot = 1, NFRAMES(vp) {
+		theta = theta0 + (rot - 1) * astep
+		uc_theta = theta	 # unit-circle for quadrant comparisons.
+		while (uc_theta >= DTWOPI)
+		    uc_theta = uc_theta - DTWOPI
+
+		# Determine line endpoints intersecting the image boundary for
+		# central projection line.
+
+		orx1 = xc - phalf * cos (uc_theta)
+		orx2 = xc + phalf * cos (uc_theta)
+		ory1 = yc - phalf * sin (uc_theta)
+		ory2 = yc + phalf * sin (uc_theta)
+
+		# Offset central projection line to hit the bottom image line of
+		# the projection plane (won't necessarily pass through image).
+
+		pdx = phalf * sin (uc_theta)
+		pdy = phalf * cos (uc_theta)
+		pstep_dx = sin (uc_theta)
+		pstep_dy = cos (uc_theta)
+		orx1 = orx1 + pdx
+		ory1 = ory1 - pdy
+		orx2 = orx2 + pdx
+		ory2 = ory2 - pdy
+		rx1 = orx1
+		ry1 = ory1
+		rx2 = orx2
+		ry2 = ory2
+
+		do oline = 1, plines {
+
+		    # Get voxel indices in p-plane contributing to central ray.
+		    call vgetincrem (rx1,ry1, rx2,ry2, pnx,pny,plines, nvox,
+			Memi[vox_x], Memi[vox_y])
+
+		    # Initialize output line.
+		    vs[2] = oline
+		    ve[2] = oline
+		    vs[3] = oband
+		    ve[3] = oband
+
+		    # If first pass through output image, initialize output
+		    # pixels.  Otherwise, we must read existing part of output
+		    # image into output buffer.
+
+		    if (first_pass) {
+			optr = impgsr (im2, vs, ve, 3)
+
+			# If opacity model, initialize output to incident int.
+			if (PTYPE(vp) == P_ATTENUATE)
+			    call amovkr (IZERO(vp), Memr[optr], pcols)
+			else
+			    call aclrr (Memr[optr], pcols)
+		    } else {
+			ioptr = imggsr (im2, vs, ve, 3)
+			optr = impgsr (im2, vs, ve, 3)
+			call amovr (Memr[ioptr], Memr[optr], pcols)
+		    }
+
+		    # Project each contributing voxel into output image line.
+		    if (nvox > 0)
+			switch (IM_PIXTYPE (im1)) {
+			case TY_SHORT:
+			    call vtransmits (Mems[buf_in], (px2-px1+1),iy,iz,ih,
+				px1,px2, Memi[vox_y], Memi[vox_x], nvox,
+				Memr[optr], vp)
+			case TY_INT:
+			    call vtransmiti (Memi[buf_in], (px2-px1+1),iy,iz,ih,
+				px1,px2, Memi[vox_y], Memi[vox_x], nvox,
+				Memr[optr], vp)
+			case TY_LONG:
+			    call vtransmitl (Meml[buf_in], (px2-px1+1),iy,iz,ih,
+				px1,px2, Memi[vox_y], Memi[vox_x], nvox,
+				Memr[optr], vp)
+			case TY_REAL:
+			    call vtransmitr (Memr[buf_in], (px2-px1+1),iy,iz,ih,
+				px1,px2, Memi[vox_y], Memi[vox_x], nvox,
+				Memr[optr], vp)
+			case TY_DOUBLE:
+			    call vtransmitd (Memd[buf_in], (px2-px1+1),iy,iz,ih,
+				px1,px2, Memi[vox_y], Memi[vox_x], nvox,
+				Memr[optr], vp)
+			case TY_COMPLEX:
+			    call vtransmitx (Memx[buf_in], (px2-px1+1),iy,iz,ih,
+				px1,px2, Memi[vox_y], Memi[vox_x], nvox,
+				Memr[optr], vp)
+			}
+
+		    # Offset endpoints for next projection line.
+		    rx1 = orx1 - oline * pstep_dx
+		    ry1 = ory1 + oline * pstep_dy
+		    rx2 = orx2 - oline * pstep_dx
+		    ry2 = ory2 + oline * pstep_dy
+		}
+
+		# Set up for next rotation.
+		oband = oband + 1
+		if (VERBOSE(vp) == YES) {
+		    call eprintf ("...end of rotation %d, theta %7.2f\n")
+		    call pargi (rot); call pargd (DRADTODEG(theta))
+		}
+	    }
+
+	    first_pass = false
+	    call imflush (im2)
+	}
+
+	call sfree (sp)
+	call fixmem (oldsize)
+end
diff --git a/pkg/proto/vol/src/vtransmit.gx b/pkg/proto/vol/src/vtransmit.gx
new file mode 100644
index 00000000..698d73d0
--- /dev/null
+++ b/pkg/proto/vol/src/vtransmit.gx
@@ -0,0 +1,146 @@
+include <imhdr.h>
+include "pvol.h"
+
+$for (silrdx)
+
+# VTRANSMIT -- Compute the intensities of each output image pixel in the
+# current line as a function of its existing intensity plus the emission
+# and absorption from each contributing voxel.
+
+procedure vtransmit$t (inbuf,nx,ny,nz,nh, px1,px2, iline,iband,nvox, oline, vp)
+PIXEL	inbuf[nx,ny,nz,nh]	# Input data buffer for current set of yz slices
+int	nx,ny,nz,nh		# Dimensions of current input buffer
+int	px1,px2			# Range of columns in current yz slice set
+int	iline[nvox]		# Input image lines for current projection ray
+int	iband[nvox]		# Input image bands for current projection ray
+int	nvox			# Number of voxels in current projection column
+real	oline[ARB]		# output image line buffer
+pointer	vp			# Volume projection descriptor
+
+bool	use_both
+
+int	i, vox, opelem, intelem, frontvox, backvox
+real	amin, amax, vox_op, vox_int, ival, attenuate, vimin, ifac, ofac, distwt
+
+begin
+	# Dereference most frequently used structure elements.
+	amin = AMIN(vp)
+	amax = AMAX(vp)
+	vimin = VIMIN(vp)
+
+	intelem = 1
+	opelem = OPACELEM(vp)
+	if (nh > 1) {
+	    use_both = true
+	    if (opelem == 1)
+		intelem = 2
+	    else if (IS_INDEFI(opelem))
+		opelem = 2
+	} else {
+	    use_both = false
+	    opelem = 1
+	}
+
+
+	# Set up for opacity, intensity, or both.
+	ifac = (IIMAX(vp) - IIMIN(vp)) / (VIMAX(vp) - vimin)
+	if (PTYPE(vp) == P_ATTENUATE || use_both)
+	    ofac = (amax - amin) / (OMAX(vp) - OMIN(vp))
+
+	# Since we are in memory anyway, it is more convenient to traverse
+	# the columns in the outer loop and the voxels from different bands
+	# and lines in the inner loop.  This is necessary when distance
+	# weighting and the distance cutoff option is on (we need to know
+	# the range of usable voxels in a given column before projecting).
+
+	if (PTYPE(vp) == P_INVDISPOW || PTYPE(vp) == P_MODN)
+	    do i = px1, px2 {
+		if (DISCUTOFF(vp) == NO) {
+		    frontvox = nvox
+		    backvox = 1
+		} else {
+		    frontvox = 1
+		    backvox = nvox
+		    do vox = 1, nvox {
+			vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],intelem]
+			if (vimin <= vox_int && vox_int < VIMAX(vp)) {
+			    frontvox = max (frontvox, vox)
+			    backvox = min (backvox, vox)
+			}
+		    }
+		}
+		if (frontvox - backvox < 0)
+		    next
+		do vox = backvox, frontvox {
+		    distwt = (real(vox-backvox+1) /
+			real(frontvox-backvox+1)) ** DISPOWER(vp)
+
+		    # Opacity transformation function.
+		    if (use_both) {
+			vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
+			    opelem] * OSCALE(vp)
+			if (vox_op < OMIN(vp))
+			    attenuate = amax
+			else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
+			    attenuate = amax - (vox_op - OMIN(vp)) * ofac
+			else
+			    attenuate = amin
+			oline[i] = oline[i] * attenuate
+		    }
+
+		    # Intensity transformation function.
+		    vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
+			intelem]
+		    if (vox_int < vimin)
+			ival = IIMIN(vp)
+		    else if (vimin <= vox_int && vox_int < VIMAX(vp))
+			ival = IIMIN(vp) + (vox_int - vimin) * ifac
+		    else
+			ival = IIMAX(vp)
+
+		    if (PTYPE(vp) == P_INVDISPOW)
+			oline[i] = oline[i] + ival * distwt
+		    else if (mod (int (ival/(IIMAX(vp)-IIMIN(vp)) * 100.0),
+			    MODN(vp)) == 0)
+			    oline[i] = oline[i] + ival * distwt
+		}
+	    }
+	else
+	    do i = px1, px2
+		do vox = 1, nvox {
+		    # Opacity transformation function.
+		    if (PTYPE(vp) == P_ATTENUATE || use_both) {
+			vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
+			    opelem] * OSCALE(vp)
+			if (vox_op < OMIN(vp))
+			    attenuate = amax
+			else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
+			    attenuate = amax - (vox_op - OMIN(vp)) * ofac
+			else
+			    attenuate = amin
+			oline[i] = oline[i] * attenuate
+		    }
+
+		    # Intensity transformation function.
+		    if (PTYPE(vp) != P_ATTENUATE || use_both) {
+			vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
+			    intelem]
+			if (vox_int < vimin)
+			    ival = IIMIN(vp)
+			else if (vimin <= vox_int && vox_int < VIMAX(vp))
+			    ival = IIMIN(vp) + (vox_int - vimin) * ifac
+			else
+			    ival = IIMAX(vp)
+
+			if (PTYPE(vp) == P_AVERAGE)
+			    oline[i] = oline[i] + ival * 1.0 / real(nvox)
+			else if (PTYPE(vp) == P_SUM)
+			    oline[i] = oline[i] + ival
+			else if (PTYPE(vp) == P_LASTONLY)
+			    if (ival > 0.0)
+				oline[i] = ival
+		    }
+		}
+end
+
+$endfor
diff --git a/pkg/proto/vol/src/vtransmit.x b/pkg/proto/vol/src/vtransmit.x
new file mode 100644
index 00000000..99716a2d
--- /dev/null
+++ b/pkg/proto/vol/src/vtransmit.x
@@ -0,0 +1,856 @@
+include <imhdr.h>
+include "pvol.h"
+
+
+
+# VTRANSMIT -- Compute the intensities of each output image pixel in the
+# current line as a function of its existing intensity plus the emission
+# and absorption from each contributing voxel.
+
+procedure vtransmits (inbuf,nx,ny,nz,nh, px1,px2, iline,iband,nvox, oline, vp)
+short	inbuf[nx,ny,nz,nh]	# Input data buffer for current set of yz slices
+int	nx,ny,nz,nh		# Dimensions of current input buffer
+int	px1,px2			# Range of columns in current yz slice set
+int	iline[nvox]		# Input image lines for current projection ray
+int	iband[nvox]		# Input image bands for current projection ray
+int	nvox			# Number of voxels in current projection column
+real	oline[ARB]		# output image line buffer
+pointer	vp			# Volume projection descriptor
+
+bool	use_both
+
+int	i, vox, opelem, intelem, frontvox, backvox
+real	amin, amax, vox_op, vox_int, ival, attenuate, vimin, ifac, ofac, distwt
+
+begin
+	# Dereference most frequently used structure elements.
+	amin = AMIN(vp)
+	amax = AMAX(vp)
+	vimin = VIMIN(vp)
+
+	intelem = 1
+	opelem = OPACELEM(vp)
+	if (nh > 1) {
+	    use_both = true
+	    if (opelem == 1)
+		intelem = 2
+	    else if (IS_INDEFI(opelem))
+		opelem = 2
+	} else {
+	    use_both = false
+	    opelem = 1
+	}
+
+
+	# Set up for opacity, intensity, or both.
+	ifac = (IIMAX(vp) - IIMIN(vp)) / (VIMAX(vp) - vimin)
+	if (PTYPE(vp) == P_ATTENUATE || use_both)
+	    ofac = (amax - amin) / (OMAX(vp) - OMIN(vp))
+
+	# Since we are in memory anyway, it is more convenient to traverse
+	# the columns in the outer loop and the voxels from different bands
+	# and lines in the inner loop.  This is necessary when distance
+	# weighting and the distance cutoff option is on (we need to know
+	# the range of usable voxels in a given column before projecting).
+
+	if (PTYPE(vp) == P_INVDISPOW || PTYPE(vp) == P_MODN)
+	    do i = px1, px2 {
+		if (DISCUTOFF(vp) == NO) {
+		    frontvox = nvox
+		    backvox = 1
+		} else {
+		    frontvox = 1
+		    backvox = nvox
+		    do vox = 1, nvox {
+			vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],intelem]
+			if (vimin <= vox_int && vox_int < VIMAX(vp)) {
+			    frontvox = max (frontvox, vox)
+			    backvox = min (backvox, vox)
+			}
+		    }
+		}
+		if (frontvox - backvox < 0)
+		    next
+		do vox = backvox, frontvox {
+		    distwt = (real(vox-backvox+1) /
+			real(frontvox-backvox+1)) ** DISPOWER(vp)
+
+		    # Opacity transformation function.
+		    if (use_both) {
+			vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
+			    opelem] * OSCALE(vp)
+			if (vox_op < OMIN(vp))
+			    attenuate = amax
+			else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
+			    attenuate = amax - (vox_op - OMIN(vp)) * ofac
+			else
+			    attenuate = amin
+			oline[i] = oline[i] * attenuate
+		    }
+
+		    # Intensity transformation function.
+		    vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
+			intelem]
+		    if (vox_int < vimin)
+			ival = IIMIN(vp)
+		    else if (vimin <= vox_int && vox_int < VIMAX(vp))
+			ival = IIMIN(vp) + (vox_int - vimin) * ifac
+		    else
+			ival = IIMAX(vp)
+
+		    if (PTYPE(vp) == P_INVDISPOW)
+			oline[i] = oline[i] + ival * distwt
+		    else if (mod (int (ival/(IIMAX(vp)-IIMIN(vp)) * 100.0),
+			    MODN(vp)) == 0)
+			    oline[i] = oline[i] + ival * distwt
+		}
+	    }
+	else
+	    do i = px1, px2
+		do vox = 1, nvox {
+		    # Opacity transformation function.
+		    if (PTYPE(vp) == P_ATTENUATE || use_both) {
+			vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
+			    opelem] * OSCALE(vp)
+			if (vox_op < OMIN(vp))
+			    attenuate = amax
+			else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
+			    attenuate = amax - (vox_op - OMIN(vp)) * ofac
+			else
+			    attenuate = amin
+			oline[i] = oline[i] * attenuate
+		    }
+
+		    # Intensity transformation function.
+		    if (PTYPE(vp) != P_ATTENUATE || use_both) {
+			vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
+			    intelem]
+			if (vox_int < vimin)
+			    ival = IIMIN(vp)
+			else if (vimin <= vox_int && vox_int < VIMAX(vp))
+			    ival = IIMIN(vp) + (vox_int - vimin) * ifac
+			else
+			    ival = IIMAX(vp)
+
+			if (PTYPE(vp) == P_AVERAGE)
+			    oline[i] = oline[i] + ival * 1.0 / real(nvox)
+			else if (PTYPE(vp) == P_SUM)
+			    oline[i] = oline[i] + ival
+			else if (PTYPE(vp) == P_LASTONLY)
+			    if (ival > 0.0)
+				oline[i] = ival
+		    }
+		}
+end
+
+
+
+# VTRANSMIT -- Compute the intensities of each output image pixel in the
+# current line as a function of its existing intensity plus the emission
+# and absorption from each contributing voxel.
+
+procedure vtransmiti (inbuf,nx,ny,nz,nh, px1,px2, iline,iband,nvox, oline, vp)
+int	inbuf[nx,ny,nz,nh]	# Input data buffer for current set of yz slices
+int	nx,ny,nz,nh		# Dimensions of current input buffer
+int	px1,px2			# Range of columns in current yz slice set
+int	iline[nvox]		# Input image lines for current projection ray
+int	iband[nvox]		# Input image bands for current projection ray
+int	nvox			# Number of voxels in current projection column
+real	oline[ARB]		# output image line buffer
+pointer	vp			# Volume projection descriptor
+
+bool	use_both
+
+int	i, vox, opelem, intelem, frontvox, backvox
+real	amin, amax, vox_op, vox_int, ival, attenuate, vimin, ifac, ofac, distwt
+
+begin
+	# Dereference most frequently used structure elements.
+	amin = AMIN(vp)
+	amax = AMAX(vp)
+	vimin = VIMIN(vp)
+
+	intelem = 1
+	opelem = OPACELEM(vp)
+	if (nh > 1) {
+	    use_both = true
+	    if (opelem == 1)
+		intelem = 2
+	    else if (IS_INDEFI(opelem))
+		opelem = 2
+	} else {
+	    use_both = false
+	    opelem = 1
+	}
+
+
+	# Set up for opacity, intensity, or both.
+	ifac = (IIMAX(vp) - IIMIN(vp)) / (VIMAX(vp) - vimin)
+	if (PTYPE(vp) == P_ATTENUATE || use_both)
+	    ofac = (amax - amin) / (OMAX(vp) - OMIN(vp))
+
+	# Since we are in memory anyway, it is more convenient to traverse
+	# the columns in the outer loop and the voxels from different bands
+	# and lines in the inner loop.  This is necessary when distance
+	# weighting and the distance cutoff option is on (we need to know
+	# the range of usable voxels in a given column before projecting).
+
+	if (PTYPE(vp) == P_INVDISPOW || PTYPE(vp) == P_MODN)
+	    do i = px1, px2 {
+		if (DISCUTOFF(vp) == NO) {
+		    frontvox = nvox
+		    backvox = 1
+		} else {
+		    frontvox = 1
+		    backvox = nvox
+		    do vox = 1, nvox {
+			vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],intelem]
+			if (vimin <= vox_int && vox_int < VIMAX(vp)) {
+			    frontvox = max (frontvox, vox)
+			    backvox = min (backvox, vox)
+			}
+		    }
+		}
+		if (frontvox - backvox < 0)
+		    next
+		do vox = backvox, frontvox {
+		    distwt = (real(vox-backvox+1) /
+			real(frontvox-backvox+1)) ** DISPOWER(vp)
+
+		    # Opacity transformation function.
+		    if (use_both) {
+			vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
+			    opelem] * OSCALE(vp)
+			if (vox_op < OMIN(vp))
+			    attenuate = amax
+			else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
+			    attenuate = amax - (vox_op - OMIN(vp)) * ofac
+			else
+			    attenuate = amin
+			oline[i] = oline[i] * attenuate
+		    }
+
+		    # Intensity transformation function.
+		    vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
+			intelem]
+		    if (vox_int < vimin)
+			ival = IIMIN(vp)
+		    else if (vimin <= vox_int && vox_int < VIMAX(vp))
+			ival = IIMIN(vp) + (vox_int - vimin) * ifac
+		    else
+			ival = IIMAX(vp)
+
+		    if (PTYPE(vp) == P_INVDISPOW)
+			oline[i] = oline[i] + ival * distwt
+		    else if (mod (int (ival/(IIMAX(vp)-IIMIN(vp)) * 100.0),
+			    MODN(vp)) == 0)
+			    oline[i] = oline[i] + ival * distwt
+		}
+	    }
+	else
+	    do i = px1, px2
+		do vox = 1, nvox {
+		    # Opacity transformation function.
+		    if (PTYPE(vp) == P_ATTENUATE || use_both) {
+			vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
+			    opelem] * OSCALE(vp)
+			if (vox_op < OMIN(vp))
+			    attenuate = amax
+			else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
+			    attenuate = amax - (vox_op - OMIN(vp)) * ofac
+			else
+			    attenuate = amin
+			oline[i] = oline[i] * attenuate
+		    }
+
+		    # Intensity transformation function.
+		    if (PTYPE(vp) != P_ATTENUATE || use_both) {
+			vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
+			    intelem]
+			if (vox_int < vimin)
+			    ival = IIMIN(vp)
+			else if (vimin <= vox_int && vox_int < VIMAX(vp))
+			    ival = IIMIN(vp) + (vox_int - vimin) * ifac
+			else
+			    ival = IIMAX(vp)
+
+			if (PTYPE(vp) == P_AVERAGE)
+			    oline[i] = oline[i] + ival * 1.0 / real(nvox)
+			else if (PTYPE(vp) == P_SUM)
+			    oline[i] = oline[i] + ival
+			else if (PTYPE(vp) == P_LASTONLY)
+			    if (ival > 0.0)
+				oline[i] = ival
+		    }
+		}
+end
+
+
+
+# VTRANSMIT -- Compute the intensities of each output image pixel in the
+# current line as a function of its existing intensity plus the emission
+# and absorption from each contributing voxel.
+
+procedure vtransmitl (inbuf,nx,ny,nz,nh, px1,px2, iline,iband,nvox, oline, vp)
+long	inbuf[nx,ny,nz,nh]	# Input data buffer for current set of yz slices
+int	nx,ny,nz,nh		# Dimensions of current input buffer
+int	px1,px2			# Range of columns in current yz slice set
+int	iline[nvox]		# Input image lines for current projection ray
+int	iband[nvox]		# Input image bands for current projection ray
+int	nvox			# Number of voxels in current projection column
+real	oline[ARB]		# output image line buffer
+pointer	vp			# Volume projection descriptor
+
+bool	use_both
+
+int	i, vox, opelem, intelem, frontvox, backvox
+real	amin, amax, vox_op, vox_int, ival, attenuate, vimin, ifac, ofac, distwt
+
+begin
+	# Dereference most frequently used structure elements.
+	amin = AMIN(vp)
+	amax = AMAX(vp)
+	vimin = VIMIN(vp)
+
+	intelem = 1
+	opelem = OPACELEM(vp)
+	if (nh > 1) {
+	    use_both = true
+	    if (opelem == 1)
+		intelem = 2
+	    else if (IS_INDEFI(opelem))
+		opelem = 2
+	} else {
+	    use_both = false
+	    opelem = 1
+	}
+
+
+	# Set up for opacity, intensity, or both.
+	ifac = (IIMAX(vp) - IIMIN(vp)) / (VIMAX(vp) - vimin)
+	if (PTYPE(vp) == P_ATTENUATE || use_both)
+	    ofac = (amax - amin) / (OMAX(vp) - OMIN(vp))
+
+	# Since we are in memory anyway, it is more convenient to traverse
+	# the columns in the outer loop and the voxels from different bands
+	# and lines in the inner loop.  This is necessary when distance
+	# weighting and the distance cutoff option is on (we need to know
+	# the range of usable voxels in a given column before projecting).
+
+	if (PTYPE(vp) == P_INVDISPOW || PTYPE(vp) == P_MODN)
+	    do i = px1, px2 {
+		if (DISCUTOFF(vp) == NO) {
+		    frontvox = nvox
+		    backvox = 1
+		} else {
+		    frontvox = 1
+		    backvox = nvox
+		    do vox = 1, nvox {
+			vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],intelem]
+			if (vimin <= vox_int && vox_int < VIMAX(vp)) {
+			    frontvox = max (frontvox, vox)
+			    backvox = min (backvox, vox)
+			}
+		    }
+		}
+		if (frontvox - backvox < 0)
+		    next
+		do vox = backvox, frontvox {
+		    distwt = (real(vox-backvox+1) /
+			real(frontvox-backvox+1)) ** DISPOWER(vp)
+
+		    # Opacity transformation function.
+		    if (use_both) {
+			vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
+			    opelem] * OSCALE(vp)
+			if (vox_op < OMIN(vp))
+			    attenuate = amax
+			else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
+			    attenuate = amax - (vox_op - OMIN(vp)) * ofac
+			else
+			    attenuate = amin
+			oline[i] = oline[i] * attenuate
+		    }
+
+		    # Intensity transformation function.
+		    vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
+			intelem]
+		    if (vox_int < vimin)
+			ival = IIMIN(vp)
+		    else if (vimin <= vox_int && vox_int < VIMAX(vp))
+			ival = IIMIN(vp) + (vox_int - vimin) * ifac
+		    else
+			ival = IIMAX(vp)
+
+		    if (PTYPE(vp) == P_INVDISPOW)
+			oline[i] = oline[i] + ival * distwt
+		    else if (mod (int (ival/(IIMAX(vp)-IIMIN(vp)) * 100.0),
+			    MODN(vp)) == 0)
+			    oline[i] = oline[i] + ival * distwt
+		}
+	    }
+	else
+	    do i = px1, px2
+		do vox = 1, nvox {
+		    # Opacity transformation function.
+		    if (PTYPE(vp) == P_ATTENUATE || use_both) {
+			vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
+			    opelem] * OSCALE(vp)
+			if (vox_op < OMIN(vp))
+			    attenuate = amax
+			else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
+			    attenuate = amax - (vox_op - OMIN(vp)) * ofac
+			else
+			    attenuate = amin
+			oline[i] = oline[i] * attenuate
+		    }
+
+		    # Intensity transformation function.
+		    if (PTYPE(vp) != P_ATTENUATE || use_both) {
+			vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
+			    intelem]
+			if (vox_int < vimin)
+			    ival = IIMIN(vp)
+			else if (vimin <= vox_int && vox_int < VIMAX(vp))
+			    ival = IIMIN(vp) + (vox_int - vimin) * ifac
+			else
+			    ival = IIMAX(vp)
+
+			if (PTYPE(vp) == P_AVERAGE)
+			    oline[i] = oline[i] + ival * 1.0 / real(nvox)
+			else if (PTYPE(vp) == P_SUM)
+			    oline[i] = oline[i] + ival
+			else if (PTYPE(vp) == P_LASTONLY)
+			    if (ival > 0.0)
+				oline[i] = ival
+		    }
+		}
+end
+
+
+
+# VTRANSMIT -- Compute the intensities of each output image pixel in the
+# current line as a function of its existing intensity plus the emission
+# and absorption from each contributing voxel.
+
+procedure vtransmitr (inbuf,nx,ny,nz,nh, px1,px2, iline,iband,nvox, oline, vp)
+real	inbuf[nx,ny,nz,nh]	# Input data buffer for current set of yz slices
+int	nx,ny,nz,nh		# Dimensions of current input buffer
+int	px1,px2			# Range of columns in current yz slice set
+int	iline[nvox]		# Input image lines for current projection ray
+int	iband[nvox]		# Input image bands for current projection ray
+int	nvox			# Number of voxels in current projection column
+real	oline[ARB]		# output image line buffer
+pointer	vp			# Volume projection descriptor
+
+bool	use_both
+
+int	i, vox, opelem, intelem, frontvox, backvox
+real	amin, amax, vox_op, vox_int, ival, attenuate, vimin, ifac, ofac, distwt
+
+begin
+	# Dereference most frequently used structure elements.
+	amin = AMIN(vp)
+	amax = AMAX(vp)
+	vimin = VIMIN(vp)
+
+	intelem = 1
+	opelem = OPACELEM(vp)
+	if (nh > 1) {
+	    use_both = true
+	    if (opelem == 1)
+		intelem = 2
+	    else if (IS_INDEFI(opelem))
+		opelem = 2
+	} else {
+	    use_both = false
+	    opelem = 1
+	}
+
+
+	# Set up for opacity, intensity, or both.
+	ifac = (IIMAX(vp) - IIMIN(vp)) / (VIMAX(vp) - vimin)
+	if (PTYPE(vp) == P_ATTENUATE || use_both)
+	    ofac = (amax - amin) / (OMAX(vp) - OMIN(vp))
+
+	# Since we are in memory anyway, it is more convenient to traverse
+	# the columns in the outer loop and the voxels from different bands
+	# and lines in the inner loop.  This is necessary when distance
+	# weighting and the distance cutoff option is on (we need to know
+	# the range of usable voxels in a given column before projecting).
+
+	if (PTYPE(vp) == P_INVDISPOW || PTYPE(vp) == P_MODN)
+	    do i = px1, px2 {
+		if (DISCUTOFF(vp) == NO) {
+		    frontvox = nvox
+		    backvox = 1
+		} else {
+		    frontvox = 1
+		    backvox = nvox
+		    do vox = 1, nvox {
+			vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],intelem]
+			if (vimin <= vox_int && vox_int < VIMAX(vp)) {
+			    frontvox = max (frontvox, vox)
+			    backvox = min (backvox, vox)
+			}
+		    }
+		}
+		if (frontvox - backvox < 0)
+		    next
+		do vox = backvox, frontvox {
+		    distwt = (real(vox-backvox+1) /
+			real(frontvox-backvox+1)) ** DISPOWER(vp)
+
+		    # Opacity transformation function.
+		    if (use_both) {
+			vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
+			    opelem] * OSCALE(vp)
+			if (vox_op < OMIN(vp))
+			    attenuate = amax
+			else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
+			    attenuate = amax - (vox_op - OMIN(vp)) * ofac
+			else
+			    attenuate = amin
+			oline[i] = oline[i] * attenuate
+		    }
+
+		    # Intensity transformation function.
+		    vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
+			intelem]
+		    if (vox_int < vimin)
+			ival = IIMIN(vp)
+		    else if (vimin <= vox_int && vox_int < VIMAX(vp))
+			ival = IIMIN(vp) + (vox_int - vimin) * ifac
+		    else
+			ival = IIMAX(vp)
+
+		    if (PTYPE(vp) == P_INVDISPOW)
+			oline[i] = oline[i] + ival * distwt
+		    else if (mod (int (ival/(IIMAX(vp)-IIMIN(vp)) * 100.0),
+			    MODN(vp)) == 0)
+			    oline[i] = oline[i] + ival * distwt
+		}
+	    }
+	else
+	    do i = px1, px2
+		do vox = 1, nvox {
+		    # Opacity transformation function.
+		    if (PTYPE(vp) == P_ATTENUATE || use_both) {
+			vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
+			    opelem] * OSCALE(vp)
+			if (vox_op < OMIN(vp))
+			    attenuate = amax
+			else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
+			    attenuate = amax - (vox_op - OMIN(vp)) * ofac
+			else
+			    attenuate = amin
+			oline[i] = oline[i] * attenuate
+		    }
+
+		    # Intensity transformation function.
+		    if (PTYPE(vp) != P_ATTENUATE || use_both) {
+			vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
+			    intelem]
+			if (vox_int < vimin)
+			    ival = IIMIN(vp)
+			else if (vimin <= vox_int && vox_int < VIMAX(vp))
+			    ival = IIMIN(vp) + (vox_int - vimin) * ifac
+			else
+			    ival = IIMAX(vp)
+
+			if (PTYPE(vp) == P_AVERAGE)
+			    oline[i] = oline[i] + ival * 1.0 / real(nvox)
+			else if (PTYPE(vp) == P_SUM)
+			    oline[i] = oline[i] + ival
+			else if (PTYPE(vp) == P_LASTONLY)
+			    if (ival > 0.0)
+				oline[i] = ival
+		    }
+		}
+end
+
+
+
+# VTRANSMIT -- Compute the intensities of each output image pixel in the
+# current line as a function of its existing intensity plus the emission
+# and absorption from each contributing voxel.
+
+procedure vtransmitd (inbuf,nx,ny,nz,nh, px1,px2, iline,iband,nvox, oline, vp)
+double	inbuf[nx,ny,nz,nh]	# Input data buffer for current set of yz slices
+int	nx,ny,nz,nh		# Dimensions of current input buffer
+int	px1,px2			# Range of columns in current yz slice set
+int	iline[nvox]		# Input image lines for current projection ray
+int	iband[nvox]		# Input image bands for current projection ray
+int	nvox			# Number of voxels in current projection column
+real	oline[ARB]		# output image line buffer
+pointer	vp			# Volume projection descriptor
+
+bool	use_both
+
+int	i, vox, opelem, intelem, frontvox, backvox
+real	amin, amax, vox_op, vox_int, ival, attenuate, vimin, ifac, ofac, distwt
+
+begin
+	# Dereference most frequently used structure elements.
+	amin = AMIN(vp)
+	amax = AMAX(vp)
+	vimin = VIMIN(vp)
+
+	intelem = 1
+	opelem = OPACELEM(vp)
+	if (nh > 1) {
+	    use_both = true
+	    if (opelem == 1)
+		intelem = 2
+	    else if (IS_INDEFI(opelem))
+		opelem = 2
+	} else {
+	    use_both = false
+	    opelem = 1
+	}
+
+
+	# Set up for opacity, intensity, or both.
+	ifac = (IIMAX(vp) - IIMIN(vp)) / (VIMAX(vp) - vimin)
+	if (PTYPE(vp) == P_ATTENUATE || use_both)
+	    ofac = (amax - amin) / (OMAX(vp) - OMIN(vp))
+
+	# Since we are in memory anyway, it is more convenient to traverse
+	# the columns in the outer loop and the voxels from different bands
+	# and lines in the inner loop.  This is necessary when distance
+	# weighting and the distance cutoff option is on (we need to know
+	# the range of usable voxels in a given column before projecting).
+
+	if (PTYPE(vp) == P_INVDISPOW || PTYPE(vp) == P_MODN)
+	    do i = px1, px2 {
+		if (DISCUTOFF(vp) == NO) {
+		    frontvox = nvox
+		    backvox = 1
+		} else {
+		    frontvox = 1
+		    backvox = nvox
+		    do vox = 1, nvox {
+			vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],intelem]
+			if (vimin <= vox_int && vox_int < VIMAX(vp)) {
+			    frontvox = max (frontvox, vox)
+			    backvox = min (backvox, vox)
+			}
+		    }
+		}
+		if (frontvox - backvox < 0)
+		    next
+		do vox = backvox, frontvox {
+		    distwt = (real(vox-backvox+1) /
+			real(frontvox-backvox+1)) ** DISPOWER(vp)
+
+		    # Opacity transformation function.
+		    if (use_both) {
+			vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
+			    opelem] * OSCALE(vp)
+			if (vox_op < OMIN(vp))
+			    attenuate = amax
+			else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
+			    attenuate = amax - (vox_op - OMIN(vp)) * ofac
+			else
+			    attenuate = amin
+			oline[i] = oline[i] * attenuate
+		    }
+
+		    # Intensity transformation function.
+		    vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
+			intelem]
+		    if (vox_int < vimin)
+			ival = IIMIN(vp)
+		    else if (vimin <= vox_int && vox_int < VIMAX(vp))
+			ival = IIMIN(vp) + (vox_int - vimin) * ifac
+		    else
+			ival = IIMAX(vp)
+
+		    if (PTYPE(vp) == P_INVDISPOW)
+			oline[i] = oline[i] + ival * distwt
+		    else if (mod (int (ival/(IIMAX(vp)-IIMIN(vp)) * 100.0),
+			    MODN(vp)) == 0)
+			    oline[i] = oline[i] + ival * distwt
+		}
+	    }
+	else
+	    do i = px1, px2
+		do vox = 1, nvox {
+		    # Opacity transformation function.
+		    if (PTYPE(vp) == P_ATTENUATE || use_both) {
+			vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
+			    opelem] * OSCALE(vp)
+			if (vox_op < OMIN(vp))
+			    attenuate = amax
+			else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
+			    attenuate = amax - (vox_op - OMIN(vp)) * ofac
+			else
+			    attenuate = amin
+			oline[i] = oline[i] * attenuate
+		    }
+
+		    # Intensity transformation function.
+		    if (PTYPE(vp) != P_ATTENUATE || use_both) {
+			vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
+			    intelem]
+			if (vox_int < vimin)
+			    ival = IIMIN(vp)
+			else if (vimin <= vox_int && vox_int < VIMAX(vp))
+			    ival = IIMIN(vp) + (vox_int - vimin) * ifac
+			else
+			    ival = IIMAX(vp)
+
+			if (PTYPE(vp) == P_AVERAGE)
+			    oline[i] = oline[i] + ival * 1.0 / real(nvox)
+			else if (PTYPE(vp) == P_SUM)
+			    oline[i] = oline[i] + ival
+			else if (PTYPE(vp) == P_LASTONLY)
+			    if (ival > 0.0)
+				oline[i] = ival
+		    }
+		}
+end
+
+
+
+# VTRANSMIT -- Compute the intensities of each output image pixel in the
+# current line as a function of its existing intensity plus the emission
+# and absorption from each contributing voxel.
+
+procedure vtransmitx (inbuf,nx,ny,nz,nh, px1,px2, iline,iband,nvox, oline, vp)
+complex	inbuf[nx,ny,nz,nh]	# Input data buffer for current set of yz slices
+int	nx,ny,nz,nh		# Dimensions of current input buffer
+int	px1,px2			# Range of columns in current yz slice set
+int	iline[nvox]		# Input image lines for current projection ray
+int	iband[nvox]		# Input image bands for current projection ray
+int	nvox			# Number of voxels in current projection column
+real	oline[ARB]		# output image line buffer
+pointer	vp			# Volume projection descriptor
+
+bool	use_both
+
+int	i, vox, opelem, intelem, frontvox, backvox
+real	amin, amax, vox_op, vox_int, ival, attenuate, vimin, ifac, ofac, distwt
+
+begin
+	# Dereference most frequently used structure elements.
+	amin = AMIN(vp)
+	amax = AMAX(vp)
+	vimin = VIMIN(vp)
+
+	intelem = 1
+	opelem = OPACELEM(vp)
+	if (nh > 1) {
+	    use_both = true
+	    if (opelem == 1)
+		intelem = 2
+	    else if (IS_INDEFI(opelem))
+		opelem = 2
+	} else {
+	    use_both = false
+	    opelem = 1
+	}
+
+
+	# Set up for opacity, intensity, or both.
+	ifac = (IIMAX(vp) - IIMIN(vp)) / (VIMAX(vp) - vimin)
+	if (PTYPE(vp) == P_ATTENUATE || use_both)
+	    ofac = (amax - amin) / (OMAX(vp) - OMIN(vp))
+
+	# Since we are in memory anyway, it is more convenient to traverse
+	# the columns in the outer loop and the voxels from different bands
+	# and lines in the inner loop.  This is necessary when distance
+	# weighting and the distance cutoff option is on (we need to know
+	# the range of usable voxels in a given column before projecting).
+
+	if (PTYPE(vp) == P_INVDISPOW || PTYPE(vp) == P_MODN)
+	    do i = px1, px2 {
+		if (DISCUTOFF(vp) == NO) {
+		    frontvox = nvox
+		    backvox = 1
+		} else {
+		    frontvox = 1
+		    backvox = nvox
+		    do vox = 1, nvox {
+			vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],intelem]
+			if (vimin <= vox_int && vox_int < VIMAX(vp)) {
+			    frontvox = max (frontvox, vox)
+			    backvox = min (backvox, vox)
+			}
+		    }
+		}
+		if (frontvox - backvox < 0)
+		    next
+		do vox = backvox, frontvox {
+		    distwt = (real(vox-backvox+1) /
+			real(frontvox-backvox+1)) ** DISPOWER(vp)
+
+		    # Opacity transformation function.
+		    if (use_both) {
+			vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
+			    opelem] * OSCALE(vp)
+			if (vox_op < OMIN(vp))
+			    attenuate = amax
+			else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
+			    attenuate = amax - (vox_op - OMIN(vp)) * ofac
+			else
+			    attenuate = amin
+			oline[i] = oline[i] * attenuate
+		    }
+
+		    # Intensity transformation function.
+		    vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
+			intelem]
+		    if (vox_int < vimin)
+			ival = IIMIN(vp)
+		    else if (vimin <= vox_int && vox_int < VIMAX(vp))
+			ival = IIMIN(vp) + (vox_int - vimin) * ifac
+		    else
+			ival = IIMAX(vp)
+
+		    if (PTYPE(vp) == P_INVDISPOW)
+			oline[i] = oline[i] + ival * distwt
+		    else if (mod (int (ival/(IIMAX(vp)-IIMIN(vp)) * 100.0),
+			    MODN(vp)) == 0)
+			    oline[i] = oline[i] + ival * distwt
+		}
+	    }
+	else
+	    do i = px1, px2
+		do vox = 1, nvox {
+		    # Opacity transformation function.
+		    if (PTYPE(vp) == P_ATTENUATE || use_both) {
+			vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
+			    opelem] * OSCALE(vp)
+			if (vox_op < OMIN(vp))
+			    attenuate = amax
+			else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
+			    attenuate = amax - (vox_op - OMIN(vp)) * ofac
+			else
+			    attenuate = amin
+			oline[i] = oline[i] * attenuate
+		    }
+
+		    # Intensity transformation function.
+		    if (PTYPE(vp) != P_ATTENUATE || use_both) {
+			vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
+			    intelem]
+			if (vox_int < vimin)
+			    ival = IIMIN(vp)
+			else if (vimin <= vox_int && vox_int < VIMAX(vp))
+			    ival = IIMIN(vp) + (vox_int - vimin) * ifac
+			else
+			    ival = IIMAX(vp)
+
+			if (PTYPE(vp) == P_AVERAGE)
+			    oline[i] = oline[i] + ival * 1.0 / real(nvox)
+			else if (PTYPE(vp) == P_SUM)
+			    oline[i] = oline[i] + ival
+			else if (PTYPE(vp) == P_LASTONLY)
+			    if (ival > 0.0)
+				oline[i] = ival
+		    }
+		}
+end
+
+
diff --git a/pkg/proto/vol/src/x_vol.x b/pkg/proto/vol/src/x_vol.x
new file mode 100644
index 00000000..0d44bccb
--- /dev/null
+++ b/pkg/proto/vol/src/x_vol.x
@@ -0,0 +1,6 @@
+# X_VOL -- Volume projection and related tasks.
+
+task	pvol = t_pvol,
+	i2sun = t_i2sun,
+	im3dtran = t_im3dtran,
+	imjoin = t_imjoin