有限點位擬合平滑曲線演算法：二次拋物線、B樣條與三次參數樣條

給定有限離散點，如何繪製平滑曲線？本文提供三種經典演算法的 C 實作：二次拋物線法、二次 B 樣條法、三次參數樣條法。適用於儀表繪圖、軌跡插值等場景。

共同資料結構

// 點陣列：x[], y[] 長度 n+2（含首尾各擴展 1 點用於邊界處理）
// k：每段插值點數（曲線細緻度）

1. 二次拋物線法

每 4 點構造一段二次拋物線，確保連續性。

void paowuxian(int *x, int *y, int n, unsigned int k) {
    unsigned int i, j;
    float t1, t2, t3, t, a, b, c, d, tx, ty;

    x[0] = x[1]; y[0] = y[1];           // 擴展首點
    x[n+1] = x[n]; y[n+1] = y[n];       // 擴展尾點

    t = 0.5f / k;
    setcolor(10);
    moveto(x[1], y[1]);

    for (i = 0; i < n-1; i++) {
        for (j = 1; j < k; j++) {
            t1 = j * t;
            t2 = t1 * t1;
            t3 = t2 * t1;
            a = 4*t2 - t1 - 4*t3;
            b = 1 - 10*t2 + 12*t3;
            c = t1 + 8*t2 - 12*t3;
            d = 4*t3 - 2*t2;
            tx = a*x[i] + b*x[i+1] + c*x[i+2] + d*x[i+3];
            ty = a*y[i] + b*y[i+1] + c*y[i+2] + d*y[i+3];
            lineto(tx, ty);
        }
    }
    lineto(x[i+2], y[i+2]);
}

2. 二次 B 樣條法

局部控制特性佳，曲線不經過控制點，但平滑度高。

void byangtiao(int *x, int *y, int n, unsigned int k) {
    unsigned int i, j;
    float t, t1, t2, a, b, c, tx, ty;

    x[0] = x[1]; y[0] = y[1];
    x[n+1] = x[n]; y[n+1] = y[n];

    t = 1.0f / k;
    setcolor(13);
    moveto((x[0]+x[1])/2.0f, (y[0]+y[1])/2.0f);

    for (i = 0; i < n; i++) {
        for (j = 1; j < k; j++) {
            t1 = j * t;
            t2 = t1 * t1;
            a = (t2 - 2*t1 + 1) / 2.0f;
            b = t1 - t2 + 0.5f;
            c = t2 / 2.0f;
            tx = a*x[i] + b*x[i+1] + c*x[i+2];
            ty = a*y[i] + b*y[i+1] + c*y[i+2];
            lineto(tx, ty);
        }
    }
}

3. 三次參數樣條法

三次多項式插值，曲線經過所有控制點，一階、二階導數連續。

void parspl(int p[][2], int n, int k) {
    int i, j;
    float t1, t2, t3, t, a, b, c, d, x, y;

    p[0][0] = p[1][0]; p[0][1] = p[1][1];
    p[n+1][0] = p[n][0]; p[n+1][1] = p[n][1];

    t = 0.5f / k;
    moveto(p[1][0], p[1][1]);

    for (i = 0; i < n-1; i++) {
        for (j = 1; j < k; j++) {
            t1 = j * t;
            t2 = t1 * t1;
            t3 = t2 * t1;
            a = 4*t2 - t1 - 4*t3;
            b = 1 - 10*t2 + 12*t3;
            c = t1 + 8*t2 - 12*t3;
            d = 4*t3 - 2*t2;
            x = a*p[i][0] + b*p[i+1][0] + c*p[i+2][0] + d*p[i+3][0];
            y = a*p[i][1] + b*p[i+1][1] + c*p[i+2][1] + d*p[i+3][1];
            lineto(x, y);
        }
    }
    lineto(p[i+2][0], p[i+2][1]);
}

演算法選擇指南

現代化建議

• 使用 std::vector<Point> 取代原始陣列。
• 封裝為 CurveFitter 類別，支援多種插值策略（策略模式）。
• 搭配 OpenGL/Direct2D/GDI+ 硬體加速繪製。
• 大量點位時考慮 Catmull-Rom 樣條或 NURBS。