By Befall


2012-03-07 11:39:06 8 Comments

I need an algorithm that can give me positions around a sphere for N points (less than 20, probably) that vaguely spreads them out. There's no need for "perfection", but I just need it so none of them are bunched together.

  • This question provided good code, but I couldn't find a way to make this uniform, as this seemed 100% randomized.
  • This blog post recommended had two ways allowing input of number of points on the sphere, but the Saff and Kuijlaars algorithm is exactly in psuedocode I could transcribe, and the code example I found contained "node[k]", which I couldn't see explained and ruined that possibility. The second blog example was the Golden Section Spiral, which gave me strange, bunched up results, with no clear way to define a constant radius.
  • This algorithm from this question seems like it could possibly work, but I can't piece together what's on that page into psuedocode or anything.

A few other question threads I came across spoke of randomized uniform distribution, which adds a level of complexity I'm not concerned about. I apologize that this is such a silly question, but I wanted to show that I've truly looked hard and still come up short.

So, what I'm looking for is simple pseudocode to evenly distribute N points around a unit sphere, that either returns in spherical or Cartesian coordinates. Even better if it can even distribute with a bit of randomization (think planets around a star, decently spread out, but with room for leeway).

15 comments

@CR Drost 2017-05-24 16:36:18

The golden spiral method

You said you couldn’t get the golden spiral method to work and that’s a shame because it’s really, really good. I would like to give you a complete understanding of it so that maybe you can understand how to keep this away from being “bunched up.”

So here’s a fast, non-random way to create a lattice that is approximately correct; as discussed above, no lattice will be perfect, but this may be good enough. It is compared to other methods e.g. at BendWavy.org but it just has a nice and pretty look as well as a guarantee about even spacing in the limit.

Primer: sunflower spirals on the unit disk

To understand this algorithm, I first invite you to look at the 2D sunflower spiral algorithm. This is based on the fact that the most irrational number is the golden ratio (1 + sqrt(5))/2 and if one emits points by the approach “stand at the center, turn a golden ratio of whole turns, then emit another point in that direction,” one naturally constructs a spiral which, as you get to higher and higher numbers of points, nevertheless refuses to have well-defined ‘bars’ that the points line up on.(Note 1.)

The algorithm for even spacing on a disk is,

from numpy import pi, cos, sin, sqrt, arange
import matplotlib.pyplot as pp

num_pts = 100
indices = arange(0, num_pts, dtype=float) + 0.5

r = sqrt(indices/num_pts)
theta = pi * (1 + 5**0.5) * indices

pp.scatter(r*cos(theta), r*sin(theta))
pp.show()

and it produces results that look like (n=100 and n=1000):

enter image description here

Spacing the points radially

The key strange thing is the formula r = sqrt(indices / num_pts); how did I come to that one? (Note 2.)

Well, I am using the square root here because I want these to have even-area spacing around the disk. That is the same as saying that in the limit of large N I want a little region R ∈ (r, r + dr), Θ ∈ (θ, θ + dθ) to contain a number of points proportional to its area, which is r dr dθ. Now if we pretend that we are talking about a random variable here, this has a straightforward interpretation as saying that the joint probability density for (R, Θ) is just c r for some constant c. Normalization on the unit disk would then force c = 1/π.

Now let me introduce a trick. It comes from probability theory where it’s known as sampling the inverse CDF: suppose you wanted to generate a random variable with a probability density f(z) and you have a random variable U ~ Uniform(0, 1), just like comes out of random() in most programming languages. How do you do this?

  1. First, turn your density into a cumulative distribution function or CDF, which we will call F(z). A CDF, remember, increases monotonically from 0 to 1 with derivative f(z).
  2. Then calculate the CDF’s inverse function F-1(z).
  3. You will find that Z = F-1(U) is distributed according to the target density. (Note 3).

Now the golden-ratio spiral trick spaces the points out in a nicely even pattern for θ so let’s integrate that out; for the unit disk we are left with F(r) = r2. So the inverse function is F-1(u) = u1/2, and therefore we would generate random points on the disk in polar coordinates with r = sqrt(random()); theta = 2 * pi * random().

Now instead of randomly sampling this inverse function we’re uniformly sampling it, and the nice thing about uniform sampling is that our results about how points are spread out in the limit of large N will behave as if we had randomly sampled it. This combination is the trick. Instead of random() we use (arange(0, num_pts, dtype=float) + 0.5)/num_pts, so that, say, if we want to sample 10 points they are r = 0.05, 0.15, 0.25, ... 0.95. We uniformly sample r to get equal-area spacing, and we use the sunflower increment to avoid awful “bars” of points in the output.

Now doing the sunflower on a sphere

