History

Ridiculously complicated minimization code.

Andrzej Kapolka committed 12 years ago

Commit 8ebb2a71 parent aa0b3a11 file · +134 −15Browse files

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

Unified · Split

interface/src/MetavoxelSystem.cpp+134−15
@@ -1657,7 +1657,7 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 1656 unmodified lines
16571657 crossing.color = colorX[0];
16581658 crossing.material = materialBase[0];
16591659 }
1660− crossing.point = glm::vec3(qAlpha(hermite), 0.0f, 0.0f);
1660+ crossing.point = glm::vec3(qAlpha(hermite) * EIGHT_BIT_MAXIMUM_RECIPROCAL, 0.0f, 0.0f);
16611661 }
16621662 if (middleY) {
16631663 if (alpha1 != alpha3) {
@@ -1671,7 +1671,7 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 7 unmodified lines
16711671 crossing.color = colorX[1];
16721672 crossing.material = materialBase[1];
16731673 }
1674− crossing.point = glm::vec3(EIGHT_BIT_MAXIMUM, qAlpha(hermite), 0.0f);
1674+ crossing.point = glm::vec3(1.0f, qAlpha(hermite) * EIGHT_BIT_MAXIMUM_RECIPROCAL, 0.0f);
16751675 }
16761676 if (alpha2 != alpha3) {
16771677 QRgb hermite = hermiteBase[hermiteStride];
@@ -1684,7 +1684,7 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 6 unmodified lines
16841684 crossing.color = colorX[size];
16851685 crossing.material = materialBase[size];
16861686 }
1687− crossing.point = glm::vec3(qAlpha(hermite), EIGHT_BIT_MAXIMUM, 0.0f);
1687+ crossing.point = glm::vec3(qAlpha(hermite) * EIGHT_BIT_MAXIMUM_RECIPROCAL, 1.0f, 0.0f);
16881688 }
16891689 if (middleZ) {
16901690 if (alpha3 != alpha7) {
@@ -1698,7 +1698,7 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 7 unmodified lines
16981698 crossing.color = colorX[offset3];
16991699 crossing.material = materialBase[offset3];
17001700 }
1701− crossing.point = glm::vec3(EIGHT_BIT_MAXIMUM, EIGHT_BIT_MAXIMUM, qAlpha(hermite));
1701+ crossing.point = glm::vec3(1.0f, 1.0f, qAlpha(hermite) * EIGHT_BIT_MAXIMUM_RECIPROCAL);
17021702 }
17031703 if (alpha5 != alpha7) {
17041704 QRgb hermite = hermiteBase[hermiteArea + VoxelHermiteData::EDGE_COUNT + 1];
@@ -1711,7 +1711,7 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 6 unmodified lines
17111711 crossing.color = colorX[offset5];
17121712 crossing.material = materialBase[offset5];
17131713 }
1714− crossing.point = glm::vec3(EIGHT_BIT_MAXIMUM, qAlpha(hermite), EIGHT_BIT_MAXIMUM);
1714+ crossing.point = glm::vec3(1.0f, qAlpha(hermite) * EIGHT_BIT_MAXIMUM_RECIPROCAL, 1.0f);
17151715 }
17161716 if (alpha6 != alpha7) {
17171717 QRgb hermite = hermiteBase[hermiteArea + hermiteStride];
@@ -1724,7 +1724,7 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 6 unmodified lines
17241724 crossing.color = colorX[offset6];
17251725 crossing.material = materialBase[offset6];
17261726 }
1727− crossing.point = glm::vec3(qAlpha(hermite), EIGHT_BIT_MAXIMUM, EIGHT_BIT_MAXIMUM);
1727+ crossing.point = glm::vec3(qAlpha(hermite) * EIGHT_BIT_MAXIMUM_RECIPROCAL, 1.0f, 1.0f);
17281728 }
17291729 }
17301730 }
@@ -1740,7 +1740,7 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 9 unmodified lines
17401740 crossing.color = colorX[1];
17411741 crossing.material = materialBase[1];
17421742 }
1743− crossing.point = glm::vec3(EIGHT_BIT_MAXIMUM, 0.0f, qAlpha(hermite));
1743+ crossing.point = glm::vec3(1.0f, 0.0f, qAlpha(hermite) * EIGHT_BIT_MAXIMUM_RECIPROCAL);
17441744 }
17451745 if (alpha4 != alpha5) {
17461746 QRgb hermite = hermiteBase[hermiteArea];
@@ -1753,7 +1753,7 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 6 unmodified lines
17531753 crossing.color = colorX[area];
17541754 crossing.material = materialBase[area];
17551755 }
1756− crossing.point = glm::vec3(qAlpha(hermite), 0.0f, EIGHT_BIT_MAXIMUM);
1756+ crossing.point = glm::vec3(qAlpha(hermite) * EIGHT_BIT_MAXIMUM_RECIPROCAL, 0.0f, 1.0f);
17571757 }
17581758 }
17591759 }
@@ -1769,7 +1769,7 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 9 unmodified lines
17691769 crossing.color = colorX[0];
17701770 crossing.material = materialBase[0];
17711771 }
1772− crossing.point = glm::vec3(0.0f, qAlpha(hermite), 0.0f);
1772+ crossing.point = glm::vec3(0.0f, qAlpha(hermite) * EIGHT_BIT_MAXIMUM_RECIPROCAL, 0.0f);
17731773 }
17741774 if (middleZ) {
17751775 if (alpha2 != alpha6) {
@@ -1783,7 +1783,7 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 7 unmodified lines
17831783 crossing.color = colorX[size];
17841784 crossing.material = materialBase[size];
17851785 }
1786− crossing.point = glm::vec3(0.0f, EIGHT_BIT_MAXIMUM, qAlpha(hermite));
1786+ crossing.point = glm::vec3(0.0f, 1.0f, qAlpha(hermite) * EIGHT_BIT_MAXIMUM_RECIPROCAL);
17871787 }
17881788 if (alpha4 != alpha6) {
17891789 QRgb hermite = hermiteBase[hermiteArea + 1];
@@ -1796,7 +1796,7 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 6 unmodified lines
17961796 crossing.color = colorX[area];
17971797 crossing.material = materialBase[area];
17981798 }
1799− crossing.point = glm::vec3(0.0f, qAlpha(hermite), EIGHT_BIT_MAXIMUM);
1799+ crossing.point = glm::vec3(0.0f, qAlpha(hermite) * EIGHT_BIT_MAXIMUM_RECIPROCAL, 1.0f);
18001800 }
18011801 }
18021802 }
@@ -1811,7 +1811,7 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 8 unmodified lines
18111811 crossing.color = colorX[0];
18121812 crossing.material = materialBase[0];
18131813 }
1814− crossing.point = glm::vec3(0.0f, 0.0f, qAlpha(hermite));
1814+ crossing.point = glm::vec3(0.0f, 0.0f, qAlpha(hermite) * EIGHT_BIT_MAXIMUM_RECIPROCAL);
18151815 }
18161816 // at present, we simply average the properties of each crossing as opposed to finding the vertex that
18171817 // minimizes the quadratic error function as described in the reference paper
@@ -1857,7 +1857,7 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 39 unmodified lines
18571857 glm::vec4 bottom;
18581858 for (int i = 0; i < crossingCount; i++) {
18591859 const EdgeCrossing& crossing = crossings[i];
1860− bottom = glm::vec4(crossing.normal, glm::dot(crossing.normal, crossing.point));
1860+ bottom = glm::vec4(crossing.normal, glm::dot(crossing.normal, crossing.point - center));
18611861
18621862 for (int j = 0; j < 4; j++) {
18631863 float angle = glm::atan(-bottom[j], r[j][j]);
@@ -1872,11 +1872,130 @@ int VoxelAugmentVisitor::visit(MetavoxelInfo& info) { · 8 unmodified lines
18721872 }
18731873 }
18741874
1875+ // extract the submatrices, form ata
1876+ glm::mat3 a(r);
1877+ glm::vec3 b(r[3]);
1878+ glm::mat3 atrans = glm::transpose(a);
1879+ glm::mat3 ata = atrans * a;
1880+
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]);
1889+ }
1890+ if (eigenvalues[1] < eigenvalues[0]) {
1891+ qSwap(eigenvalues[1], eigenvalues[0]);
1892+ }
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];
1914+ }
1915+
1916+ // form the singular matrix from the eigenvalues
1917+ glm::mat3 d;
1918+ 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+ }
1987+
1988+ // compute the pseudo-inverse, ataplus
1989+ glm::mat3 ataplus = glm::transpose(u) * d * u;
1990+
1991+ glm::vec3 solution = (ataplus * atrans * b) + center;
1992+
1993+ center = glm::clamp(solution, 0.0f, 1.0f);
1994+
18751995 if (totalWeight > 0.0f) {
18761996 materialWeights *= (EIGHT_BIT_MAXIMUM / totalWeight);
18771997 }
1878− VoxelPoint point = { info.minimum + (glm::vec3(clampedX, clampedY, clampedZ) +
1879− center * EIGHT_BIT_MAXIMUM_RECIPROCAL) * scale,
1998+ VoxelPoint point = { info.minimum + (glm::vec3(clampedX, clampedY, clampedZ) + center) * scale,
18801999 { (quint8)(red / crossingCount), (quint8)(green / crossingCount), (quint8)(blue / crossingCount) },
18812000 { (char)(normal.x * 127.0f), (char)(normal.y * 127.0f), (char)(normal.z * 127.0f) },
18822001 { materials[0], materials[1], materials[2], materials[3] },