🔒 Locked 🎮 Unity

Fast way to generate random point on hemisphere?

Started by stuh505 Nov 15, 2006 at 9:02 AM 14 replies 14.5k views

Original Post

Nov 15, 2006 at 9:02 AM

I'm looking for a fast -- the fastest-- way to generate a uniformly distributed point on an a hemisphere centered about an arbitrary vector. I know a fast way to generate points on an axis-aligned hemisphere, but it relies on the cartesian axis and therefore isn't generalizable for an arbitrary hemisphere. I could use do a quaternion rotation (this is what I'm currently doing) but I'm thinking there may be a faster, more direct, way.

Zipster

Nov 15, 2006 at 9:27 AM

The method you use for an axis-aligned hemisphere can easily be generalized to work with any set of axis. Instead of the usual i = (1,0,0), j = (0,1,0), and k = (0,0,1), you use some other set a, b, and c, such that a x b = c, b x c = a, and c x a = b. Essentially, you work entirely in the rotated frame, instead of working in the world frame and then rotating.

stuh505

Nov 15, 2006 at 9:58 AM

Oh yeah, duh :)

But it can still be faster. Rather than generating an orthonormal basis (which would require doing the first step in Gram-Schmidt I think) we could just generate a pt on unit sphere and then take the dot product with the normal and negate the vector if its negative

Bob Janova

Nov 15, 2006 at 8:35 PM

Yeah, I'd say that's the fastest (create random points on a sphere, flip those in the wrong half then scale and translate). So the question becomes: do you have a fast way to generate points on a sphere?

deks

Nov 15, 2006 at 8:40 PM

Quote:
Original post by Bob Janova
Yeah, I'd say that's the fastest (create random points on a sphere, flip those in the wrong half then scale and translate). So the question becomes: do you have a fast way to generate points on a sphere?

Simply generate a random unit length vector, multiply it by the radius of your sphere then add the center of the sphere. It's then guaranteed to be on the sphere.

If this is not fast enough, you could create a 2d lookup-table of vectors uniformly distributed on the sphere, pickup a random float number then interpolate between the 4 vectors you fall in-between. I doubt this would be faster than normalizing then scaling the vector though, depends on the platform. What's interesting with this approach is that it may be tweaked so your emisphere is simply a range in the look-up table...

[Edited by - deks on November 15, 2006 9:40:02 PM]

alvaro

Nov 16, 2006 at 8:09 AM

This seems to be reasonably fast. If you have a fast InvSqrt, you can do faster.

#include <iostream>unsigned my_rand(void){  static unsigned next1=1151752134u, next2=2070363486u;  next1 = next1 * 1701532575u + 571550083u;  next2 = next2 * 3145804233u + 4178903934u;  return (next1<<16)^next2;}struct Vector3D {  float x, y, z;    Vector3D(float x, float y, float z):x(x),y(y),z(z){  }    float dot(Vector3D const &v){    return x*v.x + y*v.y + z*v.z;  }    Vector3D & operator+=(Vector3D const &v){    x += v.x;    y += v.y;    z += v.z;        return *this;  }};std::ostream & operator<<(std::ostream &o, Vector3D const &v){  return o << '(' << v.x << ',' << v.y << ',' << v.z << ')';}float U_m1_p1(){  return float(my_rand())*(1.0f/2147483648.0f) - 1.0f;}Vector3D pick_random_point_in_sphere(){  float x0,x1,x2,x3,d2;  do{    x0=U_m1_p1();    x1=U_m1_p1();    x2=U_m1_p1();    x3=U_m1_p1();    d2=x0*x0+x1*x1+x2*x2+x3*x3;  }while(d2>1.0f);  float scale = 1.0f/d2;  return Vector3D(2*(x1*x3+x0*x2)*scale,                  2*(x2*x3+x0*x1)*scale,                  (x0*x0+x3*x3-x1*x1-x2*x2)*scale);}Vector3D pick_random_point_in_semisphere(Vector3D const &v){  Vector3D result=pick_random_point_in_sphere();  if(result.dot(v)<0){    result.x=-result.x;    result.y=-result.y;    result.z=-result.z;  }  return result;}int main(){  Vector3D dir(1,0,0), sum(0,0,0);    for(int i=0;i<10000000;++i){    sum += pick_random_point_in_semisphere(dir);  }  std::cout << sum << std::endl;}

TheAdmiral

Nov 16, 2006 at 2:55 PM

Picking an arbitrary vector (uniform in 3-space) and normalising it will not give a uniform distribution over the sphere. The density at (1, 1, 1) will be more than that at (1, 0, 0) by a factor sqrt(3).

Edit: Removed lies [rolleyes].

Regards
Admiral

[Edited by - TheAdmiral on November 18, 2006 1:55:43 PM]

Airo

Nov 16, 2006 at 3:15 PM

I've used this, I think it is ok.

