Super-ellipse Point Picking

Question

https://en.wikipedia.org/wiki/Superellipse

I have read the SO questions on how to point-pick from a circle and an ellipse.

How would one uniformly select random points from the interior of a super-ellipse?

More generally, how would one uniformly select random points from the interior of the curve described by an arbitrary super-formula?

https://en.wikipedia.org/wiki/Superformula

The discarding method is not considered a solution, as it is mathematically unenlightening.

Answer 1

In order to sample the superellipse, let's assume without loss of generality that a = b = 1. The general case can be then obtained by rescaling the corresponding axis.

The points in the first quadrant (positive x-coordinate and positive y-coordinate) can be then parametrized as:

x = r * ( cos(t) )^(2/n)
y = r * ( sin(t) )^(2/n)

with 0 <= r <= 1 and 0 <= t <= pi/2:

Now, we need to sample in r, t so that the sampling transformed into x, y is uniform. To this end, let's calculate the Jacobian of this transform:

dx*dy = (2/n) * r * (sin(2*t)/2)^(2/n - 1) dr*dt
      = (1/n) * d(r^2) * d(f(t))

Here, we see that as for the variable r, it is sufficient to sample uniformly the value of r^2 and then transform back with a square root. The dependency on t is a bit more complicated. However, with some effort, one gets

f(t) = -(n/2) * 2F1(1/n, (n-1)/n, 1 + 1/n, cos(t)^2) * cos(t)^(2/n)

where 2F1 is the hypergeometric function.

In order to obtain uniform sampling in x,y, we need now to sample uniformly the range of f(t) for t in [0, pi/2] and then find the t which corresponds to this sampled value, i.e., to solve for t the equation u = f(t) where u is a uniform random variable sampled from [f(0), f(pi/2)]. This is essentially the same method as for r, nevertheless in that case one can calculate the inverse directly.

One small issue with this approach is that the function f is not that well-behaved near zero - the infinite slope makes it quite challenging to find a root of u = f(t). To circumvent this, we can sample only the "upper part" of the first quadrant (i.e., area between lines x=y and x=0) and then obtain all the other points by symmetry (not only in the first quadrant but also for all the other ones).

An implementation of this method in Python could look like:

import numpy as np
from numpy.random import uniform, randint, seed
from scipy.optimize import brenth, ridder, bisect, newton
from scipy.special import gamma, hyp2f1

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

seed(100)

def superellipse_area(n):
    #https://en.wikipedia.org/wiki/Superellipse#Mathematical_properties
    inv_n = 1. / n
    return 4 * ( gamma(1 + inv_n)**2 ) / gamma(1 + 2*inv_n)

def sample_superellipse(n, num_of_points = 2000):

    def f(n, x):
        inv_n = 1. / n
        return -(n/2)*hyp2f1(inv_n, 1 - inv_n, 1 + inv_n, x)*(x**inv_n)

    lb = f(n, 0.5)
    ub = f(n, 0.0)

    points = [None for idx in range(num_of_points)]
    for idx in range(num_of_points):
        r = np.sqrt(uniform())

        v = uniform(lb, ub)

        w = bisect(lambda w: f(n, w**n) - v, 0.0, 0.5**(1/n))
        z = w**n

        x = r * z**(1/n)
        y = r * (1 - z)**(1/n)

        if uniform(-1, 1) < 0:
            y, x = x, y
        x = (2*randint(0, 2) - 1)*x
        y = (2*randint(0, 2) - 1)*y

        points[idx] = [x, y]
    return points

def plot_superellipse(ax, n, points):
    coords_x = [p[0] for p in points]
    coords_y = [p[1] for p in points]

    ax.set_xlim(-1.25, 1.25)
    ax.set_ylim(-1.25, 1.25)
    ax.text(-1.1, 1, '{n:.1f}'.format(n = n), fontsize = 12)
    ax.scatter(coords_x, coords_y, s = 0.6)

params = np.array([[0.5, 1], [2, 4]])

fig = plt.figure(figsize = (6, 6))
gs = gridspec.GridSpec(*params.shape, wspace = 1/32., hspace = 1/32.)

n_rows, n_cols = params.shape
for i in range(n_rows):
    for j in range(n_cols):
        n = params[i, j]
        ax = plt.subplot(gs[i, j])

        if i == n_rows-1:
            ax.set_xticks([-1, 0, 1])
        else:
            ax.set_xticks([])

        if j == 0:
            ax.set_yticks([-1, 0, 1])
        else:
            ax.set_yticks([])

        #ensure that the ellipses have similar point density
        num_of_points = int(superellipse_area(n) / superellipse_area(2) * 4000)

        points = sample_superellipse(n, num_of_points)
        plot_superellipse(ax, n, points)

fig.savefig('fig.png')

This produces: