快速分段3次样条曲线拟合和折线重采样算法实现

快速分段3次样条曲线拟合和折线重采样算法实现

cheungmine 

（保留所有权利。本文可以在互联网上转载，但不允许出版和印刷）

 

本文采用3次样条函数，用分段插值的快速计算方法，实现了用鼠标在屏幕上绘制任意光滑的曲线，并同时使用折线重采样的拟合方法，去除多余的插值点。本文所叙述的算法，可以用来绘制等高线等光滑曲线，并且由于采用了折线的重采样，以最小的数据量保证了绘图的精确度。

 

本文以作者2002年根据《计算方法》（同济大学出版社）一书所述的算法为基础修改而得，因为时间久远，这里不再叙述什么是张力样条的具体公式。仅仅对张力样条做简单介绍。

 

1 问题的提出​

 

在曲线数字化，或绘制曲线的应用中，通常要求曲线是光滑的，并且通过控制点。这在高等数学里，有专门的术语描述--导数的连续性。我们这里做个形象的比喻：拿一根弹性强的硬钢丝，通过固定在木板上的一系列点，就得到一条满足我们要求的光滑的曲线。当我们在钢丝的2端，施加巨大的拉力，可以看到，钢丝在固定点处的光滑程度越来越小，直到变为折线段。这个力就是张力。我们用一个系数表示，就是张力系数f。张力系数是张力的倒数，当f=0时，张力最大，曲线退化为折线段。当f=1，基本类似抛物线。比较合适的取值是：

f=0.1666667

 

 

2 样条函数​

 

       关于样条函数，我实在忘记原理了。请读者自行学习原理的部分。实际当中，通常是用分段拟合的方法：

       给定木板上（p0，p1，p2，p3）4个固定点，要求样条钢丝依次通过这4个点，我们就可以确定p1，p2之间的样条函数公式。根据这个公式和一定的步长step，就可以插值出p1，p2之间的一系列点q1,q2,...,qn。依次用线段连接这些点，就得到用折线段按拟合出来的样条曲线。步长越小，拟合出的折点越多，越接近曲线的真实值。然而数据量就越多。

 

 

 

 

3 分段拟合方法​

 

       给点一系列固定点p0，p1，p2，...，pn-3，pn-2，pn-1，依次采样上面的分段拟合方法：

1）  拟合p0,p1,p2,p3得到p1,p2间的点：q11,q12,q13,...,q1m。

2）  拟合q1m,p2,p3,p4得到p2,p3间的点：q21,q22,q23,...,q2n。

3）  拟合q2n,p3,p4,p5得到p3,p4间的点：q31,q32,q33,...,q3w。

4）  ...

按上面的方法，一直拟合出：pn-3，pn-2间的点。

 

上面这些拟合点和原来的控制点，按次序就构成了整条函数曲线的拟合折线。采样这种拟合方法，计算量小，极易编制程序。

 

4 拟合的精度​

但是，这种拟合得出的精度受步长的限制，我们一方面想提高精度Precision，只有把step设置的很小。又不想增加很大的数据量。这是矛盾的：

Precision <-> 1/step

 

       这种等间距step的拟合方式，无论在曲线的平缓处，还是转弯处，都采用了一样的拟合步长。而通常我们只希望在曲线的平缓处只要少量的几个点就可以达到我们的要求。而在转弯处，用密集些的点拟合。

 

       因此引入了重采样的方法，去除一些不必要的点。重采样的依据就是：

       设折线段上一系列点p0,p1,p2,...,pm，pn。分别计算p1,p2,...,pm到直线p0,pn的距离。如果这些距离的最大值pk，小于某个限值dist，则曲线上的点p1,p2,...,pm都不是必要的，我们只要保留p0,pn即可。（式1）

如果pk大于我们规定的限值dist，则分别以（p0,p1,p2,...,pk-2,pk-1,pk）和（pk,pk+1,pk+2,...,pm,pn）作为2条新的折线重新计算上面的（式1）。

 

在编制程序的时候，使用递归可以很容易解决这个问题。在我提供的测试程序中，可以看到2种情况（重采样和不重采样）下的点数是相差很大的。

 

 

 

5 结论​

 

       到此，我们近乎完美地解决问题：采样分段样条插值，并且采用重采样，减小数据量，同时保证精度不受损失。

 

6 下面就是完整的代码。我采样C语言写的：

/*==============================================================================
     cheungmine - All rights reserved
     April 5, 2009
==============================================================================*/
#include <stdio.h>
#include <float.h>
#include <assert.h>
#include <math.h>
 
