# Normal of Catmull-Rom Spline

I've implemented Catmull-Rom in python and now I want to be able to get the normal of points so I'm first finding the tangent, to look something like this. For Catmull Rom I have 4 points P0, P1, P2 and P3. Using cfh's answer I've found the tangents of the two end points, P1 and P2

  #m1 represents tangent at starting point P1
m1 = (t2-t1)*(((P1-P0)/(t1-t0))-((P2-P0)/(t2-t0))+((P2-P1)/(t2-t1)))
#m2 epresents tangent at ending point P2
m1 = (t2-t1)*(((P2-P1)/(t2-t1))-((P3-P1)/(t3-t1))+((P3-P2)/(t3-t2)))


I've then put these into the standard formula for a cubic spline following this post

cubic_spline = (2*t**3 - 3*t**2 + 1) * P1 + (t**3 - 2*t**2 + t) * m1 + (-2*t**3 + 3*t**2) * P2 + (t**3 - t**2) * m2
cubic_spline_deriv = (6*t**2 - 6*t)*P1 + (3*t**2 - 4*t + 1)*m1 + (-6*t**2 + 6*t)*P2 + (3*t**2 - 2*t)*m2


I think I'm massively misunderstanding the maths though.

I think that the output of cubic_spline_deriv is a tangent vector. However, when I try to plot it with plt.quiver I'm way off finding a tangent. I'd also expect cubic_spline to give the same output as C, but it doesn't.

Would it be at all possible for anyone to please take a look and see where I'm massively misunderstanding please?

My whole code is below:

    #https://stackoverflow.com/questions/34894837/how-to-get-connect-two-part-of-curve-and-get-the-points-position-of-connecting-c
##https://stackoverflow.com/questions/37214786/emulating-excels-scatter-with-smooth-curve-spline-function-in-matplotlib-for
def CatmullRomSpline(P0, P1, P2, P3, nPoints=100):
"""
P0, P1, P2, and P3 should be (x,y) point pairs that define the Catmull-Rom spline.
nPoints is the number of points to include in this curve segment.
"""
# Convert the points to numpy so that we can do array multiplication
P0, P1, P2, P3 = map(np.array, [P0, P1, P2, P3])

# Calculate t0 to t4
alpha = 0.5
def tj(ti, Pi, Pj):
xi, yi = Pi
xj, yj = Pj
return ( ( (xj-xi)**2 + (yj-yi)**2 )**0.5 )**alpha + ti

t0 = 0
t1 = tj(t0, P0, P1)
t2 = tj(t1, P1, P2)
t3 = tj(t2, P2, P3)

# Only calculate points between P1 and P2
t = np.linspace(t1,t2,nPoints)

# Reshape so that we can multiply by the points P0 to P3
# and get a point for each value of t.
t = t.reshape(len(t),1)
#print(t)
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
#print(A1)
#print(A2)
#print(A3)
B1 = (t2-t)/(t2-t0)*A1 + (t-t0)/(t2-t0)*A2
B2 = (t3-t)/(t3-t1)*A2 + (t-t1)/(t3-t1)*A3

C  = (t2-t)/(t2-t1)*B1 + (t-t1)/(t2-t1)*B2

#m1 represents tangent at starting point P1
m1 = (t2-t1)*(((P1-P0)/(t1-t0))-((P2-P0)/(t2-t0))+((P2-P1)/(t2-t1)))
#m2 epresents tangent at ending point P2
m2 = (t2-t1)*(((P2-P1)/(t2-t1))-((P3-P1)/(t3-t1))+((P3-P2)/(t3-t2)))

cubic_spline = (2*t**3 - 3*t**2 + 1) * P1 + (t**3 - 2*t**2 + t) * m1 + (-2*t**3 + 3*t**2) * P2 + (t**3 - t**2) * m2
cubic_spline_deriv = (6*t**2 - 6*t)*P1 + (3*t**2 - 4*t + 1)*m1 + (-6*t**2 + 6*t)*P2 + (3*t**2 - 2*t)*m2

return C, cubic_spline, cubic_spline_deriv

def CatmullRomChain(P):
"""
Calculate Catmull Rom for a chain of points and return the combined curve.
"""
sz = len(P)

# The curve C will contain an array of (x,y) points.
#C is a list of x and y coordinates
C = []
Q = []
Qd = []
for i in range(sz-3):
c, cubic_spline, cubic_spline_deriv = CatmullRomSpline(P[i], P[i+1], P[i+2], P[i+3])
C.extend(c)
Q.extend(cubic_spline)
Qd.extend(cubic_spline_deriv)

return C, Q, Qd

# Define a set of points for curve to go through
Points = [[0,1.5],[2,2],[3,1],[4,0.5],[5,1],[6,2],[7,3]]

# Calculate the Catmull-Rom splines through the points
c, q, qd = CatmullRomChain(Points)

u = qd[0][0]
v = qd[0][1]
print(c)
print(q)
print(len(qd))

# Convert the Catmull-Rom curve points into x and y arrays and plot
x,y = zip(*c)
plt.plot(x,y)

# Plot the control points
px, py = zip(*Points)
plt.plot(px,py,'or')
plt.quiver(2, 2, u, v)

plt.show()


Catmull Rom curve

• Welcome to Math.SE, @Emily. I'm not sure that here is the best place to post a request to code review. Perhaps another site around here, like Code Review, or Computer Science or one of the ones you cited posts from like Stack Overflow or Game Development. Commented May 5, 2019 at 18:53
• Thank you. I think it's more the maths I'm struggling with though rather than the code. Would it help if I posted in plain text rather than code format? Commented May 5, 2019 at 18:55
• I think that using the code format it's the best way to check code. Commented May 5, 2019 at 18:57
• I did not check your code, but I took a look at the comments on chf's post and one said "plasmacel: The approach from your link is simply wrong for the non-uniform case, it just happens to reduce to the correct solution if all time intervals are 1. And that's why you cannot match it with Barry-Goldman. The left and right time interval must appear separately in the equation, otherwise it can never work correctly! – Matthias Jul 5 '18 at 9:32" Is it relevant for your case? Commented May 5, 2019 at 22:33
• Thank you for your help. I think plasmacel is not the original commentator and was providing another parameterisation in the comments so cfh, the OP's, answer still stands. I rescaled the tangent's by multiplying by (t2-t1). I'm just confused because it's managed to be done here but with the tangent being calculated as M[k] = (P[k+1] - P[k-1])/2 . I think I don't understand how a tangent can have two outputs, an x and y? Commented May 6, 2019 at 7:27

if catmull-rom have v0 v1 v2 v3..

float t2 = t * t;

float m0 = (v2 - v0) * .5f;

float m1 = (v3 - v1) * .5f;

tangent = (t2 - t) * 6f * v1 + (3f * t2 - 4f * t + 1f) * m0 + (-6f * t2 + 6f * t) * v2 + (3f * t2 - 2f * t) * m1;

if you want normalized >> tangent.normalized

if you want normal vector >> rotate right 90 degree

if tangent is Vector2, Very Simple >> normal = (tangent.y, -tangent.x)

Have a nice day!

• This worked, thank you! Commented Dec 19, 2023 at 22:31