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

[j. van den hoff]

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.

is this clearer now?

>
> Duncan Murdoch
>
>>
>> j.
>>
>>
>>
>> --
>>
>> ______________________________________________
>> R-help at 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.
>>
>


--