I would like to fit a 2D array by an elliptic function: (x / a)² + (y / b)² = 1 ----> (and so get the a and b)
And then, be able to replot it on my graph. I found many examples on internet, but no one with this simple Cartesian equation. I probably have searched badly ! I think a basic solution for this problem could help many people.
Here is an example of the data:
Sadly, I can not put the values... So let's assume that I have an X,Y arrays defining the coordinates of each of those points.
This can be solved directly using least squares. You can frame this as minimizing the sum of squares of quantity (alpha * x_i^2 + beta * y_i^2 - 1) where alpha is 1/a^2 and beta is 1/b^2. You have all the x_i's in X and the y_i's in Y so you can find the minimizer of ||Ax - b||^2 where A is an Nx2 matrix (i.e. [X^2, Y^2]), x is the column vector [alpha; beta] and b is column vector of all ones.
The following code solves the more general problem for an ellipse of the form Ax^2 + Bxy + Cy^2 + Dx +Ey = 1 though the idea is exactly the same. The print statement gives 0.0776x^2 + 0.0315xy+0.125y^2+0.00457x+0.00314y = 1 and the image of the ellipse generated is also below
import numpy as np import matplotlib.pyplot as plt alpha = 5 beta = 3 N = 500 DIM = 2 np.random.seed(2) # Generate random points on the unit circle by sampling uniform angles theta = np.random.uniform(0, 2*np.pi, (N,1)) eps_noise = 0.2 * np.random.normal(size=[N,1]) circle = np.hstack([np.cos(theta), np.sin(theta)]) # Stretch and rotate circle to an ellipse with random linear tranformation B = np.random.randint(-3, 3, (DIM, DIM)) noisy_ellipse = circle.dot(B) + eps_noise # Extract x coords and y coords of the ellipse as column vectors X = noisy_ellipse[:,0:1] Y = noisy_ellipse[:,1:] # Formulate and solve the least squares problem ||Ax - b ||^2 A = np.hstack([X**2, X * Y, Y**2, X, Y]) b = np.ones_like(X) x = np.linalg.lstsq(A, b)[0].squeeze() # Print the equation of the ellipse in standard form print('The ellipse is given by {0:.3}x^2 + {1:.3}xy+{2:.3}y^2+{3:.3}x+{4:.3}y = 1'.format(x[0], x[1],x[2],x[3],x[4])) # Plot the noisy data plt.scatter(X, Y, label='Data Points') # Plot the original ellipse from which the data was generated phi = np.linspace(0, 2*np.pi, 1000).reshape((1000,1)) c = np.hstack([np.cos(phi), np.sin(phi)]) ground_truth_ellipse = c.dot(B) plt.plot(ground_truth_ellipse[:,0], ground_truth_ellipse[:,1], 'k--', label='Generating Ellipse') # Plot the least squares ellipse x_coord = np.linspace(-5,5,300) y_coord = np.linspace(-5,5,300) X_coord, Y_coord = np.meshgrid(x_coord, y_coord) Z_coord = x[0] * X_coord ** 2 + x[1] * X_coord * Y_coord + x[2] * Y_coord**2 + x[3] * X_coord + x[4] * Y_coord plt.contour(X_coord, Y_coord, Z_coord, levels=[1], colors=('r'), linewidths=2) plt.legend() plt.xlabel('X') plt.ylabel('Y') plt.show() 
Following the suggestion by ErroriSalvo, here is the complete process of fitting an ellipse using the SVD. The arrays x, y are coordinates of the given points, let's say there are N points. Then U, S, V are obtained from the SVD of the centered coordinate array of shape (2, N). So, U is a 2 by 2 orthogonal matrix (rotation), S is a vector of length 2 (singular values), and V, which we do not need, is an N by N orthogonal matrix.
The linear map transforming the unit circle to the ellipse of best fit is
sqrt(2/N) * U * diag(S) where diag(S) is the diagonal matrix with singular values on the diagonal. To see why the factor of sqrt(2/N) is needed, imagine that the points x, y are taken uniformly from the unit circle. Then sum(x**2) + sum(y**2) is N, and so the coordinate matrix consists of two orthogonal rows of length sqrt(N/2), hence its norm (the largest singular value) is sqrt(N/2). We need to bring this down to 1 to have the unit circle.
N = 300 t = np.linspace(0, 2*np.pi, N) x = 5*np.cos(t) + 0.2*np.random.normal(size=N) + 1 y = 4*np.sin(t+0.5) + 0.2*np.random.normal(size=N) plt.plot(x, y, '.') # given points xmean, ymean = x.mean(), y.mean() x -= xmean y -= ymean U, S, V = np.linalg.svd(np.stack((x, y))) tt = np.linspace(0, 2*np.pi, 1000) circle = np.stack((np.cos(tt), np.sin(tt))) # unit circle transform = np.sqrt(2/N) * U.dot(np.diag(S)) # transformation matrix fit = transform.dot(circle) + np.array([[xmean], [ymean]]) plt.plot(fit[0, :], fit[1, :], 'r') plt.show() 
But if you assume that there is no rotation, then np.sqrt(2/N) * S is all you need; these are a and b in the equation of the ellipse.
You could try a Singular Value Decomposition of the data matrix. https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.linalg.svd.html
First center the data by subtracting mean values of X,Y from each column respectively.
X=X-np.mean(X) Y=Y-np.mean(Y) D=np.vstack(X,Y) Then, apply SVD and extract -eigenvalues (members of s) -> axis length -eigenvectors(U) -> axis orientation
U, s, V = np.linalg.svd(D, full_matrices=True) This should be a least-squares fit. Of course, things can get more complicated than this, please see https://www.emis.de/journals/BBMS/Bulletin/sup962/gander.pdf
Using OpenCV fitEllipse:
import numpy as np import cv2 import matplotlib.pyplot as plt def get_xy_Ellipse(x_center, y_center, x_axis, y_axis, angle): # Genrate X and Y points based on the fitted Ellipse ellipse_xy = np.array([x_axis/2*np.cos(t) , y_axis/2*np.sin(t), np.ones_like(t)]) # Rotation and translation rot_trans_matrix = np.array([[np.cos(angle) , -np.sin(angle), x_center], [np.sin(angle) , np.cos(angle), y_center]]) ellipse_xy_rotated = np.matmul(rot_trans_matrix, ellipse_xy) return ellipse_xy_rotated N = 100 t = np.linspace(0, 2*np.pi, N) x_center_gt = 120 y_center_gt = 133 x_axis_gt = 50 y_axis_gt = 40 angle_gt = np.deg2rad(30) ellipse_xy_rotated_gt = get_xy_Ellipse(x_center_gt, y_center_gt, x_axis_gt, y_axis_gt, angle_gt) # Add some noise ellipse_xy_rotated_gt_noisy = ellipse_xy_rotated_gt + np.random.normal(size=ellipse_xy_rotated_gt.shape) ################################################ # Fit an ellipse to the points using OpenCV # The function calculates the ellipse that fits (in a least-squares sense) a set of 2D points best of all. ellipse = cv2.fitEllipse(ellipse_xy_rotated_gt_noisy.T.astype(np.float32)) # Extract ellipse parameters (center, axes, orientation) = ellipse (x_center, y_center) = center (x_axis, y_axis) = axes angle = np.deg2rad(orientation) ################################################ ellipse_xy_rotated = get_xy_Ellipse(x_center, y_center, x_axis, y_axis, angle) plt.figure() plt.scatter(ellipse_xy_rotated_gt[0,:], ellipse_xy_rotated_gt[1,:], label="GT Points") plt.scatter(ellipse_xy_rotated_gt_noisy[0,:], ellipse_xy_rotated_gt_noisy[1,:], label="Noisy Points") plt.scatter( ellipse_xy_rotated[0,:] , ellipse_xy_rotated[1,:], label="Fitted Points") plt.axis('equal') plt.xlabel('X') plt.ylabel('Y') plt.legend() 
See Ellipse wikipedia for equation details.
I found the SVD approaches weren't working for me, maybe because my points are not uniformly spaced?
I settled on a function to construct an ellipse (parametric, given timepoints 0 to 2 pi), an error function (finding the distance between my points and points on the ellipse), and scipy.optimize.minimize.
def ellipse(t,xc,yc,a,b,theta): # works for a list of t points x,y=a*np.cos(t),b*np.sin(t) # start with a scrunched circle c=np.cos(theta) ; s=np.sin(theta) ; R=np.asarray([[c,-s],[s,c]]) x,y=np.matmul(R,[x,y]) # apply rotation matrix return x+xc,y+yc # shift by center position def findEllipse(xs,ys): def dz(args): # error function xc,yc,a,b,theta = args ts = np.linspace(0,2*np.pi,360*3,endpoint=False) x,y=ellipse(ts,xc,yc,a,b,theta) # points for the ellipse for args passed # distance from all given points (xs,ys) to all ellipse points (x,y) distances=np.sqrt( (xs[:,None]-x[None,:])**2+(ys[:,None]-y[None,:])**2 ) # collapse to find each xs,ys points' closest point on ellipse distances = np.amin(distances,axis=1) return np.sqrt(np.sum(distances**2)) # use MSE distance as our error metric # guesses: center in x,y, width and height, zero angle to start x0 = ( np.mean(xs) , np.mean(ys) , np.ptp(xs)/2 , np.ptp(ys)/2 , 0 ) res = minimize(dz,x0) return res.x cx,cy,a,b,theta = findEllipse( xs , ys ) here's an example of my result (I am using the blue pixels to define the "edge" around which i would like to draw the ellipse:
where I found the green pixels with a simple:
mask = np.zeros(data.shape)
mask[data > np.mean(data) + np.std(data)]
and found the blue pixels via:
rolled = mask+np.roll(mask,1,axis=0)+np.roll(mask,-1,axis=0)+\ np.roll(mask,1,axis=1)+np.roll(mask,-1,axis=1) # counts how many neighbors are selected in the mask border = np.where(rolled == 3) # pixels selected by the mask(1) + 2 neighbors 