Functions for plotting 3D volumetric data.
Description
slice3D
plots a 3D dataset with a color variable as slices or on surfaces.
slicecont3D
plots a 3D dataset with a color variable as contours on slices.
isosurf3D
plots isosurfaces from a 3D dataset.
voxel3D
plots isosurfaces as scatterpoints.
createisosurf
create the isosurfaces (triangulations)
from volumetric data. Its output can be plotted with triangle3D
.
createvoxel
creates voxels (x, y, z) points from volumetric data.
Its output can be plotted with scatter3D.
Usage
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42  slice3D (x, y, z, colvar, ..., phi = 40, theta = 40,
xs = min(x), ys = max(y), zs = min(z),
col = NULL, NAcol = "white", breaks = NULL,
border = NA, facets = TRUE, colkey = NULL,
panel.first = NULL, clim = NULL,
clab = NULL, bty = "b",
lighting = FALSE, shade = NA, ltheta = 135, lphi = 0,
add = FALSE, plot = TRUE)
slicecont3D (x, y, z, colvar, ..., phi = 40, theta = 40,
xs = NULL, ys = NULL, zs = NULL, level = NULL,
col = NULL, NAcol = "white", breaks = NULL,
border = NA, facets = TRUE,
colkey = NULL, panel.first = NULL,
clim = NULL, clab = NULL, bty = "b",
dDepth = 0, add = FALSE, plot = TRUE)
isosurf3D (x, y, z, colvar, ..., phi = 40, theta = 40,
level = mean(colvar, na.rm = TRUE), isofunc = createisosurf,
col = NULL, border = NA, facets = TRUE,
colkey = NULL, panel.first = NULL,
clab = NULL, bty = "b",
lighting = FALSE, shade = 0.5, ltheta = 135, lphi = 0,
add = FALSE, plot = TRUE)
voxel3D (x, y, z, colvar, ..., phi = 40, theta = 40,
level = mean(colvar, na.rm = TRUE), eps = 0.01, operator = "=",
col = NULL, NAcol = "white", breaks = NULL, colkey = FALSE,
panel.first = NULL, bty = "b", add = FALSE, plot = TRUE)
triangle3D (tri, colvar = NULL, ..., phi = 40, theta = 40,
col = NULL, NAcol = "white", breaks = NULL,
border = NA, facets = TRUE,
colkey = NULL, panel.first = NULL,
lighting = FALSE, shade = 0.5, ltheta = 135, lphi = 0,
clim = NULL, clab = NULL,
bty = "b", add = FALSE, plot = TRUE)
createisosurf (x, y, z, colvar, level = mean(colvar, na.rm = TRUE))
createvoxel (x, y, z, colvar, level = mean(colvar, na.rm = TRUE), eps = 0.01,
operator = "=")

Arguments
x, y, z 
Vectors with x, y and zvalues.
They should be of length equal to the first, second and
third dimension of 
colvar 
The variable used for coloring.
It should be an array of dimension equal to

tri 
A threecolumned matrix (x, y, z) with triangle coordinates. A triangle is defined by three consecutive rows. 
isofunc 
A function defined as 
theta, phi 
the angles defining the viewing direction.

