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 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 ? What probability distribution of do we sample from? What if we considered the full two dimensional plane, how would we pick and 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.
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 , the point of intersection on the plane is calculated by considering the equation of a straight line – we position our light source at and without loss of generality, with the gradient of each ray being , and find a given ray intercepts at
We'll express this for convenience as the action of the projection , projecting an angle to a point , which is a stereographic projection of the -Sphere onto the reals.
The inversion of this map, , takes some point on and gives us the angle that a ray would need to intersect at . Though is already obvious, we'll make a quick detour to discuss, what Michael Spivak calls, the sneakiest substitution in the world:
Let , then
with , and hence
These formulae are trivial to derive; the genius is the definition of . 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 , which is approached as . Examining some values of and their corresponding point on the unit circle:
The points are initially spaced out angularly, and then asymptotically approach as becomes very large in magnitude.
To make use of this in our case, we want to move the unit circle to be centered on and , and re-orient the parameterization such that is orientated on the axis (i.e. flip the and components). In this manner, corresponds to a light ray traveling perpendicular to the paper, and increasing magnitudes of project points further afield. Our modified Weierstrass projection, which we'll denote , is
Note that the point will always lie on our circle; we have here a stereographic projection from to a point on the -Sphere.
Similarly to eq. (1), we'll solve for the gradient of the ray intersecting and a point given by eq. (4), namely
and therefore, we may identify
Which is precisely . If is uniformly distributed, the set of angles we may calculate will ensure the projections are uniform as well (in this case trivially, since the 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.
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 we can calculate the intersection with our circle with , and the angle of that ray with . We'll uniformly distribute :
# 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 -Sphere, shown as boxes. Now let's consider the -Sphere and points on a two dimensional plane.
Before we examine the three dimensional setup, there are a few things that are useful to cover:
the Jacobian in the context of distributions,
parameterizing non-uniform distributions, and
methods for generating uniform points on spheres.
We'll define our coordinate transformation from Cartesian to spherical polar as
with being the elevation angle, and the azimuthal angle, following the convention that corresponds to the positive axis, and is the positive axis.
As an aside, a temptation for generating random points on the surface of the sphere is to sample and from uniform distributions between and 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 on the surface of the sphere:
If we draw 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.
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 where each matrix element is the partial derivative
where are the new coordinates, and 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,
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 , this is
and therefore, if we sample from a distribution, we will have uniform points on a sphere. But how do we generate non-uniform distributions such as ?
For some random uniformly sampled , and a desired distribution with , suppose we want to use to generate a which is distributed according to . In other words, we want the cumulative density functions (CDF) of the two to be equal
such that when and , these integrals equate to unity (in practice, this often requires a normalization factor for ). The LHS is trivial, and allows us to make good use the inverse CDF method, or Smirnov transform:
Let be invertible, which maps with , and be a normalized distribution. Defining the forward action
implies the reverse yields some from .
The key observation that makes all of this useful is that for large , our random is uniform, and therefore these probabilistic results may also be used for general uniform input.
We can use this to generate points under a distribution (c.f. eq. (10)); we restrict and then define such that it is normalized. Then
Note that to change the domain of , we change the normalization factor in , which propagates to . From some uniform , we now have a method of generating uniform points on a sphere. Expressing this as a couple of Julia functions, we use our 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, θ, ϕ) end
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 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 distribution with and 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 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 -Sphere, which we can later use to test our projections.
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 be positive, then the angles are generated by
We want the number of points in an annulus to to be proportional to the area to engender constant point density, where the area element (derived from the Jacobian in two dimensions) is
Using the Smirnov transform from eq. (12), we calculate a distribution for , setting the normalized area distribution . Consequently we find
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 . It maps
The division by ensures that radial distribution is fixed between and , 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 via the golden spiral map, and via the 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 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.
For the -sphere, we now consider as projecting some and onto the - plane, which we may derive through simple geometry:
The normalizing factor is to ensure we cover the full range of . For the golden spiral projection eq. (19), we desire
and therefore . Given our understanding so far, we can immediately infer, either using the Weierstrass projection, or simply by re-arranging eq. (20), that
is a parameterization for generating uniform, equidistant projections, leveraging the golden spiral map and the Weierstrass substitution.
This map isn't quite , since it does not map , to , , but it is precisely the map we want to be able to generate our projections. Well, almost – if we analyze the limits of , we find
limiting our projection to equator of our sphere, corresponding to a 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) )
elevation_factor lets us now adjust the limits of our cone. As , we should be able to theoretically cover the whole of the - plane in points. But for now, let's set the
t_max / 3and visualize a 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 , and calculate angles which would project uniformly onto the - 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 , in this case mapping , a 4-vector describing the point of intersection with the accretion disc. In our software, we sample some 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 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.