Marching Cubes, or well Tetrahedrons

General
1

Hi. Yesterday evening I looked at this site about marching cubes. http://paulbourke.net/geometry/polygonise/

But I decided to implement the tetrahedron version of it instead, since it was a bit simpler.

https://gist.github.com/tnlogy/9073640

Maybe you have some idea on how to optimize it, currently using if-statements. I tried to use a lookup table instead, to implement the switch-case, but it slowed it down rather than sped it up. And I guess I could be a bit smarter in the x-y-z loop calculating the values. :slight_smile:

Not showing much interesting graphics for now, but could be used for many fun things :slight_smile:

2

incredible link!
how do you adapt the number of facets?

3

it is the step argument to volume() that defines the granularity in my example, but it is quite flexible

4

@tnlogy:

Very nice. I have three ideas if you want to speed it up:

(1) Use buffer()s instead of tables.

(2) Replace the chained if-elsif… chain by a divide-and-conquer approach.

(3) Replace m = m+8 by m = 8-m

Regarding (3), note that you can get to cases m == 8..0xE only if you go through m = m+8. However, those cases have the same effect as their 16-m counterparts I believe (e.g., 0xE and 1 have the same effect, etc.). So if you replace that line m = m+8 by m = 8-m (i.e., 16-(m+8)), you would only have to check for m == 1..7 (double-check me on this; I haven’t worked my way carefully through it).

Combined with (2), I think that would lead to just three comparisons for any m value.

5

A few more thoughts.

Exit early on the m == 0 case (which I think is the most common, by far).

I didn’t pay attention to the interpolation function (ip) in my first reading, but it’s actually critical in the sense that it’s the most called. So it’s worth caching math.abs (globally for sure; locally probably also) as well as iso-v1 and p2-p1. (There may be other optimizations possible here if you examine the caller pattern; not sure.)

Every iteration of the x loop could reuse half (4 out of 8) of the samples determined in the previous iteration. For an expensive field computation (function fn) that could be a significant win.

It would be nice to compute normals for your triangles, and to orient them correctly. For the orientation, I think that in the m = 8-m case the triangles should be reversed (i.e., introduce a flag for that). For the normals, I’m guessing the best way is to compute them from the field values themselves (fns[] array). I’d have to look up how to do that. (I’m thinking a cross product of tangents, which would be derivatives approximated by finite differences. That may require additional sampling, though, which isn’t cheap.)

6

Darn, you’ve got me with that problem. I’m having a hard time focusing on work instead of this :stuck_out_tongue: I suspect I’ll end up writing a version of this myself at some point.

Here is another idea: Work by 3D “blocks” (say, 10x10x10 points in each block). So instead of 3 nested loops, you’ll have 6, but the outer three will stride 10 times as fast (and the inner ones just go to 10). Precompute the fn values in each block, but subtract the desired iso value up front: This will avoid most repeated fn evaluations, and allows you to compare and interpolate wrt. zero afterward. You can also use the precomputed fn(p)-iso values to compute field gradients (which are normals!) using finite differences (central differences are probably the best compromise). It’s probably also worth precomputing a buffer of p vectors, even though their trivial. In addition, for every cube whose precomputed values all have the same sign (which should be the majority for many interesting fn functions) you can skip the six poligonise calls altogether.

Now, I just have to get this problem out of my head until I have a chance to actually experiment with it… 8-}

7

haha, thanks for the comments! I will try them out. Would be nice to see an alternative solution, so go for it. :wink: or maybe build a marching cubes algorithm?

my current is quite naive, just wanted to try it out. enjoy reading the suggestions!

8

will maybe use it to make this old project better, which just used billboards instead of polygons.

http://www.youtube.com/watch?v=SOtriszSUac

9
  1. what do you mean by buffer? shader buffer instead vertices, hmm, need some more explanation I think. :slight_smile:
  2. what kind of divide and thinking about? I tried writing it like this instead, but it slowed it down, maybe due to the loops on the lookup values instead of inline?
