submit

src/fd_grad.cpp

@@ -1,4 +1,6 @@
 #include "fd_grad.h"
+#include "fd_partial_derivative.h"
+#include <igl/cat.h>
 
 void fd_grad(
   const int nx,
@@ -7,7 +9,10 @@ void fd_grad(
   const double h,
   Eigen::SparseMatrix<double> & G)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+  Eigen::SparseMatrix<double> G1,G2,G3;
+  fd_partial_derivative(nx,ny,nz,h,0,G1);
+  fd_partial_derivative(nx,ny,nz,h,1,G2);
+  fd_partial_derivative(nx,ny,nz,h,2,G3);
+  G=igl::cat(1,G1,G2);
+  G=igl::cat(1,G,G3);
 }

src/fd_interpolate.cpp

@@ -1,5 +1,8 @@
 #include "fd_interpolate.h"
-
+#include <math.h>
+#include <vector>
+#include <iostream>
+typedef Eigen::Triplet<double> T;
 void fd_interpolate(
   const int nx,
   const int ny,
@@ -9,7 +12,31 @@ void fd_interpolate(
   const Eigen::MatrixXd & P,
   Eigen::SparseMatrix<double> & W)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+  int n=P.rows();
+  double x,y,z,dx,dy,dz;
+  std::vector<T> tripletList;
+  tripletList.reserve(n*8);
+  W.resize(n,nx*ny*nz);
+  for (int i = 0; i< n; i++){
+    x=floor((P(i,0)-corner(0))/h); 
+    y=floor((P(i,1)-corner(1))/h);
+    z=floor((P(i, 2)-corner(2))/h);
+    dx=(P(i,0)-corner(0))/h-x;
+    dy=(P(i,1)-corner(1))/h-y;
+    dz=(P(i,2)-corner(2))/h-z;
+
+    tripletList.push_back(T(i, x + y *nx + z * nx*ny, (1-dx) * (1-dy) * (1-dz)));
+    tripletList.push_back(T(i, (x+1) + y *nx + z * nx*ny, dx * (1-dy) * (1-dz)));
+
+    tripletList.push_back(T(i, x + (y+1) *nx + z * nx*ny, (1-dx) * dy * (1-dz)));
+    tripletList.push_back(T(i, x + y * nx + (z+1) * nx*ny, (1-dx) * (1-dy) * dz));
+
+    tripletList.push_back(T(i, (x+1) + (y+1) *nx + z * nx*ny, dx * dy * (1-dz)));
+    tripletList.push_back(T(i, (x+1) + y *nx + (z+1) * nx*ny, dx  * (1-dy) * dz));
+
+    tripletList.push_back(T(i, x + (y+1) *nx + (z+1) * nx*ny, (1-dx) * dy * dz));
+    tripletList.push_back(T(i, (x+1) + (y+1) * nx + (z+1) * nx*ny, dx * dy * dz));
+
+  }
+  W.setFromTriplets(tripletList.begin(),tripletList.end());
 }

src/fd_partial_derivative.cpp

@@ -8,7 +8,38 @@ void fd_partial_derivative(
   const int dir,
   Eigen::SparseMatrix<double> & D)
 {
-  ////////////////////////////////////////////////////////////////////////////
-  // Add your code here
-  ////////////////////////////////////////////////////////////////////////////
+  double h1=1/h;
+  if (dir==0){
+    D.resize((nx-1)*ny*nz,nx*ny*nz);
+    for (int i=0; i<nx-1; i++){
+      //printf("%d ",i);
+      for (int j=0; j<ny; j++)
+        for (int k=0; k<nz; k++){
+          D.insert(i+j*(nx-1)+k*(nx-1)*ny,i+j*nx+k*nx*ny)=-1;
+          D.insert(i+j*(nx-1)+k*(nx-1)*ny,i+1+j*nx+k*nx*ny)=1;
+        }
+    }
+  }
+  if (dir==1){
+    D.resize(nx*(ny-1)*nz,nx*ny*nz);
+    for (int i=0; i<nx; i++){
+      //printf("%d ",i);
+      for (int j=0; j<ny-1; j++)
+        for (int k=0; k<nz; k++){
+          D.insert(i+j*nx+k*nx*(ny-1),i+j*nx+k*nx*ny)=-1;
+          D.insert(i+j*nx+k*nx*(ny-1),i+(j+1)*nx+k*nx*ny)=1;
+        }
+    }
+  }
+  if (dir==2){
+    D.resize(nx*ny*(nz-1),nx*ny*nz);
+    for (int i=0; i<nx; i++){
+      //printf("%d ",i);
+      for (int j=0; j<ny; j++)
+        for (int k=0; k<nz-1; k++){
+          D.insert(i+j*nx+k*nx*ny,i+j*nx+k*nx*ny)=-1;
+          D.insert(i+j*nx+k*nx*ny,i+j*nx+(k+1)*nx*ny)=1;
+        }
+    }
+  }
 }

src/poisson_surface_reconstruction.cpp

@@ -1,6 +1,10 @@
 #include "poisson_surface_reconstruction.h"
 #include <igl/copyleft/marching_cubes.h>
 #include <algorithm>
+#include <igl/cat.h>
+#include "fd_interpolate.h"
+#include "fd_grad.h"
+#include <fstream>
 
 void poisson_surface_reconstruction(
     const Eigen::MatrixXd & P,
@@ -43,7 +47,30 @@ void poisson_surface_reconstruction(
     }
   }
   Eigen::VectorXd g = Eigen::VectorXd::Zero(nx*ny*nz);
+  //printf("%d %d %d %d ",n,nx,ny,nz);
+  Eigen::RowVector3d c1(h/2,0,0),c2(0,h/2,0),c3(0,0,h/2);
+  Eigen::SparseMatrix<double> W1,W2,W3,W;
+  fd_interpolate(nx-1,ny,nz,h,corner+c1,P,W1);
+  fd_interpolate(nx,ny-1,nz,h,corner+c2,P,W2);
+  fd_interpolate(nx,ny,nz-1,h,corner+c3,P,W3);
+  Eigen::VectorXd V1((nx-1)*ny*nz),V2(nx*(ny-1)*nz),V3(nx*ny*(nz-1)),vv((nx-1)*ny*nz+nx*(ny-1)*nz+nx*ny*(nz-1));
+  V1=W1.transpose()*N.col(0);
+  V2=W2.transpose()*N.col(1);
+  V3=W3.transpose()*N.col(2);
+  vv << V1,V2,V3;
 
+  Eigen::SparseMatrix<double> gg;
+  fd_grad(nx,ny,nz,h,gg);
+
+  Eigen::BiCGSTAB<Eigen::SparseMatrix<double>> solver;
+  solver.compute(gg.transpose()*gg);
+  printf("%d %d ",gg.rows(),gg.cols());
+  Eigen::MatrixXd b=gg.transpose()*vv;
+  g=solver.solve(b);
+
+  fd_interpolate(nx,ny,nz,h,corner,P,W);
+  double sigma=(Eigen::MatrixXd::Ones(1,n)*W*g).value()/n;
+  g=g-Eigen::VectorXd::Ones(nx*ny*nz)*sigma;
   ////////////////////////////////////////////////////////////////////////////
   // Add your code here
   ////////////////////////////////////////////////////////////////////////////