分段三階貝塞爾曲線的光滑連續性條件

1. 貝塞爾曲線簡介 [1]

計算機圖形學中有一類很常用的曲線，俗稱貝塞爾曲線。1962年，法國數學家Pierre Bézier 第一個研究了這種矢量繪制曲線的方法，并給出了詳細的計算公式，是以按照這樣的公式繪制出來的曲線就用他的姓氏來命名是為貝塞爾曲線。很多程式語言都有實作貝塞爾曲線的API，而該曲線本身也擁有強大的近似其它曲線的能力，即使一條不能夠勝任，那麼分段的多條貝塞爾曲線也足夠用來近似我們想繪制的曲線。

這裡給出貝塞爾曲線的一般公式：

P ( t ) = ∑ i = 0 n B i , n ( t ) P i P(t)=\sum_{i=0}^n{B_{i,n}(t)\mathbf{P_i}} P(t)=i=0∑nBi,n(t)Pi

其中：

B i , n ( t ) = n ! i ! ( n − i ) ! t i ( 1 − t ) n − i , t ∈ [ 0 , 1 ] B_{i,n}(t)=\frac{n!}{i!(n-i)!}t^i(1-t)^{n-i},t\in[0,1] Bi,n(t)=i!(n−i)!n!ti(1−t)n−i,t∈[0,1]

對于以$ P_{i-1} 和和和 P_{i} 作為起點和終點的三階貝塞爾曲線段，具有兩個控制點作為起點和終點的三階貝塞爾曲線段，具有兩個控制點作為起點和終點的三階貝塞爾曲線段，具有兩個控制點 C_{i,1} 、、、 C_{i,2} $，方程為：

P i − 1 , i ( t ) = P i − 1 + ( − 3 P i − 1 + 3 C i , 1 ) t + ( 3 P i − 1 − 6 C i , 1 + 3 C i , 2 ) t 2 + ( − P i − 1 + 3 C i , 1 − 3 C i , 2 + P i ) t 3 P_{i-1,i}(t)= P_{i-1} +(- 3P_{i-1} + 3C_{i,1} )t +( 3P_{i-1} - 6C_{i,1} + 3C_{i,2} )t^2 +(- P_{i-1} + 3C_{i,1} - 3C_{i,2} + P_{i})t^3 Pi−1,i(t)=Pi−1+(−3Pi−1+3Ci,1)t+(3Pi−1−6Ci,1+3Ci,2)t2+(−Pi−1+3Ci,1−3Ci,2+Pi)t3

對應的方程一階導數即為：

d P i − 1 , i ( t ) d t = ( − 3 P i − 1 + 3 C i , 1 ) + ( 6 P i − 1 − 12 C i , 1 + 6 C i , 2 ) t + ( − 3 P i − 1 + 9 C i , 1 − 9 C i , 2 + 3 P i ) t 2 \frac{dP_{i-1,i}(t)}{dt} = (- 3P_{i-1} + 3C_{i,1} ) +( 6P_{i-1} -12C_{i,1} + 6C_{i,2} )t +(- 3P_{i-1} + 9C_{i,1} - 9C_{i,2} + 3P_{i})t^2 dtdPi−1,i(t)=(−3Pi−1+3Ci,1)+(6Pi−1−12Ci,1+6Ci,2)t+(−3Pi−1+9Ci,1−9Ci,2+3Pi)t2

2. 光滑連續性條件

要求曲線段$ P_{i-1,i}(t) 和和和 P_{i,i+1}(t) 光滑連續，就是要求兩曲線段在點光滑連續，就是要求兩曲線段在點光滑連續，就是要求兩曲線段在點 P_{i} $處可導，則一充分非必要條件為：

d P i − 1 , i ( t ) d t ∣ t = 1 = d P i , i + 1 ( t ) d t ∣ t = 0 \frac{dP_{i-1,i}(t)}{dt}|_{t=1}=\frac{dP_{i,i+1}(t)}{dt}|_{t=0} dtdPi−1,i(t)∣t=1=dtdPi,i+1(t)∣t=0

化簡為：

P i − C i , 2 = C i + 1 , 1 − P i \mathbf{P_{i}}-\mathbf{C_{i,2}}=\mathbf{C_{i+1,1}}-\mathbf{P_{i}} Pi−Ci,2=Ci+1,1−Pi

