Uniform projections of spherical distributions

Imagine a light source above a plane of photosensitive paper that emits light into any given direction, such that every time it emits a ray of light, a single point on the paper lights up. In two dimensions, if each direction has equal probability of being chosen, we might picture something like this:


Here, the angle between any two rays is the same. Notice that the majority of the points lit up on the paper are directly beneath the light source, and few further away.

Now imagine we can control the angle θ\theta at which the source emits light. What if we wanted the distribution of points on the paper to be uniformly spaced? How should we pick θ\theta? What probability distribution of P(θ)P(\theta) do we sample θ\theta from? What if we considered the full two dimensional plane, how would we pick θ\theta and ϕ\phi together?

The aim of this blog post is to calculate angular distributions that create uniform projections on a plane, and learn some useful tricks along the way.

Stereographic projections

We'll start by picturing a one dimensional plane, i.e. a line, representing our photosensitive paper as in the above figure. For a ray at a given angle θ\theta, the point of intersection on the plane is calculated by considering the equation of a straight line – we position our light source at x=0x=0 and y=1y = 1 without loss of generality, with the gradient of each ray being tan(θ)-\tan (\theta), and find a given ray intercepts y=0y=0 at

0=tan(θ) x+1 x=1tan(θ)=cot(θ). 0 = -\tan(\theta)\, x + 1 \quad \implies \quad x = \frac{1}{\tan(\theta)} = \cot (\theta).

We'll express this for convenience as the action of the projection P:RR\mathcal{P} : \mathbb{R} \rightarrow \mathbb{R}, projecting an angle θ\theta to a point xx, which is a stereographic projection of the 11-Sphere onto the reals.

The inversion of this map, P1\mathcal{P}^{-1}, takes some point tt on R\mathbb{R} and gives us the angle θ\theta that a ray would need to intersect y=0y=0 at tt. Though P1\mathcal{P}^{-1} is already obvious, we'll make a quick detour to discuss, what Michael Spivak calls, the sneakiest substitution in the world:

Weierstrass substitution:

Let t=tan(θ2)t = \tan \left( \frac{\theta}{2} \right), then

sin(θ2)=t1+t2,andcos(θ2)=11+t2, \sin \left( \frac{\theta}{2} \right) = \frac{t}{1+t^2}, \quad \text{and} \quad \cos \left( \frac{\theta}{2} \right) = \frac{1}{1+t^2},

with θ(π,π)\theta \in \left( -\pi, \pi \right), and hence

sin(θ)=2t1+t2,andcos(θ)=1t21+t2. \sin \left( \theta \right) = \frac{2t}{1+t^2}, \quad \text{and} \quad \cos \left( \theta \right) = \frac{1-t^2}{1+t^2}.

These formulae are trivial to derive; the genius is the definition of tt. This substitution is used frequently in integral calculus to simplify trigonometric integrations, but it also sees use in stereographic projections. The Weierstrass substitution parametrizes the whole unit circle minus the point S=(1,0)S = (-1, 0), which is approached as t±t \rightarrow \pm \infty. Examining some values of tt and their corresponding point on the unit circle:


The points are initially spaced out angularly, and then asymptotically approach SS as tt becomes very large in magnitude.

To make use of this in our case, we want to move the unit circle to be centered on x=0x=0 and y=12y= \frac{1}{2}, and re-orient the parameterization such that SS is orientated on the yy axis (i.e. flip the xx and yy components). In this manner, t=0t=0 corresponds to a light ray traveling perpendicular to the paper, and increasing magnitudes of tt project points further afield. Our modified Weierstrass projection, which we'll denote W\mathcal{W}, is

t(x,y)=(t1+t2, 11+t2+1). t \mapsto (x, y) = \left( \frac{t}{1 + t^2}, \ \frac{-1}{1 + t^2} + 1 \right).

Note that the point (x,y)(x, y) will always lie on our circle; we have here a stereographic projection from R\mathbb{R} to a point on the 11-Sphere.