The changes that we need to make to dot the sphere with points merely involve switching out the polar coordinates for spherical coordinates. The radial coordinate of course doesn't enter into this because we're on a unit sphere. To keep things a little more consistent here, even though I was trained as a physicist I'll use mathematicians' coordinates where 0 ≤ φ ≤ π is latitude coming down from the pole and 0 ≤ θ ≤ 2π is longitude. So the difference from above is that we are basically replacing the variable r with φ.

Our area element, which was r dr dθ, now becomes the not-much-more-complicated sin(φ) dφ dθ. So our joint density for uniform spacing is sin(φ)/4π. Integrating out θ, we find f(φ) = sin(φ)/2, thus F(φ) = (1 − cos(φ))/2. Inverting this we can see that a uniform random variable would look like acos(1 - 2 u), but we sample uniformly instead of randomly, so we instead use φk = acos(1 − 2 (k + 0.5)/N). And the rest of the algorithm is just projecting this onto the x, y, and z coordinates:

from numpy import pi, cos, sin, arccos, arange
import mpl_toolkits.mplot3d
import matplotlib.pyplot as pp

num_pts = 1000
indices = arange(0, num_pts, dtype=float) + 0.5

phi = arccos(1 - 2*indices/num_pts)
theta = pi * (1 + 5**0.5) * indices

x, y, z = cos(theta) * sin(phi), sin(theta) * sin(phi), cos(phi);

pp.figure().add_subplot(111, projection='3d').scatter(x, y, z);
pp.show()

Again for n=100 and n=1000 the results look like: enter image description here enter image description here

Further research

I wanted to give a shout out to Martin Roberts’s blog. Note that above I created an offset of my indices by adding 0.5 to each index. This was just visually appealing to me, but it turns out that the choice of offset matters a lot and is not constant over the interval and can mean getting as much as 8% better accuracy in packing if chosen correctly. There should also be a way to get his R2 sequence to cover a sphere and it would be interesting to see if this also produced a nice even covering, perhaps as-is but perhaps needing to be, say, taken from only a half of the unit square cut diagonally or so and stretched around to get a circle.

Notes

  1. Those “bars” are formed by rational approximations to a number, and the best rational approximations to a number come from its continued fraction expression, z + 1/(n_1 + 1/(n_2 + 1/(n_3 + ...))) where z is an integer and n_1, n_2, n_3, ... is either a finite or infinite sequence of positive integers:

    def continued_fraction(r):
        while r != 0:
            n = floor(r)
            yield n
            r = 1/(r - n)
    

    Since the fraction part 1/(...) is always between zero and one, a large integer in the continued fraction allows for a particularly good rational approximation: “one divided by something between 100 and 101” is better than “one divided by something between 1 and 2.” The most irrational number is therefore the one which is 1 + 1/(1 + 1/(1 + ...)) and has no particularly good rational approximations; one can solve φ = 1 + 1/φ by multiplying through by φ to get the formula for the golden ratio.

  2. For folks who are not so familiar with NumPy -- all of the functions are “vectorized,” so that sqrt(array) is the same as what other languages might write map(sqrt, array). So this is a component-by-component sqrt application. The same also holds for division by a scalar or addition with scalars -- those apply to all components in parallel.

  3. The proof is simple once you know that this is the result. If you ask what's the probability that z < Z < z + dz, this is the same as asking what's the probability that z < F-1(U) < z + dz, apply F to all three expressions noting that it is a monotonically increasing function, hence F(z) < U < F(z + dz), expand the right hand side out to find F(z) + f(z) dz, and since U is uniform this probability is just f(z) dz as promised.

@whn 2017-11-30 15:46:35

I'm not sure why this is so far down, this is by far the best fast method to do this.

@CR Drost 2017-12-10 17:33:58

@snb thank you for the kind words! it is so far down in part because it is much, much younger than all the rest of the answers here. I am surprised that it is even doing as well as it has been.

@Felix D. 2018-05-29 22:43:26

One question that remains for me is: How many points n do i need to distribute for a given maximum distance between any two points?

@CR Drost 2018-06-01 12:39:56

@FelixD. That sounds like a question that could get very complicated very fast especially if you start using, say, great-circle distances rather than Euclidean distances. But maybe I can answer a simple question, if one converts the points on the sphere to their Voronoi diagram, one can describe each Voronoi cell as having approximately area 4π/N and one can convert this to a characteristic distance by pretending it's a circle rather than a rhombus, πr² = 4π/N. Then r=2/√(N).

@Felix D. 2018-06-01 21:09:15

Okay thanks, that sounds like a good estimate for sufficiently large N!

@Lanny 2018-10-25 03:09:05

Another good writeup on evenly distributing points on a sphere is available at extremelearning.com.au/evenly-distributing-points-on-a-spher‌​e

