Repatch (from linux) of OSX IRAF

1 files changed, 284 insertions, 0 deletions

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