aboutsummaryrefslogtreecommitdiff
path: root/pkg/proto/vol/src/t_pvol.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/proto/vol/src/t_pvol.x')
-rw-r--r--pkg/proto/vol/src/t_pvol.x284
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
generated by cgit (git 2.34.1) at 2025-09-20 07:11:53 -0400