[R] `level' definition in `computeContour3d' (misc3d package)

[Duncan Murdoch]

On 13-11-09 11:57 AM, j. van den hoff wrote:
> 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.
>

I think it does, but your data make the determination of its location 
ambiguous.

The definition is the estimated location where the continuous field 
sampled at v crosses level.

You have a field with a discontinuity (or two).  You have whole volumes 
of space where the field is equal to the level.  The marching cubes 
algorithm is designed to detect crossings, not solid regions.

For example, going back to one dimension, if your data looked like your 
original vector

data1 <- c(0, 0, 1, 1, 1, 1, 1, 0, 0)

then it is ambiguous where it crosses 1:  it could be at 3 and 7, or 
there could be multiple crossings in that range.  I believe the 
analogous situation in misc3d would treat this as a crossing at 3 and 7.

Duncan Murdoch