pythonareamontecarlogeometry-surface

Geometric estimation of area of cropped circle


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:

enter image description here

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

Solution

  • 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
    

    enter image description here