@dmckee --- ex-moderator kitten 2019-01-25 18:53:32

Using the sampling theorem with actually uniform instead of randomly-uniform input is one of those things that makes me say "Well, why the #$%& didn't I think of that?". Nice.

@CR Drost 2019-01-26 23:21:42

@dmckee thanks for saying so! Especially coming from someone I respect as much as you, that's a high compliment!

@Ruslan 2020-07-25 21:37:04

In Spacing the points radially section you say "sphere" twice. Did you instead mean "disk" there?

@Jonathan H 2020-09-02 17:21:12

Looks like the reference you cite divides by the golden ratio to obtain theta, but you seem to multiply instead. Is that a mistake?

@CR Drost 2020-09-03 07:27:46

Great question! I believe my answer is closer to the “reason it works” while Martin’s squeaks out an extra bit of precision. So the golden ratio by definition satisfies φ² = φ + 1, which rearranges to φ – 1 = 1/φ, multiplying by 2 π, that leading digit 1 just gets nuked by the trig functions. So in floating point, just subtracting the one would fill that 53rd bit with a 0 where a 1 would be more correct.

@Ismael Harun 2020-07-20 19:36:34

@robert king It's a really nice solution but has some sloppy bugs in it. I know it helped me a lot though, so never mind the sloppiness. :) Here is a cleaned up version....

from math import pi, asin, sin, degrees
halfpi, twopi = .5 * pi, 2 * pi
sphere_area = lambda R=1.0: 4 * pi * R ** 2

lat_dist = lambda lat, R=1.0: R*(1-sin(lat))

#A = 2*pi*R^2(1-sin(lat))
def sphere_latarea(lat, R=1.0):
    if -halfpi > lat or lat > halfpi:
        raise ValueError("lat must be between -halfpi and halfpi")
    return 2 * pi * R ** 2 * (1-sin(lat))

sphere_lonarea = lambda lon, R=1.0: \
        4 * pi * R ** 2 * lon / twopi

#A = 2*pi*R^2 |sin(lat1)-sin(lat2)| |lon1-lon2|/360
#    = (pi/180)R^2 |sin(lat1)-sin(lat2)| |lon1-lon2|
sphere_rectarea = lambda lat0, lat1, lon0, lon1, R=1.0: \
        (sphere_latarea(lat0, R)-sphere_latarea(lat1, R)) * (lon1-lon0) / twopi


def test_sphere(n_lats=10, n_lons=19, radius=540.0):
    total_area = 0.0
    for i_lons in range(n_lons):
        lon0 = twopi * float(i_lons) / n_lons
        lon1 = twopi * float(i_lons+1) / n_lons
        for i_lats in range(n_lats):
            lat0 = asin(2 * float(i_lats) / n_lats - 1)
            lat1 = asin(2 * float(i_lats+1)/n_lats - 1)
            area = sphere_rectarea(lat0, lat1, lon0, lon1, radius)
            print("{:} {:}: {:9.4f} to  {:9.4f}, {:9.4f} to  {:9.4f} => area {:10.4f}"
                    .format(i_lats, i_lons
                    , degrees(lat0), degrees(lat1)
                    , degrees(lon0), degrees(lon1)
                    , area))
            total_area += area
    print("total_area = {:10.4f} (difference of {:10.4f})"
            .format(total_area, abs(total_area) - sphere_area(radius)))

test_sphere()

@Fnord 2014-09-30 17:51:24

The Fibonacci sphere algorithm is great for this. It is fast and gives results that at a glance will easily fool the human eye. You can see an example done with processing which will show the result over time as points are added. Here's another great interactive example made by @gman. And here's a simple implementation in python.

import math


def fibonacci_sphere(samples=1):

    points = []
    phi = math.pi * (3. - math.sqrt(5.))  # golden angle in radians

    for i in range(samples):
        y = 1 - (i / float(samples - 1)) * 2  # y goes from 1 to -1
        radius = math.sqrt(1 - y * y)  # radius at y

        theta = phi * i  # golden angle increment

        x = math.cos(theta) * radius
        z = math.sin(theta) * radius

        points.append((x, y, z))

    return points

1000 samples gives you this:

enter image description here

@Andrew Staroscik 2014-11-17 12:13:29

a variable n is called when defining phi: phi = ((i + rnd) % n) * increment. Does n = samples?

@Fnord 2014-11-18 19:50:52

@AndrewStaroscik yes! When i first wrote the code i used "n" as a variable and changed the name later but didn't do due diligence. Thanks for catching that!

@naphier 2016-04-11 07:34:05

This is great, but how would you use this and input the desired distance between the points? For example, say I want to use this to pack sphere with radius of 2 around a sphere. What variable would I change to get a separation of 2 between the points? Thanks!

