Catmull-rom曲线没有尖点,没有自交叉

Question

我有以下代码来计算四个控制点之间的点,以生成catmull-rom曲线:

CGPoint interpolatedPosition(CGPoint p0, CGPoint p1, CGPoint p2, CGPoint p3, float t)
{
    float t3 = t * t * t;
    float t2 = t * t;

    float f1 = -0.5 * t3 + t2 - 0.5 * t;
    float f2 = 1.5 * t3 - 2.5 * t2 + 1.0;
    float f3 = -1.5 * t3 + 2.0 * t2 + 0.5 * t;
    float f4 = 0.5 * t3 - 0.5 * t2;

    float x = p0.x * f1 + p1.x * f2 + p2.x * f3 + p3.x * f4;
    float y = p0.y * f1 + p1.y * f2 + p2.y * f3 + p3.y * f4;

    return CGPointMake(x, y);
}

这工作正常,但我想创建一些我认为称为向心参数化的东西.这意味着曲线将没有尖点并且没有自交叉.如果我将一个控制点移动到另一个控制点附近,则曲线应该变得"更小".我试着找到一种方法来搜索我的眼睛.有人知道怎么做吗？

Answer 1

我也需要为工作实现这一点.您需要从的基本概念是,常规Catmull-Rom实现与修改版本之间的主要区别在于它们如何处理时间.

在原始Catmull-Rom实现的非参数化版本中,t从0开始,以1结束,并计算从P1到P2的曲线.在参数化时间实现中,t在P0处从0开始,并且在所有四个点上保持增加.因此在统一的情况下,它在P1处为1,在P2处为2,并且您将为插值传递1到2的值.

和弦情况显示| Pi + 1 - P | 随着时间跨度的变化.这只是意味着您可以使用每个段的点之间的直线距离来计算要使用的实际长度.向心情况仅使用稍微不同的方法来计算每个段使用的最佳时间长度.

所以现在我们只需要知道如何设计方程式,让我们插入新的时间值.典型的Catmull-Rom方程只有一个t,即您尝试计算值的时间.我找到了描述如何计算这些参数的最佳文章:http://www.cemyuksel.com/research/catmullrom_param/catmullrom.pdf.他们专注于曲线的数学评估,但其中有一个来自Barry和Goldman的关键公式.(1)

在上图中,箭头表示"乘以"箭头中给出的比率.

然后,这为我们提供了实际执行计算以获得所需结果所需的内容.X和Y是独立计算的,尽管我使用"距离"因子来根据2D距离修改时间,而不是1D距离.

检测结果:

