Cigarette Smoke - Part 3 - Houdini 19

Since a lot of people are asking me about this setup I thought it would be a good idea to update the setup for H19.

The logic is pretty much identical to part 1 and 2.  The differences are the following:

  • Pyro sim : I'm using a more modern Pyro workflow.
  • Cigarette Smoke Setup : I simplified the logic a bit. More specifically I ditched the Smooth SOP inside the Solver, and replaced it with a way faster Point Wrangle solution.
  • Render : using now Karma XPU ! It's so ... so much faster than Mantra. The volume renders just fine in Karma XPU without any shader assignment. I have very briefly and unsuccessfully tried to assign a volume MaterialX Volume shader following these instructions. If anybody succeeds please let me know !
  • Trivia : I've added a couple of fun advection experiments :
    • advecting a grid to see what happens to the UVs (spoiler : nothing good)
    • advecting rows and lines of the same grid.

Please find the hip file here and let me know if you've questions or suggestion.

Corona Virus Houdini Graph Tools

If you want to skip all the bla bla part and go to the gimme the stuff part, click HERE.
In this crazy time of uncertainty and tragedy caused by Covid-19, the real and only question to answer is really this : when is this going to end ?

I am not an epidemiologist or a data analyst , so please take all that follows with a grain of salt and keeping in mind I'm just a geek in search for answers and I worked on this little project for a couple of days, killing some of the time I've plenty of in these days.
We know that Corona virus (any virus maybe ?) grows exponentially. But what does this really mean ? I found a few answers watching these videos that I strongly recommend, this time not only to learn something new, but also to help giving some sense to immense amount of data concerning this pandemic.

How To Tell If We're Beating COVID-19

Exponential Growth and Epidemics

The first of the two videos is plotting the data in a way that you don't generally find online on Covid-19 statistics web sites: not only it plots data in a logarithmic scale, but it doesn't plot it against time.
Long story short, if I am not mistaken, to convert a regular value-time graph into a new-total (not sure how this type of graphs is names, please let me know if you do), this is how you do it:

X axis = log10 (total value X)
Y axis = log10 (total value today - total value yesterday)

Watching the video it seems like the data plotted this way show clearly when the virus is about to start switching from an exponential growth ... to a Sigmoid or Logistic Curve

More about Exponential Growth vs Logistic Curve here :
Days ago a friend of mine sent me a link to this github archive:

JSON time-series of coronavirus cases (confirmed, deaths and recovered) per country - updated daily since 2020-1-22.

Click on the actual data set below to take a look.

So I decided to create a tool to plot these curves in Houdini, grabbing the data directly from the net.

Disclaimer :
I am not responsible for the data contained in the above mentioned repository. I am just reading it and visualizing it using Houdini. It is very likely that at some point, hopefully soon (and hopefully for a good reason), this repository will stop being updated.

The code to download the data from the github and store it in points is really short. See for yourself.

Once you run this code in a Python SOP, you get ALL the data from the repository stored on attributes in points.

To be more precise, I created one point for each Json item. Each point stores the data in the following attributes:


I wrapped this python script into an HDA named ...

Rona Import Data

This guy needs no parameters. All it does is download the latest data on the github repo. All you need to do is drop the node.
Now we can play with the data.

Rona Filter

Self explanatory.
The parameter "Country (Regex)" accepts a regular expression string.

.* means 'don't filter anything', it's like don't using this node at all |:)
Italy will keep only data about Italy
(Italy|US|Spain) will keep data from Italy and US

Use this node immediately after a Rona Import SOP. Or don't use it at all if you want to keep all countries.

Rona Plot

This is where the fun starts.
And ends too, pretty much, really.
This node plots all the data incoming in the input (curves only, no axis, text etc).

The output of this node are polylines with prim attributes:
prim attr s@__rona_country ( a string indicating the country )
prim attr s@__rona_data (can be "deaths", "recovered" or "confirmed")
prim attr i@__rona_graph (can be 0 for time-value graph, and 1 for trend graph)
prim attr i@__rona_graph_scale (graph Y scale)
Data : you can choose between
  • Deaths - total number of deaths accumulated per day
  • Recovered - total number of recovered accumulated per day
  • Confirmed - total number of confirmed cases accumulated per day
Graph: you can choose between
  • Time-Linear - the good old time-value we're used to
  • Trend-Log - the graph I mentioned earlier. Will show if things are going better or not really.
  • Scale : sometimes the y values are huge. Sadly. We need to scale the Y axis to fit it in the screen.
Filter - Strength : daily values are pretty erratic. Use some filter to smooth the curves and let the meaning emerge.
Filter - Filter Quality : the higher this value, the more the smoothed curve will stay close to the original curve. 

Use this node immediately after a Rona Import SOP, or a Rona Filter SOP.