@AturSams 2016-11-02 15:23:10

@naphier You can do a binary search if I understand your question.

@Xarbrough 2017-03-07 11:14:57

Great! How do I scale the entire sphere when inputting a radius other than 1?

@Fnord 2017-03-07 17:00:21

@Xarbrough the code gives you points around a unit sphere, so just multiply each point by whatever scalar you want for radius.

@doublefelix 2017-03-20 22:20:53

If you are wondering whether to use randomize=True, I did some tests and found that while the minimum nearest neighbor distance (NND) was about the same whether randomize was on or not, the stdev of NNDs divided by their mean was significantly lower (better) when randomize = False. In addition, the "evenness" of the distribution gets better with a larger sample size when randomize = False.

@Amir 2018-02-10 01:54:33

@Fnord I wonder, how can one also make y totally random? Currently y is indifferent to the random seed

@judepereira 2018-02-13 04:45:02

If you're using this in context with a real world, where your camera is placed at 0,0,0, multiply the resulting output (x, y, z) with the radius of your sphere.

@Tintinabulator Zea 2019-01-17 08:56:54

how can you perform this on an uneven cylinder type shape? if you have the vertices of the geometry

@pikachuchameleon 2019-03-10 20:57:30

@Fnord: Can we do this for higher dimensions?

@Ferdinando Randisi 2019-11-16 17:31:26

Really cool!!! What tool did you use to generate that render?

@Demitri 2020-04-13 05:30:49

The function above internally calculates phi, then converts to and returns Cartesian coordinates. What is the definition of theta if I wanted to return spherical coordinates without going back and forth?

@BlueRaja - Danny Pflughoeft 2012-03-07 17:32:58

This is known as packing points on a sphere, and there is no (known) general, perfect solution. However, there are plenty of imperfect solutions. The three most popular seem to be:

  1. Create a simulation. Treat each point as an electron constrained to a sphere, then run a simulation for a certain number of steps. The electrons' repulsion will naturally tend the system to a more stable state, where the points are about as far away from each other as they can get.
  2. Hypercube rejection. This fancy-sounding method is actually really simple: you uniformly choose points (much more than n of them) inside of the cube surrounding the sphere, then reject the points outside of the sphere. Treat the remaining points as vectors, and normalize them. These are your "samples" - choose n of them using some method (randomly, greedy, etc).
  3. Spiral approximations. You trace a spiral around a sphere, and evenly-distribute the points around the spiral. Because of the mathematics involved, these are more complicated to understand than the simulation, but much faster (and probably involving less code). The most popular seems to be by Saff, et al.

A lot more information about this problem can be found here

@Befall 2012-03-07 21:08:47

I will be looking into the spiral tactic that andrew cooke posted below, however, could you please clarify the difference between what I want and what "uniform random distribution" is? Is that just 100% randomized placement of points on a sphere so that they are uniformly placed? Thanks for the help. :)

@BlueRaja - Danny Pflughoeft 2012-03-07 21:57:37

@Befall: "uniform random distribution" refers to the probability-distribution being uniform - it means, when choosing a random point on the sphere, every point has an equal likelihood of being chosen. It has nothing to do with the final spatial-distribution of the points, and thus has nothing to do with your question.

@Befall 2012-03-07 21:59:54

Ahhh, okay, thanks very much. Searching for my question lead to a ton of answers for both, and I couldn't really grasp which was pointless to me.

@AturSams 2016-11-09 09:15:22

To be clear, every point has zero probability of being chosen. The ratio of the probabilities that the point will belong to any two areas on the surface of the sphere, is equal to the ratio of the surfaces.

@AturSams 2016-11-09 09:17:31

The difference is of some importance because on a computer, when using floats of doubles, there is a finite number of points on the sphere, depending on your tolerance, and they don't (I think) have an equal chance of being picked.

@gota 2018-03-02 15:02:34

3. doesn't work for dimensions > 3... 2. works for all dimensions but the rejection rate becomes impossible for dimensions > 20... Does anyone know of an implementation of 1. in e.g. python?

@Felix D. 2018-05-29 21:36:55

The last link is now dead

@andrew cooke 2012-03-07 12:15:01

In this example code node[k] is just the kth node. You are generating an array N points and node[k] is the kth (from 0 to N-1). If that is all that is confusing you, hopefully you can use that now.

(in other words, k is an array of size N that is defined before the code fragment starts, and which contains a list of the points).

Alternatively, building on the other answer here (and using Python):

> cat ll.py
from math import asin
nx = 4; ny = 5
for x in range(nx):
    lon = 360 * ((x+0.5) / nx)
    for y in range(ny):                                                         
        midpt = (y+0.5) / ny                                                    
        lat = 180 * asin(2*((y+0.5)/ny-0.5))                                    
        print lon,lat                                                           
