[Numpy-discussion] Even Sphere Volume
Chris Colbert
sccolbert@gmail....
Tue Jul 7 20:54:45 CDT 2009
I figured I'd contribute something similar: the "even" distribution of
points on a sphere.
I can't take credit for the algorithm though, I got it from the following
page and just vectorized it, and tweaked it so it uses the golden ration
properly:
http://www.xsi-blog.com/archives/115
It makes use of the the golden ratio to increment the longitude thus tracing
a spiral out around the sphere.
############################################################
# A simple program to place evenly distributed points on a unit sphere
import numpy as np
from enthought.mayavi import mlab
def pointsOnSphere(N):
""" Generate even distributed N points on the unit sphere
centered at the origin. Uses the 'Golded Spiral'"""
x = np.array((N,), dtype='float')
y = np.array((N,), dtype='float')
z = np.array((N,), dtype='float')
phi = (1 + np.sqrt(5)) / 2 # the golden ratio
long_incr = 2*np.pi / phi # how much to increment the longitude
dz = 2.0 / float(N) # a unit sphere has diameter 2
bands = np.arange(0, N, 1) # each band will have one point placed on it
z = bands * dz - 1 + (dz/2) # the z location of each band/point
r = np.sqrt(1 - z*z) # the radius can be directly determined from
height
az = bands * long_incr # the azimuth where to place the point
x = r * np.cos(az)
y = r * np.sin(az)
return (x, y, z)
if __name__=='__main__':
x, y, z = pointsOnSphere(3000)
pts = mlab.points3d(x, y, z)
here's a picture of the result:
www.therealstevencolbert.com/dump/spherepts.png
Cheers,
Chris
On Mon, Jul 6, 2009 at 11:27 PM, Ian Mallett <geometrian@gmail.com> wrote:
> That didn't fix it. I messed around some more, but I couldn't get the
> spherical coordinates working. I decided to rework my first method. By
> raising the radius to the one third power, like for the other method,
> basically the same thing is accomplished. It's working now, thanks. :D
>
> vecs = numpy.random.standard_normal((size,size,3))
> magnitudes = numpy.sqrt((vecs*vecs).sum(axis=-1))
> vecs = vecs / magnitudes[...,numpy.newaxis]
> randlen = numpy.random.random((size,size))
> randlen = randlen ** (1.0/3.0)
> randpoints = vecs*randlen[...,numpy.newaxis]
> rgb = ((randpoints+1.0)/2.0)*255.0
>
> Ian
>
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/numpy-discussion/attachments/20090707/190cc767/attachment.html
More information about the NumPy-Discussion
mailing list