基于C/C++的椭圆等间距离散方法

参考文档与博客：
https://blog.csdn.net/baidu_38621657/article/details/87907829

https://blog.csdn.net/baidu_38621657/article/details/87907829

理论部分

由于椭圆弧长积分的解析解不存在，用数值分析方法近似求解。
下面以第一象限椭圆弧为例说明一种近似求解方法

（1）首先圆是椭圆的特殊形式，长短轴长相等，对于圆心角为θ的圆弧长为Rθ，当圆心角足够小时可以认为tanθ≈θ，弧长可以认为L≈Rtanθ≈R*θ，其中R为圆心到微段弧上点的距离，这个近似简化对于椭圆弧同样有效。

（2）对于椭圆，其参数方程为

x=acosθ
y=bsinθ
这里a,b分别对应椭圆外切圆与内接圆的半径（与下图对应）。

在此约定θ为逆时针方向与X轴之间的夹角

图中点A对应的参数角θ如图所示，α为其对应的扫略角，θ和α之间可以相互转换。

圆心O到点P之间距离|OP|=sqrt(x² + y²) =sqrt[(acosθ)²+(bsinθ)²)]，可得θ与|OP|的关系。

（3）当a=50,b=30时，得|OP|在不同θ时值
以θ为x轴，|OP|为y轴时，曲线图近似为折线段如图中绿色线curve1

绿色线curve1与x轴之间面积即为θ为0 ~ π/2时对应椭圆弧长，此时可将绿色线curve1近似为图中白色直线段curve2（为了更精确也可简化为连续折线段），即白色线curve2与x轴之间面积A=（a+b)*(π/4)(梯形面积公式),近似为θ为0~PI/2时对应椭圆弧长。

此时将椭圆弧按长度等分即可转化为将面积A按θ等分问题
curve2斜率为k=(b-a)/(π/2),方程为y=(2b-2a)/π*x+a,其与x轴之间面积积分为

∫ydx=∫[(2b-2a)/πx+a]dx=(b-a)/πx²+a*x+c,其中c为常数项，

令f(x)=(b-a)/πx²+ax，若将椭圆弧等分位n份，每份面积为A/n
f(θ1)-f(0)=f(θ1)=A/n，
f(θ2)-f(θ1)=A/n,
…
f(θn)-f(θn-1)=A/n

构造方程l来求解每一个点所对应的旋转角度：
((a-b)/π)x²+bx-4Aj/n
这里的j值的范围从1到离散坐标总数。

其次椭圆的周长计算参考：
https://www.mathsisfun.com/geometry/ellipse-perimeter.html

C++代码实现

求解各个点所对应的θ（求解方程）

