On 13-11-09 8:50 AM, j. van den hoff wrote:I'd very much appreciate some help here: I'm in the process of clarifyingwhether I can use `computeContour3d' to derive estimates of the surface area of a single closed isosurface (and prospectively the enclosed volume). getting the surface area from the list of triangles returned by `computeContour3d' is straightforward but I've stumbled over the precisemeaning of `level' here. looking into the package, ultimately the level isused in the namespace function `faceType' which reads:function (v, nx, ny, level, maxvol) { if (level == maxvol) p <- v >= level else p <- v > level v[p] <- 1 v[!p] <- 0 v[-nx, -ny] + 2 * v[-1, -ny] + 4 * v[-1, -1] + 8 * v[-nx, -1] } my question: is the discrimination of the special case `level == maxvol' (or rather of everything else) really desirable? I would arguethat always testing for `v >= level' would be better. if I feed data withdiscrete values (e.g. integer-valued) defined on a coarse grid into `computeContour3d' it presently makes a bigdifference whether there is a single data point (e.g.) with a value largerthan `level' or not. consider the 1D example: data1 <- c(0, 0, 1, 1, 1, 1, 1, 0, 0) data2 <- c(0, 0, 1, 2, 1, 1, 1, 0, 0) and level = 1this defines the isocontour `level = 1' to lie at pos 3 and 7 in for data1 but as lying at pos 4 in data2. actually I would like (and expect) to getthe same isosurface for `data2' with this `level' setting. in short: the meaning/definition of `level' changes depending on whether or not it is equal to `maxvol'. this is neither stated in the manpage nor is this desirable in my view. but maybe I miss something here. any clarification would be appreciated.I don't see why you'd expect the same output from those vectors, but since they aren't legal input to computeContour3d, maybe I don't know what you mean by them. Could you put together a reproducible example that shows bad contours?
it's not "bad" contours, actually. my question only concerns the different meaning
of `level' depending on whether `level = maxvol' or not.
here is a real example:
8<------------------------------------------------
library(misc3d)
dim <- 21
cnt <- (dim+1)/2
wid1 <- 5
wid2 <- 1
rng1 <- (cnt-wid1):(cnt+wid1)
rng2 <- (cnt-wid2):(cnt+wid2)
v <- array(0, rep (dim, 3))
#put 11x11x11 box of ones at center
v[rng1, rng1, rng1] <- 1
con1 <- computeContour3d(v, level = 1)
drawScene(makeTriangles(con1))
dum <- readline("CR for next plot")
#put an additional 3x3x3 box of twos at center
v[rng2, rng2, rng2] <- 2
con2 <- computeContour3d(v, level = 1)
drawScene(makeTriangles(con2))
8<------------------------------------------------
this first puts a 11x11x11 box one Ones at the center of the
zero-initalized array and computes `con1' for `level=1'. in the 2. step
it puts a further, 3x3x3 box of Twos at the center and computes the
`level=1' contour again which this time does not delineate
the box of Ones but lies somewhere between the two non-zero boxes since
now the test in `faceType' is for `> level'. this is not immediately
obvious from the plots (no scale) but obvious from looking at `con1' and
`con2': the `con2' isosurface is shrunk by 3 voxels at each
side relative to `con1' (so my initial mail was wrong here: `con2' does
not "jump" to the next "discrete" isocontour but rather to
a point about halfway between both plateaus ). I also (for my own problem
at hand) computed the total surface area which is
(not surprisingly...) 600 for `con1' and 64.87 for `con2'. so if one is
interested in such surfaces (I am) this makes a big difference in such
data.
the present behavior is not "wrong" per se but I would much prefer if the test where always for `>= level' (so that in the present example the resulting isosurface would in both cases delineate the box of Ones -- as is the case when using `level = 1-e-6' instead of `level=1').
I believe the isosurface for a given value of `level' should have an unambiguous meaning independent of what the data further "inside" are looking like.
is this clearer now?
Duncan Murdochj. -- ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
-- ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
