[R] how to derive true surface area from `computeContour3d' (misc3d package)

j. van den hoff veedeehjay at googlemail.com
Fri Nov 8 12:49:14 CET 2013

I want to compute the total surface area of an isosurface in 3D space as  
approximated with
`computeContour3d' by adding up the areas of the triangles making up the
contour surface.

the vertex matrix returned by `computeContour3d' in general seems to  
provide the vertices
not in the frame of reference in which the original data are given but  
apparently rather
after some linear transformation (scaling + translation (+ rotation?) -- or
I am having some fundamental misconception of what is going on.

I'm interested in the simplest case where the input data are provided as a  
3D array
on an equidistant grid (i.e. leaving the x,y,z arguments at their  

e.g. (slight modification of `example(computeContour3d)'):

x <- seq(-1,1,len=11)
g <- expand.grid(x = x, y = x, z = x)
v <- array(g$x^4 + g$y^4 + g$z^4, rep(length(x),3))
con <- computeContour3d(v, max(v), 1)

this is (approximately) a cube with edge length 10 (taking the grid  
spacing as the unit of length).
so the expected (approximate) surface area is 600.


apply(con, 2, range) yields

      [,1] [,2] [,3]
[1,]    1    1    1
[2,]   11   11   11

which might be interpreted as providing the vertices in coordinates
where the grid spacing is used as unit of length. however
I get an area of only about 430 instead of approx. 600 which is already a  
much much larger deviation
 from the ideal cube surface than I would have expected given the small  
amount of smoothing at the
box edges and corners (but I have to double-check whether my triangle area  
computation is right, although I believe it is).

choosing instead

x <- seq(-2,2,len=50)

however, the corresponding range of `con' is

       [,1]   [,2]   [,3]
[1,] 13.274 13.274 13.274
[2,] 37.726 37.726 37.726

which cannot be the "grid coordinates" (which should be in the range  
[1,50]). adopting this interpretation nevertheless
(vertices are given in grid coordinates)
the sum of the triangle areas only amounts to about 2600 instead of the  
expected approx. 49^2*6 = 14406

question 1:
am I making a stupid error (if so which one...)?

if not so:

question 2:
is there a linear transformation from the original grid coordinates (with  
range from 1 to dim(v)[n], n=1:3)
involved which yields the reported vertex coordinates?

question 3:
could someone please explain where to find this information (even if  
hidden in the source code of the package)
how to convert the vertex coordinates as delivered by `computeContour3D'  
to 'grid coordinates' (or true world coordinates
in general (if the x,y,z arguments are specified, too)?

for the wishlist: it would of course be nice if `computeContour3d' would  
indeed return the total surface area itself,
e.g. as an attribute of the triangles matrix.

for the devs: there is a typo in the manpage of this function:

      A matrix of three columns representing the triangles making up the
      contour surface. Each row represents a vertex and goups of three
      rows represent a triangle.

(should be `groups' instead of `goups')


More information about the R-help mailing list