Similarly to eq. (1), we'll solve for the gradient of the ray intersecting (0,1)(0, 1) and a point given by eq. (4), namely

11+t2=m(t1+t2)+1, \frac{-1}{1 + t^2} = m \left( \frac{t}{1 + t^2} \right) + 1,

and therefore, we may identify

m=1t θ=arctan(1t). m = \frac{1}{t} \implies \theta = \arctan\left( \frac{1}{t} \right).

Which is precisely P1\mathcal{P}^{-1}. If tt is uniformly distributed, the set of angles we may calculate will ensure the projections are uniform as well (in this case trivially, since the tt are also the points of intersection). What we have practically done is calculate where the projected points lie, inverted the map, and found the angle which would intersect at the given point – and since we're working in the two dimensional case, we can use this as a parametrization.

Visualizing the Weierstrass projection

Implementing the key maps we have so far:

𝒲(t) = t / (1 + t^2), -1 / (1 + t^2) + 1
𝒫(θ) = cot(θ)
𝒫_inv(t) = atan(1/t)

For any point of intersection tt we can calculate the intersection with our circle with W\mathcal{W}, and the angle of that ray with P1\mathcal{P}^{-1}. We'll uniformly distribute tt:

# go from -60 degrees to 60 degrees
ts = range(-tan(π/3), tan(π/3), 12)
θs = 𝒫_inv.(ts)
points = 𝒲.(ts)

Plotting these points:


Shown here, as in the first figure, are the light rays emitted by the light source, the intersection with the plane (denoted by the thick ticks beneath the rays), and the intersections with the projecting 11-Sphere, shown as boxes. Now let's consider the 22-Sphere and points on a two dimensional plane.

Useful concepts and maps

Before we examine the three dimensional setup, there are a few things that are useful to cover:

We'll define our coordinate transformation from Cartesian to spherical polar as

(xyz)=(rsin(θ)cos(ϕ)rsin(θ)sin(ϕ)rcos(θ)), \left( \begin{matrix} x\\ y\\ z\\ \end{matrix} \right) = \left( \begin{matrix} r \sin(\theta) \cos (\phi)\\ r \sin(\theta) \sin (\phi)\\ r \cos(\theta)\\ \end{matrix} \right),

with θ\theta being the elevation angle, and ϕ\phi the azimuthal angle, following the convention that θ=0\theta = 0 corresponds to the positive zz axis, and ϕ=0\phi = 0 is the positive xx axis.

As an aside, a temptation for generating random points on the surface of the sphere is to sample θ\theta and ϕ\phi from uniform distributions between [0,π]\left[0, \pi \right] and [0,2π)\left[0, 2 \pi \right) respectively, however, such a sampling generates the majority of the points close to the poles. We already know, from studying the two dimensional case, that this will not provide us with a uniform projection, but it is still worth using as a comparison. The intuition for this is picturing circles of constant and equally spaced θ\theta on the surface of the sphere:


If we draw NN points on each ring, as a uniform distribution would, there would be a higher point density where the rings are smaller. We would require a measure that preserves the point density on the sphere – and fortunately the Jacobian tells us exactly what such a distribution would look like.

The Jacobian and distributions

For geometric purposes, the Jacobian matrix describes how each coordinate in the new coordinate system changes with respect to each coordinate in the old coordinate system. It is a matrix J\mathbf{J} where each matrix element is the partial derivative

Jij=yixj, \mathbf{J}_{ij} = \frac{\partial y_i}{\partial x_j},

where yiy_i are the new coordinates, and xjx_j are the old coordinates. The Jacobian tells us how lengths change under transformation at any point, and the determinant of the Jacobian tells us how area and volume elements transform,

dy1dy2dyn=det[J(x1,x2,,xm)]dx1dx2dxm. \text{d} y_1 \text{d} y_2 \cdots \text{d} y_n = \det \left[ \mathbf{J}(x_1, x_2, \ldots, x_m) \right] \text{d} x_1 \text{d} x_2 \cdots \text{d} x_m .

