I’ve been working on a Swift library that implements Lindenmayer systems, but the past week has me deeply stuck on a specific 3D math problem. There’s a specific 3D rendering command that’s described in The Algorithmic Beauty of Plants – the ‘$’ character in that work, that has a special meaning that’s taken me quite a bit to sort out. The idea is that as you progress through rendering, this particular command rotates around the axis that you’re “growing in” such that the “up” vector (which is 90° pitched up from the “forward” vector) is as close to vertical as possible.
It turns out this little critter is critical to getting tree representations that match the examples of what’s in the book. I haven’t solved it – I thought I had earlier, but I managed to delude myself, so now I’m back to trying to sort the problem. The course of solving this had let to some nice side effects – such I know have a nice 3D debugging view for my rendered trees – but I haven’t yet finalized on HOW to solve the problem. In years past, if I saw someone stuck this badly on a problem, I’d advise them to ask for help – so that’s what I’m doing.
I’m using a 3D affine transform (a 4×4 matrix of floats) to represent a combination of translations and rotations in 3D space doing a sort of 3D turtle graphics kind of thing. From this state I know the heading that I’m going, and what I’d like to do is roll around this particular axis. The problem that I’m trying to solve is determining the angle (let’s call it θ) to roll that results in one of the heading vectors being as close to vertical (+Y) as possible while still on the plane of rotation that’s constrained to the “heading” axis.
My starting points for the “forward” heading is +1 on the Y axis, and a local “up” heading as +1 on the Z axis.
The way I was trying to tackle this problem was applying the built-up set of translations and rotations using the 3D Affine transform and then figuring out if there was trigonometry I could use to solve for the angle. Since I know the axis around which I want to roll (or can compute it by applying the transform to the unit vector that represents it – the (0,0,1) vector), I was looking at the Rodrigues rotational formula, but my linear algebra (matrix math) skills and understanding are fairly weak – and I couldn’t see a path to solving that equation for θ given known vectors for the heading, or vectors on the plane that can be used to define a cross-product that is the heading, as well as knowing the world vector.
I’m heading towards trying to solve this by iteratively applying various values of θ and homing in on the solution based on the the resulting value that has the Y component value. I can apply the roll as an affine transform that I multiply onto the current transform, and then test the result of a unit “up” vector – rinse and repeat to find the one that gives me the best “Y” component value.
I’d like to know If there’s a way to solve this directly – to compute the value of θ that I can use to directly do the rotation, rather than numerically iterate/solve into the solution. I wasn’t sure how active (or if I’d get a response) on the GameDev stack exchange but I tried asking:
If you’re familiar with 3D graphics, rotations, transforms, or the mathematics of solving this sort of thing, I’d love to have any input or advice on how to solve this problem, even just some concrete knowledge of if this problem is amenable to a direct solution – or if it’s the kind of thing that requires an iterative numerical solution like I’m considering.
UPDATE: Solved!
I got an answer from a friend in slack who saw these, and I’ve mostly implemented it. There’s a few numerical instability points I need to sort out, but the gist is: Yes, there’s definitely a way to explicitly calculate the angle of rotation needed. The summary of the answer that put me onto the right path is: “Look at this problem as projecting the vector you want to roll towards as a vector on the plane that’s the base of the heading. Once its projected, you can compute the angle between two vectors, and that should be what you need.”
The flow of the process is as follows:
- start with a unit vector that represents the ‘up’ vector that I want to compare (+Z in my case)
- pull the 3×3 translation matrix from the affine matrix, and use that to rotate the ‘up’ vector. We’ll be comparing against this later to get the angle. Because it’s explicitly set as a unit vector in the ‘up’ direction, we already know it’s on the plane that can be represented by the normal of that plane – which is our heading vector.
- Use a similar technique to rotate the ‘heading’ vector (what starts as +Y for me) by the rotation matrix.
(a quick double check here I did in my testing was that these two remained 90° apart after rotation – primarily to verify that I didn’t screw up the rotation calculation)
- Now that we have the heading, we use that as a normal vector for the plane upon which we want to project our two vectors, and from which we can get the angle desired. The vector (the ‘rotated up’ vector) is already on this plane. The other vector is the +Y world vector – the “as vertical as possible” component of this.
The formula for projecting a vector on a plane is:
vec_projected = vector - ( ( vector • plane_normal ) / plane_normal.length^2 ) * plane_normal
You can look at this conceptually as taking the vector you want to project and subtracting from it the portion of the vector that corresponds to the normal vector, which leaves you with just the component that’s aligned on the plane.
🎩-tip to Greg Titus, who referred me to ProjectionOfVectorOntoPlan at MapleSoft.
- Once both the vectors are projected onto that base plane, then use the equation to calculate the angle between two vectors:
acos(
dot(rotated_up_vector, projected_vector) /
( length(projected_vector) * length(rotated_up_vector)
)
)
There’s some numerical instability points I need to work out where my current code is returning ‘NaN’ from the final comparison, and the directional component isn’t included in that – so I need to sort out a way to determine which direction to rotate, but the fundamentals part of it is at least sorted.
Solved, even better
After I’d implemented most of the above details to prove to myself that it worked, I received a suggestion from DMGregory with a significantly better answer. His solution uses Vector cross-products to get define the plane that’s the base upon which we’ll rotate, and then uses the Arctangent function with two parameters to compute the angle from the dot products from those two angles.
It’s a denser solution, but he provided a helpful way to think about it that made a lot of sense to me after I sketched it all out in a notebook so I could understand it:
planeRight = normalize(cross(worldUp, heading));
planeUp = cross(heading, planeRight);
angle = atan2(dot(currentUp, planeRight), dot(currentUp, planeUp));
You can think of the dot products as getting us the
DMGregoryx&ycoordinates of our current up vector in this plane, and from that the two-argument arctangent gets us the angle of the vector from the positive x-axis in that plane.
Since I’m using SceneKit, I did have to fiddle around with the cross-products to make them legit for a right-handed coordinate scheme, but the resulting solution helpfully provides a proper direction to rotation as well as the value, which is something that the original solution didn’t have.
Rhonabwy 