void RanUnitVector(void)
{
scalar a,b,l;

do
{
a = 1.0 - 2.0*ran2();
b = 1.0 - 2.0*ran2();
l = a*a + b*b;
}while(1
scalar s = sqrt(1.0 - l);

x = 2.0*a*s;
y = 2.0*b*s;
z = 1.0 - 2.0*l;
}

HemoGloben

Nov 16, 2006 at 3:43 PM

Quote:
Original post by TheAdmiral
Picking an arbitrary vector (uniform in 3-space) and normalising it will not give a uniform distribution over the sphere. The density at (1, 1, 1) will be more than that at (1, 0, 0) by a factor sqrt(3).

....

Regards
Admiral

I don't want to derail the discussion, but I think that's really interesting. I can almost see it, but can you explain it a bit more?

alvaro

Nov 16, 2006 at 3:53 PM

Quote:
Original post by TheAdmiral
Picking an arbitrary vector (uniform in 3-space) and normalising it will not give a uniform distribution over the sphere. The density at (1, 1, 1) will be more than that at (1, 0, 0) by a factor sqrt(3).

Of course. But you can fix this by discarding points outside of the sphere:

Vector3D pick_random_point_in_sphere(){  float x1,x2,x3,d2;  do{    x1=U_m1_p1();    x2=U_m1_p1();    x3=U_m1_p1();    d2=x1*x1+x2*x2+x3*x3;  }while(d2>1.0f);  float scale = 1.0f/sqrt(d2); // or use a fast InvSqrt for this  return Vector3D(x1*scale,x2*scale,x3*scale);}

You only need an average of 6/Pi trials (about 1.91).

alvaro

Nov 16, 2006 at 4:20 PM

There is quite a lot of info on how to generate random points on a sphere in this Wikipedia page. Airo seems to be using Marsaglia's method, which I didn't know about. I don't quite understand why it works yet...

mooserman352

Nov 16, 2006 at 10:36 PM

Quote:
Original post by HemoGloben
I don't want to derail the discussion, but I think that's really interesting. I can almost see it, but can you explain it a bit more?

all the vectors that lay on a particular ray through the origin get mapped on to the same point when normalized. so, for a point A of a sphere, the chance of a random vector being normalized to A is roughly proportional to the length of the ray that goes through A. the problem is, the distribution of raw vectors might be a cube, and obviously the length of the ray from the center of a cube to a corner is longer than from the center to a face.

TheAdmiral

Nov 17, 2006 at 11:58 AM

Quote:
Original post by HemoGloben
I don't want to derail the discussion, but I think that's really interesting. I can almost see it, but can you explain it a bit more?

Mooserman pretty much explained it, but I'll clarify:

The original proposition was: Generate three random numbers, x,y,z, in the range [-A,A]. These points will be uniformly distributed through the cube defined by [[-A,-A,-A], [A,A,A]]. Now normalise this vector so that it lies on the sphere.

If you pick any point on the sphere and consider the pre-image of this point under the proposed transformation (which is simply a normalisation), you'll find it looks like a ray segment from the origin (exclusive) through this point, up as far as the cube's surface (inclusive). Now the probability of ending up at a given point on the sphere is proportional to the length of this ray segment. Application of Pythagoras's theorem shows us that the 'diagonal' rays - the ones that strike a corner of the cube - have length sqrt(3)A, whereas those along an axis have length A. This skews the distribution in favour of the areas furthest from an axis.

Alvaro's solution seems fine: clip any original vector that lies a distance greater than A from the origin. However, I still don't like it. While it's sure to work in practice, theory says that you computer could be furiously churning out random vector after random vector, never escaping the while-loop, 'cause it can't generate a small enough vector. Sure enough, this is infinitely improbable, but I guess I have my pure-maths upbringing to blame [rolleyes].

In terms of speed, I'm sure your method, alvaro, would outperform one that uses trigonometry [smile].

Regards
Admiral

Bob Janova

Nov 17, 2006 at 8:13 PM

Quote:
Picking two random half-angles, translating them and throwing the results into the spherical polar coordinate transformation probably won't do much damage to your frame rate (though you never said what you're doing this for) and you're guaranteed a uniform distribution.

No you're not, spherical polars have a much higher density near the poles.

Extrarius

Nov 18, 2006 at 12:14 AM

Since I don't see this method mentioned, I figure I might as well mention it:

The way that comes to my mind is to simply generate points on a circle then calculates the Z coordinate based on the equation for a sphere:

Angle = Random() * M_PI * 2;Dist = Radius * sqrt(Random());X = Dist * cos(Angle);Y = Dist * sin(Angle);Z = sqrt(Radius * Radius - Dist * Dist);if(Random() < 0.5){   Z = -Z;}

where Random() returns a number in the range [0, 1)
Maybe not the fastest, but it doesn't have a loop that depends on random numbers so you can count on it never stalling your program =-)

Dirge

Nov 18, 2006 at 12:34 AM

Sphere point picking works incredibly nice for this. Heres a link: http://mathworld.wolfram.com/SpherePointPicking.html

Implement the described algorithm as code, adjust for a hemisphere and you're set. I implemented this for an ambient occlusion raytracer I made a while back and starfield generation and it looked great!

g-luck!

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.