# Spherical Harmonics Rotation

## Recommended Posts

DarkChris    110
How am I supposed to rotate 3 bands of spherical harmonics (9 coefficients) into tangent space?
Are there any good papers on that topic, or any code already written, since mine won't differ from the conventional way anyway.

##### Share on other sites
Steve132    433
You usually don't have to. If you just want to get diffuse lighting, evaluate the polynomial the SH coordinates represent at the normal vector.

##### Share on other sites
DarkChris    110
I know... But in my case I have to... I have to multiplicate them with other coefficients that are in tangent space (e.g. PRT).

##### Share on other sites
Steve132    433
The SH rotation gritty details paper has a good section on how to do this.

Also, look up "Stupid SH Tricks"

##### Share on other sites
DarkChris    110
Ok, but what kind of Rotation will I need to rotate them into Tangent Space. ZYZ, ZXZ or ZYX (Yaw, Pitch, Roll)?

I guess it's depending on how I calculate alpha, beta and gamma... And how do I calculate them?

##### Share on other sites
The best and most robust implementation of SH rotations I've found has been in the paper "Algorithms for Spherical Harmonic Lighting" by Ian Lisle and Tracy Huang.

-= Dave

##### Share on other sites
DarkChris    110
I solved it... I just rotated about Theta (Y) and Phi (Z). And now it's working perfectly fine.

Edit: Well. Or not.

Am I right that I don't need the Binormals nor the Tangents for the rotation into the tangent space, when I'm only using unshadowed radiance transfer?

And most importantly how can a scalar product over 2 pairs of 9 coefficients be the same as the integral? Let's just take a look at (1,-1). It's part of the linear band meaning that it describes if the axis it's describing is more positive or negative. Let's just say it's describing the x axis. If we take unshadowed radiance transfer, we need to parameterize the lighting and max(0, cos(theta)). Since the cosine is neither more right nor more left, the final coefficient will be 0. In fact, it will be non-zero on all (L,0) coefficients. It basically just says that it always blocks off all incident light the same independent of the x-axis. So if the actual incident light is stronger from the right, it will result in an coefficient being non-zero. Now if we multiply those together (since the integral is the same as the scalar product), we get zero since x * 0 = 0. But that's just not correct. There is more incident light from the right, than from the left.

[Edited by - DarkChris on November 16, 2010 12:02:54 PM]

##### Share on other sites
Quote:
 Original post by DarkChrisAnd most importantly how can a scalar product over 2 pairs of 9 coefficients be the same as the integral?

The integral disappears because of the orthonormal nature of the basis. As you may know from linear algebra the dot product of two vector's in an orthonormal basis are 1 when they are the same and 0 when they are not. In a functional basis the integral of the same basis is 1 when they are the same and 0 when they are different. A functional basis has infinite components and an integral is essentially a dot product. More generally, in a functional basis we call it the inner product.

Let's start off with the integral of two functions parameterized over the hemisphere (theta and phi).
$\int_{\Omega}A(\theta, \phi)B(\theta, \phi)d\omega$

Then we project both functions onto a spherical harmonic basis Y and reconstruct using its coefficients.
$A(\theta, \phi)=\sum_{l,m}A_{l,m}Y_{l,m}(\theta, \phi)$
$B(\theta, \phi)=\sum_{l,m}B_{l,m}Y_{l,m}(\theta, \phi)$
Now we have an alternate version of A and B so we can substitute them back into the original integral.
$\int_{\Omega}(\sum_{l,m}A_{l,m}Y_{l,m}(\theta, \phi))(\sum_{l,m}B_{l,m}Y_{l,m}(\theta, \phi))d\omega$
$=\int_{\Omega}\sum_{l,m}\sum_{l,m}A_{l,m}B_{l,m}Y_{l,m}(\theta, \phi)Y_{l,m}(\theta, \phi)d\omega$

