jme3/math/Eigen3f.java

/*
 * 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;
    }
}