> python2.7 ll.py                                                      
45.0 -166.91313924                                                              
45.0 -74.0730322921                                                             
45.0 0.0                                                                        
45.0 74.0730322921                                                              
45.0 166.91313924                                                               
135.0 -166.91313924                                                             
135.0 -74.0730322921                                                            
135.0 0.0                                                                       
135.0 74.0730322921                                                             
135.0 166.91313924                                                              
225.0 -166.91313924                                                             
225.0 -74.0730322921                                                            
225.0 0.0                                                                       
225.0 74.0730322921                                                             
225.0 166.91313924
315.0 -166.91313924
315.0 -74.0730322921
315.0 0.0
315.0 74.0730322921
315.0 166.91313924

If you plot that, you'll see that the vertical spacing is larger near the poles so that each point is situated in about the same total area of space (near the poles there's less space "horizontally", so it gives more "vertically").

This isn't the same as all points having about the same distance to their neighbours (which is what I think your links are talking about), but it may be sufficient for what you want and improves on simply making a uniform lat/lon grid.

@robert king 2012-03-07 12:23:10

nice, it's good to see a mathematical solution. I was thinking of using a helix and arc length separation. I'm still not certain on how to get the optimal solution which is an interesting problem.

@andrew cooke 2012-03-07 12:51:48

did you see that i edited my answer to include an explanation of node[k] at the top? i think that may be all you need...

@Befall 2012-03-07 21:06:26

Wonderful, thanks for the explanation. I'll try it out later, as I haven't time currently, but thank you so much for helping me out. I'll let you know how it ends up working for my purposes. ^^

@Befall 2012-03-08 01:22:56

Using the Spiral method fits my needs perfectly, thanks so much for the help and clarification. :)

@Scheintod 2014-12-20 20:55:32

The link seems dead.

@Ismael Harun 2020-07-20 18:48:49

Your latitude conversion to degrees seem incorrect. Shouldn't you divide by pi also?

@Ismael Harun 2020-07-20 19:26:49

Also your latitude calculation ranges from 0 to 2 exclusively, that can't be right. Perhaps you meant to subtract 1 from that.

@Andrew Wagner 2017-06-16 23:12:47

Healpix solves a closely related problem (pixelating the sphere with equal area pixels):

http://healpix.sourceforge.net/

It's probably overkill, but maybe after looking at it you'll realize some of it's other nice properties are interesting to you. It's way more than just a function that outputs a point cloud.

I landed here trying to find it again; the name "healpix" doesn't exactly evoke spheres...

@Edward Doolittle 2015-04-16 01:09:00

What you are looking for is called a spherical covering. The spherical covering problem is very hard and solutions are unknown except for small numbers of points. One thing that is known for sure is that given n points on a sphere, there always exist two points of distance d = (4-csc^2(\pi n/6(n-2)))^(1/2) or closer.

If you want a probabilistic method for generating points uniformly distributed on a sphere, it's easy: generate points in space uniformly by Gaussian distribution (it's built into Java, not hard to find the code for other languages). So in 3-dimensional space, you need something like

Random r = new Random();
double[] p = { r.nextGaussian(), r.nextGaussian(), r.nextGaussian() };

Then project the point onto the sphere by normalizing its distance from the origin

double norm = Math.sqrt( (p[0])^2 + (p[1])^2 + (p[2])^2 ); 
double[] sphereRandomPoint = { p[0]/norm, p[1]/norm, p[2]/norm };

The Gaussian distribution in n dimensions is spherically symmetric so the projection onto the sphere is uniform.

Of course, there's no guarantee that the distance between any two points in a collection of uniformly generated points will be bounded below, so you can use rejection to enforce any such conditions that you might have: probably it's best to generate the whole collection and then reject the whole collection if necessary. (Or use "early rejection" to reject the whole collection you've generated so far; just don't keep some points and drop others.) You can use the formula for d given above, minus some slack, to determine the min distance between points below which you will reject a set of points. You'll have to calculate n choose 2 distances, and the probability of rejection will depend on the slack; it's hard to say how, so run a simulation to get a feel for the relevant statistics.

@dmckee --- ex-moderator kitten 2019-04-11 17:25:42

Upvoted for the minimum maximum distance expressions. Useful for putting limits on the number of points you want to use. A reference to an authoritative source for that would be nice, though.

@aliential 2013-12-31 05:04:45

Try:

function sphere ( N:float,k:int):Vector3 {
    var inc =  Mathf.PI  * (3 - Mathf.Sqrt(5));
    var off = 2 / N;
    var y = k * off - 1 + (off / 2);
    var r = Mathf.Sqrt(1 - y*y);
    var phi = k * inc;
    return Vector3((Mathf.Cos(phi)*r), y, Mathf.Sin(phi)*r); 
};

