Hilariously Fast Volume Computation with the Divergence Theorem

(alyssarosenzweig.ca)

164 points | by luu 6 hours ago ago

36 comments

  • physicsguy 3 hours ago ago

    This is one of those when you go "Huh, this is amazing!" or "Huh, I thought this trick was really well known!" depending on your background ;)

    Here's a similar impl from 1980 written in Fortran that also computes other properties like centroid: https://calgo.acm.org/550.zip Algorithm 550: Solid Polyhedron Measures A. M. Messner and G. Q. Taylor ACM Trans. Math. Softw., 6(1), Mar 1980, pp.121--130 Keywords: polyhedron, graphics, numerical integration Language: Fortran 66/77; Shar Index: Z; Gams: P File size: 19.1 KB;

    But Messner published it first in: A. M. Messner, "A surface Integral method for computer calculation of mass properties", Paper No. 852, 29TH ANNUAL CONF. OF THE SOCIETY OF AERONAUTICAL WEIGHT ENGINEERS, Washington, D.C., May 1970.

    I think

  • eterevsky 5 hours ago ago

    Isn't the same as just taking every triangle from the mesh, calculating the volume of a prism-like polytope between it and its projection on one the planes, and then taking it with a + sign if its projection is oriented in one direction, and with a - sign if it's oriented in another? This kind of formula works based on the basic geometry.

    • gloveone 13 minutes ago ago

      Yes, the algorithm and its derivations are elementary and rather obvious for anyone decent at undergraduate level mathematics. But still, I am glad to see more people enjoying math!

    • Sharlin 4 hours ago ago

      Yes, this is essentially what the author derived (by means of calculus rather than geometric argument but the result is unsurprisingly the same). The 2D analog is easy to grok: to compute the area of a polygon, find the sum of the signed areas of each of the trapezoids formed by an edge and its projection on the x-axis. Turns out the negative areas of the right-to-left trapezoids cancel precisely out any excess area of the left-to-right trapezoids (or in the case of edges below the x-axis, add precisely the "missing" area).

    • cgadski 4 hours ago ago

      Yep. Using the same kind of calculus ideas, I can also think about a vector field that has a Dirac mass of divergence at some point and zero divergence everywhere else. Then you get an expression that you can sum over faces to determine if a polyhedron contains some point. Again, for the right vector field there is a simple geometric interpretation, namely the solid angle that a face makes with respect to the point.

    • xigoi 4 hours ago ago

      I wonder if this could be reversed to give an intuitive “proof” of the divergence theorem.

      • meindnoch 3 hours ago ago

        The divergence theorem can be intuitively summarized in one sentence: "what comes out is whatever went in, plus whatever was produced inside"

    • joelthelion 3 hours ago ago

      > calculating the volume of a prism-like polytope between it and its projection

      There is the key insight that you don't need to explicitly compute this projection.

    • aaa_aaa 4 hours ago ago

      Yes I remember doing something like that in 90s for a survey/map engineering cad application. After delaunay triangulation, calculating approximate voulume is easy. But this probably is a more general solution

      • srean 2 hours ago ago

        Yes it goes by a couple of names. Surveyor's formula, shoelace formula.

        The Surveyor’s Area Formula Bart Braden The College Mathematics Journal, September 1986, Volume 17, Number 4.

        https://web.archive.org/web/20150406152731if_/http://www.maa...

        • sigbottle a few seconds ago ago

          Oh, it is the shoelace formula, but for 3D?

          OP's description vaguely triggered memories of the shoelace formula from a decade ago, but geometry was never my strong suit. All I remembered was positive/negative triangles (or was it trapezoids)? make magic happen for area calculation.

    • diabllicseagull 3 hours ago ago

      reminds me of that 1994 paper that reinvented the trapezoidal rule

  • srean 3 hours ago ago

    On the other hand, if you want to compute the area of a polygon that have vertices at lattice points, you can count the number of interior points I, the number of boundary points B. Then the area A is

        A = I + B/2 - 1
    
    This is Pick's theorem

    https://en.wikipedia.org/wiki/Pick's_theorem

    one of my favorite results. It does not generalize as nicely to higher dimensions unfortunately.

    If like the post you want the volume of a polyhedron you can use the three dimensional analogue of the shoelace formula (essentially equivalent).

    Let Va, Vb and Vc be the vertices of a triangle ∆ of a triangulation of the surface. You need to name the vertices in a consistent order/orientation wrt the origin.

    Then the volume V is the sum over all such triangles of the signed volumes

       V_∆ = 1/6 Va ^ Vb ^ Vc.
    
    That's the beauty of signed areas and volumes, determinants and exterior algebra.

    To understand why this is so there's this beautiful short video

    https://youtu.be/Sv7VseMsOQc

    • sebastianmestre 2 hours ago ago

      From a computational standpoint, Pick's theorem seems more useful to find the number of interior points via

          I = 2 (A - B + 1)
      
      Where area would be calculated using the sum of signed areas of triangles.
      • srean 2 hours ago ago

        Indeed.

        One of my off by one errors is a stupid hacky Monte Carlo intution for Picks theorem.

        I count the number of points inside. Now about the boundary points I must assign some fractional weight because they are not fully inside. What's a stupid fraction I can use? Well, half seems about right. Voila,

            A = I + B/2.
  • elikoga 5 hours ago ago

    My belly says the naive formula is summing the triangle pyramid volumes to the origin with sign in orientation. It looks like that's what they derived. Which is a generalization of 2d polygon area calculated by summing triangle areas for each edge, I was taught this in a math camp where we calculated map polygon areas on gis data. I remember math knowledge being hard to get pre AI era but I didn't remember it being this hard.

    No idea what the author means by "which are equivalent to rendering the mesh and then sampling the render".

    • less_less 3 hours ago ago

      > My belly says the naive formula is summing the triangle pyramid volumes to the origin with sign in orientation.

      Yeah, that would also work but it's a slightly slower formula, sum(det(v1,v2,v3))/6. This one is summing sort of prism+pyramid shapes made by projecting each triangle to the yz plane.

  • ahaferburg 2 hours ago ago

    The emphasis here is on the mesh being simple and closed. Make sure to validate these preconditions before relying on the output.

    Similar formulas exist for moments, to compute the inertia matrix for a rigid body.

  • meindnoch 3 hours ago ago

    Don't really need vector calculus for this. Geometric intuition is sufficient. It is simply the summation of signed volumes of triangular columns/prisms parallel to the X axis.

    Visualization: https://jsfiddle.net/L7r1hwca/

    I don't know what they could possibly mean by the naïve algorithms with rendering and sampling (???).

    • nyeah 3 hours ago ago

      Yeah, agreed. But words are cheap. It's one thing to say we don't need vector calculus, and another to develop that claim through the actual vector calculus, step by step.

  • arn3n 5 hours ago ago

    I love these kinds of posts. Simple, fast, AI-free, and I learn something new.

    • alex_suzuki 3 hours ago ago

      (2018)

      • IAmBroom an hour ago ago

        OMG, that absolute dummy didn't know something that a human did in 2018!

        Are they not reading the entire internet every morning, when they wake up???

        • alex_suzuki 41 minutes ago ago

          I was merely pointing out that the post was published in 2018 so AI-free is kind of implied.

  • bob1029 3 hours ago ago

    > For a ballpark number, if volume needs to be calculated every frame in a high-performance 60 frames per second application, without the aid of a GPU, only using the CPU capabilities of a $35 Raspberry Pi, around 30 million triangles could be measured every frame.

    If knowing the volume of a mesh is important, we could pre-calculate it (even using this exact technique) and store it as an attribute on the object. Lots of things in game dev that are modeled as an integral over three+ dimensions tend to work better as a baked setup rather than real time. We kickstarted an entire AI industry trying to chase real time lighting.

  • gurkwart 4 hours ago ago

    There's a really elegant solution using Geometric Algebra, that to this day is one of the most satisfying things I've ever learnt. Steven de Keninck outlines it in his 2019 Siggraph talk [1].

    [1] https://youtu.be/tX4H_ctggYo?t=4795

  • Isofarro 2 hours ago ago

    I'm sorry, English is my first language. What does "Hilariously" mean in this context? Or is there a maths specific meaning/interpretation?

    • hallgrim 2 hours ago ago

      The author is just expressing amusement at the surprising simplicity of the resulting algorithm

  • lern_too_spel 30 minutes ago ago

    You might be interested in the shoelace formula and its generalization to n dimensions. https://en.wikipedia.org/wiki/Shoelace_formula

  • gigatexal 5 hours ago ago

    Did they also work on the graphics stack for the Asahi project?

    • gigatexal 17 minutes ago ago

      I used to watch their vtuber videos on YouTube!

    • StilesCrisis 4 hours ago ago

      Yup!

    • unkeen 4 hours ago ago

      Sadly, there is no way to find out, f.ex. by a quick Google search.

  • N_Lens 5 hours ago ago

    I'll accept any kind of jocularity in the current climate!

  • FartyMcFarter 3 hours ago ago

    > (No, there won’t be jokes.)

    I must be missing something here, inside joke or something in the title?

    • Agentlien 3 hours ago ago

      I assume it's because the title contains "hilariously" and the author felt it necessary to state the content wasn't meant as humorous.