I'm trying to calculate the tight bounding box of a Bézier curve. The curve is in 3D space and can have any amount of points.
These articles are pretty much the only resources on the internet: https://pomax.github.io/bezierinfo/#aligning https://snoozetime.github.io/2018/05/22/bezier-curve-bounding-box.html
However they all use curves in 2D space.
Both of them firstly align the curve so that it starts at (0,0) and so that the last point lies on the x-axis. This is the code used: https://pomax.github.io/bezierinfo/chapters/tightbounds/tightbounds.js
translatePoints(points) {
// translate to (0,0)
let m = points[0];
return points.map(v => {
return {
x: v.x - m.x,
y: v.y - m.y
}
});
}
rotatePoints(points) {
// rotate so that last point is (...,0)
let last = points.length - 1;
let dx = points[last].x;
let dy = points[last].y;
let a = atan2(dy, dx);
return points.map(v => {
return {
a: a,
x: v.x * cos(-a) - v.y * sin(-a),
y: v.x * sin(-a) + v.y * cos(-a)
};
});
}
I am having trouble converting the rotation function from 2D to 3D. Could anyone help me?
Let's look at this mathematically, using basic vector algebra.
Note that this answer somewhat exceeds what OP is strictly asking for, but that is because I believe that for us implementors, it is crucial to understand what we are implementing, instead of just copying a given solution and hoping it works. Feel free to downvote if that irks you.
I use notation $\vec{a} = ( a_x, a_y, a_z)$ for vectors in general, and $\hat{b} = (b_x, b_y, b_z)$ for unit vectors. Also, $\lVert \vec{a} \rVert = \sqrt{\vec{a} \cdot \vec{a}} = \sqrt{a_x^2 + a_y^2 + a_z^2}$ is the magnitude or Euclidean length of vector $\vec{a}$. All vectors whose length is $1$ are unit vectors, so $\left\lVert \hat{b} \right\rVert = 1$ by definition.
By Bézier curve, OP and the linked articles are referring to piecewise cubic Bézier curves, where each segment ("piece") can be defined via the vector-valued function $\vec{p}(t)$, $$\vec{p}(t) = (1-t)^3 ~ \vec{P}_0 + 3 (1-t)^2 t ~ \vec{P}_1 + 3 (1-t) t^2 ~ \vec{P}_2 + t^3 ~ \vec{P}_3 \quad 0 \le t \le 1 \tag{1}\label{1}$$ where $\vec{P}_0$ is the starting point ($t = 0$) of the segment, $\vec{P}_3$ is the ending point ($t = 1$) of the segment, and $\vec{P}_1$ and $\vec{P}_2$ are the control points for the segment.
The method relies on the fact that when the curve is translated and rotated so that the starting point is at origin, and ending point is on the $x$ axis, the curve reaches its extremum (minimum or maximum) coordinate on the other axis/axes either at the starting or ending point, or when the derivative of that coordinate with respect to $t$ is zero. (This is true for any curve defined using a polynomial function of form $\eqref{1}$, which is everywhere continuous and differentiable.)
If we write $\eqref{1}$ in polynomial form, with $\vec{P}_0 = \vec{0}$ (starting point at origin), we have $$\begin{aligned} \vec{p}(t) &= \left( \vec{P}_3 - 3 \vec{P}_2 + 3 \vec{P}_1 \right) ~ t^3 \\ ~ &= \left( 3 \vec{P}_2 - 6 \vec{P}_1 \right) ~ t^2 \\ ~ &= \left( 3 \vec{P}_1 \right) ~ t \\ \end{aligned} \tag{2a}\label{2a}$$ and its derivative with respect to $t$ is $$\begin{aligned} \frac{d ~ \vec{p}(t)}{d t} &= \left( 3 \vec{P}_3 - 9 \vec{P}_2 + 9 \vec{P}_1 \right) ~ t^2 \\ ~ &= \left( 6 \vec{P}_2 - 12 \vec{P}_1 \right) ~ t \\ ~ &= \left( 3 \vec{P}_1 \right) \\ \end{aligned} \tag{2b}\label{2b}$$ Note that for Cartesian coordinates, these really are three equations each, one for each Cartesian coordinate axis ($x$, $y$, and $z$) in 3D.
Solving the roots of $\eqref{2b}$ (i.e. where left side is zero) for each coordinate axis gives you zero, one, or two real roots $t$. If a root is between $0$ and $1$, then it specifies the value of $t$ where the curve reaches an extremum (minimum or maximum) on that specific coordinate axis. Evaluating the coordinate at $t = 0$, $t = 1$, and the roots between $0$ and $1$, if any, yields the range of coordinates the curve segment spans on that specific axis.
Thus, this method yields a bounding box around the curve, with the box "axis" aligned with the line between the curve start and end points. The rotation of the box around that axis depends on how the curve segment is rotated, and is thus basically arbitrary.
If the curve is translated and rotated so that $\vec{P}_0$ is at origin, and $\vec{P}_3$ is on the positive $x$ axis, then $r = \max\left( \lvert y_\max \rvert, \lvert z_\max \rvert \right)$, i.e. the larger magnitude of the $y$ or $z$ extrema for this curve segment, is very useful: Then, $r$ is the radius of the enclosing right cylinder whose axis corresponds to the line between the curve segment start and end points. This cylinder fully contains the 3D curve.
To translate the curve so that $\vec{P}_0$ is at origin, we just translate all four $\vec{P}_k$ by $-\vec{P}_0$.
To rotate the curve so that the translated $\vec{P}_3$ is on the positive $x$ axis, we have two options:
We can pick the axis that yields the smallest rotation angle to bring translated $\vec{P}_3$ to the positive $x$ axis.
This is done (using original, untranslated coordinates) by rotating the translated points by $\theta$ around axis $\hat{n}$: $$ \begin{aligned} \hat{u} &= \frac{\vec{P}_3 - \vec{P}_0}{\left\lVert \vec{P}_3 - \vec{P}_0 \right\rVert} \\ \hat{v} &= (1, 0, 0) \\ \hat{n} &= \frac{ \hat{u} \times \hat{v}}{\left\lVert \hat{u} \times \hat{v} \right\rVert} \\ \cos \theta &= \hat{u} \cdot \hat{v} \\ \end{aligned}$$ and the rotation matrix $\mathbf{R}$ (for rotating points via $\vec{v}^\prime = \mathbf{R}\vec{v}$) is $$\mathbf{R} = \left[\begin{matrix} n_x^2 d + c & n_x n_y d - n_z c & n_x n_z d + n_y c \\ n_y n_x d + n_z c & n_y^2 d + c & n_y n_z d - n_x c \\ n_z n_x d - n_y c & n_z n_y d + n_x c & n_z^2 d + c \\ \end{matrix}\right]$$where $c = \cos\theta$ and $d = 1 - c$.
Note that you do not need any trigonometric functions (like $\cos$ or $\arccos$) here at all. Addition, subtraction, multiplication, and square root suffice.
If you are only interested in $r$, i.e. are constructing the bounding right cylinder for the curve, it is less work to just solve $\lvert y_\max \rvert$ and $\lvert z_\max \rvert$ from the roots of the Cartesian coordinate derivatives and then calculate $r = \max\left( \lvert y_\max \rvert, \lvert z_\max \rvert \right)$, than to try and rotate the curve so that the larger magnitude extremum would always correspond to say $y$ axis.
Also note that the axis of the cylinder, or the cylinder length, is not limited to the line segment between $\vec{P}_3 - \vec{P}_0$. Given $x_\min$ and $x_\max$ in the rotated coordinates, the base $\vec{b}$ and top $\vec{q}$ are actually $$\begin{aligned} w_\min &= \frac{x_\min}{\left\lVert \vec{P}_3 - \vec{P}_0 \right\rVert} \\ w_\max &= \frac{x_\max}{\left\lVert \vec{P}_3 - \vec{P}_0 \right\rVert} \\ \vec{b} &= (1 - w_\min) \vec{P}_0 + w_\min \vec{P}_3 \\ \vec{q} &= (1 - w_\max) \vec{P}_0 + w_\max \vec{P}_3 \\ \end{aligned}$$
We can pick the rotation such that the $y$ and $z$ components are in meaningful directions for us. We can use $\vec{Y}$ as the direction for the rotated $y$ axis –– the rotated $z$ axis will be perpendicular to both that and the line between $\vec{P}_0$ and $\vec{P}_3$, as long as $\vec{Y}$ is not parallel to the line between $\vec{P}_0$ and $\vec{P}_3$.
The idea is that scaling $\vec{P}_3 - \vec{P}_0$ to unit length gives us the new unit $x$ axis vector. We then use one step of Gram–Schmidt process to orthonormalize $\vec{Y}$ –– that is, we remove the component that is not perpendicular to $x$ from $Y$ to make it perpendicular to $\vec{P}_3 - \vec{P}_0$ to get the new unit $y$ axis vector. Finally, a vector cross product ($\times$) between the two will get us the third unit axis vector. These three unit vectors form an orthonormal basis, and directly form the rotation matrix we need.
In other words: $$\begin{aligned} \hat{u} &= \frac{\vec{P}_3 - \vec{P}_0}{\left\lVert \vec{P}_3 - \vec{P}_0 \right\rVert} \\ \hat{v} &= \frac{\vec{Y} - \left(\vec{Y} \cdot \hat{u}\right) \hat{u}}{\left\lVert \vec{Y} - \left(\vec{Y} \cdot \hat{u}\right) \hat{u} \right\rVert} \\ \hat{w} &= \hat{u} \times \hat{v} \\ \mathbf{R} &= \left[ \begin{matrix} \hat{u} & \hat{v} & \hat{w} \end{matrix} \right] = \left[ \begin{matrix} u_x & v_x & w_x \\ u_y & v_y & w_y \\ u_z & v_z & w_z \\ \end{matrix} \right] \\ \end{aligned}$$ but note that if $\vec{Y} \parallel \left(\vec{P}_3 - \vec{P}_0\right)$, then the denominator for $\hat{v}$ is zero, and the process fails.
One option is to have two unit direction vectors, $\hat{Y}$ and $\hat{Z}$, for rotated $y$ and $z$ respectively. Then, if $\left\lvert \hat{Y}\cdot \left(\vec{P}_3 - \vec{P}_0\right) \right\rvert \le \left\lvert \hat{Z} \cdot \left(\vec{P}_3 - \vec{P}_0\right) \right\rvert$, you do the above, otherwise $$\begin{aligned} \hat{w} &= \frac{\vec{Z} - \left(\vec{Z} \cdot \hat{u}\right) \hat{u}}{\left\lVert \vec{Z} - \left(\vec{Z} \cdot \hat{u}\right) \hat{u} \right\rVert} \\ \hat{v} &= \hat{w} \times \hat{u} \\ \mathbf{R} &= \left[ \begin{matrix} \hat{u} & \hat{v} & \hat{w} \end{matrix} \right] = \left[ \begin{matrix} u_x & v_x & w_x \\ u_y & v_y & w_y \\ u_z & v_z & w_z \\ \end{matrix} \right] \\ \end{aligned}$$ This means that you pick $\hat{Y}$ or $\hat{Z}$ based on which of them makes the larger angle with the line between $\vec{P}_0$ and $\vec{P}_3$.
With no rotation, only applying the translation, and ignoring the fact that the end point is no longer on the $x$ axis, we still obtain the axis-aligned bounding box for the curve segment. This is particularly useful when $\vec{P}_3 = \vec{P}_0$, i.e. the curve is a simple loop, as in that case neither of the two above cases apply (as $\hat{u}$ will be undefined due to division by zero).
Why is the 2D case so much simpler?
The key is that when you have an unit vector $\hat{u} = (u_x, u_y)$, $u_x^2 + u_y^2 = 1$, the other unit vector is always $\hat{v} = (-u_y, u_x)$ (assuming a right-handed coordinate system). In 3D, we don't have that, and instead need some other way to define the two other unit axis vectors.
If we are looking for axis-aligned bounding box for each curve segment, then no rotation or translation is necessary or desired. If we write $$\vec{P}_0 = \left[ \begin{matrix} X_0 \\ Y_0 \\ Z_0 \end{matrix} \right], \vec{P}_1 = \left[ \begin{matrix} X_1 \\ Y_1 \\ Z_1 \end{matrix} \right], \vec{P}_2 = \left[ \begin{matrix} X_2 \\ Y_2 \\ Z_2 \end{matrix} \right], \vec{P}_3 = \left[ \begin{matrix} X_3 \\ Y_3 \\ Z_3 \end{matrix} \right], \vec{p}(t) = \left[ \begin{matrix} X(t) \\ Y(t) \\ Z(t) \end{matrix} \right] \tag{3a}\label{3a}$$then the three equations for the curve segment, analogous to $\eqref{2b}$ except with $\vec{P}_0 \ne \vec{0}$, are $$\begin{aligned} X(t) &= (1 - t)^3 X_3 + 3 (1 - t)^2 t X_2 + 3 (1 - t) t^2 X_1 + t^3 X_0 \\ ~ &= \left( X_3 - 3 X_2 + 3 X_1 - X_0 \right) t^3 + \left( 3 X_2 - 6 X_1 + 3 X_0 \right) t^2 + \left( 3 X_1 - 3 X_0 \right) t + X_0 \\ \end{aligned}$$and similarly for $Y(t)$ and $Z(t)$; and the equations whose roots $t$ we want are $$\begin{aligned} \frac{d X(t)}{d t} = 0 &= \left( 3 X_3 - 9 X_2 + 9 X_1 - 3 X_0 \right) t^2 + \left( 6 X_2 - 12 X_1 + 6 X_0 \right) t + \left( 3 X_1 - 3 X_0 \right) \\ \end{aligned}$$ and similarly for $y$ and $z$ axes. Indeed, if this has real roots $t$, those roots are $$t = \frac{-X_2 + 2 X_1 - X_0 \pm \sqrt{ X_2 (X_2 - X_1) + X_1 (X_1 - X_3) + X_0 (X_3 - X_2) }}{X_3 - 3 X_2 + 3 X_1 - X_0}\tag{4a}\label{4a}$$ except $$t = \frac{X_0 - X_1}{2 X_2 - 4 X_1 + 2 X_0}, \quad \text{when} \quad X_3 - 3 X_2 + 3 X_1 - X_0 = 0 \tag{4b}\label{4b}$$ Only roots $0 \le t \le 1$ refer to an extremum on the $x$ axis. Similarly for $y$ and $z$, of course.
Here is an example cubic Bézier segment axis-aligned bounding box function in Python:
Each point must be a list or a tuple (which means it works in any number of dimensions).
For example, if $\vec{P}_0 = \vec{0}$, $\vec{P}_1 = (1, 2, -1)$, $\vec{P}_2 = (2, 2, 3)$, and $\vec{P}_3 = (3, 0, 2)$, the axis-aligned bounding box is $(0, 0, -0.162) - (3, 1.500, 2.162)$. (Run
print(Bezier_AABB([0,0,0], [1,2,-1], [2,2,3], [3,0,2]))to verify.)Another example, a particularly intriguing one, a 2D droplet: $\vec{P}_0 = \vec{P}_3 = \vec{0}$, $\vec{P}_1 = (-1, 1)$, $\vec{P}_2 = (1,1)$ is contained in axis-aligned bounding box $(-0.289, 0.000) - (0.289, 0.750)$. (Run
print(Bezier_AABB([0,0], [-1,1], [1,1], [0,0]))to verify.)In my personal opinion as someone who extensively uses cubic Bézier curves, it is this axis-aligned bounding box approach (that involves no rotation or translation, just examination of the ranges of coordinates along each Cartesian coordinate axis that the curve spans) that is most useful.
In 3D, the cylinder approach (solving for $r$) is also useful, with the limitation that the cylinder axis is always the line between the start and end points of the curve segment. When the start and end points coincide, you can calculate e.g. the bounding sphere instead.