The above function should run in loop with N loop total and k loop current iteration.

It is based on a sunflower seeds pattern, except the sunflower seeds are curved around into a half dome, and again into a sphere.

Here is a picture, except I put the camera half way inside the sphere so it looks 2d instead of 3d because the camera is same distance from all points. http://3.bp.blogspot.com/-9lbPHLccQHA/USXf88_bvVI/AAAAAAAAADY/j7qhQsSZsA8/s640/sphere.jpg

@user19371 2014-02-22 03:26:41

OR... to place 20 points, compute the centers of the icosahedronal faces. For 12 points, find the vertices of the icosahedron. For 30 points, the mid point of the edges of the icosahedron. you can do the same thing with the tetrahedron, cube, dodecahedron and octahedrons: one set of points is on the vertices, another on the center of the face and another on the center of the edges. They cannot be mixed, however.

@The Guy with The Hat 2016-01-21 06:28:50

A good idea, but it only works for 4, 6, 8, 12, 20, 24, or 30 points.

@chessofnerd 2016-04-19 06:02:54

If you want to cheat, you can use the center of faces and verticies. They will not be equi-spaced but a decent approximation. This is nice because it it deterministic.

@Matt S. 2013-04-21 06:07:42

This answer is based on the same 'theory' that is outlined well by this answer

I'm adding this answer as:
-- None of the other options fit the 'uniformity' need 'spot-on' (or not obviously-clearly so). (Noting to get the planet like distribution looking behavior particurally wanted in the original ask, you just reject from the finite list of the k uniformly created points at random (random wrt the index count in the k items back).)
--The closest other impl forced you to decide the 'N' by 'angular axis', vs. just 'one value of N' across both angular axis values ( which at low counts of N is very tricky to know what may, or may not matter (e.g. you want '5' points -- have fun ) )
--Furthermore, it's very hard to 'grok' how to differentiate between the other options without any imagery, so here's what this option looks like (below), and the ready-to-run implementation that goes with it.

with N at 20:

enter image description here
and then N at 80: enter image description here


here's the ready-to-run python3 code, where the emulation is that same source: " http://web.archive.org/web/20120421191837/http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere " found by others. ( The plotting I've included, that fires when run as 'main,' is taken from: http://www.scipy.org/Cookbook/Matplotlib/mplot3D )

from math import cos, sin, pi, sqrt

def GetPointsEquiAngularlyDistancedOnSphere(numberOfPoints=45):
    """ each point you get will be of form 'x, y, z'; in cartesian coordinates
        eg. the 'l2 distance' from the origion [0., 0., 0.] for each point will be 1.0 
        ------------
        converted from:  http://web.archive.org/web/20120421191837/http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere ) 
    """
    dlong = pi*(3.0-sqrt(5.0))  # ~2.39996323 
    dz   =  2.0/numberOfPoints
    long =  0.0
    z    =  1.0 - dz/2.0
    ptsOnSphere =[]
    for k in range( 0, numberOfPoints): 
        r    = sqrt(1.0-z*z)
        ptNew = (cos(long)*r, sin(long)*r, z)
        ptsOnSphere.append( ptNew )
        z    = z - dz
        long = long + dlong
    return ptsOnSphere

if __name__ == '__main__':                
    ptsOnSphere = GetPointsEquiAngularlyDistancedOnSphere( 80)    

    #toggle True/False to print them
    if( True ):    
        for pt in ptsOnSphere:  print( pt)

    #toggle True/False to plot them
    if(True):
        from numpy import *
        import pylab as p
        import mpl_toolkits.mplot3d.axes3d as p3

        fig=p.figure()
        ax = p3.Axes3D(fig)

        x_s=[];y_s=[]; z_s=[]

        for pt in ptsOnSphere:
            x_s.append( pt[0]); y_s.append( pt[1]); z_s.append( pt[2])

        ax.scatter3D( array( x_s), array( y_s), array( z_s) )                
        ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('Z')
        p.show()
        #end

tested at low counts (N in 2, 5, 7, 13, etc) and seems to work 'nice'

@Arthur Flower 2013-01-04 00:48:26

