Group points inside camera view in Houdini with VEX (camera culling)
TL;DR – Scroll down to the end of this page to find a VEX script made for culling points outside of camera view. In this post I’m just explaining how the script is made.
Camera Culling
There are many ways to do camera culling in Houdini but most of the setups I’ve seen are based around creating a pyramid shape from the camera view and grouping points based on this geometry and deleting the rest. I wanted to try out a method of creating this group with VEX script in a single point wrangle node while pushing my skills in mathematics and programming.
The problem
The goal of this setup is to group points that are inside the camera view. Somehow I needed to iterate over the points in the scene and find a way to check if a point is inside the camera view. I was thinking about the approach for a couple of days and eventually I chose the following path:
Measure depth along the optical axis to a plane that’s in the same depth than the iterated point. Call a position on the plane in the center of the camera view \(pc\).
Calculate a vector \( \vec{r} \) pointing to the iterated point position from \(pc\).
Calculate the horizontal and vertical components of \( \vec{r} \) relative to the camera orientation.
Calculate maximum allowed horizontal and vertical positions on the camera plane relative to the camera orientation.
Measure distance of \( \vec{r} \) horizontal component to \(pc\) and compare it with the distance of the maximum allowed horizontal position to \(pc\). Do the same for the vertical components.
If both of the distances are less than the allowed distances, the iterated point is inside the camera view.
I’m writing this post because when I try to explain my setups it helps me understand fundamentally what I’ve done. I only have high-school mathematics background so bear with me. I’m trying to explain this in a way that would be easy to understand, but I won’t be going through the simplification process of some of the mathematical formulas.
Great free resources for maths used in this post:
Solving the steps
1.
Measure depth along the optical axis to a plane that’s in the same depth than the iterated point. Call a position on the plane in the center of the camera view \(pc\).
To find the position of \(pc\) we need to know the optical axis of the camera. Then find a position on this axis that is at the depth of the iterated point.
The optical axis of the camera is the inverted z local axis of the camera transform matrix. The local z axis consists of the first three horizontal components on the third row of the matrix.
\(\require{enclose}xform = \begin{pmatrix}1 & 0 & 0 & 0\\\ 0 & 1 & 0 & 0 \\\ \enclose{circle}{0} & \enclose{circle}{0} & \enclose{circle}{1} & 0 \\\ 0 & 0 & 0 & 1 \end{pmatrix}\)
We also need the local x and y axes from this transform matrix which are found on the first and second row. These are used later on.
The formula for the position of \(pc\) is:
\(pc = cameraPosition + opticalAxis\vert \vec{d} \vert\cos(\theta)\)
Where:
\(\vec{d}\) = vector pointing to the iterated point position from the camera position. \(\vec{d} = P-cameraPosition\)
\(\theta\) = angle between optical axis and \(\vec{d}\)
The angle \(\theta\) can be calculated with this formula:
\(\theta = cos^{-1}\Biggl({{opticalAxis \bullet \vec{d}} \over {\vert opticalAxis \vert\vert \vec{d} \vert}}\Biggr)\)
So including this in the equation:
\(pc = cameraPosition + opticalAxis\vert \vec{d} \vert\cos\biggl(cos^{-1}\Biggl({{opticalAxis \bullet \vec{d}} \over {\vert opticalAxis \vert\vert \vec{d} \vert}}\Biggr)\biggr)\)
Since \(\cos(\cos^{-1}(x)) = x\), the cosine functions cancel each other out. We’re left with:
\(pc = cameraPosition + opticalAxis\vert \vec{d} \vert\Biggl({{opticalAxis \bullet \vec{d}} \over {\vert opticalAxis \vert\vert \vec{d} \vert}}\Biggr)\)
2.
Calculate a vector \( \vec{r} \) pointing to the iterated point position from \(pc\).
Now that we know where the position \(pc\) is, the vector \(\vec{r}\) pointing to the iterated point position is formed by simply subtracting the tail point from the head point. So:
\(\vec{r} = P-pc\)
3.
Calculate the horizontal and vertical components of \( \vec{r} \) relative to the camera orientation.
For this part we need the camera local x and y axes.
The camera local x axis consists of the first three horizontal components on the first row of the camera’s transform matrix.
\(xform = \begin{pmatrix}\enclose{circle}{1} & \enclose{circle}{0} & \enclose{circle}{0} & 0\\\ 0 & 1 & 0 & 0 \\\ 0 & 0 & 1 & 0 \\\ 0 & 0 & 0 & 1 \end{pmatrix}\)
And the local y axis of the first three horizontal components on the second row.
\(xform = \begin{pmatrix}1 & 0 & 0 & 0\\\ \enclose{circle}{0} & \enclose{circle}{1} & \enclose{circle}{0} & 0 \\\ 0 & 0 & 1 & 0 \\\ 0 & 0 & 0 & 1 \end{pmatrix}\)
The component of the vector \(\vec{r}\) that’s aligned with the camera local y axis can be calculated:
\(yComponent = pc + (\vert\vec{r}\vert\cos(\varphi))\cdot localYaxis\)
The angle \(\varphi\) is:
\(\varphi = \cos^{-1}\biggl({\vec{r}\bullet localYaxis \over \vert\vec{r}\vert \vert localYaxis \vert}\biggr)\)
Including this in the equation:
\(yComponent = pc + (\vert\vec{r}\vert\cos(\cos^{-1}\biggl({\vec{r}\bullet localYaxis \over \vert\vec{r}\vert \vert localYaxis \vert}\biggr)))\cdot localYaxis\)
Since \(\cos(\cos^{-1}(x))=x\) and \({\vert\vec{r}\vert \over \vert\vec{r}\vert}\) cancel each other out we’re left with:
$$ yComponent = pc + \biggr( {(\vec{r}\bullet localYaxis)\cdot localYaxis \over \vert localYaxis \vert}\biggl) $$
That’s the formula I used for the y component, and for the x component I used similar approach with some differences:
\(xComponent = pc + (\vert\vec{r}\vert\sin(\varphi))\cdot localXaxis\)
This time I used different formula to define \(\varphi\) to cancel out the sine functions.
\(\varphi = \sin^{-1}\biggr({\vert\vec{r}\times localYaxis\vert \over \vert\vec{r}\vert\vert localYaxis \vert}\biggl) \)
So for the x component we’re left with:
$$ xComponent = pc + \biggr( { \vert \vec{r}\times localYaxis \vert \cdot localXaxis \over \vert localYaxis \vert}\biggl) $$
4.
Calculate maximum allowed horizontal and vertical positions on the camera plane relative to the camera orientation.
The maximum allowed position along the camera local x axis can be calculated with this formula:
\(Xmax = pc + (\vert pc-cameraPosition\vert)\cdot\tan(\frac{angleOfView}{2})\cdot localXaxis\)
For the maximum position along the camera local y axis the formula is similar, but the result needs to be scaled down to fit the camera aspect ratio.
\(Ymax = pc + (\vert pc-cameraPosition\vert)\cdot\tan(\frac{angleOfView}{2})\cdot localYaxis \cdot (\frac{verticalRes}{horizontalRes})\)
I’ve made a separate post on the formula to calculate the angle of view:
5.
Measure distance of \( \vec{r} \) horizontal component to \(pc\) and compare it with the distance of the maximum allowed horizontal position to \(pc\). Do the same for the vertical components.
For this part I just used the distance() VEX function to find the distances.
6.
If both of the distances are less than the allowed distances, the iterated point is inside the camera view.
In VEX code:
if ( i_dist_y < max_dist_y && i_dist_x < max_dist_x ) {
i@group_inside = 1;
}Final Script
Here’s the final script. Copy and paste the code to a point wrangle. Wire in the geometry or point cloud that you want to process. Click the “Create Spare Parameters“ button next to the snippet to create the id box and use the box to select your camera node. For example: /obj/cam1
//path to camera
string op_path = chsop("id");
//initialize variables for evaluating parameters
int op_id;
int parm_index;
int vector_index;
//get transform matrix and axes from camera
matrix xform = optransform(op_path);
float a = getcomp(xform,2,0);
float b = getcomp(xform,2,1);
float c = getcomp(xform,2,2);
vector o_axis = set(a,b,c)*-1; //camera optical axis
a = getcomp(xform,1,0);
b = getcomp(xform,1,1);
c = getcomp(xform,1,2);
vector c_y = set(a,b,c); //camera local y axis
a = getcomp(xform,0,0);
b = getcomp(xform,0,1);
c = getcomp(xform,0,2);
vector c_x = set(a,b,c); //camera local x axis
float x = getcomp(xform,3,0);
float y = getcomp(xform,3,1);
float z = getcomp(xform,3,2);
vector c_pos = set(x,y,z); //camera pos
chid(op_path, "aperture", op_id, parm_index, vector_index);
float sensor_width = ch(op_id, parm_index, vector_index); //camera aperture
chid(op_path, "focal", op_id, parm_index, vector_index);
float focal_length = ch(op_id, parm_index, vector_index); //camera focal length
float aov = 2*atan( sensor_width / ( 2*focal_length ) ); //angle of view
chid(op_path, "resy", op_id, parm_index, vector_index);
float res_y = ch(op_id, parm_index, vector_index); //camera y resolution
chid(op_path, "resx", op_id, parm_index, vector_index);
float res_x = ch(op_id, parm_index, vector_index); //camera x resolution
//iterated point variables
vector i_pos = v@P; //iterated point position
vector d = i_pos-c_pos; //direction from camera to iterated point
//camera projection plane at iterated point depth
vector pc = c_pos + (o_axis*(length(d)*(dot(o_axis,d)) / (length(o_axis)*length(d)))); //center of camera plane at depth of iterated point
vector r = i_pos-pc; //vector pointing to the iterated point from the plane center
//x and y components of the vector pointing to the iterated point position from plane center
//parallel to the camera local x and y axes
vector i_pos_y = pc + ( dot(r,c_y)*c_y/ length(c_y) );
vector i_pos_x = pc + ( length(cross(r,c_y))*c_x / length(c_y) );
//maximum x and y positions on the camera plane at iterated point depth
vector max_pos_y = pc + ( distance(c_pos,pc)*(tan(aov/2)*(c_y)) * (res_y/res_x));
vector max_pos_x = pc + ( distance(c_pos,pc)*(tan(aov/2)*(c_x)));
//distance from the camera projection plane center to the iterated point positions x and y components
float i_dist_y = distance(pc,i_pos_y);
float i_dist_x = distance(pc,i_pos_x);
float pad = chf("padding");
float max_dist_y = distance(pc,max_pos_y)+pad;
float max_dist_x = distance(pc,max_pos_x)+pad;
if ( i_dist_y < max_dist_y && i_dist_x < max_dist_x ) {
i@group_inside = 1;
}