May 302017

I’m going to try to make a nice easy introduction to my two favorite functions in Houdini VEX (besides fit01 and chramp of course): xyzdist and primuv. These functions are at the core of a lot of really useful and cool tricks in Houdini, including rivets, the attributeInterpolate SOP, the old “droplets falling down a soda can” effect, and some really awesome stuff with volume shaders. I’ll do a little example of each as a way of showing off what you can do with these clever little tools.

First, let’s take a look at the VEX definition (the third overload here is the most frequently used):
float xyzdist(string geometry, vector pt, int &prim, vector &uv, float maxdist)

At its most basic, xyzdist will return the distance from the sample point pt to the nearest point on the surface geometry. Note that this doesn’t mean the nearest actual point, but the interpolated surface in between those points.

Those little “&” symbols mean that this function will write to those parameters, rather than just read from them. So if we feed this function an integer and a vector, in addition to the distance to the surface, it will also give us the primitive number prim and the parametric UVs on that primitive uv. Note that parametric UVs are not the same as regular UVs… this just means the normalized position relative to the individual primitive we found.

So, what can we do with this? Click below to find out… Continue reading »

Dec 292016

My good friend and motion graphics bromance, Eli Guerron, asked me to help him create a procedural system of branching circuits that would look something like a schematic drawing. I stupidly thought this would be an easy trick with particles, but to get the right look it actually took quite a lot of experimenting.


The reference video he showed me showed branching structures that would occasionally turn or split apart at precise angles. This part was easy enough to figure out, but the real trick was preventing collisions between branches. I first thought that I could use a modified version of the space colonization algorithm to create the non-intersecting patterns, but even after writing a function to restrict the growth directions to increments of 45 degrees (explained below), I couldn’t get the patterns to look right. Part of the reason the reference image looks “schematic-y” is because the traces will run alongside each other in parallel, sometimes splitting apart but generally running in groups. Space colonization just doesn’t grow that way, so I had to throw that method out.

The method that worked the best (it’s still not perfect) is just growing polylines in a SOP Solver, completely manually in VEX snippets. Each point has a starting velocity vector (this is left over from when I originally tried to solve this with particles), and on each timestep the endpoints will duplicate themselves, then add their velocity to their point position, and draw a polyline between the original and the duplicate. This forms the basis of the traces. The code for this is pretty straightforward:

// create new point, inheriting attrs from this point
int newpt = addpoint(0,@ptnum);
// assign new position to new point
vector newv = normalize(v@v) * speed;
vector newpos = @P + (newv);
// add primitive to connect new point to current point
int newprim = addprim(0,"polyline");
// remove old point from "growth" group...
// this can be done by just setting a point attribute

The next step is to get the traces to randomly zigzag at 45 degree angles. Since I’m dealing with vectors here, rotating vectors is most easily done (IMO) via a matrix. The math really isn’t too bad, but there’s a few steps that need to happen to make this work.

First, I need to create a generic identity matrix. This is just a matrix for the point that, when multiplied against the point position (or any other vector attribute, such as velocity), returns the same vector right back. Creating an identity matrix is easy:

matrix3 m = ident();

Next, I want to generate an angle of 45 degrees. When dealing with vector math in general, you typically want radians, not degrees. 45 degrees is equivalent to π/4 radians, which I can then randomly negate to create either a positive or negative angle.

