贝塞尔曲线拟合

接上一篇Blog,这里用贝塞尔曲线来平滑多个点。

和样条插值不同,在多个点上作贝塞尔曲线的时候,曲线只穿过首尾两个点,中间的点都是作为控制点。

移动控制点,曲线也随之形变,可以造成一种拉扯的效果。在各种作图工具中,经常使用贝塞尔曲线来画曲线。一般的操作都是先画一条线段,然后可以通过拖动一个控制点来调整线段的弯曲程度。

作多点贝塞尔曲线只需要一个公式。所有的点的X值,被归一化到[0,1]区间内。
具体理论,可以参考这个页面Bézier curves。89年创建的,可有年头了。

这里还是贴代码吧。

首先需要得到X区间的总长度。

CGPoint startPt = [[_points objectAtIndex:0] CGPointValue];
CGPoint endPt = [[_points objectAtIndex:(self.pointCount - 1)] CGPointValue];
float amount = endPt.x - startPt.x;

然后就是曲线方程了,这个比样条插值要简单不少。
rank是指总的阶数,也就是实际的点数。这个函数表示n个点的贝塞尔曲线在x处的值。
这里的ux属于区间[0,1]

float (^bezierSpline)(int rank, float ux) = ^(int rank, float ux) {
        
    float p = 0.0f;
        
    for (int i = 0; i < rank; i++)
    {
        CGPoint pt = [[_points objectAtIndex:i] CGPointValue];
            
        p += pt.y * powf((1 - ux), (rank - i - 1)) * powf(ux, i) * (factorial(rank - 1) 
                  / (factorial(i) * factorial(rank - i - 1)));
    }
    
    return p;
};

下面很容易了,画图的时候步长为1,求得ux之后,带入方程得到点的y值,画曲线。

    [path moveToPoint:startPt];
    
    for (float curX = startPt.x; (curX - endPt.x) < 1e-5; curX += 1.0f)
    {
        float u = (curX - startPt.x) / amount;
        [path addLineToPoint:
            CGPointMake(curX, bezierSpline(self.pointCount, u))];
    }
    
    CGContextSetLineWidth(context, 1.0f);
    CGContextSetStrokeColorWithColor(context, [UIColor blackColor].CGColor);
    CGContextAddPath(context, path.CGPath);
    CGContextStrokePath(context);
}

和三次样条相比,由于不经过中间点,丢失细节比较多,平滑度更好。

平滑前
no bezier

三次样条
no bezier

贝塞尔曲线
no bezier

示例代码: Sample-CurveFit
和上个例子放在一个repo里。

参考: Wiki

三次样条插值和曲线拟合

很多东西不在手上用着就容易忘,尤其是书本知识。就弄这么个类别,叫作“书到用时方恨少”,来记录这些知识。

曲线拟合是一个“数值计算“中的一个基本内容。在实际的项目中,使用拟合的目的就是从有限个点得到一条平滑曲线。曲线本身也是由点构成的,所以如何从有限个点得到曲线上的其它点,就是插值所关注的内容。

插值的方法有很多,把这些个点逐个用直线段连起来也是一种插值。

样条插值是一种工业设计中常用的、得到平滑曲线的一种插值方法,三次样条又是其中用的较为广泛的一种。

Google 三次样条插值可以看得到不少材料,这里就不罗列公式了,直接看看在代码里,我们怎么做。

首先我们需要各个点的坐标,以x,y表示。

        const int len = [_points count];
        float x[len];
        float y[len];
        for (int i = 0; i < len; i++)
        {
            CGPoint p = [[_points objectAtIndex:i] CGPointValue];
            x[i] = p.x;
            y[i] = p.y;
        }

取变量x,y

从算法中可以得知,我们的目标是样条插值函数,这是一个分段函数,x最高次数为三次,在各个点二次连续可导以保证最终函数曲线的光滑性。
我们每两个点求一个三次函数,我们有n个点,那么这里就需要4(n-1)个方程。
目前我们有n个点的坐标,有n-2个连接点,有n个函数两次连续可导,这里有n+n-2+2*(n-2)共4n-6个方程,还差两个条件。
这里一般有三种处理方法,最方便的,也是我们这里使用的是自然三次样条,也就是在首尾两个点上二次导为0。

