diff options
Diffstat (limited to 'pkg/proto/vol')
62 files changed, 6537 insertions, 0 deletions
diff --git a/pkg/proto/vol/README b/pkg/proto/vol/README new file mode 100644 index 00000000..b92c2da6 --- /dev/null +++ b/pkg/proto/vol/README @@ -0,0 +1,26 @@ +The VOLumes package is a possibly temporary collection of tasks related to +manipulating and viewing 3d or in some cases 4d "volume" images, and a few +other things. + +IMJOIN joins sets of N-dimensional images together along a specified axis. +IM3DTRAN performs 3d image transposes; if appropriate [*,-*,*] type image +sections are given as input, it also accomplishes rotates. Tasks such as +these may later be integrated into a standard IRAF package. + +PVOL projects through volume images, casting rays onto a set of output +2d images distributed along a great circle around the volume image. When +the output images are displayed or recorded onto video and played back, +the volume image appears to rotate. Various translucency and opacity +algorithms are employed. + +I2SUN is a temporary task for converting IRAF images into Sun rasterfiles, +primarily to take advantage of a Sun-specific MOVIE utility for viewing +digital movies on a workstation screen; it will no longer be necessary when +the IRAF image display servers can display movies. + +# Not currently distributed: +#VIDRECORD is an NOAO-specific IRAF script that is used to display successive +#images from template or datacube on a display device connected to a recorder. +#The recorder is assumed to be connected via rs-232 to a port on some IRAF node. +#Other sites may wish to modify the associated host-level task that sends the +#device-control commands to the rs-232 port. diff --git a/pkg/proto/vol/README.install b/pkg/proto/vol/README.install new file mode 100644 index 00000000..c3f23147 --- /dev/null +++ b/pkg/proto/vol/README.install @@ -0,0 +1,107 @@ + + Installation Instructions for the VOL Addon Package + +The volume-images package, VOL, is being distributed separately +from the IRAF V2.8 distribution. This package is being made available for +user testing on a user-beware basis. The installation instructions +that follow assume that you have copied the tar format VOL archive onto +your host machine. The method you use to copy the file (or remotely access +the tar file) is OS dependent and is not discussed in this document. If your +IRAF system has been stripped of the IRAF libraries ('mkpkg stripall') you +will not be able to build the VOL executable as described here. You must +either reload the required libraries or request a binary distribution +of VOL for your operating system. If you have any questions, please send +electronic mail to the IRAF project at + +Internet (ARPAnet, MILnet, NSFnet, etc.) iraf@noao.edu +BITnet iraf@noao.edu * +SPAN/HEPnet (DECnet) noao::iraf or 5355::iraf +UUCP/Usenet {arizona,decvax,ncar}!noao!iraf + uunet!noao.edu!iraf +* = through a user-specified gateway + + [IRAF Hotline: (602) 323-4160] + +For discussion of algorithms or ongoing bug fixes etc., contact the author, +Steve Rooke, at 602-325-9399 (or rooke@noao.edu or 5355::rooke). + +This package is being distributed as an external package making use of the +layered software enhancements in IRAF V2.8. + +[1] This discussion assumes you intend to install VOL as an independent + external layered package. You may choose to merge it into your custom + "local" package instead, but would need to edit several files not + discussed herein. + +[2] Log into the CL from the IRAF installation account. This insures you + have the proper permissions and the files will have the proper owner. + + Create a directory for vol, preferably external to the IRAF directory + tree to simplify future IRAF updates, then set an IRAF environment + variable to point to it. + + UNIX example: + cl> reset vol /local/iraftest/vol/ + VMS example: + cl> reset vol usr\$1:[localiraf.vol] + +[3] Change directory to vol and unpack the tar archive. You must + load the softools package before executing rtar. Note that the <ARCHIVE> + name must be given with the appropriate path name at the host OS level, + or may be specified as a tape device if that is how you received the + package. + + cl> cd vol + cl> softools + so> rtar -xrf <ARCHIVE> # Read the archive + so> cd vol # Change to VOL subdirectory + +[4] When the archive has been unpacked, build the VOL package executable: + + 1) Delete any bin directory symbolic links you find + 2) Create a directory called bin.ARCH where ARCH is something + like "f68881", "ffpa", "sparc", or "vax"; see your IRAF + installation guide. + 3) Then issue the following commands: + + so> mkpkg ARCH + so> mkpkg -p vol update >& spool & + + The first command sets the bin directory to be properly configured + and the second recompiles the package. For Sun multiple architecture + support, please refer to the "Sun/IRAF Site Manager's Guide". + +[5] The spool file(s) should be reviewed upon completion to make sure + there were no errors. When you are confident the installation was + successful, delete the spool file (and the archive file if you wish). + +[6] Still logged in as iraf, edit the one file in the core iraf system + that is necessary to install an external layered package: + + so> cd hlib + so> edit extern.pkg + + UNIX example, leaving out other external packages: + reset vol = /local/iraftest/vol/ + task vol.pkg = vol$vol.cl + + reset helpdb = "lib$helpdb.mip\ + [ ... other packages ... ] + ,vol$lib/helpdb.mip\ + " + VMS example, leaving out other external packages: + reset vol = usr\$2:[localiraf.vol] + task vol.pkg = vol$vol.cl + + reset helpdb = "lib$helpdb.mip\ + ,vol$lib/helpdb.mip\ + " + +[7] Finally, build the help database. Still logged in as iraf: + + so> mkhelpdb vol$lib/root.hd vol$lib/helpdb.mip + +That's it ... you should be ready to go. If you have any questions or +problems, please do not hesitate to send email to iraf@noao.edu or call the +IRAF HOTLINE at 602-323-4160 or the author at 602-325-9399. + diff --git a/pkg/proto/vol/Revisions b/pkg/proto/vol/Revisions new file mode 100644 index 00000000..c3e53369 --- /dev/null +++ b/pkg/proto/vol/Revisions @@ -0,0 +1,12 @@ +.help revisions Sep89 vol +.nf +vol$src/vtransmit.gx + Was attempting to calculate opacity factor even when mode did not + include opacity (task default parameters include INDEFs). Generated + a floating exception in iraf 2.8 not present before. (7/5/89 SRo) + +==== +V2.8 +==== + +.endhelp diff --git a/pkg/proto/vol/lib/helpdb.mip b/pkg/proto/vol/lib/helpdb.mip Binary files differnew file mode 100644 index 00000000..18d59be5 --- /dev/null +++ b/pkg/proto/vol/lib/helpdb.mip diff --git a/pkg/proto/vol/lib/mkpkg.inc b/pkg/proto/vol/lib/mkpkg.inc new file mode 100644 index 00000000..de1341ab --- /dev/null +++ b/pkg/proto/vol/lib/mkpkg.inc @@ -0,0 +1,7 @@ +# Global MKPKG definitions for the VOL package. + +$set XFLAGS = "$(XFLAGS) -p vol" + +# Special file lists, not needed at present + +#$include "vol$lib/mkpkg.sf.SUN3" diff --git a/pkg/proto/vol/lib/root.hd b/pkg/proto/vol/lib/root.hd new file mode 100644 index 00000000..db09387f --- /dev/null +++ b/pkg/proto/vol/lib/root.hd @@ -0,0 +1,5 @@ +# Root help directory for the VOL packages. This dummy package is necessary +# in order to have `vol' appear as a module in some package, so that the user +# can type "help vol" (with the `vol' given as a task). + +_vol pkg = vol$lib/rootvol.hd diff --git a/pkg/proto/vol/lib/rootvol.hd b/pkg/proto/vol/lib/rootvol.hd new file mode 100644 index 00000000..ceba8059 --- /dev/null +++ b/pkg/proto/vol/lib/rootvol.hd @@ -0,0 +1,8 @@ +# Root task entry for the VOL package help tree. Defines `vol' +# as both a task and a package in the vol help database. + +vol men = vol$vol.men, + hlp = vol$vol.men, + sys = vol$vol.hlp, + pkg = vol$vol.hd, + src = vol$vol.cl diff --git a/pkg/proto/vol/lib/strip.vol b/pkg/proto/vol/lib/strip.vol new file mode 100644 index 00000000..12157035 --- /dev/null +++ b/pkg/proto/vol/lib/strip.vol @@ -0,0 +1,12 @@ +# STRIP.VOL -- Rmfiles command script, used to strip the VOL directory +# of all files not required for ordinary runtime use of the system. + +src -allbut .hlp .hd .men .cl .par .key .dat .mip + +# Sun/IRAF only. +# --------------- +-file bin.f68881/OBJS.arc +-file bin.ffpa/OBJS.arc +-file bin.sparc/OBJS.arc +-file bin.generic/OBJS.arc +-file bin.pg/OBJS.arc diff --git a/pkg/proto/vol/lib/zzsetenv.def b/pkg/proto/vol/lib/zzsetenv.def new file mode 100644 index 00000000..637e5015 --- /dev/null +++ b/pkg/proto/vol/lib/zzsetenv.def @@ -0,0 +1,7 @@ +# Global environment definitions for the VOL package. + +set vollib = "vol$lib/" +set volsrc = "vol$src/" +set volbin = "vol$bin(arch)/" + +keep diff --git a/pkg/proto/vol/mkpkg b/pkg/proto/vol/mkpkg new file mode 100644 index 00000000..bcd5e283 --- /dev/null +++ b/pkg/proto/vol/mkpkg @@ -0,0 +1,21 @@ +# MKPKG file for the VOL Package + +$call update +$exit + +update: + $call update@src + + $ifeq (HOSTID, vms) $purge [...] $endif + ; + +relink: + $call relink@src + + $ifeq (HOSTID, vms) $purge [...] $endif + ; + +install: + $call install@src + ; + 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 diff --git a/pkg/proto/vol/vol.cl b/pkg/proto/vol/vol.cl new file mode 100644 index 00000000..d3a6a5b1 --- /dev/null +++ b/pkg/proto/vol/vol.cl @@ -0,0 +1,22 @@ +#{ VOL -- Volume Images Package + +print(" ") +print("This package contains tasks for viewing and manipulating 3d images.") +print("It is a pre-release version, and does not reflect the ultimate") +print("partitioning of n-dimensional image tasks within IRAF") +print(" ") + +# Load some needed packages now +# (None are needed yet) # if images$tv not loaded, load for vidrecord. + +package vol + +task i2sun, + im3dtran, + imjoin, + pvol = "vol$src/x_vol.e" + +#task vidrecord = "vol$src/vidrecord.cl" + +clbye() + diff --git a/pkg/proto/vol/vol.hd b/pkg/proto/vol/vol.hd new file mode 100644 index 00000000..2dcd0e39 --- /dev/null +++ b/pkg/proto/vol/vol.hd @@ -0,0 +1,10 @@ +# Help directory for the VOL package + +$doc = "vol$src/doc/" +$im3dtran = "vol$src/im3dtran/" +$i2sun = "vol$src/i2sun/" + +i2sun hlp=doc$i2sun.hlp, src=i2sun$t_i2sun.x +im3dtran hlp=doc$im3dtran.hlp, src=im3dtran$t_im3dtran.x +imjoin hlp=doc$imjoin.hlp, src=vol$src/t_imjoin.x +pvol hlp=doc$pvol.hlp, src=vol$src/t_pvol.x diff --git a/pkg/proto/vol/vol.men b/pkg/proto/vol/vol.men new file mode 100644 index 00000000..4184f44e --- /dev/null +++ b/pkg/proto/vol/vol.men @@ -0,0 +1,4 @@ + i2sun - Convert IRAF images to Sun rasterfiles + im3dtran - 3d image transpose (used for rotates as well) + imjoin - N-dimensional image join along arbitrary axis + pvol - Project volume image (generates 'rotating' volume images) diff --git a/pkg/proto/vol/vol.par b/pkg/proto/vol/vol.par new file mode 100644 index 00000000..da7b65e6 --- /dev/null +++ b/pkg/proto/vol/vol.par @@ -0,0 +1,3 @@ +# VOL package parameter file + +version,s,h,"May89" |