即是矢量$ \mathbf{C_{i,2}}\mathbf{P_{i}} $ 與矢量 $ \mathbf{P_{i}}\mathbf{C_{i+1,1}} $相等。

3. 通過點的分段三階貝塞爾曲線繪制方案

3.1 基本繪制方法

由于得到的光滑連續性條件為矢量$ \mathbf{C_{i,2}}\mathbf{P_{i}} $ 與矢量 $ \mathbf{P_{i}}\mathbf{C_{i+1,1}} $相等，但其大小和方向未知，可以采取以下方案繪制光滑連續的分段三階貝塞爾曲線：

計算相鄰三點$ P_{i-1} 、、、 P_{i} 和和和 P_{i+1} 的中點的中點的中點 P_{i-1,i}^{avg} 和和和 P_{i,i+1}^{avg} $
計算中點的中點$ P_{i-1,i,i+1}^{avg} $
将點$ P_{i-1,i}^{avg} 和和和 P_{i,i+1}^{avg} 沿矢量沿矢量沿矢量 \mathbf{P_{i-1,i,i+1}^{avg} P_{i}}$進行平移
得到的兩點就作為前一段曲線的第二控制點$ C_{i,2} 和後一段曲線的第一控制點和後一段曲線的第一控制點和後一段曲線的第一控制點 C_{i+1,1} $

通過上述方式可以得到：

C i , 2 = P i − 1 4 ( P i + 1 − P i − 1 ) C_{i,2}=P_{i}-\frac{1}{4}(P_{i+1}-P_{i-1}) Ci,2=Pi−41(Pi+1−Pi−1)

C i + 1 , 1 = P i + 1 4 ( P i + 1 − P i − 1 ) C_{i+1,1}=P_{i}+\frac{1}{4}(P_{i+1}-P_{i-1}) Ci+1,1=Pi+41(Pi+1−Pi−1)

3.2 控制點伸縮方法

對$ P_{i} 點周圍的控制點點周圍的控制點點周圍的控制點 C_{i,2} 和和和 C_{i+1,1} 進行内縮或外擴，可以通過将點沿矢量進行内縮或外擴，可以通過将點沿矢量進行内縮或外擴，可以通過将點沿矢量\mathbf{ P_{i-1,i}^{avg} P_{i,i+1}^{avg}} 左右平移得到。要使控制點線段左右平移得到。要使控制點線段左右平移得到。要使控制點線段 C_{i,2}C_{i+1,1} 為原來的為原來的為原來的k 倍，若正方向指向點倍，若正方向指向點倍，若正方向指向點 P_{i} $，有左側控制點平移矢量：

V l e f t = 1 − k 4 ( P i + 1 − P i − 1 ) V_{left}=\frac{1-k}{4}(P_{i+1}-P_{i-1}) Vleft=41−k(Pi+1−Pi−1)

右側控制點平移矢量：

V r i g h t = − 1 ⋅ 1 − k 4 ( P i + 1 − P i − 1 ) V_{right}=-1\cdot\frac{1-k}{4}(P_{i+1}-P_{i-1}) Vright=−1⋅41−k(Pi+1−Pi−1)

則：

C i , 2 = C i , 2 + V l e f t = P i − k 4 ( P i + 1 − P i − 1 ) C_{i,2}=C_{i,2}+V_{left}=P_{i}-\frac{k}{4}(P_{i+1}-P_{i-1}) Ci,2=Ci,2+Vleft=Pi−4k(Pi+1−Pi−1)

C i + 1 , 1 = C i + 1 , 1 + V r i g h t = P i + k 4 ( P i + 1 − P i − 1 ) C_{i+1,1}=C_{i+1,1}+V_{right}=P_{i}+\frac{k}{4}(P_{i+1}-P_{i-1}) Ci+1,1=Ci+1,1+Vright=Pi+4k(Pi+1−Pi−1)

其中

當 k ∈ ( 0 , 1 ] k\in(0,1] k∈(0,1]時為内縮，

當 k ∈ [ 1 , + ∞ ) k\in[1,+\infty) k∈[1,+∞)時為外擴。

一般 k k k取 [ 0.4 , 0.6 ] [0.4,0.6] [0.4,0.6]較為合适。

3.3 非封閉的分段三階貝塞爾曲線

