I have the situation described below:

v₀, v₁ … v_m 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_x is connected to 6 neighbouring points:

v_x-1, v_x+1, v_x-wh, v_x+wh, v_x+w and v_x-w

The set of points p_a,b lie on the edges connecting pairs of vertices v_a and 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 v_wh to v_wh+1 and v_wh to v_2wh 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_0,1 has alpha equal to

|v₀p_0,1| / |v₀v₁|

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_a,b to be in a new location, 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 p_a,b by moving the vertices v_x. 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_a,b still divides the edge v_av_b in the ratio alpha:(1-alpha).

I want to work out the ‘best' translations for v₀, … v_m 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_data minimises the squared difference between the desired coordinates of p’_a,b and the computed values of p’_a,b which are given by v’_a + alpha(v'_b -v’_a)

E_data = Σ_ab ( v’_a + alpha_ab(v'_b - v’_a)) - p’_a,b )²

The sum is across all edges and so each v’_a and v’_b may be involved in multiple terms.

The regularisation term E_reg will aim to minimise the differences between vertex transforms and will be

E_reg = Σ_i=0..m Σ_{j in N(i)} ( | t_i | - | t_j | )2

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 lsqnonlin 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_data in terms of v₀…v_m 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