Ellipse fitting

taby said:
Any glaring errors in logic that I'm missing?

I guess your code assumes the ellipse to be axis aligned?
If you would rotate your input points 45 degrees, would you get the same output but rotated 45 degrees as well?

Here's what i would do:
Calculate the average of all points and make it the center of the ellipse.
Calculate the line best fitting the point distribution and make it the longer axis of the ellipse.
The shorter axis is simply the perpendicular line.
Knowing both axis and thus the orientation of the ellipse, we can now measure the radi by summing up point distribution along those axis.

In 2D i would use complex numbers raised to the power of two to fit the line.
In 3D i would use covariance matrix and singular value decomposition to find orthogonal axis.

However, my proposal assumes your points are distributed with uniform density across the whole surface (or volume) of the ellipse / ellipsoid.
If your input covers only a small path of a larger ellipse for example, my proposals would fail and you need something else.

I think i saw related tools here: https://github.com/HRI-EU/WildMagic5p17/tree/master/LibMathematics/Approximation

I got it to work. Yes it takes rotation into account. I answered the question on Physics StackExchange.

P.S. Grok, Claude, and ChatGPT all had errors in the code. AI is shit.

Arbitrary rotation works… sort of ok.

taby said:
P.S. Grok, Claude, and ChatGPT all had errors in the code. AI is shit.

I still have not used such tools and can't tell.
I see lots of ‘X, but made in the 50's’ video suggestions on YT. Quite impressive. Once this can generate assets for games, well… Nintendo will change their mind.
And they are very close:

(Reminds me a bit on my own work - i knew remeshing volume data would work well for AI)

I also took the bait with AI music. My favorite:

I did not know this is AI. I started to speculate pretty quickly, but anyway: I liked it.
Then i found there is countless of similar stuff, and yes - it's all AI. Imaginary Band names, years, and album covers.

I've created thousands of album covers i guess. And i can tell you: I have seen covers much worse than those, including my own.

And i know enough about music to judge that too. It's low quality, imperfect, sometimes just wrong.
But it's the best ‘progressive rock’ album i've heard in years. The compositions are creative and interesting, and it sounds fresh and like something new.

So, you don't get away with just saying ‘AI is shit’. It's not. It has obvious skills. Maybe those skills are never exactly what the ‘creator’ wants, but without doubt it shows creativity and talent.
It can extract those skills from the training data. It sees other patterns than humans do.

And then it spits something out which can be useful, e.g. to serve the growing content addiction of a future generation of humans, which then are all ‘creators’. They can no longer play the guitar or paint an image themselves, but they can express their desire with words, and the machine makes it real for them, ready to be consumed.

However, to keep the ball rolling, the machinery needs food. It needs new data. It can not recycle itself, because then quality decreases, and the consumers are no longer happy. They might revolt against the system. This shall not happen.
To avoid the collapse, the system rewards an elite of human artists. They become idols, super rich, and they feed their creativity directly into the machine. They have mechanical hands with 20 fingers to play guitar faster than Yngwie, and shiny headsets to extract all their thoughts.

Hehehe, cyberpunk is now. And nothing really changes. : )

I am being half-joking when I call AI a bad name. It's impressive how well it works in general. However, it is also impressive how well the different models all stole the same shit code that didn't work at all. LOL

JoeJ – I've had a change of heart. I'm in love with the AIs, now that I've spent a week using them to fool around with ellipse fitting. I've tried Copilot, ChatGPT, Grok, and Claude – my favourite.

There's a lot of example code that solves for the linear ellipse coefficients, using Jacobi SVD for example. I managed to stumble upon a different method, using gradient descent, but I'm still having small problems though.

It works, almost. The main problem is that it erroneously places the orbiting body at one of the foci.

Any glaring problems with the code? It's gotta be something simple, but a solution eludes me for now.

struct Point {
	double x, y, vx, vy;
};

// Helper function for distance calculation
double distance(double x1, double y1, double x2, double y2) {
    return sqrt(pow(x2 - x1, 2) + pow(y2 - y1, 2));
}