對于封閉的分段三階貝塞爾曲線通過點集$ \mathbf{P_i}={i|i=0,1,2,…,n} ，由于其首尾相連，可以看作在點，由于其首尾相連，可以看作在點，由于其首尾相連，可以看作在點P_0 前有點前有點前有點P_n ，在點，在點，在點P_n 後有點後有點後有點P_0$，進而确定第 1 1 1個點的第一控制點 C 0 , 1 C_{0,1} C0,1和第 n n n個點的第二控制點 C n , 2 C_{n,2} Cn,2。

而對于非封閉的分段三階貝塞爾曲線通過點集$ \mathbf{P_i}={i|i=0,1,2,…,n}$，無法确定第 1 1 1個點的第一控制點 C 0 , 1 C_{0,1} C0,1和第 n n n個點的第二控制點 C n , 2 C_{n,2} Cn,2，暫時的解決方案是：

令 C 0 , 1 = C 0 , 2 , C n , 1 = C n , 2 C_{0,1}=C_{0,2},C_{n,1}=C_{n,2} C0,1=C0,2,Cn,1=Cn,2。

3.4 三階貝塞爾曲線段的邊界求解

求解三階貝塞爾曲線段的邊界，就是求解方程的一階導數為 0 0 0的點：

( − 3 P i − 1 + 3 C i , 1 ) + ( 6 P i − 1 − 12 C i , 1 + 6 C i , 2 ) t + ( − 3 P i − 1 + 9 C i , 1 − 9 C i , 2 + 3 P i ) t 2 = 0 (- 3P_{i-1} + 3C_{i,1} ) +( 6P_{i-1} -12C_{i,1} + 6C_{i,2} )t +(- 3P_{i-1} + 9C_{i,1} - 9C_{i,2} + 3P_{i})t^2 =0 (−3Pi−1+3Ci,1)+(6Pi−1−12Ci,1+6Ci,2)t+(−3Pi−1+9Ci,1−9Ci,2+3Pi)t2=0

值得注意：上下邊界用點$ P_{i} 的的的 y 值求解，左右邊界用點值求解，左右邊界用點值求解，左右邊界用點 P_{i} 的的的 x $值求解；方程存在二次項為 0 0 0的情況，是以需要判斷二次項是否為 0 0 0。

下面呈現具體計算代碼：

//首先定義一個邊界類，用以友善計算曲線邊界
public class Boundary
{
    public Boundary(double xMax = Double.NaN, double xMin = Double.NaN, double yMax = Double.NaN, double yMin = Double.NaN)
    {
        XMax = xMax;
        XMin = xMin;
        YMax = yMax;
        YMin = yMin;
    }

    public Boundary(Point point1, Point point2)
    {
        if (point1.X >= point2.X)
        {
            XMax = point1.X;
            XMin = point2.X;
        }
        else
        {
            XMax = point2.X;
            XMin = point1.X;
        }
        if (point1.Y >= point2.Y)
        {
            YMax = point1.Y;
            YMin = point2.Y;
        }
        else
        {
            YMax = point2.Y;
            YMin = point1.Y;
        }
    }

    public double XMax { get; set; }
    public double XMin { get; set; }
    public double YMax { get; set; }
    public double YMin { get; set; }

    public void Resize(double xMax, double xMin, double yMax, double yMin)
    {
        XMax = xMax;
        XMin = xMin;
        YMax = yMax;
        YMin = yMin;
    }

    public void TryInclude(params Point[] points)
    {
        var length = points.Length;
        if (length <= 0 || points == null)
        {
            return;
        }

        double xmax;
        double xmin;
        double ymax;
        double ymin;
        bool isValueInvalid = Boundary.IsValueInvalid(this);
        if (isValueInvalid)
        {
            xmax = xmin = points[0].X;
            ymax = ymin = points[0].Y;
        }
        else
        {
            xmax = this.XMax;
            xmin = this.XMin;
            ymax = this.YMax;
            ymin = this.YMin;
        }

        double x;
        double y;
        for (int i = 0; i < length; i++)
        {
            x = points[i].X;
            y = points[i].Y;
            if (x > xmax)
            {
                xmax = x;
            }
            if (x < xmin)
            {
                xmin = x;
            }
            if (y > ymax)
            {
                ymax = y;
            }
            if (y < ymin)
            {
                ymin = y;
            }
        }

        this.XMax = xmax;
        this.XMin = xmin;
        this.YMax = ymax;
        this.YMin = ymin;
    }