This works and it's deadly simple. As many points as you want:

    private function moveTweets():void {


        var newScale:Number=Scale(meshes.length,50,500,6,2);
        trace("new scale:"+newScale);


        var l:Number=this.meshes.length;
        var tweetMeshInstance:TweetMesh;
        var destx:Number;
        var desty:Number;
        var destz:Number;
        for (var i:Number=0;i<this.meshes.length;i++){

            tweetMeshInstance=meshes[i];

            var phi:Number = Math.acos( -1 + ( 2 * i ) / l );
            var theta:Number = Math.sqrt( l * Math.PI ) * phi;

            tweetMeshInstance.origX = (sphereRadius+5) * Math.cos( theta ) * Math.sin( phi );
            tweetMeshInstance.origY= (sphereRadius+5) * Math.sin( theta ) * Math.sin( phi );
            tweetMeshInstance.origZ = (sphereRadius+5) * Math.cos( phi );

            destx=sphereRadius * Math.cos( theta ) * Math.sin( phi );
            desty=sphereRadius * Math.sin( theta ) * Math.sin( phi );
            destz=sphereRadius * Math.cos( phi );

            tweetMeshInstance.lookAt(new Vector3D());


            TweenMax.to(tweetMeshInstance, 1, {scaleX:newScale,scaleY:newScale,x:destx,y:desty,z:destz,onUpdate:onLookAtTween, onUpdateParams:[tweetMeshInstance]});

        }

    }
    private function onLookAtTween(theMesh:TweetMesh):void {
        theMesh.lookAt(new Vector3D());
    }

@ksmith 2012-09-25 18:09:57

# create uniform spiral grid
numOfPoints = varargin[0]
vxyz = zeros((numOfPoints,3),dtype=float)
sq0 = 0.00033333333**2
sq2 = 0.9999998**2
sumsq = 2*sq0 + sq2
vxyz[numOfPoints -1] = array([(sqrt(sq0/sumsq)), 
                              (sqrt(sq0/sumsq)), 
                              (-sqrt(sq2/sumsq))])
vxyz[0] = -vxyz[numOfPoints -1] 
phi2 = sqrt(5)*0.5 + 2.5
rootCnt = sqrt(numOfPoints)
prevLongitude = 0
for index in arange(1, (numOfPoints -1), 1, dtype=float):
  zInc = (2*index)/(numOfPoints) -1
  radius = sqrt(1-zInc**2)

  longitude = phi2/(rootCnt*radius)
  longitude = longitude + prevLongitude
  while (longitude > 2*pi): 
    longitude = longitude - 2*pi

  prevLongitude = longitude
  if (longitude > pi):
    longitude = longitude - 2*pi

  latitude = arccos(zInc) - pi/2
  vxyz[index] = array([ (cos(latitude) * cos(longitude)) ,
                        (cos(latitude) * sin(longitude)), 
                        sin(latitude)])

@hcarver 2012-09-26 06:53:59

It'd be helpful if you wrote some text explaining what this is meant to do, so the OP doesn't have to take it on faith that it will just work.

@ninjagecko 2012-03-07 13:29:14

edit: This does not answer the question the OP meant to ask, leaving it here in case people find it useful somehow.

We use the multiplication rule of probability, combined with infinitessimals. This results in 2 lines of code to achieve your desired result:

longitude: φ = uniform([0,2pi))
azimuth:   θ = -arcsin(1 - 2*uniform([0,1]))

(defined in the following coordinate system:)

enter image description here

Your language typically has a uniform random number primitive. For example in python you can use random.random() to return a number in the range [0,1). You can multiply this number by k to get a random number in the range [0,k). Thus in python, uniform([0,2pi)) would mean random.random()*2*math.pi.


Proof

Now we can't assign θ uniformly, otherwise we'd get clumping at the poles. We wish to assign probabilities proportional to the surface area of the spherical wedge (the θ in this diagram is actually φ):

enter image description here

An angular displacement dφ at the equator will result in a displacement of dφ*r. What will that displacement be at an arbitrary azimuth θ? Well, the radius from the z-axis is r*sin(θ), so the arclength of that "latitude" intersecting the wedge is dφ * r*sin(θ). Thus we calculate the cumulative distribution of the area to sample from it, by integrating the area of the slice from the south pole to the north pole.

enter image description here (where stuff=dφ*r)

We will now attempt to get the inverse of the CDF to sample from it: http://en.wikipedia.org/wiki/Inverse_transform_sampling

First we normalize by dividing our almost-CDF by its maximum value. This has the side-effect of cancelling out the dφ and r.

azimuthalCDF: cumProb = (sin(θ)+1)/2 from -pi/2 to pi/2

inverseCDF: θ = -sin^(-1)(1 - 2*cumProb)

Thus:

let x by a random float in range [0,1]
θ = -arcsin(1-2*x)

@andrew cooke 2012-03-07 15:22:49

isn't this equivalent to the option he discarded as being "100% randomized"? my understanding is that he wants them to be more evenly spaced than a uniform random distribution.

@ninjagecko 2012-03-07 23:05:55

@BlueRaja-DannyPflughoeft: Hmm, fair enough. I guess I didn't read the question as carefully as I should have. I leave this here anyway in case others find it useful. Thanks for pointing this out.