// Objective function for optimization (least squares)
double objectiveFunction(const VectorXd& params, const vector<Point>& points, const Point& focus) {
    double h = params[0], k = params[1], a = params[2], b = params[3];
    double error = 0;
  
	for (const auto& p : points)
	{
        double dist1 = distance(p.x, p.y, h, k);
        double dist2 = distance(p.x, p.y, focus.x, focus.y);
        error += pow(dist1 + dist2 - 2 * a, 2); // Distance condition
        
        // Since we're axis-aligned, we simplify velocity condition:
        // Velocity should be more in line with the axis of the ellipse
        double velError = 0;

        if (abs(p.vx) > abs(p.vy)) 
		{ // Suggesting a is along x
			velError = pow(p.vy / p.vx - (k - p.y) / (h - p.x), 2); // Check alignment with y
        }
		else
		{
			velError = pow(p.vx / p.vy - (h - p.x) / (k - p.y), 2); // Check alignment with x
        }

        error += velError;
    }

    return error;
}

// Simple solver function using gradient descent (for demonstration)
VectorXd solveEllipseParameters(const vector<Point>& points, const Point& focus) 
{
    VectorXd params(4); // h, k, a, b

	vector<double> mvec;
	mvec.push_back(max(abs(points[0].x), abs(points[0].y)));
	mvec.push_back(max(abs(points[1].x), abs(points[1].y)));
	mvec.push_back(max(abs(points[2].x), abs(points[2].y)));
	mvec.push_back(max(abs(points[3].x), abs(points[3].y)));
	mvec.push_back(max(abs(points[4].x), abs(points[4].y)));

	sort(mvec.begin(), mvec.end());

	double m = mvec[4];
	
    params << m*0.5, m * 0.5, m * 0.5, m * 0.5; // Initial guess

    int iterations = 100000;
    double stepSize = 0.00001;

    for (int i = 0; i < iterations; ++i) {
        VectorXd gradient = VectorXd::Zero(4);
        for (int j = 0; j < 4; ++j) {
            VectorXd paramsPlus = params;
            paramsPlus[j] += stepSize;
            VectorXd paramsMinus = params;
            paramsMinus[j] -= stepSize;
            
            gradient[j] = (objectiveFunction(paramsPlus, points, focus) - objectiveFunction(paramsMinus, points, focus)) / (2 * stepSize);
        }
        params -= stepSize * gradient;
    }

    return params;
}

…

vector<Point> points;
	
Point point0(orbit_points[0].x, orbit_points[0].y, orbit_velocities[0].x, orbit_velocities[0].y);
Point point1(orbit_points[1].x, orbit_points[1].y, orbit_velocities[1].x, orbit_velocities[1].y);
Point point2(orbit_points[2].x, orbit_points[2].y, orbit_velocities[2].x, orbit_velocities[2].y);
Point point3(orbit_points[3].x, orbit_points[3].y, orbit_velocities[3].x, orbit_velocities[3].y);
Point point4(orbit_points[4].x, orbit_points[4].y, orbit_velocities[4].x, orbit_velocities[4].y);
		
points.push_back(point0);
points.push_back(point1);
points.push_back(point2);
points.push_back(point3);
points.push_back(point4);

Point focus = { 0, 0 };

VectorXd params = solveEllipseParameters(points, focus);

double h = params[0], k = params[1], a = params[2], b = params[3];

global_ep.angle = 0;
global_ep.centerX = 0;
global_ep.centerY = k;
global_ep.semiMajor = a;
global_ep.semiMinor = b;

cout << global_ep.angle << endl;
cout << global_ep.centerX << endl;
cout << global_ep.centerY << endl;
cout << global_ep.semiMajor << endl;
cout << global_ep.semiMinor << endl;

taby said:
It works, almost. The main problem is that it erroneously places the orbiting body at one of the foci.

In the picture the green points are either a line, or a small segment of a large ellipse.
For such data i would hope Dave Eberlys codebase from the link has a fitting function which works.

But in the first picture you have posted the points are more like a dsitribution of density, and my proposal of using SVD would match an ellpsoid to represent this distribution, similar to the inertia tensor used in physics simulation.

My SVD idea is somewhat different than the one you use with Eigen library. I see you give it per point data wich reminds on leat squares methods, but i don't need such data. In 3D space, i would accumulate 3x3 covariance matrix from all points, then doing SVD on that.
However, the result of my method would look very similar to what you show in the first 3 images.

I don't know much about curve fitting, but surely you would need to pick different methods, depending on what your given points typically look like.

A picture of what i mean:

