How to draw random planes

Question

I am using the following code to draw random planes in 3d that go through the origin.

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#Number of hyperplanes
n = 20
#Dimension of space
d = 3

plt3d = plt.figure().gca(projection='3d')
for i in xrange(n):
    #Create random point on unit sphere
    v = np.random.normal(size = d)
    v = v/np.sqrt(np.sum(v**2))
    # create x,y
    xx, yy = np.meshgrid(range(-5,5), range(-5,5))
    z = (-v[0] * xx - v[1] * yy)/v[2]
    # plot the surface
    plt3d.plot_surface(xx, yy, z, alpha = 0.5)
plt.show()

But looking at the picture I don't believe they are uniformly chosen. What am I doing wrong?

Answer 1

Your code is generating planes with randomly distributed normals. They just don't look that way because the z-scale is so much larger than the x- and y-scales.

You can generate a better looking image by generating points which are evenly distributed on the plane. To do that, parametrize the plane in terms of new coordinates (u, v), and then sample the plane on a uniformly spaced grid of (u,v) points. Then convert those (u,v) points into points in (x,y,z)-space.

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import math
import itertools as IT

def points_on_sphere(dim, N, norm=np.random.normal):
    """
    http://en.wikipedia.org/wiki/N-sphere#Generating_random_points
    """
    normal_deviates = norm(size=(N, dim))
    radius = np.sqrt((normal_deviates ** 2).sum(axis=0))
    points = normal_deviates / radius
    return points

# Number of hyperplanes
n = 10
# Dimension of space
d = 3

fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))
points = points_on_sphere(n, d).T
uu, vv = np.meshgrid([-5, 5], [-5, 5], sparse=True)
colors = np.linspace(0, 1, len(points))
cmap = plt.get_cmap('jet')
for nhat, c in IT.izip(points, colors):
    u = (0, 1, 0) if np.allclose(nhat, (1, 0, 0)) else np.cross(nhat, (1, 0, 0))
    u /= math.sqrt((u ** 2).sum())
    v = np.cross(nhat, u)
    u = u[:, np.newaxis, np.newaxis]
    v = v[:, np.newaxis, np.newaxis]
    xx, yy, zz = u * uu + v * vv
    ax.plot_surface(xx, yy, zz, alpha=0.5, color=cmap(c))
ax.set_xlim3d([-5,5])
ax.set_ylim3d([-5,5])
ax.set_zlim3d([-5,5])        
plt.show()

Alternatively, you could avoid the hairy math by using Till Hoffmann's pathpatch_2d_to_3d utility function:

for nhat, c in IT.izip(points, colors):
    p = patches.Rectangle((-2.5, -2.5), 5, 5, color=cmap(c), alpha=0.5)
    ax.add_patch(p)
    pathpatch_2d_to_3d(p, z=0, normal=nhat)

ax.set_xlim3d([-5,5])
ax.set_ylim3d([-5,5])
ax.set_zlim3d([-5,5])        
plt.show()

Answer 2

I suggest you check your axes. Your calculation makes the Z axis way too large, which means that you have an absurdly biased point of view.

First check that your normals are evenly distributed on the circle:

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#Number of hyperplanes
n = 1000
#Dimension of space
d = 3

plt3d = plt.figure().gca(projection='3d')
for i in xrange(n):
    #Create random point on unit sphere
    v = np.random.normal(size = d)
    v = v/np.sqrt(np.sum(v**2))
    v *= 10

    plt3d.scatter(v[0], v[1], v[2])

plt3d.set_aspect(1)
plt3d.set_xlim(-10, 10)
plt3d.set_ylim(-10, 10)
plt3d.set_zlim(-10, 10)

plt.show()

Then check that your plane is being created correctly:

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#Number of hyperplanes
n = 1
#Dimension of space
d = 3

plt3d = plt.figure().gca(projection='3d')
for i in xrange(n):
    #Create random point on unit sphere
    v = np.random.normal(size = d)
    v = v/np.sqrt(np.sum(v**2))
    v *= 10

    # create x,y
    xx, yy = np.meshgrid(np.arange(-5,5,0.3), np.arange(-5,5,0.3))
    xx = xx.flatten()
    yy = yy.flatten()
    z = (-v[0] * xx - v[1] * yy)/v[2]

    # Hack to keep the plane small
    filter = xx**2 + yy**2 + z**2 < 5**2
    xx = xx[filter]
    yy = yy[filter]
    z = z[filter]

    # plot the surface
    plt3d.scatter(xx, yy, z, alpha = 0.5)

    for i in np.arange(0.1, 1, 0.1):
        plt3d.scatter(i*v[0], i*v[1], i*v[2])

plt3d.set_aspect(1)
plt3d.set_xlim(-10, 10)
plt3d.set_ylim(-10, 10)
plt3d.set_zlim(-10, 10)

plt.show()

Then you can see that you've actually got good results!

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#Number of hyperplanes
n = 100
#Dimension of space
d = 3

plt3d = plt.figure().gca(projection='3d')
for i in xrange(n):
    #Create random point on unit sphere
    v = np.random.normal(size = d)
    v = v/np.sqrt(np.sum(v**2))
    v *= 10

    # create x,y
    xx, yy = np.meshgrid(np.arange(-5,5,0.3), np.arange(-5,5,0.3))
    xx = xx.flatten()
    yy = yy.flatten()
    z = (-v[0] * xx - v[1] * yy)/v[2]

    # Hack to keep the plane small
    filter = xx**2 + yy**2 + z**2 < 5**2
    xx = xx[filter]
    yy = yy[filter]
    z = z[filter]

    # plot the surface
    plt3d.scatter(xx, yy, z, alpha = 0.5)

plt3d.set_aspect(1)
plt3d.set_xlim(-10, 10)
plt3d.set_ylim(-10, 10)
plt3d.set_zlim(-10, 10)

plt.show()

Answer 3

Looking isn't everything. You better measure next time :-]. It seems to be non-randomly distributed because you did not fixed the axis. Therefore you see one major plane, that went up to the sky and the rest, because of the scale, look very similar and not randomly distributed.

How about this code:

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#Number of hyperplanes
n = 20
#Dimension of space
d = 3

plt3d = plt.figure().gca(projection='3d')
for i in xrange(n):
    #Create random point on unit sphere
    v = np.random.normal(size = d)
    v = v/np.sqrt(np.sum(v**2))
    # create x,y
    xx, yy = np.meshgrid(range(-1,1), range(-1,1))
    z = (-v[0] * xx - v[1] * yy)/v[2]
    # plot the surface
    plt3d.plot_surface(xx, yy, z, alpha = 0.5)

plt3d.set_xlim3d([-1,1])
plt3d.set_ylim3d([-1,1])
plt3d.set_zlim3d([-1,1])
plt.show()

It isn't perfect, but it seems much more random now...

Answer 4

I tried this one, maybe this is a better way to create uniform planes. I randomly take two different angles for spherical coordinate system and convert it to cartesian coordinates to get the normal vector of the plane. Also when you plot, you should be aware that middle point of your plane is not on the origin.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = Axes3D(fig)

for i in range(20):
    theta = 2*np.pi*np.random.uniform(-1,1)    
    psi = 2*np.pi*np.random.uniform(-1,1)
    normal = np.array([np.sin(theta)*np.cos(psi),np.sin(theta)*np.sin(psi),
                       np.cos(theta)])
    xx, yy = np.meshgrid(np.arange(-1,1), np.arange(-1,1))
    z = (-normal[0] * xx - normal[1] * yy)/normal[2]
    ax.plot_surface(xx, yy, z, alpha=0.5)