matrix3 m = ident();
float PI = 3.14159;
float angle = (PI * 0.25);
if(rand(@ptnum*@Time+seed) < 0.5) {
    angle *= -1;

Now I need to define an axis to rotate my angle around. Since my schematic is growing on the XZ plane, the axis to bend around would be the Y axis, or the vector {0,1,0}. Once I have the axis, the angle, and a matrix, I can use VEX’s rotate() function to rotate my identity matrix. All you have to do at that point is multiply any vector by that matrix, and it will rotate to the specified angle. Keep in mind that for predictable results, your vector should be normalized before multiplying it against a rotation (3×3) matrix. You can always multiply the resulting vector against the length of the original vector, once you’re done.

matrix3 m = ident();
float PI = 3.14159;
float angle = (PI * 0.25);
if(rand(@ptnum*@Time+seed) < 0.5) {
    angle *= -1;
vector up = {0,1,0};
newv = (normalize(newv) * m) * speed;

The branching step is almost exactly the same. I just randomly generate a new point and rotate its velocity by +/- 45 degrees, but don’t remove the original point from the growth group. (In my example HIP, the “tips” growth group is defined by being red, while the static group is green. The “tips” group is redefined at the beginning of each step based on color.)

This setup alone gets some pretty cool patterns, but it’s messy without collision detection:

The collision detection is where things get tricky. I wanted growth points to look ahead of themselves along their velocity vector, and if they detected a collision, to try to match the collision surface’s velocity. This process repeats in a loop until the point can verify that nothing is blocking future growth. If the loop exceeds a maximum number of tries, it kills the branch by removing the point from the growth group.

The easiest way to “look ahead” and grab the velocity from collided polylines was to use the intersect() function. This returns a primitive number (or -1 if there was no collision), and parametric UV coordinates at the collision site. The primitive number and UV coordinates can then be fed into the primuv() function, which can be used to grab any attribute value at those exact coordinates. The velocity value at the collision point can then be assigned to the growth point so it moves in the same direction.

// find intersection along velocity vector.
// if an intersection is found, inherit that point's
// velocity attribute.
// repeat this up to "maxtries" to resolve collisions.
// if there is still a collision when maxtries is reached,
// remove this particle from the tips group (turn green).

vector interP;
vector forward = {1,0,0};
vector up = {0,1,0};
float PI = 3.14159;
float interU;
float interV;
int maxtries = chi("max_tries");
int try = 0;

while(try <= maxtries) {
    int primnum = intersect(1,@P,v@v*4,interP,interU,interV);
    if(primnum != -1) {
        // we collided
        // get velocity of collision point
        vector uv = set(interU,interV,0);
        vector collideV = primuv(1,"v",primnum,uv);
        v@v = collideV;
        // if this was the last try, this line has to stop
        if(try==maxtries) {
            // turn green to remove from tips on next step
            @Cd = {0,1,0};

This works… okay, but because the points sometimes collide in between points that are at angles to each other, every once in a while the primuv() function will return angles that are a little “fuzzy,” meaning not at 45-degree angles. It doesn’t happen often, but when it does happen, you start seeing some weird curvy lines. So as a last step, the velocity vectors need to be locked to 45-degree angles, just in case.

First, I define a forward and up vector. Forward is +X (1,0,0), up is +Y (0,1,0). The forward vector is to have an axis to compare the velocity against… since this schematic is being drawn on the XZ plane, the +X axis is just a convenient starting point for figuring out rotations. To rotate angles on this XZ plane, I use the Y axis (up).

Figuring out the angle between two vectors is straightforward. The formula is this:

float theta = acos( dot(V1,V2) );

This returns the angle (typically called theta) between vectors V1 and V2, in radians. Because of the way this calculation works, it’s only going to properly return angles for vectors pointing towards +Z, so if the velocity has a negative Z component, we just multiply the angle by -1.

float angle = acos( dot( normalize(v@v),forward) );
if(v@v.z < 0) {
    angle *= -1;

Next, we need to figure out how far away our angle is from a 45-degree angle. We can do this by using the modulo operator to figure out the remainder when we divide our angle by (π/4).

float angle = acos( dot( normalize(v@v),forward) );
float rem = angle % (PI * 0.25);
angle -= rem;
if(v@v.z < 0) {
    angle *= -1;

Now we have a nice 45-degree-divisible angle relative to +X (or 0 radians). We need to turn this into a unit vector, which we can then multiply against our original velocity vector’s length to get the final velocity.

Check out this circle diagram to see what we’re doing here:

Substitute Y with Z and you can already see the exact formula we’re going to use. Thanks to the unit circle, we can easily plot out where on the circle we’ll be with a little more trig. Then we just have to multiply that against the original velocity’s magnitude (length).

float angle = acos( dot( normalize(v@v),forward) );
float rem = angle % (PI * 0.25);
angle -= rem;
if(v@v.z < 0) {
    angle *= -1;
vector newV = set(cos(angle),0,sin(angle)) * length(v@v);
v@v = newV;

The final result is pretty cool! I still don’t like how some of the branches get too close to each other, and points can still intersect each other when two tips collide with each other on the same point, but it’s pretty close to the original setup. I’d like to eventually figure out a better way to keep the spacing between traces a little more consistent (right now they tend to sometimes bunch up tighter than I’d like), but that will have to happen later on.

Here’s the HIP file, if you want to play with it!


Apr 012016

As part of the pipeline I’m currently writing for Timber, I’ve been working on an input/output system that can pass data back and forth between Maya and Houdini. The goal is to have a Maya artist send an Alembic scene to Houdini, process the scene in Houdini, then send the modified scene back to Maya. The problem is when you want the scene to be recognizable at all, and also easy to work with in Houdini.

There are several steps to making this process work, and you’re better off scripting as much as you can because it’s annoying and repetitive. Also, the method I’m using will preserve the scene hierarchy, but it will not preserve pivots or transform information (everything becomes baked to world space). Hopefully there will be a better way to manage that soon.

First off, when you import the Alembic, use an Alembic Archive geometry container. Disable “Build Hierarchy Using Subnetworks” and then build the hierarchy. Instead of packing each child transform inside subnetworks for the whole Alembic hierarhcy, you should get a flattened hierarchy inside your Alembic Archive that looks like Fig. 1.

Fig. 1. A flattened Alembic subnet. All the geometry exists inside this subnetwork, with no child subnetworks.

Fig. 1. A flattened Alembic subnet. All the geometry exists inside this subnetwork, with no child subnetworks.

The reason you want this flattened hierarchy is so that you can merge everything into another network easily. I prefer to import all my data in one place, then Object Merge it into another network to modify it in order to keep everything nice and readable. Since everything is flat, on the Object Merge node you can just merge in /obj/alembicarchive1/* and everything should appear nicely. In order to preserve the original hierarchy for export, though, we’re going to want to enable “Create Per-Primitive Path” so that a primitive attribute objname appears on every packed Alembic primitive we’re reading into the scene.

This objname parameter isn’t storing exactly what we want, though. It’s storing the path to the SOP in Houdini’s terms, not the path the object had in the original export. If you look inside the Alembic Archive you originally created and click on any of the Alembic nodes inside, you can see that there’s an objectPath parameter that is storing the original Alembic path. We just need to get this value for each primitive. Since we have objname already stored, it’s not too hard to get objectPath from each packed primitive. We just need to use a little VEX (don’t worry).

Drop down an Attribute Wrangle after your Object Merge, and try the following code:

string objname;
objname = prim(geoself(),"objname",@primnum);
// get objectPath attribute from the parent node of objname
string channame = sprintf(objname + "/objectPath");
string objpath = chs(channame);
s@abcPath = objpath;

First we fetch the Houdini object path (stored as objname) on each primitive. Next, we generate a string for the channel name that we want to query, which is objname + "/objectPath", the channel name that contains the original Alembic path of the object. We then use chs() to grab the value of that channel, and store it as a primitive string attribute called abcPath. Not too bad. If you look at the object now in the Geometry Spreadsheet, you can see that each primitive now has a string attribute pointing to the Houdini path as well as the original Alembic path of each object.

There’s another little catch here. We’re probably going to want to convert this geometry to regular Houdini geometry, since there’s not a lot we can do with packed primitives (or polygon soups, depending on what you’re trying to do). So we’ll need to append an Unpack SOP. When you unpack, though, you’ll notice your primitive attributes are gone. To transfer the primitive attribute onto the newly unpacked polygon soups, just enter abcPath into the Transfer Attributes parameter on the Unpack SOP. Now every primitive should have the correct attribute, and you can go ahead and drop down a Convert SOP to convert everything into standard Houdini geometry.

You can go ahead and smash up the geometry to your heart’s content now, as long as you maintain those primitive attributes. When you’re done and it’s time to export an Alembic, all you have to do on the Alembic ROP is select “Build Hierarchy From Attribute” and type in abcPath as your Path attribute. If you were to load the exported .ABC file back into Maya, you should see that the object names and their hierarchical relationship to each other should be unchanged (although your pivots will all be at the origin, and the transforms frozen). Now your Maya team won’t freak out when everything in your Houdini scene comes in as one giant unrecognizable object!

Oct 152015

I realize this effect was done years and years ago by smarter people than me (JT Nimoy, Adam Swaab, among others) but I’m a big nerd for procedural art and wanted to see how I could create the Quorra Heart effect described on Nimoy’s blog here. The effect itself isn’t terribly complicated to do in Houdini, but it took me a while to figure out the best approach (thanks Ray).

Originally I wrote a function in VEX to sample a volume in cross-sections. The volume was just a sphere multiplied against a 3D Perlin noise function. I stepped through the volume grid in a loop, checking to see if the volume’s value was close to the isosurface value I wanted, and then added points where it was. Here’s the code:

float minX = -.5;
float maxX = .5;
float minY = -.5;
float maxY = .5;
float minZ = -.5;
float maxZ = .5;
float interval = 0.004;
float isovalue = 0.001;

int zslices = 15;

// calculate z interval
float zint = abs(minZ - maxZ) / zslices;

for(float x=minX; x<maxX; x+=interval) {
    for(float y=minY; y<maxY; y+=interval) {
        for(float z=minZ; z<maxZ; z+= zint) {
            vector pos;
            pos.x = x;
            pos.y = y;
            pos.z = z;
            float s = volumesample(1,"density",pos);
            if(abs(s) < isovalue) {
                int newpt = addpoint(geoself(), pos);
                setpointattrib(geoself(), "sdfvalue", newpt, s);

Nothing too fancy there. Since the volume here is an SDF, I’m establishing a boundary condition around the isosurface value (the value of an SDF volume at the surface is zero) and checking to see if the current voxel is really really close to zero. If it is, I draw a point. That’s it.

This method outlined the shape in an interesting way, but it wasn’t even close to the original effect, and the points weren’t connected to each other which was problematic. I was looking at the problem the wrong way; I didn’t want to sample the volume by planes, I wanted to sample the volume by spheres!

I decided to scrap the VEX approach and do this visually. Houdini’s Cookie SOP can extract curves from the areas where two objects intersect. I just had to create a bunch of densely-packed concentric spheres, extract the intersection curves using the Cookie SOP with the volume converted to a polygon mesh, and then it was done! Very simple effect when it comes down to it.

For the final animation below, I had both concentric spheres and a stacked group of planes intersecting with the volume. The animation is just the Flow parameter of the Flow Noise VOP keyframed from zero to one… this means the animation can loop seamlessly. There’s also a little bit of VOPs work to add color to the intersection points; it’s just a distance lookup between the point’s position and the center of the sphere, fit to a 0-1 range and used to drive a ramp. You could optionally drive the point color based on the point’s distance to the SDF surface by using a Volume Sample From File VOP in a similar manner (the closer your volume sample is to zero, the closer you are to the surface).

Here’s the animation below.

isosurface intersections test 2 from Toadstorm Inc on Vimeo.

I’d like to eventually try to get this working on the GPU in TouchDesigner, but I’m still a total beginner to GLSL so it’ll be a while before I have the chops. If anyone has any tips there I’d love to hear them.

Jan 212015

I’m working on a production right now that involves sending absolutely enormous animated meshes (6.5 million polygons on average with a changing point count every frame) out of Houdini and into Maya for rendering. Normally it would be best to just render everything directly out of Houdini, but sometimes you don’t have as many Houdini licenses (or artists) as you’d like so you have to make do with what you have.

If the mesh were smaller, or at least not animated, I’d consider using the Alembic format since V-Ray can render it directly as a VRayProxy, but for a mesh this dense animating over a long sequence, the file sizes would be impossibly huge since Alembics can’t be output as sequences. Trying to load an 0.5 TB Alembic over a small file server between 20 machines trying to render the sequence and you can see why Alembic might not be ideal for this kind of situation.

Hit the jump below to see the solution.
Continue reading »

Sep 262014

Side FX added a new Cloud FX toolkit a version or two ago, and I recently had a chance to mess around with it. Their Cloud Rig shelf tool is pretty great out of the box… pick a shape, turn it into a cloud, done. The Cloud Noise and Cloud Light SOPs that are built into the rig setup can get some pretty good results, but it’s not exactly what I wanted… I was looking for a solution that would be a little less dependent on the volume container resolution, and more based on textures instead. It’s not necessarily the greatest solution for swirling, dynamic simulations, but for more-or-less static clouds, it gets resolution-independent nice results that are quick to generate. Maya’s fluid shader supports textures by default, but the Cloud Shader that Houdini uses doesn’t have much in the way of textural control, so I had to make some customizations, and that’s what this post is about. Check the link below to keep reading…

Continue reading »

Sep 112014

One of the bigger challenges with rendering liquids is that it can be difficult to get good UVs on them for texturing. Getting a displacement map on a liquid sim can make all the difference when you need some added detail without grinding out a multimillion-particle simulation. Unfortunately, liquid simulations have the annoying habit of stretching your projected UVs out after just a few seconds of movement, especially in more turbulent flows.

In Houdini smoke and pyro simulations, there’s an option to create a “dual rest field” that acts as an anchor point for texturing so that textures can be somewhat accurately applied to the fluid and they will advect through the velocity field. The trick with dual rest fields is that they will regenerate every N seconds, offset from each other by N/2 seconds. A couple of detail parameters called “rest_ratio” and “rest2_ratio” are created, which are basically just sine waves at opposite phases to each other, used as blending weights between each rest field. When it’s time for the first rest field to regenerate, its blend weight is at zero while the rest2 field is at full strength, and vice versa.

It’s great that these are built into the smoke and pyro solvers, but of course nothing in Houdini can be that easy, so for FLIP simulations we’ll have to do this manually. Rather than dig into the FLIP solver and deal with microsolvers and fields, I’ll do this using SOPs and SOP Solvers in order to simplify things and avoid as many DOPs nightmares as possible.

Here’s the basic approach: Create two point-based UV projections from the most convenient angle (XZ-axis in my case) and call them uv1 and uv2. As point attributes, they’ll automatically be advected through the FLIP solver. Then reproject each UV map at staggered intervals, so that uv2 always reprojects halfway between uv1 reprojections. We’ll also create a detail attribute to act as the rest_ratio which will always be 0 when uv1 is reprojecting, and 1 when uv2 is reprojecting. It all sounds more complicated than it really is. Here goes…

Continue reading »

Sep 082014

I ran into a problem recently where I was trying to make some nice-looking embers in houdini, complete with nice motion-blurred trails. Typically with a particle system you use the velocity attribute to handle motion blur, but geometry velocity blur is always linear, so your motion trails will always be perfectly straight even if you have nice squiggly motions with your embers.

Deformation motion blur looks great, but in most simulations particles are being born and dying all the time, and deformation motion blur doesn’t work with a changing point count.

The solution is to force a constant point count. This can be problematic when your particles need to have a lifespan, so there are a few little tricks you’re going to have to pull in order to make this work…

Continue reading »

Nov 192013

I just wanted to post a little bit about the Creep POP in Houdini because I’ve had such a hell of a time getting it to do anything useful.

I was trying to create the effect of raindrops sliding down glass. In Maya, this is a pretty simple thing- enable Needs Parent UV on the emitter, create Goal UV attributes on the particles, and set goalU=parentU and goalV=parentV on the creation expression. The particles stick to the object, and then you can work from there.

Houdini, as it usually does, makes things a little more complicated. The Creep POP looks for two important attributes: POSPRIM and POSUV. It wants to know what primitive to stick to, and what UV coordinates on that primitive to stick to. Of course, when emitting a particle randomly on a surface, you rarely have any idea what those numbers are, so we’ll have to use some expressions to create them manually.

Here goes. On the Source POP you’re emitting from, set up a birth group (leave Preserve Group off) and call it something like “justBorn.” We’ll create the attributes only on the brand new particles so they know what their starting position should be for the Creep. Next we’ll create two Attribute POPs. The first one we’ll use to create POSPRIM, so we know what primitive on the source geometry the particle was emitted from. On the first Attribute POP, set the Source Group to “justBorn”, the type to Integer, and the name to “posprim.” The expression for the value looks like this:


Here’s what this is doing. xyzdist() can return a number of things about an object based on a position in space; in this case, we’re querying based on the position of each new particle, $TX $TY $TZ. We want to get values based on the object we’re emitting from, which is ../path_to_emitter. The -1 means that we just want to return values based on the nearest primitive we can find. If we used any other value, it would return values specifically from one primitive. The 3 means that we just want to return a primitive number, and nothing else.

If you check the Details View now and run your simulation a bit, you should see that each new particle now remembers what primitive it was emitted from via our new POSPRIM attribute.

Now we have to take care of POSUV, which is done in a very similar manner. On your second Attribute POP, again set the Source Group to “justBorn.” This time the type is Float and the size is 2 (since we’re dealing with a UV coordinate) and the name is “posuv.”

Here’s the expressions for this one. The first expression goes into the first Value slot, and the second expression goes to the second Value.



It’s almost exactly the same expression every time… the difference between these expressions and the first one is that now we’re looking for U and V coordinates for a specific primitive based on our location. Since we have POSPRIM to tell us which primitive we should be sampling for every point, we can use the values 1 and 2 at the end of the expression to return U and V coordinates, respectively, based on our position. Now every particle knows what primitive it was emitted from, and the UV coordinates it was emitted from exactly. Check the expression help for xyzdist() to get a better idea of what kinds of things you can return from this expression.

Now you can just drop down a Creep POP and use the default local variables $POSPRIM and $POSUVU/V for the Prim Number and Prim UV. If you want to use forces to drive the particles, set the Behavior to Slide and use whatever force you want to push the particles along the surface. It’s not so bad once you do it a few times; it’s just too bad this behavior doesn’t have a quicker setup.

Apr 152013

My latest gigantic Houdini project is finally live! I was the Technical Director for this one, which also means Houdini tube effects, lighting, shading, rendering, etc. Thanks to my fellow Houdini artist, Alvaro Segura, for handling the inky effects, as well as for helping me to refine the tube generator OTL I spent so long constructing.