Because of the linearality of integrals we can pull the summation outside of the integral and we can pull the constant coefficient multiplies out as well.
$=\sum_{l,m}\sum_{l,m}A_{l,m}B_{l,m}\int_{\Omega}Y_{l,m}(\theta, \phi)Y_{l,m}(\theta, \phi)d\omega$

Now the magic of orthogonality means that if Y_{l,m} != Y_{l,m} we get 0 when we integrate Y against itself. That means if l and m are not equal for both A and B we get 0. This simplifies it down to the integral.
$=\sum_{l,m}A_{l,m}B_{l,m}\int_{\Omega}Y_{l,m}(\theta, \phi)Y_{l,m}(\theta, \phi)d\omega$

Also, because $\int_{\Omega}Y_{l,m}(\theta, \phi)Y_{l,m}(\theta, \phi)d\omega=1$ because the basis are always the same we simplify the original equation down to
$=\sum_{l,m}A_{l,m}B_{l,m}*1$
$=\sum_{l,m}A_{l,m}B_{l,m}$

[Edited by - David Neubelt on November 16, 2010 4:07:27 PM]

##### Share on other sites
Quote:
 Original post by DarkChrisLet's just take a look at (1,-1). It's part of the linear band meaning that it describes if the axis it's describing is more positive or negative. Let's just say it's describing the x axis. If we take unshadowed radiance transfer, we need to parameterize the lighting and max(0, cos(theta)). Since the cosine is neither more right nor more left, the final coefficient will be 0.

No, that's not correct. Cosine has no dependence on its azimuthal angle which means it has the same coefficient no matter the angle of phi. The coefficient for 1,-1 will be same coefficient as it is for 1,0 (because phi doesn't matter!)

The first 9 coefficients of max(cosine(\theta), 0) are the following,

// constant band
A_{0, 0} = 3.141593f

// linear band
A_{1,-1} = 2.095395f
A_{1, 0} = 2.095395f
A_{1, 1} = 2.095395f

A_{2,-2} = 0.785398f
A_{2,-1} = 0.785398f
A_{2, 0} = 0.785398f
A_{2, 1} = 0.785398f
A_{2, 2} = 0.785398f

EDIT: The above is wrong and I will leave for reference.
* You do want to take those 9 coefficients and multiply them by your 9 SH coeffients.
* What I was wrong about was saying that the coefficients would be the same for a clamped cos but in fact if the coefficients had represented clamped cos(\theta) then they would be 0 for any band where m != 0 but these coefficients do not reflect a clamped cos(\theta) read my post further down for details on what the coefficients reflect.

[Edited by - David Neubelt on November 17, 2010 3:18:10 AM]

##### Share on other sites
Quote:
 Original post by DarkChrisAm I right that I don't need the Binormals nor the Tangents for the rotation into the tangent space, when I'm only using unshadowed radiance transfer?

I'm not quite sure how your rendering is setup but when I implemented PRT I never needed to rotate my lighting. The visibility and indirect lighting is implemented by simply scaling the incoming lighting.

##### Share on other sites
macnihilist    377
Quote:
 The first 9 coefficients of max(cosine(\theta), 0) are the following,[...]

That doesn't seem right. Shouldn't a clamped cosine centered at +z have at least two of the coefficients in the second band set to zero?

Also, I wanted to add that (f.f) == 1 is only true for an orthonormal basis (which is true for SH). An orthogonal basis only guarantees (f.g)==0 if f!=g. (Although you can easily get an orthonormal basis from an orthogonal one by normalizing the functions.)

##### Share on other sites
DarkChris    110
But my PRT is in Tangent Space while the Environment Coefficients are obviously in World Space.

I just figured out that it doesn't depend on phi. But somehow this gives me even worse results xD

If I use unshadowed radiance transfer... Just max[0, cos theta]... Why is it giving me absolutely wrong results when doing a scalar product between the coefficients you just send me and the lighting coefficients?

Edit: Ok, I know what's wrong. I finally need the binormals.

