FitEllipsoid { public static Object[] yuryPetrov(final double[][] points) { final int nPoints = points.length; if (nPoints < 9) { throw new IllegalArgumentException( "Too few points; need at least 9 to calculate a unique ellipsoid"); } final double[][] d = new double[nPoints][9]; for (int i = 0; i < nPoints; i++) { final double x = points[i][0]; final double y = points[i][1]; final double z = points[i][2]; d[i][0] = x * x; d[i][1] = y * y; d[i][2] = z * z; d[i][3] = 2 * x * y; d[i][4] = 2 * x * z; d[i][5] = 2 * y * z; d[i][6] = 2 * x; d[i][7] = 2 * y; d[i][8] = 2 * z; } final Matrix D = new Matrix(d); final Matrix ones = new Matrix(nPoints, 1, 1.0); final Matrix V = ((D.transpose().times(D)).inverse()).times(D.transpose() .times(ones)); final double[] v = V.getColumnPackedCopy(); final Object[] matrices = matrixFromEquation(v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8]); final EigenvalueDecomposition E = (EigenvalueDecomposition) matrices[3]; final Matrix eVal = E.getD(); final Matrix diagonal = MatrixUtils.diag(eVal); final int nEvals = diagonal.getRowDimension(); final double[] radii = new double[nEvals]; for (int i = 0; i < nEvals; i++) { final double eigenvalue = diagonal.get(i, 0); if (eigenvalue <= 0) { throw new IllegalArgumentException("Ellipsoid eigenvalues must be greater than 0 " + "(positive definite matrix)"); } radii[i] = Math.sqrt(1 / eigenvalue); } final double[] centre = (double[]) matrices[0]; final double[][] eigenVectors = (double[][]) matrices[2]; return new Object[] { centre, radii, eigenVectors, v, E }; } private  FitEllipsoid(); static double[][] testEllipsoid(final double a, final double b,
		final double c, final double angle, final double xCentre,
		final double yCentre, final double zCentre, final double noise,
		final int nPoints, final boolean random); static Object[] yuryPetrov(final double[][] points);  }
@Test public void testYuryPetrov() { final double[][] points = FitEllipsoid.testEllipsoid(a, b, c, r, x, y, z, 0.0, 1000, true); final double[][] noisyPoints = FitEllipsoid.testEllipsoid(a, b, c, r, x, y, z, 0.01, 1000, true); final double[] centre = { x, y, z }; final double[] radii = { c, b, a }; Object[] result = FitEllipsoid.yuryPetrov(points); assertArrayEquals(centre, (double[]) result[0], 1e-10); assertArrayEquals(radii, (double[]) result[1], 1e-10); result = FitEllipsoid.yuryPetrov(noisyPoints); assertArrayEquals(centre, (double[]) result[0], 1e-2); assertArrayEquals(radii, (double[]) result[1], 1e-2); }
@Test(expected = IllegalArgumentException.class) public void testYuryPetrovCylinder() { final int r = 50; final int nPoints = 100; final double[][] points = new double[(int)(2 * r * nPoints)][3]; final double radialIncrement = 2 * Math.PI / nPoints; int p = 0; for (double z = -r; z < r; z += 1) { for (double theta = -Math.PI; theta < Math.PI; theta += radialIncrement) { theta += radialIncrement; final double x = Math.cos(theta) * r; final double y = Math.sin(theta) * r; final double[] point = {x, y, z}; points[p] = point; p++; } } FitEllipsoid.yuryPetrov(points); }