Following the method suggested here, I've written a Python script that generates random points along the perimeter of an ellipse. I wanted to compare a "naive" method (which fails to produce an even distribution) with the correct method, which works like this:
- Numerically compute P, the perimeter of the ellipse.
- Pick a random number s on [0,P].
- Determine the angle t associated with that arc length s.
- Compute the x- and y-coordinates associated with t from the parametric form of the ellipse.
Here is the code I am hoping to receive feedback on. I ran it through
import numpy as np import matplotlib.pyplot as plt def ellipse_arc(a, b, theta, n): """Cumulative arc length of ellipse with given dimensions""" # Divide the interval [0 , theta] into n steps at regular angles t = np.linspace(0, theta, n) # Using parametric form of ellipse, compute ellipse coord for each t x, y = np.array([a * np.cos(t), b * np.sin(t)]) # Compute vector distance between each successive point x_diffs, y_diffs = x[1:] - x[:-1], y[1:] - y[:-1] cumulative_distance =  c = 0 # Iterate over the vector distances, cumulating the full arc for xd, yd in zip(x_diffs, y_diffs): c += np.sqrt(xd**2 + yd**2) cumulative_distance.append(c) cumulative_distance = np.array(cumulative_distance) # Return theta-values, distance cumulated at each theta, # and total arc length for convenience return t, cumulative_distance, c def theta_from_arc_length_constructor(a, b, theta=2*np.pi, n=100): """ Inverse arc length function: constructs a function that returns the angle associated with a given cumulative arc length for given ellipse.""" # Get arc length data for this ellipse t, cumulative_distance, total_distance = ellipse_arc(a, b, theta, n) # Construct the function def f(s): assert np.all(s <= total_distance), "s out of range" # Can invert through interpolation since monotonic increasing return np.interp(s, cumulative_distance, t) # return f and its domain return f, total_distance def rand_ellipse(a=2, b=0.5, size=50, precision=100): """ Returns uniformly distributed random points from perimeter of ellipse. """ theta_from_arc_length, domain = theta_from_arc_length_constructor(a, b, theta=2*np.pi, n=precision) s = np.random.rand(size) * domain t = theta_from_arc_length(s) x, y = np.array([a * np.cos(t), b * np.sin(t)]) return x, y def rand_ellipse_bad(a, b, n): """ Incorrect method of generating points evenly spaced along ellipse perimeter. Points cluster around major axis. """ t = np.random.rand(n) * 2 * np.pi return np.array([a * np.cos(t), b * np.sin(t)])
And some test visualizations:
np.random.seed(4987) x1, y1 = rand_ellipse_bad(2, .5, 1000) x2, y2 = rand_ellipse(2, .5, 1000, 1000) fig, ax = plt.subplots(2, 1, figsize=(13, 7), sharex=True, sharey=True) fig.suptitle('Generating random points on perimeter of ellipse', size=18) ax.set_aspect('equal') ax.set_aspect('equal') ax.scatter(x1, y1, marker="+", alpha=0.5, color="crimson") ax.scatter(x2, y2, marker="+", alpha=0.5, color="forestgreen") ax.set_title("Bad method: Points clustered along major axis") ax.set_title("Correct method: Evenly distributed points")
# Plot arc length as function of theta theta_from_arc_length, domain = theta_from_arc_length_constructor(2, .5, theta=2*np.pi, n=100) s_plot = np.linspace(0, domain, 100) t_plot = theta_from_arc_length(s_plot) fig, ax = plt.subplots(figsize=(7,7), sharex=True, sharey=True) ax.plot(t_plot, s_plot) ax.set_xlabel(r'$\theta$') ax.set_ylabel(r'cumulative arc length')
I would appreciate a general review, but I also have a few specific questions:
- How are my comments? In my code, I know many of my comments are explaining "how" rather than "why," and my inclination is simply to remove them, since I think I have made my variable and function names clear. But I would appreciate an example of how a comment like
# Compute vector distance between each successive pointcould be rewritten in "why" terms.
- Have I computed the arc length in the most efficient way possible? I started by generating the points in two lists of x- and y-coordinates, then iterated over these lists to get the distance made at each step, cumulating those distances as I went around.
- Is it there anything unconventional or inefficient about using
theta_from_arc_length_constructorto create the inverse arc-length function? My thinking was that I need to evaluate the inverse of the arc-length function for every s, so I should go ahead and "prepare" this function rather than recalculate the total arc length each time. But doesn't
return np.interp(s, cumulative_distance, t)mean that numpy has to do the interpolation every time I call f ? Or does it get "automatically" vectorized because I feed it s as an array later on when creating the graphs?