Rona Decorators

This node takes care of adding axis, text, etc to the curves. Cause without metrics nothing makes sense.

I lied. The TREND graph actually doesn't really need metrics, since the meaningful part is the shape of the curve.
And since I didn't really know what units to write in there.
And I have been a little lazy, I admit.

You can merge together many curves upstream before feeding this node, just make sure the curves are coming from the same graph typs (linear or log) and they all have the same scale.



Plots Examples

These are all contained in the example hip file (link at the end of this article).
Filters Italy, Spain and France and plot the deceased data over time.

Filters the 5 countries with the largest number of deceased and plots the data over time.

Filters the 8 countries with the largest number of deceased, plots the data and creates a grid layout.

Trend plot.

If you're interested you can download the OTLs and an example Hip file from HERE.

Leave comments , opinions and suggestions.
Thank you for reading.

OpenCL Notes

If you're interested in OpenCL please please watch Houdini 16.5 Masterclass | OpenCL by Jeff Lait
After that, I'm sure you'll find most of the stuff I'll write in this post redundant. I'm writing this post mostly as a note to myself with code snippets (cause my memory sucks) and I thought you guys might find it useful as well.

On top of that, I'd like to have your opinion. Please correct me if I say or write something wrong, or you know a better way. I'm really looking forward to learn.

I'll keep updating this same post with new code snippets periodically, so stay tuned.


Ok let's start with an example:
I want to move all the points in my geometry up , in the positive y direction , by 1 unit.

This is how you do it in VEX ...
  • drop a Point Wrangle node
  • connect your geometry to input 1
  • write the following code in the parameter VEXpression:

... now in  OpenCL:

  • drop an OpenCL node
  • go in the bindings Tab and set it up as follows. The purpose of this tab is make geometry attributes available to the GPU. Differently from Wrangle nodes, in order to access the attributes you've to bind whatever attribute you want to read / write to a variable. In this case we bound the attribute P to the openCL variable P ( I could have chosen a different name, like for instance 'ciccio').
    It's important to note that the OpenCL node ny default will loop over every element of the first parameter appearing in the Binding Tab (as specified in the 'Options' Tab) . So , in this case, will loop over P (just like a Point Wrangle) !
  • go in the Kernel Tab and check the parameter 'Use Code Snippet'. This will allow you to write code directly on this node, opposed to point the node to a file containing the code.
  • now press the button "Generate Kernel". A window will pop up with some code. Copy the content of this window, close it and paste it in the parameter 'Kernel Code'.
  • now modify the code as follows:

You should see your geometry moving up by 1 unit.

Let's translate the OpenCL code above in English:

kernel void KernelName ( int P_length, global float * P )
Dear Houdini, please define a kernel function that I'll submit to the GPU if I press CTRL+Enter. Along with this function I'll deliver two pieces of data : one is a long list of numbers (float *) that is coming from the input connection (global) and it's called P (between me and you, this is a list of 3 numbers at a time, cause it's a list of vectors, but you don't have to worry about that for now). The other one is just an integer number (int), named P_length and will contain the length of the list P, so we know when we run out of elements.

