j. van den hoff
2013-Nov-08 11:49 UTC
[R] how to derive true surface area from `computeContour3d' (misc3d package)
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 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: 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') --