November 28th, 2016 Here I will only consider pseudo-random point generation. Well-spaced point generation is a different topic. The code examples are intended for clarity. Uniform points on the unit disc One method to generate a uniform random point $p$ on the unit disc ($\mathbb{D}$) is as follows. Generate two uniform values: $a \in \left[ 0,1 \right] $ , $ \theta \in \left[ -\pi, \pi \right) $ and apply1: \[p = \left(\sqrt{a}\cos\left(\theta\right), ~\sqrt{a}\sin\left(\theta\right)\right)\] Using symmetry and trig identities we can drop this down to one trig approximation and pair of square roots. Rejection method Choose a 2D point $p$ uniformly in each dimension: $p_x,p_y \in \left[ -1,1 \right] $ and reject values when outside the unit disc. The average number of iterations is $\frac{4}{\pi} \approx 1.27324$. Although previously known the first CS publication of this technique seems to be von Neumann 19512. An example implementation that returns $p \cdot p$ as well as the generated point: float uniform_disc(vec2_t* p) { float d,x,y; do { // lose one or both remaps (2*rng_f32()-1) // for uniform half or quarter disc x = 2.f*rng_f32()-1.f; // [-1,1) y = 2.f*rng_f32()-1.f; // [-1,1) d = x*x + y*y; // p.p } while(d >= 1.f); // branch prob: 1-pi/4 ~= 0.214602 p->x = x; p->y = y; return d; // p.p } A better method of producing range $[-1,1)$ is described here. By dropping one of the range remapping of $x$ or $y$ gives a uniform distribution on the half disc and dropping both the quarter disc. The efficiency of the these are identical to the full disc case. 125 250 500 1000 2000 Ad-hoc non-uniform If we rescale the component by some function of the dot product $d = (p \cdot p)$ we can get radial ad-hoc non-uniform distributions. Directly scaling by the dot-product: 125 250 500 1000 2000 Scaling by $\frac{d}{1+d}$: 125 250 500 1000 2000 Uniform points on unit circle or uniform 2D rotation Like the disc case we can generate a uniform point on the unit circle (1-sphere or $\mathbb{S}^1$) using trigonometry and symmetry/identities to reduce runtime complexity. \[\left(\cos\left(\theta\right), ~\sin\left(\theta\right)\right)\] In the same paper2 as the disc rejection method John von Neumann presented a pair of techinques to transform uniform points on the unit disc to uniform points on the circle. If we have a uniform point $p = \left(x,y \right)$ (considered as a complex number) on the unit disc, then the first method is to simply normalize: \[\left\Vert p \right\Vert = \frac{p}{\sqrt{p \cdot p}} = \frac{1}{\sqrt{x^2+y^2}} \left(x, y \right)\] The second method observes that squaring drops the square root and since it doubles the angle maintains uniformity. \[\left\Vert p \right\Vert^2 = \frac{p^2}{p \cdot p} = \frac{1}{x^2+y^2} \left(x^2-y^2, 2xy \right)\] The second method is widely known since it trades a square root for two products and two additions which historically was an assured win. With the introduction of fast inverse square root approximations the lesser known normalization method becomes interesting. Also choosing between zero, one or two NR iterations give a knob of speed vs. bias of the result. The remaining issue with a practical implementation is the degenerate case. I’ll assume that uniform float values are generated by “equidistance” method (pretty much the only one used). This method will generates values: $\frac{u}{2^{24}}$ where $u$ is a uniform integer on $\left[ 0, 2^{24} \right)$. The only problem result is then when both are zero. We can correct inverse square root by introducing a bias to the dot product which is sufficiently small that only contributes when both are zero such as $2^{-72}$. We also need to project that input to unit circle which requires choosing the branch point for the origin. Let’s send it to $(1,0)$ which needs adding a bias to the $x$ coordinate. This gives us: #define BIAS (1.f/68719476736.f) // 0x1p-36f void uniform_s1(vec2_t* p) { float d = uniform_disc(p); float s = rsqrt(d+(BIAS*BIAS)); // 1/sqrt approx + X NR steps p->x += BIAS; // map (0,0) -> (1,0) p->x *= s; p->y *= s; } Since we’re simply normalizing we can create half and quarter circles as per same disc above. 25 50 100 200 Ad-hoc non-uniform Taking the uniform variants and applying some simple transforms can yield some potentially interesting non-uniform distributions. If we take a uniform disc, translate it by one in x and then normalize we get a distribution that is roughly normal (standard deviation ~0.615113) on $\left[-\frac{\pi}{2},\frac{\pi}{2}\right] $ Using the angle-doubling version, well, doubles the angle output range. Uniform points on the 3D unit sphere In 1972 George Marsaglia3 introduced a method to transform a uniform point on the unit disc to a uniform point on the 3D sphere $\left(\mathbb{S}^2\right)$: \[\left(2x\sqrt{1-\left(x^2+y^2\right)},~2y\sqrt{1-\left(x^2+y^2\right)},~1-2\left(x^2+y^2\right)\right)\] I walked through to this equation in the post “Implied normals (unit bivectors as rotations)”4 where it’s described as the sequence: orthogonal project disc to half-sphere followed by a double-angle transform. Some bullet points are: an area-preserving map from disc to sphere (equivalent to Lamberts) radial lines (through origin) on disc are mapped to longitude circles on sphere (angle measure maintained) circles are mapped to latitudes (see below) $\left( 0, ~0 \right) $ is mapped to the north pole the disc: $r \in \left[ 0, ~\frac{1}{\sqrt{2}} \right] $ is mapped to the positive half sphere $\left( \mathbf{z} \geq 0 \right) $ the annulus: $r \in \left( \frac{1}{\sqrt{2}}, ~1 \right) $ is mapped to the negative half sphere $\left( \mathbf{z} < 0 \right )$ all points on the unit circle are mapped to south pole void uniform_s2(vec3_t* p) { float d,s; vec2_t v; d = uniform_disc(&v); s = 2.f*sqrtf(1.f-d); p->x = s*v.x; p->y = s*v.y; p->z = 1.f-2.f*d; } Given how points are mapped from the disc to the sphere, if we knew what “what circle maps to which latitude” we could form rotationally symettric distributions on the sphere directly from those on the disc. Dropping the angle-doubling map we end up with a cosine distribution on the half-sphere. Spherical caps As with non-uniform distributions, given the way the disc is mapped to the sphere we can directly formulate uniform points on spherical caps. The areas of the disc, sphere and spherical cap (height $h$) are respectively: \[A_d = \pi r^2 \\ A_s = 4\pi r^2 \\ A_c = 2\pi r h\] So the mapping from disc to sphere involves a constant factor of four and the area of the cap is linear with height. So we can form a point set on the unit disc, apply a uniform scaling and apply the disc to sphere mapping. The cap has a radius of one, divide out the constant factor of the disc to sphere map and relate the areas: \[\frac{2\pi h}{4} = \pi r^2\] Solving for $r$ give us: \[r = \sqrt{\frac{h}{2}}\] Plugging this scale factor into the map gives us the direct formulation: \[\left(x\sqrt{h\left(2-h\left(x^2+y^2\right)\right)},~y\sqrt{h\left(2-h\left(x^2+y^2\right)\right)},~1-h\left(x^2+y^2\right)\right)\] void uniform_s2_cap(vec3_t* p, float h) { float k,s; vec2_t v; k = h*uniform_disc(&v); // h(x^2+y^2) s = sqrtf(h*(2.f-k)); // sqrt( h(2-h(x^2+y^2)) ) p->x = s*v.x; p->y = s*v.y; p->z = 1.f-k; } Extra examples Cosine distribution on half-circle Generating a uniform line and projecting to circle results in $P(\theta) = \frac{1}{2}\cos(\theta)$ with $\theta \in \left[\frac{\pi}{2}, \frac{\pi}{2} \right]$ void cosine_s1(vec3_t* p) { float y = 2.f*rng_f32()-1.f; // (-1,1) float x = sqrtf(1.f-y*y); p->x = x; p->y = y; } Squaring this result converts the half-circle to circle with: $P(\theta) = \frac{1}{4}\cos(\frac{\theta}{2})$ Cosine distribution on half-sphere Seemingly the most common fragment shader version is as follows: vec3 cosDistHalfSphere(in vec3 n , in float seed) { float d2 = rng(seed++); float a = 2.0*PI*rng(seed++)); // generate a uniform point on disc (trig method) // p = ra(rx,ry) float ra = sqrt(d2); float rx = cos(a); float ry = sin(a); // project to sphere (Z) float rz = sqrt(1.0-d2); // an orthonormal basis (complement of 'n') vec3 xp = normalize(cross(n, vec3(0,1,1))); vec3 yp = cross(xp, n); // map 'p' to generated basis: ra*(rx*xp + ry*yp) // map 'z' to 'n': rz*n return normalize(ra*(rx*uu + ry*vv) + rz*n); } Note that the basis computation is degenerate as ‘n’ approaches $\pm \frac{1}{2}\left(0, \sqrt{2}, \sqrt{2}\right)$ and this requires the otherwise unneeded normalization of the final result to reduce the degenerate zone. So one option is use a method without a degenerate case5: vec3 cosDistHalfSphere(in vec3 n , in float seed) { float d2 = hash(seed++); float t = 2.0*PI*hash(seed++); // assuming hash result on [0,1) float ra = sqrt(d2); float rx = cos(t); float ry = sin(t); float rz = sqrt(1.0-d2); // some alternate non-degenerate method // see post including links float sz = n.z >= 0.0 ? 1.0 : -1.0; float a = n.y/(n.z+sz); float b = n.y*a; float c = n.x*a; vec3 xp = vec3(-n.z-b, c, n.x); vec3 yp = vec3(c, b-sz, n.y); return ra*(rx*xp + ry*yp) + rz*n; } Another option if the direction will never approach some fixed (say negative $z$) is the method described by Colin Barre-Brisebois and Stephen Hill6. References and Footnotes “Disk Point Picking”, Eric W. Weisstein, MathWorld ↩ “Various Techniques Used in Connection With Random Digits”, John von Neumann, 1951 (PDF) ↩ ↩2 “Choosing a Point from the Surface of a Sphere”, George Marsaglia, 1972 (PDF) ↩ “Implied normals (unit bivectors as rotations)”, (local post) ↩ “Orthonormal basis from normal via quaternion similarity”, (local post) ↩ “Blending in Detail”, Colin Barre-Brisebois and Stephen Hill, 2012 (link) ↩ Comments distribution (7) distribution (13) , random (16)