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.1415926static 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);elsea = 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");
}