minillinim / ellipsoid

Example code for enclosing points in an elipsoid
GNU General Public License v3.0
54 stars 16 forks source link

How to get a axis-parallel ellipsoid? #6

Open OvOtoQAQ opened 2 months ago

OvOtoQAQ commented 2 months ago

Hi,

Thanks for you great work! I wonder if it's possible to specify the angle of fitted ellipsoid? I want to get a axis-parallel ellipsoid/ellipse, i.e. set the angle to zero, could you please give me some advice?

All the best

minillinim commented 2 months ago

Hey there. Thanks for the nice words.

You'll have to forgive me, it's been a while since I wrote this code so I'm a bit rusty with it but I don't think it's going to be a matter of making a small adjustment to this code. The issue is that the optimisation is repeatedly rotating the ellipse and adjusting the radii to minimise the volume. If you want an ellipse that has it's axes parallel to the X, Y, Z, ... axes then this will involve modifying the optimisation code at https://github.com/minillinim/ellipsoid/blob/master/ellipsoid.py#L44-L54 to include a constraint so that it keeps the ellipse aligned and only modifies the radii.

This will result in a larger ellipse than the minimal-unaligned / unconstrained one this code currently outputs.

Unfortunately, I really can't remember how this code works anymore so I personally can't help much more than this. That said, I'd love to see how you solve it if you do work it out.

Luckily for us, we live in the future! You might have done this already, but if you paste my code and your question into something like ChatGPT then it's very likely to give you a good place to start out.

Good luck!

OvOtoQAQ commented 2 months ago

Hi,

Thanks for your kind reply! Modifying the Khachiyan Algorithm is a little too fundamental and hard for me, so I tried another iterative optimization method to get a ellipse for 2d points as follows,

from scipy.optimize import minimize

def ellipse_area(params):
    r_maj, r_min, x_c, y_c = params
    return np.pi * r_maj * r_min

def constraint(params, x, y):
    r_maj, r_min, x_c, y_c = params
    return 1 - ((x - x_c)**2 / r_maj**2 + (y - y_c)**2 / r_min**2)

def quadratic_programming(layer_points, init_a, init_b):
    x = layer_points[:, 0]
    y = layer_points[:, 1]

    r_maj_init = init_a
    r_min_init = init_b
    x_c_init = np.mean(x)
    y_c_init = np.mean(y)
    initial_guess = [r_maj_init, r_min_init, x_c_init, y_c_init]

    cons = {'type': 'ineq', 'fun': lambda params: constraint(params, x, y)}

    result = minimize(ellipse_area, initial_guess, constraints=cons)
    x0, y0, a, b = result.x[2], result.x[3], result.x[0], result.x[1]
    return x0, y0, a, b

Hope it will help:)