Approximating Catmull-Rom Curves with Quadratic Beziers

See QuadraticApproximation.pdf in this directory for the math underlying this approach.
In this notebook, I’m figuring out a simple algebraic representation of the equations in the paper that specialize the code to specific parameters that specify a centripetal Catmull-Rom spline (α = 0.5) and minimize distortion (γ = 0.5). With these parameter values, the resulting equations are simpler and we can avoid doing extra math at runtime.

Quadratic Approximation of Cubic Curves
NGHIA TRUONG, University of Utah
CEM YUKSEL, University of Utah
LARRY SEILER, Facebook Reality Labs
https://ttnghia.github.io/pdf/QuadraticApproximation.pdf

α=(1/2);(* set the alpha parameter to 0.5 for a centripetal catmull-rom spline *)
γ=(1/2); (* split point parameter quadratic approx of cubic bezier. 1/2 simplifies the math, and is the parameter value that minimizes distortion; the two middle endpoints of the quadratic curves will lie on the approximated cubic. *)

(* distances *)
dist[a_,b_]:=Norm[b-a]
(*dist[a_,b_]:=Sqrt[Dot[a-b,a-b]]*)
d1 = dist[p0,p1];
d2 = dist[p1,p2];
d3 = dist[p2,p3];

(*Section 6: Given four consecutive data points p0,p1,p2,and p3,the Bézier form of a cubic Catmull-Rom curve segment that interpolates p1 and p2 can be written using the Bézier control points: *)
b0 = p1;
b1 =p1 +(((d1^(2α))(p2-p1) + (d2^(2α))(p1-p0))/(3(d1^α)((d1^α)+(d2^α))));
b2=p2+(((d3^(2α))(p1-p2)+(d2^(2α))(p2-p3))/(3(d3^α)((d3^α)+(d2^α))));
b3=p2;

(* Using our degree reduction in Equations 13-15,we
can write the control points of the corresponding pair of quadratic segments: *)
q00=p1;
q01=p1+γ((((d1^(2α))(p2-p1)+(d2^(2α))(p1-p0))/(2(d1^α)((d1^α)+(d2^α)))))
q02=q10=(1-γ)q01-γ q11;
q11=p2+(1-γ)(((d3^(2α))(p1-p2)+(d2^(2α))(p2-p3))/(2(d3^α)((d3^α)+(d2^α))))
q12=p2;
(* The six control points of the two resulting quadratic beziers are q01..q02 and q11..q12. *)

(* To degree-lower a cubic curve to two quadratic curves, see eqns. 13, 14, and 15, and the diagram above that provides explanations for the labels. *)

Module[{b0=cp0, b1=cp1, b2=cp2, b3=cp3, q1, q2, q3},
q1=b0 + (3/2)γ(b1-b0);
q3=b3+(3/2)(1-γ)(b2-b3);
q2=(1-γ)q1+γ(q3);
{
{b0, q1, q2},
{q2,q3, b3}
}
]//Simplify
%//MatrixForm