History

Use Jacobi eigenvalue algorithm to get eigenvalues/eigenvectors.

Andrzej Kapolka committed 12 years ago

Commit 802b14d1 parent c4be35f1 file · +28 −100Browse files

Showing the changes to interface/src/MetavoxelSystem.cpp (1 of 1 file). Show all

Unified · Split

interface/src/MetavoxelSystem.cpp+28−100
@@ -1878,114 +1878,42 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 1877 unmodified lines
18781878 glm::mat3 atrans = glm::transpose(a);
18791879 glm::mat3 ata = atrans * a;
18801880
1881− // compute the SVD of ata; first, find the eigenvalues
1882− // (see http://en.wikipedia.org/wiki/Eigenvalue_algorithm#3.C3.973_matrices)
1883− glm::vec3 eigenvalues;
1884− float p1 = ata[0][1] * ata[0][1] + ata[0][2] * ata[0][2] + ata[1][2] * ata[1][2];
1885− if (p1 < EPSILON) {
1886− eigenvalues = glm::vec3(ata[0][0], ata[1][1], ata[2][2]);
1887− if (eigenvalues[2] < eigenvalues[1]) {
1888− qSwap(eigenvalues[2], eigenvalues[1]);
1881+ // find the eigenvalues and eigenvectors of ata
1882+ // (see http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm)
1883+ glm::mat3 d = ata;
1884+ glm::quat combinedRotation;
1885+ const int MAX_ITERATIONS = 20;
1886+ for (int i = 0; i < MAX_ITERATIONS; i++) {
1887+ glm::vec3 offDiagonals = glm::abs(glm::vec3(d[1][0], d[2][0], d[2][1]));
1888+ int largestIndex = (offDiagonals[0] > offDiagonals[1]) ? (offDiagonals[0] > offDiagonals[2] ? 0 : 2) :
1889+ (offDiagonals[1] > offDiagonals[2] ? 1 : 2);
1890+ const float DESIRED_PRECISION = 0.00001f;
1891+ if (offDiagonals[largestIndex] < DESIRED_PRECISION) {
1892+ break;
18891893 }
1890− if (eigenvalues[1] < eigenvalues[0]) {
1891− qSwap(eigenvalues[1], eigenvalues[0]);
1894+ int largestJ = (largestIndex == 2) ? 1 : 0;
1895+ int largestI = (largestIndex == 0) ? 1 : 2;
1896+ float sjj = d[largestJ][largestJ];
1897+ float sii = d[largestI][largestI];
1898+ float angle = (sii == sjj ? PI_OVER_TWO : glm::atan(2.0f * d[largestJ][largestI], sjj - sii)) / 2.0f;
1899+ glm::quat rotation = glm::angleAxis(angle, largestIndex == 0 ? glm::vec3(0.0f, 0.0f, -1.0f) :
1900+ (largestIndex == 1 ? glm::vec3(0.0f, 1.0f, 0.0f) : glm::vec3(-1.0f, 0.0f, 0.0f)));
1901+ if (rotation.w == 0.0f) {
1902+ break;
18921903 }
1893− if (eigenvalues[2] < eigenvalues[1]) {
1894− qSwap(eigenvalues[2], eigenvalues[1]);
1895− }
1896− } else {
1897− float q = (ata[0][0] + ata[1][1] + ata[2][2]) / 3.0f;
1898− float d1 = ata[0][0] - q, d2 = ata[1][1] - q, d3 = ata[2][2] - q;
1899− float p2 = d1 * d1 + d2 * d2 + d3 * d3 + 2.0f * p1;
1900− float p = glm::sqrt(p2 / 6.0f);
1901− glm::mat3 b = (ata - glm::mat3(q)) / p;
1902− float r = glm::determinant(b) / 2.0f;
1903− float phi;
1904− if (r <= -1.0f) {
1905− phi = PI / 3.0f;
1906− } else if (r >= 1.0f) {
1907− phi = 0.0f;
1908− } else {
1909− phi = glm::acos(r) / 3.0f;
1910− }
1911− eigenvalues[2] = q + 2.0f * p * glm::cos(phi);
1912− eigenvalues[0] = q + 2.0f * p * glm::cos(phi + (2.0f * PI / 3.0f));
1913− eigenvalues[1] = 3.0f * q - eigenvalues[0] - eigenvalues[2];
1904+ combinedRotation = glm::normalize(rotation * combinedRotation);
1905+ glm::mat3 matrix = glm::mat3_cast(combinedRotation);
1906+ d = matrix * ata * glm::transpose(matrix);
19141907 }
19151908
19161909 // form the singular matrix from the eigenvalues
1917− glm::mat3 d;
19181910 const float MIN_SINGULAR_THRESHOLD = 0.1f;
1919− d[0][0] = (eigenvalues[0] < MIN_SINGULAR_THRESHOLD) ? 0.0f : 1.0f / eigenvalues[0];
1920− d[1][1] = (eigenvalues[1] < MIN_SINGULAR_THRESHOLD) ? 0.0f : 1.0f / eigenvalues[1];
1921− d[2][2] = (eigenvalues[2] < MIN_SINGULAR_THRESHOLD) ? 0.0f : 1.0f / eigenvalues[2];
1922−
1923− glm::mat3 m[] = { ata - glm::mat3(eigenvalues[0]), ata - glm::mat3(eigenvalues[1]),
1924− ata - glm::mat3(eigenvalues[2]) };
1925−
1926− // form the orthogonal matrix from the eigenvectors
1927− // see http://www.geometrictools.com/Documentation/EigenSymmetric3x3.pdf
1928− bool same01 = glm::abs(eigenvalues[0] - eigenvalues[1]) < EPSILON;
1929− bool same12 = glm::abs(eigenvalues[1] - eigenvalues[2]) < EPSILON;
1930− glm::mat3 u;
1931− if (!(same01 && same12)) {
1932− if (same01 || same12) {
1933− int i = same01 ? 2 : 0;
1934− for (int j = 0; j < 3; j++) {
1935− glm::vec3 first = glm::vec3(m[i][0][j], m[i][1][j], m[i][2][j]);
1936− int j2 = (j + 1) % 3;
1937− glm::vec3 second = glm::vec3(m[i][0][j2], m[i][1][j2], m[i][2][j2]);
1938− glm::vec3 cross = glm::cross(first, second);
1939− float length = glm::length(cross);
1940− if (length > EPSILON) {
1941− u[0][i] = cross[0] / length;
1942− u[1][i] = cross[1] / length;
1943− u[2][i] = cross[2] / length;
1944− break;
1945− }
1946− }
1947− i = (i + 1) % 3;
1948− for (int j = 0; j < 3; j++) {
1949− glm::vec3 first = glm::vec3(m[i][0][j], m[i][1][j], m[i][2][j]);
1950− float length = glm::length(first);
1951− if (length > EPSILON) {
1952− glm::vec3 second = glm::cross(first, glm::vec3(1.0f, 0.0f, 0.0f));
1953− length = glm::length(second);
1954− if (length < EPSILON) {
1955− second = glm::cross(first, glm::vec3(0.0f, 1.0f, 0.0f));
1956− length = glm::length(second);
1957− }
1958− u[0][i] = second[0] / length;
1959− u[1][i] = second[1] / length;
1960− u[2][i] = second[2] / length;
1961− second = glm::normalize(glm::cross(second, first));
1962− i = (i + 1) % 3;
1963− u[0][i] = second[0];
1964− u[1][i] = second[1];
1965− u[2][i] = second[2];
1966− break;
1967− }
1968− }
1969− } else {
1970− for (int i = 0; i < 3; i++) {
1971− for (int j = 0; j < 3; j++) {
1972− glm::vec3 first = glm::vec3(m[i][0][j], m[i][1][j], m[i][2][j]);
1973− int j2 = (j + 1) % 3;
1974− glm::vec3 second = glm::vec3(m[i][0][j2], m[i][1][j2], m[i][2][j2]);
1975− glm::vec3 cross = glm::cross(first, second);
1976− float length = glm::length(cross);
1977− if (length > EPSILON) {
1978− u[0][i] = cross[0] / length;
1979− u[1][i] = cross[1] / length;
1980− u[2][i] = cross[2] / length;
1981− break;
1982− }
1983− }
1984− }
1985− }
1986− }
1911+ d[0][0] = (d[0][0] < MIN_SINGULAR_THRESHOLD) ? 0.0f : 1.0f / d[0][0];
1912+ d[1][1] = (d[1][1] < MIN_SINGULAR_THRESHOLD) ? 0.0f : 1.0f / d[1][1];
1913+ d[2][2] = (d[2][2] < MIN_SINGULAR_THRESHOLD) ? 0.0f : 1.0f / d[2][2];
19871914
19881915 // compute the pseudo-inverse, ataplus, and use to find the minimizing solution
1916+ glm::mat3 u = glm::mat3_cast(combinedRotation);
19891917 glm::mat3 ataplus = glm::transpose(u) * d * u;
19901918 glm::vec3 solution = (ataplus * atrans * b) + center;
19911919