(1)PJ Barry和RN Goldman.一类catmull-rom样条的递归评估算法.SIGGRAPH Computer Graphics,22(4):199 {204,1988.

我在Java中最终实现的源代码如下所示:

/**
 * This method will calculate the Catmull-Rom interpolation curve, returning
 * it as a list of Coord coordinate objects.  This method in particular
 * adds the first and last control points which are not visible, but required
 * for calculating the spline.
 *
 * @param coordinates The list of original straight line points to calculate
 * an interpolation from.
 * @param pointsPerSegment The integer number of equally spaced points to
 * return along each curve.  The actual distance between each
 * point will depend on the spacing between the control points.
 * @return The list of interpolated coordinates.
 * @param curveType Chordal (stiff), Uniform(floppy), or Centripetal(medium)
 * @throws gov.ca.water.shapelite.analysis.CatmullRomException if
 * pointsPerSegment is less than 2.
 */
public static List<Coord> interpolate(List<Coord> coordinates, int pointsPerSegment, CatmullRomType curveType)
        throws CatmullRomException {
    List<Coord> vertices = new ArrayList<>();
    for (Coord c : coordinates) {
        vertices.add(c.copy());
    }
    if (pointsPerSegment < 2) {
        throw new CatmullRomException("The pointsPerSegment parameter must be greater than 2, since 2 points is just the linear segment.");
    }

    // Cannot interpolate curves given only two points.  Two points
    // is best represented as a simple line segment.
    if (vertices.size() < 3) {
        return vertices;
    }

    // Test whether the shape is open or closed by checking to see if
    // the first point intersects with the last point.  M and Z are ignored.
    boolean isClosed = vertices.get(0).intersects2D(vertices.get(vertices.size() - 1));
    if (isClosed) {
        // Use the second and second from last points as control points.
        // get the second point.
        Coord p2 = vertices.get(1).copy();
        // get the point before the last point
        Coord pn1 = vertices.get(vertices.size() - 2).copy();

        // insert the second from the last point as the first point in the list
        // because when the shape is closed it keeps wrapping around to
        // the second point.
        vertices.add(0, pn1);
        // add the second point to the end.
        vertices.add(p2);
    } else {
        // The shape is open, so use control points that simply extend
        // the first and last segments

        // Get the change in x and y between the first and second coordinates.
        double dx = vertices.get(1).X - vertices.get(0).X;
        double dy = vertices.get(1).Y - vertices.get(0).Y;

        // Then using the change, extrapolate backwards to find a control point.
        double x1 = vertices.get(0).X - dx;
        double y1 = vertices.get(0).Y - dy;

        // Actaully create the start point from the extrapolated values.
        Coord start = new Coord(x1, y1, vertices.get(0).Z);

        // Repeat for the end control point.
        int n = vertices.size() - 1;
        dx = vertices.get(n).X - vertices.get(n - 1).X;
        dy = vertices.get(n).Y - vertices.get(n - 1).Y;
        double xn = vertices.get(n).X + dx;
        double yn = vertices.get(n).Y + dy;
        Coord end = new Coord(xn, yn, vertices.get(n).Z);

        // insert the start control point at the start of the vertices list.
        vertices.add(0, start);

        // append the end control ponit to the end of the vertices list.
        vertices.add(end);
    }

    // Dimension a result list of coordinates. 
    List<Coord> result = new ArrayList<>();
    // When looping, remember that each cycle requires 4 points, starting
    // with i and ending with i+3.  So we don't loop through all the points.
    for (int i = 0; i < vertices.size() - 3; i++) {

        // Actually calculate the Catmull-Rom curve for one segment.
        List<Coord> points = interpolate(vertices, i, pointsPerSegment, curveType);
        // Since the middle points are added twice, once for each bordering
        // segment, we only add the 0 index result point for the first
        // segment.  Otherwise we will have duplicate points.
        if (result.size() > 0) {
            points.remove(0);
        }

        // Add the coordinates for the segment to the result list.
        result.addAll(points);
    }
    return result;

}

/**
 * Given a list of control points, this will create a list of pointsPerSegment
 * points spaced uniformly along the resulting Catmull-Rom curve.
 *
 * @param points The list of control points, leading and ending with a 
 * coordinate that is only used for controling the spline and is not visualized.
 * @param index The index of control point p0, where p0, p1, p2, and p3 are
 * used in order to create a curve between p1 and p2.
 * @param pointsPerSegment The total number of uniformly spaced interpolated
 * points to calculate for each segment. The larger this number, the
 * smoother the resulting curve.
 * @param curveType Clarifies whether the curve should use uniform, chordal
 * or centripetal curve types. Uniform can produce loops, chordal can
 * produce large distortions from the original lines, and centripetal is an
 * optimal balance without spaces.
 * @return the list of coordinates that define the CatmullRom curve
 * between the points defined by index+1 and index+2.
 */
public static List<Coord> interpolate(List<Coord> points, int index, int pointsPerSegment, CatmullRomType curveType) {
    List<Coord> result = new ArrayList<>();
    double[] x = new double[4];
    double[] y = new double[4];
    double[] time = new double[4];
    for (int i = 0; i < 4; i++) {
        x[i] = points.get(index + i).X;
        y[i] = points.get(index + i).Y;
        time[i] = i;
    }

    double tstart = 1;
    double tend = 2;
    if (!curveType.equals(CatmullRomType.Uniform)) {
        double total = 0;
        for (int i = 1; i < 4; i++) {
            double dx = x[i] - x[i - 1];
            double dy = y[i] - y[i - 1];
            if (curveType.equals(CatmullRomType.Centripetal)) {
                total += Math.pow(dx * dx + dy * dy, .25);
            } else {
                total += Math.pow(dx * dx + dy * dy, .5);
            }
            time[i] = total;
        }
        tstart = time[1];
        tend = time[2];
    }
    double z1 = 0.0;
    double z2 = 0.0;
    if (!Double.isNaN(points.get(index + 1).Z)) {
        z1 = points.get(index + 1).Z;
    }
    if (!Double.isNaN(points.get(index + 2).Z)) {
        z2 = points.get(index + 2).Z;
    }
    double dz = z2 - z1;
    int segments = pointsPerSegment - 1;
    result.add(points.get(index + 1));
    for (int i = 1; i < segments; i++) {
        double xi = interpolate(x, time, tstart + (i * (tend - tstart)) / segments);
        double yi = interpolate(y, time, tstart + (i * (tend - tstart)) / segments);
        double zi = z1 + (dz * i) / segments;
        result.add(new Coord(xi, yi, zi));
    }
    result.add(points.get(index + 2));
    return result;
}

/**
 * Unlike the other implementation here, which uses the default "uniform"
 * treatment of t, this computation is used to calculate the same values but
 * introduces the ability to "parameterize" the t values used in the
 * calculation. This is based on Figure 3 from
 * http://www.cemyuksel.com/research/catmullrom_param/catmullrom.pdf
 *
 * @param p An array of double values of length 4, where interpolation
 * occurs from p1 to p2.
 * @param time An array of time measures of length 4, corresponding to each
 * p value.
 * @param t the actual interpolation ratio from 0 to 1 representing the
 * position between p1 and p2 to interpolate the value.
 * @return
 */
public static double interpolate(double[] p, double[] time, double t) {
    double L01 = p[0] * (time[1] - t) / (time[1] - time[0]) + p[1] * (t - time[0]) / (time[1] - time[0]);
    double L12 = p[1] * (time[2] - t) / (time[2] - time[1]) + p[2] * (t - time[1]) / (time[2] - time[1]);
    double L23 = p[2] * (time[3] - t) / (time[3] - time[2]) + p[3] * (t - time[2]) / (time[3] - time[2]);
    double L012 = L01 * (time[2] - t) / (time[2] - time[0]) + L12 * (t - time[0]) / (time[2] - time[0]);
    double L123 = L12 * (time[3] - t) / (time[3] - time[1]) + L23 * (t - time[1]) / (time[3] - time[1]);
    double C12 = L012 * (time[2] - t) / (time[2] - time[1]) + L123 * (t - time[1]) / (time[2] - time[1]);
    return C12;
}   

Answer 2

有一种更简单,更有效的方法来实现它,只需要您使用不同的公式计算切线,而无需实现Barry和Goldman的递归评估算法.

如果您采用Barry-Goldman参数化(在Ted的答案中引用)C(t)作为节点(t0,t1,t2,t3)和控制点(P0,P1,P2,P3),其闭合形式相当复杂,但是当你将它约束到区间(t1,t2)时,它仍然是t中的三次多项式.所以我们需要完全描述的是两个端点t1和t2的值和切线.如果我们计算出这些值(我在Mathematica中做到了这一点),我们发现

C(t1)  = P1
C(t2)  = P2
C'(t1) = (P1 - P0) / (t1 - t0) - (P2 - P0) / (t2 - t0) + (P2 - P1) / (t2 - t1)
C'(t2) = (P2 - P1) / (t2 - t1) - (P3 - P1) / (t3 - t1) + (P3 - P2) / (t3 - t2)

我们可以简单地将其插入到标准公式中,用于计算在端点处具有给定值和切线的三次样条,并且我们具有非均匀的Catmull-Rom样条.需要注意的是,上面的切线是针对区间(t1,t2)计算的,因此如果要在标准区间(0,1)中评估曲线,只需通过将切线与系数(t2-t1)相乘来重新缩放切线. ).

