This is a quick extension to a previous post on uniform points on disc, circle, sphere and spherical caps1. We can also move away from trig or the rejection method by approximing equal-area disc2. I’ll assume basic understanding of representing rotations by quaternions. For a matrix simply follow through the deriviation by applying the similarity transform to an orthogonal basis.
Informally there are three general types of random rotations in 3D…one, two and three degrees of freedom.
About fixed axis one degree of freedom
The simplest general case is generating a random rotation about a fixed axis (considered as unit bivector
void uniform_quat_about_u(quat_t* q, vec3_t* u)
{
vec2_t p;
float d = uniform_circle_px(&p);
quat_set(q, p.y*u->x, p.y*u->y, p.y*u->z, p.x);
}
For non-uniform, if the distribution on the half-circle has angle probability function
With a fixed reference direction two degrees of freedom
In axis-angle speak moving to two degress of freedom fixes the axis-of-rotation to a plane. We can choose a fixed reference direction which is orthogonal to that plane (say positive
First generate a uniform point on the disc:
Find the quaternion that rotates
Composing the previous operations results in:
Breaking this down we have:
- the direction of the axis of rotation:
- CCW orthogonal to point on disc (which is logically in
-plane) - magnitude is logically
- CCW orthogonal to point on disc (which is logically in
- the scalar part is logically
We could have directly generated this result but requires figuring out how to set up a Jacobian…and I hate doing that. Since our values on the disc are random, we can also simply reorder the bivector part and drop the sign (not rotate) without any impact, giving us:
void uniform_quat_from_z(quat_t* q)
{
vec2_t p;
float d = uniform_disc(&p);
float s = sqrtf(1-d); // root on (0,1] can approx + fix-up
quat_set(q, p.x, p.y, 0.f s);
}
For non-uniform (with the reordering and sign drop) then if the probability on the disc is
The Real Deal (arbitrary) three degrees of freedom
George Marsaglia in same paper4 as uniform points on the 3D sphere
Given two uniform points on the unit disc:
Making an arbitrary choice for the scalar part (
void uniform_quat(quat_t* q)
{
vec2_t p0,p1;
float d1 = uniform_disc(&p1) + EPS;
float s1 = rsqrt(d1);
float d0 = uniform_disc(&p0); // or positive in 'x' since -Q & Q are equivalent
float s0 = sqrtf(1.f-d0);
float s = s0*s1;
quat_set(q, p0.y, s*p1.x, s*p1.y, p0.x);
}
References and Footnotes
-
“Uniform points on disc, circle, sphere and caps”, (local post) ↩ ↩2
-
“Square/Disc mappings”, (local post) ↩
-
“Implied normals (unit bivectors as rotations)”, (local post) ↩
-
“Choosing a Point from the Surface of a Sphere”, George Marsaglia, 1972 (PDF) ↩
-
Per Vognsen provides a simple justification in this gist: (link) ↩