Generate hemisphere around normal

(L) [2008/10/08] [jogshy] [Generate hemisphere around normal] Wayback!

Anybody can tell me how to implement a uniform and cosine-sampled hemisphere around a world-space normal with angle delimitation(a cone) in a robust way, pls? ( like Renderman's "gather" function ).
Currently i'm precomputing a hemisphere ( uniform and cosine distributed... through the specified angle interval ) aligned with Z up. Then, all the points are aligned to a specific normal using a Möller-Hughes quaternion-rotate-around-vector routine ( [LINK http://jgt.akpeters.com/papers/MollerHughes99/code.html] )... The problem is that the points(which are normalized vectors really) turn around suddenly when the normal is aligned with any world axis ( X=(1,0,0), Y=(0,1,0) or Z=(0,0,1) )... or near poles... which can cause strange artifacts.
I also tried to generate a sphere... then use dot products to discriminate the vectors those don't lie in the specified angle interval... but it's not good because I need a fixed amount of rays gathered around the normal... and the dot products gives me only an approximation to the desired # of rays.
Before you ask... yes, this is for ambient occlusion [SMILEY :lol:]
thanks.

(L) [2008/10/08] [greenhybrid] [Generate hemisphere around normal] Wayback!

How about sampling a disk with radius=1 uniformly, and then pushing that generated point-on-disk p(x,0,z) onto the sphere p(x,sqrt(1-x*x-z*z),z) (according to the Global Illumination Compendium this gives you a cosine distribution), and taking that point as a vector which to transform with a matrix you can create out of a single vector (pbrt has code for that, it looks like: Code: [LINK # Select all]void computeCoordinateSystem (Vector3d &v1, Vector3d &v2, Vector3d &v3) const {
v1 = computeNormal();
// from PBRT
if (fabs (v1.m[0]) > fabs (v1.m[1])) {
real ilen = 1.0 / sqrt (v1.m[0]*v1.m[0] + v1.m[2]*v1.m[2]);
v2 = Vector3d (-v1.m[2] * ilen, 0.0, v1.m[0] * ilen);
} else {
real ilen = 1.0 / sqrt (v1.m[1]*v1.m[1] + v1.m[2]*v1.m[2]);
v2 = Vector3d (0.0, v1.m[2] * ilen, -v1.m[1] * ilen);
}
v3 = v1.computeCross (v2);
}, where v1 would be your surface normal, note that the other axes are orthogonal but quasi random, say, not usable for anisotropic models)

(L) [2008/10/08] [davepermen] [Generate hemisphere around normal] Wayback!

depending on the distribution it doesn't work. but normally, just randomly create a vector length one (a point on the sphere surface). then make a dotproduct of that vector with the normal. multiply by the sign of the dotproduct to flip the vector if needed, to point in the correct direction.
that way, you can generate an exact amount of vectors, no dropouts.
now, to generate random points on a sphere, there are two ways. one is the russian roulete one: create a random vector with [-1,1] for each component, while length > 1, do it again. normalize the final vector.
the other one is algorithmically and is somewhere on flipcode afaik.
but this wouldn't work if you need a certain distribution function...

(L) [2008/10/08] [jogshy] [Generate hemisphere around normal] Wayback!

Thanks for links, they're very interesting. There's a thing I cannot understand tho:
Code: [LINK # Select all]void computeCoordinateSystem (Vector3d &v1, Vector3d &v2, Vector3d &v3) const
{
v1 = computeNormal();

// from PBRT
if (fabs (v1.m[0]) > fabs (v1.m[1]))
{
real ilen = 1.0 / sqrt (v1.m[0]*v1.m[0] + v1.m[2]*v1.m[2]);
v2 = Vector3d (-v1.m[2] * ilen, 0.0, v1.m[0] * ilen);
}
else
{
real ilen = 1.0 / sqrt (v1.m[1]*v1.m[1] + v1.m[2]*v1.m[2]);
v2 = Vector3d (0.0, v1.m[2] * ilen, -v1.m[1] * ilen);
}
v3 = v1.computeCross (v2);
}

That's not very robust... because "sqrt (v1.m[0]*v1.m[0] + v1.m[2]*v1.m[2])" can be near zero and then a divBy0 will occur.
[LINK http://jgt.akpeters.com/papers/MollerHughes99/code.html] seems to be more stable.
I'm using a similar mechamism but the points move suddenly near the world axis (XYZ) producing some artifacts due to numerical stability. At some point, the fabs (v1.m[0]) > fabs (v1.m[1]) can result in the branch above ( for example v1=(0.1,0.99) and the next v1=(0.1,0.101)... and the points will be suddenly rotated in a similar way than, without the pertinent detection, a vector camera that looks down turns suddenly around ). That will cause a heavy discontinuity for sure.
Perhaps it's better to gather the hemisphere directions directly around a world-space normal instead of pre-computing the hemisphere points and then rotate them using a normal. I wonder if there is any method to compute it in that way.
On the other hand, almost all these methods are based on random sampling(Monte Carlo Integration). Won't be better to sample uniformly + a small jitter the hemisphere so we can cover all the angle surface better? I mean a "for" from 0 to pi using small increments(+small jitter) or whatever.

(L) [2008/10/08] [beason] [Generate hemisphere around normal] Wayback!

>> jogshy wrote:Thanks for links, they're very interesting. There's a thing I cannot understand tho:
Code: [LINK # Select all]void computeCoordinateSystem (Vector3d &v1, Vector3d &v2, Vector3d &v3) const
{
v1 = computeNormal();

// from PBRT
if (fabs (v1.m[0]) > fabs (v1.m[1]))
{
real ilen = 1.0 / sqrt (v1.m[0]*v1.m[0] + v1.m[2]*v1.m[2]);
v2 = Vector3d (-v1.m[2] * ilen, 0.0, v1.m[0] * ilen);
}
else
{
real ilen = 1.0 / sqrt (v1.m[1]*v1.m[1] + v1.m[2]*v1.m[2]);
v2 = Vector3d (0.0, v1.m[2] * ilen, -v1.m[1] * ilen);
}
v3 = v1.computeCross (v2);
}

That's not very robust... because "sqrt (v1.m[0]*v1.m[0] + v1.m[2]*v1.m[2])" can be near zero and then a divBy0 will occur.

Are you sure? That divide is conditioned on v1.m[0]>v1.m[1]. Remember v1 is unit length so not all three components can be near zero. I think these conditions ensure you use the maximum magnitude component, so guaranteed no divBy0.

(L) [2008/10/08] [jogshy] [Generate hemisphere around normal] Wayback!

>> beason wrote:Are you sure? That divide is conditioned on v1.m[0]>v1.m[1]. Remember v1 is unit length so not all three components can be near zero. I think these conditions ensure you use the maximum magnitude component, so guaranteed no divBy0.
Well, if v1 has a length of one it's ok, but in case of very small triangles with a very small normal modulus perhaps it's a problem... after the computeNormal() an assert or any other comparison to test the normal's length will be good. To normalize very small vectors can be problematic, specially if you're using single precision.
What preocupies me is the sudden turn of the "if ( fabs > fabs )" brach.. .because it will suddenly rotate all the points. If you use the same point distribution for each pixel and you just align it using the surface normal then a discontinuity could appear.

(L) [2008/10/08] [beason] [Generate hemisphere around normal] Wayback!

If v1 was a vertex coordinate, then I would agree with you. v1 could be [0,0,0] and you would divide by zero.
However:
Code: [LINK # Select all]v1 = computeNormal();
...and I see you already know that. And you worry about the accuracy of the result of this computation. Well, doesn't your Moller technique that you prefer require a "to" vector, which you compute using the same technique? If so then it suffers from the same problems.
EDIT: You don't like like the Moller technique either. Sorry... okay I have to think about this some more [SMILEY :)]

(L) [2008/10/09] [jogshy] [Generate hemisphere around normal] Wayback!

>> beason wrote:EDIT: You don't like like the Moller technique either. Sorry... okay I have to think about this some more
Yep yep, I don't like any alignment solution because it causes discontinuities. Neither Monte Carlo sampling because generates too much noise for a low number of samples.
I think the solution is just to use a shared evenly-distributed sphere, then perform a dot product with the desired pixel normal and voila. The uniform sphere must contain 360*nSamplesPerPixel/coneAngle samples so you'll get a fixed # of AO dirs per pixel.
For the cosine-weighted version just maintain the evenly distribution but weight each resulting AO direction using the the dot of the desired normal and the AO dir.
And... if you like noise then just rotate randomly around a world axis the evenly distributed sphere for each pixel.

(L) [2008/10/09] [greenhybrid] [Generate hemisphere around normal] Wayback!

Another method lycium made me aware of once goes like:
Code: [LINK # Select all]v1 = computeNormal(); // if length is <> 1
tmp = randomUnitVector();
// assert that v1 and tmp are not parallel {fabs(v1*tmp-1)<epsilon} or
// accept that there is only a small probability that v1*tmp=1[+-]epsilon
v2 = v1.computeCross (tmp);
v3 = v1.computeCross (v2);

This produces a kind of noisy continuity [SMILEY :)]
I am currently writing a sky-map tool, where both methods aren't good to transform heightmap-normals to the hemisphere, so I am intersecting three rays with the sphere:
Code: [LINK # Select all]{pseudo code}
X = intersectSphere (A)
U = (intersectSphere (B) - X).computeNormal()
V = (intersectSphere (C) - X).computeNormal()

Where the rays A, B, C are generated out of a method f(p,q) -> (ray(x,y,z),(u,v,w)), like you do for your primary rays. Then A = f (x,y), B = f (x+1,y), C = f (x,y+1).
This doesn't give you a transformation where the spanning-vectors are orthogonal, of course, but it suits my needs to overcome the discontinuitites you mention.
I don't know how much this helps you, though.

(L) [2008/10/09] [raymo] [Generate hemisphere around normal] Wayback!

Another solution is to change the branch so it brances very rarely in the scene and should make the discontinuity a lot less noticeable.
i.e Change
Code: [LINK # Select all]If fabs(v1.m[0]) < fabs(v1.m[1] )
to
Code: [LINK # Select all]if (fabs(v1.m[0]) < 0.999f )
The division by zero case stills doesn't occurr and the discontinuity is only on very coplanar polygons ( which you probably have none with smooth normals ).

Another option is reversing the process. In my own tests I have found that image based lighting using a bunch of directional lights is roughly equivilant in speed to ambient occlusion with long rays ( where the rays are roughly half the scene size or greater), but the visual result is a lot better. Since the directional lights are all in world space you don't have to transform anything.

(L) [2008/10/09] [jogshy] [Generate hemisphere around normal] Wayback!

>> raymo wrote: In my own tests I have found that image based lighting using a bunch of directional lights is roughly equivilant in speed to ambient occlusion with long rays ( where the rays are roughly half the scene size or greater), but the visual result is a lot better. Since the directional lights are all in world space you don't have to transform anything.
Hmmm, that sounds interesting. Do you have any screenshot showing the difference?

Generate hemisphere around normal back

Board: Board index ‹ Ray tracing ‹ General Development