Functions for plotting 3-D volumetric data.

Share:

Description

slice3D plots a 3-D dataset with a color variable as slices or on surfaces.

slicecont3D plots a 3-D dataset with a color variable as contours on slices.

isosurf3D plots isosurfaces from a 3-D 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 z-values. They should be of length equal to the first, second and third dimension of colvar respectively.

colvar

The variable used for coloring. It should be an array of dimension equal to c(length(x), length(y), length(z)). For triangle3D, colvar should be of length = nrow(tri) / 3. It must be present.

tri

A three-columned matrix (x, y, z) with triangle coordinates. A triangle is defined by three consecutive rows.

isofunc

A function defined as function(x, y, z, colvar, level), and that returns the three-columned matrix with triangle coordinates. The default, createisosurf uses function computeContour3d from package misc3d.

theta, phi

the angles defining the viewing direction. theta gives the azimuthal direction and phi the colatitude. see persp.

col

Colors to be used for coloring the colvar variable. If col is NULL then a red-yellow-blue colorscheme (jet.col) will be used.

NAcol

Colors to be used for colvar values that are NA.

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, NA, will disable the drawing of borders.

facets

If TRUE, then col denotes the color of the surface facets. If FALSE, then the surface facets are colored “white” and the border (if NA) will be colored as specified by col. If NA then the facets will be transparent. It is usually faster to draw with facets = FALSE.

colkey

A logical, NULL (default), or a list with parameters for the color key (legend). List parameters should be one of side, plot, length, width, dist, shift, addlines, col.clab, cex.clab, side.clab, line.clab, adj.clab, font.clab and the axis parameters at, labels, tick, line, pos, outer, font, lty, lwd, lwd.ticks, col.box, col.axis, col.ticks, hadj, padj, cex.axis, mgp, tck, tcl, las. The defaults for the parameters are side = 4, plot = TRUE, length = 1, width = 1, dist = 0, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL) See colkey.

The default is to draw the color key on side = 4, i.e. in the right margin. If colkey = NULL then a color key will be added only if col is a vector. Setting colkey = list(plot = FALSE) will create room for the color key without drawing it. if colkey = FALSE, no color key legend will be added.

panel.first

A function to be evaluated after the plot axes are set up but before any plotting takes place. This can be useful for drawing background grids or scatterplot smooths. The function should have as argument the transformation matrix, e.g. it should be defined as function(pmat). See last example and example of persp3D.

clab

Only if colkey is not NULL or FALSE, the label to be written on top of the color key. The label will be written at the same level as the main title. To lower it, clab can be made a vector, with the first values empty strings.

clim

Only if colvar is specified, the range of the color variable, used for the color key. Values of colvar that extend the range will be put to NA.

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 colvar will be projected on these slices. Matrices specify a surface on which the colvar will be projected.

level

The level(s) at which the contour will be generated or the isosurfaces generated.

There can be more than one level, but for slicecont3D too many will give a crowded view, and one is often best. For isosurf3D, the use of multiple values may need transparent colors to visualise. For voxel3D, level should either be one number (if operator equals '=', '<', '>') or two numbers (for operator = '<>').

lighting

If not FALSE the facets will be illuminated, and colors may appear more bright. To switch on lighting, the argument lighting should be either set to TRUE (using default settings) or it can be a list with specifications of one of the following: ambient, diffuse, specular, exponent, sr and alpha.

Will overrule shade not equal to NA.

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 ltheta and lphi, the surface is shaded as though it was being illuminated from the direction specified by azimuth ltheta and colatitude lphi. See persp.

bty

The type of the box, the default only draws background panels. Only effective if the persp argument (box) equals TRUE (this is the default). See perspbox.

eps

The voxel precision, only used when operator = "=". A point is selected if it closer than eps*diff(range(colvar)) to the required level.

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 dDepth times the persp argument expand (usually = 1), to make them appear in front of the polygons. Too large values of dDepth may create visible artifacts. See contour3D.

add

Logical. If TRUE, then the slices, voxels or surfaces will be added to the current plot. If FALSE a new plot is started.

plot

Logical. If TRUE (default), a plot is created, otherwise the viewing transformation matrix is returned (as invisible).

...

additional arguments passed to the plotting methods.

The following persp arguments can be specified: xlim, ylim, zlim, xlab, ylab, zlab, main, sub, r, d, scale, expand, box, axes, nticks, ticktype. The arguments xlim, ylim, zlim only affect the axes. All objects will be plotted, including those that fall out of these ranges. To select objects only within the axis limits, use plotdev.

In addition, the perspbox arguments col.axis, col.panel, lwd.panel, col.grid, lwd.grid can also be given a value.

alpha can be given a value inbetween 0 and 1 to make colors transparent.

For all functions, the arguments lty, lwd can be specified.

The arguments after ... must be matched exactly.

Value

The plotting functions return the viewing transformation matrix, See trans3D.

Function createisosurf returns a three-columned 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 z-variable. 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 163-169 (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+1e-3)

## 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)