Left is expected result from density SVD, right is curve fitting.
I guess you need to decide for the better method depending on your data.

I found a solution. I was using the ellipse centre instead of the 2nd focus when calculating the distance equation. So, I calculate the 2nd focus now, and I use it.

   struct EllipseParameters2 {
       double semi_major_axis;
       double eccentricity;
       double angle;
       Eigen::Vector2d center;
   };

   // Helper function for distance calculation
   double distance(double x1, double y1, double x2, double y2) 
   {
       return sqrt(pow(x2 - x1, 2) + pow(y2 - y1, 2));
   }
   
   // Objective function for optimization (least squares)
   double objectiveFunction(
       const VectorXd& params, 
       const vector<cartesian_point>& points, 
       const vector<cartesian_point>& velocities, 
       const cartesian_point& focus) 
   {
       double h = params[0], k = params[1], a = params[2], b = params[3];
       double error = 0;
   
       EllipseParameters2 ep;
       ep.angle = 0; 
       ep.center(0) = 0;
       ep.center(1) = k;
       ep.semi_major_axis = a;
       ep.eccentricity = sqrt(1 - (b * b) / (a * a));
   
       Eigen::Vector2d focus1;
       focus1(0) = focus.x;
       focus1(1) = focus.y;
   
       Eigen::Vector2d focus2 = calculateSecondFocus(ep, focus1);
    
       for(size_t i = 0; i < points.size(); i++)
       {
           cartesian_point p = points[i];
           cartesian_point v = velocities[i];
   
           double dist1 = distance(p.x, p.y,  focus1(0),  focus1(1));
           double dist2 = distance(p.x, p.y,  focus2(0),  focus2(1));
   
           error += pow(dist1 + dist2 - 2 * a, 2);
   
           // Since we're axis-aligned, we simplify velocity condition:
           // Velocity should be more in line with the axis of the ellipse
           double velError = 0;
   
           if (abs(v.x) > abs(v.y)) // Suggesting a is along x
               velError = pow(v.y / v.x - (k - p.y) / (h - p.x), 2); // Check alignment with y
           else
               velError = pow(v.x / v.y - (h - p.x) / (k - p.y), 2); // Check alignment with x
   
           error += velError;
       }
   
       return error;
   }

Also, I tweaked how the gradient descent initial parameters are set.

   // Simple solver function using gradient descent (for demonstration)
   VectorXd solveEllipseParameters(const vector<cartesian_point>& points, const vector<cartesian_point>& velocities, const cartesian_point& focus)
   {
       // Get max distance data
       vector<double> mvec;
   
       mvec.push_back(points[0].length());
       mvec.push_back(points[1].length());
       mvec.push_back(points[2].length());
       mvec.push_back(points[3].length());
       mvec.push_back(points[4].length());
   
       sort(mvec.begin(), mvec.end());
   
       // Use the maximum distance data
       const double m = mvec[4];
   
       const double d = (mvec[4] - mvec[0]) / mvec[4];
       
       VectorXd params(4); // h, k, a, b
   
       cout << "d: " << d << endl;
   
       if(d < 0.1)
       params << 1, 1, m, m; // Initial guess
       else
       params << 1, 1, m*0.125, m*0.125; // Initial guess
   
   
       int iterations = 100000;
       double stepSize = 0.0001;
   
       for (int i = 0; i < iterations; i++) 
       {
           VectorXd gradient = VectorXd::Zero(4);
   
           for (int j = 0; j < 4; j++) 
           {
               VectorXd paramsPlus = params;
               paramsPlus[j] += stepSize;
               VectorXd paramsMinus = params;
               paramsMinus[j] -= stepSize;
               
               gradient[j] = (objectiveFunction(paramsPlus, points, velocities, focus) - objectiveFunction(paramsMinus, points, velocities, focus)) / (2 * stepSize);
           }
   
           params -= stepSize * gradient;
       }
   
       return params;
   }

The results are not quite perfect, but they're closer. Any suggestions on how to better go about calculating the initial parameters for the gradient descent are very welcome.

I fixed the problematic setting of initial guesses for the gradient descent. Works great!