##### Share on other sites
Quote:
Original post by macnihilist
Quote:
 The first 9 coefficients of max(cosine(\theta), 0) are the following,[...]

That doesn't seem right. Shouldn't a clamped cosine centered at +z have at least two of the coefficients in the second band set to zero?

You are not going to represent a clamped cosine perfectly and remember the signal is a sum of the positive and negative bands. There is ringing that occurs on the bottom half.

Quote:
 Also, I wanted to add that (f.f) == 1 is only true for an orthonormal basis (which is true for SH). An orthogonal basis only guarantees (f.g)==0 if f!=g. (Although you can easily get an orthonormal basis from an orthogonal one by normalizing the functions.)

You are absolutely correct and I will go edit my post and update my terminology.

-= Dave

##### Share on other sites
Quote:
 Original post by DarkChrisBut my PRT is in Tangent Space while the Environment Coefficients are obviously in World Space.

You can compute your PRT in world space and just rotate your vectors as opposed to rotating your SH at runtime.

-= Dave

##### Share on other sites
DarkChris    110
In my case I can't... And which vectors? I thought we are talking about a simple scalar product?

I haven't fully rotated into tangent space, since I thought that irradiance doesn't care about the directions anyway. I didn't use any binormals or tangents at all, I just rotated back so that the normal would be pointing along Z. But is it possible that this might be causing the problems?

[Edited by - DarkChris on November 16, 2010 4:06:03 PM]

##### Share on other sites
macnihilist    377
Quote:
Original post by David Neubelt
Quote:
Original post by macnihilist
Quote:
 The first 9 coefficients of max(cosine(\theta), 0) are the following,[...]

That doesn't seem right. Shouldn't a clamped cosine centered at +z have at least two of the coefficients in the second band set to zero?

You are not going to represent a clamped cosine perfectly and remember the signal is a sum of the positive and negative bands. There is ringing that occurs on the bottom half.

True, but still a cosine in direction of +z is symmetric wrt x and y. So the projection of the cosine onto the basis functions that are just linear functions in x and y, respectively, should be zero.

