两个圆盘相交区域的均匀采样

8

如果我们有一个二维均匀变量,我们可以像这里讨论的那样在单位圆中生成均匀分布。

我的问题类似,我希望均匀采样两个相交圆的交集区域,其中一个圆总是单位圆,而另一个圆可以自由移动和调整大小,就像这里一样。

enter image description here

我试图将该区域分成两个区域(如上所示),并根据各自的圆盘对每个区域进行采样。我的方法基于上述均匀圆盘算法。要对中心线右侧的第一个区域进行采样,我会限制θ在两个交点之间。接下来,需要基于该θ投影r,以使点被推入我们的中线和圆盘半径之间的区域。Python示例代码可在此处找到。
u = unifrom2D()
A;B; // Intersection points
for p in allPoints
    theta = u.x * (getTheta(A) - getTheta(B)) + getTheta(B)
    r = sqrt(u.y + (1- u.y)*length2(lineIntersection(theta)))  
    p = (r * cos(theta), r * sin(theta))

然而,这种方法相当昂贵,而且无法保持一致性。只是为了澄清,我不想使用拒绝抽样。

2
理论上,您可以找到相交区域中每条线的长度并进行积分,创建一个“概率分布”,并希望它具有易于计算的反函数。但是,我会使用拒绝抽样,即使我知道您明确拒绝它。 - user354134
3
为什么你不想使用拒绝采样?如果你使用一个合适的边界框来采样,它甚至可以非常高效。 - pjs
1
我想在光线追踪中使用这种方法,根据之前选择的图像样本位置来选择镜头样本点,参见此处。我目前正在blender-cycles-renderer中实现此方法,此处。当前使用的拒绝抽样会导致图像的某些区域接收较少的样本,因此噪声更多。此外,所有样本一次性发送到gpu,为被拒绝的样本请求新样本将非常复杂和昂贵。 - Knork
@barrycarter 我的回答计算了这个“概率分布” - 由于角度和三角函数混合在其中,因此没有解析逆。但是可以通过解析计算导数,因此可以快速进行数值反演。 - coproc
1个回答

4
我不确定这是否比拒绝抽样更好,但是这里有一个解决方案,用于均匀采样圆弧段(中心角度<= pi),涉及计算逆函数的数值。 (两个圆的交点的均匀采样可以由分段、扇形和三角形的采样组成-具体取决于如何将交点分割成更简单的图形。)
首先,我们需要知道如何生成具有给定分布F的随机值Z,即我们想要。
P(Z < x) = F(x)                   <=>  (x = F^-1(y))
P(Z < F^-1(y)) = F(F^-1(y)) = y   <=>  (F is monotonous)
P(F(Z) < y) = y

这意味着:如果 Z 具有所请求的分布 F,那么 F(Z) 是均匀分布的。反过来也是一样的:
Z = F^-1(Y), 

其中Y在区间[0,1]上均匀分布,满足所需的分布。

如果F的形式为

       / 0,                             x < a
F(x) = | (F0(x)-F0(a)) / (F0(b)-F0(a)), a <= x <= b
       \ 1,                             b < x

然后,我们可以在区间[F(a), F(b)]中均匀选择一个Y0,并设置Z = F0^-1(Y0)
我们选择用(theta, r)来参数化该线段,其中中心角度theta是从一个线段侧面测量的。当线段的中心角度为alpha时,与从线段开始的角度为theta的扇形相交的线段面积为(对于单位圆,theta[0, alpha/2]).
F0_theta(theta) = 0.5*(theta - d*(s - d*tan(alpha/2-theta)))

在这里输入图片描述 其中 s = AB/2 = sin(alpha/2) 以及 d = dist(M,AB) = cos(alpha/2) (圆心到线段的距离)。(本文不考虑alpha/2 <= theta <= alpha 的对称情况)。 我们需要一个随机变量theta,满足 P(theta < x) = F_theta(x)。无法通过符号计算得出 F_theta 的逆 - 必须通过优化算法(例如牛顿迭代)确定它。 一旦确定了theta,我们需要在以下范围内随机生成半径r

[r_min, 1], r_min = d/cos(alpha/2-theta).

对于x[0, 1-r_min]范围内,分布必须是

F0_r(x) = (x+r_min)^2 - r_min^2 = x^2 + 2*x*r_min.

