Month: August 2018
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.