# [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.

problem:
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
defaults).

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

library(misc3d)
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)
drawScene(makeTriangles(con))

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.

indeed,

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
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
(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:
Value:

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

--

```