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?

@Veedrac Can't you define a plane that goes through the origin by a single normal vector? See http://math.stackexchange.com/questions/952525/select-a-random-hyperplane . Also the method I am using is meant to be the one described in "Another easy way to pick a random point[..]" at http://mathworld.wolfram.com/SpherePointPicking.html — Simd, Oct 01 '14 at 19:36

score 4 · Answer 1 · answered Oct 01 '14 at 20:30

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

A sphere of points around the normal

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

A satellite dish... sort of.

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

It's a sphere made of spherically bound planes!

score 4 · Accepted Answer · edited May 23 '17 at 10:33

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

enter image description here

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

enter image description here

I really like that you linked to pathpatch_2d_to_3d , thank you! — Simd, Oct 02 '14 at 09:28

Jendas · Answer 3 · 2014-10-01T20:50:10.650

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...

edited Oct 01 '14 at 20:50

answered Oct 01 '14 at 19:24

Jendas

3,359
3
27
55

The method I am using to pick the random point on the sphere is at the end of http://mathworld.wolfram.com/SpherePointPicking.html . This then defines the plane as in http://math.stackexchange.com/questions/952525/select-a-random-hyperplane Clearly I am doing something wrong, however. – Simd Oct 01 '14 at 19:39
Well it makes sense, but to tell the truth, I do not really see how it becomes uniform distribution... – Jendas Oct 01 '14 at 19:44
@Jendas Haha! It looks like you got the answer in an hour before me, not 3 minutes! Well, the resulting output from this is [a bit of a mess.](http://imgur.com/G8JyArX) Is this what you get too? – Veedrac Oct 01 '14 at 20:36
Yea, I just kept editing the deleted answer and finally got here. I also get the mess, but I focused on keeping the original code untouched as much as possible for asker to better understand... Plus he questioned the uniform distribution of planes and that is what mattered. – Jendas Oct 01 '14 at 20:41

score 1 · Answer 4 · answered Oct 01 '14 at 20:23

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)

How to draw random planes

4 Answers4