I'm trying to implement Barnes-Hut algorithm, with a binary tree. My initial conditions are a uniform mass distribution in a sphere with radius $R$. How can I create uniform dots distribution in a sphere?

  • 2
    $\begingroup$ Do you want to distribute a totally arbitrary number of points $N$, or do you want to just have a grid of evenly distributed points (with some more implicit control parameter)? Do you distribute points only on the surface of the sphere or throughout the volume? $\endgroup$
    – Anton Menshov
    Aug 1, 2018 at 12:24
  • $\begingroup$ Just to make sure (question is a bit ambiguous), do you need to generate points on the sphere or in the ball ? $\endgroup$
    – BrunoLevy
    Aug 9, 2018 at 8:29
  • $\begingroup$ It sounds a reasonable question to me, I do not understand the downvotes (but I agree the question could be more precise...). $\endgroup$
    – BrunoLevy
    Aug 9, 2018 at 9:26
  • 1
    $\begingroup$ uniform mass distribution in a sphere ... sounds to me that the OP wants to generate points inside a ball of radius $R$ $\endgroup$
    – GoHokies
    Aug 9, 2018 at 16:51

2 Answers 2


To generate a specified number of random points in a 3D ball, one possibility is to use rejection sampling to have an initial point distribution, then improve the sampling with Lloyd relaxation.

Examples Examples with 30, 300 and 10000 points

Rejection sampling: Generate x,y,z coordinates using a uniform random number generator, ignore points that are outside the ball, and continue until you get the required number of points.

Lloyd relaxation It is an iterative method that makes the distribution of points more regular. Each iteration works as follow: compute the Voronoi diagram of the points, then move each point to the barycenter of its Voronoi cell. See [1] for more detailed explanations. There is an implementation in my Geogram software library [2]. With Geogram, you need to use a tetrahedralized version of the sphere.

In volume / on surface There is also a version of the algorithm that works on a 3D surface (sphere or arbitrary surface), called Constrained Centroidal Voronoi Tesselation [3], also implemented in Geogram [2].

You may also use my software Graphite [4] (graphic user interface around geogram). To sample a sphere in Graphite:

Scene->Create object->OGF::MeshGrob
Surface->Shapes->Create sphere
Points->Sample volume or Points->Sample surface


minitutorial: sampling a ball in Graphite

[1] https://en.wikipedia.org/wiki/Lloyd%27s_algorithm

[2] http://alice.loria.fr/software/geogram/doc/html/index.html

[3] Constrained Centroidal Voronoi Tesselation, Du, Gunzburger and Ju, SIAM J. on computing, 2003

[4] http://alice.loria.fr/software/graphite/doc/html/


There is a nice document:

which basically summarizes the two easiest methods to generate a distribution of points on a surface of a sphere:

  1. Uniform random distribution (generate $N$ random $z$ and $\varphi$ in $[-R, R]$ and $[0,2\pi]$, respectively and calculate $(x,y)$ as $(\sqrt{r^2-z^2}\cos{\varphi},\sqrt{r^2-z^2}\sin{\varphi})$. However, with small enough $N$ the outpput does not look too uniform.

  2. Generating circles of latitude with constant intervals with points at almost the same distance from each other (see the paper above). Does not generate EXACTLY the requested number of points, but tries to do it close enough.

A sample Python script I quickly wrote:

import numpy as np 

n = 5000
r = 1
generateRandom = True
generateFixed = True

if (generateRandom):
    print("Generating randomly %d points on a sphere centered at the origin" % (n))
    theta = np.random.uniform(0.0,2*np.pi,n)
    z = np.random.uniform(-1.0,1.0,n)
    for i in range (0,n):
        zp = z[i]
        xp = np.sqrt(r*r - zp*zp)* np.cos(theta[i])
        yp = np.sqrt(r*r - zp*zp)* np.sin(theta[i])
if (generateFixed):
    print("Generating fixed %d points on a sphere centered at the origin" % (n))
    alpha = 4.0*np.pi*r*r/n
    d = np.sqrt(alpha)
    m_nu = int(np.round(np.pi/d))
    d_nu = np.pi/m_nu
    d_phi = alpha/d_nu
    count = 0
    for m in range (0,m_nu):
        nu = np.pi*(m+0.5)/m_nu
        m_phi = int(np.round(2*np.pi*np.sin(nu)/d_phi))
        for n in range (0,m_phi):
            phi = 2*np.pi*n/m_phi
            xp = r*np.sin(nu)*np.cos(phi)
            yp = r*np.sin(nu)*np.sin(phi)
            zp = r*np.cos(nu)
            count = count +1

enter image description here

  • $\begingroup$ I couldn't understand your code. I wrote another version (here), let me know if my interpretation is correct. $\endgroup$
    – nicoguaro
    Aug 7, 2018 at 21:54
  • $\begingroup$ @nicoguaro I don't think the second version of the generation is correct. The logic is a bit more complicated to generate it uniformly: you cannot just have a fixed number of $\phi$ and $\theta$. Take a look at the reference I provided. $\endgroup$
    – Anton Menshov
    Aug 7, 2018 at 22:22
  • $\begingroup$ Well, I didn't check the reference. But in your code there is just a triplet (x,y,z) at the end. $\endgroup$
    – nicoguaro
    Aug 7, 2018 at 22:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.