c++ - Bounding Ellipse Implementation in Java -

i'm scratching head implementation of algorithm have been led believe calculate coefficients equation give me ellipse bound set of points. given don't know how algorithm supposed do, i'm stumped why algorithm isn't working think should. sure can have implemented algorithm accurately. realise have stuffed somewhere.

my implementation modelled this implementation in c++ because found easier work pseudo code given here. op of c++ implementation says based on same pseudo code.

here implementation:

// 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 = 0; < 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);     // m matrix formed mdata cells on     // diagonal populated values greater zero, row     // , column values identical, , equal     // place maximum value occurs in mdata. matrix m     // never used again in algorithm, , hence creation of     // matrix m 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 = 0; < 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 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 = (1.0 / (double) d) * (p * u * p.transpose() - (p * u) * (p * u).transpose()).inverse(); // broken down 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 = (1.0 / (double) d) * (p * u * p.transpose() - (p * u) * (p * u).transpose()).inverse() realmatrix afinal=seven.scalarmultiply(eight);  // matrixxd c = p * u; has been calculated in sub-step (1) above , stored cfinal.  system.out.println(); system.out.println("the coefficients of ellipse's equation given 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 axis shifts given 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(); }  // centre of set of points, centre of // ellipse. part not included in c++ // implementation. guess op considered 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 ellipse given follows:"); system.out.printf(" x axis %3.8f.\n",x); system.out.printf(" y axis %3.8f.\n",y);  system.out.println(); system.out.println("the algorithm completed ["+count+"] iterations of while loop.");  // code constructs , displays yellow ellipse black border.  arraylist<integer>pointsx=new arraylist<>(); arraylist<integer>pointsy=new arraylist<>(); (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()); 

this program generates following output:

the coefficients of ellipse's equation given follows:   0.00085538  0.00050693   0.00050693  0.00093474  axis shifts given follows:   54.31114965   55.60647648  centre of ellipse given follows:  x axis 72.00000000.  y axis 47.00000000.  algorithm completed [23] iterations of while loop. 

i find strange entries of 2x2 matrix small. led believe these entries coefficients used describe orientation of ellipse, while second 2x1 matrix describes size of ellipse's x , y axes.

i understand equations used obtain points referred parametric equations. have form quoted here.

the location of centre of ellipse , computation of these values has been added me. not appear in c++ implementation, , after married output of algorithm parametric equations used describe ellipse, lead believe op of c++ implementation gave wrong impression 2x1 matrix described ellipse's centre. admit impression formed wrong, because if 1 assumes right, centre (half-way between lowest , highest values of both axes) appears wrong; less size of y axis, take radial measure.

when plug these values parametric equations ellipse generate points use create polygon, generated shape occupies single pixel. considering values given in 2x2 matrix describing orientation, expect.

hence, seems me there problem in how generate 2x2 matrix produces orientation.

i have done best concise , provide relevant facts, relevant impressions have formed whether right or wrong. hope can provide answer question.

unfortunately, haven't been able find question.

however, able find compromise solution involving implementations in multiple languages enclosing circle here. i'll leave question else answer if better solution can offered.