#ifndef CG_BIPI
#define CG_BIPI        6.28318530717958647692
#endif
 
#ifndef IN
#define IN
#endif
 
#ifndef OUT
#define OUT
#endif
 
#ifndef INOUT
#define INOUT
#endif
 
/* Vertex structure */
typedef struct                       
{
     double  x, y;
} cg_vertex_t;
 
/**
 * CG_vertex_get_dist_to_line
 *  计算点到直线的距离和垂足ft
 */
double CG_vertex_get_dist_to_line (
              IN const cg_vertex_t  *pt, 
              IN const cg_vertex_t  *start, 
              IN const cg_vertex_t  *end, 
              OUT cg_vertex_t       *ft)
{
     double a1 = pt->x - start->x;
     double a2 = pt->y - start->y;
     double s = sqrt(a1*a1 + a2*a2);
 
     a1 = atan2(a2, a1);
     if(a1 < 0) a1 += CG_BIPI;
 
     a2 = atan2((end->y-start->y), (end->x-start->x));
     if(a2 < 0) a2 += CG_BIPI;
 
     a1 = fabs(sin(a1-a2)) * s;
 
     if (ft){
         s *= cos(a1-a2);
         ft->x = start->x + s*cos(a2);
         ft->y = start->y + s*sin(a2);
     }
     return a1;
}
 
 
 
/*===========================================================================
                                 Spline3 Functions
===========================================================================*/
/* 返回pl的点到直线[start,end]的最大距离点的索引imax和最大距离smax */
static int max_dist_at(
                   const   cg_vertex_t* start,
                   const   cg_vertex_t* end,
                   const   cg_vertex_t* pl, int np,
                   double *smax)
{
     int     i, imax = 0;
     double  s;
     *smax = 0;
 
     // 计算全部点到直线的距离, 找出最大的点
     for(i=0; i<np; i++){
         s = CG_vertex_get_dist_to_line(&pl[i], start, end, 0);
         if (s > *smax){
              *smax = s;
              imax = i;
         }
     }
     return imax;
}
 
/* 递归重采样方法*/
static void resample_recursive(
                   const  cg_vertex_t* start,
                   const  cg_vertex_t* end,
                   const  cg_vertex_t* pl,
                   int    np,
                   cg_vertex_t *pout,
                   int   *nout,
                   const double *dist)
{
     double smax;
     int    imax;
     if (np==0) return;
     imax = max_dist_at(start, end, pl, np, &smax);
     if (imax==0||imax==np-1) return;
     if (smax<*dist) return;
 
     resample_recursive( start, &pl[imax], pl, imax, pout, nout, dist);
 
     pout[(*nout)++] = pl[imax];
 
     resample_recursive( &pl[imax], end, &pl[imax+1], np-imax-1, pout, nout, dist);
}
 
/**
 * CG_vertices_fit_resample
 * 曲线重采样拟合, 返回拟合之后的pout数组的点数
 */
int CGAL_CALL CG_vertices_fit_resample (
                       IN const  cg_vertex_t* start, /* 折线或多边形的起点*/
                       IN const  cg_vertex_t* end,   /* 折线或多边形的终点*/
                       IN const  cg_vertex_t* mids,  /* 折线或多边形的中间点数组*/
                       IN      int            nmid,  /* 折线或多边形的中间点数组的点数*/
                       OUT   cg_vertex_t*     pout,  /* 返回拟合后的中间点数组, 点数组由客户分配, 
                                                             至少与输入的中间点数组的点数相同*/
                       IN   double            dist   /* 拟合的最小距离, 必须指定>0的有效值*/                    
)
{
     int  nout = 0;
 
     resample_recursive(start, end, mids, nmid, pout, &nout, &dist);
 
     return nout;
}
 
/**
 * CG_vertex_insert_pt_spline
 * 计算张力样条插值点: p1, p2之间的点, 返回插值点数目
 * p0-----p1........p2------p3
 * 本程序为减少数据复制, 采用points数组返回插值点或拟合后的值, 
 * 因此调用者必须分配足够的空间以存储中间计算结果和返回的值
 * 必须满足:
 *            size > (|p1,p2|/step)*2
 * 如果设置的step太小, 系统则自动根据size/2来调整step, 确保
 * 计算能正确进行
 */
