`level' definition in `computeContour3d' (misc3d package)

On Sat, 09 Nov 2013 17:16:28 +0100, Duncan Murdoch
<murdoch.duncan at gmail.com> wrote:

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
clarifying
whether 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 precise
meaning of `level' here. looking into the package, ultimately the level
is
used 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 argue
that always testing for `v >= level' would be better. if I feed data
with
discrete values (e.g. integer-valued) defined
on a coarse grid into `computeContour3d' it presently makes a big
difference whether there is a single data point (e.g.) with a value
larger
than `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 = 1

this 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
get
the 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.