Analytical Foam

comments 17
Free Tutorials

As position based simulation techniques are very fast theses days, it is possible to use them to confrom intersecting spheres to each other, to model foam. Nevertheless there are other possibilities, too.

In this tutorial, Manuel shows you how to come up with an analytical operator that removes the intersections from a dense sphere packing by deforming the intersecting points. This is even faster than PBD and can be used without any simulation at all. The video deals with implementing the setup. The particle effect from the demo is included in the hip file though.

Paul Bourke’s site on the intersection of circles

Download the hip file here

Liked it? Take a second to support Manuel on Patreon!
Become a patron at Patreon!


  1. brilliant, i’ve not got the brains to do this myself but have been looking at sops waays for years! so easy when you are sitting next door! thanks

  2. Very enjoyable tutorial. Thanks very much for sharing your work!

    I think it would be interesting to extend this so that the bubbles’ centers are scattered not only on the ground plane, but in a volume. They could inflate as a function of their height as well, from bottom up. In this case, the intersection plane wouldn’t be perpendicular to the ground plane, so the transform in and out of “local space” would be a little more involved.

    (and also, to cut off the bottoms of all bubbles, with their intersections with the ground plane)

    Again, a great tutorial, a good build on your last one, too. Thanks!

    • Hi Maury. That is supported. Just try it. The calculations are not restricted to a plane. Itβ€˜ll work in a volume, too. Cheers mnu

      • You’re right, you can fill VDBs with foam!

        I noticed that some intersections were getting small “collars”, if sectdist was negative. This occurs when a smaller sphere’s center is within it’s larger neighbor’s radius. (You have to look pretty closely to see this, it’s usually quite a small thing, though occasionally it makes largeish “disk” shapes)

        Here is one possible way to solve this (very small) issue:

        In the “gather data” wrangle:

        int neighbours[] = nearpoints(0,@P,maxdist*2);
        pop(neighbours,0); //eliminate current point from array
        int intersecting=0;
        float sectdist[];
        matrix xform[];
        //**** new variable:
        float intersection_radius[];
        foreach(int i; int neighbour; neighbours)

        float n_pscale=point(0,”pscale”,neighbour);
        vector n_P=point(0,”P”,neighbour);
        vector between = n_P-@P;
        float dist = length(between);
        float clear = n_pscale+@pscale;

        if (dist<clear)
        float mysectdist = (pow(@pscale,2) – pow(n_pscale,2)+pow(dist,2))/(2*dist);
        append(sectdist, mysectdist);

        // calculate intersection radius, to remove "collars"
        //myintersection_radius corresponds to "h" on the drawing from paul bourke's website:
        //h=sqrt(r^2-a^2) where r is pscale, a is mysectdist
        float myintersection_radius = pow( pow(@pscale,2)-pow(mysectdist,2),0.5);

        vector zaxis = normalize(between);
        vector tmp= normalize(cross(zaxis,{0,1,0}));
        vector yaxis=cross(zaxis,tmp);
        vector xaxis=cross(zaxis,yaxis);

        matrix myxform=set(xaxis, yaxis, zaxis, @P);
        myxform.xa = 0;
        myxform.ya = 0; = 0;

        f[]@sectdist = sectdist;
        4[]@xform = xform;
        //******* new attribute:
        f[]@intersection_radius = intersection_radius;


        And in the "deform" wrangle:


        float sectdists[]=point(1,"sectdist",@id);
        matrix xforms[] = point(1,"xform",@id);
        //************ new var:
        float intersection_radii[]=point(1,"intersection_radius",@id);

        foreach (int i; matrix xform; xforms)
        //*********optionally highlight collar problem instances
        // there’s an intersection
        //local_P=set(local_P.x,local_P.y, sectdists[i]);
        //************************ new code:
        if (sectdists[i]intersection_radii[i]) //this is a “collar” point, as it extends past the intersection radius
        local_P*=(intersection_radii[i]/dist); // scale this point back to the intersection radius
        local_P.z=sectdists[i]; //move this back to the intersection plane.

        //********* if you’d like a flattened bottom, like a ground plane:


        • thanx for this addition! Now it’s even more flexible!

        • I’m still fairly new to coding, so my problem solving skills in this area aren’t the best. I’ve been trying to institute this new code in my scene to remove the collars and flattened spheres, but I can’t get it to work. I feel like there’s something missing or maybe a typo, but I can’t figure out what it is. I do think it’s something in the Deform wrangle, though.

  3. Another great one! However, those of us who get an anxious feeling when we see the typo that the instructor has made and the feeling doesn’t go away until they correct it… minute 29 to 30 was a minor form of torture πŸ™‚

    • I know exactly what you mean. I feel the same when someone does this and I apologize for having so many typos in the code πŸ™‚

      • Hah, no need to apologize. It was a fun and funny part of the tutorial. To be honest, it might have actually cured me of “typo anxiety” as I believe that minute desensitized me to it :p

  4. Hello

    Thank you for a incredible yet simple lecture. Can you maybe recommened me some books with which I can start exploring the backgroung behing your process? I have good understanding of basic geometrie, vextors and medium level mathematics.


  5. invert matrix and local_p= @P*invert matrix please help me about why is it

  6. Hi Manuel, thanks for another great tutorial, but I do think the same can be achieved a little easier by directly working with the intersection planes between spheres. This removes the need for matrix constructions and inversions and it will also be much faster, especially when using larger numbers of spheres.

    The idea is to use the so-called “normal form” of a plane, i.e. a 4D vector storing the plane normal and the plane’s distance from the world origin and then later moving/forcing/projecting impacted points back onto this plane, just like you did via the matrix approach. The nice thing is that this only requires a 2 dot products (one in gather, the other in deform), rather than all this matrix malarkey… πŸ™‚

    For reference:

    // gather code
    float rp = @pscale;
    int neighbors[] = nearpoints(0, @P, chf(“../set_pscale/max”) * 2);
    vector4 planes[];

    for(int i=1, n = len(neighbors); i<n; i++) {
    vector q = point(0, "P", neighbors[i]);
    float rq = point(0, "pscale", neighbors[i]);
    vector dir = q – @P;
    float dist = length(dir);
    if (dist 0) {
    // project back onto plane
    @P -= normalize(n) * d;

    Another tiny thing: I’d recommend updating the attrib settings in “copy to points” – by default, all these array attribs you’ve created in “gather” will be copied to each point of each sphere, which is very wasteful, since all you need is the “id” attrib to look these up in “deform” from the original seed points. Not nitpicking, just trying to be helpful πŸ™‚

  7. I’m new to houdini so forgive me if this is a silly question but, in the attribute VOP that sets the pscale. You say the ramp is used to remap the values between 0 and 1. How come the first fit doesn’t output values between 0 and 1? I don’t really get why the ramp is necessary.

  8. hi, I can’t download this ‘properly’, when I unzip it thearchive, the .hip file is missing from the resulting folder structure

Leave a Reply