以前学数值计算经常要拟合曲线,昨天突然收到一医院朋友的求助, 要将颈椎椎体的样点拟合为圆弧曲线,如下图, 具体问题如下: 1、测量方法：从C2椎体齿状突后上缘A点至C7椎体后下缘B点，将每一颈椎椎体后上缘和后下缘坐标定点，记录相应的像素位置，然后输入专门设计的软件程序连接各标点并优化成弧线c？，连接A、B两点成直线a，测量弦a和弧c的长度。 2、计算方法：将数据输入公式，解二元一次方程c＝2πr/360，sin（α/2）＝a/2/r（r为模拟圆的半径，α为弧c所对的圆心角），自动求得人体颈椎正常的生理曲度α。 开始求解: 首先,就想到了最小二乘法,  根据题意, 该圆必然通过两个端点,圆心位置就在大弦的垂直平分线上, 为了节省搜索时间,采用下述方法: 1.首先求除端点外剩余点的中心点, 2.然后,通过两端点和中心点确定一圆 3.求大弦的垂直平分线,并求圆心到大弦中点的距离H 4.从圆心出发,沿着大弦的垂直平分线往两边以0.01H为单位往两边搜索,利用最小二乘原理确定最合理值。 程序附后,结果如下,验证合理: 数据点个数为 12相应坐标为：   1:     6.440   3.600   2:     6.380   2.610   3:     6.400   2.340   4:     6.480   2.020   5:     6.540   1.800   6:     6.670   1.450   7:     6.720   1.260   8:     6.920   0.990   9:     7.060   0.820  10:     7.390   0.490  11:     7.600   0.360  12:     7.930   0.000计算结果为：圆心 X: 9.8070  Y: 2.8852弧长c:  4.1415弦长a:  3.8962生理曲度α: 68.9396 程序代码: void Calculate(){ //m_dX,m_dY分别存储样本点坐标，TWE为点个数 double x0 = m_dX[0]; double y0 =m_dY[0]; double x1 = m_dX[TWE-1]; double y1 =m_dY[TWE-1];     //获取中间点集的中点    double x2 = 0,y2 =0; for( size_t i=1;i< TWE-1;i++) {   x2+=m_dX[i];  y2+=m_dY[i]; } x2/=(TWE-2); y2/=(TWE-2);     //分别求两个弦的中点和垂直平分线  double X1=(x0+x1)/2.; double Y1=(y0+y1)/2.; double K1 = -(x1-x0)/(y1-y0); double X2=(x1+x2)/2.; double Y2=(y1+y2)/2.; double K2 = -(x2-x1)/(y2-y1);    //求取圆心和半径 m_X = (Y2-Y1+K1*X1-K2*X2)/(K1-K2); m_Y = Y1+K1*(m_X-X1); r = sqrt((x0-m_X)*(x0-m_X) + (y0-m_Y)*(y0-m_Y));    //求圆心到大弦中点的距离 double H = sqrt((m_X-X1)*(m_X-X1) + (m_Y-Y1)*(m_Y-Y1));    //求误差平方和 eps = 0.; for(  i=1;i< TWE-1;i++) {   double len=r - sqrt((m_dX[i]-m_X)*(m_dX[i]-m_X) + (m_dY[i]-m_Y)*(m_dY[i]-m_Y));  eps+=len*len; }    eps/=(TWE-2);     //采用最小二乘法搜索最合理圆心位置 double realX = m_X; double realY = m_Y; //虚拟圆心坐标 virX,virY virr vireps; double virX,virY,virr,vireps;  for(int j=-100;j<=100;j++) {  virX = realX + 0.01*H*j*1./(sqrt(K1*K1+1));  virY = realY + 0.01*H*j*K1/(sqrt(K1*K1+1));  virr = sqrt((x0-virX)*(x0-virX) + (y0-virY)*(y0-virY));    vireps =0.;  for(  i=1;i< TWE-1;i++)  {    double len=virr - sqrt((m_dX[i]-virX)*(m_dX[i]-virX) + (m_dY[i]-virY)*(m_dY[i]-virY));   vireps+=len*len;  }  vireps/=(TWE-2);    if(vireps<eps)  {   eps=vireps;   m_X=virX;   m_Y=virY;   r=virr;  } } //求相关参数 a = sqrt( (m_dX[0]-m_dX[TWE-1])*(m_dX[0]-m_dX[TWE-1]) +  (m_dY[0]-m_dY[TWE-1])*(m_dY[0]-m_dY[TWE-1])); a1= 2.*asin(a/2./r)*180./3.1415926; c=r*2.*asin(a/2./r);}  

评论