VectorXd solveEllipseParameters(const vector<cartesian_point>& points, const vector<cartesian_point>& velocities, const cartesian_point& focus)
{
	// Get max distance data
	vector<double> mvec;

	mvec.push_back(points[0].length());
	mvec.push_back(points[1].length());
	mvec.push_back(points[2].length());
	mvec.push_back(points[3].length());
	mvec.push_back(points[4].length());

	sort(mvec.begin(), mvec.end());

	// Use the maximum distance data	
	const double m = mvec[4];
		
	double d = (mvec[4] - mvec[0]) / mvec[4];
	d = pow(1 - d, 20.0);

	VectorXd params(4); // h, k, a, b
	params << 1, 1, m * d, m * d; // Initial guess

	int iterations = 1000;
	double stepSize = 0.0001;

	for (int i = 0; i < iterations; i++)
	{
		VectorXd gradient = VectorXd::Zero(4);

		for (int j = 0; j < 4; j++)
		{
			VectorXd paramsPlus = params;
			paramsPlus[j] += stepSize;
			VectorXd paramsMinus = params;
			paramsMinus[j] -= stepSize;

			gradient[j] = (objectiveFunction(paramsPlus, points, velocities, focus) - objectiveFunction(paramsMinus, points, velocities, focus)) / (2 * stepSize);
		}

		params -= stepSize * gradient;
	}

	return params;
}

I’ve been using ChatGPT, Grok, Copilot, and Claude. I carefully studied their output, when it came to solving for the ellipse parameters using both gradient descent, and Jacobi SVD. The ellipse that it solves for is not really similar to the ellipse gotten via numerical integration of position and velocity using 4th-order symplectic solver.

The AI is amazing in so many ways, but in the end, it fails.

taby said:
The AI is amazing in so many ways, but in the end, it fails.

Not anymore.

You did post your solution on githib, no? ; )

I did. LOL There's a lot of good code in there though, like converting an angle-time coordinate into Cartesian coordinates and velocity → https://github.com/sjhalayka/ellipse_fitting

But in the end, Jacobi SVD and gradient descent are both numerical solutions that don't really work. If we're going to use a numerical solution, then we might as well just track the orbit based on the input position and velocity, using symplectic integration. Gathering the angle, semi-major axis and semi-minor axis, and the centre from the orbit data would be pretty straightforward. If we demand that we use a numerical method, then the orbit tracking is the best bang for the buck.

I finally found an analytical solution to the ellipse fitting problem. Basically, the person who asked the question on Physics Stack Exchange said that they had already tried the AIs to generate the necessary code. Well, I finally found the answer using Meta AI. The instructions that I gave were “how do i use gauss's method with three input points for orbit ellipse center, semi-major axis, and semi minor axis determination. please use c++. use heliocentric position. use equatorial plane for orbit. don't forget v2.”

So I wasted all that time on numerical solutions. :)

The answer and GitHub are both up-to-date:

https://physics.stackexchange.com/a/839919/142234

https://github.com/sjhalayka/ellipse_fitting

A point belonging to an ellipse can be described as a linear combination of 1, x, y, x^2, xy and y^2 being zero. You can fix the coefficient for 1 to 1 and fit the others, using least squares regression. This should be fast and robust.

I tried solving for the ellipse coefficients (A, B, C, D, E, F) and / or parameters (centre, axes, angle) using Jacobi SVD and gradient descent. Both looked nice, but neither was an actual close answer. I now know a lot more about linear algebra, even though the project failed. For what it's worth, when solving Mx = b, I had to set b's elements to -1, not zero.

taby said:
For what it's worth, when solving Ax = b, I had to set b's elements to -1, not zero.

That's definitely why it failed. When setting up the system it should look like this:

A: numPoints * numCoefficients
x: numCoefficients * 1
b: numPoints * 1

A:
[ 1 x0 y0 x0*x0 x0*y0 y0*y0 ]
[ 1 x1 y1 x1*x1 x1*y1 y1*y1 ]
[ .... ]
[ 1 xN yN xN*xN xN*yN yN*yN ]

b:
[ 0 ]
[ ... ]
[ 0 ]

You also will need at least as many points as coefficients, or else the system is underdetermined. Once you set up the system as above, it should just be a simple least squares solve, no gradient descent or iterative algorithms required.

https://en.wikipedia.org/wiki/Ellipse#General_ellipse

No, I’ve made you misunderstand. It worked insomuch that it found an ellipse based on the 5 input points and 5 input velocities and one focus.. It looks close on the screen, but It’s just not close enough. It’s a failure.