1 files changed, 92 insertions, 0 deletions
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