From bbaba4051817d4e9e9a96d7e9c515bbe8d38c390 Mon Sep 17 00:00:00 2001 From: Minhua Zhu Date: Fri, 28 Sep 2018 18:32:04 -0400 Subject: [PATCH] Minhua Zhu HW1 --- src/fd_grad.cpp | 11 ++++++ src/fd_interpolate.cpp | 38 ++++++++++++++++++ src/fd_partial_derivative.cpp | 53 +++++++++++++++++++++++++- src/poisson_surface_reconstruction.cpp | 48 +++++++++++++++++++++++ 4 files changed, 149 insertions(+), 1 deletion(-) diff --git a/src/fd_grad.cpp b/src/fd_grad.cpp index 9206181..4e5f39d 100644 --- a/src/fd_grad.cpp +++ b/src/fd_grad.cpp @@ -1,4 +1,6 @@ #include "fd_grad.h" +#include "fd_partial_derivative.h" +#include void fd_grad( const int nx, @@ -10,4 +12,13 @@ void fd_grad( //////////////////////////////////////////////////////////////////////////// // Add your code here //////////////////////////////////////////////////////////////////////////// + Eigen::SparseMatrix Dx, Dy, Dz, Dxy; + + 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); + + igl::cat(1,Dx,Dy,Dxy); + igl::cat(1,Dxy,Dz,G); + } diff --git a/src/fd_interpolate.cpp b/src/fd_interpolate.cpp index 2a32675..aa3c2f1 100644 --- a/src/fd_interpolate.cpp +++ b/src/fd_interpolate.cpp @@ -1,4 +1,6 @@ #include "fd_interpolate.h" +#include +#include void fd_interpolate( const int nx, @@ -12,4 +14,40 @@ void fd_interpolate( //////////////////////////////////////////////////////////////////////////// // Add your code here //////////////////////////////////////////////////////////////////////////// + + W.resize(P.rows(),nx * ny * nz); + typedef Eigen::Triplet T; + std::vector< T > tripletList; + tripletList.reserve(8*P.rows()); + + for (int i = 0; i < P.rows(); i++) { + // distance + double dx = P(i,0) - corner[0]; + double dy = P(i,1) - corner[1]; + double dz = P(i,2) - corner[2]; + + // position of grids + int gx = floor(dx / h); + int gy = floor(dy / h); + int gz = floor(dz / h); + + // relative distance + double rx = (dx - gx*h) / h; + double ry = (dy - gy*h) / h ; + double rz = (dz - gz*h) / h; + + // eight point of one cube + tripletList.push_back(T(i, gx+nx*(gy+gz*ny), (1-rx)*(1-ry)*(1-rz))); + tripletList.push_back(T(i, gx+1+nx*(gy+gz*ny), rx*(1-ry)*(1-rz))); + tripletList.push_back(T(i, gx+nx*(gy+1+gz*ny), (1-rx)*ry*(1-rz))); + + tripletList.push_back(T(i, gx+nx*(gy+(gz+1)*ny), (1-rx)*(1-ry)*rz)); + tripletList.push_back(T(i, gx+1+nx*(gy+1+gz*ny), rx*ry*(1-rz))); + tripletList.push_back(T(i, gx+1+nx*(gy+(gz+1)*ny), rx*(1-ry)*rz)); + + tripletList.push_back(T(i, gx+nx*(gy+1+(gz+1)*ny), (1-rx)*ry*rz)); + tripletList.push_back(T(i, gx+1+nx*(gy+1+(gz+1)*ny), rx*ry*rz)); + + } + W.setFromTriplets(tripletList.begin(), tripletList.end()); } diff --git a/src/fd_partial_derivative.cpp b/src/fd_partial_derivative.cpp index c3c4deb..7ae335a 100644 --- a/src/fd_partial_derivative.cpp +++ b/src/fd_partial_derivative.cpp @@ -11,4 +11,55 @@ void fd_partial_derivative( //////////////////////////////////////////////////////////////////////////// // Add your code here //////////////////////////////////////////////////////////////////////////// -} + int m = 0; + if (dir == 0) { // nx-1 + int change_x = nx - 1; + m = change_x*ny*nz; + D.resize(m, nx*ny*nz); + for (int i=0; i #include +#include "fd_partial_derivative.h" +#include "fd_grad.h" +#include "fd_interpolate.h" +#include +#include +#include void poisson_surface_reconstruction( const Eigen::MatrixXd & P, @@ -44,6 +50,8 @@ void poisson_surface_reconstruction( } Eigen::VectorXd g = Eigen::VectorXd::Zero(nx*ny*nz); + + //////////////////////////////////////////////////////////////////////////// // Add your code here //////////////////////////////////////////////////////////////////////////// @@ -52,5 +60,45 @@ void poisson_surface_reconstruction( // Run black box algorithm to compute mesh from implicit function: this // function always extracts g=0, so "pre-shift" your g values by -sigma //////////////////////////////////////////////////////////////////////////// + + Eigen::SparseMatrix Weightx(n,(nx-1)*ny*nz); + Eigen::RowVector3d cornerx = corner; + cornerx[0] = corner[0] + h/2; + fd_interpolate(nx-1,ny,nz,h,cornerx,P,Weightx); + + Eigen::SparseMatrix Weighty(n,nx*(ny-1)*nz); + Eigen::RowVector3d cornery = corner; + cornery[1] = corner[1] + h/2; + fd_interpolate(nx,ny-1,nz,h,corner,P,Weighty); + + Eigen::SparseMatrix Weightz(n,nx*ny*(nz-1)); + Eigen::RowVector3d cornerz = corner; + cornerz[2] = corner[2] + h/2; + fd_interpolate(nx,ny,nz-1,h,corner,P,Weightz); + + Eigen::SparseMatrix W(n,nx*ny*nz); + fd_interpolate(nx,ny,nz,h,corner,P, W); + + Eigen::SparseMatrix Gradient; + fd_grad(nx,ny,nz,h,Gradient); + + Eigen::VectorXd vx(n,1), vy(n,1), vz(n,1); + vx=Weightx.transpose()*N.col(0); + vy=Weighty.transpose()*N.col(1); + vz=Weightz.transpose()*N.col(2); + Eigen::VectorXd B((nx-1)*ny*nz+ nx*(ny-1)*nz+ nx*ny*(nz-1)); + B << vx,vy,vz; + + Eigen::SparseMatrix A = Gradient.transpose()*Gradient; + Eigen::VectorXd b = Gradient.transpose()*B; + //Eigen::BiCGSTAB> bicg; + Eigen::ConjugateGradient, Eigen::UpLoType::Lower | Eigen::UpLoType::Upper> bicg; + g = bicg.compute(A).solve(b); + + Eigen::VectorXd i(g.size()); + i.setOnes(); + double sigma = (1/(n*1.0)) * (W*g).sum(); + g = g - sigma*i; + igl::copyleft::marching_cubes(g, x, nx, ny, nz, V, F); }