I have a rectangular frame and a circle with center and radius randomly generated. The center is always located within the limits of the frame, as shown:
I need to estimate the area of the fraction of the circle that is located within the frame. Currently I employ a simple Monte Carlo estimate that works ok, but I'd like to compare this with an exact geometric estimation of this area.
Is there a library and/or method to do this? I'm open to pretty much anything that can be installed with conda
or pip
.
import numpy as np
import matplotlib.pyplot as plt
def circFrac(cx, cy, rad, x0, x1, y0, y1, N_tot=100000):
"""
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'.
"""
# Source: https://stackoverflow.com/a/50746409/1391441
r = rad * np.sqrt(np.random.uniform(0., 1., N_tot))
theta = np.random.uniform(0., 1., N_tot) * 2 * np.pi
xr = cx + r * np.cos(theta)
yr = cy + r * np.sin(theta)
# Points within the circle that are within the frame.
msk_xy = (xr > x0) & (xr < x1) & (yr > y0) & (yr < y1)
# The area is the points within circle and frame over the points within
# circle.
return msk_xy.sum() / N_tot
for _ in range(10):
# Random (x, y) limits of the frame
x0, y0 = np.random.uniform(0., 500., 2)
x1, y1 = np.random.uniform(500., 1000., 2)
# Random center coordinates *always* within the frame
cx = np.random.uniform(x0, x1)
cy = np.random.uniform(y0, y1)
# Random radius
rad = np.random.uniform(10., 500)
frac = circFrac(cx, cy, rad, x0, x1, y0, y1)
plt.xlim(x0, x1)
plt.ylim(y0, y1)
circle = plt.Circle((cx, cy), rad, fill=False)
plt.gca().add_artist(circle)
plt.scatter(
cx, cy, marker='x', c='r', label="({:.0f}, {:.0f}), r={:.0f}".format(
cx, cy, rad))
plt.legend()
plt.title("Fraction of circle inside frame: {:.2f}".format(frac))
plt.axes().set_aspect('equal')
plt.show()
You can use shapely for that:
import shapely.geometry as g
xr,yr,r = 0,0,5
circle = g.Point(xr,yr).buffer(r)
x1,y1,x2,y2 = -1,-2,3,5
rectangle = g.Polygon([(x1,y1),(x1,y2),(x2,y2),(x2,y1)])
intersection = rectangle.intersection(circle)
intersection.area