As point density is merely the number of points per unit volume, the determinant of the Jacobian also tells us the shape of distributions that preserve point density, as it tells us how we may preserve volumes. For the case of our spherical coordinates, where we have fixed the radius to r=1r=1, this is

det[J(r,θ,ϕ)]=xrxθxϕyryθyϕzrzθzϕr=1=sin(θ), \det \left[ \mathbf{J} \left( r, \theta, \phi \right) \right] = \left\lvert \begin{matrix} \frac{\partial x}{\partial r} & \frac{\partial x}{\partial \theta} & \frac{\partial x}{\partial \phi} \\ \frac{\partial y}{\partial r} & \frac{\partial y}{\partial \theta} & \frac{\partial y}{\partial \phi} \\ \frac{\partial z}{\partial r} & \frac{\partial z}{\partial \theta} & \frac{\partial z}{\partial \phi} \\ \end{matrix} \right\rvert_{r=1} = \sin (\theta),

and therefore, if we sample θ\theta from a sin\sin distribution, we will have uniform points on a sphere. But how do we generate non-uniform distributions such as sin\sin?

Parameterising non-uniform distributions

For some random uniformly sampled xinU(0,1)x_\text{in} \sim \mathcal{U}(0, 1), and a desired distribution P(θ)P(\theta) with θ[θ1,θ2)\theta \in \left[\theta_1, \theta_2\right), suppose we want to use xinx_\text{in} to generate a θout\theta_\text{out} which is distributed according to P(θ)P(\theta). In other words, we want the cumulative density functions (CDF) of the two to be equal

0xindx=θ1θoutP(θ)dθ, \int_{0}^{x_\text{in}} \text{d} x = \int_{\theta_1}^{\theta_\text{out}} P(\theta) \text{d} \theta,

such that when xin=1x_\text{in} = 1 and θout=θ2\theta_\text{out} = \theta_2, these integrals equate to unity (in practice, this often requires a normalization factor for P(θ)P(\theta)). The LHS is trivial, and allows us to make good use the inverse CDF method, or Smirnov transform:

Smirnov transform:

Let Q:RRQ : \mathbb{R} \rightarrow \mathbb{R} be invertible, which maps θoutQ(θout)=xin\theta_\text{out} \mapsto Q(\theta_\text{out}) = x_\text{in} with xinU(0,1)x_\text{in} \sim \mathcal{U}(0,1), and P(θ)P(\theta) be a normalized distribution. Defining the forward action

Q(θout):=θ1θoutP(θ)dθ, Q(\theta_\text{out}) := \int_{\theta_1}^{\theta_\text{out}} P(\theta) \text{d} \theta,

implies the reverse θout=Q1(xin)\theta_\text{out} = Q^{-1}(x_\text{in}) yields some θoutP(θ)\theta_\text{out} \sim P(\theta) from xinx_\text{in}.

The key observation that makes all of this useful is that for large NN, our random xinx_\text{in} is uniform, and therefore these probabilistic results may also be used for general uniform input.

We can use this to generate points under a sinθ\sin \theta distribution (c.f. eq. (10)); we restrict θ[0,π)\theta \in \left[0, \pi \right) and then define P(θ)=12sinθ P(\theta) = \frac{1}{2} \sin \theta such that it is normalized. Then

xin=120θoutsinθdθ=12(1cosθout), x_\text{in} = \frac{1}{2} \int_0^{\theta_\text{out}} \sin \theta \text{d} \theta = \frac{1}{2} \left( 1 - \cos \theta_\text{out} \right),

and therefore

Q1(x)=arccos(12x). Q^{-1}(x) = \arccos\left( 1 - 2x \right).

Note that to change the domain of θ\theta, we change the normalization factor AA in P(θ)P(\theta), which propagates to arccos(1xA)\arccos \left(1 - \frac{x}{A} \right). From some uniform xx, we now have a method of generating uniform points on a sphere. Expressing this as a couple of Julia functions, we use our Q1Q^{-1} and then transform our angles to Cartesian coordinates:

