Fitting an ellipsoid to 3D data points

Question

I have a large set of 3D data points to which I want to fit to an ellipsoid.

My maths is pretty poor, so I'm having trouble implementing the least squares method without any math libraries.

Does anyone know of or have a piece of code that can fit an ellipsoid to data which I can plug straight into my project? In C would be best, but it should be no problem for me to convert from C++, Java, C#, python etc.

EDIT: Just being able to find the centre would be a huge help too. Note that the points aren't evenly spaced so taking the mean won't result in the centre.

Answer 1

Least Squares data fitting is probably a good methodology give the nature of the data you describe. The GNU Scientific Library contains linear and non-linear least squares data fitting routines. In your case, you may be able to transform your data into a linear space and use linear least-squares, but that would depend on your actual use case. Otherwise, you'll need to use non-linear methods.

Answer 2

here you go:

This paper describes fitting an ellipsoid to multiple dimensions AS WELL AS finding the center for the ellipois. Hope this helps,

http://www.physics.smu.edu/~scalise/SMUpreprints/SMU-HEP-10-14.pdf

(btw, I'm assuming this answer is a bit late, but I figured I would add this solution for anyone who stumbles across your question in search for the same thing :)

Answer 3

If you want the minimum-volume enclosing ellipsoid, check out this SO answer for a bounding ellipsoid.

If you want the best fitting ellipse in a least-squares sense, check out this MATLAB code for error ellipsoids where you find the covariance matrix of your mean-shifted 3D points and use that to construct the ellipsoid.

Answer 4

I could not find a good Java based algorithm for fitting an ellipsoid, so I ended up writing it myself. There were some good algorithms for an ellipse with 2D points, but not for an ellipsoid with 3D points. I experimented with a few different MATLAB scripts and eventually settled on Yury Petrov's Ellipsoid Fit. It fits an ellipsoid to the polynomial Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz = 1. It doesn't use any constraints to force an ellipsoid, so you have to have a fairly large number of points to prevent a random quardic from being fit instead of the ellipsoid. Other than that, it works really well. I wrote a small Java library using Apache Commons Math that implements Yury Petrov's script in Java. The GIT repository can be found at https://github.com/BokiSoft/EllipsoidFit.

Answer 5

We developed a set of Matlab and Java codes to fit ellipsoids here: https://github.com/pierre-weiss

You can also check our open-source Icy plugin. The following tutorial can be helpful: https://www.youtube.com/endscreen?video_referrer=watch&v=nXnPOG_YCxw

Note: most of the existing codes fit a generic quadric and do not impose an ellipsoidal shape. To get more robustness, you need to go to convex programming rather than just linear algebra. This is what is done in the indicated sources.

Cheers, Pierre

Answer 6

I have an idea. Approximately solution, not the best but will keep points inside. In XY plane find the radius R1 that will obtain all points. Same do for the XZ plane (R2) and YZ plane (R3). Then use the maximums on each axes. A=max(R1,R2), B=max(R1,R3) and C=max(R2,R3). But, first of all find the average (center) of all points and align it to origin.

Answer 7

I have just gone through the same process. Here is a python module which is based on work by Nima Moshtagh. Referenced in many places but also in this question about a Bounding ellipse

This module also handles plotting of the final ellipsoid. Enjoy!

https://github.com/minillinim/ellipsoid/blob/master/ellipsoid.py

Answer 8

I ported Yury Petrov's least-squares Matlab fitter to Java some time ago, it only needs JAMA: https://github.com/mdoube/BoneJ/blob/master/src/org/doube/geometry/FitEllipsoid.java

Answer 9

Here is very simple method to find foci's, based on random search. No linear algebra used. So if you are fine with non-perfect solution:

3DPoints - Array of some Amount of points
vecCenter - average of 3DPoints
AngleCosine - cos of angle between two vectors

RandomOrder(3DPoints)

vecFocus := (0, 0, 0)
for i := 0 to Amount:
    vecRadius := 3DPoints[i] - vecCenter

    // Change vecRadius direction to parallel
    if AngleCosine(vecFocus, vecRadius) < 0 then
        vecRadius *= -1
    vecFocus += (vecRadius - vecFocus) / Amount

Result can be improved by passing array 2nd time in backward direction Foci's placed at vecCenter +/- vecFocus