Processing math: 0%
\newcommand{\n}{\hat{n}}\newcommand{\thetai}{\theta_\mathrm{i}}\newcommand{\thetao}{\theta_\mathrm{o}}\newcommand{\d}[1]{\mathrm{d}#1}\newcommand{\w}{\hat{\omega}}\newcommand{\wi}{\w_\mathrm{i}}\newcommand{\wo}{\w_\mathrm{o}}\newcommand{\wh}{\w_\mathrm{h}}\newcommand{\Li}{L_\mathrm{i}}\newcommand{\Lo}{L_\mathrm{o}}\newcommand{\Le}{L_\mathrm{e}}\newcommand{\Lr}{L_\mathrm{r}}\newcommand{\Lt}{L_\mathrm{t}}\newcommand{\O}{\mathrm{O}}\newcommand{\degrees}{{^{\large\circ}}}\newcommand{\T}{\mathsf{T}}\newcommand{\mathset}[1]{\mathbb{#1}}\newcommand{\Real}{\mathset{R}}\newcommand{\Integer}{\mathset{Z}}\newcommand{\Boolean}{\mathset{B}}\newcommand{\Complex}{\mathset{C}}\newcommand{\un}[1]{\,\mathrm{#1}}

VEX treats positions and vectors as row vectors

VEX treats positions and vectors as row vectors
Intro · How HOM does it · How numpy does it · How VEX does it · Conclusion · Future work
2022-11-03

   

Intro

You can do vector and matrix multiplications in various places in Houdini: HDK, HOM and VEX. HDK and HOM docs state that vectors, positions are treated as row vectors during matrix multiplication. In case of VEX, I couldn't really find this information in the docs and VEX's behavior when testing it seemed odd. That motivated me to investigate it. I was wondering - does VEX treat vectors as column vectors?

Usually I just find the convention for a given library or an application and follow it: column vector (4 \times 1 matrix) is post-multiplied (vector is on the right), row vector (1 \times 4 matrix) is pre-multiplied (vector is on the left). It's not possible to use the wrong order, because then the two matrices would not have compatible shapes. VEX takes another approach though.

Here's the starting point that will be transformed in the examples below.

Starting point, cube at origin, top view

   

How HOM does it

Create a Python SOP with the following snippet. I will go over it below.

import hou node = hou.pwd() geo = node.geometry() # Initialize T, R matrices T = hou.hmath.buildTranslate(1, 0, 0) # 4x4, move 1 unit along +X axis R = hou.hmath.buildRotate(0, 45, 0) # 4x4, rotate 45 degrees along Y (up) axis for point in geo.points(): # Iterate over all points P = point.position() # Get current position new_P = P #new_P = P * T # P is treated as row vector #new_P = T * P # This does not work - as expected #new_P = P * T * R # Translate, rotate #new_P = R * T * P # This does not work - as expected point.setPosition(new_P) # Set new position

At first we import hou module, store reference to processed geometry and create T, R matrices for translation and rotation. See hmath docs, T, R are of type hou.Matrix4 and are in row-major form.

import hou node = hou.pwd() geo = node.geometry() # Initialize T, R matrices T = hou.hmath.buildTranslate(1, 0, 0) # 4x4, move 1 unit along +X axis R = hou.hmath.buildRotate(0, 45, 0) # 4x4, rotate 45 degrees along Y (up) axis

Then we iterate over each point, modifying its position.

for point in geo.points(): # Iterate over all points P = point.position() # Get current position new_P = P # ... point.setPosition(new_P) # Set new position

Now take a look at the following 4 lines. Try uncommenting each of them individually to see results or errors.

#new_P = P * T # P is treated as row vector #new_P = T * P # This does not work - as expected #new_P = P * T * R # Translate, rotate #new_P = R * T * P # This does not work - as expected

We know from the documentation that vectors are treated as row vectors, so this behaves as expected: new_P = P * T is the correct way to translate P by T. new_P = T * P does not work because P is row vector, not column one.

P * T

To combine translation, rotation transformations, we can do the following: new_P = P * T * R. new_P = R * T * P does not work again, nothing new here.

P * T * R

What I've done might seem simple and as expected. But hold on, it will be contrasted to VEX later in this post.

Note some details about the type of P, new_P variables. They are actually of hou.Vector3 type. When we multiply P * T, Python calls P.__mul__(T), which treats P automatically as position, or a hou.Vector4 with 4th component set to 1. This implementation detail allows us to multiply hou.Vector3 with hou.Matrix4.

I generally like this approach, although the implicit treatment of hou.Vector3s as positions is not ideal in my opinion.

   

How numpy does it

This is how I'd do the same transformations, but with numpy.

Create a Python SOP again with the following snippet. I will go over it too.

import hou import numpy as np node = hou.pwd() geo = node.geometry() # Initialize T_hou, R_hou matrices with the help from hou.hmath T_hou = hou.hmath.buildTranslate(1, 0, 0) # 4x4, move 1 unit along +X axis R_hou = hou.hmath.buildRotate(0, 45, 0) # 4x4, rotate 45 degrees along Y (up) axis # Create numpy arrays from our hou matrices # We know that T, R are row-major matrices T = np.array(T_hou.asTupleOfTuples()) R = np.array(R_hou.asTupleOfTuples()) for point in geo.points(): # Iterate over all points P = np.array(point.position()) # Get current position, store it in numpy array P = np.append(P, 1.0) # Append 1.0 as the 4th element, effectively treating it as a position P = P.reshape(1, 4) # Reshape it into a row vector, or 1x4 matrix new_P = P #new_P = P @ T # P is treated as row vector #new_P = T @ P # This does not work - as expected #new_P = T.T @ P.T # But we are free to treat P as a column vector and T as a column major matrix with transpose #new_P = P @ T @ R # Translate, rotate #new_P = R @ T @ P # Again, this does not work - as expected #new_P = R.T @ T.T @ P.T # But with transposes we can do it without issues point.setPosition(new_P.reshape(4)) # Set new position from numpy array

You should see similarities with the HOM example above. There is however some boilerplate code for converting between hou and numpy types. I build numpy matrices T, R with the help of hou.hmath functions, which create T_hou, R_hou matrices. When iterating over points, I reshape point's position into a row vector and later set it back to the iterated point.

import hou import numpy as np node = hou.pwd() geo = node.geometry() # Initialize T_hou, R_hou matrices with the help from hou.hmath T_hou = hou.hmath.buildTranslate(1, 0, 0) # 4x4, move 1 unit along +X axis R_hou = hou.hmath.buildRotate(0, 45, 0) # 4x4, rotate 45 degrees along Y (up) axis # Create numpy arrays from our hou matrices # We know that T, R are row-major matrices T = np.array(T_hou.asTupleOfTuples()) R = np.array(R_hou.asTupleOfTuples()) for point in geo.points(): # Iterate over all points P = np.array(point.position()) # Get current position, store it in numpy array P = np.append(P, 1.0) # Append 1.0 as the 4th element, effectively treating it as a position P = P.reshape(1, 4) # Reshape it into a row vector, or 1x4 matrix new_P = P # ... point.setPosition(new_P.reshape(4)) # Set new position from numpy array

Now let's take a look at the following 6 lines of code. Numpy uses the @ operator for matrix multiplication.

#new_P = P @ T # P is treated as row vector #new_P = T @ P # This does not work - as expected #new_P = T.T @ P.T # But we are free to treat P as a column vector and T as a column major matrix with transpose #new_P = P @ T @ R # Translate, rotate #new_P = R @ T @ P # Again, this does not work - as expected #new_P = R.T @ T.T @ P.T # But with transposes we can do it without issues

Since we have our position as a row vector and translation as row-major matrix, we can simply do this multiplication: P @ T. If we feel like using column-major form, then just transpose everything and swap multiplication order, resulting in the same transformation: T.T @ P.T

Trying T @ P does not make much sense and results in an error.

P @ T, or T.T @ P.T

To combine translation and rotation, we can either do P @ T @ R, or R.T @ T.T @ P.T. Both result in the following transformation.

P @ T @ R, or R.T @ T.T @ P.T

I like this approach the most because it's explicit and doesn't make any assumptions. Numpy just multiplies matrices and leaves it on us to decide what they represent. Of course as long as we have compatible matrix shapes. We can treat our vectors and matrices in a row-major or column-major forms, with identical results and without much hassle.

To better visualize what's going on, I've written the actual values of our matrices here. P is the currently processed point, P' is the transformed point. T, R are translation, rotation matrices.

Row-major form:

P' = P T R

P' = \begin{bmatrix} x & y & z & 1 \\ \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} 0.71 & 0 & -0.71 & 0 \\ 0 & 1 & 0 & 0 \\ 0.71 & 0 & 0.71 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}

Column-major form, R^T is the transpose of R and so on:

P^{'T} = R^T T^T P^T

P^{'T} = \begin{bmatrix} 0.71 & 0 & 0.71 & 0 \\ 0 & 1 & 0 & 0 \\ -0.71 & 0 & 0.71 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} x \\ y \\ z \\ 1 \\ \end{bmatrix}

   

How VEX does it

Having seen HOM, numpy examples above, let's have a look at VEX.

The whole snippet for Attribute Wrangle SOP is here, I will dissect it below.

// Initialize T, R matrices matrix T = ident(); // 4x4 matrix R = ident(); // 4x4 // Mutate T, R with respective transforms translate(T, {1, 0, 0}); // Move 1 unit along +X axis rotate(R, radians(45), {0,1,0}); // Rotate 45 degrees along Y (up) axis // Both cases below work in VEX, producing identical result //v@P = T * v@P; // P is treated as a column vector? (or 4x1 matrix), is post-multiplied //v@P = v@P * T; // P is treated as a row vector? (or 1x4 matrix), is pre-multiplied //v@P = T * v@P * R; // Hmm, interesting //v@P = (T * v@P) * R; // The same result as above: translate, then rotate //v@P = T * (v@P * R); // A different result: rotate, then translate //v@P = v@P * T * R; // Result: translate, rotate //v@P = R * T * v@P; // Result: rotate, translate

This section creates 2 matrices stored in variables T, R for translation and rotation.

// Initialize T, R matrices matrix T = ident(); // 4x4 matrix R = ident(); // 4x4 // Mutate T, R with respective transforms translate(T, {1, 0, 0}); // Move 1 unit along +X axis rotate(R, radians(45), {0,1,0}); // Rotate 45 degrees along Y (up) axis

Now comes the confusing part. Uncomment either of the two lines below and you will see the same result.

// Both cases below work in VEX, producing identical result //v@P = T * v@P; // P is treated as a column vector? (or 4x1 matrix), is post-multiplied //v@P = v@P * T; // P is treated as a row vector? (or 1x4 matrix), is pre-multiplied

Translated cube, T * v@P, or v@P * T

I'd expect one line to result in an error. This would reveal how's vectors treated in VEX. But it didn't.

Then I've tried this: v@P should now decide whether it's a row or a column and one multiplication should fail. Again, try uncommenting the lines below (only one at the time).

//v@P = T * v@P * R; // Hmm, interesting //v@P = (T * v@P) * R; // The same result as above: translate, then rotate //v@P = T * (v@P * R); // A different result: rotate, then translate

None of the lines above resulted in an error.

Result of T * v@P * R and (T * v@P) * R

T * (v@P * R)

Now it looks like in T * v@P * R, position gets first translated (as a column vector) and then rotated (as a row vector). Still unsure what's going on here. Let's try the following:

//v@P = v@P * T * R; // Result: translate, rotate //v@P = R * T * v@P; // Result: rotate, translate

v@P * T * R, translate, rotate

R * T * v@P, rotate, translate

Based on this it seems like VEX's evaluation rules are like this:

   

Conclusion

So perhaps unsurprisingly it seems that VEX treats positions, vectors as row vectors - consistently with HDK, HOM. But it allows for multiplication in cases where it should not be possible. I find doing so really confusing and wouldn't recommend using such style.

I suspect SideFX has made this feature to automatically handle common user errors, but in my opinion at the cost of clarity and transparency of how things work internally.

And the big question again: does VEX treat vectors as column vectors?

Kind of

That's it for this post :)

   

Future work

A similar confusion is in the fact that both v@P, v@N attributes are treated as 3-vectors in both VEX and Geometry Spreadsheet. Thus hiding the differentiation between the two, which would be in the 4th element in a 4-vector type. Perhaps Houdini differentiates between the two cases with a different mechanism, e.g. to save memory? But I am digressing now and that would be for another post anyway.