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>]
Clearly, we do not get a uniform distribution of points in the cirle this way. The points are clustered near the center. Lets take a step back and see why this went wrong. Lets say we already have a uniform sampling inside the circle, then on integrating the point density over $\theta$ we would find the the density along $r$ is not constant (i.e not uniform). In fact, the density at any $r$ would be proportional to the value of $r$ at that point($=2r$ specifically). This means drawing $r$ from a uniform distribution was not right.
We need to draw $r$ from a distribution that starts at 0 and increases linearly till $r=r_0$. We can use the Inverse transform sampling procedure to draw from the correct distribution. For this to work, we need to integrate the $PDF$ to get the $CDF$ and then invert the $CDF$. Applying the inverse of the $CDF$ on a uniformly sampled random variable will give us draws from the $PDF$ that we need. Given a continuous uniform variable $ U$ in $ [0,1)$ and an invertible cumulative distribution function $F_{X}$, the random variable $X=F_{X}^{-1}(U)$ has distribution $F_{X}$.[wikipedia page].
In this particular case the correct distributions are $PDF(r)=2r$ and $CDF=r^2$ and so $CDF^{-1}(x)=x^{\frac{1}{2}}$. So, for the radius we draw uniform samples again but this time instead of using it directly, we take square-root. Alternative way to think about this is to imagine that we are taking uniform samples for $r^2$ from $[0,r_0^2]$.
r = r0*np.sqrt(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))
ax.add_artist(plt.Circle((0, 0), r0,alpha=0.25))
ax.plot(x,y,'ro', alpha=0.5, mec='none')
ax.set_aspect('equal')
This looks uniform now. In my opinion, this is the most elegant method. It needed a total of $2N$ samples. This is not necessarily the most efficient method. Rejection sampling is faster in low dimensions. For circles rejection method was about 3.5-4 times faster than this method.
There is an alternative way to sample the "linear" probability distribution, without going the inverse transform sampling route. This needs two independent draws from a uniform distribution for the radius (i.e. three draws total). First, we add the two draws. The resulting points look like they have been drawn from a distribution that looks like a hat function with a peak at $r_0$ and domain from $0$ to $2r_0$. 'Folding back' the portion in $[r_0,2r_0]$ onto $[0,r_0]$, we get our linearly increasing distribution. (Somewhat distantly related: Central limit theorem)
x1 = uniform(0.0,r0,N)
x2 = uniform(0.0,r0,N)
r = x1+x2
r = r0-np.abs(r-r0) #same as r=2*r0-r if r>r0 else r
########################################
x,y= r*np.cos(t),r*np.sin(t)
fig,ax=plt.subplots(figsize=(6,6))
ax.add_artist(plt.Circle((0, 0), r0,alpha=0.25))
ax.plot(x,y,'ro', alpha=0.5, mec='none')
ax.set_aspect('equal')
This looks alright as well. We needed to draw $3N$ samples.
First, uniform samples in the square region $[-r_0,r_0]\times[-r_0,r_0]$ are drawn. Then, only the samples that fall inside the circle are kept. This is repeated until we have the necessary number of samples.
x=[]
y=[]
r0square=r0**2
ns=0
k=0
while(ns<N):
k+=1
xp = uniform(-r0,r0,1)
yp = uniform(-r0,r0,1)
if (xp**2+yp**2) < r0square:#acceptable
x.append(xp)
y.append(yp)
ns+=1
########################################
fig,ax=plt.subplots(figsize=(6,6))
ax.add_artist(plt.Circle((0, 0), r0,alpha=0.25))
ax.plot(x,y,'ro', alpha=0.5, mec='none')
ax.set_aspect('equal')
print('Total number of samples drawn:',2*k)
Total number of samples drawn: 2552
Unlike the previous methods, number of draws needed for this method is not known apriori. This is very simple to program easily extendable to 3D and different shapes like ellipses and ellipsoids for instance.
This is a very elegant method described in Barthe et. al, 2005 and generalizes to arbitrary dimensions.
N=1000
u=standard_normal(N)
v=standard_normal(N)
e=-2*np.log(uniform(size=N))
pts=np.c_[u,v]/((u**2+v**2 +e)**0.5)[:,np.newaxis]
fig,ax=plt.subplots(figsize=(6,6))
ax.add_artist(plt.Circle((0, 0), 1.0,alpha=0.25))
_=ax.plot(pts[:,0],pts[:,1],'ro', alpha=0.5, mec='none')
This is another interesting general method. The first step is to get N uniformly random points on a 4-sphere. This can be done by drawing coordinates from a standard normal distribution and normalizing. Then, we have $(x_1,x_2,x_3,x_4)$ such that $x_1^2 + x_2^2 + x_3^2 + x_4^2 = 1$. Now, drop two the last two coordinates (any two really) from all the points. The resulting points $(x_1,x_2)$ are uniformly random inside a circle. The N points formed by the dropped coordinates are also uniformly random inside a circle, but they are not independent of the first N due to the normalization. I have plotted both sets of points below in different colors, but one only set should be used as mentioned.
N=500
u=standard_normal((N,4))
#normalize each row
u=u/np.linalg.norm(u,axis=1)[:,np.newaxis]
fig,ax=plt.subplots(figsize=(6,6))
_=ax.plot(u[:,0],u[:,1],'ro', alpha=0.5, mec='none')
#dont use these points. these are not independent of the previous points!!!
_=ax.plot(u[:,2],u[:,3],'go', alpha=0.5, mec='none')
_=ax.add_artist(plt.Circle((0, 0), 1.0,alpha=0.25))
The math used here leads to apparently weird results. Consider points on the surface of a sphere given by $(x_i,y_i,z_i)$ where $x_i^2+y_i^2+z_i^2=1$. According to the rule used, each of $x_i$, $y_i$ and $z_i$ are uniform over $(-1,1)$. This is a bit non-intuitive. If we had drawn $x_i$, $y_i$ and $z_i$ uniformly from $(-1,1)$ then we would have uniformly sampled inside the cube $(-1,1)\times(-1,1)\times(-1,1)$. The discrepancy is due to $x,y,z$ being non-independent when coming from the surface of a sphere.
Sampling 'on a circle' is trivial. We can draw from a uniform distribution in $[0,2\pi)$. This is the same as generating random directions or randomly oriented vectors in 2D world. The same in 3D is less clear. Generating a uniform sampling on the surface of a sphere is the same as generating uniformly randomly oriented vectors in 3D space. (Note: This is different from randomly rotating a 3D object. A vector does not change when rotated about itself, whereas a general 3D object would). Again, I start with a incorrect naive implementation based on the spherical coordinate parametrization and then show methods that generate the correct distribution.
The surface of a sphere can be parametrized in $\theta$ and $\phi$. Lets see what we get when we naively sample these two from a uniform distribution. We'll take it to be a unit sphere.
N=800
phi = np.random.uniform(0,2*np.pi,size=N)
theta = np.random.uniform(0,2*np.pi,size=N)
st=np.sin(theta)
ct=np.cos(theta)
sp=np.sin(phi)
cp=np.cos(phi)
v=np.c_[st*cp,st*sp,ct]
fig,ax=plt.subplots(figsize=(10,10))
ax = plt.axes(projection='3d')
ax.plot(v[:,0],v[:,1],v[:,2],'ro', alpha=0.5, mec='none')
# ax.set_aspect('equal')
# Make spherical surface
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
# Plot the surface
collection=ax.plot_surface(x, y, z,alpha=0.25)
It is perhaps a little hard to make out but you should see points clustering at the poles. There are more points at the north and south poles that at the equator. This is similar to the problem seen in the earlier section with the naive sampling in the circle. To address this issue we can again do a inverse transform sampling of $\theta$ as shown in the next section. Notice that there is no problem in the $\phi$ direction. So $\phi$ was okay to be drawn from a uniform distribution.
When a hollow sphere is sliced along the latitudes(i.e. with equidistant parallel planes), the surface area of the slice (referring to surface area of a sliced hollow sphere and not the "cross-section" area of a sliced solid sphere) does not depend on the latitude, i.e. the $z$ coordinate of uniformly random points on a sphere are uniformly drawn from $(-1,1)$. Applying the arc cosine function to uniformly drawn samples from $(-1,1)$ to give the right spread for $\theta$. $\phi\in[0,2\pi)$ can be picked uniformly.
N=800
phi = uniform(0,2*np.pi,size=N)
theta = np.arccos(uniform(-1,1,size=N))
st=np.sin(theta)
ct=np.cos(theta)
sp=np.sin(phi)
cp=np.cos(phi)
pts=np.c_[st*cp,st*sp,ct]
fig,ax=plt.subplots(figsize=(10,10))
ax = plt.axes(projection='3d')
ax.plot(pts[:,0],pts[:,1],pts[:,2],'ro')
# For Visualization, make spherical surface
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
# Plot the surface
collection=ax.plot_surface(x, y, z,alpha=0.25)
This is an intriguing method and its not immediately apparent why it works. I will first demonstrate it in 2D i.e. sampling on a circle for illustration (of course, sampling on a circle is just uniformly sampling from $[0,2\pi)$).
For sampling on a circle, we generate (x,y) pairs independently from a standard normal distribution. Lets first see what that gives us.
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111, aspect='equal', autoscale_on=True)
N=1000
pts=standard_normal(size=[N,2])
ax.add_artist(plt.Circle((0, 0), 1.0,fc='none', ec='g',lw=1.0 ))
_=ax.plot(pts[:,0],pts[:,1],'ro', alpha=0.5, mec='none')
We notice that the distribution of the points is isotropic. It exhibits radial symmetry. This phenomenon makes the sampling procedure possible. Since they are radially symmetric, normalizing the points will cause them to lie on the boundary of the circle and be uniformly randomly distrubuted, as shown by the next code block.
pts/=np.sqrt(np.sum(pts**2,axis=-1)[:,np.newaxis])
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111, aspect='equal', autoscale_on=True)
ax.add_artist(plt.Circle((0, 0), 1.0,alpha=0.8, fc='none', ec='g', lw=1.0))
_=ax.plot(pts[:,0],pts[:,1],'ro', alpha=0.5, mec='none')
The extension to sampling on spherical surfaces is routine. In fact, this method generalizes to arbitrary dimensions. It is perhaps not the most efficient method in low dimensions but I like the generality.
For surface of sphere, we first get points (in 3D space) by independently drawing the x,y,z coordinates for each point from a normal distribution. The resulting points will be spherically symmetric. We then normalize the position vectors so that the points lie on a sphere. The resulting distribution will be uniform on a sphere.
N=2000
pts=standard_normal(size=[N,3]) #get 3 N-length values
pts/=np.sqrt(np.sum(pts**2,axis=-1)[:,np.newaxis]) #normalize: v=v/|v|
fig,ax=plt.subplots(figsize=(10,10))
ax = plt.axes(projection='3d')
ax.plot(pts[:,0],pts[:,1],pts[:,2],'ro', mec='none')
# For Visualization, make spherical surface
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
# Plot the surface
collection=ax.plot_surface(x, y, z,alpha=0.25)
Choosing a point from the surface of sphere, George Marsaglia 1972, The Annals of Mathematical Statistics, vol. 43, No. 2, 645 - 646.
From the paper:
Generate $V_1$ and $V_2$ independent uniform on (-1,1) until $S= V_1^2+ V_2^2<1$ then form $(~2V_1(1-S)^{1/2}~,~2V_2(1-S)^{1/2}~,~1-2S~)$
Essentially we first generate uniform random samples on disk($\equiv$ inside a circle) and then transform them as prescribed to get points on a sphere. To generate the points on a disk, instead of rejection sampling as described, we will use the inverese transform sampling method described earlier.
N=500
r = np.sqrt(uniform(0.0,1.0,N))
t = uniform(0,2*np.pi,N)
v1,v2 = r*np.cos(t),r*np.sin(t)
s = v1**2+v2**2
pts = np.c_[2*v1*(np.sqrt(1-s)),2*v2*(np.sqrt(1-s)),1-2*s]
fig,ax=plt.subplots(figsize=(10,10))
ax = plt.axes(projection='3d')
ax.plot(pts[:,0],pts[:,1],pts[:,2],'ro')
# For visualization, make spherical surface
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
# Plot the surface
collection=ax.plot_surface(x, y, z,alpha=0.25)
A Quaternion is a set of 4 numbers represented as $(w,x,y,z)$ or $w+x~i+y~j+z~k$. This is an extension of complex numbers to 3D space. $i,j,k$ are imaginary numbers such that $i^2=j^2=k^2=ijk=-1$. The $i,j,k$ do not commute($ijk\neq ikj$). Quaternions can be used to rotate objects in 3D space. A vector (==axi-symmetric object) in 3D space can be rotated to a general position by rotating it about each of three axis in succession. This operation does not commute either. However, it turns out that any rotated position can be obtained by rotating the original vector by a particular angle about a particular axis. This is the axis-angle specification of a rotation. Quaternion is an alternative to the axis-angle method to represent rotations.
Choose three points $u_1,u_2,u_3 \in [0,1]$ uniformly at random. A uniform, random quaternion is given by the simple expression $h = (\sqrt{1-u_1}\sin 2 \pi u_2, \sqrt{1-u_1}\cos 2 \pi u_2, \sqrt{u_1}\sin 2 \pi u_3, \sqrt{u_1}\cos 2 \pi u_3)$. (K. Shoemake. Uniform random rotations. D. Kirk, editor, Graphics Gems III, pages 124-132. Academic, New York, 1992)
u1 = uniform(0.0,1.0,N)
u2 = uniform(0.0,1.0,N)
u3 = uniform(0.0,1.0,N)
rnd_quats=np.c_[np.sqrt(1-u1)*np.sin(2*np.pi*u2), np.sqrt(1-u1)*np.cos(2*np.pi*u2),
np.sqrt( u1 )*np.sin(2*np.pi*u3), np.sqrt( u1 )*np.cos(2*np.pi*u3) ]
Now that we have the random quaternions, we can generate the needed points by applying the quaternions to a reference position vector. I will choose $[1,0,0]$. To maintain the uniformity of the distribution each quaternion has to be applied to the same reference vector, which in this case is $[1,0,0]$. The actual rotation is done by multiplying the vector with a transformation matrix that is generated from the quaternion. The function below gives us the rotation matrix for a given quaternion.
def rotation_matrix_q(q):
"""Returns rotation matrix from quaternion.
"""
q *= np.sqrt(2.0 / (q @ q))
q = np.outer(q, q)
return np.asarray([
[1.0-q[2, 2]-q[3, 3], q[1, 2]-q[3, 0], q[1, 3]+q[2, 0]],
[ q[1, 2]+q[3, 0], 1.0-q[1, 1]-q[3, 3], q[2, 3]-q[1, 0]],
[ q[1, 3]-q[2, 0], q[2, 3]+q[1, 0], 1.0-q[1, 1]-q[2, 2]],])
Now we can rotate the the reference vector by applying the rotation matrix on the vector as shown below
ref_v=np.asarray([1.0,0.0,0.0])
pts=np.empty([N,3])
for i,q in enumerate(rnd_quats):
rot_mat = rotation_matrix_q(q)
pts[i] = rot_mat.dot(ref_v.T).T
fig,ax=plt.subplots(figsize=(10,10))
ax = plt.axes(projection='3d')
ax.plot(pts[:,0],pts[:,1],pts[:,2],'ro')
# For Visualization, make spherical surface
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
# Plot the surface
collection=ax.plot_surface(x, y, z,alpha=0.25)
The following function calculates the rotation matrix from an axis-angle specification of a rotational transformation
def rotation_matrix_aa(axis,angle):
"""Returns matrix for rotation by 'angle' radians about 'axis' through the origin
"""
axis/=np.sqrt(np.sum(axis**2))
l,m,n = axis
sina,cosa=np.sin(angle),np.cos(angle)
return np.asarray([
[ l*l*(1-cosa)+ cosa , l*m*(1-cosa)+n*sina , l*n*(1-cosa)-m*sina ],
[ m*l*(1-cosa)-n*sina , m*m*(1-cosa)+ cosa , m*n*(1-cosa)+l*sina ],
[ n*l*(1-cosa)+m*sina , n*m*(1-cosa)-l*sina , n*n*(1-cosa)+ cosa ]])
This function is tested below. Say, we want to get the transformation matrix that rotates one vector v1 to another vector v2. Two of the randomly generated vectors earlier are used for the demonstration (pts[4] and pts[5]). We rotate v1 to v2 and check if the rotation worked properly.
v1=pts[4]
v2=pts[5]
angle=np.arccos(np.dot(v1,v2)) # these are unit vectors so dot product is cosine of angle
axis=np.cross(v2,v1) # order matters
rot_mat=rotation_matrix_aa(axis,angle)
v1_rot=rot_mat.dot(v1.T).T
print(f"Are v1 and v2 equal? : {np.allclose(v1_rot,v2)}")
Are v1 and v2 equal? : True