Here is the coefficient for the SH that oscillates only in x (=cos(phi)). (Let's say its number is (1,-1) and N is the normalization factor):
$c_{1,-1} = \int_{\Omega}{Y_{1,-1}(\theta, \phi) \cos{\theta} \mathrm{d}\omega} = \int_{0}^{2\pi}{ \int_{0}^{\pi/2}{ N_{1,-1} \cos{\phi} \cos{\theta} \sin{\theta} \mathrm{d}\theta} \mathrm{d}\phi} = N_{1,-1} \int_{0}^{2\pi}{ \cos{\phi} \left( \int_{0}^{\pi/2}{ \cos{\theta} \sin{\theta} \mathrm{d}\theta} \right) \mathrm{d}\phi} = N_{1,-1} I_{\cos{\theta}\sin{\theta}} \int_{0}^{2\pi}{ \cos{\phi} \mathrm{d}\phi} = 0$

I'm pretty sure this is correct (apart from typos).

EDIT: Removed duplicate quote.

[Edited by - macnihilist on November 16, 2010 5:45:48 PM]

##### Share on other sites
Quote:
 Original post by macnihilistTrue, but still a cosine in direction of +z is symmetric wrt x and y. So the projection of the cosine onto the basis functions that are just linear functions in x and y, respectively, should be zero.

You are correct but let me clear up the confusion. I've misstated some facts and I've seem to make this problem more confusing.

First let me correct my mistake. The first 9 coefficients I wrote are not the projection of max(cos(\theta),0) onto a SH basis. You are correct a diffuse transfer function will have zero coefficients in any part of the basis where m != 0. Additionally, all odd functions l > 1 will be zero. The coefficients I wrote before actually represent the transfer function max(cos(\theta, 0) and the rotational operator to bring the lighting into the local space of the transfer function. You do want to multiply those 9 numbers I posted above by your incoming radiance.

The derivation I showed above is just a simple fact about inner product of orthonormal basis functions. When doing PRT you simply take the dot product of the coefficients vectors. When sampling irradiance you must scale by the 9 coefficients I posted above and then multiply by your basis functions. It doesn't seem to follow on why in one formula its a simple dot product and in another formula you have to evaluate the SH function.

Let me show you the original formula for irradiance to help explain where the numbers come from. Assuming a simple diffuse transfer function max(cos(\theta), 0) then we write irradiance as a parameterization of its surface normal,
$E(\alpha, \beta) = \int_{\Omega}L(\theta_i, \phi_i)A(\theta'_i)d\Omega$

The ' primes indicate local coordinates. The alpha and beta are the normal's orientation. L is radiance and A is the cos transfer function. As you see lighting get its radiance in a global coordinate frame but the transfer function is being evaluated in local coordinates. The transfer function cos is naturally evaluated in local coordinates so we will rotate our incoming lighting into a local coordinate frame.

We make a rotational operator R parameterized by our surface normal and rotate our global lighting into a local coordinate frame

$L(\theta_i, \phi_i) =L(R^{\alpha, \beta}(\theta'_i, \phi'_i))$

and substitute it into our original irradiance formula

$E(\alpha, \beta)= \int_{\Omega'}L(R^{\alpha,\beta}(\theta'_i,\phi'_i))A(\theta'_i)d\Omega'$

I'm going to skip a few steps a head and simply write all the coefficients those 3 functions end up creating. The steps I skipped are not important for understanding why you need to multiply each band by A_l and why you need to sample the SH basis Y. In the following equation you'll have you're rotational operator's coefficients D, lighting coefficients L and transfer coefficients A.
$E(\alpha, \beta)= \sum_{n=0}^{\infty}\sum_{l=0}^{\infty}\sum_{m=-l}^{\l}\sum_{m'=-l}^{l}L_{l,m}A_nD_{m,m'}^l(\alpha)\times exp(Im'\beta)exp(Im'\gamma)\int_{\phi'_i=0}^{2\pi}\int_{\theta'_i=0}^{2\pi}Y_{l,m'}(\theta'_i,\phi'_i)Y_{n,0}(\theta'_i, phi'_i)\times sin \theta'_i d\theta'_i d\phi'_i$

This equation because of orthogonality will cancel terms for anything that doesn't have n=l, m'=0. This gives us the equation
$E(\alpha, \beta)= \sum_{l=0}^{\infty}\sum_{m=-l}^{\l}L_{l,m}A_nD_{m,0}^l(\alpha) exp(Im'\beta)$

and finally we can make a substitution to get rid of D's coefficients and the exponential with the Legendre polynomials and a normalization term

$E(\alpha, \beta)= \sum_{l=0}^{\infty}\sum_{m=-l}^{\l}\left(\frac{4\pi}{2l + 1}\right)^\frac{1}{2}L_{l,m}A_nY_{l,m}(\alpha,\beta)$

A lot of what I posted is slight of hand mathematics where I skip a lot of steps but it should give you a high level view on why PRT is different then irradiance and why you need to multiply
A_0, A_1, A_2 by all your coefficients in each band and multiply by the sh basis in the direction of your normal.

It was also confusing for me at first on why PRT is a dot product but irradiance is formulated differently. At one point I dug into the derivation and it made sense but if you can remember the rules then there is no point going through the derivation yourself.

##### Share on other sites
Dark Chris, I'll post pseudo code on my PRT simulator later and hopefully it will help you out. I still don't understand why you need to rotate your SH.

##### Share on other sites
DarkChris    110
It seems like there are two possibilities to do PRT.

The tangent space one:
- Rotate L coefficients into tangent space
- E = L00 * A00 + L1m1 * A1m1 + L10 * A10 + L11 * A11 + ...
- Everything except Al,0 is always 0

The world space one:
- E(theta, phi) = L00 * A0 * Y00(theta, phi) + L1m1 * A1 * Y1m1(theta, phi) + L10 * A1 * Y10(theta, phi) + L11 * A1 * Y11(theta, phi) + ...
- A is only defined for l; all m's are the same

##### Share on other sites
Quote:
 Original post by DarkChrisIt seems like there are two possibilities to do PRT.The tangent space one: - Rotate L coefficients into tangent space - E = L00 * A00 + L1m1 * A1m1 + L10 * A10 + L11 * A11 + ... - Everything except Al,0 is always 0The world space one: - E(theta, phi) = L00 * A0 * Y00(theta, phi) + L1m1 * A1 * Y1m1(theta, phi) + L10 * A1 * Y10(theta, phi) + L11 * A1 * Y11(theta, phi) + ... - A is only defined for l; all m's are the same

Ok, I need to understand what you are doing.

Are you at minimum calculating a visibility term in SH for all your geometry in your scene?

If you are then the PRT SH visibility coefficients should already be in world space and a simple dot product against your irradiance env map projected into SH is all you need.

If you explain to me more what you are doing I can help you better.
-= Dave

##### Share on other sites
DarkChris    110
I don't want to "leak" too much at the moment, but trust me... I'm only able to calculate the transfer coefficients in tangent space. But I believe I can just replace the already existing A coefficients and than calculate the World Space PRT lighting, just like you tried to explain it to me.

##### Share on other sites
macnihilist    377
Quote:
 Original post by David Neubeltlet me clear up the confusion.

Yes, as I awoke this morning the suspicion entered my mind that the A's you wrote down are used to project an oriented clamped cosine into SH, while I was thinking about the actual SH coefficients of a fixed lobe aligned with the z axis. So the last equation you wrote down in the post above is basically a on-the-fly projection of the oriented cosine combined with the approximation of the double-product integral of the cosine factor and the incident radiance.
Thanks for the explanation and sorry for the interference.

EDIT: fixed typo

##### Share on other sites
Let's imagine a few different scenarios,

1) Combining normal mapping with world space PRT
If you are doing PRT and you want to do it in tangent space because you want to do normal mapping or evaluate a BRDF then you do the following.

Calculate the visibility, diffuse transfer, etc. in world space. Call that set of coefficients T for transfer.

Calculate the incident lighting (env map for example) and call that I.

Then you have your BRDF C, which is the A's, and your basis function y(n)

The final equation is,

y(n)^TCRMI.

This equation is described in more detail in the siggraph 2005 spherical harmonics course. The best description of the equation is in the "Normal Mapping for Precomputed Radiance Transfer" power point presentation by peter pike sloan. Please notice that C is a diagnol matrix with the cosine term which means the A coefficients get multiplied by all the incident lighting coefficients there are no zero's in the equation that cancel it out.

2) Computing the visibility term into a texture for runtime visibility normal mapping.

"Normal Mapping with Low-Frequency Precomputed Visibility"

In this case they store the visibility term in a texture and then store the incident lighting coefficients for every normal direction (in a cube map). Of course, even in this case you'd still want to computer your visibility in world space.

3) Store incident radiance in a sh lightmap and sample at runtime to combine with normal mapping (pretty much what halo 3 does and march of the froblins demo).

In this case you store your incident radiance in worldspace and use the formula from ravi's irradiance environment maps.

4) Basic PRT (such as the paper by peter pike sloan or robin green).
In this case you store PRT in world space and incident lighting in world space and do a dot product between the PRT and incident lighting SH.

Anyway, in almost all cases I've ever seen you would never multiply any part of the bands by 0. If you want further help you'll have to describe your situation, sorry.

-= Dave

##### Share on other sites
DarkChris    110
I guess, that's all the help I needed. I'm pretty confident that it'll work after the changes and will even be significantly faster since their won't be any rotation involved anymore.

I'll post an update if it actually worked tomorrow. But like always. I'm so glad that you helped me. Big thanks to both of you.