Courbe Catmull-rom sans cusps et sans auto-intersections

il y a une façon beaucoup plus simple et efficace de mettre en œuvre ce qui ne nécessite que vous pour calculer vos tangentes en utilisant une formule différente, sans avoir à mettre en œuvre l'algorithme d'évaluation récursive de Barry et Goldman.

si vous prenez la paramétrisation Barry-Goldman (référencée dans la réponse de Ted) C( t) pour les noeuds (t0, t1, t2, t3) et les points de contrôle( P0,P1, P2, P3), sa forme fermée est assez compliquée, mais à la fin c'est toujours un polynôme cubique en t quand vous limitez-le à l'intervalle (t1,t2). Tout ce qu'il nous faut donc pour le décrire pleinement, ce sont les valeurs et les tangentes aux deux extrémités t1 et t2. Si nous calculons ces valeurs (Je l'ai fait dans Mathematica), nous trouvons

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)

nous pouvons simplement brancher ceci dans le formule standard pour le calcul d'une spline cubique avec des valeurs données et tangentes aux extrémités et nous avons notre spline Catmull-Rom non uniforme. Une mise en garde est que les tangentes ci-dessus sont calculées pour l'intervalle (t1,t2), donc si vous voulez évaluer la courbe dans l'intervalle standard (0,1), il suffit de redimensionner les tangentes en les multipliant par le facteur (t2-t1).

j'ai mis un exemple C++ sur Ideone: http://ideone.com/NoEbVM

je vais aussi coller le code ci-dessous.

#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;
}

Voici une version iOS du code de Ted. J'ai exclu les parties en "z".

.h

typedef enum {
    CatmullRomTypeUniform,
    CatmullRomTypeChordal,
    CatmullRomTypeCentripetal
} CatmullRomType ;

.m

-(NSMutableArray *)interpolate:(NSArray *)coordinates withPointsPerSegment:(NSInteger)pointsPerSegment andType:(CatmullRomType)curveType {

    NSMutableArray *vertices = [[NSMutableArray alloc] initWithArray:coordinates copyItems:YES];

    if (pointsPerSegment < 3)
        return vertices;

    //start point
    CGPoint pt1 = [vertices[0] CGPointValue];
    CGPoint pt2 = [vertices[1] CGPointValue];

    double dx = pt2.x - pt1.x;
    double dy = pt2.y - pt1.y;

    double x1 = pt1.x - dx;
    double y1 = pt1.y - dy;

    CGPoint start = CGPointMake(x1*.5, y1);

    //end point
    pt2 = [vertices[vertices.count-1] CGPointValue];
    pt1 = [vertices[vertices.count-2] CGPointValue];

    dx = pt2.x - pt1.x;
    dy = pt2.y - pt1.y;

    x1 = pt2.x + dx;
    y1 = pt2.y + dy;

    CGPoint end = CGPointMake(x1, y1);

    [vertices insertObject:[NSValue valueWithCGPoint:start] atIndex:0];
    [vertices addObject:[NSValue valueWithCGPoint:end]];

    NSMutableArray *result = [[NSMutableArray alloc] init];

    for (int i = 0; i < vertices.count - 3; i++) {
        NSMutableArray *points = [self interpolate:vertices forIndex:i withPointsPerSegment:pointsPerSegment andType:curveType];

        if ([points count] > 0)
            [points removeObjectAtIndex:0];

        [result addObjectsFromArray:points];
    }

    return result;
}

-(double)interpolate:(double*)p  time:(double*)time t:(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;
    }

-(NSMutableArray*)interpolate:(NSArray *)points forIndex:(NSInteger)index withPointsPerSegment:(NSInteger)pointsPerSegment andType:(CatmullRomType)curveType {
    NSMutableArray *result = [[NSMutableArray alloc] init];

    double x[4];
    double y[4];
    double time[4];

    for (int i=0; i < 4; i++) {
        x[i] = [points[index+i] CGPointValue].x;
        y[i] = [points[index+i] CGPointValue].y;
        time[i] = i;
    }

    double tstart = 1;
    double tend = 2;

    if (curveType != CatmullRomTypeUniform) {
        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 == CatmullRomTypeCentripetal) {
                total += pow(dx * dx + dy * dy, 0.25);
            }
            else {
                total += pow(dx * dx + dy * dy, 0.5); //sqrt
            }
            time[i] = total;
        }
        tstart = time[1];
        tend = time[2];
    }

    int segments = pointsPerSegment - 1;

    [result addObject:points[index+1]];

    for (int i =1; i < segments; i++) {
        double xi = [self interpolate:x time:time t:tstart + (i * (tend - tstart)) / segments];
        double yi = [self interpolate:y time:time t:tstart + (i * (tend - tstart)) / segments];
        NSLog(@"(%f,%f)",xi,yi);
        [result addObject:[NSValue valueWithCGPoint:CGPointMake(xi, yi)]];
    }
    [result addObject:points[index+2]];

    return result;
}

aussi, Voici une méthode pour transformer un tableau de points en un chemin plus Bézier pour le dessin, en utilisant le

