Geometry

Converting between quaternions and rotation matrices

Posted on by John

In the previous post I wrote about representing rotations with quaternions. This representation has several advantages, such as making it clear how rotations compose. Rotations are often represented as matrices, and so it’s useful to be able to go between the two representations.

A unit-length quaternion (q0, q1, q2, q3) represents a rotation by an angle θ around an axis in the direction of (q1, q2, q3) where cos(θ/2) = q0. The corresponding rotation matrix is given below.

R = \begin{pmatrix} 2(q_0^2 + q_1^2) - 1 & 2(q_1 q_2 - q_0 q_3) & 2(q_1 q_3 + q_0 q_2) \\ 2(q_1 q_2 + q_0 q_3) & 2(q_0^2 + q_2^2) - 1 & 2(q_1 q_3 - q_0 q_1) \\ 2(q_1 q_3 - q_0 q_2) & 2(q_2 q_3 + q_0 q_1) & 2(q_0^2 + q_3^2) - 1 \end{pmatrix}

Going the other way around, inferring a quaternion representation from a rotation matrix, is harder. Here is a mathematically correct but numerically suboptimal method known [1] as the Chiaverini-Siciliano method.

\begin{align*} q_0 &= \frac{1}{2} \sqrt{1 + r_{11} + r_{22} + r_{33}} \\ q_1 &= \frac{1}{2} \sqrt{1 + r_{11} - r_{22} - r_{33}} \text{ sgn}(r_{32} - r_{32}) \\ q_2 &= \frac{1}{2} \sqrt{1 - r_{11} + r_{22} - r_{33}} \text{ sgn}(r_{13} - r_{31}) \\ q_3 &= \frac{1}{2} \sqrt{1 - r_{11} - r_{22} + r_{33}} \text{ sgn}(r_{21} - r_{12}) \end{align*}

Here sgn is the sign function; sgn(x) equals 1 if x is positive and −1 if x is negative. Note that the components only depend on the diagonal of the rotation matrix, aside from the sign terms. Better numerical algorithms make more use of the off-diagonal elements.

Accounting for degrees of freedom

Something seems a little suspicious here. Quaternions contain four real numbers, and 3 by 3 matrices contain nine. How can four numbers determine nine numbers? And going the other way, out of the nine, we essentially choose three that determine the four components of a quaternion.

Quaterions have four degrees of freedom, but we’re using unit quaternions, so there are basically three degrees of freedom. Likewise orthogonal matrices have three degrees of freedom. An axis of rotation is a point on a sphere, so that has two degrees of freedom, and the degree of rotation is the third degree of freedom.

In topological terms, the unit quaternions and the set of 3 by 3 orthogonal matrices are both three dimensional manifolds, and the former is a double cover of the latter. It is a double cover because a unit quaternion q corresponds to the same rotation as −q.

Python code

Implementing the equations above is straightforward.

import numpy as np

def quaternion_to_rotation_matrix(q):
    q0, q1, q2, q3 = q
    return np.array([
        [2*(q0**2 + q1**2) - 1, 2*(q1*q2 - q0*q3), 2*(q1*q3 + q0*q2)],
        [2*(q1*q2 + q0*q3), 2*(q0**2 + q2**2) - 1, 2*(q2*q3 - q0*q1)],
        [2*(q1*q3 - q0*q2), 2*(q2*q3 + q0*q1), 2*(q0**2 + q3**2) - 1]
    ]) 

def rotation_matrix_to_quaternion(R):
    r11, r12, r13 = R[0, 0], R[0, 1], R[0, 2]
    r21, r22, r23 = R[1, 0], R[1, 1], R[1, 2]
    r31, r32, r33 = R[2, 0], R[2, 1], R[2, 2]
    
    # Calculate quaternion components
    q0 = 0.5 * np.sqrt(1 + r11 + r22 + r33)
    q1 = 0.5 * np.sqrt(1 + r11 - r22 - r33) * np.sign(r32 - r23)
    q2 = 0.5 * np.sqrt(1 - r11 + r22 - r33) * np.sign(r13 - r31)
    q3 = 0.5 * np.sqrt(1 - r11 - r22 + r33) * np.sign(r21 - r12)
    
    return np.array([q0, q1, q2, q3])