    public void Merge(Boundary boundary)
    {
        bool isValueInvalid = Boundary.IsValueInvalid(this);
        if (!isValueInvalid)
        {
            if (Boundary.IsValueInvalid(boundary))
            {
                return;
            }
            if (boundary.XMax > this.XMax)
            {
                this.XMax = boundary.XMax;
            }
            if (boundary.XMin < this.XMin)
            {
                this.XMin = boundary.XMin;
            }
            if (boundary.YMax > this.YMax)
            {
                this.YMax = boundary.YMax;
            }
            if (boundary.YMin < this.YMin)
            {
                this.YMin = boundary.YMin;
            }
        }
        else
        {
            this.XMax = boundary.XMax;
            this.XMin = boundary.XMin;
            this.YMax = boundary.YMax;
            this.YMin = boundary.YMin;
        }
    }
    
    public static bool IsValueInvalid(Boundary boundary)
    {
        if (Double.IsNaN(boundary.XMax))
        {
            return true;
        }
        if (Double.IsNaN(boundary.XMin))
        {
            return true;
        }
        if (Double.IsNaN(boundary.YMax))
        {
            return true;
        }
        if (Double.IsNaN(boundary.YMin))
        {
            return true;
        }
        return false;
    }
}
//接下來給出具體貝塞爾曲線邊界的計算方法，該方法用以傳回某段曲線的邊界
private Boundary CalculateBezierSegmentBoundary(Point p0, Point c1, Point c2, Point p1)
{
    Boundary boundaryBS = new Boundary(p0, p1);
    double AY = 3 * c1.Y - 3 * c2.Y + p1.Y - p0.Y;
    double BY = 3 * (p0.Y + c2.Y - 2 * c1.Y);
    double CY = 3 * (c1.Y - p0.Y);
    double DY = p0.Y;

    double AX = 3 * c1.X - 3 * c2.X + p1.X - p0.X;
    double BX = 3 * (p0.X + c2.X - 2 * c1.X);
    double CX = 3 * (c1.X - p0.X);
    double DX = p0.X;

    double ay = 3 * AY;
    double by = 2 * BY;
    double cy = CY;
    double deltay = Math.Pow(by, 2) - 4 * ay * cy;

    double ax = 3 * AX;
    double bx = 2 * BX;
    double cx = CX;
    double deltax = Math.Pow(bx, 2) - 4 * ax * cx;

    bool IsAYApproximateZero = (ay > -1 * 1e-10) && (ay < 1e-10);
    bool IsAXApproximateZero = (ax > -1 * 1e-10) && (ax < 1e-10);
    if (IsAYApproximateZero || IsAXApproximateZero)
    {
        if (IsAYApproximateZero)
        {
            double ty = -1 * cy / by;
            if (ty > 0 && ty < 1)
            {
                double x = AX * ty * ty * ty + BX * ty * ty + CX * ty + DX;
                double y = AY * ty * ty * ty + BY * ty * ty + CY * ty + DY;
                boundaryBS.TryInclude(new Point(x, y));
            }
        }
        if (IsAXApproximateZero)
        {
            double tx = -1 * cx / bx;
            if (tx > 0 && tx < 1)
            {
                double x = AX * tx * tx * tx + BX * tx * tx + CX * tx + DX;
                double y = AY * tx * tx * tx + BY * tx * tx + CY * tx + DY;
                boundaryBS.TryInclude(new Point(x, y));
            }
        }
        return boundaryBS;
    }
    bool IsDeltaYGreaterThanZero = deltay > 0;
    bool IsDeltaXGreaterThanZero = deltax > 0;
    bool IsDeltaYApproximateZero = (deltay > -1 * 1e-10) && (deltay < 1e-10);
    bool IsDeltaXApproximateZero = (deltax > -1 * 1e-10) && (deltax < 1e-10);
    if (IsDeltaYGreaterThanZero || IsDeltaXGreaterThanZero)
    {
        if (!IsAYApproximateZero)
        {
            if (IsDeltaYGreaterThanZero)
            {
                double t1y = (-1 * by + Math.Pow(deltay, 0.5)) / (2 * ay);
                double t2y = (-1 * by - Math.Pow(deltay, 0.5)) / (2 * ay);
                if (t1y > 0 && t1y < 1)
                {
                    double x = AX * t1y * t1y * t1y + BX * t1y * t1y + CX * t1y + DX;
                    double y = AY * t1y * t1y * t1y + BY * t1y * t1y + CY * t1y + DY;
                    boundaryBS.TryInclude(new Point(x, y));
                }
                if (t2y > 0 && t2y < 1)
                {
                    double x = AX * t2y * t2y * t2y + BX * t2y * t2y + CX * t2y + DX;
                    double y = AY * t2y * t2y * t2y + BY * t2y * t2y + CY * t2y + DY;
                    boundaryBS.TryInclude(new Point(x, y));
                }
            }
        }
        if (!IsAXApproximateZero)
        {
            if (IsDeltaXGreaterThanZero)
            {
                double t1x = (-1 * bx + Math.Pow(deltax, 0.5)) / (2 * ax);
                double t2x = (-1 * bx - Math.Pow(deltax, 0.5)) / (2 * ax);
                if (t1x > 0 && t1x < 1)
                {
                    double x = AX * t1x * t1x * t1x + BX * t1x * t1x + CX * t1x + DX;
                    double y = AY * t1x * t1x * t1x + BY * t1x * t1x + CY * t1x + DY;
                    boundaryBS.TryInclude(new Point(x, y));
                }
                if (t2x > 0 && t2x < 1)
                {
                    double x = AX * t2x * t2x * t2x + BX * t2x * t2x + CX * t2x + DX;
                    double y = AY * t2x * t2x * t2x + BY * t2x * t2x + CY * t2x + DY;
                    boundaryBS.TryInclude(new Point(x, y));
                }
            }
        }
        return boundaryBS;
    }
    else if (IsDeltaYApproximateZero || IsDeltaXApproximateZero)
    {
        if (!IsAYApproximateZero)
        {
            if (IsDeltaYApproximateZero)
            {
                double ty = -1 * by / (2 * ay);
                if (ty > 0 && ty < 1)
                {
                    double x = AX * ty * ty * ty + BX * ty * ty + CX * ty + DX;
                    double y = AY * ty * ty * ty + BY * ty * ty + CY * ty + DY;
                    boundaryBS.TryInclude(new Point(x, y));
                }
            }
        }
        if (!IsAXApproximateZero)
        {
            if (IsDeltaXApproximateZero)
            {
                double tx = -1 * bx / (2 * ax);
                if (tx > 0 && tx < 1)
                {
                    double x = AX * tx * tx * tx + BX * tx * tx + CX * tx + DX;
                    double y = AY * tx * tx * tx + BY * tx * tx + CY * tx + DY;
                    boundaryBS.TryInclude(new Point(x, y));
                }
            }
        }
        return boundaryBS;
    }
    else
    {
        return boundaryBS;
    }
}
//最後可以使用Boundary類中的Merge方法合并邊界以計算所有分段三階貝塞爾曲線邊界