@High Performance Mark 2012-03-07 12:03:15

Take the two largest factors of your N, if N==20 then the two largest factors are {5,4}, or, more generally {a,b}. Calculate

dlat  = 180/(a+1)
dlong = 360/(b+1})

Put your first point at {90-dlat/2,(dlong/2)-180}, your second at {90-dlat/2,(3*dlong/2)-180}, your 3rd at {90-dlat/2,(5*dlong/2)-180}, until you've tripped round the world once, by which time you've got to about {75,150} when you go next to {90-3*dlat/2,(dlong/2)-180}.

Obviously I'm working this in degrees on the surface of the spherical earth, with the usual conventions for translating +/- to N/S or E/W. And obviously this gives you a completely non-random distribution, but it is uniform and the points are not bunched together.

To add some degree of randomness, you could generate 2 normally-distributed (with mean 0 and std dev of {dlat/3, dlong/3} as appropriate) and add them to your uniformly distributed points.

@andrew cooke 2012-03-07 12:06:10

that would look a lot better if you worked in sin(lat) rather than lat. as it is, you will get a lot of bunching near the poles.

@robert king 2012-03-07 12:20:11

with small numbers of points you could run a simulation:

from random import random,randint
r = 10
n = 20
best_closest_d = 0
best_points = []
points = [(r,0,0) for i in range(n)]
for simulation in range(10000):
    x = random()*r
    y = random()*r
    z = r-(x**2+y**2)**0.5
    if randint(0,1):
        x = -x
    if randint(0,1):
        y = -y
    if randint(0,1):
        z = -z
    closest_dist = (2*r)**2
    closest_index = None
    for i in range(n):
        for j in range(n):
            if i==j:
                continue
            p1,p2 = points[i],points[j]
            x1,y1,z1 = p1
            x2,y2,z2 = p2
            d = (x1-x2)**2+(y1-y2)**2+(z1-z2)**2
            if d < closest_dist:
                closest_dist = d
                closest_index = i
    if simulation % 100 == 0:
        print simulation,closest_dist
    if closest_dist > best_closest_d:
        best_closest_d = closest_dist
        best_points = points[:]
    points[closest_index]=(x,y,z)


print best_points
>>> best_points
[(9.921692138442777, -9.930808529773849, 4.037839326088124),
 (5.141893371460546, 1.7274947332807744, -4.575674650522637),
 (-4.917695758662436, -1.090127967097737, -4.9629263893193745),
 (3.6164803265540666, 7.004158551438312, -2.1172868271109184),
 (-9.550655088997003, -9.580386054762917, 3.5277052594769422),
 (-0.062238110294250415, 6.803105171979587, 3.1966101417463655),
 (-9.600996012203195, 9.488067284474834, -3.498242301168819),
 (-8.601522086624803, 4.519484132245867, -0.2834204048792728),
 (-1.1198210500791472, -2.2916581379035694, 7.44937337008726),
 (7.981831370440529, 8.539378431788634, 1.6889099589074377),
 (0.513546008372332, -2.974333486904779, -6.981657873262494),
 (-4.13615438946178, -6.707488383678717, 2.1197605651446807),
 (2.2859494919024326, -8.14336582650039, 1.5418694699275672),
 (-7.241410895247996, 9.907335206038226, 2.271647103735541),
 (-9.433349952523232, -7.999106443463781, -2.3682575660694347),
 (3.704772125650199, 1.0526567864085812, 6.148581714099761),
 (-3.5710511242327048, 5.512552040316693, -3.4318468250897647),
 (-7.483466337225052, -1.506434920354559, 2.36641535124918),
 (7.73363824231576, -8.460241422163824, -1.4623228616326003),
 (10, 0, 0)]

@robert king 2012-03-07 12:26:36

to improve my answer you should change closest_index = i to closest_index = randchoice(i,j)

Related Questions

Sponsored Content

30 Answered Questions

[SOLVED] Is floating point math broken?

42 Answered Questions

27 Answered Questions

[SOLVED] Limiting floats to two decimal points

64 Answered Questions

[SOLVED] How do you split a list into evenly sized chunks?

5 Answered Questions

3 Answered Questions

[SOLVED] Distributing points evenly spaced within a sphere

  • 2019-01-24 01:52:43
  • Some1Else
  • 652 View
  • 1 Score
  • 3 Answer
  • Tags:   math

9 Answered Questions

3 Answered Questions

[SOLVED] Find evenly distributed random points on a spherical cap

  • 2019-01-19 10:38:30
  • Qwertie
  • 203 View
  • 1 Score
  • 3 Answer
  • Tags:   algorithm math geo

7 Answered Questions

9 Answered Questions

[SOLVED] Distribute points on a circle as evenly as possible

Sponsored Content