Random testing

We’d like to test the code above by generating random quaternions, converting the quaternions to rotation matrices, then back to quaternions to verify that the round trip puts us back essentially where we started. Then we’d like to go the other way around, starting with randomly generated rotation matrices.

To generate a random unit quaternion, we generate a vector of four independent normal random values, then normalize by dividing by its length. (See this recent post.)

To generate a random rotation matrix, we use a generator that is part of SciPy.

Here’s the test code:

def randomq():
    q = norm.rvs(size=4)
    return q/np.linalg.norm(q)

def randomR():
    return special_ortho_group.rvs(dim=3)

np.random.seed(20250507)
N = 10

for _ in range(N):
    q = randomq()
    R = quaternion_to_rotation_matrix(q)
    t = rotation_matrix_to_quaternion(R)
    print(np.linalg.norm(q - t))
    
for _ in range(N):
    R = randomR()
    q = rotation_matrix_to_quaternion(R)
    T = quaternion_to_rotation_matrix(q)
    print(np.linalg.norm(R - T))

The first test utterly fails, returning six 2s, i.e. the round trip vector is as far as possible from the vector we started with. How could that happen? It must be returning the negative of the original vector. Now go back to the discussion above about double covers: q and −q correspond to the same rotation.

If we go back and add the line

    q *= np.sign(q[0])

then we standardize our random vectors to have a positive first component, just like the vectors returned by rotation_matrix_to_quaternion.

Now our tests all return norms on the order of 10−16 to 10−14. There’s a little room to improve the accuracy, but the results are good.

Update: I did some more random testing, and found errors on the order of 10−10. Then I was able to create a test case where rotation_matrix_to_quaternion threw an exception because one of the square roots had a negative argument. In [1] the authors get around this problem by evaluating two theoretically equivalent expressions for each of the square root arguments. The expressions are complementary in the sense that both should not lead to numerical difficulties at the same time.

[1] See “Accurate Computation of Quaternions from Rotation Matrices” by Soheil Sarabandi and Federico Thomas for a better numerical algorithm. See also the article “A Survey on the Computation of Quaternions From Rotation Matrices” by the same authors.

Composing rotations using quaternions

Posted on by John

Every rotation in 3D fixes an axis [1]. This is Euler’s rotation theorem from 1775. Another way to state the theorem is that no matter how you rotate a sphere about its center, two points stay fixed.

The composition of two rotations is a rotation. So the first rotation fixes some axis, the second rotation fixes some other axis, and the composition fixes some third axis. It’s easy to see what these axes are if we work with quaternions. (Quaternions were discovered decades after Euler proved his rotation theorem.)

A rotation by θ about the axis given by a unit vector v = (v1, v2, v3) corresponds to the quaternion

q = (cos(θ/2), sin(θ/2)v1, sin(θ/2)v2, sin(θ/2)v3).

To rotate a point p = (p1, p2, p3) by an angle θ about the axis v, first embed p as a quaternion by setting its first coordinate to 0:

p → (0, p1, p2, p3)

and multiply the quaternion p on the left by q and on the right by the conjugate of q, written q*. That is, the rotation takes p to

p′ = qpq*.

This gives us a quaternion p′, not a 3D vector. We recover the vector by undoing the embedding, i.e. chopping off the first coordinate.

Making things clearer

Since q has unit length, the conjugate of q is also its inverse: q* = q−1. Usually rotations are described as above: multiply on the left by q and on the right by q*. In my opinion it’s clearer to say

p′ = qpq−1.

Presumably sources say q* instead of q−1 because it’s obvious how to compute q* from q; it’s not quite as obvious that this also gives the inverse of q.

Another thing about the presentation above that, while standard, could be made clearer is the role Euclidean space and quaternions. It’s common to speak of the real and quaternion representations of a vector, but we could make this more explicit by framing this as an embedding E from ℝ³ to the quaternions ℍ and a projection P from ℍ back to ℝ³ [3].

