Implémentation d'Ellipse Bounding en Java

Question

Je me gratte la tête avec la mise en œuvre d'un algorithme que l'on m'a laissé supposer calculer les coefficients d'une équation qui me donnera une ellipse qui va lier un ensemble de points. Étant donné que je ne sais pas comment l'algorithme fait réellement ce qu'il est censé faire, je suis perplexe quant à la raison pour laquelle l'algorithme ne fonctionne pas comme je le pense. Je suis aussi sûr que je peux être que j'ai implémenté l'algorithme avec précision. Je me rends compte que je me suis probablement bourré quelque part.

Ma mise en œuvre a été modélisée à partir de this implementation in C++ car je l'ai trouvé plus facile à utiliser que le pseudo code given here. L'OP de l'implémentation C++ dit qu'il est basé sur le même pseudo code.

Voici mon implémentation:

// 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()); 

Ce programme génère la sortie suivante:

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. 

Je trouve un peu étrange que les entrées de la matrice 2x2 sont très petites. Je suis amené à croire que ces entrées sont les coefficients utilisés pour décrire l'orientation de l'ellipse, tandis que la deuxième matrice 2x1 décrit la taille des axes x et y de l'ellipse. Je comprends que les équations utilisées pour obtenir les points sont appelées équations paramétriques. Ils ont un formulaire qui est coté here.

L'emplacement du centre de l'ellipse et le calcul de ces valeurs a été ajouté par moi. Ils n'apparaissent pas dans l'implémentation C++, et après avoir épousé la sortie de cet algorithme aux équations paramétriques utilisées pour décrire l'ellipse, j'ai été amené à croire que l'OP de l'implémentation C++ donnait la fausse impression que cette matrice 2x1 décrivait centre de l'ellipse. J'avoue que l'impression que j'ai formée pourrait être fausse, parce que si l'on suppose que j'ai raison, alors le centre (à mi-chemin entre les valeurs les plus basses et les plus hautes des deux axes) semble être faux; il est inférieur à la taille de l'axe des y, que je prends est une mesure radiale.

Lorsque je branche ces valeurs dans les équations paramétriques de l'ellipse pour générer les points que j'utilise ensuite pour créer un Polygon, la forme générée occupe un seul pixel. Considérant les valeurs données dans la matrice 2x2 décrivant l'orientation, c'est ce à quoi je m'attendais.

Par conséquent, il me semble qu'il y a un problème dans la façon dont je génère la matrice 2x2 qui produit l'orientation.

J'ai fait de mon mieux pour être concis et pour fournir tous les faits pertinents, ainsi que toutes les impressions pertinentes que j'ai formées, qu'elles soient justes ou fausses. J'espère que quelqu'un peut fournir une réponse à cette question.

Answer 1

Malheureusement, je n'ai pas trouvé d'aide pour cette application.

Cependant, j'ai été capable de trouver une solution de compromis impliquant des implémentations dans plusieurs langues pour un cercle englobant here. Je vais laisser cette question à quelqu'un d'autre pour répondre si une meilleure solution peut être offerte.