bezier曲线拟合椭圆弧线

椭圆弧线用bezier曲线拟合。
先计算出椭圆中心起始角度旋转角度

$Step 1: Compute(x'_1, y'_1)$
$\begin{pmatrix} x'_1 \\ \\ y'_1 \end{pmatrix} = \begin{pmatrix} \cos\varphi & \sin\varphi \\ \\ -\sin\varphi & \cos\varphi \end{pmatrix} \cdot \begin{pmatrix} \dfrac{ x_1-x_2}{2} \\ \\ \dfrac{ y_1-y_2}{2} \end{pmatrix}$

$Step 2: Compute(c'_x, c'_y)$
$\begin{pmatrix} c'_x \\ \\ c'_y \end{pmatrix} = \pm \sqrt { \dfrac{ r^2_x y^2_y - r^2_x(y'_1)^2 - r^2_y (x'_1)^2} {r^2_x (y'_1)^2 + r^2_y(x'_1)^2 } } \begin{pmatrix} \dfrac{ r_xy'_1}{r_y} \\ \\ -\dfrac{ r_y x'_1}{r_x} \end{pmatrix}$

在这里插入图片描述

然后拟合
在这里插入图片描述

在这里插入图片描述

附代码

#include <math.h>
#include <stdio.h>
#include <errno.h>
#include <fcntl.h>
#include <unistd.h>
#include <stdlib.h>

#define MAX(a, b) (((a) > (b)) ? (a) : (b))
#define MIN(a, b) (((a) < (b)) ? (a) : (b))
#define PI 3.1415926

static float angle(float ux, float uy, float vx, float vy) {   
  float res = (ux*vx + uy*vy)/sqrtf((ux*ux+uy*uy)*(vx*vx+vy*vy));
  res = MIN(MAX(res,-1),1);  //-1 <= res <= 1;

  float a ; 
  if ((ux*vy - uy*vx) < 0)
   a = -acosf(res);
  else
   a = acosf(res);

  return a;
}

static int calcEllipse(float rx, float ry, float ph, int fa, int fs, 
        float x1, float y1, float x2, float y2, 
        float* pcx, float* pcy, float* pth, float* pdth) 
{
    if (!rx || !ry) return -1; 
    rx = (rx > 0) ? rx : -rx;
    ry = (ry > 0) ? ry : -ry;

    float xp, yp; 
    {   
        float tx = (x1 - x2) / 2;
        float ty = (y1 - y2) / 2;

        xp = cosf(ph)*tx + sinf(ph)*ty;
        yp = -sinf(ph)*tx + cosf(ph)*ty;

        float lambda = xp*xp/(rx*rx) + yp*yp/(ry*ry);

        if (lambda > 1) {
            lambda = sqrtf(lambda);

            rx = lambda*rx;
            ry = lambda*ry;
        }
    }

    float cx,cy;
    float th, dth;
    {
        float pr;
        {
            float tp = rx*rx*yp*yp + ry*ry*xp*xp;
            pr = sqrtf((rx*rx*ry*ry - tp)/tp);
        }

        float cxp,cyp;

        if (fa != fs) {
            cxp =  pr * rx * yp / ry;
            cyp = -pr * ry * xp / rx;
        } else {
            cxp = -pr * rx * yp / ry;
            cyp =  pr * ry * xp / rx;
        }
        cx = cosf(ph)*cxp - sinf(ph)*cyp + (x1+x2)/2;
        cy = sinf(ph)*cxp + cosf(ph)*cyp + (y1+y2)/2;

        th = angle(1, 0, (xp - cxp) / rx, (yp - cyp) / ry);
        dth = angle((xp-cxp)/rx,(yp-cyp)/ry,(-xp-cxp)/rx,(-yp-cyp)/ry);
        dth = fmod(dth,2*PI);
    }

    if ( fs==0 && dth > 0)
        dth -= 2*PI;
    if ( fs==1 && dth < 0)
        dth += 2*PI;

    *pcx = cx;
    *pcy = cy;
    *pth = th;
    *pdth = dth;
    return 0;
}
//  EFAULT pathname points outside your accessible address space
static int adrNotValid(void* p)
{
    int fd = open(p, 0, 0);
    int e = errno;
    if (fd == -1 && e == EFAULT)
        return 1;
    else if (fd != -1)
        close(fd);
    return 0;
}

static int _bezierEllipse(float cx, float cy, float rx, float ry,float th, float sth, float eth,float* res)
{
    if(!res || adrNotValid(res) || adrNotValid(&res[5]))
        return -1;

    float x1 = cx + rx * cosf(th)*cosf(sth) - ry * sinf(th)*sinf(sth);
    float y1 = cy + rx * sinf(th)*cosf(sth) + ry * cosf(th)*sinf(sth);

    float x2 = cx + rx * cosf(th)*cosf(eth) - ry * sinf(th)*sinf(eth);
    float y2 = cy + rx * sinf(th)*cosf(eth) + ry * cosf(th)*sinf(eth);

    float dx1 = -1 * rx * cosf(th) * sinf(sth) - ry * sinf(th) * cosf(sth);
    float dy1 = -1 * rx * sinf(th) * sinf(sth) + ry * cosf(th) * cosf(sth);

    float dx2 = -1 * rx * cosf(th) * sinf(eth) - ry * sinf(th) * cosf(eth);
    float dy2 = -1 * rx * sinf(th) * sinf(eth) + ry * cosf(th) * cosf(eth);
    float tmp = tan((eth-sth)/2);
    float alpha = sinf(eth - sth)*(sqrtf(4+3*tmp*tmp)-1)/3;

 //   *res++ = x1;
 //   *res++ = y1;
    *res++ = x1 + alpha*dx1; // p1x
    *res++ = y1 + alpha*dy1; // p1y
    *res++ = x2 - alpha*dx2; // p2x
    *res++ = y2 - alpha*dy2; // p2y
    *res++ = x2;
    *res++ = y2;
    return 0;
}

int bezierEllipse(float cx, float cy, float rx, float ry,float th, float sth, float eth, int n,float* res)
{
    if (n <= 0 || sth == eth) return -1;

    float step = (eth - sth)/n;
    float tmp = sth;
    for (int i = 0; i < n; i++) {
        int e = _bezierEllipse(cx, cy, rx, ry,th, tmp, tmp+step,res);
        if (e) return -1;
        tmp += step;
        res += 6;
    }
    return 0;
}

int main()
{
 float rx = 50;
 float ry = 40;
 float phi = -120.0/180 * PI ;
 float x1 = 100;
 float y1 = 80;
 float x2 = 120;
 float y2 = 60;

 float cx,cy,th,dth;

 calcEllipse( rx, ry, phi, 1, 0, x1, y1, x2, y2,&cx,&cy,&th,&dth);
 printf("%f %f %f %f\n",cx,cy,th,dth);

 int dn = ceilf(((dth > 0)? dth : -dth)  * 4 / PI) ;

 if (dn < 1)
     return 0;

 int size = 2+6*dn;
 float* data = (float*)malloc(size*sizeof(float));

 data[0] = x1;
 data[1] = y1;

 bezierEllipse(cx, cy, rx, ry,phi, th, th+dth, dn,data+2);

 data[size-2] = x2;
 data[size-1] = y2;
 printf("M %f %f ",x1,y1);
 for(int i = 0; i < dn; i++) {
     printf(" C ");
     for(int j = 0; j < 6; j++)
     printf(" %f ",data[2+6*i+j]);
 }
 printf("\n");
}