From 6b53b0bb551c0b764a746cd8ea83ca695cd710f9 Mon Sep 17 00:00:00 2001 From: Ziheng Liang Date: Sun, 28 Jan 2018 20:17:05 -0500 Subject: [PATCH 1/4] compile version --- src/fd_grad.cpp | 16 ++++++++++++++++ src/fd_interpolate.cpp | 25 +++++++++++++++++++++++++ src/fd_partial_derivative.cpp | 26 ++++++++++++++++++++++++++ src/poisson_surface_reconstruction.cpp | 21 +++++++++++++++++++++ 4 files changed, 88 insertions(+) diff --git a/src/fd_grad.cpp b/src/fd_grad.cpp index 9206181..0f43f9c 100644 --- a/src/fd_grad.cpp +++ b/src/fd_grad.cpp @@ -1,4 +1,5 @@ #include "fd_grad.h" +#include "fd_partial_derivative.h" void fd_grad( const int nx, @@ -9,5 +10,20 @@ void fd_grad( { //////////////////////////////////////////////////////////////////////////// // Add your code here + int dx = (nx-1)*ny*nz; + int dy = nx*(ny-1)*nz; + int dz = nx*ny*(nz-1); + + Eigen::SparseMatrix Dx = Eigen::SparseMatrix(dx, nx*ny*nz); + Eigen::SparseMatrix Dy = Eigen::SparseMatrix(dy, nx*ny*nz); + Eigen::SparseMatrix Dz = Eigen::SparseMatrix(dz, nx*ny*nz); + + fd_partial_derivative(nx, ny, nz, h, 0, Dx); + fd_partial_derivative(nx, ny, nz, h, 1, Dy); + fd_partial_derivative(nx, ny, nz, h, 2, Dz); + + // G.block(0, 0, dx, nx*ny*nz) = Dx; + // G.block(dx, nx*ny*nz, dx+dy, nx*ny*nz) = Dy; + // G.block(dx+dy, nx*ny*nz, dx+dy+dz, nx*ny*nz) = Dz; //////////////////////////////////////////////////////////////////////////// } diff --git a/src/fd_interpolate.cpp b/src/fd_interpolate.cpp index 2a32675..0eea826 100644 --- a/src/fd_interpolate.cpp +++ b/src/fd_interpolate.cpp @@ -11,5 +11,30 @@ void fd_interpolate( { //////////////////////////////////////////////////////////////////////////// // Add your code here + int i_int, j_int, k_int; + double i, j, k; + double w_x, w_y, w_z; + for (int row=0; row((nx-x_diff)*(ny-y_diff)*(nz-z_diff), nx*ny*nz); + for(int i = 0; i < nx - x_diff; i++) + { + for(int j = 0; j < ny - y_diff; j++) + { + for(int k = 0; k < nz - z_diff; k++) + { + // Convert subscript to index + const auto ind = i + nx*(j + k * ny); + const auto ind_next = (x_diff + nx*(y_diff + z_diff * ny)); + D.insert(ind, ind) = -1; + D.insert(ind, ind_next) = 1; + } + } + } //////////////////////////////////////////////////////////////////////////// } diff --git a/src/poisson_surface_reconstruction.cpp b/src/poisson_surface_reconstruction.cpp index 4fe7c1e..ec79ba4 100644 --- a/src/poisson_surface_reconstruction.cpp +++ b/src/poisson_surface_reconstruction.cpp @@ -1,6 +1,10 @@ #include "poisson_surface_reconstruction.h" #include #include +#include +#include "fd_grad.h" +#include "fd_interpolate.h" +#include void poisson_surface_reconstruction( const Eigen::MatrixXd & P, @@ -47,7 +51,24 @@ void poisson_surface_reconstruction( //////////////////////////////////////////////////////////////////////////// // Add your code here //////////////////////////////////////////////////////////////////////////// + int dx = (nx-1)*ny*nz; + int dy = nx*(ny-1)*nz; + int dz = nx*ny*(nz-1); + Eigen::SparseMatrix W(P.rows(), nx*ny*nz); + Eigen::SparseMatrix G(dx+dy+dz, nx*ny*nz); + fd_interpolate(nx, ny, nz, h, corner, P, W); + fd_grad(nx, ny, nz, h, G); + Eigen::VectorXd v; + Eigen::MatrixXd temp = W * N; + v << temp.col(0), temp.col(1), temp.col(2); + // Eigen::Map v((W * N).data(), (W * N).size()); + + Eigen::SimplicialLDLT> solver; + solver.compute(G.transpose() * G); + g = solver.solve(G.transpose() * v); + + // g -= W * g; //////////////////////////////////////////////////////////////////////////// // Run black box algorithm to compute mesh from implicit function: this // function always extracts g=0, so "pre-shift" your g values by -sigma From 52f5ce57389f92a232b2077420d1feed4a51a5c8 Mon Sep 17 00:00:00 2001 From: Ziheng Liang Date: Sun, 28 Jan 2018 20:31:53 -0500 Subject: [PATCH 2/4] change .block() to concatenate --- src/fd_grad.cpp | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/src/fd_grad.cpp b/src/fd_grad.cpp index 0f43f9c..2e61371 100644 --- a/src/fd_grad.cpp +++ b/src/fd_grad.cpp @@ -22,8 +22,18 @@ void fd_grad( fd_partial_derivative(nx, ny, nz, h, 1, Dy); fd_partial_derivative(nx, ny, nz, h, 2, Dz); - // G.block(0, 0, dx, nx*ny*nz) = Dx; - // G.block(dx, nx*ny*nz, dx+dy, nx*ny*nz) = Dy; - // G.block(dx+dy, nx*ny*nz, dx+dy+dz, nx*ny*nz) = Dz; + //concatenate method copy from stackoverflow + G.reserve(Dx.nonZeros() + Dy.nonZeros() + Dz.nonZeros()); + for(int c=0; c::InnerIterator itDx(Dx, c); itDx; ++itDx) + G.insertBack(itDx.row(), c) = itDx.value(); + for(Eigen::SparseMatrix::InnerIterator itDy(Dy, c); itDy; ++itDy) + G.insertBack(itDy.row(), c) = itDy.value(); + for(Eigen::SparseMatrix::InnerIterator itDz(Dz, c); itDz; ++itDz) + G.insertBack(itDz.row(), c) = itDz.value(); + } + + G.finalize(); + //////////////////////////////////////////////////////////////////////////// } From 77abbd0a98754f85590f93326e1e2d16536db32e Mon Sep 17 00:00:00 2001 From: Ziheng Liang Date: Tue, 30 Jan 2018 00:10:38 -0500 Subject: [PATCH 3/4] save --- include/fd_grad.h | 6 +++ src/fd_grad.cpp | 40 ++++++++++------ src/fd_interpolate.cpp | 65 ++++++++++++++++++++------ src/fd_partial_derivative.cpp | 10 ++-- src/poisson_surface_reconstruction.cpp | 51 ++++++++++++++++---- 5 files changed, 129 insertions(+), 43 deletions(-) diff --git a/include/fd_grad.h b/include/fd_grad.h index fa0df87..5ed1b36 100644 --- a/include/fd_grad.h +++ b/include/fd_grad.h @@ -19,4 +19,10 @@ void fd_grad( const int nz, const double h, Eigen::SparseMatrix & G); + +void stack_sparse_matric( + Eigen::SparseMatrix &A, + Eigen::SparseMatrix &B, + Eigen::SparseMatrix &C, + Eigen::SparseMatrix &out); #endif diff --git a/src/fd_grad.cpp b/src/fd_grad.cpp index 2e61371..3d2a5b6 100644 --- a/src/fd_grad.cpp +++ b/src/fd_grad.cpp @@ -1,5 +1,6 @@ #include "fd_grad.h" #include "fd_partial_derivative.h" +#include void fd_grad( const int nx, @@ -21,19 +22,32 @@ void fd_grad( fd_partial_derivative(nx, ny, nz, h, 0, Dx); fd_partial_derivative(nx, ny, nz, h, 1, Dy); fd_partial_derivative(nx, ny, nz, h, 2, Dz); + stack_sparse_matric(Dx, Dy, Dz, G); + //////////////////////////////////////////////////////////////////////////// +} - //concatenate method copy from stackoverflow - G.reserve(Dx.nonZeros() + Dy.nonZeros() + Dz.nonZeros()); - for(int c=0; c::InnerIterator itDx(Dx, c); itDx; ++itDx) - G.insertBack(itDx.row(), c) = itDx.value(); - for(Eigen::SparseMatrix::InnerIterator itDy(Dy, c); itDy; ++itDy) - G.insertBack(itDy.row(), c) = itDy.value(); - for(Eigen::SparseMatrix::InnerIterator itDz(Dz, c); itDz; ++itDz) - G.insertBack(itDz.row(), c) = itDz.value(); +//concatenate method copy from stackoverflow +void stack_sparse_matric( + Eigen::SparseMatrix &A, + Eigen::SparseMatrix &B, + Eigen::SparseMatrix &C, + Eigen::SparseMatrix &out) +{ + std::vector > tripletList; + for (int k = 0; k < A.outerSize(); ++k) { + for (Eigen::SparseMatrix::InnerIterator it(A, k); it; ++it) { + tripletList.push_back(Eigen::Triplet(it.row(), it.col(), it.value())); + } } - - G.finalize(); - - //////////////////////////////////////////////////////////////////////////// + for (int k = 0; k < B.outerSize(); ++k) { + for (Eigen::SparseMatrix::InnerIterator it(B, k); it; ++it) { + tripletList.push_back(Eigen::Triplet(A.rows() + it.row(), it.col(), it.value())); + } + } + for (int k = 0; k < C.outerSize(); ++k) { + for (Eigen::SparseMatrix::InnerIterator it(C, k); it; ++it) { + tripletList.push_back(Eigen::Triplet(A.rows() + B.rows() + it.row(), it.col(), it.value())); + } + } + out.setFromTriplets(tripletList.begin(), tripletList.end()); } diff --git a/src/fd_interpolate.cpp b/src/fd_interpolate.cpp index 0eea826..fb54c37 100644 --- a/src/fd_interpolate.cpp +++ b/src/fd_interpolate.cpp @@ -1,4 +1,7 @@ #include "fd_interpolate.h" +#include +#include +#include void fd_interpolate( const int nx, @@ -14,27 +17,59 @@ void fd_interpolate( int i_int, j_int, k_int; double i, j, k; double w_x, w_y, w_z; - for (int row=0; row (); + // Eigen::MatrixXd P_frac = temp - P_int.cast (); + + std::vector > tripletList; + + for (int row=0; row(row, col, P_frac(row,0) * P_frac(row,1) * P_frac(row,2))); + // tripletList.push_back(Eigen::Triplet(row, col+1, (1-P_frac(row,0)) * P_frac(row,1) * P_frac(row,2))); + // tripletList.push_back(Eigen::Triplet(row, col+nx, P_frac(row,0) * (1-P_frac(row,1)) * P_frac(row,2))); + // tripletList.push_back(Eigen::Triplet(row, col+nx*ny, P_frac(row,0) * P_frac(row,1) * (1-P_frac(row,2)))); + // tripletList.push_back(Eigen::Triplet(row, col+1+nx, (1-P_frac(row,0)) * (1-P_frac(row,1)) * P_frac(row,2))); + // tripletList.push_back(Eigen::Triplet(row, col+1+nx*ny, (1-P_frac(row,0)) * P_frac(row,1) * (1-P_frac(row,2)))); + // tripletList.push_back(Eigen::Triplet(row, col+nx+nx*ny, (P_frac(row,0) * (1-P_frac(row,1)) * (1-P_frac(row,2))))); + // tripletList.push_back(Eigen::Triplet(row, col+1+nx+nx*ny, (1-P_frac(row,0)) * (1-P_frac(row,1)) * (1-P_frac(row,2)))); + + + tripletList.push_back(Eigen::Triplet(row, col, (1-w_x) * (1-w_y) * (1-w_z))); + tripletList.push_back(Eigen::Triplet(row, col+1, w_x * (1-w_y) * (1-w_z))); + tripletList.push_back(Eigen::Triplet(row, col+nx, (1-w_x) * w_y * (1-w_z))); + tripletList.push_back(Eigen::Triplet(row, col+nx*ny, (1-w_x) * (1-w_y) * w_z)); + tripletList.push_back(Eigen::Triplet(row, col+1+nx, w_x * w_y * (1-w_z))); + tripletList.push_back(Eigen::Triplet(row, col+1+nx*ny, w_x * (1-w_y) * w_z)); + tripletList.push_back(Eigen::Triplet(row, col+nx+nx*ny, (1-w_x) * w_y * w_z)); + tripletList.push_back(Eigen::Triplet(row, col+1+nx+nx*ny, w_x * w_y * w_z)); + } + W.setFromTriplets(tripletList.begin(), tripletList.end()); //////////////////////////////////////////////////////////////////////////// } diff --git a/src/fd_partial_derivative.cpp b/src/fd_partial_derivative.cpp index 30d7a61..ff0780e 100644 --- a/src/fd_partial_derivative.cpp +++ b/src/fd_partial_derivative.cpp @@ -13,7 +13,6 @@ void fd_partial_derivative( int x_diff = 0, y_diff = 0, z_diff = 0; if (dir == 0) { x_diff = 1; - } else if (dir == 1) { y_diff = 1; @@ -21,7 +20,9 @@ void fd_partial_derivative( else if (dir == 2) { z_diff = 1; } + D = Eigen::SparseMatrix((nx-x_diff)*(ny-y_diff)*(nz-z_diff), nx*ny*nz); + const auto offset = (x_diff + nx*(y_diff + z_diff * ny)); for(int i = 0; i < nx - x_diff; i++) { for(int j = 0; j < ny - y_diff; j++) @@ -29,10 +30,9 @@ void fd_partial_derivative( for(int k = 0; k < nz - z_diff; k++) { // Convert subscript to index - const auto ind = i + nx*(j + k * ny); - const auto ind_next = (x_diff + nx*(y_diff + z_diff * ny)); - D.insert(ind, ind) = -1; - D.insert(ind, ind_next) = 1; + const auto ind = i + (nx-x_diff)*(j + k * (ny-y_diff)); + D.insert(ind, ind) = -1/h; + D.insert(ind, ind + offset) = 1/h; } } } diff --git a/src/poisson_surface_reconstruction.cpp b/src/poisson_surface_reconstruction.cpp index ec79ba4..ac2cf2b 100644 --- a/src/poisson_surface_reconstruction.cpp +++ b/src/poisson_surface_reconstruction.cpp @@ -4,7 +4,7 @@ #include #include "fd_grad.h" #include "fd_interpolate.h" -#include +#include void poisson_surface_reconstruction( const Eigen::MatrixXd & P, @@ -55,20 +55,51 @@ void poisson_surface_reconstruction( int dy = nx*(ny-1)*nz; int dz = nx*ny*(nz-1); + Eigen::SparseMatrix Wx(P.rows(), dx); + Eigen::SparseMatrix Wy(P.rows(), dy); + Eigen::SparseMatrix Wz(P.rows(), dz); Eigen::SparseMatrix W(P.rows(), nx*ny*nz); Eigen::SparseMatrix G(dx+dy+dz, nx*ny*nz); - fd_interpolate(nx, ny, nz, h, corner, P, W); + std::cout << "test1" << std::endl; fd_grad(nx, ny, nz, h, G); - Eigen::VectorXd v; - Eigen::MatrixXd temp = W * N; - v << temp.col(0), temp.col(1), temp.col(2); - // Eigen::Map v((W * N).data(), (W * N).size()); + std::cout << "test2" << std::endl; + fd_interpolate(nx-1, ny, nz, h, corner + Eigen::RowVector3d(h/2, 0, 0), P, Wx); + std::cout << "test2.1" << std::endl; + fd_interpolate(nx, ny-1, nz, h, corner + Eigen::RowVector3d(0, h/2, 0), P, Wy); + std::cout << "test2.2" << std::endl; + fd_interpolate(nx, ny, nz-1, h, corner + Eigen::RowVector3d(0, 0, h/2), P, Wz); + fd_interpolate(nx, ny, nz, h, corner, P, W); + + std::cout << W * x - P << std::endl; + std::cout << (W * x - P).norm() << std::endl; + std::cout << "test3" << std::endl; + Eigen::VectorXd vx, vy, vz; + vx = Wx.transpose() * N.col(0); + vy = Wy.transpose() * N.col(1); + vz = Wz.transpose() * N.col(2); + + + Eigen::VectorXd v(dx+dy+dz); + v << vx, vy, vz; - Eigen::SimplicialLDLT> solver; - solver.compute(G.transpose() * G); - g = solver.solve(G.transpose() * v); + // Eigen::Map v(temp.data(), temp.size()); + // v << temp.col(0), temp.col(1), temp.col(2); + std::cout << "test3.5" << std::endl; + std::cout << "test4" << std::endl; + std::cout << W.rows() << "|" << W.cols() << std::endl; + std::cout << g.rows() << std::endl; + + std::cout << "test5" << std::endl; + Eigen::ConjugateGradient> solver; + Eigen::SparseMatrix G_square = G.transpose() * G; + Eigen::VectorXd Gv = G.transpose() * v; + solver.compute(G_square); + g = solver.solve(Gv); + std::cout << "test6" << std::endl; + + double sigma = (W * g).sum() / n; + g = (g.array() - sigma).matrix(); - // g -= W * g; //////////////////////////////////////////////////////////////////////////// // Run black box algorithm to compute mesh from implicit function: this // function always extracts g=0, so "pre-shift" your g values by -sigma From 0338e10481aeb2aaeb072365acaadc1f8072b13e Mon Sep 17 00:00:00 2001 From: Ziheng Liang Date: Tue, 30 Jan 2018 00:40:58 -0500 Subject: [PATCH 4/4] debug && cleanup --- src/fd_partial_derivative.cpp | 7 ++++--- src/poisson_surface_reconstruction.cpp | 19 +------------------ 2 files changed, 5 insertions(+), 21 deletions(-) diff --git a/src/fd_partial_derivative.cpp b/src/fd_partial_derivative.cpp index ff0780e..f9d7a7e 100644 --- a/src/fd_partial_derivative.cpp +++ b/src/fd_partial_derivative.cpp @@ -30,9 +30,10 @@ void fd_partial_derivative( for(int k = 0; k < nz - z_diff; k++) { // Convert subscript to index - const auto ind = i + (nx-x_diff)*(j + k * (ny-y_diff)); - D.insert(ind, ind) = -1/h; - D.insert(ind, ind + offset) = 1/h; + const auto stagger_idx = i + j * (nx-x_diff) + k * (ny-y_diff) * (nx-x_diff); + const auto original_idx = i + j * nx + k * ny * nx; + D.insert(stagger_idx, original_idx) = -1/h; + D.insert(stagger_idx, original_idx + offset) = 1/h; } } } diff --git a/src/poisson_surface_reconstruction.cpp b/src/poisson_surface_reconstruction.cpp index ac2cf2b..b4e2d28 100644 --- a/src/poisson_surface_reconstruction.cpp +++ b/src/poisson_surface_reconstruction.cpp @@ -60,42 +60,25 @@ void poisson_surface_reconstruction( Eigen::SparseMatrix Wz(P.rows(), dz); Eigen::SparseMatrix W(P.rows(), nx*ny*nz); Eigen::SparseMatrix G(dx+dy+dz, nx*ny*nz); - std::cout << "test1" << std::endl; + fd_grad(nx, ny, nz, h, G); - std::cout << "test2" << std::endl; fd_interpolate(nx-1, ny, nz, h, corner + Eigen::RowVector3d(h/2, 0, 0), P, Wx); - std::cout << "test2.1" << std::endl; fd_interpolate(nx, ny-1, nz, h, corner + Eigen::RowVector3d(0, h/2, 0), P, Wy); - std::cout << "test2.2" << std::endl; fd_interpolate(nx, ny, nz-1, h, corner + Eigen::RowVector3d(0, 0, h/2), P, Wz); fd_interpolate(nx, ny, nz, h, corner, P, W); - std::cout << W * x - P << std::endl; - std::cout << (W * x - P).norm() << std::endl; - std::cout << "test3" << std::endl; Eigen::VectorXd vx, vy, vz; vx = Wx.transpose() * N.col(0); vy = Wy.transpose() * N.col(1); vz = Wz.transpose() * N.col(2); - Eigen::VectorXd v(dx+dy+dz); v << vx, vy, vz; - - // Eigen::Map v(temp.data(), temp.size()); - // v << temp.col(0), temp.col(1), temp.col(2); - std::cout << "test3.5" << std::endl; - std::cout << "test4" << std::endl; - std::cout << W.rows() << "|" << W.cols() << std::endl; - std::cout << g.rows() << std::endl; - - std::cout << "test5" << std::endl; Eigen::ConjugateGradient> solver; Eigen::SparseMatrix G_square = G.transpose() * G; Eigen::VectorXd Gv = G.transpose() * v; solver.compute(G_square); g = solver.solve(Gv); - std::cout << "test6" << std::endl; double sigma = (W * g).sum() / n; g = (g.array() - sigma).matrix();