[SciPy-user] How to fit a surface from a list of measured 3D points ?

[Zachary Pincus]

>> If not, then you'll have to:
>> a) fit a surface to all of the data in 3D (something done a lot by
>> computer graphics and robotics people, who get point clouds as return
>> data from LIDAR scanners and similar, and then try to fit the points
>> to 3D surfaces for visualization / navigation)
>>
>> b) Find locally-smooth patches and fit surfaces to these individually
>> (the manifold-learning folks do this, e.g. "Hessian LLE"). Say you're
>> interested in curvature around a given data point (x, y, z)... you
>> could take the points within some neighborhood and then either fit
>> them to a simple 3d surface (like some kind of paraboloid), or figure
>> out (with e.g. PCA) the best projection of those data points to a
>> plane, and then fit a surface to f(x, y) -> z for the transformed  
>> data.
>>
>> or perhaps even c) just calculate what you need from the data points
>> directly. If you just need very local curvature data, you could
>> probably calculate that from a point and its nearest neighbors. (This
>> is really just a degenerate case of b...)
>>
>> Lots of tools for these tasks are in scipy, but nothing off-the-shelf
>> that I know if.
>
> The method c) seems the simplest at first sight but I see two issues
> for this local approach :
>  - the measured data are noisy. Using the nearest neighbor could give
> a noisy result two, especially when looking at a radius of curvature
>  - I don't see how to use this approach to plot the variation of
> radius of curvature along the surface, It can give me an array of
> radius of curvature, but as my data are not regularly spaced, it won't
> be easy to handle.
>
> Th method b) seems very fuzzy to me : I don't have any knowledge in
> manifold-learning and I would have the second issue of the method c)
> too.
>
> The method a) is what I had initially in mind, but I didn't see how to
> do this in scipy :-(
> I believed at first that I could make a sort of parametric bispline
> fit with the functions available in scipy.interpolate, but I didn't
> succeed in.

Right. A bivariate spline is to fit data of the form f(x, y) -> z,  
whereas your data are just scattered (x, y, z) points. That's not  
going to work unless you can reliably transform them to a function in  
2D. (Perhaps a PCA step to find two axes of maximum variance in the  
data, and then a projection down to these axes will be a general  
method to transform the data to something a bi-spline can fit. On the  
other hand, if the data are, e.g., helical, there would be no  
transform you could do to render them as something that could be fit  
by a function f(x, y) -> z...)

So you'll have to choose some parameterized family of surfaces that  
you feel will fit the data sufficiently, and then find the best-fit  
parameters to give a curve that approximates the data. This is a very  
general problem (nonlinear regression, or just function optimization,  
depending on how you cast it) for which there are many tools in scipy.  
But for figuring out the right family of surfaces to use, you'll have  
to dig into the graphics or robotics literature. Quite likely there's  
some generalized spline surface you could use, but there's nothing of  
the shelf in scipy to do that, other than the basic numerical tools.

Zach