diff --git a/src/modules/voxel/CMakeLists.txt b/src/modules/voxel/CMakeLists.txt index 3fb41a9704..5396ab8b2d 100644 --- a/src/modules/voxel/CMakeLists.txt +++ b/src/modules/voxel/CMakeLists.txt @@ -3,10 +3,12 @@ set(SRCS private/CubicSurfaceExtractor.h private/CubicSurfaceExtractor.cpp private/MarchingCubesSurfaceExtractor.h private/MarchingCubesSurfaceExtractor.cpp private/MarchingCubesTables.h + private/DualContouringSurfaceExtractor.h private/DualContouringSurfaceExtractor.cpp SurfaceExtractor.h SurfaceExtractor.cpp ChunkMesh.h Face.h Face.cpp + QEF.h QEF.cpp MaterialColor.h MaterialColor.cpp Mesh.h Mesh.cpp RawVolume.h RawVolume.cpp diff --git a/src/modules/voxel/QEF.cpp b/src/modules/voxel/QEF.cpp new file mode 100644 index 0000000000..6bcd46700e --- /dev/null +++ b/src/modules/voxel/QEF.cpp @@ -0,0 +1,801 @@ +//---------------------------------------------------------------------------- +// ThreeD Quadric Error Function +//---------------------------------------------------------------------------- + +#include "QEF.h" +#include "core/StandardLib.h" +#include +#include +#include + +namespace voxel { + +//---------------------------------------------------------------------------- + +#define MAXROWS 12 +#define EPSILON 1e-5 + +//---------------------------------------------------------------------------- + +glm::vec3 evaluateQEF(double mat[][3], double *vec, int rows) { + // perform singular value decomposition on matrix mat + // into u, v and d. + // u is a matrix of rows x 3 (same as mat); + // v is a square matrix 3 x 3 (for 3 columns in mat); + // d is vector of 3 values representing the diagonal + // matrix 3 x 3 (for 3 colums in mat). + double u[MAXROWS][3], v[3][3], d[3]; + computeSVD(mat, u, v, d, rows); + + // solve linear system given by mat and vec using the + // singular value decomposition of mat into u, v and d. + if (d[2] < 0.1) + d[2] = 0.0; + if (d[1] < 0.1) + d[1] = 0.0; + if (d[0] < 0.1) + d[0] = 0.0; + + double x[3]; + solveSVD(u, v, d, vec, x, rows); + + return {static_cast(x[0]), static_cast(x[1]), static_cast(x[2])}; +} + +//---------------------------------------------------------------------------- + +void computeSVD(double mat[][3], // matrix (rows x 3) + double u[][3], // matrix (rows x 3) + double v[3][3], // matrix (3x3) + double d[3], // vector (1x3) + int rows) { + core_memcpy(u, mat, rows * 3 * sizeof(double)); + + double *tau_u = d; + double tau_v[2]; + + factorize(u, tau_u, tau_v, rows); + + unpack(u, v, tau_u, tau_v, rows); + + diagonalize(u, v, tau_u, tau_v, rows); + + singularize(u, v, tau_u, rows); +} + +//---------------------------------------------------------------------------- + +void factorize(double mat[][3], // matrix (rows x 3) + double tau_u[3], // vector, (1x3) + double tau_v[2], // vector, (1x2) + int rows) { + int y; + + // bidiagonal factorization of (rows x 3) matrix into :- + // tau_u, a vector of 1x3 (for 3 columns in the matrix) + // tau_v, a vector of 1x2 (one less column than the matrix) + + for (int i = 0; i < 3; ++i) { + + // set up a vector to reference into the matrix + // from mat(i,i) to mat(m,i), that is, from the + // i'th column of the i'th row and down all the way + // through that column + double *ptrs[MAXROWS]; + int num_ptrs = rows - i; + for (int q = 0; q < num_ptrs; ++q) + ptrs[q] = &mat[q + i][i]; + + // perform householder transformation on this vector + double tau = factorize_hh(ptrs, num_ptrs); + tau_u[i] = tau; + + // all computations below this point are performed + // only for the first two columns: i=0 or i=1 + if (i + 1 < 3) { + + // perform householder transformation on the matrix + // mat(i,i+1) to mat(m,n), that is, on the sub-matrix + // that begins in the (i+1)'th column of the i'th + // row and extends to the end of the matrix at (m,n) + if (tau != 0.0) { + for (int x = i + 1; x < 3; ++x) { + double wx = mat[i][x]; + for (y = i + 1; y < rows; ++y) + wx += mat[y][x] * (*ptrs[y - i]); + double tau_wx = tau * wx; + mat[i][x] -= tau_wx; + for (y = i + 1; y < rows; ++y) + mat[y][x] -= tau_wx * (*ptrs[y - i]); + } + } + + // perform householder transformation on i'th row + // (remember at this point, i is either 0 or 1) + + // set up a vector to reference into the matrix + // from mat(i,i+1) to mat(i,n), that is, from the + // (i+1)'th column of the i'th row and all the way + // through to the end of that row + ptrs[0] = &mat[i][i + 1]; + if (i == 0) { + ptrs[1] = &mat[i][i + 2]; + num_ptrs = 2; + } else // i == 1 + num_ptrs = 1; + + // perform householder transformation on this vector + tau = factorize_hh(ptrs, num_ptrs); + tau_v[i] = tau; + + // perform householder transformation on the sub-matrix + // mat(i+1,i+1) to mat(m,n), that is, on the sub-matrix + // that begins in the (i+1)'th column of the (i+1)'th + // row and extends to the end of the matrix at (m,n) + if (tau != 0.0) { + for (y = i + 1; y < rows; ++y) { + double wy = mat[y][i + 1]; + if (i == 0) + wy += mat[y][i + 2] * (*ptrs[1]); + double tau_wy = tau * wy; + mat[y][i + 1] -= tau_wy; + if (i == 0) + mat[y][i + 2] -= tau_wy * (*ptrs[1]); + } + } + + } // if (i + 1 < 3) + } +} + +//---------------------------------------------------------------------------- + +double factorize_hh(double *ptrs[], int n) { + double tau = 0.0; + + if (n > 1) { + double xnorm; + if (n == 2) + xnorm = glm::abs(*ptrs[1]); + else { + double scl = 0.0; + double ssq = 1.0; + for (int i = 1; i < n; ++i) { + double x = glm::abs(*ptrs[i]); + if (x != 0.0) { + if (scl < x) { + ssq = 1.0 + ssq * (scl / x) * (scl / x); + scl = x; + } else + ssq += (x / scl) * (x / scl); + } + } + xnorm = scl * sqrt(ssq); + } + + if (xnorm != 0.0) { + double alpha = *ptrs[0]; + double beta = sqrt(alpha * alpha + xnorm * xnorm); + if (alpha >= 0.0) + beta = -beta; + tau = (beta - alpha) / beta; + + double scl = 1.0 / (alpha - beta); + *ptrs[0] = beta; + for (int i = 1; i < n; ++i) + *ptrs[i] *= scl; + } + } + + return tau; +} + +//---------------------------------------------------------------------------- + +void unpack(double u[][3], // matrix (rows x 3) + double v[3][3], // matrix (3x3) + double tau_u[3], // vector, (1x3) + double tau_v[2], // vector, (1x2) + int rows) { + int i, y; + + // reset v to the identity matrix + v[0][0] = v[1][1] = v[2][2] = 1.0; + v[0][1] = v[0][2] = v[1][0] = v[1][2] = v[2][0] = v[2][1] = 0.0; + + for (i = 1; i >= 0; --i) { + double tau = tau_v[i]; + + // perform householder transformation on the sub-matrix + // v(i+1,i+1) to v(m,n), that is, on the sub-matrix of v + // that begins in the (i+1)'th column of the (i+1)'th row + // and extends to the end of the matrix at (m,n). the + // householder vector used to perform this is the vector + // from u(i,i+1) to u(i,n) + if (tau != 0.0) { + for (int x = i + 1; x < 3; ++x) { + double wx = v[i + 1][x]; + for (y = i + 1 + 1; y < 3; ++y) + wx += v[y][x] * u[i][y]; + double tau_wx = tau * wx; + v[i + 1][x] -= tau_wx; + for (y = i + 1 + 1; y < 3; ++y) + v[y][x] -= tau_wx * u[i][y]; + } + } + } + + // copy superdiagonal of u into tau_v + for (i = 0; i < 2; ++i) + tau_v[i] = u[i][i + 1]; + + // below, same idea for u: householder transformations + // and the superdiagonal copy + + for (i = 2; i >= 0; --i) { + // copy superdiagonal of u into tau_u + double tau = tau_u[i]; + tau_u[i] = u[i][i]; + + // perform householder transformation on the sub-matrix + // u(i,i) to u(m,n), that is, on the sub-matrix of u that + // begins in the i'th column of the i'th row and extends + // to the end of the matrix at (m,n). the householder + // vector used to perform this is the i'th column of u, + // that is, u(0,i) to u(m,i) + if (tau == 0.0) { + u[i][i] = 1.0; + if (i < 2) { + u[i][2] = 0.0; + if (i < 1) + u[i][1] = 0.0; + } + for (y = i + 1; y < rows; ++y) + u[y][i] = 0.0; + } else { + for (int x = i + 1; x < 3; ++x) { + double wx = 0.0; + for (y = i + 1; y < rows; ++y) + wx += u[y][x] * u[y][i]; + double tau_wx = tau * wx; + u[i][x] = -tau_wx; + for (y = i + 1; y < rows; ++y) + u[y][x] -= tau_wx * u[y][i]; + } + for (y = i + 1; y < rows; ++y) + u[y][i] = u[y][i] * -tau; + u[i][i] = 1.0 - tau; + } + } +} + +//---------------------------------------------------------------------------- + +void diagonalize(double u[][3], // matrix (rows x 3) + double v[3][3], // matrix (3x3) + double tau_u[3], // vector, (1x3) + double tau_v[2], // vector, (1x2) + int rows) { + int i, j; + + chop(tau_u, tau_v, 3); + + // progressively reduce the matrices into diagonal form + + int b = 3 - 1; + while (b > 0) { + if (tau_v[b - 1] == 0.0) + --b; + else { + int a = b - 1; + while (a > 0 && tau_v[a - 1] != 0.0) + --a; + int n = b - a + 1; + + double u1[MAXROWS][3]; + double v1[3][3]; + for (j = a; j <= b; ++j) { + for (i = 0; i < rows; ++i) + u1[i][j - a] = u[i][j]; + for (i = 0; i < 3; ++i) + v1[i][j - a] = v[i][j]; + } + + qrstep(u1, v1, &tau_u[a], &tau_v[a], rows, n); + + for (j = a; j <= b; ++j) { + for (i = 0; i < rows; ++i) + u[i][j] = u1[i][j - a]; + for (i = 0; i < 3; ++i) + v[i][j] = v1[i][j - a]; + } + + chop(&tau_u[a], &tau_v[a], n); + } + } +} + +//---------------------------------------------------------------------------- + +void chop(double *a, double *b, int n) { + double ai = a[0]; + for (int i = 0; i < n - 1; ++i) { + double bi = b[i]; + double ai1 = a[i + 1]; + if (glm::abs(bi) < EPSILON * (glm::abs(ai) + glm::abs(ai1))) + b[i] = 0.0; + ai = ai1; + } +} + +//---------------------------------------------------------------------------- + +void qrstep(double u[][3], // matrix (rows x cols) + double v[][3], // matrix (3 x cols) + double tau_u[], // vector (1 x cols) + double tau_v[], // vector (1 x cols - 1) + int rows, int cols) { + int i; + + if (cols == 2) { + qrstep_cols2(u, v, tau_u, tau_v, rows); + return; + } + + // handle zeros on the diagonal or at its end + for (i = 0; i < cols - 1; ++i) + if (tau_u[i] == 0.0) { + qrstep_middle(u, tau_u, tau_v, rows, cols, i); + return; + } + if (tau_u[cols - 1] == 0.0) { + qrstep_end(v, tau_u, tau_v, cols); + return; + } + + // perform qr reduction on the diagonal and off-diagonal + + double mu = qrstep_eigenvalue(tau_u, tau_v, cols); + double y = tau_u[0] * tau_u[0] - mu; + double z = tau_u[0] * tau_v[0]; + + double ak = 0.0; + double bk = 0.0; + double zk; + double ap = tau_u[0]; + double bp = tau_v[0]; + double aq = tau_u[1]; + //double bq = tau_v[1]; + + for (int k = 0; k < cols - 1; ++k) { + double c, s; + + // perform Givens rotation on V + + computeGivens(y, z, &c, &s); + + for (i = 0; i < 3; ++i) { + double vip = v[i][k]; + double viq = v[i][k + 1]; + v[i][k] = vip * c - viq * s; + v[i][k + 1] = vip * s + viq * c; + } + + // perform Givens rotation on B + + double bk1 = bk * c - z * s; + double ap1 = ap * c - bp * s; + double bp1 = ap * s + bp * c; + double zp1 = aq * -s; + double aq1 = aq * c; + + if (k > 0) + tau_v[k - 1] = bk1; + + ak = ap1; + bk = bp1; + zk = zp1; + ap = aq1; + + if (k < cols - 2) + bp = tau_v[k + 1]; + else + bp = 0.0; + + y = ak; + z = zk; + + // perform Givens rotation on U + + computeGivens(y, z, &c, &s); + + for (i = 0; i < rows; ++i) { + double uip = u[i][k]; + double uiq = u[i][k + 1]; + u[i][k] = uip * c - uiq * s; + u[i][k + 1] = uip * s + uiq * c; + } + + // perform Givens rotation on B + + double ak1 = ak * c - zk * s; + bk1 = bk * c - ap * s; + double zk1 = bp * -s; + + ap1 = bk * s + ap * c; + bp1 = bp * c; + + tau_u[k] = ak1; + + ak = ak1; + bk = bk1; + zk = zk1; + ap = ap1; + bp = bp1; + + if (k < cols - 2) + aq = tau_u[k + 2]; + else + aq = 0.0; + + y = bk; + z = zk; + } + + tau_v[cols - 2] = bk; + tau_u[cols - 1] = ap; +} + +//---------------------------------------------------------------------------- + +void qrstep_middle(double u[][3], // matrix (rows x cols) + double tau_u[], // vector (1 x cols) + double tau_v[], // vector (1 x cols - 1) + int rows, int cols, int col) { + double x = tau_v[col]; + double y = tau_u[col + 1]; + for (int j = col; j < cols - 1; ++j) { + double c, s; + + // perform Givens rotation on U + + computeGivens(y, -x, &c, &s); + for (int i = 0; i < rows; ++i) { + double uip = u[i][col]; + double uiq = u[i][j + 1]; + u[i][col] = uip * c - uiq * s; + u[i][j + 1] = uip * s + uiq * c; + } + + // perform transposed Givens rotation on B + + tau_u[j + 1] = x * s + y * c; + if (j == col) + tau_v[j] = x * c - y * s; + + if (j < cols - 2) { + double z = tau_v[j + 1]; + tau_v[j + 1] *= c; + x = z * -s; + y = tau_u[j + 2]; + } + } +} + +//---------------------------------------------------------------------------- + +void qrstep_end(double v[][3], // matrix (3 x 3) + double tau_u[], // vector (1 x 3) + double tau_v[], // vector (1 x 2) + int cols) { + double x = tau_u[1]; + double y = tau_v[1]; + + for (int k = 1; k >= 0; --k) { + double c, s; + + // perform Givens rotation on V + + computeGivens(x, y, &c, &s); + + for (int i = 0; i < 3; ++i) { + double vip = v[i][k]; + double viq = v[i][2]; + v[i][k] = vip * c - viq * s; + v[i][2] = vip * s + viq * c; + } + + // perform Givens rotation on B + + tau_u[k] = x * c - y * s; + if (k == 1) + tau_v[k] = x * s + y * c; + if (k > 0) { + double z = tau_v[k - 1]; + tau_v[k - 1] *= c; + + x = tau_u[k - 1]; + y = z * s; + } + } +} + +//---------------------------------------------------------------------------- + +double qrstep_eigenvalue(double tau_u[], // vector (1 x 3) + double tau_v[], // vector (1 x 2) + int cols) { + double ta = tau_u[1] * tau_u[1] + tau_v[0] * tau_v[0]; + double tb = tau_u[2] * tau_u[2] + tau_v[1] * tau_v[1]; + double tab = tau_u[1] * tau_v[1]; + double dt = (ta - tb) / 2.0; + double mu; + if (dt >= 0.0) + mu = tb - (tab * tab) / (dt + sqrt(dt * dt + tab * tab)); + else + mu = tb + (tab * tab) / (sqrt(dt * dt + tab * tab) - dt); + return mu; +} + +//---------------------------------------------------------------------------- + +void qrstep_cols2(double u[][3], // matrix (rows x 2) + double v[][3], // matrix (3 x 2) + double tau_u[], // vector (1 x 2) + double tau_v[], // vector (1 x 1) + int rows) { + int i; + double tmp; + + // eliminate off-diagonal element in [ 0 tau_v0 ] + // [ 0 tau_u1 ] + // to make [ tau_u[0] 0 ] + // [ 0 0 ] + + if (tau_u[0] == 0.0) { + double c, s; + + // perform transposed Givens rotation on B + // multiplied by X = [ 0 1 ] + // [ 1 0 ] + + computeGivens(tau_v[0], tau_u[1], &c, &s); + + tau_u[0] = tau_v[0] * c - tau_u[1] * s; + tau_v[0] = tau_v[0] * s + tau_u[1] * c; + tau_u[1] = 0.0; + + // perform Givens rotation on U + + for (i = 0; i < rows; ++i) { + double uip = u[i][0]; + double uiq = u[i][1]; + u[i][0] = uip * c - uiq * s; + u[i][1] = uip * s + uiq * c; + } + + // multiply V by X, effectively swapping first two columns + + for (i = 0; i < 3; ++i) { + tmp = v[i][0]; + v[i][0] = v[i][1]; + v[i][1] = tmp; + } + } + + // eliminate off-diagonal element in [ tau_u0 tau_v0 ] + // [ 0 0 ] + + else if (tau_u[1] == 0.0) { + double c, s; + + // perform Givens rotation on B + + computeGivens(tau_u[0], tau_v[0], &c, &s); + + tau_u[0] = tau_u[0] * c - tau_v[0] * s; + tau_v[0] = 0.0; + + // perform Givens rotation on V + + for (i = 0; i < 3; ++i) { + double vip = v[i][0]; + double viq = v[i][1]; + v[i][0] = vip * c - viq * s; + v[i][1] = vip * s + viq * c; + } + } + + // make colums orthogonal, + + else { + double c, s; + + // perform Schur rotation on B + + computeSchur(tau_u[0], tau_v[0], tau_u[1], &c, &s); + + double a11 = tau_u[0] * c - tau_v[0] * s; + double a21 = -tau_u[1] * s; + double a12 = tau_u[0] * s + tau_v[0] * c; + double a22 = tau_u[1] * c; + + // perform Schur rotation on V + + for (i = 0; i < 3; ++i) { + double vip = v[i][0]; + double viq = v[i][1]; + v[i][0] = vip * c - viq * s; + v[i][1] = vip * s + viq * c; + } + + // eliminate off diagonal elements + + if ((a11 * a11 + a21 * a21) < (a12 * a12 + a22 * a22)) { + + // multiply B by X + + tmp = a11; + a11 = a12; + a12 = tmp; + tmp = a21; + a21 = a22; + a22 = tmp; + + // multiply V by X, effectively swapping first + // two columns + + for (i = 0; i < 3; ++i) { + tmp = v[i][0]; + v[i][0] = v[i][1]; + v[i][1] = tmp; + } + } + + // perform transposed Givens rotation on B + + computeGivens(a11, a21, &c, &s); + + tau_u[0] = a11 * c - a21 * s; + tau_v[0] = a12 * c - a22 * s; + tau_u[1] = a12 * s + a22 * c; + + // perform Givens rotation on U + + for (i = 0; i < rows; ++i) { + double uip = u[i][0]; + double uiq = u[i][1]; + u[i][0] = uip * c - uiq * s; + u[i][1] = uip * s + uiq * c; + } + } +} + +//---------------------------------------------------------------------------- + +void computeGivens(double a, double b, double *c, double *s) { + if (b == 0.0) { + *c = 1.0; + *s = 0.0; + } else if (glm::abs(b) > glm::abs(a)) { + double t = -a / b; + double s1 = 1.0 / sqrt(1 + t * t); + *s = s1; + *c = s1 * t; + } else { + double t = -b / a; + double c1 = 1.0 / sqrt(1 + t * t); + *c = c1; + *s = c1 * t; + } +} + +//---------------------------------------------------------------------------- + +void computeSchur(double a1, double a2, double a3, double *c, double *s) { + double apq = a1 * a2 * 2.0; + + if (apq == 0.0) { + *c = 1.0; + *s = 0.0; + } else { + double t; + double tau = (a2 * a2 + (a3 + a1) * (a3 - a1)) / apq; + if (tau >= 0.0) + t = 1.0 / (tau + sqrt(1.0 + tau * tau)); + else + t = -1.0 / (sqrt(1.0 + tau * tau) - tau); + *c = 1.0 / sqrt(1.0 + t * t); + *s = t * (*c); + } +} + +//---------------------------------------------------------------------------- + +void singularize(double u[][3], // matrix (rows x 3) + double v[3][3], // matrix (3x3) + double d[3], // vector, (1x3) + int rows) { + int i, j, y; + + // make singularize values positive + + for (j = 0; j < 3; ++j) + if (d[j] < 0.0) { + for (i = 0; i < 3; ++i) + v[i][j] = -v[i][j]; + d[j] = -d[j]; + } + + // sort singular values in decreasing order + + for (i = 0; i < 3; ++i) { + double d_max = d[i]; + int i_max = i; + for (j = i + 1; j < 3; ++j) + if (d[j] > d_max) { + d_max = d[j]; + i_max = j; + } + + if (i_max != i) { + // swap eigenvalues + double tmp = d[i]; + d[i] = d[i_max]; + d[i_max] = tmp; + + // swap eigenvectors + for (y = 0; y < rows; ++y) { + tmp = u[y][i]; + u[y][i] = u[y][i_max]; + u[y][i_max] = tmp; + } + for (y = 0; y < 3; ++y) { + tmp = v[y][i]; + v[y][i] = v[y][i_max]; + v[y][i_max] = tmp; + } + } + } +} + +//---------------------------------------------------------------------------- + +void solveSVD(double u[][3], // matrix (rows x 3) + double v[3][3], // matrix (3x3) + double d[3], // vector (1x3) + double b[], // vector (1 x rows) + double x[3], // vector (1x3) + int rows) { + static double zeroes[3] = {0.0, 0.0, 0.0}; + + int i, j; + + // compute vector w = U^T * b + + double w[3]; + core_memcpy(w, zeroes, sizeof(w)); + for (i = 0; i < rows; ++i) { + if (b[i] != 0.0) + for (j = 0; j < 3; ++j) + w[j] += b[i] * u[i][j]; + } + + // introduce non-zero singular values in d into w + + for (i = 0; i < 3; ++i) { + if (d[i] != 0.0) + w[i] /= d[i]; + } + + // compute result vector x = V * w + + for (i = 0; i < 3; ++i) { + double tmp = 0.0; + for (j = 0; j < 3; ++j) + tmp += w[j] * v[i][j]; + x[i] = tmp; + } +} + +} // namespace voxel diff --git a/src/modules/voxel/QEF.h b/src/modules/voxel/QEF.h new file mode 100644 index 0000000000..eb2cda59e3 --- /dev/null +++ b/src/modules/voxel/QEF.h @@ -0,0 +1,99 @@ +//---------------------------------------------------------------------------- +// ThreeD Quadric Error Function +//---------------------------------------------------------------------------- +#pragma once + +#include + +namespace voxel { + +/** + * QEF, implementing the quadric error function + * E[x] = P - Ni . Pi + * + * Given at least three points Pi, each with its respective + * normal vector Ni, that describe at least two planes, + * the QEF evalulates to the point x. + */ +glm::vec3 evaluateQEF(double mat[][3], double *vec, int rows); // TODO make these STL containers + +// compute svd + +void computeSVD(double mat[][3], // matrix (rows x 3) + double u[][3], // matrix (rows x 3) + double v[3][3], // matrix (3x3) + double d[3], // vector (1x3) + int rows); + +// factorize + +void factorize(double mat[][3], // matrix (rows x 3) + double tau_u[3], // vector (1x3) + double tau_v[2], // vectors, (1x2) + int rows); + +double factorize_hh(double *ptrs[], int n); + +// unpack + +void unpack(double u[][3], // matrix (rows x 3) + double v[3][3], // matrix (3x3) + double tau_u[3], // vector, (1x3) + double tau_v[2], // vector, (1x2) + int rows); + +// diagonalize + +void diagonalize(double u[][3], // matrix (rows x 3) + double v[3][3], // matrix (3x3) + double tau_u[3], // vector, (1x3) + double tau_v[2], // vector, (1x2) + int rows); + +void chop(double *a, double *b, int n); + +void qrstep(double u[][3], // matrix (rows x cols) + double v[][3], // matrix (3 x cols) + double tau_u[], // vector (1 x cols) + double tau_v[], // vector (1 x cols - 1) + int rows, int cols); + +void qrstep_middle(double u[][3], // matrix (rows x cols) + double tau_u[], // vector (1 x cols) + double tau_v[], // vector (1 x cols - 1) + int rows, int cols, int col); + +void qrstep_end(double v[][3], // matrix (3 x 3) + double tau_u[], // vector (1 x 3) + double tau_v[], // vector (1 x 2) + int cols); + +double qrstep_eigenvalue(double tau_u[], // vector (1 x 3) + double tau_v[], // vector (1 x 2) + int cols); + +void qrstep_cols2(double u[][3], // matrix (rows x 2) + double v[][3], // matrix (3 x 2) + double tau_u[], // vector (1 x 2) + double tau_v[], // vector (1 x 1) + int rows); + +void computeGivens(double a, double b, double *c, double *s); + +void computeSchur(double a1, double a2, double a3, double *c, double *s); + +// singularize + +void singularize(double u[][3], // matrix (rows x 3) + double v[3][3], // matrix (3x3) + double d[3], // vector, (1x3) + int rows); + +// solve svd +void solveSVD(double u[][3], // matrix (rows x 3) + double v[3][3], // matrix (3x3) + double d[3], // vector (1x3) + double b[], // vector (1 x rows) + double x[3], // vector (1x3) + int rows); +} // namespace voxel diff --git a/src/modules/voxel/SurfaceExtractor.cpp b/src/modules/voxel/SurfaceExtractor.cpp index c49ff2e105..71bdddf688 100644 --- a/src/modules/voxel/SurfaceExtractor.cpp +++ b/src/modules/voxel/SurfaceExtractor.cpp @@ -6,6 +6,7 @@ #include "voxel/MaterialColor.h" #include "voxel/Region.h" #include "voxel/private/CubicSurfaceExtractor.h" +#include "voxel/private/DualContouringSurfaceExtractor.h" #include "voxel/private/MarchingCubesSurfaceExtractor.h" namespace voxel { @@ -23,11 +24,21 @@ SurfaceExtractionContext buildMarchingCubesContext(const RawVolume *volume, cons false, false, false); } +SurfaceExtractionContext buildDualContouringContext(const RawVolume *volume, const Region ®ion, ChunkMesh &mesh, + const palette::Palette &palette) { + return SurfaceExtractionContext(volume, palette, region, mesh, glm::ivec3(0), SurfaceExtractionType::DualContouring, + false, false, false); +} + void extractSurface(SurfaceExtractionContext &ctx) { if (ctx.type == SurfaceExtractionType::MarchingCubes) { voxel::Region extractRegion = ctx.region; extractRegion.shrink(-1); voxel::extractMarchingCubesMesh(ctx.volume, ctx.palette, extractRegion, &ctx.mesh); + } else if (ctx.type == SurfaceExtractionType::DualContouring) { + voxel::Region extractRegion = ctx.region; + extractRegion.shrink(-1); + voxel::extractDualContouringMesh(ctx.volume, ctx.palette, extractRegion, &ctx.mesh); } else { voxel::Region extractRegion = ctx.region; extractRegion.shiftUpperCorner(1, 1, 1); diff --git a/src/modules/voxel/SurfaceExtractor.h b/src/modules/voxel/SurfaceExtractor.h index d62788f308..4d12b79581 100644 --- a/src/modules/voxel/SurfaceExtractor.h +++ b/src/modules/voxel/SurfaceExtractor.h @@ -15,7 +15,7 @@ class RawVolume; class Region; struct ChunkMesh; -enum class SurfaceExtractionType { Cubic, MarchingCubes, Max }; +enum class SurfaceExtractionType { Cubic, MarchingCubes, DualContouring, Max }; struct SurfaceExtractionContext { SurfaceExtractionContext(const RawVolume *_volume, const palette::Palette &_palette, const Region &_region, @@ -40,7 +40,8 @@ SurfaceExtractionContext buildCubicContext(const RawVolume *volume, const Region bool reuseVertices = true, bool ambientOcclusion = true); SurfaceExtractionContext buildMarchingCubesContext(const RawVolume *volume, const Region ®ion, ChunkMesh &mesh, const palette::Palette &palette); - +SurfaceExtractionContext buildDualContouringContext(const RawVolume *volume, const Region ®ion, ChunkMesh &mesh, + const palette::Palette &palette); void extractSurface(SurfaceExtractionContext &ctx); } // namespace voxel diff --git a/src/modules/voxel/private/DualContouringSurfaceExtractor.cpp b/src/modules/voxel/private/DualContouringSurfaceExtractor.cpp new file mode 100644 index 0000000000..d8cc54fe76 --- /dev/null +++ b/src/modules/voxel/private/DualContouringSurfaceExtractor.cpp @@ -0,0 +1,376 @@ +/** + * @file + */ + +#include "core/Assert.h" +#include "core/Pair.h" +#include "palette/Palette.h" +#include "voxel/Mesh.h" +#include "voxel/QEF.h" +#include "voxel/RawVolume.h" +#include "voxel/Region.h" +#include "voxel/Voxel.h" +#include "voxel/VoxelVertex.h" +#include "voxel/ChunkMesh.h" +#ifndef GLM_ENABLE_EXPERIMENTAL +#define GLM_ENABLE_EXPERIMENTAL +#endif +#include +#include +#include +#include + +// BUG We will get duplication of edges if the surface is along region boundaries + +namespace voxel { +namespace { + +struct PositionNormal { + glm::vec3 position; + glm::vec3 normal; +}; + +struct EdgeData { + glm::vec3 normal{0.0f}; + /** fraction (0.0-1.0) along the edge in the positive direction that the intersection happens */ + float fraction = 0.0f; + bool intersects = false; +}; + +struct CellData { + EdgeData edges[3]; + IndexType vertexIndex; +}; + +using ThresholdType = float; + +static const float DualContouringMaxDensity = 255.0f; + +static inline float convertToDensity(const Voxel &voxel) { + return isAir(voxel.getMaterial()) ? 0.0f : DualContouringMaxDensity; +} + +static inline EdgeData calculateEdge(const float &vA, const float &vB, const glm::vec3 &gA, const glm::vec3 &gB, + const ThresholdType &threshold) { + EdgeData edge; + + const float divisor = (vA - vB); + if (divisor == 0.0f) { + edge.fraction = 0.0f; + } else { + edge.fraction = (vA - threshold) / divisor; + } + + if (glm::min(vA, vB) <= threshold && glm::max(vA, vB) > threshold) { + edge.intersects = true; + } else { + edge.intersects = false; + return edge; + } + + edge.normal = (gA * edge.fraction + gB * (1.0f - edge.fraction)); + const float v = glm::length2(edge.normal); + if (v != 0.0f) { + edge.normal *= glm::inversesqrt(v); + } + + return edge; +} + +static inline PositionNormal computeVertex(EdgeData *edges[12]) { + glm::vec3 massPoint{0.0f}; // The average of the intersection vertices + + glm::vec3 vertices[12]; + + vertices[0] = {edges[0]->fraction, 0, 0}; + vertices[1] = {0, edges[1]->fraction, 0}; + vertices[2] = {0, 0, edges[2]->fraction}; + vertices[3] = {1, edges[3]->fraction, 0}; + vertices[4] = {1, 0, edges[4]->fraction}; + vertices[5] = {0, 1, edges[5]->fraction}; + vertices[6] = {edges[6]->fraction, 1, 0}; + vertices[7] = {edges[7]->fraction, 0, 1}; + vertices[8] = {0, edges[8]->fraction, 1}; + vertices[9] = {1, 1, edges[9]->fraction}; + vertices[10] = {1, edges[10]->fraction, 1}; + vertices[11] = {edges[11]->fraction, 1, 1}; + + int numIntersections = 0; + for (int i = 0; i < 12; ++i) { + if (!edges[i]->intersects) { + continue; + } + + ++numIntersections; + massPoint += vertices[i]; + } + + massPoint /= numIntersections; // Make the average + + glm::vec3 cellVertexNormal{0, 0, 0}; + + double matrix[12][3]; + double vector[12]; + int rows = 0; + + for (int i = 0; i < 12; ++i) { + if (!edges[i]->intersects) { + continue; + } + + glm::vec3 normal = edges[i]->normal; + matrix[rows][0] = normal.x; + matrix[rows][1] = normal.y; + matrix[rows][2] = normal.z; + + const glm::vec3 product = normal * (vertices[i] - massPoint); + + vector[rows] = product.x + product.y + product.z; + + cellVertexNormal += normal; + + ++rows; + } + + const glm::vec3 &vertexPosition = evaluateQEF(matrix, vector, rows) + massPoint; + + core_assert_msg(vertexPosition.x > -0.01 && vertexPosition.y > -0.01 && vertexPosition.z > -0.01 && + vertexPosition.x < 1.01 && vertexPosition.y < 1.01 && vertexPosition.z < 1.01, + "Vertex is outside unit cell %f:%f:%f", vertexPosition.x, vertexPosition.y, vertexPosition.z); + + const float v = glm::length2(cellVertexNormal); + if (v != 0.0f) { + cellVertexNormal *= glm::inversesqrt(v); + } + + return {vertexPosition, cellVertexNormal}; +} + +static inline uint32_t convert(uint32_t x, uint32_t y, uint32_t z, uint32_t width, uint32_t height) { + return z * height * width + y * width + x; +} + +} // namespace + +void extractDualContouringMesh(const voxel::RawVolume *volData, const palette::Palette &palette, const Region ®ion, + ChunkMesh *result) { + const float threshold = 127.5f; // controller.getThreshold(); + + const int regionXDimension = region.getDimensionsInVoxels().x; + const int regionYDimension = region.getDimensionsInVoxels().y; + const int regionZDimension = region.getDimensionsInVoxels().z; + + const auto gradientRegionXDimension = regionXDimension + 2; + const auto gradientRegionYDimension = regionYDimension + 2; + const auto gradientRegionZDimension = regionZDimension + 2; + + const size_t gradientsSize = + (size_t)gradientRegionXDimension * (size_t)gradientRegionYDimension * (size_t)gradientRegionZDimension; + core::DynamicArray> gradients; + gradients.reserve(gradientsSize); + + RawVolume::Sampler volSampler{volData}; + volSampler.setPosition(region.getLowerCorner() - 1); + + const auto lowerCornerX = region.getLowerCorner().x; + const auto lowerCornerY = region.getLowerCorner().y; + const auto lowerCornerZ = region.getLowerCorner().z; + + for (int32_t z = 0; z < gradientRegionZDimension; z++) { + volSampler.setPosition(lowerCornerX - 1, lowerCornerY - 1, + lowerCornerZ + z - 1); // Reset x and y and increment z + for (int32_t y = 0; y < gradientRegionYDimension; y++) { + volSampler.setPosition(lowerCornerX - 1, lowerCornerY + y - 1, + lowerCornerZ + z - 1); // Reset x and increment y (z remains the same) + for (int32_t x = 0; x < gradientRegionXDimension; x++) { + volSampler.movePositiveX(); // Increment x + + const float voxel = convertToDensity(volSampler.voxel()); + const float voxel1px = convertToDensity(volSampler.peekVoxel1px0py0pz()); + const float voxel1py = convertToDensity(volSampler.peekVoxel0px1py0pz()); + const float voxel1pz = convertToDensity(volSampler.peekVoxel0px0py1pz()); + const float voxel1nx = convertToDensity(volSampler.peekVoxel1nx0py0pz()); + const float voxel1ny = convertToDensity(volSampler.peekVoxel0px1ny0pz()); + const float voxel1nz = convertToDensity(volSampler.peekVoxel0px0py1nz()); + + gradients.emplace_back(voxel, glm::vec3(voxel1nx - voxel1px, voxel1ny - voxel1py, voxel1nz - voxel1pz)); + } + } + } + + const int32_t cellRegionXDimension = regionXDimension + 2; + const int32_t cellRegionYDimension = regionYDimension + 2; + const int32_t cellRegionZDimension = regionZDimension + 2; + + const size_t cellSize = (size_t)cellRegionXDimension * (size_t)cellRegionYDimension * (size_t)cellRegionZDimension; + core::DynamicArray cells; + cells.reserve(cellSize); + + for (int32_t cellZ = 0; cellZ < cellRegionZDimension; cellZ++) { + for (int32_t cellY = 0; cellY < cellRegionYDimension; cellY++) { + for (int32_t cellX = 0; cellX < cellRegionXDimension; cellX++) { + // For each cell, calculate the edge intersection points and normals + const auto &g000 = gradients[convert(cellX, cellY, cellZ, cellRegionXDimension, cellRegionYDimension)]; + + // For the last columns/rows, only calculate the interior edge + if (cellX < cellRegionXDimension - 1 && cellY < cellRegionYDimension - 1 && + cellZ < cellRegionZDimension - 1) // This is the main bulk + { + const auto &g100 = + gradients[convert(cellX + 1, cellY, cellZ, cellRegionXDimension, cellRegionYDimension)]; + const auto &g010 = + gradients[convert(cellX, cellY + 1, cellZ, cellRegionXDimension, cellRegionYDimension)]; + const auto &g001 = + gradients[convert(cellX, cellY, cellZ + 1, cellRegionXDimension, cellRegionYDimension)]; + cells.push_back({{calculateEdge(g000.first, g100.first, g000.second, g100.second, threshold), + calculateEdge(g000.first, g010.first, g000.second, g010.second, threshold), + calculateEdge(g000.first, g001.first, g000.second, g001.second, threshold)}, + 0}); + } else if (cellX == cellRegionXDimension - 1 || cellY == cellRegionYDimension - 1 || + cellZ == cellRegionZDimension - 1) // This is the three far edges and the far corner + { + cells.push_back({}); // Default and empty + } else if (cellX == cellRegionXDimension - 1) // Far x side + { + const auto &g100 = + gradients[convert(cellX + 1, cellY, cellZ, cellRegionXDimension, cellRegionYDimension)]; + cells.push_back({{calculateEdge(g000.first, g100.first, g000.second, g100.second, threshold), + EdgeData(), EdgeData()}, + 0}); + } else if (cellY == cellRegionYDimension - 1) // Far y side + { + const auto &g010 = + gradients[convert(cellX + 1, cellY, cellZ, cellRegionXDimension, cellRegionYDimension)]; + cells.push_back( + {{EdgeData(), calculateEdge(g000.first, g010.first, g000.second, g010.second, threshold), + EdgeData()}, + 0}); + } else if (cellZ == cellRegionZDimension - 1) // Far z side + { + const auto &g001 = + gradients[convert(cellX + 1, cellY, cellZ, cellRegionXDimension, cellRegionYDimension)]; + cells.push_back({{EdgeData(), EdgeData(), + calculateEdge(g000.first, g001.first, g000.second, g001.second, threshold)}, + 0}); + } + } + } + } + + EdgeData *edges[12]; // Create this now but it will be overwritten for each cell + + for (int32_t cellZ = 0; cellZ < cellRegionZDimension; cellZ++) { + for (int32_t cellY = 0; cellY < cellRegionYDimension; cellY++) { + for (int32_t cellX = 0; cellX < cellRegionXDimension; cellX++) { + if (cellZ >= 1 && cellY >= 1 && cellX >= 1) { + // After the first rows and columns are done, start calculating vertex positions + const int32_t cellXVertex = cellX - 1; + const int32_t cellYVertex = cellY - 1; + const int32_t cellZVertex = cellZ - 1; + + auto &cell = cells[convert(cellXVertex, cellYVertex, cellZVertex, cellRegionXDimension, + cellRegionYDimension)]; + + edges[0] = &cell.edges[0]; + edges[1] = &cell.edges[1]; + edges[2] = &cell.edges[2]; + + edges[3] = &cells[convert(cellXVertex + 1, cellYVertex, cellZVertex, cellRegionXDimension, + cellRegionYDimension)] + .edges[1]; + edges[4] = &cells[convert(cellXVertex + 1, cellYVertex, cellZVertex, cellRegionXDimension, + cellRegionYDimension)] + .edges[2]; + + edges[5] = &cells[convert(cellXVertex, cellYVertex + 1, cellZVertex, cellRegionXDimension, + cellRegionYDimension)] + .edges[2]; + edges[6] = &cells[convert(cellXVertex, cellYVertex + 1, cellZVertex, cellRegionXDimension, + cellRegionYDimension)] + .edges[0]; + + edges[7] = &cells[convert(cellXVertex, cellYVertex, cellZVertex + 1, cellRegionXDimension, + cellRegionYDimension)] + .edges[0]; + edges[8] = &cells[convert(cellXVertex, cellYVertex, cellZVertex + 1, cellRegionXDimension, + cellRegionYDimension)] + .edges[1]; + + edges[9] = &cells[convert(cellXVertex + 1, cellYVertex + 1, cellZVertex, cellRegionXDimension, + cellRegionYDimension)] + .edges[2]; + + edges[10] = &cells[convert(cellXVertex + 1, cellYVertex, cellZVertex + 1, cellRegionXDimension, + cellRegionYDimension)] + .edges[1]; + + edges[11] = &cells[convert(cellXVertex, cellYVertex + 1, cellZVertex + 1, cellRegionXDimension, + cellRegionYDimension)] + .edges[0]; + + if (edges[0]->intersects || edges[1]->intersects || edges[2]->intersects || edges[3]->intersects || + edges[4]->intersects || edges[5]->intersects || edges[6]->intersects || edges[7]->intersects || + edges[8]->intersects || edges[9]->intersects || edges[10]->intersects || + edges[11]->intersects) //'if' Maybe not needed? + { + const PositionNormal &vertex = computeVertex(edges); + + VoxelVertex v; + v.info = 0; + + // TODO: use the sampler + const int x = region.getLowerCorner().x + cellXVertex; + const int y = region.getLowerCorner().y + cellYVertex; + const int z = region.getLowerCorner().z + cellZVertex; + v.colorIndex = volData->voxel(x, y, z).getColor(); + v.position.x = vertex.position.x + (float)cellXVertex; + v.position.y = vertex.position.y + (float)cellYVertex; + v.position.z = vertex.position.z + (float)cellZVertex; + + cell.vertexIndex = result->mesh[0].addVertex(v); + result->mesh[0].setNormal(cell.vertexIndex, vertex.normal); + + // TODO: transparency mesh + if (cellZVertex >= 1 && cellYVertex >= 1 && cellXVertex >= 1) { + // Once the second rows and colums are done, start connecting up edges + if (cell.edges[0].intersects) { + const auto &v1 = cells[convert(cellXVertex, cellYVertex - 1, cellZVertex, + cellRegionXDimension, cellRegionYDimension)]; + const auto &v2 = cells[convert(cellXVertex, cellYVertex, cellZVertex - 1, + cellRegionXDimension, cellRegionYDimension)]; + const auto &v3 = cells[convert(cellXVertex, cellYVertex - 1, cellZVertex - 1, + cellRegionXDimension, cellRegionYDimension)]; + result->mesh[0].addTriangle(cell.vertexIndex, v1.vertexIndex, v2.vertexIndex); + result->mesh[0].addTriangle(v3.vertexIndex, v2.vertexIndex, v1.vertexIndex); + } + + if (cell.edges[1].intersects) { + const auto &v1 = cells[convert(cellXVertex - 1, cellYVertex, cellZVertex, + cellRegionXDimension, cellRegionYDimension)]; + const auto &v2 = cells[convert(cellXVertex, cellYVertex, cellZVertex - 1, + cellRegionXDimension, cellRegionYDimension)]; + const auto &v3 = cells[convert(cellXVertex - 1, cellYVertex, cellZVertex - 1, + cellRegionXDimension, cellRegionYDimension)]; + result->mesh[0].addTriangle(cell.vertexIndex, v1.vertexIndex, v2.vertexIndex); + result->mesh[0].addTriangle(v3.vertexIndex, v2.vertexIndex, v1.vertexIndex); + } + + if (cell.edges[2].intersects) { + const auto &v1 = cells[convert(cellXVertex - 1, cellYVertex, cellZVertex, + cellRegionXDimension, cellRegionYDimension)]; + const auto &v2 = cells[convert(cellXVertex, cellYVertex - 1, cellZVertex, + cellRegionXDimension, cellRegionYDimension)]; + const auto &v3 = cells[convert(cellXVertex - 1, cellYVertex - 1, cellZVertex, + cellRegionXDimension, cellRegionYDimension)]; + result->mesh[0].addTriangle(cell.vertexIndex, v1.vertexIndex, v2.vertexIndex); + result->mesh[0].addTriangle(v3.vertexIndex, v2.vertexIndex, v1.vertexIndex); + } + } + } + } + } + } + } +} + +} // namespace voxel diff --git a/src/modules/voxel/private/DualContouringSurfaceExtractor.h b/src/modules/voxel/private/DualContouringSurfaceExtractor.h new file mode 100644 index 0000000000..0cd339a690 --- /dev/null +++ b/src/modules/voxel/private/DualContouringSurfaceExtractor.h @@ -0,0 +1,17 @@ +/** + * @file + */ + +#pragma once + +namespace palette { +class Palette; +} + +namespace voxel { +struct ChunkMesh; +class RawVolume; +class Region; +void extractDualContouringMesh(const voxel::RawVolume *volData, const palette::Palette &palette, const Region ®ion, ChunkMesh *mesh); + +} // namespace voxel diff --git a/src/modules/voxelformat/FormatConfig.cpp b/src/modules/voxelformat/FormatConfig.cpp index d57e0fd213..d43247af9d 100644 --- a/src/modules/voxelformat/FormatConfig.cpp +++ b/src/modules/voxelformat/FormatConfig.cpp @@ -23,7 +23,7 @@ bool FormatConfig::init() { core::Var::get(cfg::VoxformatMergequads, "true", core::CV_NOPERSIST, "Merge similar quads to optimize the mesh", core::Var::boolValidator); core::Var::get(cfg::VoxelMeshMode, core::string::toString((int)voxel::SurfaceExtractionType::Cubic), - core::CV_SHADER, "0 = cubes, 1 = marching cubes", + core::CV_SHADER, "0 = cubes, 1 = marching cubes, 2 = dual contouring", core::Var::minMaxValidator<(int)voxel::SurfaceExtractionType::Cubic, (int)voxel::SurfaceExtractionType::Max - 1>); core::Var::get(cfg::VoxformatReusevertices, "true", core::CV_NOPERSIST, "Reuse vertices or always create new ones", diff --git a/src/tools/voxedit/modules/voxedit-ui/MenuBar.cpp b/src/tools/voxedit/modules/voxedit-ui/MenuBar.cpp index ecc4cfa574..9be8cec6da 100644 --- a/src/tools/voxedit/modules/voxedit-ui/MenuBar.cpp +++ b/src/tools/voxedit/modules/voxedit-ui/MenuBar.cpp @@ -116,7 +116,7 @@ bool MenuBar::update(ui::IMGUIApp *app, command::CommandExecutionListener &liste ui::metricOption(); static const core::Array meshModes = { - _("Cubes"), _("Marching cubes")}; + _("Cubes"), _("Marching cubes"), _("Dual contouring")}; ImGui::ComboVar(_("Mesh mode"), cfg::VoxelMeshMode, meshModes); ImGui::InputVarInt(_("Model animation speed"), cfg::VoxEditAnimationSpeed); ImGui::InputVarInt(_("Autosave delay in seconds"), cfg::VoxEditAutoSaveSeconds);