标签：yi dim 椭圆 mat int double xi 拟合 乘法
转自：https://blog.csdn.net/weixin_39591047/article/details/87542496
参考文章：
最小二乘法拟合椭圆——MATLAB和Qt-C++实现
https://blog.csdn.net/sinat_21107433/article/details/80877758
以上文章中，C++代码有问题。因此参考如下文章，得到正确的结果。
矩阵求逆-高斯消元法介绍及其实现
https://blog.csdn.net/qithon/article/details/80100029
代码如下：
1 import android.util.Log;
2
3 public class EclipseFitting {
4
5 /**
6 *
7 * 功能说明： 对平面上的一些列点给出最小二乘的椭圆拟合，利用消元法
8 * 解得最小二乘解作为椭圆参数。
9 *
10 *调用形式： EclipseFitting2（arrayx,arrayy,Result）；
11 *参数说明： arrayx: arrayx[n] 为第n个点的x坐标
12 arrayy: arrayy[n]为第n点的y坐标
13 n : 点的个数
14 Result : Result[0][],椭圆的五个参数，分别为center.x,center.y,2a,2b,xtheta
15 Result[1][],椭圆方程的的五个系数（A,B,C,D,E），X^2+A*XY+B*Y^2+C*X+D*Y+E=0
16 esp: 解精度，通常取1e-6,这个是解方程用的说
17 *
18 *
19 * @param arrayx
20 * @param arrayy
21 * @param n
22 * @param Result
23 */
24
25 void EclipseFitting2(int [] arrayx, int [] arrayy,int N,double [][] Result)
26 {
27 double A = 0.00,B = 0.00,C = 0.00,D = 0.00,E = 0.00;
28 double x2y2=0.0,x1y3=0.0,x2y1=0.0,x1y2=0.0,x1y1=0.0,yyy4=0.0,yyy3=0.0,yyy2=0.0,xxx2=0.0,xxx1=0.0,yyy1=0.0,x3y1=0.0,xxx3=0.0;
29 //平面上任意位置的一个椭圆，其中心坐标为（x0，y0），半长轴a，半短轴b，长轴偏角为θ，方程通式为
30 //方程为：X^2+A*XY+B*Y^2+C*X+D*Y+E=0
31
32 for(int i = 0; i < N; i++)
33 {
34 double xi = arrayx[i], yi = arrayy[i];
35
36 x2y2 += xi*xi*yi*yi;
37 x1y3 += xi*yi*yi*yi;
38 x2y1 += xi*xi*yi;
39 x1y2 += xi*yi*yi;
40 x1y1 += xi*yi;
41 yyy4 += yi*yi*yi*yi;
42 yyy3 += yi*yi*yi;
43 yyy2 += yi*yi;
44 xxx2 += xi*xi;
45 xxx1 += xi;
46 yyy1 += yi;
47 x3y1 += xi*xi*xi*yi;
48 xxx3 += xi*xi*xi;
49 }
50
51 //5*5 matrix
52
53
54 double[][] matrix={{x2y2,x1y3,x2y1,x1y2,x1y1},
55 {x1y3,yyy4,x1y2,yyy3,yyy2},
56 {x2y1,x1y2,xxx2,x1y1,xxx1},
57 {x1y2,yyy3,x1y1,yyy2,yyy1},
58 {x1y1,yyy2,xxx1,yyy1,N}
59 };
60
61 double[] matrix2={x3y1,x2y2,xxx3,x2y1,xxx2};
62 double[] matrix3 = {A,B,C,D,E};
63
64 Log.i("EclipseFitting2","系数矩阵");
65 LogMatrix(matrix,5);
66
67
68 求矩阵matrix的逆，结果为InverseMatrix
69
70 ///
71 double[][] mInverseMatrix=InverseMatrix(matrix,5);
72
73 ///求参数A,B,C,D,E
74 for(int i=0;i<5;i++)
75 {
76 for(int j=0;j<5;j++)
77 {
78 matrix3[i] += mInverseMatrix[i][j]*(-matrix2[j]);
79 }
80 }
81 A = matrix3[0];
82 B = matrix3[1];
83 C = matrix3[2];
84 D = matrix3[3];
85 E = matrix3[4];
86
87 ///求拟合结果重要参数
88
89 double Xc = (2*B*C-A*D)/(A*A-4*B);
90 double Yc = (2*D-A*D)/(A*A-4*B);
91 double a = Math.sqrt(Math.abs(2*(A*C*D-B*C*C-D*D+4*B*E-A*A*E)/((A*A-4*B)*(B-Math.sqrt(A*A+(1-B)*(1-B))+1))));
92 double b = Math.sqrt(Math.abs(2*(A*C*D-B*C*C-D*D+4*B*E-A*A*E)/((A*A-4*B)*(B+Math.sqrt(A*A+(1-B)*(1-B))+1))));
93 double theta = Math.atan2(a*a-b*b*B,a*a*B-b*b);
94
95 Result[0][0]=Xc; Result[0][1]=Yc;
96 Result[0][2]=a; Result[0][3]=b;
97 Result[0][4]=theta*180/3.1415926;
98
99 Result[1][0]=A; Result[1][1]=B; Result[1][2]=C;
100 Result[1][3]=D; Result[1][4]=E;
101
102 }
103
104 /**
105 matrix为所求矩阵
106 dim为矩阵的维度
107 */
108 public static double[][] InverseMatrix(double[][] matrix, int dim)
109 {
110 double eps = 1e-6;
111 double[][] mat = new double [dim][ dim * 2];
112 //构造增广矩阵W
113 for(int i = 0;i < dim; i++)
114 {
115 for(int j = 0;j < 2 * dim; j++)
116 {
117 if(j < dim)
118 {
119 mat[i][j] = matrix[i][j];
120 }
121 else
122 {
123 mat[i][j] = (j - dim == i) ? 1 : 0;
124 }
125 }
126 }
127
128
129
130 //从上往下消元
131 for(int i = 0;i < dim; i++)
132 {
133 if(Math.abs(mat[i][i]) < eps)
134 //若mat[i,i]为0，则往下找一个不为0的行，加到第i行
135 {
136 int j;
137 for ( j = i + 1; j < dim; j++)
138 {
139 if (Math.abs(mat[j][ i]) > eps) break;
140 }
141 if (j == dim) return null;
142 for(int r = i; r < 2 * dim; r++)
143 {
144 mat[i][ r] += mat[j][ r];
145 }
146 }
147
148
149 double ep = mat[i][ i];
150 //将mat[i][i]变为1
151 for (int r = i; r < 2 * dim; r++)
152 {
153 mat[i][ r] /= ep;
154 }
155
156 for(int j = i + 1; j < dim; j++)
157 {
158 double e = -1 * (mat[j][ i] / mat[i][ i]);
159 for(int r = i; r < 2 * dim; r++)
160 {
161 mat[j][ r] += e * mat[i][ r];
162 }
163 }
164 }
165
166 LogMatrix(mat,10);
167
168 //从下往上消元
169 for(int i = dim - 1; i >= 0; i--)
170 {
171 for (int j = i - 1; j >= 0; j--)
172 {
173 double e = -1 * (mat[j][ i] / mat[i][ i]);
174 for (int r = i; r < 2 * dim; r++)
175 {
176 mat[j][ r] += e * mat[i][ r];
177 }
178 }
179 }
180
181
182 LogMatrix(mat,10);
183
184
185 double[][] result = new double[dim][ dim];
186 for(int i = 0;i < dim; i++)
187 {
188 for(int r = dim; r < 2 * dim; r++)
189 {
190 result[i][ r - dim] = mat[i][ r];
191 }
192 }
193
194 return result;
195 }
196
197
198 // 原文：https://blog.csdn.net/qithon/article/details/80100029
199
200
201
202
203
204 private static void LogMatrix(double[][] Matrix, int width){
205
206 for(int i=0;i<5;i++){
207 String information=" ";
208
209 for(int j=0;j<width;j++){
210 information+=Matrix[i][j]+",";
211 }
212
213
214 Log.i("EclipseFitting2",information);
215
216 }
217
218 }
219
220
221
222 // 原文：https://blog.csdn.net/sinat_21107433/article/details/80877758
223
224
225 }
标签：yi,dim,椭圆,mat,int,double,xi,拟合,乘法 来源： https://www.cnblogs.com/peifx/p/16610507.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。