A note about the rejection method: in high dimensions, this becomes very expensive. This is because the volume of a ball inscribed in a unit cube goes to zero (very quickly) as the dimension goes to infinity, so you start rejecting points with high probability.
The graphics are nice, but this article really could have used some complexity analysis and some probability theory.<p>The author neither discusses the asymptomatic complexity of the algorithms, nor the run-time of his implementation (which is more pertinent here), nor gives any proofs of why some of the algorithms sample uniformly from the unit ball and why some of them don't.<p>Also, it would have been really nice to generalize this problem to n dimensions. I assume that there is some small value of n where naive rejection sampling is worse in practice than one of the more sophisticated methods.
"about 48% of the points chosen are discarded. This may not be a huge problem if the number of points that need to be generated is small"<p>... or if the number of points isn't small -- it's only 2x the cost of the simple linear case without a loss of uniformity or adding any costly trig functions. That's still pretty good!
Correct uniform sampling in things like this comes up a lot in Monte Carlo radiation transport problems and is covered very well in your basic first year graduate Monte Carlo class in a good nuclear engineering program. If you want to go deeper, check out the literature around that.
That (last method) is a lot of sines and cosines, cube root and an acos... I bet for 3D the rejection method is not that much slower, in wall clock time.<p>OK so I coded up the first one and the last one in C and compiled on my mac using whatever gcc links to. I found 1e7 points in the sphere:<p>first routine: 0.521s<p>last routine: 0.850s<p>[edit]: C code here: <a href="https://pastebin.com/zH0BehSM" rel="nofollow">https://pastebin.com/zH0BehSM</a>
This problem is equivalent to generating random orientations, which can be easily done generating random unit quaternions (see, for example, [1]). Each quaternion found will correspond to a point in the surface of a 4D sphere, which can be mapped to a point inside a 3D sphere.<p>[1] <a href="http://mathproofs.blogspot.com/2005/05/uniformly-distributed-random-unit.html" rel="nofollow">http://mathproofs.blogspot.com/2005/05/uniformly-distributed...</a>
Interestingly, it seems to be all about bending some distribution and then unbending it again. Orthogonal dimensions are independent, so choosing three random numbers generates uniform in-cube distribution. But once you bend it via polar coordinates or vectors, you “wrap” many points of polar/vector abstract space to near the cartesian origin. And then you have to bend back, to make it uniform again. That raises interesting question on whether it relates to topology problems. Since generating N randoms is a fastest primitive operation, is there a quick and uniform translation from a cube to the sphere, such that all points are still uniform? Can we apply space-filling curve magic here (like Hilbert’s)? I feel like there is much more geometry power under that, but sadly I’m just a guy with no math bg.
The main issue with the second method, which isn't discussed at all, is that it's more likely to pick a direction towards the corners of the unit cube. You can even see that bias in the point cloud.<p>The problem with choosing d uniformly is also better explained by pointing out that shells of larger radius have larger surface, thus they'll be underrepresented.
Interestingly, the method in which you determine your "random point" can greatly affect the end result, so swapping out the methods (e.g. "generate three coordinates and discard" method for the "generate a direction and a length" method) may produce wildly different results.<p>An interesting illustration of this is Bertrand's Paradox, which asks the question: "if you draw a random chord inside a circle, what is the probability that its length is greater than sqrt(3)?". The correct answer is 1/2. The correct answer is also 1/3. Oh, and the correct answer is also 1/4. It depends on three equally valid (and uniformly distributed) ways of defining your random chord.<p>See: <a href="http://web.mit.edu/tee/www/bertrand/" rel="nofollow">http://web.mit.edu/tee/www/bertrand/</a>
When the author mentions "normally distributed numbers" it sounds like he is sampling a normal (Gaussian) distribution. In fact, he is sampling from a uniform distribution.
For people who like the simplicity of the rejection algorithm, it might be worth checking out Ziggurat sampling which is based on a similar idea, but shapes the sampling region to reduce rejections: <a href="https://en.m.wikipedia.org/wiki/Ziggurat_algorithm" rel="nofollow">https://en.m.wikipedia.org/wiki/Ziggurat_algorithm</a>
Related: selecting points on the surface of a sphere <a href="https://www.jasondavies.com/maps/random-points/" rel="nofollow">https://www.jasondavies.com/maps/random-points/</a><p>Poisson disc sampling can be quite useful.
For the second example, sampling the radius uniformly leads to an error with a higher point density in the centre of sphere. The radius, r, has a weight proportional to r^2 or the surface area. So, one can take a uniform random sample, s between 0 and R^2, and take the square root of s. That will give higher probability density around larger radius. Then, the second picture should look like the first one.
Similar problem which ends up not being quite as simple in the general case as you might think: sampling from a unit simplex: <a href="http://www.cs.cmu.edu/~nasmith/papers/smith+tromble.tr04.pdf" rel="nofollow">http://www.cs.cmu.edu/~nasmith/papers/smith+tromble.tr04.pdf</a>
I think the section "Using normally distributed random numbers" deserves an additional explanation. The code may look simple, but like most of the other methods mentioned on the page, generating normally (Gaussian) distributed random numbers still requires the use of rejection sampling and sine/cosine functions. See <a href="https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform" rel="nofollow">https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform</a>
I thought this was kinda fun, and I was trying to figure out why you couldn't just take a random x and use it to find a y that preserves x²+y²=r² (ie by setting y to a random number between 0 and ±√(r²-x²)). In the end I figured out that you need to generate x the same way, but if you have, say: n=rand(r), x=rand(±√(r²-n²)), y=rand(±√(r²-x²)) it seems to work.<p>Here's a version I extended to three dimensions, designed to be run in the JS console of the original article (it'll show up underneath the first sphere):<p><pre><code> new SphereSimulation(document.getElementById("spheres1"), () => {
const rand = Math.random
const root = (x) => rand() < 0.5 ? Math.sqrt(x) : -Math.sqrt(x)
const invdist = (a=0, b=0) => root(RADIUS*RADIUS - a*a - b*b)
const x1 = rand() * invdist()
const x2 = rand() * invdist(x1)
const y2 = rand() * invdist(x2)
const x3 = rand() * invdist(x2, y2)
const y3 = rand() * invdist(x3, y2)
const z3 = rand() * invdist(x3, y3)
return new THREE.Vector3(x3, y3, z3)
})
</code></pre>
It looks okay to me - am I missing something? Is this a reasonable approach?
So what ends up being the general purpose method for uniformly distributed random numbers in some arbitrary shape? (other than the rejection method of course)<p>If arbitrary shapes always requires the rejection method, what are some "classes" that have easier algorithms (the sphere is obviously included here)
This is also very important in many global illumination techniques which requires sampling the hemisphere around a point.
As case in point a simple Monte Carlo path tracer uses this to sample random directions, see e.g. [0]. In addition to this there are many 'tricks' one can use to reduce the variance/noise using a differently weighted sampling (to help the converge of MC).
A nice technique to get a cosine-weighted hemisphere sampling is to generate a random point on the unit disc and then project them onto the (hemi)sphere.
Seems somewhat relevant here.<p>[0] <a href="https://en.wikipedia.org/wiki/Path_tracing#Algorithm" rel="nofollow">https://en.wikipedia.org/wiki/Path_tracing#Algorithm</a>
The method I immediately think of is just based on the relation.<p>Sqrt(x^2 + y^2 +z^2) <= 1
reduces to
x^2 + y^2 +z^2 <= 1<p>Generate the random value for x^2, any random value <= 1<p>Now we have y^2 +z^2 <= 1-x^2<p>Generate a random value for y^2, any random value <= 1-x^2<p>Finally generate a random value for z^2, any random value <= 1-x^2-y^2<p>Take the square roots of each number we generated to get the final coordinates.<p>But wait! Theres another step..the most important one! There are 6 orderings of (x,y,z) We must choose a random one because otherwise x will be biased towards higher values compared to y and z.<p>So I kind of lied. All we needed was 3 random values via the subtraction methods above. Then we need to randomly choose their ordering as (#,#,#)<p>This seems like the easiest most trivial solution that would be uniformly distributed.<p>Let me know if I am wrong.
Section 3 of this paper explains this topic pretty well in N-dimensions:<p><a href="http://jmlr.org/papers/volume17/blaser16a/blaser16a.pdf" rel="nofollow">http://jmlr.org/papers/volume17/blaser16a/blaser16a.pdf</a><p>To get a random point <i>within</i> the sphere (rather than on the surface of the sphere, as described in the paper) you also need to pick a random distance from the origin to randomly scale the point on the surface.
This is similar to the method used for generating random numbers using Box-Muller transform.<p><a href="https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform" rel="nofollow">https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform</a>
Maybe a mathematician can correct me. But could you not also model 3D vertices of uniform density into a 1D space and then uniformly sample the 1D space. Of course this may have issues with highly dense objects or very large ones.
Notice how the rejection method becomes useless for higher dimensions.<p><a href="https://en.m.wikipedia.org/wiki/Curse_of_dimensionality" rel="nofollow">https://en.m.wikipedia.org/wiki/Curse_of_dimensionality</a>
I don't agree with the way he normalizes the radius in the last 2 methods. I think it should be a square root, not a cubic one.
This is because the volume diferential r^2•sin(phi)•d_phi•d_theta•dr is proportional to the square of r.
This mistake is also made in a source he cites: <a href="https://math.stackexchange.com/questions/87230/picking-random-points-in-the-volume-of-sphere-with-uniform-probability/2872878" rel="nofollow">https://math.stackexchange.com/questions/87230/picking-rando...</a>
and I tried to point it out there as well.