col 
Colors to be used for coloring the 
NAcol 
Colors to be used for 
breaks 
a set of finite numeric breakpoints for the colors; must have one more breakpoint than color and be in increasing order. Unsorted vectors will be sorted, with a warning. 
border 
The color of the lines drawn around the surface facets.
The default, 
facets 
If 
colkey 
A logical, The default is to draw the color key on side = 4, i.e. in the right margin.
If 
panel.first 
A 
clab 
Only if 
clim 
Only if 
xs, ys, zs 
Vectors or matrices.
Vectors specify the positions in x, y or z where the slices (planes) are to be drawn.
The values of 
level 
The level(s) at which the There can be more than one level, but for 
lighting 
If not Will overrule See examples in jet.col. 
shade 
the degree of shading of the surface facets. Values of shade close to one yield shading similar to a point light source model and values close to zero produce no shading. Values in the range 0.5 to 0.75 provide an approximation to daylight illumination. See persp. 
ltheta, lphi 
if finite values are specified for 
bty 
The type of the box, the default only draws background panels.
Only effective if the persp
argument ( 
eps 
The voxel precision, only used when 
operator 
One of '=', '<', '>', '<>' for selection of points ‘equal’ (within precision), larger or smaller than the required level or to be within an interval. 
dDepth 
When a contour is added on an image, the image polygons may
hide some contour segments. To avoid that, the viewing depth of the segments
can be artificially decreased with the factor 
add 
Logical. If 
plot 
Logical. If 
... 
additional arguments passed to the plotting methods. The following persp arguments can be specified:
In addition, the perspbox arguments
For all functions, the arguments The arguments after ... must be matched exactly. 
Value
The plotting functions return the viewing transformation matrix, See trans3D.
Function createisosurf
returns a threecolumned matrix (x, y, z) with
triangle coordinates. One triangle is defined by three consecutive rows.
It can be plotted with triangle3D
.
Function createvoxel
returns a list with the elements x, y, z
defining the points that are at a distance of less than
eps*diff(range(colvar))
from the required level
.
Its output can be plotted with scatter3D.
Note
The isosurf3D
function uses function computeContour3d
,
from package misc3d
, which is based on the marching cubes algorithm.
Please cite the package misc3d
(Feng & Tierney, 2008) when using isosurf3D
.
For voxel3D
, coloring is always according to the zvariable. A more flexible
coloration can be achieved by using createvoxel
, followed by scatter3D.
See examples.
Author(s)
Karline Soetaert <karline.soetaert@nioz.nl>
References
Lorensen, W.E. and Cline, H.E., Marching Cubes: a high resolution 3D surface reconstruction algorithm, Computer Graphics, Vol. 21, No. 4, pp 163169 (Proc. of SIGGRAPH), 1987.
Dai Feng, Luke Tierney, Computing and Displaying Isosurfaces in R, Journal of Statistical Software 28(1), 2008. URL http://www.jstatsoft.org/v28/i01/.
See Also
Oxsat for another example of slice3D
.
plotdev for zooming, rescaling, rotating a plot.
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173  # save plotting parameters
pm < par("mfrow")
pmar < par("mar")
## =======================================================================
## Simple slice3D examples
## =======================================================================
par(mfrow = c(2, 2))
x < y < z < seq(1, 1, by = 0.1)
grid < mesh(x, y, z)
colvar < with(grid, x*exp(x^2  y^2  z^2))
# default is just the panels
slice3D (x, y, z, colvar = colvar, theta = 60)
# contour slices
slicecont3D (x, y, z, ys = seq(1, 1, by = 0.5), colvar = colvar,
theta = 60, border = "black")
slice3D (x, y, z, xs = c(1, 0.5, 0.5), ys = c(1, 0, 1),
zs = c(1, 0), colvar = colvar,
theta = 60, phi = 40)
## =======================================================================
## coloring on a surface
## =======================================================================
XY < mesh(x, y)
ZZ < XY$x*XY$y
slice3D (x, y, z, xs = XY$x, ys = XY$y, zs = ZZ, colvar = colvar,
lighting = TRUE, lphi = 90, ltheta = 0)
## =======================================================================
## Specifying transparent colors
## =======================================================================
par(mfrow = c(1, 1))
x < y < z < seq(4, 4, by = 0.2)
M < mesh(x, y, z)
R < with (M, sqrt(x^2 + y^2 + z^2))
p < sin(2*R) /(R+1e3)
## Not run:
# This is very slow  alpha = 0.5 makes it transparent
slice3D(x, y, z, colvar = p, col = jet.col(alpha = 0.5),
xs = 0, ys = c(4, 0, 4), zs = NULL, d = 2)
## End(Not run)
slice3D(x, y, z, colvar = p, d = 2, theta = 60, border = "black",
xs = c(4, 0), ys = c(4, 0, 4), zs = c(4, 0))
## =======================================================================
## A section along a transect
## =======================================================================
data(Oxsat)
Ox < Oxsat$val[, Oxsat$lat >  5 & Oxsat$lat < 5, ]
slice3D(x = Oxsat$lon, z = Oxsat$depth, y = 1:5, colvar = Ox,
ys = 1:5, zs = NULL, NAcol = "black",
expand = 0.4, theta = 45, phi = 45)
## =======================================================================
## isosurf3D example  rather slow
## =======================================================================
par(mfrow = c(2, 2), mar = c(2, 2, 2, 2))
x < y < z < seq(2, 2, length.out = 15)
xyz < mesh(x, y, z)
F < with(xyz, log(x^2 + y^2 + z^2 +
10*(x^2 + y^2) * (y^2 + z^2) ^2))
# use shading for level = 1  show triangulation with border
isosurf3D(x, y, z, F, level = 1, shade = 0.9,
col = "yellow", border = "orange")
# lighting for level  2
isosurf3D(x, y, z, F, level = 2, lighting = TRUE,
lphi = 0, ltheta = 0, col = "blue", shade = NA)
# three levels, transparency added
isosurf3D(x, y, z, F, level = seq(0, 4, by = 2),
col = c("red", "blue", "yellow"),
clab = "F", alpha = 0.2, theta = 0, lighting = TRUE)
# transparency can also be added afterwards with plotdev()
## Not run:
isosurf3D(x, y, z, F, level = seq(0, 4, by = 2),
col = c("red", "blue", "yellow"),
shade = NA, plot = FALSE, clab = "F")
plotdev(lighting = TRUE, alpha = 0.2, theta = 0)
## End(Not run)
# use of creatisosurf
iso < createisosurf(x, y, z, F, level = 2)
head(iso)
triangle3D(iso, col = "green", shade = 0.3)
## Not run:
# higher resolution
x < y < z < seq(2, 2, length.out = 50)
xyz < mesh(x, y, z)
F < with(xyz, log(x^2 + y^2 + z^2 +
10*(x^2 + y^2) * (y^2 + z^2) ^2))
# three levels
isosurf3D(x, y, z, F, level = seq(0, 4, by = 2),
col = c("red", "blue", "yellow"),
shade = NA, plot = FALSE, clab = "F")
plotdev(lighting = TRUE, alpha = 0.2, theta = 0)
## End(Not run)
## =======================================================================
## voxel3D example
## =======================================================================
par(mfrow = c(2, 2), mar = c(2, 2, 2, 2))
# fast but needs high resolution grid
x < y < z < seq(2, 2, length.out = 70)
xyz < mesh(x, y, z)
F < with(xyz, log(x^2 + y^2 + z^2 +
10*(x^2 + y^2) * (y^2 + z^2) ^2))
voxel3D(x, y, z, F, level = 4, pch = ".", cex = 5)
## =======================================================================
## rotation
## =======================================================================
plotdev(theta = 45, phi = 0)
plotdev(theta = 90, phi = 10)
# same using createvoxel  more flexible for coloring
vox < createvoxel(x, y, z, F, level = 4)
scatter3D(vox$x, vox$y, vox$z, colvar = vox$y,
bty = "g", colkey = FALSE)
## =======================================================================
## voxel3D to show hypox sites
## =======================================================================
par(mfrow = c(1, 1), mar = c(2, 2, 2, 2))
Hypox < createvoxel(Oxsat$lon, Oxsat$lat, Oxsat$depth[1:19],
Oxsat$val[,,1:19], level = 40, operator = "<")
panel < function(pmat) { # an image at the bottom
Nx < length(Oxsat$lon)
Ny < length(Oxsat$lat)
M < mesh(Oxsat$lon, Oxsat$lat)
xy < trans3D(pmat = pmat, x = as.vector(M$x), y = as.vector(M$y),
z = rep(1000, length.out = Nx*Ny))
x < matrix(nrow = Nx, ncol = Ny, data = xy$x)
y < matrix(nrow = Nx, ncol = Ny, data = xy$y)
Bat < Oxsat$val[,,1]; Bat[!is.na(Bat)] < 1
image2D(x = x, y = y, z = Bat, NAcol = "black", col = "grey",
add = TRUE, colkey = FALSE)
}
scatter3D(Hypox$x, Hypox$y, Hypox$z, colvar = Hypox$cv,
panel.first = panel, pch = ".", bty = "b",
theta = 30, phi = 20, ticktype = "detailed",
zlim = c(1000,0), xlim = range(Oxsat$lon),
ylim = range(Oxsat$lat) )
# restore plotting parameters
par(mfrow = pm)
par(mar = pmar)