int CGAL_CALL CG_vertex_insert_pt_spline (
                       IN   double         step,    /* 插值的步长, 必须指定有效的步长*/
                       IN   cg_vertex_t    p0,      /* 起点*/
                       IN   cg_vertex_t    p1,      /* 插值段的起点*/
                       IN   cg_vertex_t    p2,      /* 插值段的终点*/
                       IN   cg_vertex_t    p3,      /* 终点*/
                       IN OUT cg_vertex_t *points,  /* 返回的插值点数组, 数组由调用者分配*/
                       IN   int            size,    /* 插值点数组的最大尺寸, 此尺寸必须足够容纳2倍的插值点数*/
                       IN   double         ratio,  /* 张力系数[0,1], 最佳0.1666667. -1 为取默认.1666667 */
                       IN   double         dist     /* 拟合的最小距离, 必须>0. <=0 不拟合*/
                   )
{
     int     i, count, at;
     double  x, y, ca, sa, s12, h1, h2, h3, u1, u2, v1, v3, d1, d2, D, M1, M2;
     cg_vertex_t  end = {p2.x, p2.y};
 
     // 必须计算p1,p2距离, 确保不会溢出
     d1 = p2.x-p1.x;
     d2 = p2.y-p1.y;
     s12 = sqrt(d1*d1 + d2*d2);
 
     // p1,p2距离太小, 返回
     if (dist>0 && s12 < dist)
         return 0;
 
     // 步距离太长, 计算p1与p2中点即可
     if (step > s12/2 ){         
         points[0].x = (p2.x+p1.x)/2;
         points[0].y = (p2.y+p1.y)/2;
         return 1;
     }
 
     // 平移原点至p1
     p0.x -= p1.x;
     p0.y -= p1.y; 
     p2.x = d1;
     p2.y = d2;
     p3.x -= p1.x;
     p3.y -= p1.y;
 
     // 旋转至p1->p2为x轴方向
     ca = p2.x / s12;   // cos(a) 
     sa = p2.y / s12;   // sin(a)
 
     x = p0.x * ca + p0.y * sa;
     p0.y = -p0.x * sa + p0.y * ca;
     p0.x = x;
     p2.x = s12;
     p2.y = 0;
     x = p3.x * ca + p3.y * sa;
     p3.y = -p3.x * sa + p3.y * ca;
     p3.x = x;
 
     // 判断是否是单值函数
     if ( p0.x * p2.x >= 0 || p2.x * (p2.x-p3.x) >= 0 )
          return 0;
 
     // 计算系数, 这里为清晰起见, 没有优化代码
     h1 = -p0.x;
     h2 = p2.x;
     h3 = p3.x - p2.x;
     u1 = h2/(h1+h2);   // u1=p2.x/(p2.x-p0.x);
     u2 = h2/(h2+h3);   // u2=p2.x/(p3.x-p1.x);
 
     v1 = -p0.y;
     v3 = p3.y;
 
     d1 = -6.0*p0.y/((h1+h2)*p0.x);
     d2 = 6.0*p3.y/(h3*p3.x);
     D  = 4.0-u1*u2;
 
     if(D < 0.001)
         return 0;
 
     M1 = (2*d1-u1*d2)/D;
     M2 = (2*d2-u2*d1)/D;
 
     // 计算插值点数, 如果点数超过size/2, 则自动增大step
     count = (int) floor((s12+step/2)/step);
     if (count > size/2){
         count = size/2;
         step=s12/(count+1);
     }
 
     // 张力系数, 最佳=0.16666667
     if(ratio<0) ratio = 0.16666667;
 
     x = 0.0;
     y = 0.0;
     at = (dist>0 && count<=size/2)? count:0;  
     for(i=0; i < count; i++){
         x += step;
         y =  ratio*x*(M1*(h2-x)*(x-2*h2)+M2*(x-h2)*(x+h2))/h2;
         points[at+i].x = x*ca - y*sa + p1.x;
         points[at+i].y = x*sa + y*ca + p1.y;
     }
 
     // 需要重采样: p1 as start, p2 as end
     if (dist>0)
         count = CG_vertices_fit_resample(&p1, &end, points+at, count, points, dist);
 
     return count;
}

7 我还提供了一个MFC的测试程序。你只要用鼠标在上面点即可。测试程序演示如何调用上面的函数，实现绘制连续样条曲线。测试程序在下面的链接可以下载得到

 

程序没实现封闭的样条曲线。在封闭部分有：

pn-3~~~~qw~pn-2=====pn-1=====p0=====p1~q1~~~~p2

其中为做曲线计算的部分以=====表示。读者可以多次调用CG_vertex_insert_pt_spline以分别计算拟合pn-2=====pn-1、pn-1=====p0和p0=====p1段的曲线。

 

 

2009-4-5 by win32 -cheungmine -sdk