find maxima on a cardinal (catmull rom) spline I am looking to find the maximum deviation from the straight line joining control points on a cardinal spline. This link gives a polynomial equation for the spline:
$$q(t) = 0.5 \cdot(   (2 \cdot P_1) + (-P_0 + P_2) \cdot t + (2\cdot P_0 - 5\cdot P_1 + 4\cdot P_2 - P_3)  \cdot t^2 + (-P_0 + 3\cdot P_1- 3 \cdot P_2 + P_3) \cdot t^3)$$
The only input I have is the 4 points in 3 dimensions.
My idea was to compute the tangent by taking the first derivative of the above equation so that I get the slope at any point on the spline. Then I could compare it with the slope of the line joining the two points. But I am kind of confused how to do it in 3D. How do I compare the slopes in 3D? I can get a unit vector in the direction of the straight line joining the points, but how can I get a unit vector at the tangent. Is there a simpler way to do this?
 A: I don't think your idea will work in 3D. At the point that's furthest away from the straight line, the curve tangent might not be parallel to the line. Counterexample: take the curve with Bezier control points $\mathbf{P}_0 = (-1,0,0)$, $\mathbf{P}_1 = (0,1,1)$, $\mathbf{P}_2 = (0,1,-1)$, $\mathbf{P}_3 = (1,0,0)$. The chord $\mathbf{P}_0\mathbf{P}_3$ is the $x$-axis. By symmetry, the maximum deviation of the curve from the $x$-axis must occur at $t = \tfrac12$, but the tangent vector at this point is not parallel to the $x$-axis.
You can write down an expression that gives the distance $d(t)$ between the line and a point $\mathbf{Q}(t)$ on your curve. See here if you need help with that. Then you need to find where $d(t)$ is a maximum, which you can do using derivatives, as you described. 
A: Let's look at an approach that works for all cubic curves:
$$\vec{p}(t) = \vec{p}_0 + t \vec{p}_1 + t^2 \vec{p}_2 + t^3 \vec{p}_3 \tag{1}\label{1}$$
with $0 \le t \le 1$, $t \in \mathbb{R}$.
The straight line is from $ \vec{p}(0) = \vec{p}_0 $ to $ \vec{p}(1) = \vec{p}_0 + \vec{p}_1 + \vec{p}_2 + \vec{p}_3 $.
Let the direction vector for the line be $(x_E , y_E , z_E ) = \vec{p}(1) - \vec{p}(0) = \vec{p}_1 + \vec{p}_2 + \vec{p}_3$.
First, we translate and rotate the coordinate system so that this line will be on the $x$ axis. The translation is by $-\vec{p}_0$, and the orthonormal rotation matrix is
$$\mathbf{R} = \left [ \begin{array}{ccc}
r_{11} & r_{12} & r_{13} \\
r_{21} & r_{22} & r_{23} \\
r_{31} & r_{32} & r_{33} \end{array} \right ]$$
where$$r_{11} = \frac{ x_E }{\sqrt{ x_E^2 + y_E^2 + z_E^2 }}$$
$$r_{12} = \frac{ y_E }{\sqrt{ x_E^2 + y_E^2 + z_E^2 }}$$
$$r_{13} = \frac{ z_E }{\sqrt{ x_E^2 + y_E^2 + z_E^2 }}$$
We have an underdetermined situation: we do not care how the other two axes end up, as long as the matrix itself is valid (orthonormal). So, use the Gram-Schmidt process to find the second row vector ($r_{21}$, $r_{22}$, and $r_{23}$) starting with any initial vector not parallel to $(x_E , y_E , z_E )$; remember to normalize the result to unit length. The third row vector is the cross product of the two first row vectors.
Translate and rotate the control points, to get new ones:
$$\begin{cases}
\vec{c}_0 = (0, 0, 0) \\
\vec{c}_1 = \mathbf{R} \left ( \vec{p}_1 - \vec{p}_0 \right ) \\
\vec{c}_2 = \mathbf{R} \left ( \vec{p}_2 - \vec{p}_0 \right ) \\
\vec{c}_3 = \mathbf{R} \left ( \vec{p}_3 - \vec{p}_0 \right ) \end{cases}$$
In the transformed coordinates, the curve is now
$$\vec{c} = t \vec{c}_1 + t^2 \vec{c}_2 + t^3 \vec{c}_3 \tag{2}\label{2}$$
Using Cartesian coordinates, we can write $\eqref{2}$ as
$$\begin{cases}
x(t) = t x_1 + t^2 x_2 + t^3 x_3 \\
y(t) = t y_1 + t^2 y_2 + t^3 y_3 \\
z(t) = t z_1 + t^2 z_2 + t^3 z_3 \end{cases}$$
where $\vec{c}_1 = ( x_1 , y_1 , z_1 )$, $\vec{c}_2 = ( x_2 , y_2 , z_2 )$, and $\vec{c}_3 = ( x_3 , y_3 , z_3 )$.
Remember, $0 \le t \le 1$. We also know that $\vec{c}(0) = (0, 0, 0)$, and that $\vec{c}(1) = (L, 0, 0)$, where $L = \sqrt{x_E^2 + y_E^2 + z_E^2}$. This means that $x(t)$ is the distance from origin along the line.
The distance to from point $\vec{c}(t)$ to the line ($x$-axis) is
$$r(t) = \sqrt{ y(t)^2 + z(t)^2 } \tag{3}\label{3}$$
and we are looking for the extrema of $\eqref{3}$ within $0 \le t \le 1$. Because $r(t)$ is a nonnegative function, it reaches its extrema where its square does. Therefore, we can actually look for the extremas of
$$R(t) = r(t)^2 = y(t)^2 + z(t)^2$$
which occur when $t = 0$, $t = 1$, or the derivative is zero;
$$\begin{array}{rl}
\frac{d \, R(t)}{d \, t} = 0 = & t^5 ( 6 y_3^2 + 6 z_3^2 ) + \\
\; & t^4 ( 10 y_2 y_3 + 10 z_2 z_3 ) + \\
\; & t^3 ( 4 y_2^2 + 4 z_2^2 + 8 y_1 y_3 + 8 z_1 z_3 ) + \\
\; & t^2 ( 6 y_1 y_2 + 6 z_1 z_2 ) + \\
\; & t ( 2 y_1^2 + z_1^2 ) \end{array}\tag{4}\label{4}$$
Note that $t = 0$ is always a root; $t = 1$ is mathematically also a root (because of how the coefficients have been defined; but not in general for the above equation!), but rounding errors may perturb the root at $1$ slightly.
$\eqref{4}$ does have analytical solutions, but you'll probably want to use Sage, Maple, Mathematica, or something similar to massage the solutions into more useful form; the solutions I found using Sage are just too long to paste here.
Let's say you find roots $t_i$, $0 \le t_i \le 1$. Pick $t_{max} = t_k$ where $R(t_k) \ge R(t_i)$ for all $i$; i.e. pick the $t$ that yields the largest $R(t)$. Then:


*

*The point where the curve is furthest from the straight line occurs at $t_{max}$.

*In original coordinates, the curve is farthest from the straight line at $\vec{p}(t_{max})$.

*The distance is then $r(t_{max}) = \sqrt{R(t_{max})}$.

*Distance from $\vec{p}(0)$ to the point on the line that is perpendicular to, and closest to the furthest point, is $x(t_{max})$.
