java 边界椭圆

声明:本页面是StackOverFlow热门问题的中英对照翻译,遵循CC BY-SA 4.0协议,如果您需要使用它,必须同样遵循CC BY-SA许可,注明原文地址和作者信息,同时你必须将它归于原作者(不是我):StackOverFlow 原文地址: http://stackoverflow.com/questions/1768197/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me): StackOverFlow

提示:将鼠标放在中文语句上可以显示对应的英文。显示中英文
时间:2020-10-29 17:53:01  来源:igfitidea点击:

Bounding ellipse

javageometrycollision-detectionbounding

提问by Martin

I have been given an assignement for a graphics module, one part of which is to calculate the minimum bounding ellipse of a set of arbitrary shapes. The ellipse doesn't have to be axis aligned.

我得到了一个图形模块的作业,其中一部分是计算一组任意形状的最小边界椭圆。椭圆不必与轴对齐。

This is working in java (euch) using the AWT shapes, so I can use all the tools shape provides for checking containment/intersection of objects.

这是使用 AWT 形状在 java (euch) 中工作的,因此我可以使用 shape 提供的所有工具来检查对象的包含/相交。

回答by Jacob

You're looking for the Minimum Volume Enclosing Ellipsoid, or in your 2D case, the minimum area. This optimization problem is convex and can be solved efficiently. Check out the MATLAB code in the link I've included - the implementation is trivial and doesn't require anything more complex than a matrix inversion.

您正在寻找最小体积封闭椭圆体,或者在您的 2D 情况下,最小面积。这个优化问题是凸的,可以有效地解决。查看我包含的链接中的 MATLAB 代码 - 实现很简单,不需要比矩阵求逆更复杂的东西。

Anyone interested in the math should read this document.

任何对数学感兴趣的人都应该阅读这份文件

Also, plotting the ellipse is also simple - this can be found here, but you'll need a MATLAB-specific function to generate points on the ellipse.

此外,绘制椭圆也很简单 - 可以在此处找到,但您需要一个特定于 MATLAB 的函数来在椭圆上生成点。

Butsince the algorithm returns the equation of the ellipse in the matrix form,

由于算法以矩阵形式返回椭圆的方程,

matrix form

matrix form

you can use this code to see how you can convert the equation to the canonical form,

您可以使用此代码查看如何将方程转换为规范形式,

canonical

canonical

using Singular Value Decomposition (SVD). And then it's quite easy to plot the ellipse using the canonical form.

使用奇异值分解(SVD)。然后很容易使用规范形式绘制椭圆。

Here's the result of the MATLAB code on a set of 10 random 2D points (blue). results

这是 MATLAB 代码在一组 10 个随机二维点(蓝色)上的结果。 results

Other methods like PCAdoes not guarantee that the ellipse obtained from the decomposition (eigen/singular value) will be minimum bounding ellipse since points outside the ellipse is an indication of the variance.

PCA等其他方法不保证从分解(特征值/奇异值)中获得的椭圆将是最小边界椭圆,因为椭圆外的点是方差的指示。

EDIT:

编辑:

So if anyone read the document, there are two ways to go about this in 2D: here's the pseudocode of the optimal algorithm - the suboptimal algorithm is clearly explained in the document:

因此,如果有人阅读了该文档,那么在 2D 中有两种方法可以解决此问题:这是最佳算法的伪代码 - 文档中清楚地解释了次优算法:

Optimal algorithm:

最优算法:

Input: A 2x10 matrix P storing 10 2D points 
       and tolerance = tolerance for error.
Output: The equation of the ellipse in the matrix form, 
        i.e. a 2x2 matrix A and a 2x1 vector C representing 
        the center of the ellipse.

// Dimension of the points
d = 2;   
// Number of points
N = 10;  

// Add a row of 1s to the 2xN matrix P - so Q is 3xN now.
Q = [P;ones(1,N)]  

// Initialize
count = 1;
err = 1;
//u is an Nx1 vector where each element is 1/N
u = (1/N) * ones(N,1)       

