I pared it down to only needing 3 points.

EllipseParameters extractEllipseParameters(const Eigen::VectorXd& coefficients)
{
    //A(i, 0) = x * x;
    //A(i, 1) = 0;// x* y;
    //A(i, 2) = x * x;
    //A(i, 3) = 0;// x;
    //A(i, 4) = y;

    double a = coefficients(0);
    double b = 0;
    double c = coefficients(1);
    double d = 0;
    double e = coefficients(2);
    double f = 1;

    // Calculate center
    double centerX = (2 * c * d - b * e) / (b * b - 4 * a * c);
    double centerY = (2 * a * e - b * d) / (b * b - 4 * a * c);

    // Calculate rotation angle
    double theta = 0.5 * atan2(b, (a - c));

    // Calculate semi-axes
    double ct = cos(theta);
    double st = sin(theta);
    double ct2 = ct * ct;
    double st2 = st * st;
    double a2 = a * ct2 + b * ct * st + c * st2;
    double c2 = a * st2 - b * ct * st + c * ct2;

    // Calculate constants
    double term = 2 * (a * centerX * centerX + b * centerX * centerY +
        c * centerY * centerY + d * centerX + e * centerY + f);

    double semiMajor = sqrt(abs(term / (2 * std::min(a2, c2))));
    double semiMinor = sqrt(abs(term / (2 * std::max(a2, c2))));

    if (a2 > c2) 
    {
        cout << "Rotating ellipse" << endl;
        std::swap(semiMajor, semiMinor);
        theta += pi / 2;
    }

    EllipseParameters params;
    params.centerX = centerX;
    params.centerY = centerY;
    params.semiMajor = semiMinor;
    params.semiMinor = semiMajor;
    params.angle = theta;

    return params;
}

EllipseParameters fitEllipse(const std::vector<cartesian_point>& points)
{
    if (points.size() != 3 ) {
        std::cerr << "Error: 3 points are required.\n";
        return EllipseParameters();
    }

    Eigen::MatrixXd A(3, 3);
    Eigen::VectorXd b(3);

    // Fill the matrix A and vector b with the equations from the points
    for (size_t i = 0; i < 3; ++i)
    {
        double x = points[i].x;
        double y = points[i].y;

        A(i, 0) = x * x;
        A(i, 1) = y * y;
        A(i, 2) = y;

        b(i) = -1;
    }

    // Solve for the ellipse parameters
    Eigen::VectorXd ellipseParams = A.jacobiSvd(Eigen::ComputeFullU | Eigen::ComputeFullV).solve(b);


    // Compute center of ellipse
    EllipseParameters ep = extractEllipseParameters(ellipseParams);

    global_ep.angle = ep.angle;
    global_ep.centerX = ep.centerX;
    global_ep.centerY = ep.centerY;
    global_ep.semiMajor = ep.semiMajor;
    global_ep.semiMinor = ep.semiMinor;

    return ep;
}

taby said:
I pared it down to only needing 3 points.

From 3 points we can always form a circle, so there can't be enough information to describe an ellipse with that.

I’m sorry, but that’s untrue. Check out the code:

https://github.com/sjhalayka/ellipse_fitting/blob/850496d6aa5530735e0775ad8df3ccaa936fb823/main.cpp#L139

P.S. all I did was disable rotation and translation on the x axis. That still counts as an ellipse.

taby said:
I’m sorry, but that’s untrue. Check out the code:

The code uses a linear algebra math library. Sadly, that's my beyond background and education.

But it's irrelevant. I'm not wrong. Draw 3 random points on paper. You can always connct them with a circle. Always.

Ofc. it's possible to find infinitely many ellipses which connect the points as well. And your code calculates on such random solution. But that's not what you really want, i guess?

If that were true, then Gauss’ method wouldn’t work either.

We have an equation given in https://en.wikipedia.org/wiki/Ellipse#General_ellipse

The equation is $A x^2 + B x y + C y^2 + D x + E y + F = 0$.

The $x y$ term is rotation, and the $x$ term is translation on the x axis. I cut them out, since I didn't need them. The equation is now $A x^2 + C y^2 + E y + 1 = 0$. We only need 3 points to define the three coefficients of an axis-aligned ellipse. We then fill out a matrix of a system with the system of equations, and then use the solver to arrive at the solution for $x$ in the equation $M x = b$, where $b$ is an array full of $-1$. Does this even make sense? I'm trying. :)

… it’s a special ellipse, not general.

taby said:
The xy term is rotation, and the x term is translation on the x axis. I cut them out, since I didn't need them.

Aha, so that's the precious parts you were hiding from me.
I guess cutting them out means they are zero, which is additional information beside the given 3 points.

Does this even make sense? I'm trying. :)

Idk, but sounds reasonable to me. Though, the -1 feels confusing.

Some years back i have used Eigen too, also for some kind of Mx=b problem. It was about diffusing energy over a mesh so it spreads equally fast in all directions, independent of triangle area and angles.

I as very proud i could get it to work. This Mx=b seems so… important? :)
But sadly i have not really learned anything. I ended up using my own custom iterative solver instead. It was slower for big problems, but i understood how it works. : /

That’s deadly awesome man. Yes, the Jacobi SVD algorithm is not clear to me. If you still have it, please post your solver code. Other than GLSL, I know quite little about linear algebra. I couldn’t pass a university class on the subject, that’s for sure.