Skip to content

Surface Reconstruction

Link to paper

Note that the notes is a simplified version of the paper's implementation. We use fixed unit grid node instead of adaptive grid node. Also, we use the trilinear interpolation instead of more complex weighted methods.

Problem Setup

Given a set of scanned point samples PP on the surface of an object and the estimated point normals NN, we want to find an explicit continuous surface representation, i.e. a triangle mesh.

Voxel-based Implicit Surface

However, if we directly convert the point cloud to a mesh by connecting points together, it's difficult to ensure the certain topological postconditions, i.e. a closed manifold with a small number of holes. Instead we will first convert the point cloud sampling representation into an implicit surface representation. Let the unknown surface be the level set of some function g:R3Rg:\mathbb R^3\rightarrow\mathbb R, i.e.

S={xR3g(x)=c}\partial S = \{x\in\mathbb R^3\mid g(x) = c\}

Then, we discretize an implicit function by defining a regular 3D grid of voxels containing at least the bounding box of SS. At each node in the grid xi,j,kx_{i,j,k} we store the value of the implicit function $g(x_{i,j,k}) and then we can get the value everywhere by trilinear interpolation over the grid cell.

Matching Cubes Algorithm

We want the implicit surface representation so that we can contoured into a triangle mesh via Marching Cubes Algorithm.

Characteristic Function of Solids

Assume that our set of point PP lies on the surface S\partial S. This solid object SS must have some non-trivial volume that we can calculate as the integral of unit density over the solid

SdA=R3xS(x)dA=R3I(xS)dA\int_S dA = \int_{\mathbb R^3}\mathcal x_S(x)dA = \int_{\mathbb R^3}\mathbb I(x\in S)dA

From the Poisson Surface Reconstruction [Kazhdan et al. 2006] , the gradient of a infinitesimally mollified characteristic function

  • points in the direction of the normal near the surface S\partial S
  • is zero everywhere else

Therefore, using points PP and normals NN, we can optimize an implicit function gg over a regular grid, so that g\nabla g meets the two properties above, and gg will be an approximation of the mollified characteristic function.

Poisson Surface Reconstruction

Starting from some unrealistic assumptions (note that these assumptions are not to be made but to make the description more intuitive)

  • We know each of g(xi,j,k)\nabla g(x_{i,j,k})
  • All of input points lies perfectly at grid nodes: pP.xi,j,k=p\forall p\in P. \exists x_{i,j,k} = p

With these assumptions, we simply have

g(xi,j,k):R3R3:=vi,j,k=I(xi,j,kpl=xi,j,k)nl\nabla g(x_{i,j,k}):\mathbb R^3\rightarrow\mathbb R^3:= v_{i,j,k} = \mathbb I(x_{i,j,k}\mid \exists p_l = x_{i,j,k}) \:n_l

So that gg is defined via a simple set of linear equations

Since our system is over determined, this can be turned into an optimization problem.

mingijk12g(xi,j,k)vi,j,k2\min_g \sum_i\sum_j\sum_k \frac12\|g(x_{i,j,k}) - v_{i,j,k}\|^2

where gg is a vector of unknown grid nodes values, gi,j,k=g(xi,j,k)g_{i,j,k} = g(x_{i,j,k})

With the assumptions, we can compute approximations of the x,y,z components of g\nabla g via a matrix multiplication of a "gradient matrix" GG and vector unknown grid value gg, so that the optimization problem becomes

ming12Ggv2=ming12gTGTGggTGTv+vTv\min_g \frac12\|Gg - v\|^2 = \min_g\frac12 g^TG^TGg - g^TG^Tv + v^Tv

This is a quadratic function, so that we can set the gradient being 0

g12gTGTGggTGTv+vTv=GTGgGTv=0\frac{\partial}{\partial g} \frac12 g^TG^TGg - g^TG^Tv + v^Tv = G^TGg - G^Tv = 0

Gradient estimated from regular grid

Since we have a grid, it's straightforward to estimate the gradient as a finite difference. For example, gx(xi1/2,j,k)=1h(gi,j,kgi1,j,k)\frac{\partial g}{\partial x}(x_{i-1/2, j, k}) = \frac1h(g_{i,j,k} - g_{i-1,j,k}).
With this observation, we can construct a partial difference matrix DxR(nx1)nynz×nxnynzD^x \in \mathbb R^{(n_x -1)n_yn_z\times n_xn_yn_z}, where each row Di1/2,j,kxD^x_{i-1/2,j,k} computes the partial derivative at xi1/2,j,kx_{i-1/2,j,k}, i.e.

Di1/2,j,kx(l)={1l=i11l=10otherwiseD^x_{i-1/2,j,k}(l) = \begin{cases}-1 &l=i-1\\1 &l=1\\0&\text{otherwise}\end{cases}

We can do similar things on Dy,DzD^y, D^z and obtain the desired

G=(DxDyDz)R(nx1)nynz+nx(ny1)nz+nxny(nz1)×nxnynzG = \begin{pmatrix}D^x\\D^y\\D^z\end{pmatrix} \in \mathbb R^{(n_x-1)n_yn_z + n_x(n_y-1)n_z + n_xn_y(n_z-1)\times n_xn_yn_z}

Of course, when we implement this, we cannot have indexing as 0.50.5, so we can shift the indexing down by 0.50.5, and do the same for the vector vv.

Estimate v

Note that our normals NN does not lie on the grid node. In order to obtain vv on the ideal grid, we can distribute each n\vec n to its 8 neighboring grid nodes via trilinear interpolation weights, for example

nx=w12,0,0(x1,14,12) v12,0,0x+w32,0,0(x1,14,12) v112,0,0x+w32,1,1(x1,14,12) v112,1,1x.\begin{align*} n_x =\:& w_{ \frac12 ,0,0}( \mathbf{x}_{1,\frac14 ,\frac12 } ) \ v^x_{ \frac12 ,0,0} + \\ &w_{\frac32 ,0,0}( \mathbf{x}_{1,\frac14 ,\frac12 } ) \ v^x_{1\frac12 ,0,0} + \\ &\vdots \\ &w_{\frac32 ,1,1}( \mathbf{x}_{1,\frac14 ,\frac12 } )\ v^x_{1\frac12 ,1,1}. \end{align*}

we can do this for each input normal and assmeble a parse matrix Wxn×(nx1)nynzW^x\in n\times (n_x-1)n_yn_z Then, we can have vx=(Wx)TNxv^x = (W^x)^TN^x. Doing similar for vy,vzv^y, v^z and stack them together, we can obtain vv

Poisson Equation

Consider some discrete energy minimization problem

E(g)=ΩgV2dAE(g) = \int_{\Omega}\|\nabla g - V\|^2 dA

the Poisson's equation states that the minimizers on Ω\Omega will satisfy

g=v\nabla\cdot \nabla g = \nabla\cdot v