Hermite 曲线的子曲线不适合原始曲线
Sub-curve of a Hermite Curve does not fit original curve
我正在尝试在一个项目中使用 Hermite 曲线，但我对数学的实际运作方式的理解有限，而且我 运行 遇到了一些我不理解的行为。我已经用下面的最小代码示例证明了我的困惑，但基本上我希望沿着埃尔米特曲线的子曲线（即使用原始曲线上的点和切线定义的子曲线）来拟合原始曲线，但这似乎是假的。
以下 c# 代码定义了一条 Hermite 曲线 class，该曲线提供了用于计算沿曲线以某种比率的点的位置和切线的函数。我 copy/pasted 来自互联网上其他地方的两个函数的数学运算。
然后，一个小型测试工具执行我期望成功但没有成功的测试。我不清楚我的代码中是否存在错误、我的数学错误，或者我是否误解了 Hermite 曲线的工作原理以及该测试实际上不应通过。
如有任何见解，我们将不胜感激。
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);
}
}
}
程序输出：
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
您使用的是指针相等而不是对象相等？
简短的回答是，数学根本无法按照您想要的方式运行。
我已经很久没有玩弄多项式曲线了，所以为了好玩，我拼凑了一个 Python 程序来计算“分裂”厄米曲线的控制点，以及你的“错误”曲线。实际上，使用 de Casteljau 算法可能会更好。
这个实现可能在很多方面都很糟糕，但至少它似乎产生了正确的结果。 :-)
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()
您代码中的问题是您从原始 Hermite 曲线中获取一阶导数向量并直接使用它们来创建 sub-curve。这样创建的 sub-curve 实际上不会与原始 Hermite 曲线的 [0.5, 1] 部分相同。您可以尝试绘制曲线，您会发现您的原始曲线与 sub-curve 不匹配。
要使 sub-curve 匹配原始曲线的 [0.5, 1.0] 部分，您必须使用“m0prime/2”和“m1/2”来创建您的“曲线”。执行此操作后，您会发现计算出的“curvePoint”和“curvePrimePoint”将相同，并且“curvePrimeTangent”将是“curveTangent”的一半。
简而言之，三次Hermite曲线和三次贝塞尔曲线本质上是同一个东西（三次多项式曲线），只是表现形式不同。三次 Hermite 曲线的定义包含点和向量，而三次贝塞尔曲线的定义仅包含点。所以，三次 Hermite 曲线确实有时会引起混淆。但是三次贝塞尔曲线也有它的缺点，中间两个控制点不在曲线上，这使得使用三次贝塞尔曲线时“点插值”问题更难解决。
我正在尝试在一个项目中使用 Hermite 曲线，但我对数学的实际运作方式的理解有限，而且我 运行 遇到了一些我不理解的行为。我已经用下面的最小代码示例证明了我的困惑，但基本上我希望沿着埃尔米特曲线的子曲线（即使用原始曲线上的点和切线定义的子曲线）来拟合原始曲线，但这似乎是假的。
以下 c# 代码定义了一条 Hermite 曲线 class，该曲线提供了用于计算沿曲线以某种比率的点的位置和切线的函数。我 copy/pasted 来自互联网上其他地方的两个函数的数学运算。
然后，一个小型测试工具执行我期望成功但没有成功的测试。我不清楚我的代码中是否存在错误、我的数学错误，或者我是否误解了 Hermite 曲线的工作原理以及该测试实际上不应通过。
如有任何见解，我们将不胜感激。
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);
}
}
}
程序输出：
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
您使用的是指针相等而不是对象相等？
简短的回答是，数学根本无法按照您想要的方式运行。
我已经很久没有玩弄多项式曲线了，所以为了好玩，我拼凑了一个 Python 程序来计算“分裂”厄米曲线的控制点，以及你的“错误”曲线。实际上，使用 de Casteljau 算法可能会更好。
这个实现可能在很多方面都很糟糕，但至少它似乎产生了正确的结果。 :-)
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()
您代码中的问题是您从原始 Hermite 曲线中获取一阶导数向量并直接使用它们来创建 sub-curve。这样创建的 sub-curve 实际上不会与原始 Hermite 曲线的 [0.5, 1] 部分相同。您可以尝试绘制曲线，您会发现您的原始曲线与 sub-curve 不匹配。
要使 sub-curve 匹配原始曲线的 [0.5, 1.0] 部分，您必须使用“m0prime/2”和“m1/2”来创建您的“曲线”。执行此操作后，您会发现计算出的“curvePoint”和“curvePrimePoint”将相同，并且“curvePrimeTangent”将是“curveTangent”的一半。
简而言之，三次Hermite曲线和三次贝塞尔曲线本质上是同一个东西（三次多项式曲线），只是表现形式不同。三次 Hermite 曲线的定义包含点和向量，而三次贝塞尔曲线的定义仅包含点。所以，三次 Hermite 曲线确实有时会引起混淆。但是三次贝塞尔曲线也有它的缺点，中间两个控制点不在曲线上，这使得使用三次贝塞尔曲线时“点插值”问题更难解决。