# Analytical Foam

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

1. harry

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. Maury

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!

• Manuel

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

• Maury

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:
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

//*************************
intersecting++;
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;
myxform.za = 0;
append(xform,myxform);
}

}
i@intersecting=intersecting;
f[]@sectdist = sectdist;
i@id=@ptnum;
4[]@xform = xform;
//******* new attribute:

//**********************

And in the "deform" wrangle:

//*******************
//*******************

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

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

}
//********* if you’d like a flattened bottom, like a ground plane:
/*
if(@P.y<0)
{
//@P.y=0;
}
*/

//******************************

• Manuel

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

• Brian

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. Kevin

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 π

• Manuel

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 π

• Kevin

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. Edis

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.

Cheers

5. ibrahim

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… π

// 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. Nizzac

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. Giorgio

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