I'm using a Truncated Signed Distance Function to perform 3D reconstruction from depth images (http://graphics.stanford.edu/papers/volrange/volrange.pdfTruncated Signed Distance Function to perform 3D reconstruction) from depth images.
Essentially I have a large voxel grid where each voxel contains the signed distance to an implicit surface.
The The surface can be extracted form the Voxel Grid using a marching cubes algorithm (See Wikipedia: Marching_Cubes, Apologies but I can't publish more than 2 links)marching cubes algorithm to generate a mesh. The vertices of this generated mesh lie on the edges of the cubes formed by joining the central points of 8 adjacent voxels.
v0, v1 … vm are$v_0, v_1, \cdots, v_m$ are points in 3D space. They initially lie on a regular lattice with width w$w$, height h$h$ and depth d$d$ so that m = whd$m = whd$.
Other than at the boundaries of the grid, each vx$v_x$ is connected to 6 neighbouring points:
vx-1, vx+1, vx-wh, vx+wh, vx+w and vx-w$v_{x-1}, v_{x+1}, v_{x-wh}, v_{x+wh}, v_{x+h}$ and $v_{x-w}$.
The set of points pa,b$p_{a,b}$ lie on the edges connecting pairs of vertices va and$v_a$ and vb$v_b$. I have at most one point per edge that I care about but for some edges, I have no point. For example, in the image, the edges connecting vwh$v_{wh}$ to vwh+1$v_{wh+1}$ and vwh$v_{wh}$ to v2wh$v_{2wh}$ have no point defined.
For each defined point, we can define alpha$\alpha$ (green in the image), which is the proportion of the way along the edge from the first vertex to the second. So for example, p0,1$p_{0,1}$ has alpha$\alpha$ equal to
|v0p0,1| / |v0v1|$$\alpha = \frac{v_0 p_{0,1}}{v_0 v_1}$$
At time 1, I require each of the each of the points pa,b$p_{a,b}$ to be in a new location, p’a,b$p'_{a,b}$ (known) which requires each of these points to undergo a translation. The overall transformation of the point cloud is not rigid, i.e. each point will be translated by a different amount.
I can only translate the points pa,b$p_{a,b}$ by moving the vertices vx$v_x$. When I do this, the distances between neighbouring vertices can change arbitrarily however the values of alpha$\alpha$ computed after remain constant, i.e. the point pa,b$p_{a,b}$ still divides the edge vavb$v_{a} v_{b}$ in the ratio alpha:(1-alpha)$\alpha:(1 - \alpha)$.
I want to work out the ‘best' translations for v0, … vm which$v_0, \cdots, v_m$ which will :
- Get each of the points p$p$ as close as possible to its desired position; and
- Minimise the differences between the translations applied to neighbouring v’s$v$s
The data term, Edata$E_\text{data}$ minimises the squared difference between the desired coordinates of p’a,b$p'_{a,b}$ and the computed values of p’a,b$p'_{a,b}$ which are given by v’a + alpha(v'b -v’a)$v_a + \alpha v_b$
Edata = Σab ( v’a + alphaab(v'b - v’a))$$E_\text{data} = \sum_{a, b} (v'_a + \alpha_{ab}(v'_b - v'_a) - p_{a,b})^2$$
The sum is across all edges and so each -$v'_a$ and p’a,b$v
_b$ may be involved in multiple terms.
The regularisation term )2
The sum is across all edges and so each v’a and v’b may be involved in multiple terms.
The regularisation term Ereg will aim to minimise the differences between vertex transforms and will be
Ereg = Σi=0..m Σj in N(i) ( | ti | - | tj | )2
$E_\text{reg}$ will aim to minimise the differences between vertex transforms and will be$$E_\text{reg} = \sum_{i=0}^m \sum_{j\in N(i)} (|t_i| - |t_j|)^2$$
Where N(i)$N(i)$ is the neighbourhood of i and comprises all of the adjacent vertices.
Note that there are more regularisation terms (m$m$) than there are data terms (because there’s not one p$p$ per v$v$).
I’m thinking of using non-linear least squares to solve this. I’ve written code which does this for a simpler problem in MatlabMATLAB and it works fine. The simpler problem uses :
- Only 2 dimensions
- A very small grid
- The built in MatlabMATLAB lsqnonlin function.
- It seems like using Newton’s method should work however to compute the gradient vector Hessian matrix it feels like I need the same number of data and regularisation terms so I need to express Edata$E_\text{data}$ in terms of v0…vm but$v_0, \cdots, v_m$ but am unsure how to do this. Is this required ?
- Given that m$m$ is of the order of tens, if not hundreds of thousands, computing the Hessian and it’s inverse seems like it will take a very long time. The Hessian should be sparse but is there a better way to solve this ?
- It feels like avoiding constructing matrices altogether would be beneficial but is there a simple (to understand) way of achieving this ?
- I feel that using CUDA to parallelize the problem would help with performance but again, I'm not sure how I would go about crafting the algorithm to use.
Thanks
References
- Curless, Brian, and Marc Levoy. "A volumetric method for building complex models from range images." Proceedings of the 23rd annual conference on Computer graphics and interactive techniques. ACM, 1996.
