Uniform sampling of the area of ​​intersection of two disks

Question

For a two-dimensional homogeneous variable, we can generate a uniform distribution in a unit disk, as discussed here .

My problem is similar in that I want to evenly try the intersection area of ​​two intersecting disks, where one disk is always a single disk, and the other can be freely moved and changed, for example here

I tried to divide an area into two areas (as shown above) and try each individual region based on a reputable disk. My approach is based on the single disc algorithm given above. To try the first area to the right of the center line, I would limit theta to two intersection points. Next you should project r based on the theta being so that the points are pressed in the area between our middle line and the radius of the disk. The python sample code can be found here .

u = unifrom2D() A;B; // Intersection points for p in allPoints theta = ux * (getTheta(A) - getTheta(B)) + getTheta(B) r = sqrt(uy + (1- uy)*length2(lineIntersection(theta))) p = (r * cos(theta), r * sin(theta)) 

However, this approach is quite expensive and does not allow to maintain uniformity. To clarify, I do not want to use sampling.

Answer 1

I'm not sure if this is better than reject sampling, but here is a solution for uniformly sampling a circle segment (with a central angle <= pi) using a numerical calculation of the inverse function. (A single sample of the intersection of two circles can then be composed of a sample of segments, sectors, and triangles, depending on how the intersection can be divided into simpler shapes.)

First, we need to know how to create a random value Z with a given distribution F , that is, we want

 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 

This means: if Z has the requested distribution F , then F(Z) is evenly distributed. Another way:

 Z = F^-1(Y), 

where Y evenly distributed over [0,1] , has the requested distribution.

If F has the form

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

then we can choose a Y0 uniformly in [F(a),F(b)] and set Z = F0^-1(Y0) .

We choose the parameterization of the segment at (theta,r) , where the central angle of theta is measured on one side of the segment. When the center angle of the segment is alpha , the area of ​​the segment intersects with the sector of the angle theta , starting at the beginning of the segment (for the unit circle theta in [0,alpha/2] )

 F0_theta(theta) = 0.5*(theta - d*(s - d*tan(alpha/2-theta))) 

where s = AB/2 = sin(alpha/2) and d = dist(M,AB) = cos(alpha/2) (the distance from the center of the circle to the segment). (The case alpha/2 <= theta <= alpha symmetric and is not considered here.) We need a random theta with P(theta < x) = F_theta(x) . The inverse to F_theta cannot be calculated symbolically - it must be determined by some optimization algorithm (for example, Newton-Raphson).

Once theta fixed, we need a random radius r in the range

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

For x in [0, 1-r_min] distribution should be

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

Here the inverse can be calculated symbolically:

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

Here is a Python implementation to prove the concept:

 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 # (ie 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) 

The visualization, apparently, checks the uniform distribution (of the upper half of the segment with the central angle pi/2 ):

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