Sub-curve of a Hermite Curve does not fit original curve

Question

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

Answer 1

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()

Answer 2

you're using pointer equality instead of object equality?

Answer 3

The problem in your codes is that you are taking the first derivative vectors from original Hermite curve and use them directly to create the sub-curve. The sub-curve created this way will not actually be the same as the [0.5, 1] portion of the original Hermite curve. You can try to draw the curves and you will see that your original curve and the sub-curve do not match.

To make the sub-curve match the [0.5, 1.0] portion of the original curve, you will have to use "m0prime/2" and "m1/2" to create your "curvePrime". Once you do that, you will find that the computed "curvePoint" and "curvePrimePoint" will be identical and the "curvePrimeTangent" will be half of "curveTangent".

In a nutshell, cubic Hermite curve and cubic Bezier curve are essentially the same thing (cubic polynomial curve) but with different representation. Cubic Hermite curve's definition contains points and vectors while cubic Bezier curve's definition contains points only. So, indeed cubic Hermite curve will sometime cause confusions. But cubic Bezier curves have its downside as well as the two middle control points do not lie on the curve, which makes the "point interpolation" problem harder to solve when using cubic Bezier curve.