超椭圆点拾取

Question

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

我已阅读有关如何从圆和椭圆中进行点选的 SO 问题。

如何从超椭圆内部均匀地选择随机点？

更一般地说，如何从任意超公式描述的曲线内部均匀地选择随机点？

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

丢弃方法不被视为解决方案，因为它在数学上没有启发性。

Answer 1

为了对超椭圆进行采样，我们不失一般性地假设a = b = 1。然后通过重新缩放相应的轴可以得到一般情况。

第一象限中的点（正x坐标和正y坐标）则可以parametrized为：

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

与0 <= r <= 1 和0 <= t <= pi/2：

现在，我们需要在r, t 中采样，以便转换为x, y 的采样是统一的。为此，我们来计算一下这个变换的雅可比：

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

在这里，我们看到对于变量r，对r^2的值进行均匀采样，然后用平方根变换回来就足够了。对t 的依赖有点复杂。然而，通过一些努力，一个人得到了

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

其中2F1 是hypergeometric 函数。

为了得到x,y中的均匀采样，我们现在需要在[0, pi/2]中对t的f(t)的范围进行均匀采样，然后找到这个采样值对应的t，即求解t 方程u = f(t)，其中u 是从[f(0), f(pi/2)] 采样的均匀随机变量。这与r 的方法基本相同，但在这种情况下可以直接计算逆。

这种方法的一个小问题是函数f 在接近零时表现不佳 - 无限斜率使得找到u = f(t) 的根非常具有挑战性。为了避免这种情况，我们可以只采样第一象限的“上部”（即线x=y和x=0之间的区域），然后通过对称获得所有其他点（不仅在第一象限，而且对于所有其他的）。

此方法在 Python 中的实现可能如下所示：

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')

这会产生：