我在Ideone上放了一个有效的C++示例:http://ideone.com/NoEbVM

我还会粘贴下面的代码.

#include <iostream>
#include <cmath>

using namespace std;

struct CubicPoly
{
    float c0, c1, c2, c3;

    float eval(float t)
    {
        float t2 = t*t;
        float t3 = t2 * t;
        return c0 + c1*t + c2*t2 + c3*t3;
    }
};

/*
 * Compute coefficients for a cubic polynomial
 *   p(s) = c0 + c1*s + c2*s^2 + c3*s^3
 * such that
 *   p(0) = x0, p(1) = x1
 *  and
 *   p'(0) = t0, p'(1) = t1.
 */
void InitCubicPoly(float x0, float x1, float t0, float t1, CubicPoly &p)
{
    p.c0 = x0;
    p.c1 = t0;
    p.c2 = -3*x0 + 3*x1 - 2*t0 - t1;
    p.c3 = 2*x0 - 2*x1 + t0 + t1;
}

// standard Catmull-Rom spline: interpolate between x1 and x2 with previous/following points x0/x3
// (we don't need this here, but it's for illustration)
void InitCatmullRom(float x0, float x1, float x2, float x3, CubicPoly &p)
{
    // Catmull-Rom with tension 0.5
    InitCubicPoly(x1, x2, 0.5f*(x2-x0), 0.5f*(x3-x1), p);
}