# returns θ, ϕ
angles_on_sphere(x1, x2) = (acos(1 - 2x1), 2π * x2)
# spherical to cartesian
spher_to_cart(r, θ, ϕ) = (r * sin(θ) * cos(ϕ), r * sin(θ) * sin(ϕ), r * cos(θ))

function point_on_sphere(x1, x2; R = 1.0)
    θ, ϕ = angles_on_sphere(x1, x2)
    spher_to_cart(R, θ, ϕ)

Note that x1 == x2 is ideal, however for illustrative purpose I have decoupled these inputs so we can draw more azimuthal points per ring. We use these functions to map rings onto the sphere

# six rings, only go to 0.5 for upper hemisphere
x1_range = range(1-cos(π/12), 0.5, 6)

# fine range to ensure good resolution 
# using the adjoint to broadcast with tensor products
x2_range = adjoint(0.0:0.01:1)

rings = point_on_sphere.(x1_range, x2_range)

Visualizing these rings on the upper hemisphere:


Left panel: the dashed lines are when θ\theta is equidistant, to contrast how the spacing changes; we see the difference between the rings is initial much greater, and converges with the equidistant lines the by the time we reach the equator. Right panel: the sin\sin distribution with θ\theta and ϕ\phi generating points on the upper hemisphere.

This is great and all, as for large samples the point density across the surface is constant, but there's still a problem – the method for generating ϕ\phi gives the same angles for each ring, resulting in slices down the great circle, like the pieces of an orange. We ideally desire some mechanism for generating a sprinkling of points on the surface of the 22-Sphere, which we can later use to test our projections.

Golden spirals

The golden spirals method generates discrete uniform points in polar coordinates. It works by placing a point, rotating by the golden ratio and moving outwards a little, then placing another point and repeating. This ensures no two points are ever directly intersected by the same line through the origin.

In other words, let tR+t \in \mathbb{R}^+ be positive, then the angles are generated by

ϕ=t π(1+5). \phi = t \, \pi \left( 1 + \sqrt{5} \right) .

We want the number of points in an annulus rr to r+drr + \text{d} r to be proportional to the area to engender constant point density, where the area element (derived from the Jacobian in two dimensions) is

dA=02πrdrdθ=2πdr. \text{d} A = \int_0^{2\pi} r \text{d} r \text{d} \theta = 2 \pi \text{d} r.

Using the Smirnov transform from eq. (12), we calculate a distribution for r[0,1]r \in \left[ 0, 1 \right], setting the normalized area distribution P(r)=2rP(r) = 2 r. Consequently we find

Q(r)=r2, Q(r) = r^2,

and therefore

Q1(t)=t. Q^{-1}(t) = \sqrt{t}.

The full action of the golden spiral map may now be be expressed:

Golden spiral map:

The golden spiral map creates uniform polar points within the unit circle, and is denoted G:R+R2\mathcal{G} : \mathbb{R}^+ \rightarrow \mathbb{R}^2. It maps

t(r,ϕ)=(ttmax, t π(1+5)), t \mapsto (r, \phi) = \left( \sqrt{\frac{t}{t_\text{max}}},\ t\, \pi\left( 1 + \sqrt{5} \right) \right),

for tR+t \in \mathbb{R}^+.

The division by tmaxt_\text{max} ensures that radial distribution is fixed between 00 and 11, i.e. the unit circle. In code:

golden_spiral(t, tmax) = [√(t/tmax), (t * π)*(1 + √5)]

# generate 300 points
ts = 0:300
points = map(t -> golden_spiral(t, last(ts)), ts)

This can be generalized to a sphere by calculated ϕ\phi via the golden spiral map, and θ\theta via the sin\sin distribution in eq. (14). Plotting these:


Left panel: plot of the 300 generated points in polar coordinates, with a bounding unit circle. If you look closely, you can see how the spiral paths seem to appear in the distribution. Right panel: using the golden spiral method to generate ϕ\phi for our uniform distribution of spherical coordinates.

