View source for Centripetal Catmull–Rom spline

{{Short description|Variant form of the Catmull-Rom spine}}
In [[computer graphics]], the '''centripetal Catmull–Rom spline''' is a variant form of the [[Catmull–Rom spline]], originally formulated by [[Edwin Catmull]] and [[Raphael Rom]],<ref>{{Cite book |title=Computer Aided Geometric Design |last1=Catmull |first1=Edwin |last2=Rom |first2=Raphael |date=1974 |isbn=978-0-12-079050-0 |editor-last=Barnhill |editor-first=Robert E. |pages=317–326 |chapter=A class of local interpolating splines |doi=10.1016/B978-0-12-079050-0.50020-5 |author-link=Edwin Catmull |author-link2=Raphael Rom |editor-last2=Riesenfeld |editor-first2=Richard F.}}</ref> which can be evaluated using a recursive algorithm proposed by Barry and Goldman.<ref>{{Cite conference|last1=Barry |first1=Phillip J. |last2=Goldman |first2=Ronald N. |date=August 1988|title=A recursive evaluation algorithm for a class of Catmull–Rom splines |conference=Proceedings of the 15st Annual Conference on Computer Graphics and Interactive Techniques, [[SIGGRAPH]] 1988|doi=10.1145/378456.378511|publisher=[[Association for Computing Machinery]]|url=https://www.researchgate.net/publication/220720141 |volume=22 |issue=4 |pages=199–204}}</ref> It is a type of interpolating spline (a curve that goes through its control points) defined by four control points <math>\mathbf{P}_0, \mathbf{P}_1, \mathbf{P}_2, \mathbf{P}_3</math>, with the curve drawn only from <math>\mathbf{P}_1</math> to <math>\mathbf{P}_2</math>.

[[File:Catmull-Rom Spline.png|thumb|Catmull–Rom spline interpolation with four points]]

==Definition==
[[File:Cubic Catmull-Rom formulation.png|thumb|Barry and Goldman's pyramidal formulation]]
[[File:Catmull-Rom Parameterized Time.png|thumb|Knot parameterization for the Catmull–Rom algorithm]]

Let <math>\mathbf{P}_i = [x_i \quad y_i]^T</math> denote a point. For a curve segment <math>\mathbf{C}</math> defined by points <math>\mathbf{P}_0, \mathbf{P}_1, \mathbf{P}_2, \mathbf{P}_3</math> and knot sequence <math>t_0, t_1, t_2, t_3</math>, the centripetal Catmull–Rom spline can be produced by:

: <math>\mathbf{C} = \frac{t_2-t}{t_2-t_1}\mathbf{B}_1+\frac{t-t_1}{t_2-t_1}\mathbf{B}_2</math>

where

: <math>\mathbf{B}_1 = \frac{t_{2}-t}{t_{2}-t_0}\mathbf{A}_1+\frac{t-t_0}{t_{2}-t_0}\mathbf{A}_2</math>
: <math>\mathbf{B}_2 = \frac{t_{3}-t}{t_{3}-t_1}\mathbf{A}_2+\frac{t-t_1}{t_{3}-t_1}\mathbf{A}_3</math>
: <math>\mathbf{A}_1 = \frac{t_{1}-t}{t_{1}-t_0}\mathbf{P}_0+\frac{t-t_0}{t_{1}-t_0}\mathbf{P}_1</math>
: <math>\mathbf{A}_2 = \frac{t_{2}-t}{t_{2}-t_1}\mathbf{P}_1+\frac{t-t_1}{t_{2}-t_1}\mathbf{P}_2</math>
: <math>\mathbf{A}_3 = \frac{t_{3}-t}{t_{3}-t_2}\mathbf{P}_2+\frac{t-t_2}{t_{3}-t_2}\mathbf{P}_3</math>

and

:<math>t_{i+1} = \left[\sqrt{(x_{i+1}-x_i)^2+(y_{i+1}-y_i)^2}\right]^{\alpha} + t_i</math>

in which <math>\alpha</math> ranges from 0 to 1 for knot parameterization, and <math>i = 0,1,2,3</math> with <math>t_0 = 0 </math>. For centripetal Catmull–Rom spline, the value of <math>\alpha</math> is <math>0.5</math>. When <math>\alpha = 0</math>, the resulting curve is the standard [[uniform Catmull–Rom spline]]; when <math>\alpha = 1</math>, the result is a [[chordal Catmull–Rom spline]].
[[File:Uniform, Centripetal and Chordal Parameterization for Catmull-Rom Spline.gif|thumb|Gif animation for ''uniform'', ''centripetal'' and ''chordal'' parameterization of Catmull–Rom spline depending on the ''α'' value]]
Plugging <math>t = t_1</math> into the spline equations <math> \mathbf{A}_1, \mathbf{A}_2, \mathbf{A}_3, \mathbf{B}_1, \mathbf{B}_2,</math> and <math> \mathbf{C}</math> shows that the value of the spline curve at <math>t_1</math> is <math>\mathbf{C} = \mathbf{P}_1</math>. Similarly, substituting <math>t = t_2</math> into the spline equations shows that <math> \mathbf{C} = \mathbf{P}_2 </math> at <math>t_2</math>. This is true independent of the value of <math>\alpha</math> since the equation for <math>t_{i+1}</math> is not needed to calculate the value of <math>\mathbf{C}</math> at points <math>t_1</math> and <math>t_2</math>.
[[File:3D Centripetal Catmull-Rom Spline segment.png|alt=|thumb|3D ''centripetal'' Catmull-Rom spline segment.]]
The extension to 3D points is simply achieved by considering <math>\mathbf{P}_i = [x_i \quad y_i \quad z_i]^T</math>a generic 3D point <math>\mathbf{P}_i</math> and

:<math>t_{i+1} = \left[\sqrt{(x_{i+1}-x_i)^2+(y_{i+1}-y_i)^2+(z_{i+1}-z_i)^2}\right]^{\alpha} + t_i</math>

==Advantages==
Centripetal Catmull–Rom splines have several desirable mathematical properties compared to the original and the other types of Catmull-Rom formulation.<ref name="Yuksel">{{Cite journal |last1=Yuksel |first1=Cem |last2=Schaefer |first2=Scott |last3=Keyser |first3=John |date=July 2011 |title=Parameterization and applications of Catmull-Rom curves |url=http://www.cemyuksel.com/research/catmullrom_param/ |journal=Computer-Aided Design |volume=43 |issue=7 |pages=747–755 |doi=10.1016/j.cad.2010.08.008|citeseerx=10.1.1.359.9148 }}</ref> First, loops and self-intersections do not occur within a curve segment. Second, [[cusp_(singularity)|cusps]] will never occur within a curve segment. Third, it follows the control points more tightly.<ref>{{Cite web|last=Yuksel; Schaefer; Keyser|first=Cem; Scott; John|title=On the Parameterization of Catmull-Rom Curves|url=http://www.cemyuksel.com/research/catmullrom_param/catmullrom.pdf}}</ref>{{vague|date=September 2016}}

[[File:Catmull-Rom examples with parameters..png|thumb|In this figure, there is a self-intersection/loop on the uniform Catmull-Rom spline (green), whereas for chordal Catmull-Rom spline (red), the curve does not follow tightly through the control points.]]

==Other uses==
In [[computer vision]], centripetal Catmull-Rom spline has been used to formulate an active model for segmentation. The method is termed '''active spline model'''.<ref>{{Cite journal |last1=Jen Hong |first1=Tan |last2=Acharya |first2=U. Rajendra |date=2014 |title=Active spline model: A shape based model-interactive segmentation |url=https://doppiomovimento.files.wordpress.com/2015/01/active-spline-model-s.pdf |journal=[[Digital Signal Processing (journal)|Digital Signal Processing]] |volume=35 |pages=64–74 |arxiv=1402.6387 |doi=10.1016/j.dsp.2014.09.002|s2cid=6953844 }}</ref> The model is devised on the basis of [[active shape model]], but uses centripetal Catmull-Rom spline to join two successive points (active shape model uses simple straight line), so that the total number of points necessary to depict a shape is less. The use of centripetal Catmull-Rom spline makes the training of a shape model much simpler, and it enables a better way to edit a contour after segmentation.

==Code example in Python==
The following is an implementation of the Catmull–Rom spline in [[Python (programming language)|Python]] that produces the plot shown beneath.

<syntaxhighlight lang="numpy">
import numpy
import matplotlib.pyplot as plt


QUADRUPLE_SIZE: int = 4


def num_segments(point_chain: tuple) -> int:
    # There is 1 segment per 4 points, so we must subtract 3 from the number of points  
    return len(point_chain) - (QUADRUPLE_SIZE - 1)


def flatten(list_of_lists) -> list:
    # E.g. mapping [[1, 2], [3], [4, 5]] to  [1, 2, 3, 4, 5] 
    return [elem for lst in list_of_lists for elem in lst]


def catmull_rom_spline(
    P0: tuple,
    P1: tuple,
    P2: tuple,
    P3: tuple,
    num_points: int,
    alpha: float = 0.5,
):
    """
    Compute the points in the spline segment
    :param P0, P1, P2, and P3: The (x,y) point pairs that define the Catmull-Rom spline
    :param num_points: The number of points to include in the resulting curve segment
    :param alpha: 0.5 for the centripetal spline, 0.0 for the uniform spline, 1.0 for the chordal spline.
    :return: The points
    """

    # Calculate t0 to t4. Then only calculate points between P1 and P2.
    # Reshape linspace so that we can multiply by the points P0 to P3
    # and get a point for each value of t.
    def tj(ti: float, pi: tuple, pj: tuple) -> float:
        xi, yi = pi
        xj, yj = pj
        dx, dy = xj - xi, yj - yi
        l = (dx ** 2 + dy ** 2) ** 0.5
        return ti + l ** alpha

    t0: float = 0.0
    t1: float = tj(t0, P0, P1)
    t2: float = tj(t1, P1, P2)
    t3: float = tj(t2, P2, P3)
    t = numpy.linspace(t1, t2, num_points).reshape(num_points, 1)

    A1 = (t1 - t) / (t1 - t0) * P0 + (t - t0) / (t1 - t0) * P1
    A2 = (t2 - t) / (t2 - t1) * P1 + (t - t1) / (t2 - t1) * P2
    A3 = (t3 - t) / (t3 - t2) * P2 + (t - t2) / (t3 - t2) * P3
    B1 = (t2 - t) / (t2 - t0) * A1 + (t - t0) / (t2 - t0) * A2
    B2 = (t3 - t) / (t3 - t1) * A2 + (t - t1) / (t3 - t1) * A3
    points = (t2 - t) / (t2 - t1) * B1 + (t - t1) / (t2 - t1) * B2
    return points


def catmull_rom_chain(points: tuple, num_points: int) -> list:
    """
    Calculate Catmull-Rom for a sequence of initial points and return the combined curve.
    :param points: Base points from which the quadruples for the algorithm are taken
    :param num_points: The number of points to include in each curve segment
    :return: The chain of all points (points of all segments)
    """
    point_quadruples = (  # Prepare function inputs
        (points[idx_segment_start + d] for d in range(QUADRUPLE_SIZE))
        for idx_segment_start in range(num_segments(points))
    )
    all_splines = (catmull_rom_spline(*pq, num_points) for pq in point_quadruples)
    return flatten(all_splines)


if __name__ == "__main__":
    POINTS: tuple = ((0, 1.5), (2, 2), (3, 1), (4, 0.5), (5, 1), (6, 2), (7, 3))  # Red points
    NUM_POINTS: int = 100  # Density of blue chain points between two red points

    chain_points: list = catmull_rom_chain(POINTS, NUM_POINTS)
    assert len(chain_points) == num_segments(POINTS) * NUM_POINTS  # 400 blue points for this example

    plt.plot(*zip(*chain_points), c="blue")
    plt.plot(*zip(*POINTS), c="red", linestyle="none", marker="o")
    plt.show()

</syntaxhighlight>
[[File:Catmull–Rom-spline as per the Python implementation in English Wikipedia Catmull–Rom-spline article.svg|thumb|left|Plot obtained by the Python example code given above]]
{{Clear}}

==Code example in Unity C#==
<syntaxhighlight lang="csharp">
using UnityEngine;

// a single catmull-rom curve
public struct CatmullRomCurve
{
	public Vector2 p0, p1, p2, p3;
	public float alpha;

	public CatmullRomCurve(Vector2 p0, Vector2 p1, Vector2 p2, Vector2 p3, float alpha)
   {
		(this.p0, this.p1, this.p2, this.p3) = (p0, p1, p2, p3);
		this.alpha = alpha;
	}

	// Evaluates a point at the given t-value from 0 to 1
	public Vector2 GetPoint(float t)
    {
		// calculate knots
		const float k0 = 0;
		float k1 = GetKnotInterval(p0, p1);
		float k2 = GetKnotInterval(p1, p2) + k1;
		float k3 = GetKnotInterval(p2, p3) + k2;

		// evaluate the point
		float u = Mathf.LerpUnclamped(k1, k2, t);
		Vector2 A1 = Remap(k0, k1, p0, p1, u);
		Vector2 A2 = Remap(k1, k2, p1, p2, u);
		Vector2 A3 = Remap(k2, k3, p2, p3, u);
		Vector2 B1 = Remap(k0, k2, A1, A2, u);
		Vector2 B2 = Remap(k1, k3, A2, A3, u);
		return Remap(k1, k2, B1, B2, u);
	}

	static Vector2 Remap(float a, float b, Vector2 c, Vector2 d, float u)
    {
		return Vector2.LerpUnclamped(c, d, (u - a) / (b - a));
	}

	float GetKnotInterval(Vector2 a, Vector2 b)
    {
		return Mathf.Pow(Vector2.SqrMagnitude(a - b), 0.5f * alpha);
	}
}
</syntaxhighlight>


<syntaxhighlight lang="csharp">
using UnityEngine;

// Draws a catmull-rom spline in the scene view,
// along the child objects of the transform of this component
public class CatmullRomSpline : MonoBehaviour
{
	[Range(0, 1)]
    public float alpha = 0.5f;
	int PointCount => transform.childCount;
	int SegmentCount => PointCount - 3;
	Vector2 GetPoint(int i) => transform.GetChild(i).position;

	CatmullRomCurve GetCurve(int i)
    {
		return new CatmullRomCurve(GetPoint(i), GetPoint(i+1), GetPoint(i+2), GetPoint(i+3), alpha);
	}

	void OnDrawGizmos()
    {
		for (int i = 0; i < SegmentCount; i++)
			DrawCurveSegment(GetCurve(i));
	}

	void DrawCurveSegment(CatmullRomCurve curve)
    {
		const int detail = 32;
		Vector2 prev = curve.p1;

		for (int i = 1; i < detail; i++)
        {
			float t = i / (detail - 1f);
			Vector2 pt = curve.GetPoint(t);
			Gizmos.DrawLine(prev, pt);
			prev = pt;
		}
	}
}
</syntaxhighlight>

==Code example in Unreal C++==
<syntaxhighlight lang="c++">

float GetT( float t, float alpha, const FVector& p0, const FVector& p1 )
{
    auto d  = p1 - p0;
    float a = d | d; // Dot product
    float b = FMath::Pow( a, alpha*.5f );
    return (b + t);
}

FVector CatmullRom( const FVector& p0, const FVector& p1, const FVector& p2, const FVector& p3, float t /* between 0 and 1 */, float alpha=.5f /* between 0 and 1 */ )
{
    float t0 = 0.0f;
    float t1 = GetT( t0, alpha, p0, p1 );
    float t2 = GetT( t1, alpha, p1, p2 );
    float t3 = GetT( t2, alpha, p2, p3 );
    t = FMath::Lerp( t1, t2, t );
    FVector A1 = ( t1-t )/( t1-t0 )*p0 + ( t-t0 )/( t1-t0 )*p1;
    FVector A2 = ( t2-t )/( t2-t1 )*p1 + ( t-t1 )/( t2-t1 )*p2;
    FVector A3 = ( t3-t )/( t3-t2 )*p2 + ( t-t2 )/( t3-t2 )*p3;
    FVector B1 = ( t2-t )/( t2-t0 )*A1 + ( t-t0 )/( t2-t0 )*A2;
    FVector B2 = ( t3-t )/( t3-t1 )*A2 + ( t-t1 )/( t3-t1 )*A3;
    FVector C  = ( t2-t )/( t2-t1 )*B1 + ( t-t1 )/( t2-t1 )*B2;
    return C;
}
</syntaxhighlight>

==See also==
* [[Catmull–Rom spline#Catmull–Rom spline|Cubic Hermite splines]]

==References==
{{Reflist}}

==External links==
* [https://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections/19283471#19283471  Catmull-Rom curve with no cusps and no self-intersections]{{snd}} implementation in Java
* [https://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections/23980479#23980479 Catmull-Rom curve with no cusps and no self-intersections]{{snd}} simplified implementation in C++
* [http://nbviewer.jupyter.org/github/empet/geom_modeling/blob/master/Catmull-Rom-splines.ipynb Catmull-Rom splines]{{snd}} interactive generation via Python, in a Jupyter notebook
* [https://qroph.github.io/2018/07/30/smooth-paths-using-catmull-rom-splines.html Smooth Paths Using Catmull-Rom Splines]{{snd}} another versatile implementation in C++ including centripetal CR splines


{{DEFAULTSORT:Centripetal Catmull-Rom spline}}
[[Category:Splines (mathematics)]]
[[Category:Articles with example C Sharp code]]
[[Category:Articles with example Python (programming language) code]]