## In which Ellipses are Dealt with Handily

I’m making a user interface where people can draw a lasso select on an image, and I want to see if certain points in the image are within that lasso select. Now that I think about it, I could probably have done some form of polygon membership test (e.g. raycasting) for each of the points that I wanted to check, with the points of the lasso forming the polygon, but that has a few problems. One of them is that the ends of the lasso can overlap, but that can be handled by e.g. a convex hull. In fact, I may end up doing that, because the current approach has the problem that it can’t handle an arbitrarily-shaped lasso selection.

What I ended up doing as a first cut, however, has some other useful properties. I had the software find a minimum enclosing ellipse of the points for the lasso selection, and then tested the points I wanted checked to see if they were inside the ellipse.

```def min_bounding_oval(points, tolerance=0.01):
#Create a numpy 2xN matrix from a list of 2D points
P = np.matrix(points).T

# Pseudocode from https:#stackoverflow.com/questions/1768197/bounding-ellipse
# Pythonificiation by me
# Input: A 2xN matrix P storing N 2D points
#        and tolerance = tolerance for error.
# Output: The equation of the ellipse in the matrix form,
#         i.e. a 2x2 matrix A and a 2x1 vector C representing
#         the center of the ellipse.

# Dimension of the points
d = 2;
# Number of points
N = len(points);

# Add a row of 1s to the 2xN matrix P - so Q is 3xN now.
Q = np.vstack([P,np.ones((1,N))])

# Initialize
count = 1;
err = 1;
#u is an Nx1 vector where each element is 1/N
u = (1.0/N) * np.ones((N,1))

# Khachiyan Algorithm
while err > tolerance:
# Matrix multiplication:
X = Q * np.diagflat(u) * Q.T

M = np.diagonal(Q.T * X.I * Q)

# Find the value and location of the maximum element in the vector M
maximum = M.max()
j = np.argmax(M);

# Calculate the step size for the ascent
step_size = (maximum - d -1)/((d+1)*(maximum-1));

# Calculate the new_u:
# Take the vector u, and multiply all the elements in it by (1-step_size)
new_u = (1 - step_size) * u ;

# Increment the jth element of new_u by step_size
new_u[j] = new_u[j] + step_size;

# Store the error by taking finding the square root of the SSD
# between new_u and u
err = math.sqrt(((new_u - u)**2).sum());

# Increment count and replace u
count = count + 1;
u = new_u;

# Put the elements of the vector u into the diagonal of a matrix
# U with the rest of the elements as 0
U = np.diagflat(u);

# Compute the A-matrix
A = (1.0/d) * (P * U * P.T - (P * u)*(P*u).T ).I

# And the center,
c = P * u

return [A, c]```

That gets me A, a 2×2 matrix representing the ellipse and c, a 2×1 matrix representing the center of the ellipse.

To check if a point is in the ellipse, some fiddling has to be done to A.

```
#Convert A to a whitening matrix (W = A**-1/2)
w, v = np.linalg.eig(A)
D = np.diagflat(w)
W = v * np.sqrt(D) * v.I

#Subtract the center of the ellipse and use the whitening matrix
#tag_x and tag_y are the x and y positions of the point to check
p = np.matrix([[tag_x],[tag_y]])
p_center = p - c
p_white = W * p_center

#Check if the whitened point is in the ellipse
if np.linalg.norm(p_white) <= 1:
return True #or whatever you want to do with it```

Note that I haven’t proven this code correct or anything, just run it and had it work for me.