Abhilash Reddy M, Originally written sometime in 2016. updated on June 03 2022
This notebook describes methods for uniformly random sampling over some common shapes. This is a python notebook that can be found here. There is a follow-up piece here. Use the table of contents below to jump to a particular section.
We start with a small amount of boiler plate code to import NumPy and Matplotlib modules. We are using two distributions functions from the library here: uniform() and standard_normal(). If your programming language does not come with a standard normal distribution function, you can make use of the Box-Muller transform to generate standard normal distributed samples from a uniform distribution. So, just a uniform random generator is sufficient to use any of the methods described. In C codes, one finds u=((double)rand())/RAND_MAX; as a means of generating uniform random samples. This is not the best uniform distribution one can get with 8-byte reals, but it should be more than enough for most cases. C++ 2011 introduced the random header that has more random number generators.
#set up the notebook and the environment
%matplotlib inline
from jupyterthemes import jtplot
jtplot.style() # set the style of plots to match the notebook theme
import numpy as np
from numpy.random import uniform,standard_normal
import matplotlib.pyplot as plt
I will assume that we have access to a 1D uniform distribution sampling function and a 1D normal distribution sampling function. So, uniform sampling on a line becomes trivial by mapping the line to $[0,1)$. Then uniformly random draws from $[0,1)$ can be mapped back to the line.
N=1000
verts=np.asarray([[0,0],[0.7,1],[1,0.5]])
fig,ax=plt.subplots(figsize=[6,6])
_=ax.add_artist(plt.Polygon(verts,alpha=0.3,lw=2,ec='k'))
p1 = np.sqrt(uniform(0,1,N)) # taking sqrt here itself
p2 = uniform(0,1,N)
ts=np.empty([N,2])
ts[:,0] = (1-p1)*verts[0,0] + p1*(1-p2)*verts[1,0] + p2*p1*verts[2,0]
ts[:,1] = (1-p1)*verts[0,1] + p1*(1-p2)*verts[1,1] + p2*p1*verts[2,1]
_=ax.plot(ts[:,0],ts[:,1],'ro')
This works for any triangle.
An alternative solution is to form a parallellogram using any two sides of the triangle. And generate random points in the parallellogram (shown below). Some points will fall outside the main triangle and in the second triangle. These points can be "folded" back, across the edge, on to the main triangle and used as a uniform sample.
This is pretty straightforward. Each point is made of two independent random draws. This works for rectangles and cubes and cuboids just as well.
N=1000
verts=np.asarray([[0,0],[0,1],[1,1],[1,0]])
fig,ax=plt.subplots(figsize=[6,6])
_=ax.add_artist(plt.Polygon(verts,alpha=0.3,lw=2,ec='k'))
u1 = uniform(0,1,N)
u2 = uniform(0,1,N)
_=ax.plot(u1,u2,'ro')
For a general square that is not aligned with the coordinate axis, we could calculate a transformation that maps it to a reference unit square. We can then sample the points on the reference geometry and use the inverse transformation to get the distribution on the actual geometry. However, there is an easier way shown below that works for any arbitrary parallellogram. Any point inside a parallellogram can be written in terms of the two edge vectors that form the parallellogram. This will not work for a general quadrilateral, which will be dealt with further down.
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111, aspect='equal', autoscale_on=True)
verts=2.5+np.asarray([[0,0],[0.2,0.8],[1,1],[0.8,0.2]])
_=ax.add_artist(plt.Polygon(verts,alpha=0.3,lw=2,ec='k'))
N=4000
u1 = uniform(0,1,N)
u2 = uniform(0,1,N)
e1= verts[1]-verts[0]
e2= verts[3]-verts[0]
p=verts[0] + u1[:,np.newaxis]*e1+ u2[:,np.newaxis]*e2
_=ax.plot(p[:,0],p[:,1],'ro', alpha=0.5 , mfc='r', mec='none')
For a general quadrilateral, there is not a method analogous to what we did for a triangle. The mapping from any triangle to any other triangle is an affine transformation. The same mapping for quadrilaterals is non-linear and not affine in general. Parallellograms are an affine transformation of squares so we could deal with parallellograms easily.
We can use the triangle sampler from above to generate points uniformly inside a polygon in the following way:
We can make use of the general rejection sampling method.
This is a very general method. The following demonstrates this with a general polygon. There is one additional task here: testing if a point is inside a polygon. There are multiple ways to check this. If the polygon is convex, there is a very fast test that can be used:
def pinp_c(p,q):
"""
Determine if point p is inside convex polygon q
Returns True if p is inside q, False otherwise
we loop over all edges of q and check if p is on the same side of all the edges
If yes, then p is inside q
Input:
p: point to be tested, array of size 2
q: vertices of convex polygon, Nx2 array, have to be ordered CW or CCW
"""
d=np.cross(np.roll(q,-1,axis=0)-q,p-q)
return np.all(d>0) or np.all(d<0)
The following function uses a ray casting method and works for compex polygons including non-convex polygons, polygons with holes, and inverted polygons (based on an algorithm in Paul Bourke's site attributed to Dr. Randolph Franklin)
def pinp_nc(p,q):
"""
Determine if point p is inside general non-convex polygon q
Returns True if p is inside q, False otherwise
This uses the ray casting method. based on a Paul Bourke article
"""
c=False
j=len(q)-1
for i in range(len(q)):
t1=(q[i,1]<=p[1] and p[1]<q[j,1]) or (q[j,1]<=p[1] and p[1]<q[i,1])
t2=p[0]<(q[j,0]-q[i,0])*(p[1]-q[i,1])/(q[j,1]-q[i,1])+q[i,0]
if t1 and t2:
c=not c
j=i
return c
Following code shows generation of points in arbitrary polygons using rejection based on the function pinp_nc. Try changing the polygon defnition and the test function to pinp_nc and see what happens. There are other methods like the winding number method or the angle subtended method also which might be used.
If you can guarantee that the polygon will be convex then the convex-only test should be used because it is very fast.
fig,ax=plt.subplots(figsize=[6,6])
ax.set_aspect('equal')
poly=2.5 + np.asarray([[-0.1,0],[1.4,0.1],[1.6,0.7], [1.0,0.6] ,[0.3,0.8], [0.1,1.2]])
_=ax.add_artist(plt.Polygon(poly,alpha=0.3,lw=2,ec='k'))
bbox=np.array([poly.min(axis=0),poly.max(axis=0)])
N=1000
ns=0
while ns<N:
u1=uniform()
u2=uniform()
p=bbox[0]+np.array([u1*(bbox[1,0]-bbox[0,0]),u2*(bbox[1,1]-bbox[0,1])])
if pinp_nc(p,poly):
ns+=1
_=ax.plot(p[0],p[1],'ro', alpha=0.8, mec='none')
else:
_=ax.plot(p[0],p[1],'go', alpha=0.2, mec='none')
Finally, if a lot of points are needed inside a large complicated polygon (imagine the boundary of a meshed region), it will be faster to use an auxiliary background grid and use a distance function based method to test the points.
Now, we move on to sampling inside circles. in a disk. in the 2-d ball. First a naive implementation based on the polar parametrization of the circle is demonstrated to produce wrong distribution. Then, several procedures that generate proper samples are described.
We want to generate $N$ random points inside the circle, such that they are distributed with an equal density. Using polar coordinates, $r$ and $\theta$, any point inside a circle or radius $r_0$ can be represented as
$$(x,y)=(r~cos(\theta),r~sin(\theta))$$for $$r\in[0,r_0),~\theta\in[0,2\pi) $$
The naive approach is to randomly pick $N$ values each of $r\in[0,r_0]$ and a $\theta\in[0,2\pi]$ and calculate $x$ and $y$. As we shall see this does not give us a uniform distribution inside the circle. Lets say $N=1000$ and $r_0=2$.
r0 = 2.0
N = 1000
r = r0*uniform(0,1,N)
t = uniform(0,2*np.pi,N)
x,y= r*np.cos(t),r*np.sin(t)
fig,ax=plt.subplots(figsize=(6,6), subplot_kw=dict(aspect='equal'))
ax.add_artist(plt.Circle((0, 0), r0,alpha=0.25))
ax.plot(x,y,'ro', alpha=0.5, mec='none')
[<matplotlib.lines.Line2D at 0x2934edb7100>]