Moving least squares fitting for point displacements having issues

Question

Explanation of the problem: I have points with (x,y,z) coordinates at two+ distinct times. For convenience, they can be imagined as irregularly spaced points along the surface of an inverted paraboloid. There is some minimal thickness to the paraboloid. The paraboloid changes shape slightly as time proceeds (like a balloon inflating) and when it does so, all of the points move. By substracting the coordinates at time2 - time1, I can get the displacement vectors at each point.

It is important to note (and I suspect this might be the source of the problem) that at the first time point, the x and y coordinates range from 0 to 2000, and the z coordinates are all within a narrower range - say 350 to 450. During the deformation, each point has an x component of displacement, y component, and z component. The x and y components are small (~50 at most), while the z component is the largest (goes up to 400 near the center, much less near the edges).

Using weighted moving least squares at the location of each point, I am trying to fit the components of displacements to a second degree polynomial surface in terms of the original x,y,z coordinates of the point: eg. x component of

displacement = ax^2 + bxy + cx + dy^2.. + hz^2 + iz + j

I use the lsqr function in MATLAB,like so, looping through each point for each time interval:

Ux = displacements{k,1}(:,1);
Cx = lsqr((adjust_B_matrix'*W*adjust_B_matrix),(adjust_B_matrix'*W*Ux),1e-7,10000); 

W is the weight matrix, and adjust_B_matrix is the matrix of all (x,y,z) coordinates at time 1, shifted so that they're all centered around the point at which I'm trying to fit the function.

What is going wrong? It's just not working -- once I have the functions, they're re-centered around the actual coordinates of the points. But once I plot the resulting points (initial pointx + displacementx, initial pointy + displacementy, initial pointz + displacementz) by plugging in the coordinates at time 1 into the now-discovered functions, it just spits out a surface that looks just like the surface at time 1.

What might be going wrong? Things I have tried:

1. It's not an issue with the code itself- I generated 'fake' data using a grid of points and it worked perfectly. The predicted locations were superimposed with the actual coordinates and I was able to get back the function I started with. But in my trial example, I used x,y,z from 0 to 5, evenly spaced.
2. Global fitting works (but I need local fitting...). I tried MATLAB's curve fitting toolbox and just tried to fit one of the displacements to only x and y coordinates, globally. It worked perfectly.
3. I think I shouldn't have a singular matrix issue because I use a large radius (around 75-80) points in the calculations, somewhat dispersed in 3D space.

Suspicions: I think it has to do with the uneven distribution of initial (x,y,z) coordinates, but I don't know why or how to fix the issue, or even what method I can use.

If you read this far, thank you so much. Any advice would be greatly appreciated.

Figure for reference: green = predicted points at time 2. Overlapping mostly with red, the actual coordinates of the points at time 1. blue is the correct coordinates of points at time 2 (this is where the green ones should be if things were working).

image

Updated link for files: http://a.tmp.ninja/eWfkNmFZyTFk.zip Contents - code, sample data (please load the .mat files).

Answer 1

I can't actually access the code you posted, so here's some general suggestions.

It does look like the curve fitting toolbox has tools that do exactly what you are looking for, checkout the bottom of this page: https://www.mathworks.com/help/curvefit/polynomial.html#bt9ykh.

It looks like for whatever your learned function for the displacement is just very small or zero everywhere. I suspect the issue is just a minor typo/error on your part somewhere in your pipeline, possibly translating what you have to work with the fit function will reveal the issue.

This really shouldn't be the issue, but in the future if you had much more unbalanced data you could normalize it all before fitting (x_norm = (x - x_mu)/x_std).

Also, I don't think this is your problem either, but you can check if your matrix is close to singular by checking the condition number using the cord() function. So you could check cond(adjust_B_matrix'Wadjust_B_matrix). Second, If you check the documentation for lsqr there is an option to get a debug return flag, that is worth checking too.