问Java中的边界椭圆实现

EN

// 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.