Skip to content

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

6 messages · j. van den hoff, Duncan Murdoch

j. van den hoff
# 
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.

j.



--
Duncan Murdoch
#
On 13-11-09 8:50 AM, j. van den hoff wrote:
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?

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

            

      
      
      
    

        
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
    

    

      
      #
    

  


  

      
On 13-11-09 11:57 AM, j. van den hoff wrote:

            

      
      
      
    

        
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

  
  


  


      

    

    

      
      

        


  

        

      j. van den hoff
    

    

      
      #
    

  


  

      
On Sat, 09 Nov 2013 18:18:23 +0100, Duncan Murdoch

        
<murdoch.duncan at gmail.com> wrote:

        

            

      
      
      
    

        
I was imprecise: what I meant is: the isosurface should not change in my  
example between both cases.

            

      
      
      
    

        
understood/agreed.

            

      
      
      
    

        
yes, it does that. and it is clear that due to your interpretation it  
selects
about point 3.5 and 6.5 (?) for

data2 <- c(0, 0, 1, 1, 2, 1, 1, 0, 0).

still, depending on application I would maintain that it can make (more)  
sense
to keep the isosurface at 3,7 in this case.

I believe the problem maybe is not so much "discrete vs. continuous" but  
whether there is
a constant "plateau" in the data: even for the underlying continuous field  
it is a matter
of convention, then, where to put the level=1 contour: at the first  
crossing, in the middle
of the plateau, or at the second crossing. I understand `computeContour3d'  
essentially puts
the contour in the middle of the plateau.

I do not want to claim present behavior is a bug. it just is not  
necessarily what
is needed. an additional argument/flag to `computeContour3d', e.g., to  
select behaviour in this
sort of "degenerate" cases (exhibiting strictly constant plateaus) would  
be great.

 from a purely practical point of view: as explained I want to get the  
"outer isosurface" of such
preprocessed discretized data and quantify the surface area. the question  
here is "when do the data
first cross the threshold" and that's where the present behaviour causes a  
problem for me.

but anyway thanks for bothering. if it is deemed undesirable to change  
present behaviour (or to add
a further flag for controling the behaviour of `level') I can fix it  
locally for my needs by
just changing the test in `faceType' (or so I presume).

joerg

            

      
      
      
    

        
--

  
  


  


      

    

    

      
      

        


  

        

      Duncan Murdoch
    

    

      
      #
    

  


  

      
On 13-11-09 12:53 PM, j. van den hoff wrote:

            

      
      
      
    

        
I don't think it does.  I think it picks 4 and 6.  In your 3d example, 
the smaller cube runs from 9 to 13 in each coordinate (though it misses 
the corners).  You can see this if you plot it using the "rgl" engine, 
then call rgl::decorate3d() to add axes.

            

      
      
      
    

        
That makes just as much sense, but not more.  Anywhere from 3 to 4 is a 
sensible left end, anywhere from 6 to 7 is fine for the right end.

            

      
      
      
    

        
No, it doesn't.  As you've seen, it handles plateaus inconsistently 
depending on whether they are the max of the field or not.  For the case 
where it is an interior plateau, you get your contour at height 
level+epsilon in an approximation to the field.  For the max of the 
field, you get level.

Try comparing the contour you get at level 1 with v and at level -1 with 
-v.  They are not the same.

            

      
      
      
    

        
You can suggest this to the maintainer of the package (I am not author 
or maintainer of it).

            

      
      
      
    

        
You may need to write your own function for this.

Duncan Murdoch