\[\begin{tikzcd} {\mathbb{H}} && {\mathbb{H}} \\ \\ {\mathbb{R}^3} && {\mathbb{R}^3} \arrow[

The commutative diagram says we end up with the same result regardless of which of two paths we take: we can do the rotation directly from ℝ³ to ℝ³, or we could project into ℍ, multiply by q on the left and divide by q on the right, and project back down to ℝ³.

Composition

Composing rotations represented by quaternions is simple. Rotating by a quaterions q and then by a quaterion r is the same as rotating by rq. Proof:

r(qpq−1)r−1 = (rq)p(q−1r−1) = (rq)p(rq)−1.

See the next post for how to convert between quaternion and matrix representations of a rotation.

Related posts

[1] Note that this is not the case in two dimensions, nor is it true in higher even dimensions.

[2] We assumed v had unit length as a vector in ℝ³. So

||q||² = cos²(θ/2) + sin²(θ/2) = 1.

[3] Why ℍ for quaternions? ℚ for rationals was already taken, so we use ℍ for William Rowan Hamilton.

How many ways can you triangulate a regular polygon?

Posted on by John

In this post we want to count the number of ways to divide a regular polygon [1] into triangles by connecting vertices with straight lines that do not cross.

Squares

For a square, there are two possibilities: we either connect the NW and SE corners,

or we connect the SW and NE corners.

Pentagons

For a pentagon, we pick one vertex and connect it to both non-adjacent vertices.

We can do this for any vertex, so there are five possible triangulations. All five triangulations are rotations of the same triangulation. What if we consider these rotations as equivalent? We’ll get to that later.

Hexagons

For a hexagon, things are more interesting. We can again pick any vertex and connect it to all non-adjacent vertices, giving six triangulations.

But there are more possibilities. We could connect every other vertex, creating an equilateral triangle inside. We can do this two ways, connecting either the even-numbered vertices or the odd-numbered vertices. Either triangulation is a rotation of the other.

We can also connect the vertices in a zig-zag pattern, creating an N-shaped pattern inside. We could also rotate this triangulation one or two turns. (Three turns gives us the same pattern again.)

Finally, we could also connect the vertices creating a backward N pattern.

General case

So to recap, we have 2 ways to triangulate a square, 5 ways to triangulate a pentagon, and 6 + 2 + 3 + 3 = 14 ways to triangulate a hexagon. Also, there is only 1 way to triangulate a triangle: do nothing.

Let Cn be the number of ways to triangulate a regular (n + 2)-gon. Then we have C1 = 1, C2 = 2, C3 = 5, and C4 = 14.

In general,

C_n = \frac{1}{n+1}\binom{2n}{n}

which is the nth Catalan number.

Catalan numbers are the answers to a large number of questions. For example, Cn is also the number of ways to fully parenthesize a product of n + 1 terms, and the number of full binary trees with n + 1 nodes.

The Catalan numbers have been very well studied, and we know that asymptotically

C \sim \frac{4^n}{n^{3/2} \sqrt{\pi}}

so we can estimate Cn for large n. For example, we could use the formula above to estimate the number of ways to triangulate a 100-gon to be 5.84 ×1055. The 98th Catalan number is closer to 5.77 ×1055. Two takeaways: Catalan numbers grow very quickly, and we can estimate them within an order of magnitude using the asymptotic formula.

Equivalence classes

Now let’s go back and count the number of triangulations again, considering some variations on a triangulation to be the same triangulation.

We’ll consider rotations of the same triangulation to count only once. So, for example, we’ll say there is only one triangulation of a pentagon and four triangulations of a hexagon. If we consider mirror images to be the same triangulation, then there are three triangulations of a hexagon, counting the N pattern and the backward N pattern to be the same.

Grouping rotations

The number of equivalence classes of n-gon triangulations, grouping rotations together, is OEIS sequence A001683. Note that the sequence starts at 2.

OEIS gives a formula for this sequence:

a(n) = \frac{1}{2n}C_{n-2} + \frac{1}{4}C_{n/2-1} + \frac{1}{2} C_{\lceil (n+1)/2\rceil - 2} + \frac{1}{3} C_{n/3 - 1}
where Cx is zero when x is not an integer. So a(6) = 4, as expected.

Grouping rotations and reflections

The number of equivalence classes of n-gon triangulations, grouping rotations and reflections together, is OEIS sequence A000207. Note that the sequence starts at 3.

OEIS gives a formula for this sequence as well:

a(n) = \frac{1}{2n}C_{n-2} + \frac{1}{4}C_{n/2-1} + \frac{1}{2} C_{\lceil (n+1)/2\rceil - 2} + \frac{1}{3} C_{n/3 - 1}

As before, Cx is zero when x is not an integer. This gives a(6) = 3, as expected.

The formula on the OEIS page is a little confusing since it uses C(n) to denote Cn−2 .

Related posts

[1] Our polygons do not need to be regular, but they do need to be convex.

Erdős-Mordell triangle theorem

Posted on by John

If any field of mathematics has been thoroughly combed over, it’s Euclidean geometry. But once in a while someone will come up with a new theorem that seems it should have been discovered centuries ago.

Here’s a theorem conjectured by Paul Erdős in 1935 and proved by Louis Mordell later the same year.

If from a point O inside a given triangle ABC perpendiculars OD, OE, OF are drawn to its sides, then

OA + OB + OC ≥ 2(OD + OE + OF).

Equality holds if and only if triangle ABC is equilateral.

To put it more succinctly,

From any interior point, the distances to the vertices are at least twice the distances to the sides.

Here’s an illustration. In the figure above, the theorem says the dashed blue lines together are more than twice as long as the solid red lines.

In the units I used in drawing the figure above, the blue lines have combined length 9.5 and the red lines have combined length 4.7.

Hojoo Lee gave an elementary proof of the Erdős-Mordell theorem in 2001 that takes about one printed page [1].

In my opinion the Erdős-Mordell theorem feels like a theorem an ancient geometer could have discovered and proved. Here’s a generalization of the theorem that feels much more contemporary [2].

Let Ri be the distance from an interior point O to the ith vertex and ri the distance to the side opposite the ith vertex. Let λ1, λ2, and λ3 be any three positive real numbers. Then

\sum_{i=1}^3 \lambda_iR_i \geq 2\sqrt{\lambda_1\lambda_2\lambda_3} \sum_{i=1}^3 \frac{r_i}{\sqrt{\lambda_i}}

When all the λs equal 1, we get the original Erdős-Mordell theorem.

You could say the weighted distances to the vertices are at least twice the weighted distances to the sides, but you have to say more about the weights, and in general the weights work differently on both sides of the inequality.

[1] Hojoo Lee. Another Proof of the Erdős-Mordell Theorem. Forum Geometricorum, Volume 1 (2001) p. 7–8

[2] Seannie Dar, Shay Gueron. A Weighted Erdős-Mordell Inequality. The American Mathematical Monthly. Vol. 108, No. 2, Feb., 2001

Interior of a conic

Posted on by John

What is the interior of a circle? Obvious.

What is the interior of a parabola? Not quite as obvious.

What is the interior of a hyperbola? Not at all obvious.

Is it possible to define interior in a way that applies to all conic sections?

Circles

If you remove a circle from the plane, there are two components left. Which one is the interior and which one is the exterior?

Obviously the bounded part is the interior and the unbounded part is the exterior. But using boundedness as our criteria runs into immediate problems.

Parabolas

If you remove a parabola in the plane, which component is the interior and which is the exterior? You might say there is no interior because both components of the plane minus the parabola are unbounded. Still, if you had to label one of the components the interior, you’d probably say the smaller one.

But is the “smaller” component really smaller? Both components have infinite area. You could patch this up by taking a square centered at the origin and letting its size grow to infinity. The interior of the parabola is the component that has smaller area inside the square all along the way.

Hyperbolas

A hyperbola divides the plane into three regions. Which of these is the interior? If we try to look at area inside an expanding square, it’s not clear which component(s) will have more or less area. Seems like it may depend on the location of the center of the square relative to the position of the hyperbola.

Tangents to a circle

Here’s another way to define the interior of a circle. Look at the set of all lines that are tangent to a point on the circle. None of them go through the interior of the circle. We can define the interior of the circle as the set of points that no tangent line passes through.

This clearly works for a circle, and it’s almost as clear that it would work for an ellipse.

How do we define the exterior of a circle? We could just say it’s the part of the plane that isn’t the interior or the circle itself. But there is a more interesting definition. If the interior of the circle consists of points that tangent lines don’t pass through, the exterior of the circle consists of the set of points that tangent lines do pass thorough. Twice in fact: every point outside the circle is at the intersection of two lines tangent to the circle.

To put it another way, consider the set of all tangent lines to a circle. Every point in the plane is part of zero, one, or two of these lines. The interior of the circle is the set of points that belong to zero tangent lines. The circle is the set of points that belong to one tangent line. The exterior of the circle is the set of points that belong to two tangent lines.

Tangents to a parabola

If we apply the analogous definition to a parabola, the interior of the parabola works out to be the part we’d like to call the interior.

It’s not obvious that every point of the plane not on the parabola and not in the interior lies at the intersection of two tangent lines, but it’s true.

Tangents to a hyperbola

If we look at the hyperbola x² − y² = 1 and draw tangent lines, the interior, the portion of the plane with no crossing tangent lines, is the union of two components, one containing (−∞, −1) and one containing (1, ∞). The exterior is then the component containing the line y = x. In the image above, the pink and lavender components are the interior and the green component is the exterior.

It’s unsatisfying that the interior of the hyperbola is disconnected. Also, I believe the exterior is missing the origin. Both of these annoyances go away when we add points at infinity. In the projective plane, the complement of a conic section consists of two connected components, the interior and the exterior. The origin lies on two tangent lines: one connecting (−∞, −∞) to (∞, ∞) and one connecting (−∞, ∞) to (∞, −∞).

Do perimeter and area determine a triangle?

Posted on by John

Is the shape of a triangle determined by its perimeter and area? In other words, if two triangles have the same area and the same perimeter, are the triangles similar? [1]

It’s plausible. A triangle has three degrees of freedom: the lengths of the three sides. Specifying the area and perimeter removes two degrees of freedom. Allowing the triangles to be similar rather than congruent accounts for a third degree of freedom.

Here’s another plausibility argument. Heron’s formula computes the area of a triangle from the lengths of the sides.

A = \sqrt{s(s-a)(s-b)(s-c)}

Here s is the semi-perimeter, half of the sum of the lengths of the sides. So if the perimeter and area are known, we have a third order equation for the sides:

(a - s)(b - s)(c - s) = -\frac{A^2}{s}

If the right-hand side were 0, then we could solve for the lengths of the sides. But the right-hand side is not zero. Is it still possible that the sides are uniquely determined, up to rearranging how we label the sides?

It turns out the answer is no [2], and yet it is not simple to construct counterexamples. If all the sides of a triangle are rational numbers, it is possible to find a non-congruent triangle with the same perimeter and area, but the process of finding this triangle is a bit complicated.

One example is the triangles with sides (20, 21, 29) and (17, 25, 28). Both have perimeter 70 and area 210. But the former is a right triangle and the latter is not.

Where did our algebraic argument go wrong? How can a cubic equation have two sets of solutions? But we don’t have a cubic equation in one variable; we have an equation in three variables that is the product of three linear terms.

What third piece of information would specify a triangle uniquely? If you knew the perimeter, area, and the length of one side, then the triangle is determined. What if you specified the center of the triangle? There are many ways to define a center of a triangle; would some, along with perimeter and area, uniquely determine a triangle while others would not?

Related posts

[1] Two triangles are similar if you can transform one into the other by scaling and/or rotation.

[2] Mordechai Ben-Ari. Mathematical Surprises. Springer, 2022. The author sites this blog post as his source.

Area of a quadrilateral from the lengths of its sides

Posted on by John

Last week Heron’s formula came up in the post An Unexpected Triangle. Given the lengths of the sides of a triangle, there is a simple expression for the area of the triangle.

A = \sqrt{s(s-a)(s-b)(s-c)}

where the sides are a, b, and c and s is the semiperimeter, half the perimeter.

Is there an analogous formula for the area of a quadrilateral? Yes and no. If the quadrilateral is cyclic, meaning there exists a circle going through all four of its vertices, then Brahmagupta’s formula for the area of a quadrilateral is a direct generalization of Heron’s formula for the area of a triangle. If the sides of the cyclic quadrilateral are a, b, c, and d, then the area of the quadrilateral is

A = \sqrt{(s-a)(s-b)(s-c)(s-d)}

where again s is the semiperimeter.

But in general, the area of a quadrilateral is not determined by the length of its sides alone. There is a more general expression, Bretschneider’s formula, that expresses the area of a general quadrilateral in terms of the lengths of its sides and the sum of two opposite angles. (Either pair of opposite angles lead to the same value.)

A = \sqrt {(s-a)(s-b)(s-c)(s-d) - abcd \, \cos^2 \left(\frac{\alpha + \gamma}{2}\right)}

In a cyclic quadrilateral, the opposite angles α and γ add up to π, and so the cosine term drops out.

The contrast between the triangle and the quadrilateral touches on an area of math called distance geometry. At first this term may sound redundant. Isn’t geometry all about distances? Well, no. It is also about angles. Distance geometry seeks results, like Heron’s theorem, that only depend on distances.

Related posts

The mathematics of GPS

Posted on by John

The basic idea of GPS is that if you know the distance to several satellites, you can figure out your position. But you don’t actually know, or need to know, the distance to the satellites: you know the time (according to each satellite’s clock) when the signals were sent, and you know the time (according to your clock) when the signals arrived.

The atomic clocks on satellites are synchronized with each other to within a nanosecond, but they’re not synchronized with your clock. There is some offset t between your clock and the satellites’ clocks. Presumably t is small, but it matters.

If you observe m satellites, you have a system of m equations in 4 unknowns:

|| ai − x || = ti − t

where ai is the known position of the ith satellite in 3 dimensions, x is the observer’s position in three dimensions, and ti is the difference between the time when the signal left the ith satellite (according to its clock) and the time when the signal arrived (according to the observer’s clock). This assumes we choose units so that the speed of light is c = 1.

So we have a system of m equations in 4 unknowns. It’s plausible there could be a unique solution provided m = 4. However, this is not guaranteed.

Here’s an example to suggest why there may not be a unique solution. Suppose t is known to be 0. Then observing 3 satellites will give us 3 equations in 3 unknowns. Each ti determines a sphere of radius ti. Suppose two spheres intersect in a circle, and the third sphere intersects this circle in two points. This means we have two solutions to our system of equations.

In [1] the authors thoroughly study the solution to the GPS system of equations. They allow the satellites and the observer to be anywhere in space and look for conditions under which the system has a unique solution. In practice, GPS satellites are approximately confined to a sphere (more on that here) and the observer is also approximately confined to a sphere, namely the earth’s surface, but the authors do not take advantage of these assumptions.

The authors also assume the problem is formulated in n dimensional space, where n does not necessarily equal 3. It’s complicated to state when the system of equations has a unique solution, but allowing n to vary does not add to the complexity.

I’m curious whether there are practical uses for the GPS problem when n > 3. There are numerous practical problems involving the intersections of spheres in higher dimensions, where the dimensions are not Euclidean spacial dimensions but rather abstract degrees of freedom. But off hand I cannot think of a problem that would involve the time offset that GPS location finding has.

[1] Mireille Boutin, Gregor Kemperc. Global positioning: The uniqueness question and a new solution method. Advances in Applied Mathematics 160 (2024)

Triangle circle maximization problem

Posted on by John

Let a, b, and c be the sides of a triangle. Let r be the radius of an inscribed circle and R the radius of a circumscribed circle. Finally, let p be the perimeter. Then the previous post said that

2prR = abc.

We could rewrite this as

2rR = abc / (a + b + c)

The right hand side is maximized when a = b = c. To prove this, maximize abc subject to the constraint a + b + c = p using Lagrange multipliers. This says

[bc, ac, ab] = λ[1, 1, 1]

and so ab = bc = ac, and from there we conclude a = b = c. This means among triangles with any given perimeter, the product of the inner and outer radii is maximized for an equilateral triangle.

The inner radius for an equilateral triangle is (√3 / 6)a and the outer radius is a/√3, so the maximum product is a²/6.

Related posts