Introduction to Problem ----------------------- I'm using a Truncated Signed Distance Function to perform 3D reconstruction from depth images (http://graphics.stanford.edu/papers/volrange/volrange.pdf). Essentially I have a large voxel grid where each voxel contains the signed distance to an implicit surface. 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) 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. Now the volume that I'm reconstructing is moving in a non-rigid way. I am able to track the motion of points on the surface of the object. In this way I can determine how a particular point on the mesh moved between one frame and another. I want to use this information to generate the effective transformation (translation) of each voxel in the voxel grid which would lead to the deformation of the implicit mesh which most closely matches the observed scene. I have the situation described below: [![https://i.sstatic.net/MtAYt.png][1]][1] *v<sub>0</sub>, v<sub>1</sub> … v<sub>m</sub>* are points in 3D space. They initially lie on a regular lattice with width *w*, height *h* and depth *d* so that *m* = *whd*. Other than at the boundaries of the grid, each *v<sub>x</sub>* is connected to 6 neighbouring points: *v<sub>x-1</sub>, v<sub>x+1</sub>, v<sub>x-wh</sub>, v<sub>x+wh</sub>, v<sub>x+w</sub> and v<sub>x-w</sub>* The set of points *p<sub>a,b</sub>* lie on the edges connecting pairs of vertices *v<sub>a</sub>* and *v<sub>b</sub>*. 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 *v<sub>wh</sub>* to *v<sub>wh+1</sub>* and *v<sub>wh</sub>* to *v<sub>2wh</sub>* have no point defined. For each defined point, we can define *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, *p<sub>0,1</sub>* has *alpha* equal to |*v<sub>0</sub>p<sub>0,1</sub>*| / |*v<sub>0</sub>v<sub>1</sub>*| The problem ----------- The grid as shown represents the situation at time 0. At time 1, I require each of the each of the points *p<sub>a,b</sub>* to be in a new location, *p’<sub>a,b</sub>* (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 *p<sub>a,b</sub>* by moving the vertices *v<sub>x</sub>*. When I do this, the distances between neighbouring vertices can change arbitrarily however the values of *alpha* computed after remain constant, i.e. the point *p<sub>a,b</sub>* still divides the edge *v<sub>a</sub>v<sub>b</sub>* in the ratio *alpha*:(*1-alpha*). I want to work out the ‘best' translations for *v<sub>0</sub>, … v<sub>m</sub>* which will : * Get each of the points *p* as close as possible to its desired position; and - Minimise the differences between the translations applied to neighbouring *v*’s What I think is true --------------------- I need to construct an energy equation comprising a data term and a regularisation term and then minimise this equation. The data term, **E<sub>data</sub>** minimises the squared difference between the desired coordinates of *p’<sub>a,b</sub>* and the computed values of *p’<sub>a,b</sub>* which are given by *v’<sub>a</sub> + alpha(v'<sub>b</sub> -v’<sub>a</sub>)* **E<sub>data</sub>** = Σ<sub>*ab*</sub> ( *v’<sub>a</sub>* + *alpha<sub>ab</sub>*(*v'<sub>b</sub>* - *v’<sub>a</sub>*)) - *p’<sub>a,b</sub>* )<sup>2</sub> The sum is across all edges and so each *v’<sub>a</sub>* and *v’<sub>b</sub>* may be involved in multiple terms. The regularisation term **E<sub>reg</sub>** will aim to minimise the differences between vertex transforms and will be **E<sub>reg</sub>** = Σ<sub>*i=0..m*</sub> Σ<sub>*j in N(i)*</sub> ( | *t<sub>i</sub>* | - | *t<sub>j</sub>* | )<sup>2</sup> Where *N(i)* is the neighbourhood of i and comprises all of the adjacent vertices. **Note that there are more regularisation terms (*m*) than there are data terms (because there’s not one *p* per *v*).** My questions ------------- I’m thinking of using non-linear least squares to solve this. I’ve written code which does this for a simpler problem in Matlab and it works fine. The simpler problem uses : - Only 2 dimensions - A very small grid - The built in Matlab <code>lsqnonlin<code> function. Now I want to build my own code which can solve this for the full scale problem in 3D but I have some concerns: - 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 **E<sub>data</sub>** in terms of *v<sub>0</sub>…v<sub>m</sub>* but am unsure how to do this. Is this required ? - Given that *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. Can anyone point me at appropriate literature or suggest a good starting point for research ? Thanks [1]: https://i.sstatic.net/MtAYt.png