/*
* Copyright (c) 2009-2012 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.FINE, "{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;
}
}