local tbl = {}
tbl[0xe] = {{1,2,1,3,1,4}}
tbl[1] = tbl[0xe]
tbl[0xd] = {{2,1,2,4,2,3}}
tbl[2] = tbl[0xd]
tbl[0xc] = {{1,4,1,3,2,4}, {2,4,2,3,1,3}}
tbl[3] = tbl[0xc]
tbl[0xb] = {{3,1,3,2,3,4}}
tbl[4] = tbl[0xb]
tbl[0xa] = {{1,2,3,4,1,4}, {1,2,2,3,3,4}}
tbl[5] = tbl[0xa]
tbl[9] = {{1,2,2,4,3,4}, {1,2,1,3,3,4}}
tbl[6] = tbl[9]
tbl[7] = {{4,1,4,3,4,2}}
tbl[8] = tbl[7]

function polygonise(ps, fns, iso, i0,i1,i2,i3, vs)
    local pts = {ps[i0],ps[i1],ps[i2],ps[i3]}
    local vts = {fns[i0],fns[i1],fns[i2],fns[i3]}
    local m = 0
    if vts[1] < iso then m = m + 1 end
    if vts[2] < iso then m = m + 2 end
    if vts[3] < iso then m = m + 4 end
    if vts[4] < iso then m = m + 8 end

    local ts = tbl[m]
    if ts then
        for i,t in ipairs(ts) do
            local t0,t1,t2,t3,t4,t5 = unpack(t)
            table.insert(vs, ip(iso,pts[t0],pts[t1],vts[t0],vts[t1]))
            table.insert(vs, ip(iso,pts[t2],pts[t3],vts[t2],vts[t3]))
            table.insert(vs, ip(iso,pts[t4],pts[t5],vts[t4],vts[t5]))
        end
    end
end
10

I dont think (3) works, or at least not in my tests. I think you miss some cases then. See m as an encoding of the corners of the tetra to define if the point is within the volume or not.

11

@tnlogy:

what do you mean by buffer?

The mesh structure points to buffers, which are manually resized and efficiently indexed (a little better than tables). You can create your own like so:

    b = buffer(); b:resize(10000)
    b[1] = vec3(1, 2, 3)

(I’ve been working on some mesh generators, and I use these all the time.)

what kind of divide and thinking about?

“Divide and conquer”. Instead of writing:

  if m == 1 then
    ...
  elseif m == 2 then
    ...
  elseif m == 3 then
    ...
  elseif m == 4 then
    ...
  elseif m == 5 then
    ...
  elseif m == 6 then
    ...
  else if m == 7 then
    ...
  end

write:

  if m < 5 then
    if m < 3 then
      if m == 1 then
        ...
      else -- m == 2 case
        ...
      end
    else
      if m == 3 then
        ...
      else -- m == 4 case
        ...
      end
    end
  else
    if m == 5 then
      ...
    elseif m == 6 then
      ...
    else
      ...
    end
  end

With the latter you’ll perform less work on average.

I dont think (3) works, or at least not in my tests. I think you miss some cases then. See m as an encoding of the corners of the tetra to define if the point is within the volume or not.

I’m pretty sure there is a solution, but I’ll confirm after I code it up (maybe tonight, but probably not; professional bugs are getting priority :wink: ).

For your table approach, try making everything local. That way the Lua interpreter can keep everything in “registers” rather than having to look up the table every time. Also, use table constructor notation instead of inserting subtables one at a time (and out of order): I suspect the out-of-order insertion (e.g., tbl[0xe] before tbl[1]) might cause Lua to treat it via hashing instead of indexing (but I’m not sure what the Lua heuristic is for small indices).

12

Ah, thanks for explaining! I’ve read somewhere on the forum that buffers were slower, but good that they arent!

I stared myself blind since my if-else was based on a switch statement. Thanks!

13

@tnlogy:

I’ve read somewhere on the forum that buffers were slower, …

My mistake: I don’t know why I had concluded that buffers were faster than tables, but additional benchmarking indicates that that’s not the case.

14

Ah, maybe need a performance FAQ to keep in memory what is faster. :slight_smile:

It’s also good to know that it is faster to create a table of vertices, rather than modifying an existing array with mesh:vertex(index, …).