I developed a simple function for estimating the area of a possibly cropped circle inside a frame
It is simple, and that has value. It also has (severe) cost. It's very slow and very inaccurate. The worst case is when the intersection area is small - the relative error skyrockets.
Unfortunately, to make it fast and accurate takes a lot of work. The math isn't particularly difficult, there's just a lot of it.
This:
shows that there are many, many edge cases that need to be accommodated for; here d intersects eg instead of ea due to a narrow rectangle.
Specific to your original implementation,
- Pass in an optional random generator for test reproducibility.
- Pre-compute
xr[msk_in_circ] and yr[msk_in_circ] before using them multiple times in an expression. - Prefer count_nonzero over a simple sum.
import time from collections import defaultdict from pprint import pprint import numpy as np from matplotlib import pyplot as plt from scipy.spatial import distance def circ_frac_montecarlo( cx: float, cy: float, rad: float, x0: float, x1: float, y0: float, y1: float, n: int = 100_000, rand: np.random.Generator | None = None, ): """ Use Monte Carlo to estimate the fraction of the area of a circle centered in (cx, cy) with a radius of 'rad', that is located within the frame given by the limits 'x0, x1, y0, y1'. """ if rand is None: rand = np.random.default_rng() # Generate a square containing the circle. xmin, xmax = cx - rad, cx + rad ymin, ymax = cy - rad, cy + rad # Generate 'N_tot' uniform random points inside that square. xr = rand.uniform(xmin, xmax, n) yr = rand.uniform(ymin, ymax, n) # Obtain the distances of those points to the center of the circle. # Find the points within the circle. dist = distance.cdist([(cx, cy)], np.array([xr, yr]).T)[0] msk_in_circ = dist < rad # Find the points that are within the frame, from the points that are # within the circle. xmask = xr[msk_in_circ] ymask = yr[msk_in_circ] msk_xy = ( (xmask > x0) & (xmask < x1) & (ymask > y0) & (ymask < y1) ) # The area is the points within circle and frame over the points within circle. return np.count_nonzero(msk_xy) / np.count_nonzero(msk_in_circ) def intersect( ex: float, ey: float, # line origin cx: float, cy: float, # circular intersection ) -> float: """ Single-axis intersection of a ray (0c) and either of x=ex or y=ey Half of the callers transpose x and y. """ num = ey*cx if num == 0: return ex ix = max(ex, num/cy) return ix def quadrant_area(x0: float, x1: float, y0: float, y1: float) -> float: """ Per quadrant, calculate the intersection of a unit circle and a rectangle. The only constraints on the rectangle are x0 <= x1, y0 <= y1. In the most complex case, the intersection looks like | 1 + | theta1 y1 + a+-----)b | | / ) theta0 | | / )c | d|/ /| | /| f/ | y0 + /e+------+--+g | / / |/ / |/ / |/ / 0 +---+---------+--+-- x0 x1 1 abd: right triangle, include in area cfg: right triangle, include in area bcfed: truncated sector, include in area 0de: scalene, exclude from sector area; or if d is below e: aed*2: rectangle, include in area 0ef: scalene, exclude from sector area; or if f is above e: efg*2: rectangle, include in area """ # clip to within [0,1], [0,1] x0, x1, y0, y1 = np.clip((x0, x1, y0, y1), a_min=0, a_max=1) ex, ey = x0, y0 # intersection of x^2+y^2=1 and y=y1 bx, by = np.sqrt(1 - y1*y1), y1 if bx < x0: # intersection of x^2+y^2=1 and x=x0 bx, by = x0, np.sqrt(1 - x0*x0) if by < y0: # rectangle's upper right corner before the circle; simple rectangular area return (x1 - x0)*(y1 - y0) elif bx > x1: return (x1 - x0)*(y1 - y0) # intersection of x^2+y^2=1 and x=x1 cx, cy = x1, np.sqrt(1 - x1*x1) if cy < y0: # intersection of x^2+y^2=1 and y=y0 cx, cy = np.sqrt(1 - y0*y0), y0 if cx < x0: return (x1 - x0)*(y1 - y0) elif cy > y1: return (x1 - x0)*(y1 - y0) ax, ay = ex, by gx, gy = cx, ey dx = intersect(ex, ey, bx, by) # intersection of y=ey and y=byx/bx dy = intersect(ey, ex, by, bx) # intersection of x=ex and y=byx/bx fx = intersect(ex, ey, cx, cy) # intersection of y=ey and y=cyx/cx fy = intersect(ey, ex, cy, cx) # intersection of x=ex and y=cyx/cx aed2 = (dx - ex)*(ay - ey) efg2 = (gx - ex)*(fy - ey) abd = 0.5*(bx - dx)*(by - dy) cfg = 0.5*(cx - fx)*(cy - fy) # Signed-heron area de = dy - ey + ex - dx abs_de = np.abs(de) d0 = np.hypot(dx, dy) e0 = np.hypot(ex, ey) sp = 0.5*(d0 + e0 + abs_de) # semiperimeter de0 = np.copysign(np.sqrt(sp*(sp - d0)*(sp - e0)*(sp - abs_de)), de) ef = fx - ex + ey - fy abs_ef = np.abs(ef) f0 = np.hypot(fx, fy) sp = 0.5*(e0 + f0 + abs_ef) ef0 = np.copysign(np.sqrt(sp*(sp - e0)*(sp - f0)*(sp - abs_ef)), ef) t0 = np.atan2(cy, cx) t1 = np.atan2(by, bx) sector = 0.5*(t1 - t0) # * 1**2 area = max(0, aed2 + efg2 + abd + cfg + sector - de0 - ef0) # dirty; do not use this in production # _debug_quadrant(**locals()) return area def _debug_quadrant( x0,y0, x1,y1, ax,ay, bx,by, cx,cy, dx,dy, ex,ey, fx,fy, gx,gy, area, **kwargs, ) -> None: plt.figure() plt.plot( [x0, x0, x1, x1, x0], [y0, y1, y1, y0, y0], marker='o', ) plt.plot( [0, dx, bx, cx, fx, 0], [0, dy, by, cy, fy, 0], marker='o', ) coords = { 'a': (ax, ay), 'b': (bx, by), 'c': (cx, cy), 'd': (dx, dy), 'e': (ex, ey), 'f': (fx, fy), 'g': (gx, gy), } pprint({ **coords, 'area': area, 'x0': x0, 'x1': x1, 'y0': y0, 'y1': y1, }) by_value = defaultdict(list) for name, (x, y) in coords.items(): by_value[x, y].append(name) for loc_coords, names in by_value.items(): plt.annotate( ','.join(names), loc_coords, ) def circ_frac_analytic( cx: float, cy: float, rad: float, x0: float, x1: float, y0: float, y1: float, ) -> float: # The fractional area is invariant under affine transformation. # Eliminate cx,cy,rad by transforming to a unit circle. x0 = (x0 - cx)/rad x1 = (x1 - cx)/rad y0 = (y0 - cy)/rad y1 = (y1 - cy)/rad ''' q1 | q0 ---+--- q2 | q3 ''' q0 = quadrant_area(x0=x0, x1=x1, y0=y0, y1=y1) q1 = quadrant_area(x0=-x1, x1=-x0, y0=y0, y1=y1) q2 = quadrant_area(x0=-x1, x1=-x0, y0=-y1, y1=-y0) q3 = quadrant_area(x0=x0, x1=x1, y0=-y1, y1=-y0) intersection_area = q0 + q1 + q2 + q3 circle_area = np.pi return intersection_area/circle_area def test_100k() -> None: """OP's original test scheme: big square, small circle""" rand = np.random.default_rng(seed=0) # Define the (x, y) limits of the frame x0, x1 = 0., 1. y0, y1 = 0., 1. xyr = rand.uniform(low=(x0, y0, 0.05), high=(x1, y1, 0.5), size=(40, 3)) fracs = [ circ_frac_montecarlo(cx, cy, rad, x0, x1, y0, y1, rand=rand) for cx, cy, rad in xyr ] afracs = [ circ_frac_analytic(cx, cy, rad, x0, x1, y0, y1) for cx, cy, rad in xyr ] # allow for Monte Carlo error assert np.allclose(fracs, afracs, atol=0, rtol=6e-3) def test_1m() -> None: rand = np.random.default_rng(seed=0) axis = np.linspace(start=-2, stop=1, num=11) # Due to a needle-like sector, the worst Monte Carlo error is seen for e.g. # x0=-1.7 x1=-0.7 y0=-0.8 y1=-0.7 for x0 in axis: x1 = x0 + 1 for y0 in axis: y1 = y0 + 0.1 ref = circ_frac_montecarlo(cx=0, cy=0, rad=1, x0=x0, y0=y0, x1=x1, y1=y1, n=1_000_000, rand=rand) analytic = circ_frac_analytic(cx=0, cy=0, rad=1, x0=x0, y0=y0, x1=x1, y1=y1) assert np.allclose(ref, analytic, rtol=0, atol=4e-4) # plt.show() # if debug() is called def benchmark() -> None: rand = np.random.default_rng(seed=0) n = 1_000_000 kwargs = { 'cx': 0.1, 'cy': -0.1, 'rad': 1.1, 'x0': -0.2, 'y0': -0.3, 'x1': 0.8, 'y1': 0.7, } t0 = time.perf_counter() circ_frac_montecarlo(**kwargs, n=n, rand=rand) t1 = time.perf_counter() circ_frac_analytic(**kwargs) t2 = time.perf_counter() print(f'Monte Carlo: {1e3*(t1-t0):.1f} ms') print(f'Analytic: {1e6*(t2-t1):.1f} us') if __name__ == '__main__': test_100k() test_1m() benchmark()
Monte Carlo: 29.1 ms Analytic: 83.9 us
