Abstract

This code aims for ellipse fitting but only contains the center’s calculation.

CODE

function [x0,y0,theta,ap,bp] = fitellipseDirect(x,y)
    n = size(x,1);
    centerx = mean(x);
    centery = mean(y);
    A = zeros(n,6);
    scale =  double(0.01);
    for i = 1:n
        px =  double((x(i) - centerx)*scale);
        py =  double((y(i) - centery)*scale);
        A(i,1) =  double(px^2);
        A(i,2) =  double(px*py);
        A(i,3) =  double(py^2);
        A(i,4) =  double(px);
        A(i,5) =  double(py);
        A(i,6) = 1.0;
    end

    DM =  double(A' * A);
    DM =  double(DM/n);

    S1 =  DM(1:3,1:3);
    S2 =  DM(1:3,4:6);
    S2T =  DM(4:6,1:3);
    S3 =  DM(4:6,4:6);

    C1_inv = [0,   0,    0.5;
          0,   -1,     0;
          0.5, 0.0,  0.0];
    M =  double(C1_inv*(S1 - S2*inv(S3)*S2T));
    [v,D] = eig(M);
    index = 1;
    if D(1,1)>0
        index = 1;
    elseif D(2,2)>0
        index = 2;
    elseif D(3,3)>0
        index = 3;
    end

    vec =  [v(1,index);v(2,index);v(3,index)];
    theta_1 =  double(vec);
    theta_2 =  double(-inv(S3)*S2T*theta_1);

    %% method 1, only calculate the ellipse's center
    % b2_4ac =  theta_1(2)*theta_1(2) - 4*theta_1(1)*theta_1(3);
    % cd_be =  2*theta_1(3)*theta_2(1) - theta_1(2)*theta_2(2);
    % ae_bd =  2*theta_1(1)*theta_2(2) - theta_1(2)*theta_2(1);
	% x0 =  (cd_be/b2_4ac/scale + centerx);
	% y0 =  (ae_bd/b2_4ac/scale + centery); 

    %% method 2, calculate the ellipse's parameters
    % assumes a cartesian form ax^2 + 2bxy + cy^2 + 2dx + 2fy + g = 0.
    a = theta_1(1);
    b = theta_1(2) / 2;
    c = theta_1(3);
    d = theta_2(1) / 2;
    f = theta_2(2) / 2;
    g = theta_2(3);
    den = (b^2 - a*c);
    x0 = (c*d - b*f) / den /scale + centerx;
    y0 = (a*f - b*d) / den /scale + centery;
    
    num = 2 * (a*f^2 + c*d^2 + g*b^2 - 2*b*d*f - a*c*g);
    fac = sqrt((a - c)^2 + 4*b^2);
    ap = sqrt(num / den / (fac - a - c))/scale;
    bp = sqrt(num / den / (-fac - a - c))/scale;

    theta = 0;
    if (theta_1(2)  == 0) 
        if (theta_1(1)  < theta_1(3) )
            theta = 0;
        else 
            theta = pi/2.;
        end
    else 
        theta = pi/2. + 0.5*atan2(theta_1(2) , (theta_1(1)  - theta_1(3) ));
    end

    %% obtain the angle of bounding rectangle
    % angle = 0;
    % if( ap > bp )
    %     angle = (mod((90 + theta*180/pi),180.0)) ;
    % else 
    %     angle = (mod(theta*180/pi,180.0));
    % end

    %% ellipse formulation
    % x=a*cost*cosθ-b*sint*sinθ+X,
    % y=a*cost*sinθ+b*sint*cosθ+Y.

    x_ = [];
    y_ = [];
    for t = 0:0.02:2*pi
        x_ = [x_;ap*cos(t)*cos(theta)-bp*sin(t)*sin(theta)+x0];
        y_ = [y_;ap*cos(t)*sin(theta)+bp*sin(t)*cos(theta)+y0];
    end

    %% compare the results
    plot(x,y,'ro');
    hold on;
    plot(x_,y_,'b-');
end

使用：

[x0,y0,theta,ap,bp]  = fitellipseDirect(xcell,ycell);

xcell为x坐标数组，ycell为y坐标数组；

输出结果：