// Khachiyan Algorithm
while err > tolerance
{
    // Matrix multiplication: 
    // diag(u) : if u is a vector, places the elements of u 
    // in the diagonal of an NxN matrix of zeros
    X = Q*diag(u)*Q'; // Q' - transpose of Q    

    // inv(X) returns the matrix inverse of X
    // diag(M) when M is a matrix returns the diagonal vector of M
    M = diag(Q' * inv(X) * Q); // Q' - transpose of Q  

    // Find the value and location of the maximum element in the vector M
    maximum = max(M);
    j = find_maximum_value_location(M);

    // Calculate the step size for the ascent
    step_size = (maximum - d -1)/((d+1)*(maximum-1));

    // Calculate the new_u:
    // Take the vector u, and multiply all the elements in it by (1-step_size)
    new_u = (1 - step_size)*u ;

    // Increment the jth element of new_u by step_size
    new_u(j) = new_u(j) + step_size;

    // Store the error by taking finding the square root of the SSD 
    // between new_u and u
    // The SSD or sum-of-square-differences, takes two vectors 
    // of the same size, creates a new vector by finding the 
    // difference between corresponding elements, squaring 
    // each difference and adding them all together. 

    // So if the vectors were: a = [1 2 3] and b = [5 4 6], then:
    // SSD = (1-5)^2 + (2-4)^2 + (3-6)^2;
    // And the norm(a-b) = sqrt(SSD);
    err = norm(new_u - u);

    // Increment count and replace u
    count = count + 1;
    u = new_u;
}

// Put the elements of the vector u into the diagonal of a matrix
// U with the rest of the elements as 0
U = diag(u);

// Compute the A-matrix
A = (1/d) * inv(P * U * P' - (P * u)*(P*u)' );

// And the center,
c = P * u;

回答by Peter

Thanks to Jacob's pseudocode I was able to implement the Minimum Volume Enclosing Ellipsoid (MVEE) in Java. There are public methods to get the center point, the "A" matrix, and a method to generate a list of coordinates that can be used to render the ellipse. The latter is based on MatLab code posted by Peter Lawrence in the comments to the original MVEE code. Note that the code references a class called "Eigen" - a modified version of Jama's EigenvalueDecomposition class (I took out Matrix class dependencies). I would add it but there is a 30k character limit to answers...

感谢 Jacob 的伪代码,我能够在 Java 中实现最小体积封闭椭圆体 (MVEE)。有一些公共方法可以获取中心点、“A”矩阵,以及一种生成可用于渲染椭圆的坐标列表的方法。后者基于 Peter Lawrence 在对原始 MVEE代码的评论中发布的 MatLab代码。请注意,代码引用了一个名为“Eigen”的类——Jama 的 EigenvalueDecomposition 类的修改版本(我去掉了 Matrix 类依赖项)。我会添加它,但答案有 30k 个字符的限制......

public class Ellipse {

    private double[] center;
    private double[][] A;
    private double l1;
    private double l2;
    private double thu;

  //**************************************************************************
  //** Constructor
  //**************************************************************************
  /** @param P An array of points. Each entry in the array contains an x,y
   *  coordinate.
   */
    public Ellipse(double[][] P, double tolerance){

         // Dimension of the points
        double d = 2;

        // Number of points
        int N = P.length;

        // Rotate the array of points
        P = transpose(P);


        // Add a row of 1s to the 2xN matrix P - so Q is 3xN now.
        //Q = [P;ones(1,N)]
        double[][] Q = merge(P, ones(1,N));


        // Initialize
        int count = 1;
        double err = 1;


        //u is an Nx1 vector where each element is 1/N
        //u = (1/N) * ones(N,1)
        double[] u = new double[N];
        for (int i=0; i<u.length; i++) u[i] = (1D/(double)N);




        // Khachiyan Algorithm
        while (err > tolerance){


            // Matrix multiplication:
            // diag(u) : if u is a vector, places the elements of u
            // in the diagonal of an NxN matrix of zeros
            //X = Q*diag(u)*Q'; // Q' - transpose of Q
            double[][] X = multiply(multiply(Q,diag(u)), transpose(Q));

            // inv(X) returns the matrix inverse of X
            // diag(M) when M is a matrix returns the diagonal vector of M
            //M = diag(Q' * inv(X) * Q); // Q' - transpose of Q
            double[] M = diag(multiply(multiply(transpose(Q), inv(X)), Q));


            //Find the value and location of the maximum element in the vector M
            double maximum = max(M);
            int j = find_maximum_value_location(M, maximum);


            // Calculate the step size for the ascent
            double step_size = (maximum - d -1)/((d+1)*(maximum-1));


            // Calculate the new_u:
            // Take the vector u, and multiply all the elements in it by (1-step_size)
            double[] new_u = multiply((1 - step_size), u);


            // Increment the jth element of new_u by step_size
            new_u[j] = new_u[j] + step_size;


            // Calculate error by taking finding the square root of the SSD
            // between new_u and u
            err = Math.sqrt(ssd(new_u, u));


            // Increment count and replace u
            count = count + 1;
            u = new_u;
        }


        // Compute center point
        //c = P * u
        double[][] c = multiply(P, u);
        center = transpose(c)[0];



        // Put the elements of the vector u into the diagonal of a matrix
        // U with the rest of the elements as 0
        double[][] U = diag(u);




        // Compute the A-matrix
        //A = (1/d) * inv(P * U * P' - (P * u)*(P*u)' );
        double[][] pup = multiply(multiply(P, U) , transpose(P));
        double[][] pupu = multiply((multiply(P, u)), transpose(multiply(P, u)));
        double[][] pup_pupu = minus(pup, pupu);
        A = multiply((1/d), inv(pup_pupu));



        // Compute Eigen vectors and values
        //A=inv(A);
        //[Ve,De]=eig(A);
        Eigen eig = new Eigen(inv(A));
        double[][] Ve = eig.getV(); //eigenvalues
        double[][] De = eig.getD(); //right eigenvectors
        reorderEigenVectors(De);
        reorderEigenValues(Ve);


        //v=sqrt(diag(De));
        double[] v = sqrt(diag(De));


        //[l1,Ie] = max(v);
        l1 = max(v);
        int Ie = find_maximum_value_location(v, l1); //off by one from MatLab but I think it's ok here



        //veig=Ve(:,Ie);
        double[] veig = new double[Ve.length];
        for (int i=0; i<veig.length; i++){
            veig[i] = Ve[Ie][i];
        }


        //thu=atan2(veig(2),veig(1));
        thu = Math.atan2(veig[1], veig[0]);


        //l2=v(setdiff([1 2],Ie));
        l2 = v[setdiff(new int[]{0,1}, Ie)];
    }


  //**************************************************************************
  //** getCenter
  //**************************************************************************
  /** Returns the center point of the ellipse
   */
    public double[] getCenter(){
        double[] pt = new double[2];
        pt[0] = center[0];
        pt[1] = center[1];
        return pt;
    }


  //**************************************************************************
  //** getMatrix
  //**************************************************************************
  /** Returns a matrix containing all the information regarding the shape of
   *  the ellipsoid. To get the radii and orientation of the ellipsoid take
   *  the Singular Value Decomposition of the matrix.
   */
    public double[][] getMatrix(){
        return A;
    }


  //**************************************************************************
  //** getBoundingCoordinates
  //**************************************************************************
  /** Returns a list of coordinates that can be used to render the ellipse.
   *  @param numPoints The number of points used to represent the ellipse.
   *  The higher the number the more dense the ellipse outline, the more
   *  accurate the shape.
   */
    public double[][] getBoundingCoordinates(int numPoints){

        //tq=linspace(-pi,pi,50);
        double[] tq = linspace(-Math.PI, Math.PI, numPoints);


        //U=[cos(thu) -sin(thu);sin(thu) cos(thu)]*[l1*cos(tq);l2*sin(tq)];
        double[][] U = multiply(
            new double[][]{
                createVector(Math.cos(thu), -Math.sin(thu)),
                createVector(Math.sin(thu), Math.cos(thu))
            },
            new double[][]{
                multiply(l1, cos(tq)),
                multiply(l2, sin(tq))
            }
        );
        //System.out.println(toString(transpose(U)));



        double[][] coords = transpose(U);
        for (int i=0; i<coords.length; i++){
            double x = coords[i][0] + center[0];
            double y = coords[i][1] + center[1];

            coords[i][0] = x;
            coords[i][1] = y;
        }

        return coords;
    }


  //**************************************************************************
  //** reorderEigenVectors
  //**************************************************************************
  /** Eigen values generated from Apache Common Math and JAMA are different
   *  than MatLab. The vectors are in the reverse order than expected. This
   *  function will update the array to what we expect to see in MatLab.
   */
    private void reorderEigenVectors(double[][] De){
        rotateMatrix(De);
        rotateMatrix(De);
    }


  //**************************************************************************
  //** reorderEigenValues
  //**************************************************************************
  /** Eigen values generated from Apache Common Math and JAMA are different
   *  than MatLab. The vectors are in reverse order than expected and with an
   *  opposite sign. This function will update the array to what we expect to
   *  see in MatLab.
   */
    private void reorderEigenValues(double[][] Ve){
        rotateMatrix(Ve);
        for (int i=0; i<Ve.length; i++){
            for (int j=0; j<Ve[i].length; j++){
                Ve[i][j] = -Ve[i][j];
            }
        }
    }


  //**************************************************************************
  //** linspace
  //**************************************************************************
    private double[] linspace(double min, double max, int points) {
        double[] d = new double[points];
        for (int i = 0; i < points; i++){
            d[i] = min + i * (max - min) / (points - 1);
        }
        return d;
    }



  //**************************************************************************
  //** ssd
  //**************************************************************************
  /** Returns the sum-of-square-differences between tow arrays. Takes two
   *  vectors of the same size, creates a new vector by finding the difference
   *  between corresponding elements, squaring each difference and adding them
   *  all together. So if the vectors were: a = [1 2 3] and b = [5 4 6], then:
   *  SSD = (1-5)^2 + (2-4)^2 + (3-6)^2;
   */
    private double ssd(double[] a, double[] b){
        double ssd = 0;
        for (int i=0; i<a.length; i++){
            ssd += Math.pow(a[i]-b[i], 2);
        }
        return ssd;
    }


  //**************************************************************************
  //** ones
  //**************************************************************************
  /** Creates an array of all ones. For example, ones(2,3) returns a 2-by-3
   *  array of ones.
    <pre>
        1 1 1
        1 1 1
    </pre>
   * Reference: https://www.mathworks.com/help/matlab/ref/ones.html
   */
    private double[][] ones(int rows, int cols){
        double[][] arr = new double[rows][];
        for (int i=0; i<arr.length; i++){
            double[] row = new double[cols];
            for (int j=0; j<row.length; j++){
                row[j] = 1;
            }
            arr[i] = row;
        }
        return arr;
    }


  //**************************************************************************
  //** merge
  //**************************************************************************
  /** Used to combine two arrays into one
   */
    private double[][] merge(double[][] m1, double[][] m2) {
        int x = 0;
        double[][] out = new double[m1.length + m2.length][];
        for (int i=0; i<m1.length; i++){
            out[x] = m1[i];
            x++;
        }
        for (int i=0; i<m2.length; i++){
            out[x] = m2[i];
            x++;
        }
        return out;
    }


  //**************************************************************************
  //** multiply
  //**************************************************************************
  /** Used to multiply all the values in the vector (arr) by n. This is called
   *  scalar multiplication.
   */
    private double[] multiply(double n, double[] arr){
        double[] out = new double[arr.length];
        for (int i=0; i<arr.length; i++){
            out[i] = arr[i]*n;
        }
        return out;
    }



  //**************************************************************************
  //** multiply
  //**************************************************************************
  /** Used to multiply all the values in the matrix (arr) by n
   */
    private double[][] multiply(double n, double[][] arr){
        double[][] out = new double[arr.length][];
        for (int i=0; i<arr.length; i++){
            double[] row = arr[i];
            double[] r = new double[row.length];
            for (int j=0; j<row.length; j++){
                r[j] = row[j]*n;
            }
            out[i] = r;
        }
        return out;
    }


  //**************************************************************************
  //** multiply
  //**************************************************************************
  /** Multiply a matrix with a vector by converting the vector to a matrix
   */
    private double[][] multiply(double[][] P, double[] u){
        double[][] m2 = new double[u.length][];
        for (int i=0; i<m2.length; i++){
            double[] row = new double[1];
            row[0] = u[i];
            m2[i] = row;
        }
        return multiply(P, m2);
    }


  //**************************************************************************
  //** multiply
  //**************************************************************************
  /** Used to multiply two matrices. Credit:
   *  https://stackoverflow.com/a/23817780
   */
    private double[][] multiply(double[][] m1, double[][] m2) {
        int m1ColLength = m1[0].length; // m1 columns length
        int m2RowLength = m2.length;    // m2 rows length
        if(m1ColLength != m2RowLength) return null; // matrix multiplication is not possible
        int mRRowLength = m1.length;    // m result rows length
        int mRColLength = m2[0].length; // m result columns length
        double[][] mResult = new double[mRRowLength][mRColLength];
        for(int i = 0; i < mRRowLength; i++) {         // rows from m1
            for(int j = 0; j < mRColLength; j++) {     // columns from m2
                for(int k = 0; k < m1ColLength; k++) { // columns from m1
                    mResult[i][j] += m1[i][k] * m2[k][j];
                }
            }
        }
        return mResult;
    }


  //**************************************************************************
  //** diag
  //**************************************************************************
  /** Returns a matrix for a given vector. The values in the vector will
   *  appear diagonally in the output.
   *  Reference: https://www.mathworks.com/help/matlab/ref/diag.html
   */
    private double[][] diag(double[] arr){
        double[][] out = new double[arr.length][];
        for (int i=0; i<arr.length; i++){
            double[] row = new double[arr.length];
            for (int j=0; j<row.length; j++){
                if (j==i) row[j] = arr[i];
                else row[j] = 0;
            }
            out[i] = row;
        }
        return out;
    }


  //**************************************************************************
  //** diag
  //**************************************************************************
  /** Returns a vector representing values that appear diagonally in the given
   *  matrix.
   *  Reference: https://www.mathworks.com/help/matlab/ref/diag.html
   */
    private double[] diag(double[][] arr){
        double[] out = new double[arr.length];
        for (int i=0; i<arr.length; i++){
            out[i] = arr[i][i];
        }
        return out;
    }


  //**************************************************************************
  //** transpose
  //**************************************************************************
  /** Interchanges the row and column index for each element
   *  Reference: https://www.mathworks.com/help/matlab/ref/transpose.html
   */
    private double[][] transpose(double[][] arr){
        int rows = arr.length;
        int cols = arr[0].length;

        double[][] out = new double[cols][rows];
        for (int x = 0; x < cols; x++) {
            for (int y = 0; y < rows; y++) {
                out[x][y] = arr[y][x];
            }
        }

        return out;
    }


  //**************************************************************************
  //** inv
  //**************************************************************************
  /** Returns the inverse of a matrix. Relies on 2 different implementations.
   *  The first implementation is more accurate (passes inverse check) but
   *  has the potential to fail. If so, falls back to second method that
   *  relies on partial-pivoting Gaussian elimination.
   *  Reference: https://www.mathworks.com/help/matlab/ref/inv.html
   */
    private double[][] inv(double[][] matrix){
        try{
            return inv1(matrix);
        }
        catch(Exception e){
            try{
                return inv2(matrix);
            }
            catch(Exception ex){
                throw new RuntimeException(ex);
            }
        }
    }


  //**************************************************************************
  //** inv1
  //**************************************************************************
  /** Returns the inverse of a matrix. This implementation passes inverse
   *  check so I think it's valid but it has a tendency to fail. For example,
   *  the following matrix fails with a ArrayIndexOutOfBoundsException in the
   *  determinant method.
    <pre>
        1171.18     658.33
         658.33    1039.55
    </pre>
   *  Credit: https://github.com/rchen8/Algorithms/blob/master/Matrix.java
   */
    private double[][] inv1(double[][] matrix){
        double[][] inverse = new double[matrix.length][matrix.length];

        // minors and cofactors
        for (int i = 0; i < matrix.length; i++)
            for (int j = 0; j < matrix[i].length; j++)
                inverse[i][j] = Math.pow(-1, i + j)
                        * determinant(minor(matrix, i, j));

        // adjugate and determinant
        double det = 1.0 / determinant(matrix);
        for (int i = 0; i < inverse.length; i++) {
            for (int j = 0; j <= i; j++) {
                double temp = inverse[i][j];
                inverse[i][j] = inverse[j][i] * det;
                inverse[j][i] = temp * det;
            }
        }

        return inverse;
    }
    private static double determinant(double[][] matrix) {
        if (matrix.length != matrix[0].length)
            throw new IllegalStateException("invalid dimensions");

        if (matrix.length == 2)
            return matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];

        double det = 0;
        for (int i = 0; i < matrix[0].length; i++)
            det += Math.pow(-1, i) * matrix[0][i]
                    * determinant(minor(matrix, 0, i));
        return det;
    }
    private static double[][] minor(double[][] matrix, int row, int column) {
        double[][] minor = new double[matrix.length - 1][matrix.length - 1];

        for (int i = 0; i < matrix.length; i++)
            for (int j = 0; i != row && j < matrix[i].length; j++)
                if (j != column)
                    minor[i < row ? i : i - 1][j < column ? j : j - 1] = matrix[i][j];
        return minor;
    }


  //**************************************************************************
  //** inv2
  //**************************************************************************
  /** Returns the inverse of a matrix. This implementation successfully
   *  executes but does not pass the inverse check.
   *  Credit: https://www.sanfoundry.com/java-program-find-inverse-matrix/
   */
    public static double[][] inv2(double a[][]){

        int n = a.length;
        double x[][] = new double[n][n];
        double b[][] = new double[n][n];
        int index[] = new int[n];

        for (int i=0; i<n; ++i)
            b[i][i] = 1;


          //Transform the matrix into an upper triangle
            gaussian(a, index);


          //Update the matrix b[i][j] with the ratios stored
            for (int i=0; i<n-1; ++i){
                for (int j=i+1; j<n; ++j){
                    for (int k=0; k<n; ++k){
                        b[index[j]][k]
                             -= a[index[j]][i]*b[index[i]][k];
                    }
                }
            }


      //Perform backward substitutions
        for (int i=0; i<n; ++i){
            x[n-1][i] = b[index[n-1]][i]/a[index[n-1]][n-1];
            for (int j=n-2; j>=0; --j){
                x[j][i] = b[index[j]][i];
                for (int k=j+1; k<n; ++k){
                    x[j][i] -= a[index[j]][k]*x[k][i];
                }
                x[j][i] /= a[index[j]][j];
            }
        }

        return x;
    }



    // Method to carry out the partial-pivoting Gaussian
    // elimination.  Here index[] stores pivoting order.
    public static void gaussian(double a[][], int index[]) {

        int n = index.length;
        double c[] = new double[n];


        // Initialize the index
        for (int i=0; i<n; ++i)
            index[i] = i;



        // Find the rescaling factors, one from each row
        for (int i=0; i<n; ++i) {
            double c1 = 0;
            for (int j=0; j<n; ++j){
                double c0 = Math.abs(a[i][j]);
                if (c0 > c1) c1 = c0;
            }
            c[i] = c1;
        }



        // Search the pivoting element from each column
        int k = 0;
        for (int j=0; j<n-1; ++j){
            double pi1 = 0;
            for (int i=j; i<n; ++i){
                double pi0 = Math.abs(a[index[i]][j]);
                pi0 /= c[index[i]];
                if (pi0 > pi1) {
                    pi1 = pi0;
                    k = i;
                }
            }



            // Interchange rows according to the pivoting order
            int itmp = index[j];
            index[j] = index[k];
            index[k] = itmp;
            for (int i=j+1; i<n; ++i){
                double pj = a[index[i]][j]/a[index[j]][j];

                // Record pivoting ratios below the diagonal
                a[index[i]][j] = pj;


                // Modify other elements accordingly
                for (int l=j+1; l<n; ++l)
                    a[index[i]][l] -= pj*a[index[j]][l];
            }
        }
    }



  //**************************************************************************
  //** max
  //**************************************************************************
  /** Returns the max value in a vector
   */
    private double max(double[] arr){
        double max = arr[0];
        for (double d : arr){
            max = Math.max(d, max);
        }
        return max;
    }


  //**************************************************************************
  //** find_maximum_value_location
  //**************************************************************************
  /** Returns the index of the max value in a vector
   */
    private int find_maximum_value_location(double[] arr, double max){
        for (int i=0; i<arr.length; i++){
            if (arr[i]==max) return i;
        }
        return 0;
    }


  //**************************************************************************
  //** minus
  //**************************************************************************
  /** Used to subtract array B from array A and returns the result
   *  Reference: https://www.mathworks.com/help/matlab/ref/minus.html
   */
    private double[][] minus(double[][] a, double[][] b){
        double[][] out = new double[a.length][];
        for (int i=0; i<out.length; i++){
            double[] row = new double[a[i].length];
            for (int j=0; j<row.length; j++){
                row[j] = a[i][j]-b[i][j];
            }
            out[i] = row;
        }
        return out;
    }


  //**************************************************************************
  //** sqrt
  //**************************************************************************
  /** Returns the square root of each element in a vector
   *  Reference: https://www.mathworks.com/help/matlab/ref/sqrt.html
   */
    private double[] sqrt(double[] arr){
        double[] out = new double[arr.length];
        for (int i=0; i<out.length; i++){
            out[i] = Math.sqrt(arr[i]);
        }
        return out;
    }


    private double[] cos(double[] arr){
        double[] out = new double[arr.length];
        for (int i=0; i<out.length; i++){
            out[i] = Math.cos(arr[i]);
        }
        return out;
    }

    private double[] sin(double[] arr){
        double[] out = new double[arr.length];
        for (int i=0; i<out.length; i++){
            out[i] = Math.sin(arr[i]);
        }
        return out;
    }


  //**************************************************************************
  //** setdiff
  //**************************************************************************
  /** Partial implementation of setdiff
   */
    private int setdiff(int[] arr, int x){
        for (int i : arr){
            if (i!=x) return i;
        }
        return 0; //?
    }


  //**************************************************************************
  //** rotateMatrix
  //**************************************************************************
    private void rotateMatrix(double mat[][]) {
        int N = mat[0].length;

        // Consider all squares one by one
        for (int x = 0; x < N / 2; x++)
        {
            // Consider elements in group of 4 in
            // current square
            for (int y = x; y < N-x-1; y++)
            {
                // store current cell in temp variable
                double temp = mat[x][y];

                // move values from right to top
                mat[x][y] = mat[y][N-1-x];

                // move values from bottom to right
                mat[y][N-1-x] = mat[N-1-x][N-1-y];

                // move values from left to bottom
                mat[N-1-x][N-1-y] = mat[N-1-y][x];

                // assign temp to left
                mat[N-1-y][x] = temp;
            }
        }
    }


  //**************************************************************************
  //** createVector
  //**************************************************************************
  /** Used to generate a vector for testing purposes
   */
    private double[] createVector(double ...d){
        double[] arr = new double[d.length];
        for (int i=0; i<arr.length; i++){
            arr[i] = d[i];
        }
        return arr;
    }
}

Here an example output using 10 random points and a tolerance of 0.001. The ellipse is rendered using straight lines connecting points generated via the Ellipse.getBoundingCoordinates() method using 50 points.

这是一个使用 10 个随机点和 0.001 容差的示例输出。椭圆使用直线连接通过 Ellipse.getBoundingCoordinates() 方法生成的点使用 50 个点进行渲染。

Minimum Bounding Ellipse

Minimum Bounding Ellipse