[Numpy-discussion] Even Sphere Volume

[Chris Colbert]

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