分段三阶贝塞尔曲线的光滑连续性条件

1. 贝塞尔曲线简介 ^[1]

计算机图形学中有一类很常用的曲线，俗称贝塞尔曲线。1962年，法国数学家Pierre Bézier 第一个研究了这种矢量绘制曲线的方法，并给出了详细的计算公式，因此按照这样的公式绘制出来的曲线就用他的姓氏来命名是为贝塞尔曲线。很多程序语言都有实现贝塞尔曲线的API，而该曲线本身也拥有强大的近似其它曲线的能力，即使一条不能够胜任，那么分段的多条贝塞尔曲线也足够用来近似我们想绘制的曲线。

这里给出贝塞尔曲线的一般公式：
$P(t)=\sum_{i=0}^n{B_{i,n}(t)\mathbf{P_i}}$
其中：
$B_{i,n}(t)=\frac{n!}{i!(n-i)!}t^i(1-t)^{n-i},t\in[0,1]$
对于以 $P_{i-1}$ 和 $P_{i}$ 作为起点和终点的三阶贝塞尔曲线段，具有两个控制点 $C_{i,1}$ 、 $C_{i,2}$ ，方程为：
$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$
对应的方程一阶导数即为：
$\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$

2. 光滑连续性条件

要求曲线段 $P_{i-1,i}(t)$ 和 $P_{i,i+1}(t)$ 光滑连续，就是要求两曲线段在点 $P_{i}$ 处可导，则一充分非必要条件为：
$\frac{dP_{i-1,i}(t)}{dt}|_{t=1}=\frac{dP_{i,i+1}(t)}{dt}|_{t=0}$
化简为：
$\mathbf{P_{i}}-\mathbf{C_{i,2}}=\mathbf{C_{i+1,1}}-\mathbf{P_{i}}$
即是矢量 $\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}-\frac{1}{4}(P_{i+1}-P_{i-1})$
$C_{i+1,1}=P_{i}+\frac{1}{4}(P_{i+1}-P_{i-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_{left}=\frac{1-k}{4}(P_{i+1}-P_{i-1})$
右侧控制点平移矢量：
$V_{right}=-1\cdot\frac{1-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})$
$C_{i+1,1}=C_{i+1,1}+V_{right}=P_{i}+\frac{k}{4}(P_{i+1}-P_{i-1})$
其中
当 $k\in(0,1]$ 时为内缩，
当 $k\in[1,+\infty)$ 时为外扩。
一般 $k$ 取 $[0.4,0.6]$ 较为合适。

3.3 非封闭的分段三阶贝塞尔曲线

对于封闭的分段三阶贝塞尔曲线通过点集 $\mathbf{P_i}=\{i|i=0,1,2,...,n\}\$ ，由于其首尾相连，可以看作在点 $P_0$ 前有点 $P_n$ ，在点 $P_n$ 后有点 $P_0$ ，从而确定第 $1$ 个点的第一控制点 $C_{0,1}$ 和第 $n$ 个点的第二控制点 $C_{n,2}$ 。
而对于非封闭的分段三阶贝塞尔曲线通过点集 $\mathbf{P_i}=\{i|i=0,1,2,...,n\}\$ ，无法确定第 $1$ 个点的第一控制点 $C_{0,1}$ 和第 $n$ 个点的第二控制点 $C_{n,2}$ ，暂时的解决方案是：
令 $C_{0,1}=C_{0,2},C_{n,1}=C_{n,2}$ 。

3.4 三阶贝塞尔曲线段的边界求解

求解三阶贝塞尔曲线段的边界，就是求解方程的一阶导数为 $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$
值得注意：上下边界用点 $P_{i}$ 的 $y$ 值求解，左右边界用点 $P_{i}$ 的 $x$ 值求解；方程存在二次项为 $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应用问题的一些解决方法

分段三阶贝塞尔曲线的光滑连续性条件及在WPF应用问题的一些解决方法

分段三阶贝塞尔曲线的光滑连续性条件

1. 贝塞尔曲线简介 ^[1]

2. 光滑连续性条件