-(UIBezierPath*)bezierPathFromPoints:(NSArray *)points withGranulaity:(NSInteger)granularity
{
    UIBezierPath __block *path = [[UIBezierPath alloc] init];

    NSMutableArray *curve = [self interpolate:points withPointsPerSegment:granularity andType:CatmullRomTypeCentripetal];

    CGPoint __block p0 = [curve[0] CGPointValue];

    [path moveToPoint:p0];

    //use this loop to draw lines between all points
    for (int idx=1; idx < [curve count]; idx+=1) {
        CGPoint c1 = [curve[idx] CGPointValue];

        [path addLineToPoint:c1];
    };

    //or use this loop to use actual control points (less smooth but probably faster)
//    for (int idx=0; idx < [curve count]-3; idx+=3) {
//        CGPoint c1 = [curve[idx+1] CGPointValue];
//        CGPoint c2 = [curve[idx+2] CGPointValue];
//        CGPoint p1 = [curve[idx+3] CGPointValue];
//
//        [path addCurveToPoint:p1 controlPoint1:c1 controlPoint2:c2];
//    };

    return path;
}

j'ai codé quelque chose en Python (adaptation de la page de Wikipedia Catmull-Rom) qui compare les lignes uniformes, centripèdes, et chorales (bien que vous pouvez définir alpha à tout ce que vous voulez) en utilisant des données aléatoires (vous pouvez utiliser vos propres données et les fucntions fonctionnent bien). À noter que pour les points de terminaison j'ai juste coincé dans un rapide "hack" qui maintient la pente de la première et de la dernière 2 points, bien que la distance entre ce point et le premier/perdu point connu est arbitraire (je l'ai mis à 1% de la domaine... pour aucune raison du tout. Alors gardez cela à l'esprit avant de vous appliquer à quelque chose d'important):

# coding: utf-8

# In[1]:

import numpy
import matplotlib.pyplot as plt
get_ipython().magic(u'pylab inline')


# In[2]:

def CatmullRomSpline(P0, P1, P2, P3, a, nPoints=100):
  """
  P0, P1, P2, and P3 should be (x,y) point pairs that define the Catmull-Rom spline.
  nPoints is the number of points to include in this curve segment.
  """
  # Convert the points to numpy so that we can do array multiplication
  P0, P1, P2, P3 = map(numpy.array, [P0, P1, P2, P3])

  # Calculate t0 to t4
  alpha = a
  def tj(ti, Pi, Pj):
    xi, yi = Pi
    xj, yj = Pj
    return ( ( (xj-xi)**2 + (yj-yi)**2 )**0.5 )**alpha + ti

  t0 = 0
  t1 = tj(t0, P0, P1)
  t2 = tj(t1, P1, P2)
  t3 = tj(t2, P2, P3)

  # Only calculate points between P1 and P2
  t = numpy.linspace(t1,t2,nPoints)

  # Reshape so that we can multiply by the points P0 to P3
  # and get a point for each value of t.
  t = t.reshape(len(t),1)

  A1 = (t1-t)/(t1-t0)*P0 + (t-t0)/(t1-t0)*P1
  A2 = (t2-t)/(t2-t1)*P1 + (t-t1)/(t2-t1)*P2
  A3 = (t3-t)/(t3-t2)*P2 + (t-t2)/(t3-t2)*P3

  B1 = (t2-t)/(t2-t0)*A1 + (t-t0)/(t2-t0)*A2
  B2 = (t3-t)/(t3-t1)*A2 + (t-t1)/(t3-t1)*A3

  C  = (t2-t)/(t2-t1)*B1 + (t-t1)/(t2-t1)*B2
  return C

def CatmullRomChain(P,alpha):
  """
  Calculate Catmull Rom for a chain of points and return the combined curve.
  """
  sz = len(P)

  # The curve C will contain an array of (x,y) points.
  C = []
  for i in range(sz-3):
    c = CatmullRomSpline(P[i], P[i+1], P[i+2], P[i+3],alpha)
    C.extend(c)

  return C


# In[8]:

# Define a set of points for curve to go through
Points = numpy.random.rand(12,2)
x1=Points[0][0]
x2=Points[1][0]
y1=Points[0][1]
y2=Points[1][1]
x3=Points[-2][0]
x4=Points[-1][0]
y3=Points[-2][1]
y4=Points[-1][1]
dom=max(Points[:,0])-min(Points[:,0])
rng=max(Points[:,1])-min(Points[:,1])
prex=x1+sign(x1-x2)*dom*0.01
prey=(y1-y2)/(x1-x2)*dom*0.01+y1
endx=x4+sign(x4-x3)*dom*0.01
endy=(y4-y3)/(x4-x3)*dom*0.01+y4
print len(Points)
Points=list(Points)
Points.insert(0,array([prex,prey]))
Points.append(array([endx,endy]))
print len(Points)


# In[9]:

#Define alpha
a=0.

# Calculate the Catmull-Rom splines through the points
c = CatmullRomChain(Points,a)

# Convert the Catmull-Rom curve points into x and y arrays and plot
x,y = zip(*c)
plt.plot(x,y,c='green',zorder=10)

# Plot the control points
px, py = zip(*Points)
plt.plot(px,py,'or')

a=0.5
c = CatmullRomChain(Points,a)
x,y = zip(*c)
plt.plot(x,y,c='blue')

a=1.
c = CatmullRomChain(Points,a)
x,y = zip(*c)
plt.plot(x,y,c='red')


plt.grid(b=True)
plt.show()


# In[10]:

Points


# In[ ]:

code original:https://en.wikipedia.org/wiki/Centripetal_Catmull%E2%80%93Rom_spline