时间戳

首页 > 解决方案 > Hermite 曲线的子曲线不适合原始曲线

c# - Hermite 曲线的子曲线不适合原始曲线

问题描述

我正在尝试在项目中使用 Hermite 曲线，但对数学的实际工作方式的理解有限，而且我遇到了一些我不理解的行为。我已经用下面的最小代码示例证明了我的困惑，但基本上我希望沿着 Hermite 曲线的子曲线（即使用原始曲线上的点和切线定义的子曲线）来拟合原始曲线，但这似乎是假的。

下面的 c# 代码定义了一个 Hermite 曲线类，该类提供了用于计算沿曲线的某个比率的点的位置和切线的函数。我从互联网上的其他地方复制/粘贴了这两个函数的数学。

然后，一个小型测试工具执行我希望成功但没有成功的测试。我不清楚我的代码中是否存在错误，我的数学中是否存在错误，或者我是否误解了 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

标签： c#splinehermite

解决方案

简短的回答是，数学根本无法按照您想要的方式工作。

自从我玩弄多项式曲线以来已经有很多年了，所以只是为了好玩，我编写了一个 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()

推荐阅读