{, let's start looping over P using idx as index

int idx = get_global_id(0);  
if (idx >= P_length) 
To start, kindly put the index of the data you're currently looping over in the variable idx. I know you'll be increasing idx at every iteration, but please if you run out of elements to loop over (if (idx>=P_length)...) , just stop (return) so we don't get one of those annoying crashes.

float3 pos = vload3(idx, P);
... fetch the current vector P in a variable called pos ...

pos.y += 1;
... now take whatever number was in the component y of the vector pos , and add 1 to it. Then store it back into pos.

vstore3 (pos, idx, P);
... and finally replace the current P with the content of pos.

... let's go to the next element in P and start a new loop iteration.


I want to put the number -100 in the 3rd element of the array attribute A attached to point with @ptnum 14.


In OpenCL:

Why ?
Well, let's pretend we have this array point attribute named A :

in Vex the syntax would be 

and this is how we're used to see this data in the Geometry Spreadsheet:

in OpenCL all the arrays (one per point) are concatenated and since each array can potentially have an arbitrary length, we need a way to identify where the array for each point begins, in this long list. That's exactly what A_index is for.

A_index[ptnum] contains the beginning of the array for the point ptnum.

So, the code to find the beginning of the array for the point 14 is the following:

int mem_loc = A_index[14]

and then count 3 elements from that memory location, like so :

int element = A[mem_loc+3]

Furthermore, if you want to know the length of the array for the point ptnum, this is the code:

arraylength = A_index[ptnum+1] - A_index[ptnum]


Let's say we have a vector grid called "vel".
First off, make sure this is a Volume Grid (not a VDB). (If you find a way to pass VDBs into OpenCL please I'd love to know how).
In other words we're going to need 3 float grids (vel.x, vel.y, vel.z).
In the example below, I chose to add the Voxel Resolution, because I wanted to know weather or not I was outside of the grid boundaries, but it's not strictly necessary.

When you import a volume grid , for instance in this case it's called 'vel.x', this will result in this code in the Kernel:

int velx_stride_x
int velx_stride_y
int velx_stride_z
int velx_stride_offset
global float * velx ,

If you select "Voxel Resolution" for the grid 'vel.x', this will result in this code in the Kernel:

int velx_res_x
int velx_res_y
int velx_res_z

If you select "Volume Transform to Voxel" for the grid 'vel.x', this will result in this code in the Kernel:

float16 velx_xformtovoxel
(note this is a 4x4 matrix)

Access voxels from a world position
If you want to find the 3 voxel indices starting from a world position P, here's how (this was kindly described in the Masterclass mentioned above).

// first off we need to fetch pos from the P attribute ....
float3 pos = vload3(idx,P);

// now this suuuuper weird formula to find the volume indice
float4 voxelpos = pos.x * velx_xformtovoxel.lo.lo + pos.y * velx_xformtovoxel.lo.hi + pos.z * velx_xformtovoxel.hi.lo + 1 * velx_xformtovoxel.hi.hi;

// and of course the indices must be integers so ...
// NOTE : if you are not 100% sure P is within the grid boundaries you better modify this code to make sure you're not trying to fetch a non existing voxel (Houdini might very easily crash if you do). In this case let's pretend I'm sure P is within the volume boundaries.
int x = (int)(floor(voxelpos.x));
int y = (int)(floor(voxelpos.y));
int z = (int)(floor(voxelpos.z));

// and finally we can access the content of the voxel with indices (x,y,z)
float vx = velx[ velx_stride_offset + velx_stride_x * x + velx_stride_y * y + velz_stride_z * z ];

Detail Attributes

Important :
You need to create the detail attribute that you want to populate in advance, meaning before the OpenCL node, using Attribute Create SOP or Attribute Wrangle SOP.

detail int/float attribute

i@myattrib1 = 12;
f@myattrib2 = 23.4;


The variable myattrib_length returns the number of elements. If this was a point attribute it would be the number of points. Since detail attributes are only one attribute per geometry, for detail attributes this variable is always 1.

detail vector attribute

v@myattrib = set(12.0, 54.0, 34.0);


detail int/float array attribute

Important :
OpenCL cannot resize an array. So the array size must be pre-allocated before the OpenCL node. You can use an Attribute Create SOP or Attribute Wrangle SOP to set the size of the array in advance.

f[]@myattrib = {12.0, 34.2, -4.2, 3};

First off make sure to initialize the detail array attribute with the desired length using , for instance, an Attribute Wrangle SOP. The code below creates an array attribute with length 4:

resize(@myattrib, 4);

then you can finally append the OpenCL node:

The variable myattrib_index is created in OpenCL by default every time there's an array attribute. This variable makes more sense when the array attribute is a point , vertex, or prim attribute, cause in that case every point has his own array, and in this case myattrib_index contains the beginning of the array for each point, vertex or prim. In the case of detail attributes there's only one array, hence myattrib_index is always 0, so can be ignored.

detail string attribute

Looking at the menu that provides the different options for the variable data type, I don't see a "string" option so I must conclude it's not supported.

Sand without grains - PART 2 - VEX Implementation

I'll assume you've read PART 1 so you've a good grasp of this super simple algorithm.

Oh, and in case you're not familiar with HeightFields or Vex , I recommend to watch these two amazing and life changing Masterclasses first.

Let's quickly recap from PART 1 the ingredients:
  • the ideal data structure for this algorithm is a 2 dimensional grid.
  • we need 2 grids: one to store the height, one to store the 'toppling' direction
Now there is one data structure in Houdini that happens to be perfect for what we need :

                                                                   HeightFields (HF)

For those who don't know what HeightFields are in Houdini:
A heightfield is a grid, just instead of being a polygonal grid, it's a volumetric grid. To be more precise it's a 2d Volume. Each voxel of the grid contains a floating point number that corresponds to the height of the visual representation of that voxel.

Let's get familiar with HF first in Houdini.

In a SOP context drop a HeightField SOP node , set size to 100 , 100 and middle click on it to check the info.

You can see the node generates 2 volumes: 'height' and 'mask'.
Please note that even if there are 2 volumes here, in the view port you'll see only the 'height' volume. Mask volume is pretty much identical to Height , but is invisible and in the view port is used to color the Height volume in red for those Mask voxels bigger than 0. Mask is used to select areas of the Height volume to perform arbitrary HF operations. We will use Mask as a debugging tool.

In order to get familiar with these new concepts let's now drop a Volume Wrangle and enter the following code in the "VEXpression" field.

You'll notice that all the voxels that lay in the negative X plane (v@P.x < 0) are now 20 units taller (f@height=20). And you'll notice that all the voxels in the negative Z plane (v@P.z < 0) are now red (f@mask=1).

Remember that every Wrangle node is really a loop that iterates over all the specific elements the wrangle node is set to. In the case of a Point Wrangle the loop will iterate over all the points. In our case, Volume Wrangle , the loop will iterate over all the voxels. We don't have to write any code for the loop to happen, we only have to write the code for what happens for each iteration. So when I say 'current voxel' what I mean is 'the voxel that is currently being looped over'.

In this example we used the voxel World position (v@P) to choose where to change Height or Mask. This means that if you translate or rotate the grid before the Wrangle node, the result will 'slide'.

What if we want the result to be independent from the grid position / orientation ?
Well, we'll need to use a different way to address the voxels : indexing !
The Volume Wrangle node provides several pre-built attributes (v@P is one of those). The ones we want to use now are i@ix, and i@iy (we'll ignore i@iz since we're dealing with a 2d volume, so we'll only consider the components x and y , and always assume the z component as 0).
i@ix and i@iy contain the voxel index. Just what we need.
Ok let's drop a different Volume Wrangle in a new branch and write this code in the VEXpression field.

Since we set the grid to be 100x100 in the HeightField SOP node, with the code i@iy<50 I'm selecting only half of the voxels in the y component to change the Height to 20, and only half of the voxels in the x component to change the Mask to 1.

The result in the viewport did not change at all, but now the result of the Volume Wrangle is not affected by the grid location.
You can test the difference between the two wrangle nodes dropping a Transform SOP for each branch before the Wrangle node, and trying to move the HeightField around.

Now let's try something silly, just because I want to prove a point.

Let's attach a new Volume Wrangle to volumewrangle2 and let's try to color in red (using f@mask) all those voxels that are now higher.
It should be easy since we just set Height to 20 for certain voxels. All we need to do is select all the voxels whose world position y is bigger than, say ... 10 ?

It didn't work !
The reason is explained by the fact that we never changed the actual voxel position. We changed the data *contained* in voxels ('height'), but we never changed v@P.y (the y component of the voxel world position). Even if we wanted to change the voxel v@P.y, we wouldn't be able to do it using a Volume Wrangle node.
If you think about that, it makes perfect sense : a volume is a grid. It wouldn't make sense to be able to move voxels around. If we were able to do so , we would lose the advantages that comes from being ... a grid (we'll see this advantage in a minute).

I don't want to go too deep into HF (please make sure to watch the Masterclasses I mentioned before), but there are a couple of important facts to keep in mind:
  • Even if you see a HeightHield appearing as a 3d object in the viewport , internally it is a 2d grid.
  • The reason cause you see a 3d object, is because the 'height' volume has a special visualizer that shows the 'height' content to visually y-displace the grid in the view port.
  • The y component of a HF voxel world position *is not* 'height'.
  • Voxels cannot be moved in space by a Volume Wrangle. Volume Wrangle can only change the data (for example 'height' or 'mask') contained in the voxel.
  • Voxels can be moved in space with a Transform SOP or the HeightField SOP ('Center' parameter).
  • HF voxels can be uniquely identified by 2d coordinates (x, y). In VEX this coordinate is { i@ix , i@iy , 0 }
  • this is weird : the unit for the 'height' field is 10x bigger than Houdini Unit Length, which is 1 m by default. This means that in order to move a voxel up by 1 unit (or 1 meter) you've to set height = 0.1 (not 1).

Ok enough about HeightFields. Let's dive into the sand algorithm.

In PART 1 I named each cell 'stack'. Since we're now dealing with Volumes, I'll use the term 'voxel' to refer to the stacks.

First off we need to create an environment that will allow us to play with HeightFields. We'll use the Mask volume (generated by the HeightField SOP) as a debugging tool, as we have already done previously.

Cool, let's create a new Geo context, and name it "sandpile_VEX". Dive inside and create the following.

heightfield1 node generates a 20x20 sized grid. Grid Spacing is 1, meaning that there is 1 voxel per unit. If we put 0.5 in there we would have twice the number of voxels. This parameter is a very quick way to increase the HF resolution.

So we have a total of 400 voxels, and the volume indices are as follow:

i@ix = [ 0:19 ]
i@iy = [ 0:19 ]

heightfield_noise1 node modifies the 'height' field using a noise function. This is a quick way to create an initial sand distribution.
In the view port you should see something similar to this.

In PART 1 we split the algorithm in 2 steps. Let's work on step 1.


For each voxel (or stack), find out if any of the 4 adjacent voxels is lower by a certain amount (let's say 1).

Drop a Volume Wrangle node and write the following code in the parameter 'VEXpression'

In the first line we set the topple_threshold. If the current voxel is taller than the adjacent voxels by a height value bigger than 'topple_threshold', it means the current voxel should topple. But don't worry, we'll use this variable later.

Now, the biggest advantage of using volumes as a data structure is that I can find the adjacent voxels very very quickly. The grid position of the 'current voxel' is given by i@ix and i@iy. So if I want to access the voxel to the left, right , top or bottom of the current voxel all I need to do is add or subtract a 1 here and there.

The following 4 lines read the 'height' value from the 4 voxels adjacent to the current voxel.
We'll ignore the voxels in the corners.

Let's translate the following line in English :

Dear Houdini, please read the content of the node attached the first input (@OpInput1), which I assume is a volume. But since there are 2 volumes flowing through that connection ('height' and 'mask') I want to be more specific : the volume I am interested in is named "height", but hey ... not the whole thing. No way. I want the content of one voxel, and more specifically the one sitting at the index ( ix+1 , iy , 0 ). Once you get that data contained in that voxel, please put it in a float variable named 'hxp'.

@OpInput1 is just a fancy way of saying "node input 0". You can write directly 0 in there, and it will still work perfectly. I prefer to use the fancy form cause it's more descriptive for my lazy brain.

Ok now the variables hxp, hxn, hyp, hyn contain the height value of the 4 adjacent voxels. At the moment we're doing nothing with this data.

Let's add some more Vex code:

The purpose of the new 4 lines is to flag which of the adjacent voxels is a candidate for the current voxel to topple towards to.

Let's make an example:
Let's assume the height of the current voxel (f@height) is 3.
Let's assume the height of the adjacent voxel on the positive x (hxp) is 6.
And we know that topple_threshold is 1.
3 - 6 = -3
is -3 bigger than 1 ? False . Then txp = False (or 0).
Meaning, the current voxel cannot topple towards the adjacent voxel on the positive x, cause that voxel is actually taller ! Actually that voxel could topple towards the current, but that will be evaluated in some other iteration of the Volume Wrangle node loop.

Another example:
Let's assume the height of the current voxel (f@height) is 3.
Let's assume the height of the adjacent voxel on the negative y (hyn) is 1.
And we know that topple_threshold is 1.
3 - 1 = 2
is 2 bigger than 1 ? True . Then tyn = True (or 1).
This means that the current voxel could topple towards the adjacent voxel on the negative y, cause it's lower than the threshold value.

Once we collect all the values for txp, txn, typ, tyn we can put this into an array.

The line with hdiff[] is initializing the array to an array of 4 integers. All zeros to start.
Then I'm replacing each element of the array with txp ... tyn, that we know can be only 0 or 1.

Please note that in this moment we just established an important convention. Meaning that the first element in the array corresponds to the upper voxel (yp). The 3rd element of the array corresponds to the lower voxel (yn) and so on. We could have chosen a different convention of course, I just like starting from top and going clockwise. Feel free to go crazy with it, just remember to maintain the same convention later.

Fig 1

In the very last line I sum all the values contained in the array and store it in a new variable called 'topple_candidates'. This variable ultimately answers the most important question : "is the current stack going to topple ?"
If this variable is zero, it means that the array is composed by 4 zeroes, which in turn means that none of the adjacent voxels is a candidate for the current stack to topple.
In other words, if this variable is bigger than zero, the current stack is unstable. And it's really easy now to color code those voxels using the grid 'mask'.
All we need to do is add the following line :

f@mask = topple_candidates > 0;

Which brings the entire code to this :

Now if you look at the view port, you'll probably see this :

Yes it does look like a napkin full of blood. Let's try to ignore that.
This is a very comforting result : note how the voxels that lie on the ground area are all white, cause none of those have adjacent voxels that are much lower. Those white voxels are stable and will not topple. The red voxels instead are on steep slopes and it makes a lot of sense that they might topple down.

But there is one issue we need to resolve before moving to the next step: at the moment if you change the resolution of the grid ( heightfield1 -- Grid Spacing Parameter) you'll notice that the toppling areas (in red) shrink or grow. 

This is not what we want. We want to be able to develop our algorithm on a fast grid (say 20x20), and once we're happy with it, move to a way higher res grid without fearing that the result will change. In other words we need to make our system resolution independent.

First off let's try to understand why this is happening.
When we change the resolution of the grid , we effectively change the size of each voxel. But the 'topple_threshold' height is always 1. Somehow we need to tie grid spacing with topple_threshold.

when grid spacing becomes smaller, that means that the voxel becomes smaller, consequently even  topple_threshold should become smaller.

grid_spacing ---- voxel size ---- topple_threshold

So the parameter Grid Spacing and the variable 'topple_threshold' seem to be connected by a linear dependency. Let's just multiply those, and use the result as our new 'adaptive' topple_threshold.

First off I parametrized a new variable called 'topple_threshould_mult', and created a new parameter called 'grid_spacing' linked to the parameter 'grid_spacing' on the heightfield1 node.
At this point I can finally calculate the new 'topple_threshold' :

topple_threshold = grid_spacing * topple_threshold_mult ;

Now the red areas will stay consistent even if we change the resolution, on top of that we've an additional slider to change the toppling threshold !

Switch to a more decent resolution, something that will show some cool detail. On the node "Heightfield1" set the parameter Grid Spacing to 0.05.

Great !
Now let's get back for a moment to the variable 'topple_candidates'.
We said that if this is bigger than 0 it means the current stack will topple.
Now the question is "in which direction will it topple ?".
In PART 1 I decided to choose the toppling direction randomly. Note that I made this choice cause it's the simplest one. More elaborate solutions would be to choose the toppling direction based on a vector (which could be wind, or forces like mm .. collisions / proximity with other objects). For the sake of simplicity we'll go for the random route, but please feel free to experiment with other solutions on your own.

We need to choose one of the non-zero elements in the array randomly and store the index in an additional grid that we'll name 'dropdir'.

Let's abstract the problem here :

I have an array of 1s and 0s

example :
A = [ 0 , 1 , 0 , 1 ]
I want my code to *randomly* return 1 or 3 (the indices of the array that correspond to 1s).

another example : 
A = [ 1 , 0 , 0 , 1 ]
I want my code to *randomly* return 0 or 3 (the indices of the array that correspond to 1s).

another example : 
A = [ 0 , 0 , 1 , 0 ]
I want my code to always return 2 (the only index of the array that correspond to 1).

A really simple solution is the following (if you find a smarter and/or faster solution please please let me know):
  1. count how many 1s are present in the array (we already have this value, it's the variable 'topple_candidates')
    Example :
    if the array A = [  0 , 1 , 0 , 1 ]  , topple_candidates = 2
    if the array A = [  1 , 0 , 1 , 1 ]  , topple_candidates = 3
  2. if the array A = [  0 , 0 , 0 , 0 ]  , topple_candidates = 0
  3. find a random integer number between 0 and topple_candidates ( let's call it R )
    if topple_candidates = 2 , R = 1 or 2
    if topple_candidates = 3 , R = 1 , 2 or 3
    if topple_candidates = 1 , R = 1
  4. loop over the array until I find as many 1s as R
  5. store the array index in 'dropdir'
Ok first off we need to create an additional grid named 'dropdir'. As we just mentioned earlier each voxel of this grid will contain a number between 0 and 3. This number will give an indication of the toppling direction, according to the convention established in fig.1.

Modify your network this way:

In other words I'm duplicating the field 'height', renaming it to '_toppledir' (I like to ad a leading underscore to temporary data) and merging it back with the other 2 fields 'height' and 'mask'. The 'primitive1' node will disable the visibility of this new grid. For now I disabled this node cause we want to see what happens on this grid when we start storing values. We'll hide it after we are sure it's working properly.

In the view port you should see something like this

Now if you middle click on the merge1 node you should see 3 grids.

  • height - the bloody napkin, now more similar to a bloody mountain
  • mask - not visible as a grid, but visible as red mask on the 'height' field
  • _toppledir - a copy of the height field, but completely flat since all the voxels have a zero value.
Now modify the code in volumewrangle1 this way:

( I know , you hate me cause you can't copy and paste the code. But guess what, I did it on purpose ! 😼 I want you to write it all cause I want you to get used to code in VEX, opposed to become a copy and paste black belt. On top of that formatting code with Blogger is a real pain ... )

Let's quickly go over the new code.
First off, the whole new code is almost entirely enclosed in a big if-statement: if topple_candidates is zero, we mark that voxel as stable ( or -1 ) and we move on to the next voxel.
Otherwise ... 
first off we set a random number 'rnd' between 0.0 and 1.0.

let's focus on this line :
int r = 1 + floor ( rnd * topple_candidates )

The first thing to keep in mind is that the function floor( a ) returns the integer part of a. 
For instance :
floor(0.3) will return 0
floor(1.543) will return 1
floor(12.543) will return 12

now ...
(rnd * topple_candidates) is just a way to re-map a number that is in the range between 0.0 and 0.99999... to a number in a range between 0.0 and topple_candidates-1.

Floor will return the integer part of it, which is ... 0 - (topple_candidates-1)

But, remember that the purpose here is to create a random number between 1 and topple_candidates cause we want to use this number (r) as a counter ... so we just need to add 1. Hence ...
int r = 1 + floor ( rnd * topple_candidates )

In the for-loop below I start looping over the 4 elements of the array hdiff[]. This array will have at least one 1, otherwise we would have been kicked out of the if-statement earlier. 
In this loop , a variable 'c' will be increased every time we meet a 1 in the array. As soon as 'c' reaches the desired count 'r', we have found the index, and we can exit the loop (break statement).
All we need to do now is store the index in '_toppledir' !

Now if you look at the view port you'll see something weird ...

... and kinda beautiful.
The weird looking spiky thing is the '_toppledir' grid. This grid should contain numbers between -1 and 3. In fact, if you switch to the front view , you'll notice exactly that ! 

On top of that, if you look closely, you'll notice the following:
  • areas corresponding to white areas in the 'height' grid are flat and at level -1 
  • areas corresponding to red areas in the 'height' grid are always higher than -1
  • if you scrub the frame in the timeline, the peaks should change randomly (cause the variable 'rnd' is dependent on @Time) , but not the ones in the white areas of the 'height' field.
  • there are a few flat areas that don't correspond to white areas of the 'height' field ! Those should correspond to slopes where there can be only one toppling direction. These areas will be flat, yes, but higher than -1 level.
The good news is that Step 1 is done !
Let's rename volumewrangle1 to step1.

Oh and let's store the value calculated for 'topple_threshold' in a detail attribute, on the geometry itself. This way I don't need to calculate it again in Step 2.

 The VEX line above will use the command 'setdetailattribute' to create a detail attribute

A detail attribute is simply a per-object attribute, opposed to a per-point, per-vertex, per-prim attributes.
I generally use detail attributes to store values that exist for the whole object. They're super convenient cause they're kinda similar to parameters, but attached to the geometry itself, so they flow along the whole graph and you can access it from everywhere. On top of that, they're very light since there's only one per object, and the size is independent on how heavy is the geometry.
Why not call them 'object attributes' instead of  'detail attributes' ? MMMmmm ... well ...


Let's quickly recap.
At this stage we've 3 volumes:
  • Height : contains our sand. It's a 2d volume. Each voxel contains the height of the corresponding stack of sand.
  • _dropdir : contains the info about which stack is stable or unstable. On top of that, for those who are unstable, it even tells us in which direction they'll topple. More specifically (check Fig.1) :
    • -1 : stable
    • 0 : topple up
    • 1 : topple right
    • 2 : topple down
    • 3 : topple left 
  • Mask : used only for debugging purpose. At the moment it flags red all the unstable Height voxels.
So, well ... we know everything. All we need to do is actually do the toppling thing.

All we need to do is this:

For each voxel V in the grid 'Height':
  • find which of the 4 adjacent voxels is marked to topple towards V and increase V's Height accordingly.

Looking at V's 4 adjacent voxels, two of them will topple on V. One is the voxel on the right that contains the value -3 (which, according to out convention, means it's going topple left), and the other is the voxel down that contains the value 0 (which , according to our convention, means it's going to topple up). So a total of 2 'grains of sand' will land on V, increasing it's Height.

  • if V itself is not stable (meaning that it will topple on one of its 4 adjacent vectors) decrease V's 'Height" accordingly. In other words, if the current voxel V itself is not stable (meaning _toppledir contains any value different than -1) it means that one 'grain of sand' is going to topple , decreasing V's Height.

Let's implement Step 2 in Houdini : drop a new volumewrangle node and get the data about the 4 adjacent voxels.

The first 4 lines should be familiar. We're fetching the current voxel's adjacent voxels content for the grid _toppledir that we populated in step 1.

In line 5, I'm just fetching topple_threshold stored in step 1.
Why am I dividing it by 4 ? I'll explain it in a moment.

Cool, now we know in what direction the adjacent voxels will topple. Now we need to check if any of the adjacent stacks are going to topple on the current voxel. In that case we need to increase the current voxel's Height by the topple_threshold.
Let's remember that topple_threshold is the threshold value that decides if a voxel is unstable. So , in step 2, if all 4 the adjacent voxels are toppling towards the current voxel, it means that the current voxel's Height will be increased by topple_threshold*4.
Now, this wouldn't be 'fair' towards the adjacent voxels, because in Step 1, each one of them alone had a reason to topple on the current voxel. If when an adjacent stack topples on the current voxel I increase it by the *whole* topple_threshold amount, this alone would be enough to re-evaluate the reasons for the remaining voxels to topple on the current voxels. We need, instead, to preserve the consistency with the decisions made in Step 1. This is why I'm dividing the value by 4.0.
This way, if all the adjacent stacks topple towards the current stack , Height will never get higher than any of those adjacent stacks.

Cool, all we need to do, is to check if the adjacent stacks have any intention to topple on the current stack. In other words ...
  • if the stack on the right wants to topple on the left increase the current stack by threshold/4
  • if the stack on the left wants to topple on the right increase the current stack by threshold/4
  • if the stack up wants to topple down increase the current stack by threshold/4
  • if the stack down wants to topple up increase the current stack by threshold/4
 ... that's not it yet. We need to check if the current stack is not stable ..
  • if the current stack wants to topple (we don't care where) decrease the current stack by threshold/4
Let's write all this in VEX :

... first off I'm initializing to zero a new variable height_incr where I want to store the toppling contributions coming from all the directions.

All the if... statements are doing exactly what I described earlier.
  • if (dxp == 3) ...
    means ... if the stack on the right (xp = x positive) wants to topple towards the left (remember that according to our convention in Fig.1 the number 3 means 'left').
  • if (dxn == 1) ...
    means ... if the stack on the right (xn = x negative) wants to topple towards the right (remember that according to our convention in Fig.1 the number 1 means 'right').
  • ...
  • if (f@_toppledir!=-1)
    means ... if the current voxel is not stable (meaning _toppledir does not contain -1) then , regardless of what direction it'll topple, has to go down, decrease.
After accumulating all the sand contribution in height_incr , I'll add it to the current height (last line).

Done. Finish. Caput. Finito.


Ok to recap this is what your network should look like

This is the vex code contained in 'step1' Point Wrangle node:

And this is the code contained in 'step2' Point Wrangle node:

Now, if you check your bloody napk .. uhm your sand distribution in the viewport, not much is changed. Well, check better ... it is actually changed, but veeery very little. You need to really zoom in a detail, and if you enable and disable the node 'step2' you'll notice that something is actually changing.

This happens of course because this is only one iteration. We need to repeat this process over and over in order to see the something interesting, and of course Solver SOP is once again our friend.
All we need to do is move the nodes step1 and step2 in a Solver SOP. The Solver SOP will fetch the geometry only at the frame specified in the parameter Start Frame (typically frame 1). For all the subsequent frames, the Solver will use the result of the previous iteration. In other words, will perform a 'Feedback Loop'.

Let's do that. 

  • Create a Solver SOP node, and copy and paste the nodes step1 and step2 inside the Solver SOP network.
  • inside the Solver SOP network connect Prev_Frame input to the node step1
  • inside the Solver SOP , change the expression on the node step2 to

  • inside the Solver SOP, on the node step 1 change the parameter Topple Threshold Mult to 2. This will make sure that only the very steep parts of the sand distribution will be modified, which is more realistic.

  • outside of the Solver SOP connect the node heightfield_noise1 with input 1 of the Solver SOP.
  • outside of the Solver SOP connect the output of the Solver SOP to a delete node, and use it to delete the field '__toppledir'.
    This field was only a temporary storage to store data between step1 and step2.  Once we're out of the Solver SOP we no longer need it.
  • outside of the Solver SOP, Set the number of Sub Steps to 20.

Make sure you're on Frame 1 , and in your ViewPort you should see something like this :

Notice that the mask field is marking in red what elements of the sand distribution are not stable.
Let's press play and see what happens.

It's cool to see how the unstable sand areas get smaller and smaller , that is comforting cause it means the sand distribution is adapting to find a stable state. Let's remove the red areas, so we can actually see the result of all our work.

In the delete node, append @name=mask to the Group parameter to remove the mask feedback from the viewport.

Now the red areas should be gone and we can run a new simulation...
(for the sim below I increased Sub Steps to 100. It's a little slower per frame, but it takes less frames to see the progress.)

Note how the steep area collapse ...

.. and the sand accumulates at the base , creating an advancing leading edge ..

Furthermore, note how the process reaches an equilibrium , and little by little the sand distribution becomes stable and the result is an entirely new distribution.

Advanced ideas for a more production-ready tool :
  • use f@mask to art direct the areas that should be affected by this process
  • use a noise function to create patterns / more interesting toppling effect
  • implement a wind direction
  • emit sand (this time particles) from the non-stable areas to simulate something like this
  • use other objects as 'colliders'. For instance create a sand footstep fx where the sand gets pushed on the sides and accumulates in taller areas around the foot, when the foot sinks in the sand (I haven't implemented this one, but would be cool to see an example of this !) (Thank you for the idea Daniel !)
I hope you enjoyed this article and apologies for taking so long before releasing part-2.
Next time we'll implement this same setup using OpenCL and compare the speed with the VEX implementation.

You'll find this setup HERE.
Note : I've implemented this on H17.0.XXX., but it should work great even on 16.

Featured Post

Cigarette Smoke - Part 3 - Houdini 19

Since a lot of people are asking me about this setup I thought it would be a good idea to update the setup for H19. The logic is pretty much...

Most Popular