[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
`computeContour3d' by adding up the areas of the triangles making up the
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
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
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).
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
am I making a stupid error (if so which one...)?
if not so:
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?
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