// compute coefficients for a nonuniform Catmull-Rom spline
void InitNonuniformCatmullRom(float x0, float x1, float x2, float x3, float dt0, float dt1, float dt2, CubicPoly &p)
{
    // compute tangents when parameterized in [t1,t2]
    float t1 = (x1 - x0) / dt0 - (x2 - x0) / (dt0 + dt1) + (x2 - x1) / dt1;
    float t2 = (x2 - x1) / dt1 - (x3 - x1) / (dt1 + dt2) + (x3 - x2) / dt2;

    // rescale tangents for parametrization in [0,1]
    t1 *= dt1;
    t2 *= dt1;

    InitCubicPoly(x1, x2, t1, t2, p);
}

struct Vec2D
{
    Vec2D(float _x, float _y) : x(_x), y(_y) {}
    float x, y;
};

float VecDistSquared(const Vec2D& p, const Vec2D& q)
{
    float dx = q.x - p.x;
    float dy = q.y - p.y;
    return dx*dx + dy*dy;
}

void InitCentripetalCR(const Vec2D& p0, const Vec2D& p1, const Vec2D& p2, const Vec2D& p3,
    CubicPoly &px, CubicPoly &py)
{
    float dt0 = powf(VecDistSquared(p0, p1), 0.25f);
    float dt1 = powf(VecDistSquared(p1, p2), 0.25f);
    float dt2 = powf(VecDistSquared(p2, p3), 0.25f);

    // safety check for repeated points
    if (dt1 < 1e-4f)    dt1 = 1.0f;
    if (dt0 < 1e-4f)    dt0 = dt1;
    if (dt2 < 1e-4f)    dt2 = dt1;

    InitNonuniformCatmullRom(p0.x, p1.x, p2.x, p3.x, dt0, dt1, dt2, px);
    InitNonuniformCatmullRom(p0.y, p1.y, p2.y, p3.y, dt0, dt1, dt2, py);
}


int main()
{
    Vec2D p0(0,0), p1(1,1), p2(1.1,1), p3(2,0);
    CubicPoly px, py;
    InitCentripetalCR(p0, p1, p2, p3, px, py);
    for (int i = 0; i <= 10; ++i)
        cout << px.eval(0.1f*i) << " " << py.eval(0.1f*i) << endl;
}