分段三階貝塞爾曲線的光滑連續性條件及其WPF應用分段三階貝塞爾曲線的光滑連續性條件

分段三階貝塞爾曲線的光滑連續性條件

1. 貝塞爾曲線簡介 [1]

2. 光滑連續性條件

3. 通過點的分段三階貝塞爾曲線繪制方案

3.1 基本繪制方法

3.2 控制點伸縮方法

3.3 非封閉的分段三階貝塞爾曲線

3.4 三階貝塞爾曲線段的邊界求解

繼續閱讀

存在是“一”（分化規律）

WPF中檢測U盤插入、拔出執行個體

XAML中綁定一個枚舉

wpf連接配接資料庫

WPF學習感受WPF的XAML好的UI類庫

對CPersistentFrame的問題

閑話WinFrom與WPF(4) 架構篇——結合VS擴充,插件式開發架構

WPF 中 InitializeComponent 不存在解決方案

簡述Shadow Mapping和Shadow Volume的新方法簡介相關工作高分辨率Shadow Mapping壓縮基于聚類和度量樹的實時陰影渲染方法度量樹建構結果與讨論參考文獻

opengl 和direct3d坐标系的差別

第18章自定義元素(1)——概述

WPF實作3D粒子波浪效果

bgfx跨平台渲染庫學習筆記:概述

WPF開發帶水印的輸入框