Marine-Snow-Simulation in Houdini

Don’t lose your bearings!

One feature I wanted to add was not directly related to the simulation, but to the behaviour of the objects that are bound to the particles for rendering. It would be nice if the Marine Snow wasn’t just a mass moving up and down, but if the particles also had a slight rotation or orientation of their own.

Calculating the orientation itself was less of a problem, as there are numerous good explanations and examples on the Internet. All approaches are based on the calculation of an angle that changes over time in order to obtain a rotation. However, what is not yet taken into account here is the axis around which a particle rotates. The normal vector of a particle seems to be the most obvious choice here, except that POP particles have no normal at all. A corresponding attribute can be created with a VEX line, e.g. with the Y-axis as the normal vector:

@N = set(0,1,0);

However, this does not go far enough, as all particles would have the same normal and therefore rotation axis, which would lead to an unnaturally synchronised movement. It would therefore be appropriate to give each particle a random axis of rotation that takes all three directions XYZ into account.

In VEX, the rand() function, which returns a value between 0 and 1, is usually responsible for randomness. This suited my requirements very well, as the components of the axis vectors are also between 0 and 1. So that not every particle is always assigned the same vector, the particle number @id can be used as a seed to generate the random numbers. And to generate different values for each axis, a different seed is required for each XYZ. This is achieved by simply adding or multiplying any value to @id.

vector axis = set(rand(@id), rand(@id * 1.1), rand(@id * 1.2)

This is a good approach, but it has a weak point. The axis and angle are linked to determine the orientation. However, a random vector can distort the result. Therefore, another transformation is used that has already been used more frequently: normalise().

vector axis = normalise(set(rand(@id), rand(@id * 1.1), rand(@id * 1.2))

The problem with this formula is that the values are always positive, which creates a certain tendency in one direction. By subtracting 0.5, the values are shifted into a symmetrical range between -0.5 and 0.5. This means that the rotation axes are distributed in all directions.

vector axis = normalise(set(rand(@id), rand(@id * 1.1), rand(@id * 1.2)) - 0.5);

What I still needed was the rotation angle. This angle should change over time, i.e. become larger, but never be too high. After all, the particles should not be travelling on a merry-go-round, but should be moving slowly. Furthermore, the speed had to be slightly different for each particle, because I didn’t want a ballet either.

A good approach seemed to be to make the angle dependent on the so-called time increment @TimeInc and to introduce a new rotation_speed parameter. @TimeInc is the current time increment of the simulation, not the time itself! I have set 0.1 as the default value for rotation_speed.

float angle = @TimeInc * rotation_speed * fit01(rand(@id), 1, 3);

Quaternions are required to determine the orientation. I will save myself the explanation of this number range. However, it is important to mention that quaternions are particularly important for rotations because they bypass the so-called gimbal lock. With gimbal lock, two axes of rotation are parallel in certain orientations, which leads to errors. To obtain the orientation, I was able to use two VEX functions, which saved me tedious matrix calculations.

vector4 orient = p@orient;
vector4 q = quaternion(angle, axis);
orient = qmultiply(orient, q);

p@orient = orient;

I could only judge the result of the orientation calculation after I had bound a box to the particles. Keyword: Copy to Points SOP. A box because the angle changes can only be seen with angular objects. This preview showed me that the rotation still needed damping/slowing down. After some trial and error, I again decided in favour of an exponential function, which works in the same way as the depth dependency. It is of course very practical when formulae can be recycled.

angle *= exp(-@TimeInc * 0.5);

I then inserted this line of code between the angle definition float angle and the quaternion calculation vector 4 q.