diff options
Diffstat (limited to 'engine/src/core/com/jme3/math/Eigen3f.java')
| -rw-r--r-- | engine/src/core/com/jme3/math/Eigen3f.java | 421 |
1 files changed, 421 insertions, 0 deletions
diff --git a/engine/src/core/com/jme3/math/Eigen3f.java b/engine/src/core/com/jme3/math/Eigen3f.java new file mode 100644 index 0000000..e356057 --- /dev/null +++ b/engine/src/core/com/jme3/math/Eigen3f.java @@ -0,0 +1,421 @@ +/*
+ * Copyright (c) 2009-2010 jMonkeyEngine
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are
+ * met:
+ *
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ *
+ * * Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * * Neither the name of 'jMonkeyEngine' nor the names of its contributors
+ * may be used to endorse or promote products derived from this software
+ * without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+package com.jme3.math;
+
+import java.util.logging.Level;
+import java.util.logging.Logger;
+
+public class Eigen3f implements java.io.Serializable {
+
+ static final long serialVersionUID = 1;
+
+ private static final Logger logger = Logger.getLogger(Eigen3f.class
+ .getName());
+
+ float[] eigenValues = new float[3];
+ Vector3f[] eigenVectors = new Vector3f[3];
+
+ static final double ONE_THIRD_DOUBLE = 1.0 / 3.0;
+ static final double ROOT_THREE_DOUBLE = Math.sqrt(3.0);
+
+
+ public Eigen3f() {
+
+ }
+
+ public Eigen3f(Matrix3f data) {
+ calculateEigen(data);
+ }
+
+ public void calculateEigen(Matrix3f data) {
+ // prep work...
+ eigenVectors[0] = new Vector3f();
+ eigenVectors[1] = new Vector3f();
+ eigenVectors[2] = new Vector3f();
+
+ Matrix3f scaledData = new Matrix3f(data);
+ float maxMagnitude = scaleMatrix(scaledData);
+
+ // Compute the eigenvalues using double-precision arithmetic.
+ double roots[] = new double[3];
+ computeRoots(scaledData, roots);
+ eigenValues[0] = (float) roots[0];
+ eigenValues[1] = (float) roots[1];
+ eigenValues[2] = (float) roots[2];
+
+ float[] maxValues = new float[3];
+ Vector3f[] maxRows = new Vector3f[3];
+ maxRows[0] = new Vector3f();
+ maxRows[1] = new Vector3f();
+ maxRows[2] = new Vector3f();
+
+ for (int i = 0; i < 3; i++) {
+ Matrix3f tempMatrix = new Matrix3f(scaledData);
+ tempMatrix.m00 -= eigenValues[i];
+ tempMatrix.m11 -= eigenValues[i];
+ tempMatrix.m22 -= eigenValues[i];
+ float[] val = new float[1];
+ val[0] = maxValues[i];
+ if (!positiveRank(tempMatrix, val, maxRows[i])) {
+ // Rank was 0 (or very close to 0), so just return.
+ // Rescale back to the original size.
+ if (maxMagnitude > 1f) {
+ for (int j = 0; j < 3; j++) {
+ eigenValues[j] *= maxMagnitude;
+ }
+ }
+
+ eigenVectors[0].set(Vector3f.UNIT_X);
+ eigenVectors[1].set(Vector3f.UNIT_Y);
+ eigenVectors[2].set(Vector3f.UNIT_Z);
+ return;
+ }
+ maxValues[i] = val[0];
+ }
+
+ float maxCompare = maxValues[0];
+ int i = 0;
+ if (maxValues[1] > maxCompare) {
+ maxCompare = maxValues[1];
+ i = 1;
+ }
+ if (maxValues[2] > maxCompare) {
+ i = 2;
+ }
+
+ // use the largest row to compute and order the eigen vectors.
+ switch (i) {
+ case 0:
+ maxRows[0].normalizeLocal();
+ computeVectors(scaledData, maxRows[0], 1, 2, 0);
+ break;
+ case 1:
+ maxRows[1].normalizeLocal();
+ computeVectors(scaledData, maxRows[1], 2, 0, 1);
+ break;
+ case 2:
+ maxRows[2].normalizeLocal();
+ computeVectors(scaledData, maxRows[2], 0, 1, 2);
+ break;
+ }
+
+ // Rescale the values back to the original size.
+ if (maxMagnitude > 1f) {
+ for (i = 0; i < 3; i++) {
+ eigenValues[i] *= maxMagnitude;
+ }
+ }
+ }
+
+ /**
+ * Scale the matrix so its entries are in [-1,1]. The scaling is applied
+ * only when at least one matrix entry has magnitude larger than 1.
+ *
+ * @return the max magnitude in this matrix
+ */
+ private float scaleMatrix(Matrix3f mat) {
+
+ float max = FastMath.abs(mat.m00);
+ float abs = FastMath.abs(mat.m01);
+
+ if (abs > max) {
+ max = abs;
+ }
+ abs = FastMath.abs(mat.m02);
+ if (abs > max) {
+ max = abs;
+ }
+ abs = FastMath.abs(mat.m11);
+ if (abs > max) {
+ max = abs;
+ }
+ abs = FastMath.abs(mat.m12);
+ if (abs > max) {
+ max = abs;
+ }
+ abs = FastMath.abs(mat.m22);
+ if (abs > max) {
+ max = abs;
+ }
+
+ if (max > 1f) {
+ float fInvMax = 1f / max;
+ mat.multLocal(fInvMax);
+ }
+
+ return max;
+ }
+
+ /**
+ * Compute the eigenvectors of the given Matrix, using the
+ * @param mat
+ * @param vect
+ * @param index1
+ * @param index2
+ * @param index3
+ */
+ private void computeVectors(Matrix3f mat, Vector3f vect, int index1,
+ int index2, int index3) {
+ Vector3f vectorU = new Vector3f(), vectorV = new Vector3f();
+ Vector3f.generateComplementBasis(vectorU, vectorV, vect);
+
+ Vector3f tempVect = mat.mult(vectorU);
+ float p00 = eigenValues[index3] - vectorU.dot(tempVect);
+ float p01 = vectorV.dot(tempVect);
+ float p11 = eigenValues[index3] - vectorV.dot(mat.mult(vectorV));
+ float invLength;
+ float max = FastMath.abs(p00);
+ int row = 0;
+ float fAbs = FastMath.abs(p01);
+ if (fAbs > max) {
+ max = fAbs;
+ }
+ fAbs = FastMath.abs(p11);
+ if (fAbs > max) {
+ max = fAbs;
+ row = 1;
+ }
+
+ if (max >= FastMath.ZERO_TOLERANCE) {
+ if (row == 0) {
+ invLength = FastMath.invSqrt(p00 * p00 + p01 * p01);
+ p00 *= invLength;
+ p01 *= invLength;
+ vectorU.mult(p01, eigenVectors[index3])
+ .addLocal(vectorV.mult(p00));
+ } else {
+ invLength = FastMath.invSqrt(p11 * p11 + p01 * p01);
+ p11 *= invLength;
+ p01 *= invLength;
+ vectorU.mult(p11, eigenVectors[index3])
+ .addLocal(vectorV.mult(p01));
+ }
+ } else {
+ if (row == 0) {
+ eigenVectors[index3] = vectorV;
+ } else {
+ eigenVectors[index3] = vectorU;
+ }
+ }
+
+ Vector3f vectorS = vect.cross(eigenVectors[index3]);
+ mat.mult(vect, tempVect);
+ p00 = eigenValues[index1] - vect.dot(tempVect);
+ p01 = vectorS.dot(tempVect);
+ p11 = eigenValues[index1] - vectorS.dot(mat.mult(vectorS));
+ max = FastMath.abs(p00);
+ row = 0;
+ fAbs = FastMath.abs(p01);
+ if (fAbs > max) {
+ max = fAbs;
+ }
+ fAbs = FastMath.abs(p11);
+ if (fAbs > max) {
+ max = fAbs;
+ row = 1;
+ }
+
+ if (max >= FastMath.ZERO_TOLERANCE) {
+ if (row == 0) {
+ invLength = FastMath.invSqrt(p00 * p00 + p01 * p01);
+ p00 *= invLength;
+ p01 *= invLength;
+ eigenVectors[index1] = vect.mult(p01).add(vectorS.mult(p00));
+ } else {
+ invLength = FastMath.invSqrt(p11 * p11 + p01 * p01);
+ p11 *= invLength;
+ p01 *= invLength;
+ eigenVectors[index1] = vect.mult(p11).add(vectorS.mult(p01));
+ }
+ } else {
+ if (row == 0) {
+ eigenVectors[index1].set(vectorS);
+ } else {
+ eigenVectors[index1].set(vect);
+ }
+ }
+
+ eigenVectors[index3].cross(eigenVectors[index1], eigenVectors[index2]);
+ }
+
+ /**
+ * Check the rank of the given Matrix to determine if it is positive. While
+ * doing so, store the max magnitude entry in the given float store and the
+ * max row of the matrix in the Vector store.
+ *
+ * @param matrix
+ * the Matrix3f to analyze.
+ * @param maxMagnitudeStore
+ * a float array in which to store (in the 0th position) the max
+ * magnitude entry of the matrix.
+ * @param maxRowStore
+ * a Vector3f to store the values of the row of the matrix
+ * containing the max magnitude entry.
+ * @return true if the given matrix has a non 0 rank.
+ */
+ private boolean positiveRank(Matrix3f matrix, float[] maxMagnitudeStore, Vector3f maxRowStore) {
+ // Locate the maximum-magnitude entry of the matrix.
+ maxMagnitudeStore[0] = -1f;
+ int iRow, iCol, iMaxRow = -1;
+ for (iRow = 0; iRow < 3; iRow++) {
+ for (iCol = iRow; iCol < 3; iCol++) {
+ float fAbs = FastMath.abs(matrix.get(iRow, iCol));
+ if (fAbs > maxMagnitudeStore[0]) {
+ maxMagnitudeStore[0] = fAbs;
+ iMaxRow = iRow;
+ }
+ }
+ }
+
+ // Return the row containing the maximum, to be used for eigenvector
+ // construction.
+ maxRowStore.set(matrix.getRow(iMaxRow));
+
+ return maxMagnitudeStore[0] >= FastMath.ZERO_TOLERANCE;
+ }
+
+ /**
+ * Generate the base eigen values of the given matrix using double precision
+ * math.
+ *
+ * @param mat
+ * the Matrix3f to analyze.
+ * @param rootsStore
+ * a double array to store the results in. Must be at least
+ * length 3.
+ */
+ private void computeRoots(Matrix3f mat, double[] rootsStore) {
+ // Convert the unique matrix entries to double precision.
+ double a = mat.m00, b = mat.m01, c = mat.m02,
+ d = mat.m11, e = mat.m12,
+ f = mat.m22;
+
+ // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0. The
+ // eigenvalues are the roots to this equation, all guaranteed to be
+ // real-valued, because the matrix is symmetric.
+ double char0 = a * d * f + 2.0 * b * c * e - a
+ * e * e - d * c * c - f * b * b;
+
+ double char1 = a * d - b * b + a * f - c * c
+ + d * f - e * e;
+
+ double char2 = a + d + f;
+
+ // Construct the parameters used in classifying the roots of the
+ // equation and in solving the equation for the roots in closed form.
+ double char2Div3 = char2 * ONE_THIRD_DOUBLE;
+ double abcDiv3 = (char1 - char2 * char2Div3) * ONE_THIRD_DOUBLE;
+ if (abcDiv3 > 0.0) {
+ abcDiv3 = 0.0;
+ }
+
+ double mbDiv2 = 0.5 * (char0 + char2Div3 * (2.0 * char2Div3 * char2Div3 - char1));
+
+ double q = mbDiv2 * mbDiv2 + abcDiv3 * abcDiv3 * abcDiv3;
+ if (q > 0.0) {
+ q = 0.0;
+ }
+
+ // Compute the eigenvalues by solving for the roots of the polynomial.
+ double magnitude = Math.sqrt(-abcDiv3);
+ double angle = Math.atan2(Math.sqrt(-q), mbDiv2) * ONE_THIRD_DOUBLE;
+ double cos = Math.cos(angle);
+ double sin = Math.sin(angle);
+ double root0 = char2Div3 + 2.0 * magnitude * cos;
+ double root1 = char2Div3 - magnitude
+ * (cos + ROOT_THREE_DOUBLE * sin);
+ double root2 = char2Div3 - magnitude
+ * (cos - ROOT_THREE_DOUBLE * sin);
+
+ // Sort in increasing order.
+ if (root1 >= root0) {
+ rootsStore[0] = root0;
+ rootsStore[1] = root1;
+ } else {
+ rootsStore[0] = root1;
+ rootsStore[1] = root0;
+ }
+
+ if (root2 >= rootsStore[1]) {
+ rootsStore[2] = root2;
+ } else {
+ rootsStore[2] = rootsStore[1];
+ if (root2 >= rootsStore[0]) {
+ rootsStore[1] = root2;
+ } else {
+ rootsStore[1] = rootsStore[0];
+ rootsStore[0] = root2;
+ }
+ }
+ }
+
+ public static void main(String[] args) {
+ Matrix3f mat = new Matrix3f(2, 1, 1, 1, 2, 1, 1, 1, 2);
+ Eigen3f eigenSystem = new Eigen3f(mat);
+
+ logger.info("eigenvalues = ");
+ for (int i = 0; i < 3; i++)
+ logger.log(Level.INFO, "{0} ", eigenSystem.getEigenValue(i));
+
+ logger.info("eigenvectors = ");
+ for (int i = 0; i < 3; i++) {
+ Vector3f vector = eigenSystem.getEigenVector(i);
+ logger.info(vector.toString());
+ mat.setColumn(i, vector);
+ }
+ logger.info(mat.toString());
+
+ // eigenvalues =
+ // 1.000000 1.000000 4.000000
+ // eigenvectors =
+ // 0.411953 0.704955 0.577350
+ // 0.404533 -0.709239 0.577350
+ // -0.816485 0.004284 0.577350
+ }
+
+ public float getEigenValue(int i) {
+ return eigenValues[i];
+ }
+
+ public Vector3f getEigenVector(int i) {
+ return eigenVectors[i];
+ }
+
+ public float[] getEigenValues() {
+ return eigenValues;
+ }
+
+ public Vector3f[] getEigenVectors() {
+ return eigenVectors;
+ }
+}
|