具体计算不在此列举了,根据算法构建一个方程组求一组中间值sx,左边是一个三对角矩阵。

        float h[len];
        float u[len];
        float lam[len];
        for (int i = 0; i < len-1; i++)
        {
            h[i] = x[i+1] - x[i];
        }    

        u[0] = 0;
        lam[0] = 1;
        for (int i = 1; i < (len - 1); i++)
        {
            u[i] = h[i-1]/(h[i] + h[i-1]);
            lam[i] = h[i]/(h[i] + h[i-1]);
        }
        
        float a[len];
        float b[len];
        float c[len];
            
        float m[len][len];
        for (int i = 0; i < len; i++)
        {
            for (int j = 0; j < len; j++)
            {
                m[i][j] = 0;
            }
            if (i == 0)
            {
                m[i][0] = 2;
                m[i][1] = 1;
                
                b[0] = 2;
                c[0] = 1;
            }
            else if (i == (len - 1))
            {
                m[i][len - 2] = 1;
                m[i][len - 1] = 2;
                
                a[len-1] = 1;
                b[len-1] = 2;
            }
            else
            {
                m[i][i-1] = lam[i];
                m[i][i] = 2;
                m[i][i+1] = u[i];
                
                a[i] = lam[i];
                b[i] = 2;
                c[i] = u[i];
            }
        }

求三对角矩阵,自下而上对角线上的参数是a,b,c

当然需要得到方程组右边的值

        float g[len];
        g[0] = 3 * (y[1] - y[0])/h[0];
        g[len-1] = 3 * (y[len - 1] - y[len - 2])/h[len - 2];
        for (int i = 1; i < len - 1; i++)
        {
            g[i] = 3 * ((lam[i] * (y[i] - y[i-1])/h[i-1]) + u[i] * (y[i+1] - y[i])/h[i]);
        }

下面就是求解这个方程组了,对于三对角矩阵,使用追赶法

        //< Solve the Equations
        float p[len];
        float q[len];

        p[0] = b[0];
        for (int i = 0; i < len - 1; i++)
        {
            q[i] = c[i]/p[i];
        }

        for (int i = 1; i < len; i++)
        {
            p[i] = b[i] - a[i]*q[i-1];
        }
        
        float su[len];
        float sq[len];
        float sx[len];
        
        su[0] = c[0]/b[0];
        sq[0] = g[0]/b[0];
        for (int i = 1; i < len - 1; i++)
        {
            su[i] = c[i]/(b[i] - su[i-1]*a[i]);
        }

        for (int i = 1; i < len; i++)
        {
            sq[i] = (g[i] - sq[i-1]*a[i])/(b[i] - su[i-1]*a[i]);
        }
        
        sx[len-1] = sq[len-1];
        for (int i = len - 2; i >= 0; i--)
        {
            sx[i] = sq[i] - su[i]*sx[i+1];
        }        

求得了参数,现在就得到分段插值函数了。这种函数实际上有两个参数。
一个是所在的段位,这个是整数。另一个才是自变量。
在Objective-C里用block做一个临时的函数很是方便。

        double (^func)(int k, float vX) = ^(int k, float vX) {
            double p1 =  (ph[k] + 2.0 * (vX - px[k])) * ((vX - px[k+1]) * (vX - px[k+1])) * py[k] / (ph[k] *ph[k] * ph[k]);
            double p2 =  (ph[k] - 2 * (vX - px[k+1])) * powf((vX - px[k]), 2.0f) * py[k+1] / powf(ph[k], 3.0f);
            double p3 =  (vX - px[k]) * powf((vX - px[k+1]), 2.0f) * psx[k] / powf(ph[k], 2.0f);
            double p4 =  (vX - px[k+1]) * powf((vX - px[k]), 2.0f) * psx[k+1] / powf(ph[k], 2.0f);
            return p1 + p2 + p3 + p4;
        };

到此为止,得到了曲线的函数。当然在实际画曲线的时候,我们仍然是逐个点画,
取点的时候取得密集一些就形成曲线了。这里delta取为1。

        for (int i = 0; i < [_points count]; i++)
        {
            CGPoint pt = [[_points objectAtIndex:i] CGPointValue];
            if (i == 0)
            {
                CGPathMoveToPoint(path, NULL, pt.x, pt.y);
            }
            else
            {
                CGPoint curP = [[_points objectAtIndex:i-1] CGPointValue];
                float delta = 1.0f;
                for (float pointX = curP.x; fabs(pointX - pt.x) > 1e-5f; pointX += delta)
                {
                    float pointY = func(i-1, pointX);
                    CGPathAddLineToPoint(path, NULL, pointX, pointY);
                }
            }
        }

拟合前

拟合后

示例代码在GitHub上: Sample-CurveFit