Java 中的边界椭圆实现
Bounding Ellipse Implementation in Java
我正在为实现一个算法而挠头，我被引导相信该算法将计算方程的系数，该方程将给我一个椭圆，该椭圆将限制一组点。鉴于我不知道该算法实际上是如何完成它应该做的事情的，所以我很困惑为什么该算法实际上没有像我认为的那样工作。我非常确定我已经准确地实现了算法。我意识到我可能已经塞满了某个地方。
我的实现是根据 because I found it easier to work with than the pseudo code given here 建模的。 C++实现的OP说是基于同样的伪代码。
这是我的实现：
// double tolerance = 0.2;
// int n = 10;
// int d = 2;
double tolerance=0.02;
int n=10;
int d=2;
// MatrixXd p = MatrixXd::Random(d,n);
RealMatrix p=new BlockRealMatrix(d,n,new double[][]{{70,56,44,93,77,12,30,51,35,82,74,38,92,49,22,69,71,91,39,13}},false);
// MatrixXd q = p;
// q.conservativeResize(p.rows() + 1, p.cols());
RealMatrix q=p.createMatrix(d+1,n);
q.setSubMatrix(p.getData(),0,0);
// for(size_t i = 0; i < q.cols(); i++)
// {
// q(q.rows() - 1, i) = 1;
// }
//
// const double init_u = 1.0 / (double) n;
// MatrixXd u = MatrixXd::Constant(n, 1, init_u);
double[]ones=new double[n];
double[]uData=new double[n];
for(int i=0;i<n;i++){
ones[i]=1;
uData[i]=((double)1)/((double)n);
}
q.setRow(d,ones);
// int count = 1;
// double err = 1;
int count=0;
double err=1;
while(err>tolerance){
// MatrixXd Q_tr = q.transpose();
RealMatrix qTr=q.transpose();
// MatrixXd X = q * u.asDiagonal() * Q_tr;
RealMatrix uDiag=MatrixUtils.createRealDiagonalMatrix(uData);
RealMatrix qByuDiag=q.multiply(uDiag);
RealMatrix x=qByuDiag.multiply(qTr);
// MatrixXd M = (Q_tr * X.inverse() * q).diagonal();
RealMatrix qTrByxInverse=qTr.multiply(MatrixUtils.inverse(x));
RealMatrix qTrByxInverseByq=qTrByxInverse.multiply(q);
int r=qTrByxInverseByq.getRowDimension();
double mData[]=new double[r];
for(int i=0;i<r;i++){
mData[i]=qTrByxInverseByq.getEntry(i,i);
}
// double maximum = M.maxCoeff(&j_x, &j_y);
// As M is a matrix formed from mData where only cells on the
// diagonal are populated with values greater than zero, the row
// and column values will be identical, and will be equal to the
// place where the maximum value occurs in mData. The matrix M
// is never used again in the algorithm, and hence creation of
// the matrix M is unnecessary.
int iMax=0;
double dMax=0;
for(int i=0;i<mData.length;i++){
if(mData[i]>dMax){
dMax=mData[i];
iMax=i;
}
}
// double step_size = (maximum - d - 1) / ((d + 1) * (maximum + 1));
double stepSize=(dMax-d-1)/((d+1)*(dMax+1));
// MatrixXd new_u = (1 - step_size) * u;
double[]uDataNew=new double[n];
for(int i=0;i<n;i++){
uDataNew[i]=(((double)1)-stepSize)*uData[i];
}
// new_u(j_x, 0) += step_size;
uDataNew[iMax]+=stepSize;
// MatrixXd u_diff = new_u - u;
// for(size_t i = 0; i < u_diff.rows(); i++)
// {
// for(size_t j = 0; j < u_diff.cols(); j++)
// u_diff(i, j) *= u_diff(i, j); // Square each element of the matrix
// }
// err = sqrt(u_diff.sum());
double sum=0;
for(int i=1;i<n;i++){
double cell=uDataNew[i]-uData[i];
sum+=(cell*cell);
}
err=Math.sqrt(sum);
// count++
// u = new_u;
count++;
uData=uDataNew;
}
// MatrixXd U = u.asDiagonal();
RealMatrix uFinal=MatrixUtils.createRealDiagonalMatrix(uData);
// MatrixXd A = (1.0 / (double) d) * (p * U * p.transpose() - (p * u) * (p * u).transpose()).inverse();
// Broken down into the following 9 sub-steps:
// 1 p * u
double[][]uMatrixData=new double[1][];
uMatrixData[0]=uData;
RealMatrix u=new BlockRealMatrix(n,1,uMatrixData,false);
RealMatrix cFinal=p.multiply(u);
// 2 (p * u).transpose()
RealMatrix two=cFinal.transpose();
// 3 (p * u) * (p * u).transpose()
RealMatrix three=cFinal.multiply(two);
// 4 p * U
RealMatrix four=p.multiply(uFinal);
// 5 p * U * p.transpose()
RealMatrix five=four.multiply(p.transpose());
// 6 p * U * p.transpose() - (p * u) * (p * u).transpose()
RealMatrix six=five.subtract(three);
// 7 (p * U * p.transpose() - (p * u) * (p * u).transpose()).inverse()
RealMatrix seven=MatrixUtils.inverse(six);
// 8 1.0 / (double) d
double eight=((double)1)/d;
// 9 MatrixXd A = (1.0 / (double) d) * (p * U * p.transpose() - (p * u) * (p * u).transpose()).inverse()
RealMatrix aFinal=seven.scalarMultiply(eight);
// MatrixXd c = p * u; This has been calculated in sub-step (1) above and stored as cFinal.
System.out.println();
System.out.println("The coefficients of ellipse's equation are given as follows:");
for(int i=0;i<aFinal.getRowDimension();i++){
for(int j=0;j<aFinal.getColumnDimension();j++){
System.out.printf(" %3.8f",aFinal.getEntry(i,j));
}
System.out.println();
}
System.out.println();
System.out.println("The the axis shifts are given as follows:");
for(int i=0;i<cFinal.getRowDimension();i++){
for(int j=0;j<cFinal.getColumnDimension();j++){
System.out.printf(" %3.8f",cFinal.getEntry(i,j));
}
System.out.println();
}
// Get the centre of the set of points, which will be the centre of the
// ellipse. This part was not actually included in the C++
// implementation. I guess the OP considered it too trivial.
double xmin=0;
double xmax=0;
double ymin=0;
double ymax=0;
for(int i=0;i<p.getRowDimension();i++){
double x=p.getEntry(i,0);
double y=p.getEntry(i,1);
if(i==0){
xmin=xmax=x;
ymin=ymax=y;
}else{
if(x<xmin){
xmin=x;
}else if(x>xmax){
xmax=x;
}
if(y<ymin){
ymin=y;
}else if(y>ymax){
ymax=y;
}
}
}
double x=(xmax-xmin)/2+xmin;
double y=(ymax-ymin)/2+ymin;
System.out.println();
System.out.println("The centre of the ellipse is given as follows:");
System.out.printf(" The x axis is %3.8f.\n",x);
System.out.printf(" The y axis is %3.8f.\n",y);
System.out.println();
System.out.println("The algorithm completed ["+count+"] iterations of its while loop.");
// This code constructs and displays a yellow ellipse with a black border.
ArrayList<Integer>pointsx=new ArrayList<>();
ArrayList<Integer>pointsy=new ArrayList<>();
for (double t=0;t<2*Math.PI;t+=0.02){ // <- or different step
pointsx.add(this.getWidth()/2+(int)(cFinal.getEntry(0,0)*Math.cos(t)*aFinal.getEntry(0,0)-cFinal.getEntry(1,0)*Math.sin(t)*aFinal.getEntry(0,1)));
pointsy.add(this.getHeight()/2+(int)(cFinal.getEntry(0,0)*Math.cos(t)*aFinal.getEntry(1,0)+cFinal.getEntry(1,0)*Math.sin(t)*aFinal.getEntry(1,1)));
}
int[]xpoints=new int[pointsx.size()];
Iterator<Integer>xpit=pointsx.iterator();
for(int i=0;xpit.hasNext();i++){
xpoints[i]=xpit.next();
}
int[]ypoints=new int[pointsy.size()];
Iterator<Integer>ypit=pointsy.iterator();
for(int i=0;ypit.hasNext();i++){
ypoints[i]=ypit.next();
}
g.setColor(Color.yellow);
g.fillPolygon(xpoints,ypoints,pointsx.size());
g.setColor(Color.black);
g.drawPolygon(xpoints,ypoints,pointsx.size());
此程序生成以下输出：
The coefficients of ellipse's equation are given as follows:
0.00085538 0.00050693
0.00050693 0.00093474
The axis shifts are given as follows:
54.31114965
55.60647648
The centre of the ellipse is given as follows:
The x axis is 72.00000000.
The y axis is 47.00000000.
The algorithm completed [23] iterations of its while loop.
我觉得有点奇怪，2x2 矩阵的条目非常小。我被引导相信这些条目是用于描述椭圆方向的系数，而第二个 2x1 矩阵描述了椭圆的 x 和 y 轴的大小。
据我了解，用于获取点的方程称为参数方程。他们有一个引用 here.
的表格
椭圆中心的位置和这些值的计算是我加的。它们没有出现在 C++ 实现中，在我将这个算法的输出与用于描述椭圆的参数方程结合之后，我被引导相信 C++ 实现的 OP 给人一种错误的印象，即这个 2x1 矩阵描述了椭圆的中心。我承认我形成的印象可能是错误的，因为如果有人认为我是对的，那么中心（两个轴的最低值和最高值之间的中间位置）似乎是错误的；它小于y轴的大小，我取的是径向测量。
当我将这些值代入椭圆的参数方程以生成我随后用来创建 Polygon 的点时，生成的形状占据一个像素。考虑到描述方向的 2x2 矩阵中给出的值，这是我所期望的。
因此，在我看来，我如何生成产生方向的 2x2 矩阵存在一些问题。
我已尽力简明扼要地提供所有相关事实，以及我形成的任何相关印象，无论它们是对还是错。我希望有人能提供这个问题的答案。
很遗憾，我找不到有关此问题的帮助。
但是，我找到了一个折衷的解决方案，涉及用多种语言实现一个封闭的圆 here。如果可以提供更好的解决方案，我会把这个问题留给其他人回答。
我今天尝试使用您在问题中引用的相同伪代码来实现 Minimum Bounding Ellipse in Java。我必须承认，我不懂数学，但我知道您在 "coefficients of ellipse" 矩阵中看到的超小数字是正常的。
对于我的实现，我通过 Nima Moshtagh 移植了 MatLab 代码。在 MatLab 页面的评论部分，Peter Lawrence 的一些代码看起来像这样：
C=inv(C);
tq=linspace(-pi,pi,M); % for display purposes, M is the number of points on the ellipse
[Ve,De]=eig(C);
De=sqrt(diag(De));
[l1,Ie] = max(De);
veig=Ve(:,Ie);
thu=atan2(veig(2),veig(1));
l2=De(setdiff([1 2],Ie));
U=[cos(thu) -sin(thu);sin(thu) cos(thu)]*[l1*cos(tq);l2*sin(tq)];
plot(U(1,:)+m(1),U(2,:)+m(2))
C 是描述椭圆的 2x2 矩阵，m 是中心点。
结果 U 是一个 x,y 坐标数组。如果将坐标添加到中心点，您将得到一圈可用于渲染椭圆的点。
同样，我不是 100% 确定数学在做什么，但它有效。
我正在为实现一个算法而挠头，我被引导相信该算法将计算方程的系数，该方程将给我一个椭圆，该椭圆将限制一组点。鉴于我不知道该算法实际上是如何完成它应该做的事情的，所以我很困惑为什么该算法实际上没有像我认为的那样工作。我非常确定我已经准确地实现了算法。我意识到我可能已经塞满了某个地方。
我的实现是根据
这是我的实现：
// double tolerance = 0.2;
// int n = 10;
// int d = 2;
double tolerance=0.02;
int n=10;
int d=2;
// MatrixXd p = MatrixXd::Random(d,n);
RealMatrix p=new BlockRealMatrix(d,n,new double[][]{{70,56,44,93,77,12,30,51,35,82,74,38,92,49,22,69,71,91,39,13}},false);
// MatrixXd q = p;
// q.conservativeResize(p.rows() + 1, p.cols());
RealMatrix q=p.createMatrix(d+1,n);
q.setSubMatrix(p.getData(),0,0);
// for(size_t i = 0; i < q.cols(); i++)
// {
// q(q.rows() - 1, i) = 1;
// }
//
// const double init_u = 1.0 / (double) n;
// MatrixXd u = MatrixXd::Constant(n, 1, init_u);
double[]ones=new double[n];
double[]uData=new double[n];
for(int i=0;i<n;i++){
ones[i]=1;
uData[i]=((double)1)/((double)n);
}
q.setRow(d,ones);
// int count = 1;
// double err = 1;
int count=0;
double err=1;
while(err>tolerance){
// MatrixXd Q_tr = q.transpose();
RealMatrix qTr=q.transpose();
// MatrixXd X = q * u.asDiagonal() * Q_tr;
RealMatrix uDiag=MatrixUtils.createRealDiagonalMatrix(uData);
RealMatrix qByuDiag=q.multiply(uDiag);
RealMatrix x=qByuDiag.multiply(qTr);
// MatrixXd M = (Q_tr * X.inverse() * q).diagonal();
RealMatrix qTrByxInverse=qTr.multiply(MatrixUtils.inverse(x));
RealMatrix qTrByxInverseByq=qTrByxInverse.multiply(q);
int r=qTrByxInverseByq.getRowDimension();
double mData[]=new double[r];
for(int i=0;i<r;i++){
mData[i]=qTrByxInverseByq.getEntry(i,i);
}
// double maximum = M.maxCoeff(&j_x, &j_y);
// As M is a matrix formed from mData where only cells on the
// diagonal are populated with values greater than zero, the row
// and column values will be identical, and will be equal to the
// place where the maximum value occurs in mData. The matrix M
// is never used again in the algorithm, and hence creation of
// the matrix M is unnecessary.
int iMax=0;
double dMax=0;
for(int i=0;i<mData.length;i++){
if(mData[i]>dMax){
dMax=mData[i];
iMax=i;
}
}
// double step_size = (maximum - d - 1) / ((d + 1) * (maximum + 1));
double stepSize=(dMax-d-1)/((d+1)*(dMax+1));
// MatrixXd new_u = (1 - step_size) * u;
double[]uDataNew=new double[n];
for(int i=0;i<n;i++){
uDataNew[i]=(((double)1)-stepSize)*uData[i];
}
// new_u(j_x, 0) += step_size;
uDataNew[iMax]+=stepSize;
// MatrixXd u_diff = new_u - u;
// for(size_t i = 0; i < u_diff.rows(); i++)
// {
// for(size_t j = 0; j < u_diff.cols(); j++)
// u_diff(i, j) *= u_diff(i, j); // Square each element of the matrix
// }
// err = sqrt(u_diff.sum());
double sum=0;
for(int i=1;i<n;i++){
double cell=uDataNew[i]-uData[i];
sum+=(cell*cell);
}
err=Math.sqrt(sum);
// count++
// u = new_u;
count++;
uData=uDataNew;
}
// MatrixXd U = u.asDiagonal();
RealMatrix uFinal=MatrixUtils.createRealDiagonalMatrix(uData);
// MatrixXd A = (1.0 / (double) d) * (p * U * p.transpose() - (p * u) * (p * u).transpose()).inverse();
// Broken down into the following 9 sub-steps:
// 1 p * u
double[][]uMatrixData=new double[1][];
uMatrixData[0]=uData;
RealMatrix u=new BlockRealMatrix(n,1,uMatrixData,false);
RealMatrix cFinal=p.multiply(u);
// 2 (p * u).transpose()
RealMatrix two=cFinal.transpose();
// 3 (p * u) * (p * u).transpose()
RealMatrix three=cFinal.multiply(two);
// 4 p * U
RealMatrix four=p.multiply(uFinal);
// 5 p * U * p.transpose()
RealMatrix five=four.multiply(p.transpose());
// 6 p * U * p.transpose() - (p * u) * (p * u).transpose()
RealMatrix six=five.subtract(three);
// 7 (p * U * p.transpose() - (p * u) * (p * u).transpose()).inverse()
RealMatrix seven=MatrixUtils.inverse(six);
// 8 1.0 / (double) d
double eight=((double)1)/d;
// 9 MatrixXd A = (1.0 / (double) d) * (p * U * p.transpose() - (p * u) * (p * u).transpose()).inverse()
RealMatrix aFinal=seven.scalarMultiply(eight);
// MatrixXd c = p * u; This has been calculated in sub-step (1) above and stored as cFinal.
System.out.println();
System.out.println("The coefficients of ellipse's equation are given as follows:");
for(int i=0;i<aFinal.getRowDimension();i++){
for(int j=0;j<aFinal.getColumnDimension();j++){
System.out.printf(" %3.8f",aFinal.getEntry(i,j));
}
System.out.println();
}
System.out.println();
System.out.println("The the axis shifts are given as follows:");
for(int i=0;i<cFinal.getRowDimension();i++){
for(int j=0;j<cFinal.getColumnDimension();j++){
System.out.printf(" %3.8f",cFinal.getEntry(i,j));
}
System.out.println();
}
// Get the centre of the set of points, which will be the centre of the
// ellipse. This part was not actually included in the C++
// implementation. I guess the OP considered it too trivial.
double xmin=0;
double xmax=0;
double ymin=0;
double ymax=0;
for(int i=0;i<p.getRowDimension();i++){
double x=p.getEntry(i,0);
double y=p.getEntry(i,1);
if(i==0){
xmin=xmax=x;
ymin=ymax=y;
}else{
if(x<xmin){
xmin=x;
}else if(x>xmax){
xmax=x;
}
if(y<ymin){
ymin=y;
}else if(y>ymax){
ymax=y;
}
}
}
double x=(xmax-xmin)/2+xmin;
double y=(ymax-ymin)/2+ymin;
System.out.println();
System.out.println("The centre of the ellipse is given as follows:");
System.out.printf(" The x axis is %3.8f.\n",x);
System.out.printf(" The y axis is %3.8f.\n",y);
System.out.println();
System.out.println("The algorithm completed ["+count+"] iterations of its while loop.");
// This code constructs and displays a yellow ellipse with a black border.
ArrayList<Integer>pointsx=new ArrayList<>();
ArrayList<Integer>pointsy=new ArrayList<>();
for (double t=0;t<2*Math.PI;t+=0.02){ // <- or different step
pointsx.add(this.getWidth()/2+(int)(cFinal.getEntry(0,0)*Math.cos(t)*aFinal.getEntry(0,0)-cFinal.getEntry(1,0)*Math.sin(t)*aFinal.getEntry(0,1)));
pointsy.add(this.getHeight()/2+(int)(cFinal.getEntry(0,0)*Math.cos(t)*aFinal.getEntry(1,0)+cFinal.getEntry(1,0)*Math.sin(t)*aFinal.getEntry(1,1)));
}
int[]xpoints=new int[pointsx.size()];
Iterator<Integer>xpit=pointsx.iterator();
for(int i=0;xpit.hasNext();i++){
xpoints[i]=xpit.next();
}
int[]ypoints=new int[pointsy.size()];
Iterator<Integer>ypit=pointsy.iterator();
for(int i=0;ypit.hasNext();i++){
ypoints[i]=ypit.next();
}
g.setColor(Color.yellow);
g.fillPolygon(xpoints,ypoints,pointsx.size());
g.setColor(Color.black);
g.drawPolygon(xpoints,ypoints,pointsx.size());
此程序生成以下输出：
The coefficients of ellipse's equation are given as follows:
0.00085538 0.00050693
0.00050693 0.00093474
The axis shifts are given as follows:
54.31114965
55.60647648
The centre of the ellipse is given as follows:
The x axis is 72.00000000.
The y axis is 47.00000000.
The algorithm completed [23] iterations of its while loop.
我觉得有点奇怪，2x2 矩阵的条目非常小。我被引导相信这些条目是用于描述椭圆方向的系数，而第二个 2x1 矩阵描述了椭圆的 x 和 y 轴的大小。
据我了解，用于获取点的方程称为参数方程。他们有一个引用 here.
的表格椭圆中心的位置和这些值的计算是我加的。它们没有出现在 C++ 实现中，在我将这个算法的输出与用于描述椭圆的参数方程结合之后，我被引导相信 C++ 实现的 OP 给人一种错误的印象，即这个 2x1 矩阵描述了椭圆的中心。我承认我形成的印象可能是错误的，因为如果有人认为我是对的，那么中心（两个轴的最低值和最高值之间的中间位置）似乎是错误的；它小于y轴的大小，我取的是径向测量。
当我将这些值代入椭圆的参数方程以生成我随后用来创建 Polygon 的点时，生成的形状占据一个像素。考虑到描述方向的 2x2 矩阵中给出的值，这是我所期望的。
因此，在我看来，我如何生成产生方向的 2x2 矩阵存在一些问题。
我已尽力简明扼要地提供所有相关事实，以及我形成的任何相关印象，无论它们是对还是错。我希望有人能提供这个问题的答案。
很遗憾，我找不到有关此问题的帮助。
但是，我找到了一个折衷的解决方案，涉及用多种语言实现一个封闭的圆 here。如果可以提供更好的解决方案，我会把这个问题留给其他人回答。
我今天尝试使用您在问题中引用的相同伪代码来实现 Minimum Bounding Ellipse in Java。我必须承认，我不懂数学，但我知道您在 "coefficients of ellipse" 矩阵中看到的超小数字是正常的。
对于我的实现，我通过 Nima Moshtagh 移植了 MatLab 代码。在 MatLab 页面的评论部分，Peter Lawrence 的一些代码看起来像这样：
C=inv(C);
tq=linspace(-pi,pi,M); % for display purposes, M is the number of points on the ellipse
[Ve,De]=eig(C);
De=sqrt(diag(De));
[l1,Ie] = max(De);
veig=Ve(:,Ie);
thu=atan2(veig(2),veig(1));
l2=De(setdiff([1 2],Ie));
U=[cos(thu) -sin(thu);sin(thu) cos(thu)]*[l1*cos(tq);l2*sin(tq)];
plot(U(1,:)+m(1),U(2,:)+m(2))
C 是描述椭圆的 2x2 矩阵，m 是中心点。
结果 U 是一个 x,y 坐标数组。如果将坐标添加到中心点，您将得到一圈可用于渲染椭圆的点。
同样，我不是 100% 确定数学在做什么，但它有效。