root roots(double r1, double r2,int j, int N)
{root solu;double a, b, c;a = (r1 - r2) / pi;b = r2;double A = (r1 + r2) * pi / 4;c = -4 * A / N * (j+1);double discriminant = b * b - 4 * a * c;if (discriminant > 0) {solu.m_x1 = (-b + sqrt(discriminant)) / (2 * a);solu.m_x2 = (-b - sqrt(discriminant)) / (2 * a);}else if (discriminant == 0) {//cout << "Roots are real and same." << endl;solu.m_x1 = -b / (2 * a);cout << "x1 = x2 =" << solu.m_x1 << endl;}else {//realPart = -b / (2 * a);//imaginaryPart = sqrt(-discriminant) / (2 * a);cout << "Roots are complex and different." << endl;//cout << "x1 = " << realPart << "+" << imaginaryPart << "i" << endl;//cout << "x2 = " << realPart << "-" << imaginaryPart << "i" << endl;}return solu;
}

椭圆离散函数：

void discretizeEllipse(//椭圆离散化为多边形double x, double y, double r1, double r2, double alpha, double L)//x,y为中心点坐标，r1长轴长度，r2短轴长度，alpha椭圆旋转角度，L离散的线段长度
{double P = 2 * pi * sqrt((r1*r1 + r2*r2) / 2);double N = P / L; //计算需要划分的段数N = 4 * round((double)N / 4);//将段数转换为4的倍数'if (N == 0)N = 4;for (int j = 0; j < N / 4; j++){root rts;rts = roots(r1,r2,j,N);root_out.push_back(rts.m_x1);}for (int k = 0; k < N; k++){if (k < N / 4){p_out.push_back(root_out[k]);}else if ((k >= N / 4) && (k < N / 2 - 1)){p_out.push_back(pi - root_out[N / 2 - k - 2]);}else if (k == N / 2 - 1){p_out.push_back(3.1415);}else if ((k > N / 2 - 1) && k < N * 3 / 4){p_out.push_back(pi + root_out[k - (N / 2)]);}else if (k >= N * 3 / 4 && k < N - 1){p_out.push_back(2 * pi - root_out[N - k - 2]);}else{p_out.push_back(2 * pi);}}for (int m = 0; m < N; m++){double x1 = r1 * cos(p_out[m]);double y1 = r2 * sin(p_out[m]);double x2 = x1 * cos(alpha) - y1 * sin(alpha);double y2 = x1 * sin(alpha) + y1 * cos(alpha);x2 = x2 + x;y2 = y2 + y;vec_out.push_back(x2);vec_out.push_back(y2);}
}

完整代码

#include <algorithm>
#include <graphics.h>
#include<vector>
#include <fstream>
#include<math.h>
#include<iostream>
#define pi 3.1415926using namespace std;
class root
{
public:  //访问权限  公共的权限//属性double m_x1;//根double m_x2;
};
vector <double> root_out;//存储根
vector <double> p_out; //存储整合后的根
vector <double> vec_out;//坐标输出root roots(double r1, double r2,int j, int N)
{root solu;double a, b, c;a = (r1 - r2) / pi;b = r2;double A = (r1 + r2) * pi / 4;c = -4 * A / N * (j+1);double discriminant = b * b - 4 * a * c;if (discriminant > 0) {solu.m_x1 = (-b + sqrt(discriminant)) / (2 * a);solu.m_x2 = (-b - sqrt(discriminant)) / (2 * a);}else if (discriminant == 0) {//cout << "Roots are real and same." << endl;solu.m_x1 = -b / (2 * a);cout << "x1 = x2 =" << solu.m_x1 << endl;}else {//realPart = -b / (2 * a);//imaginaryPart = sqrt(-discriminant) / (2 * a);cout << "Roots are complex and different." << endl;//cout << "x1 = " << realPart << "+" << imaginaryPart << "i" << endl;//cout << "x2 = " << realPart << "-" << imaginaryPart << "i" << endl;}return solu;
}void discretizeEllipse(//椭圆离散化为多边形double x, double y, double r1, double r2, double alpha, double L)//x,y为中心点坐标，r1长轴长度，r2短轴长度，alpha椭圆旋转角度，L离散的线段长度
{double P = 2 * pi * sqrt((r1*r1 + r2*r2) / 2);//float P = 2 * pi * sqrtf((r1^2 + r2 ^ 2) / 2);double N = P / L;   N = 4 * round((double)N / 4);if (N == 0)N = 4;for (int j = 0; j < N / 4; j++){root rts;rts = roots(r1,r2,j,N);root_out.push_back(rts.m_x1);}for (int k = 0; k < N; k++){if (k < N / 4){p_out.push_back(root_out[k]);}else if ((k >= N / 4) && (k < N / 2 - 1)){p_out.push_back(pi - root_out[N / 2 - k - 2]);}else if (k == N / 2 - 1){p_out.push_back(3.1415);}else if ((k > N / 2 - 1) && k < N * 3 / 4){p_out.push_back(pi + root_out[k - (N / 2)]);}else if (k >= N * 3 / 4 && k < N - 1){p_out.push_back(2 * pi - root_out[N - k - 2]);}else{p_out.push_back(2 * pi);}}for (int m = 0; m < N; m++){double x1 = r1 * cos(p_out[m]);double y1 = r2 * sin(p_out[m]);double x2 = x1 * cos(alpha) - y1 * sin(alpha);double y2 = x1 * sin(alpha) + y1 * cos(alpha);x2 = x2 + x;y2 = y2 + y;vec_out.push_back(x2);vec_out.push_back(y2);}
//输出至文件ofstream ofs;ofs.open("test1.txt", ios::out);ofs << vec_out.size() << endl;for (int i = 0; i < vec_out.size(); i++){ofs << vec_out[i] << endl;}ofs.close();
}
void main()
{discretizeEllipse(0,0,50,40,0,20);system("pause");
}

matlab绘制离散结果