这里可以以符号方式计算其反函数:

F0_r^-1(y) = -r_min + sqrt(r_min^2+y)

这里有一个Python的实现,用于概念证明:

from math import sin,cos,tan,sqrt
from scipy.optimize import newton

# area of segment of unit circle
# alpha: center angle of segment (0 <= alpha <= pi)
def segmentArea(alpha):
  return 0.5*(alpha - sin(alpha))

# generate a function that gives the area of a segment of a unit circle
# intersected with a sector of given angle, where the sector starts at one end of the segment. 
# The returned function is valid for [0,alpha/2].
# For theta=alpha/2 the returned function gives half of the segment area.
# alpha: center angle of segment (0 <= alpha <= pi)
def segmentAreaByAngle_gen(alpha):
  alpha_2 = 0.5*alpha
  s,d = sin(alpha_2),cos(alpha_2)
  return lambda theta: 0.5*(theta - d*(s - d*tan(alpha_2-theta)))

# generate derivative function generated by segmentAreaByAngle_gen
def segmentAreaByAngleDeriv_gen(alpha):
  alpha_2 = 0.5*alpha
  d = cos(alpha_2)
  return lambda theta: (lambda dr = d/cos(alpha_2-theta): 0.5*(1 - dr*dr))()

# generate inverse of function generated by segmentAreaByAngle_gen
def segmentAreaByAngleInv_gen(alpha):
  x0 = sqrt(0.5*segmentArea(alpha)) # initial guess by approximating half of segment with right-angled triangle
  return lambda area: newton(lambda theta: segmentAreaByAngle_gen(alpha)(theta) - area, x0, segmentAreaByAngleDeriv_gen(alpha))

# for a segment of the unit circle in canonical position
# (i.e. symmetric to x-axis, on positive side of x-axis)
# generate uniformly distributed random point in upper half
def randomPointInSegmentHalf(alpha):
  FInv = segmentAreaByAngleInv_gen(alpha)
  areaRandom = random.uniform(0,0.5*segmentArea(alpha))
  thetaRandom = FInv(areaRandom)
  alpha_2 = 0.5*alpha
  d = cos(alpha_2)
  rMin = d/cos(alpha_2-thetaRandom)
  secAreaRandom = random.uniform(0, 1-rMin*rMin)
  rRandom = sqrt(rMin*rMin + secAreaRandom)
  return rRandom*cos(alpha_2-thetaRandom), rRandom*sin(alpha_2-thetaRandom)

这个可视化似乎验证了均匀分布(在中心角为pi/2的线段的上半部分):

import matplotlib.pyplot as plot
segmentPoints = [randomPointInSegmentHalf(pi/2) for _ in range(500)]
plot.scatter(*zip(*segmentPoints))
plot.show()

enter image description here


1
@barrycarter 将两个圆的交点分解为两个部分(每个圆一个),如果一个部分包含圆心,则进一步将该部分分解为扇形和三角形是很简单的。(总共将有两个部分或一个圆的部分,另一个圆的扇形和三角形。)相应的面积 A1、A2 等可以很容易地计算出来。要对整个交集进行采样,可以在顶部放置一个开关:Y = uniform(0, A1+A2+...);如果 (Y < A1),则获取区域 1 的采样点,否则如果 (Y-A1 < A2),则获取区域 2 的采样点... - coproc
1
好的,我想我的问题是:为什么不编写解决整个问题的代码,然后显示在两个圆的交集上看起来均匀的分布呢? - user354134
1
@Knork区域(由角度参数化)作为CDF:没问题。当您固定一个角度theta0时,您希望P(点在theta < theta0的区域内)=area(theta0)/area(theta_max)。因此,area(theta0)/area(theta_max)是CDF。 - coproc
1
@Knork 如果您想从“[0,1]”中选择均匀分布的变量,则需要进行归一化。对于编程而言,更容易省略归一化,并从适当的区间选择均匀分布的变量(请参见我在“F0”和“F”之间的区别)。 - coproc
1
@Knork 当交点包含两个圆心时:在你的例子中,你可以将交点分解为红色圆的一段和蓝色圆的一个扇形以及它们之间的三角形 - 对于所有组件,都知道如何均匀采样(请参见我在回答barrycarter的第一个问题时上面的评论)。 - coproc
显示剩余3条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接