Fit ellipsoid to (x,y,z) data

I have a 3D dataset having (x,y,z) coordinates. The x and y values are equally spaced (regular grid). How can I fit an ellipsoid of the form (x-p)^2/a^2 + (y-q)^2/b^2 + (z-r)^2/c^2 , where (p,q,r) are the coordinates of the center of the ellipsoid, and a,b,c are the radii?

[Matt J]

This ellipsoid would have no rotation (all its axes are parallel to the x,y,z axes). Is that an important constraint for you?

[Geetartha Dutta]

No, I don't have any rotation constraint. The only constraint I think I have is the fitted surface should be at most half of an ellipsoid (e.g. in case of a sphere, the points should fit at most a (upper) hemisphere). That's because I'm trying to fit topographical data (like a hill) using an ellipsoid (or at most half of it).

[Matt J]

A Gaussian lobe (which has elliptical cross sections) sounds like it would be better for that. Most hills that I've seen start to slope upward gradually at the base. A hemisphere does not.

[Geetartha Dutta]

Unfortunately I need to fit an ellipsoid (does not have to be a sphere). The fit does not have to be very good.

[Geetartha Dutta]

Attached is my xyz data. I know that there seems to be two modes in the data, so an ellipsoid would not fit the data very well. It would be fine if I get a reasonably good fit near the greater mode, so I could try to fit locally in that region.

[Matt J]

I know that there seems to be two modes in the data

Looks like a lot more than that. I can't tell which is supposed to be the "greater" mode. In any case, if you want a good fit in a particular region, you will have to prune the data to exclude the other regions.

[Geetartha Dutta]

Attached is the pruned data. It would be great if I could get a reasonably good fit to this data.

Answer 1

[Matt J]

I'm finding that a decent fitting strategy is to first fit with a Gaussian, but then use the parameters of the Gaussian to construct an ellipsoid hemisphere. For the Gaussian fitting, I used gaussfitn, which is downloadable from,

https://www.mathworks.com/matlabcentral/fileexchange/69116-gaussfitn

load xyz

[maxval,i]=max(z(:));
mu0=[x(i);y(i)];
D0=min(z(:));

opts={'FunctionTolerance',1e-14, 'OptimalityTolerance',1e-14, 'StepTolerance',1e-14};
G0={D0,maxval-D0,mu0,100*eye(2)};
LB={0,0,[],[]};
UB={D0,maxval,[],[]};

G = gaussfitn([x(:),y(:)],z(:),G0,LB,UB,opts{:});

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

%Disaply surfaces
[Zg,Ze]=getSurf(x,y,G);

surf(x,y,z,'FaceAlpha',0.5,'FaceColor','b');
surface(x,y,Ze,'FaceColor','r'); xlabel X, ylabel Y
legend('Raw Data','Fit')

function [Zg,Ze]=getSurf(x,y,G)

  [D,A,mu,sig]=deal(G{:});

  sz=size(x);

  xy=[x(:),y(:)]'-mu;

  Zg=D+A*exp(-0.5*sum( (sig\xy).*xy,1)); Zg=reshape(Zg,sz); %Gaussian Fit
 
  Ze=D+A*sqrt(1-sum( (sig\xy).*xy)); Ze=reshape(Ze,sz);   %Ellipsoid Fit

end

[Geetartha Dutta]

The fit looks good, but I am trying to figure out how I can get the radii of the fitted ellipsoid. My math is not very strong and I'm still trying to understand what you did.

[Matt J]

You would process the parameters from this line, the output of gaussfitn

[D,A,mu,sig]=deal(G{:});

For x and y:

You would simply do an eigendecomposition of sig. The square root of the eigenvalues will be the ellipsoid radii (in the xy dimensions) and the eigenvectors will give the major and minor axis directions. The ellipsoid center coordinates are mu. See also,

https://en.wikipedia.org/wiki/Ellipsoid#As_a_quadric

For z:

The radius in the z-direction is A and the center coordinate is D.

[Geetartha Dutta]

Perfect! Thanks for your help!

[Geetartha Dutta]

For another dataset (attached file 'xyz_2.mat'), I am not able to plot the ellipsoid fit as the output 'Ze' from 'getSurf' is complex valued. I am not sure why. Does this method fail to fit an ellipsoid to this dataset?

[Matt J]

You should set the complex values to NaN. They correspond to (x,y) outside the footprint of the ellipsoid.

[Geetartha Dutta]

I see, thanks!

Answer 2

[Torsten]

https://uk.mathworks.com/matlabcentral/answers/2038016-how-to-plot-the-best-fitted-ellipse-or-circle?s_tid=srchtitle

https://uk.mathworks.com/matlabcentral/fileexchange/87584-object-oriented-tools-for-fitting-conics-and-quadrics

Answer 3

[Matt J]

Using quadricFit from,

https://uk.mathworks.com/matlabcentral/fileexchange/87584-object-oriented-tools-for-fitting-conics-and-quadrics

%%%%%%%%%%%Fake input data

[X,Y,Z] = sphere;

[X,Y,Z]=deal(1+40*X, 2+20*Y,3+30*Z); %stretch into an ellipsoid

surf(X,Y,Z); axis equal

%%%%%%%%%%%  Do the fit

    XYZ=[X(:),Y(:),Z(:)]';

    [XYZ,T]=quadricFit.homogNorm(XYZ);
    
    X=XYZ(1,:).';
    Y=XYZ(2,:).';
    Z=XYZ(3,:).';
    
    e=+ones(size(X,1),1);
    M= [X.^2, [],   [],   X, ...
              Y.^2, [],   Y, ...
                    Z.^2  Z, ...
                          e];
    
   
    coeffs=quadricFit.mostnull(M);
    ABCDEFGHIJ=zeros(1,10);
    ABCDEFGHIJ([1,4,5,7:10])=coeffs;
    ABCDEFGHIJ=num2cell(ABCDEFGHIJ);
    [A,B,C,D,E,F,G,H,I,J]=deal(ABCDEFGHIJ{:});
    
    
    Q=[A, B, C; %D
       0  E, F; %G
       0  0  H];%I
                %J

    Q=Q/2+Q.'/2;            
                
    W=T.'*[Q,[D;G;I]/2;[D,G,I]/2,J]*T;
    
    Q=W(1:3,1:3);
    x0=-Q\W(1:3,end);
    
    T=eye(4); T(1:3,4)=x0;
    
    W=T.'*W*T; W=-W/W(end);

    rad=sqrt(1./diag(W(1:3,1:3)));

    [a,b,c]=deal(rad(1),rad(2),rad(3))  %ellipsoid radii

a = 40.0000

b = 20

c = 30.0000

    [p,q,r]=deal(x0(1),x0(2),x0(3))   %ellipsoid center coordinates

p = 1.0000

q = 2.0000

r = 3.0000

[Geetartha Dutta]

I tried the above code using my data, and it gives complex values for a and b. I am not sure why.

[Matt J]

Attach your xyz data in a .mat file, so it can be examined.