This looks perfectly and equally even! However, its projections unfortunately won't be:


This suffers from the same pathologies as equidistant angles in the two dimensional case.

Uniform 22-Sphere projections

For the 22-sphere, we now consider P\mathcal{P} as projecting some θ\theta and ϕ\phi onto the xx-yy plane, which we may derive through simple geometry:

(θ,ϕ)(x,y)=(cot(θ2)cos(ϕ),cot(θ2)sin(ϕ)). (\theta, \phi) \mapsto (x, y) = \left( \cot\left(\frac{\theta}{2}\right) \cos(\phi), \cot\left(\frac{\theta}{2}\right) \sin(\phi) \right).

The normalizing factor 12\frac{1}{2} is to ensure we cover the full range of cot\cot. For the golden spiral projection eq. (19), we desire

(x,y)=(ttmaxcos(ϕ),ttmaxsin(ϕ)), (x, y) = \left( \sqrt{\frac{t}{t_\text{max}}} \cos (\phi), \sqrt{\frac{t}{t_\text{max}}} \sin (\phi) \right),

and therefore cot(θ2)=t\cot(\frac{\theta}{2}) = \sqrt{t}. Given our understanding so far, we can immediately infer, either using the Weierstrass projection, or simply by re-arranging eq. (20), that

t(θ,ϕ)=(2arctan(tmaxt ), t π(1+5)), t \mapsto (\theta, \phi) = \left( 2\arctan\left( \sqrt{\frac{t_\text{max}}{t}} \ \right), \ t\, \pi \left( 1 + \sqrt{5} \right) \right),

is a parameterization for generating uniform, equidistant projections, leveraging the golden spiral map and the Weierstrass substitution.

This map isn't quite P1\mathcal{P}^{-1}, since it does not map xx, yy to θ\theta, ϕ\phi, but it is precisely the map we want to be able to generate our projections. Well, almost – if we analyze the limits of tt, we find

t=tmax θ=2arctan(1)=π2, t = t_\text{max} \implies \theta = 2\arctan\left( 1 \right) = \frac{\pi}{2},

limiting our projection to equator of our sphere, corresponding to a 4545^\circ cone. However, this does have the practical benefit that we may scale the cone by simply dividing by a constant factor in the square root. Let's implement a possible version of this

angles_on_sphere(t; elevation_factor=1) = (
    2atan(√(t/elevation_factor)), (t * π)*(1 + √5)

This elevation_factor lets us now adjust the limits of our cone. As tt \rightarrow \infty, we should be able to theoretically cover the whole of the xx-yy plane in points. But for now, let's set the elevation_factor to t_max / 3and visualize a 6060^\circ cone:


Perfect! With this, let's reconsider our initial light source above photosensitive paper.


The parameterization given in eq. (22) allows us to take some uniformly distributed positive tR+t \in \mathbb{R}^+, and calculate angles which would project uniformly onto the xx-yy plane. The result is already quite interesting, but it's applications should hopefully prove fruitful:

The lamp post model is a coronal model for accreting black holes, which consists of a super-luminous cloud of pure electrons above the spin axis of the black hole. Simulations of such a model rely on using a geodesic integrator to approximate P\mathcal{P}, in this case mapping (θ,ϕ)(x0,x1,x2,x3)(\theta, \phi) \mapsto (x^0, x^1, x^2, x^3), a 4-vector describing the point of intersection with the accretion disc. In our software, we sample some NN points this way, and then use a custom Delaunay interpolation algorithm (there's a future blog post in how we devised this) to create a smooth field of 4-vectors over the disc.

Up until now, these models would evenly distribute angles across the entire sky of the corona, which resulted in the majority of the 4-vectors we end up with being close the the source, and extremely few towards the edges of the disc (c.f. eq. (14)). This situation is not ideal, and required extremely large NN to rectify – but now, with this method for uniform projections, we anticipate much better coverage.

This not only optimizes the run-time of the method, but should allow our interpolations to more faithfully reconstruct the true underlying field.