2017-08-18 110 views
0

我在执行一个算法时正在摸索我的头脑,我相信会计算一个方程的系数,这个方程会给我一个将限定一组点的椭圆。鉴于我不知道该算法如何实际执行它应该做的事情,我很难理解为什么该算法实际上并不像我认为的那样工作。我确信我已经准确地实现了该算法。我意识到我可能已经塞满了某个地方。在Java中实现边界椭圆

我的实现是从this implementation in C++建模的,因为我发现它比伪代码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. 

我觉得有些奇怪的是,2×2矩阵的输入是非常小的。我被引导认为这些条目是用于描述椭圆的方向的系数,而第二个2×1矩阵描述了椭圆的x和y轴的大小。

我知道用于获得点的方程式称为参数方程。他们有一个表格,引用here

椭圆的中心位置和这些值的计算已经被我添加了。它们没有出现在C++实现中,并且在将此算法的输出结合到用于描述椭圆的参数方程之后,我导致相信C++实现的OP给了错误的印象,即这个2×1矩阵描述了椭圆的中心。我承认我形成的印象可能是错误的,因为如果一个人认为我是对的,那么这个中心(两个轴的最低值和最高值之间的中间值)似乎是错误的;它小于y轴的大小,我采取的是径向测量。

当我将这些值插入椭圆的参数方程来生成点I然后使用创建一个Polygon时,生成的形状占据一个像素。考虑到描述方向的2×2矩阵中给出的值,这就是我所期望的。

因此,在我看来,我如何生成产生方向的2x2矩阵存在一些问题。

我已尽力简洁,并提供所有相关事实,以及我形成的任何相关印象,无论它们是对还是错。我希望有人能提供这个问题的答案。

+0

请:编辑你的问题和工作,最后一段。关注这个问题。请注意:当人们告诉你“问题不好”时,那么你的反应不应该是“无论如何,请继续前进”。因为那你就忽略了有价值的反馈。 – GhostCat

+0

听起来像一个好主意。感谢您的忠告,我会编辑我的问题并继续前进。 – 000

回答

0

不幸的是,我一直无法找到这个问题的帮助。

但是,我找到了一个折衷解决方案,涉及多种语言的实现,用于包含圆圈here。如果可以提供更好的解决方案,我会让其他人回答这个问题。