Sub-curve of a Hermite Curve does not fit original curve
I’m trying to use a Hermite curve in a project, with an admittedly limited understanding of how the math actually works, and I’ve run into some behavior I don’t understand. I've demonstrated my confusion with a minimal code sample below, but basically I would expect points along a subcurve of a hermite curve (i.e. a subcurve defined using points and tangents on the original curve) to fit the original curve, but this seems to be false.
The following c# code defines a Hermite curve class that provides functions for computing the position and the tangent of a point at some ratio along the curve. I copy/pasted the math for both functions from other places on the internet.
A small test harness then performs the test that I would expect to succeed, but doesn’t. It is unclear to me if there is a bug in my code, a mistake in my math, or if I misunderstand something about how Hermite curves work and this test actually should not pass.
Any insight is appreciated.
using System;
using System.Numerics;
class Program
{
class HermiteCurve
{
Vector2 start;
Vector2 startTangent;
Vector2 end;
Vector2 endTangent;
public HermiteCurve(Vector2 start, Vector2 startTangent, Vector2 end, Vector2 endTangent)
{
this.start = start;
this.startTangent = startTangent;
this.end = end;
this.endTangent = endTangent;
}
public Vector2 GetPoint(float t)
{
var t2 = t * t;
var t3 = t2 * t;
return
( 2f*t3 - 3f*t2 + 1f) * start +
(-2f*t3 + 3f*t2) * end +
(t3 - 2f*t2 + t) * startTangent +
(t3 - t2) * endTangent;
}
public Vector2 GetTangent(float t)
{
var t2 = t * t;
return
(6f*t2 - 6*t) * start +
(-6f*t2 + 6*t) * end +
(3f*t2 - 4f*t + 1) * startTangent +
(3f*t2 - 2f*t) * endTangent;
}
}
static void Main(string[] args)
{
Vector2 p0 = new Vector2(0, 0);
Vector2 m0 = new Vector2(1, 0);
Vector2 p1 = new Vector2(1, 1);
Vector2 m1 = new Vector2(0, 1);
HermiteCurve curve = new HermiteCurve(p0, m0, p1, m1);
Vector2 p0prime = curve.GetPoint(0.5f);
Vector2 m0prime = curve.GetTangent(0.5f);
HermiteCurve curvePrime = new HermiteCurve(p0prime, m0prime, p1, m1);
Vector2 curvePoint = curve.GetPoint(0.75f);
Vector2 curveTangent = curve.GetTangent(0.75f);
Vector2 curvePrimePoint = curvePrime.GetPoint(0.5f);
Vector2 curvePrimeTangent = curvePrime.GetTangent(0.5f);
// Why does this check fail?
if (curvePoint != curvePrimePoint || curveTangent != curvePrimeTangent)
{
Console.WriteLine("fail");
Console.WriteLine("curvePosition - x: " + curvePoint.X + " y: " + curvePoint.Y);
Console.WriteLine("curvePrimePosition - x: " + curvePrimePoint.X + " y: " + curvePrimePoint.Y);
Console.WriteLine("curveTangent - x: " + curveTangent.X + " y: " + curveTangent.Y);
Console.WriteLine("curvePrimeTangent - x: " + curvePrimeTangent.X + " y: " + curvePrimeTangent.Y);
}
}
}
Program output:
fail
curvePosition - x: 0.890625 y: 0.703125
curvePrimePosition - x: 0.96875 y: 0.71875
curveTangent - x: 0.8125 y: 1.3125
curvePrimeTangent - x: 0.25 y: 0.375
The short answer is that the math simply does not work the way you want it to.
It's been ages since I toyed around with polynomial curves, so just for fun, I hacked together a Python program that computes control points for a "split" hermite curve, as well as your "wrong" curve. In practice, you might be better off using the de Casteljau algorithm.
This implementation is probably horrible in a gazillion ways, but at least it seems to produce correct results. :-)
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import numpy.polynomial.polynomial as npp
# Hermite basis matrix
B = np.array([[1, 0, -3, 2],
[0, 1, -2, 1],
[0, 0, 3, -2],
[0, 0, -1, 1]])
# Create the parameter space, for plotting. T has one column for each point
# that we plot, with the first row t^0, second row t^1 etc.
step = 0.01
Tt = np.arange(0.0, 1.01, step)
T = np.vstack([np.ones(Tt.shape[0]), Tt, Tt ** 2, Tt ** 3])
# Control points for first curve. One column for each point/tangent.
C1 = np.array([[0, 1, 1, 0],
[0, 0, 1, 1]])
# Coefficients for first curve. @ is matrix multiplication. A1 will be a
# matrix with two rows, one for the "x" polynomial and one for the "y"
# polynomial, and four columns, one for each power in the polynomial.
# So, x(t) = Σ A[0,k] t^k and y(t) = Σ A[1,k] t^k.
A1 = C1 @ B
# Points for first curve.
X1 = A1 @ T
# Parameter value where we will split the curve.
t0 = 0.5
# Evaluate the first curve and its tangent at t=t0. The 'polyval' function
# evaluates a polynomial; 'polyder' computes the derivative of a polynomial.
# Reshape and transpose business because I want my COordinates to be COlumns,
# but polyval/polyder wants coefficients to be columns...
p = npp.polyval(t0, A1.T).reshape(2, 1)
d = npp.polyval(t0, npp.polyder(A1.T, 1)).reshape(2, 1)
# Control points for second, "wrong", curve (last two points are same as C1)
C2 = np.hstack([p, d, C1[:,2:]])
# Coefficients for second, "wrong", curve
A2 = C2 @ B
# Points for second, "wrong", curve
X2 = A2 @ T
# We want to create a new curve, such that that its parameter
# u ∈ [t0, 1] maps to t ∈ [0,1], so we let t = k * u + m.
# To get the curve on [0, t0] instead, set k=t0, m=0.
k = 1 - t0
m = t0
k2 = k * k; k3 = k2 * k; m2 = m * m; m3 = m2 * m
# This matrix maps a polynomial from "t" space to "u" space.
KM = np.array([[1, 0, 0, 0],
[m, k, 0, 0],
[m2, 2 * k * m, k2, 0],
[m3, 3 * k * m2, 3 * k2 * m, k3]])
# A3 are the coefficients for a polynomial which gives the same values
# on the range [0, 1] as the A1 coefficients give on the range [t0, 1].
A3 = A1 @ KM
X3 = A3 @ T
# Compute the control points corresponding to the A3 coefficients, by
# multiplying by the inverse of the B matrix.
C3 = A3 @ np.linalg.inv(B)
# C3 = (np.linalg.solve(B.T, A3.T)).T # Possibly better version!
print(C3)
# Plot curves
fig, ax = plt.subplots()
ax.plot(X1[0,:], X1[1,:], linewidth=3)
ax.plot(X2[0,:], X2[1,:])
ax.plot(X3[0,:], X3[1,:])
ax.set_aspect('equal')
ax.grid()
plt.show()