Memo on Fresnel equations

Version : 1.3 – Living blog – First version was 29 April 2013

This post is a memo for me about misc thing I found related to the Fresnel equation as there is plenty of formula into the wild. I spend some times to gather these information and I was thinking it can interest others. This will not be really relevant to game rendering but still good to know. This could be useful when doing reference or when dealing with total internal reflection (frequent with multi layered BRDF). This is a memo, not a tutorial, so I won’t give basic explanation of many concepts like reflection/refraction, Snell’s law, index of refraction (IOR), total internal reflection (TIR), etc… At the end of the post I provide a Mathematica file with all the equations and graphs.

Notation for this post:

Conductor mean Metal and dielectric mean no-Metal material.

Light moves from a medium of a given IOR n_i (incoming) into a second medium with IOR n_t (transmitted).

Conductor have complex IOR with an imaginary part k_t (note that’s t “transmitted” is a bad choice for conductor but I found it more identifiable than n_1 and n_2, k_2).

Fresnel Equation basis

All equations using an index of refraction n_t can be replace with the same equation using a complex index of refraction n_t -\mathbf{i} k_t.

Snell’s law dielectric-conductor interface : \frac{\sin\theta_i}{\sin\theta_t}=\frac{n_t -\mathbf{i} k_t}{n_i}

Snell’s law dielectric-dielectric interface : \frac{\sin\theta_i}{\sin\theta_t}=\frac{n_t}{n_i}

The calculations of the reflectance R (What we are looking for when we want to calculate the percentage of reflection and transmission) depend on p- and s-polarization of the incident ray. Rs and Rp are the reflectivity for the two planes of polarization. Rs is perpendicular (s = German senkrecht) and Rp is parallel. The reflectance R for unpolarized light is the average of Rs and Rp:

R=\frac{(R_s+R_p)}{2}=\frac{(r_\perp^2+r_\parallel^2)}{2}

Following expressions use R_s and R_t or r_\perp and r_\parallel depend on cases to simplify notation.

Fresnel equation for dielectric-conductor interface is [5]:

r_\perp=\frac{n_i \cos\theta_i - (n_t -\mathbf{i} k_t)\cos\theta_t }{n_i \cos\theta_i + (n_t -\mathbf{i} k_t)\cos\theta_t }

r_\parallel=\frac{(n_t -\mathbf{i} k_t) \cos\theta_t - n_i\cos\theta_i }{(n_t -\mathbf{i} k_t)\cos\theta_t + n_i\cos\theta_i }

Fresnel equation for dielectric-dielectric interface is [2]:

r_\perp=\frac{n_i \cos\theta_i - n_t\cos\theta_t }{n_i \cos\theta_i + n_t\cos\theta_t }

r_\parallel=\frac{n_t \cos\theta_t - n_i\cos\theta_i }{n_t \cos\theta_t + n_i\cos\theta_i }

Note as the expression is square “i” and “t” can be swap.

Dielectric Fresnel’s formulas are always derive from conductor Fresnel’s formulas by settings k_t=0.

Several Fresnel’s formulas assume an Air-Dielectric/Conductor interface. Air IOR is 1.000293 and is approximate to 1. Thus using n_i=1 allow many simplification.

We define \eta=\frac{n_t}{n_i} and \eta_k=\frac{k_t}{n_i} the ratio of indices of refraction at the surface interface.
Fresnel’s formulas with n_i, n_t and k_t can always be replace by Fresnel’s formulas with \eta and \eta_k.

All Fresnel’s formulas of the rest of the post are based on derivation of the Fresnel equation above. Arrangement is done to express and optimize Fresnel’s formulas depending on \theta_i alias \theta in this post using expression derives from Snell’s law:

\cos\theta_t=\sqrt{1 - (\frac{n_i}{n_t}\sin\theta_i)^2}

\cos\theta_t= \sqrt{1 - (\frac{n_i}{n_t}^2(1-\cos^2\theta_i))}

Fresnel Equation in practice

Dielectric-Conductor interface

From [1]:

a^2=\frac{1}{2n_i^2}(\sqrt{(n_t^2-k_t^2-n_i^2\sin^2\theta)^2+4n_t^2 k_t^2}+n_t^2-k_t^2-n_i^2\sin^2\theta)

b^2=\frac{1}{2n_i^2}(\sqrt{(n_t^2-k_t^2-n_i^2\sin^2\theta)^2+4n_t^2 k_t^2}-n_t^2+k_t^2+n_i^2\sin^2\theta)

R_s=\frac{a^2+b^2-2 a\cos\theta + \cos^2\theta }{a^2 + b^2 + 2 a\cos\theta + \cos^2\theta}

R_p=R_s\frac{a^2+b^2-2 a\sin\theta \tan\theta + \sin^2\theta \tan^2\theta}{a^2+b^2+2 a\sin\theta \tan\theta + \sin^2\theta \tan^2\theta}

using a^2 and b^2 with \eta and \eta_k as follow give the same result:

a^2=\frac{1}{2}(\sqrt{(\eta^2-\eta_k^2-\sin^2\theta)^2+4\eta^2 \eta_k^2}+\eta^2-\eta_k^2-\sin^2\theta)

b^2=\frac{1}{2}(\sqrt{(\eta^2-\eta_k^2-\sin^2\theta)^2+4\eta^2 \eta_k^2}-\eta^2+\eta_k^2+\sin^2\theta)

Derivation from conductor Fresnel equation can be found in [5] (p.111) and use

(n_t -\mathbf{i} k_t)\cos\theta_t= (n_t -\mathbf{i} k_t)\sqrt{1 - (\frac{n_i}{(n_t -\mathbf{i} k_t)}\sin\theta)^2}=\sqrt{(n_t -\mathbf{i} k_t)^2 - n_i^2\sin^2\theta}

a -\mathbf{i}b=\sqrt{(n_t -\mathbf{i} k_t)^2 - n_i^2\sin^2\theta)}

In practice there is some simplification possible:

a^2 + b^2 =\sqrt{(\eta^2-\eta_k^2-\sin^2\theta)^2+4\eta^2 \eta_k^2}

\tan\theta = \frac{\sin \theta}{\cos \theta}

R_p=R_s\frac{ \cos^2 \theta (a^2+b^2)-2 a\cos\theta \sin^2\theta + \sin^4\theta}{\cos^2 \theta (a^2+b^2)+2 a\cos\theta \sin^2\theta + \sin^4\theta}

resulting in hlsl code:

float3 FresnelDieletricConductor(float3 Eta, float3 Etak, float CosTheta)
{  
   float CosTheta2 = CosTheta * CosTheta;
   float SinTheta2 = 1 - CosTheta2;
   float3 Eta2 = Eta * Eta;
   float3 Etak2 = Etak * Etak;

   float3 t0 = Eta2 - Etak2 - SinTheta2;
   float3 a2plusb2 = sqrt(t0 * t0 + 4 * Eta2 * Etak2);
   float3 t1 = a2plusb2 + CosTheta2;
   float3 a = sqrt(0.5f * (a2plusb2 + t0));
   float3 t2 = 2 * a * CosTheta;
   float3 Rs = (t1 - t2) / (t1 + t2);

   float3 t3 = CosTheta2 * a2plusb2 + SinTheta2 * SinTheta2;
   float3 t4 = t2 * SinTheta2;   
   float3 Rp = Rs * (t3 - t4) / (t3 + t4);

   return 0.5 * (Rp + Rs);
}

A more efficient approximation for Air-Conductor interface is available in [3]:

Rs=\frac{(n_t^2+k_t^2)- 2n_t \cos \theta + \cos^2 \theta }{(n_t^2+k_t^2)+ 2n_t \cos \theta + \cos^2 \theta}

Rp=\frac{(n_t^2+k_t^2)\cos^2 \theta - 2n_t \cos \theta + 1 }{(n_t^2+k_t^2)\cos^2 \theta + 2n_t \cos \theta + 1}

Calling the above expressions with \eta and \eta_k instead of n_t and k_t allow to approximate Dielectric-Conductor interface. Hsls code is:

float3 FresnelDieletricConductorApprox(float3 Eta, float3 Etak, float CosTheta)
{
    float  CosTheta2 = CosTheta * CosTheta;
    float3 TwoEtaCosTheta = 2 * Eta * CosTheta;

    float3 t0 = Eta * Eta + Etak * Etak;
    float3 t1 = t0 * CosTheta2;
    float3 Rs = (t0 - TwoEtaCosTheta + CosTheta2) / (t0 + TwoEtaCosTheta + CosTheta2);
    float3 Rp = (t1 - TwoEtaCosTheta + 1) / (t1 + TwoEtaCosTheta + 1);

    return 0.5* (Rp + Rs);
}

For k_t=0 the approximation (Blue curve) is not good, for other case, the approximation is great:

FresnelConductorApproximation

This formulation is really specialized for conductor and should not be used for Dielectric-Dielectric interface.

Dielectric-Dielectric interface

From [2]:

r_\perp=\frac{n_i \cos\theta -n_t\sqrt{1 - (\frac{n_i}{n_t}^2 sin^2\theta)}}{n_i \cos\theta +n_t\sqrt{1 - (\frac{n_i}{n_t}^2 sin^2\theta)}}=\frac{n_i \cos\theta -n_t\sqrt{1 - (\frac{n_i}{n_t}^2(1-\cos^2\theta))}}{n_i \cos\theta +n_t\sqrt{1 - (\frac{n_i}{n_t}^2(1-\cos^2\theta))}}=\frac{\cos\theta -\eta\sqrt{1 - (\frac{1}{\eta}^2(1-\cos^2\theta))}}{\cos\theta +\eta\sqrt{1 - (\frac{1}{\eta}^2(1-\cos^2\theta))}}

r_\parallel=\frac{n_i \sqrt{1 - (\frac{n_i}{n_t}^2 sin^2\theta)} -n_t\cos\theta}{n_i \sqrt{1 - (\frac{n_i}{n_t}^2 sin^2\theta)} +n_t\cos\theta}=\frac{n_i \sqrt{1 - (\frac{n_i}{n_t}^2(1-\cos^2\theta))} -n_t\cos\theta}{n_i \sqrt{1 - (\frac{n_i}{n_t}^2(1-\cos^2\theta))} +n_t\cos\theta}=\frac{\sqrt{1 - (\frac{1}{\eta}^2(1-\cos^2\theta))} -\eta\cos\theta}{\sqrt{1 - (\frac{1}{\eta}^2(1-\cos^2\theta))} +\eta\cos\theta}

In hlsl this is:

float FresnelDielectricDielectric(float Eta, float CosTheta)
{
   float SinTheta2 = 1 - CosTheta * CosTheta;

   float t0 = sqrt(1 - (SinTheta2 / (Eta * Eta)));
   float t1 = Eta * t0;
   float t2 = Eta * CosTheta;

   float rs = (CosTheta - t1) / (CosTheta + t1);
   float rp = (t0 - t2) / (t0 + t2);

   return 0.5 * (rs * rs + rp * rp);
}

An equivalent expression can be derive from the [1] formulation by setting k_t=0. Source [7]:

c=\cos \theta

g=\sqrt{\eta^2+ c^2-1}

R=\frac{(g-c)^2}{2(g+c)^2}(1+\frac{(c(g+c)-1)^2}{(c(g-c)+1)^2}

If g is imaginary, this indicates a total internal reflection and R=1.

Note that’s here we get the reflectance directly. In hlsl this is:

float FresnelDielectricDielectric(float Eta, float CosTheta)
{
    float c = CosTheta;
    float temp = Eta* Eta + c * c - 1;

   if (temp < 0)
      return 1;

   float g = sqrt(temp);
   return 0.5 * Square((g - c) / (g + c)) *
                       (1 + Square(( (g + c)  * c - 1) / ((g - c) * c+ 1)));
}

Fresnel Schlick’s approximation

The Fresnel Schlick’s approximation [9] is widely used and knows (Call FresnelSchlick later) but its properties are not.

There is two possible versions of the approximation.

The first Schlick’s approximation use R(0) , the  Fresnel reflectance at 0° (reflectance at normal incidence), also name specular color in PBR game :

R=R(0) + (1 - R(0))(1-\cos \theta)^5

in hlsl:

float3 FresnelSchlick(float3 SpecularColor, float3 E,float3 H)
{
    return SpecularColor + (1.0f - SpecularColor) * pow(1.0f - saturate(dot(E, H)), 5);
}

For conductor:

R(0)=\frac{(n_t - n_i)^2+k_t^2}{(n_t + n_i)^2+k_t^2}=\frac{(\eta - 1)^2+\eta_k^2}{(\eta + 1)^2+\eta_k^2}=\mathrm{Re}{\frac{((n_t +\mathbf{i} k_t) -n_i)((n_t -\mathbf{i} k_t) -n_i) }{((n_t +\mathbf{i} k_t) +n_i)((n_t -\mathbf{i} k_t) +n_i)}}

For dielectric:

R(0)=\frac{(n_t - n_i)^2}{(n_t + n_i)^2}=\frac{(\eta - 1)^2}{(\eta + 1)^2}

In game we always define specular color as R(0) with n_i=1 for air as we always want to deal with air-dielectric/conductor interface. But if we need to interact with another interface like water-dielectric, we need to store another specular color. In case of specular color define for air-dielectric/conductor interface, IOR can be retrieving with:

n_t=\frac{(1+\sqrt{R(0)})}{(1-\sqrt{R(0)})}

However for complex IOR there is no way to retrieve the correct n_t and k_t. So in case of multi layer BRDF prepare to move to parameters with complex IOR instead of specular color.

The FresnelSchlick (Blue) with air-dielectric interface is only accurate in the range 1.4 < n_t < 2.2. Outside this range the relative error grows [8]:

FresnelSchlick

The accurate range still covers most use material in game (water at 1.33 still ok), so no need to worry too much about this.

The other form of the Schlick approximation is without R(0). This is exactly the same as the first one, just replace R(0) by one of the equality provided above:

Example for dielectric-dielectric:

R=\frac{(n_t - n_i)^2+4 n_t n_i(1-\cos \theta)^5}{(n_t + n_i)^2}

Example for dielectric-conductor which is the formulation provide in Lazanyi et al [8].

R=\frac{(n_t - n_i)^2+4 n_t n_i(1-\cos \theta)^5+k_t^2}{(n_t+n_i)^2+k_t^2}

The FresnelSchlick (Blue) for air-conductor don’t match exactly the Fresnel conductor equation (Red):

FresnelMetal

There is a noticeable “dip” just before going to white in case of several metals. Lazanyi et al [8] try to provide a better FresnelSchlick approximation for Metals. They try to add a compensation term -a \cos \theta (1 - \cos \theta)^\alpha to take into account the error introduce by the FresnelSchlick. The term is a rational approximation of the error between Fresnel equation and FresnelSchlick equation. Here is an example with n_t=1.5, n_k=5. Blue curve is the error term, red is the rational approximation with a=3 and \alpha=7:

ErrorFresnelMetal

They suggest to use a = 2 n_t but the \alpha value is more problematic. To understand why, examine p (Red) and s-polarized (Yellow) curves of Fresnel equation for an air-dielectric interface (Left) and an air-conductor interface (Right):

BrewsterAngle

The left figure show that’s the “dip” come from the perpendicular curve. For dielectric, the  perpendicular curve crosses the 0 line at the Brewster’s angle: \theta_b = ArcTan(n_t / n_i) = ArcTan(\eta).
The reflectance R_p does not reach zero for a conductor as it does for a dielectric, but the angle for which it is a minimum is called the pseudo Brewster. However there is no simple solution to retrieve it. That’s why the \alpha is difficult to get right as it should be related to this pseudo-Brewster’s angle which define the local minimun of the error term.

Authors suggest that for most metal local minimun of the error term is somewhere between 0.1 and 0.15. I write a small “solver” in Mathematica to retrieve a good value. It retrieves the local minimum by brute force, then solves the error compensation equation with it. For the author’s case of n_t=1.5 and k_t=5 I get alpha=7.7. Note that the calculated value is not the optimal value it just match the local minimum. Look at the Mathematica file if interested. The equation provided in the paper follow the second Schlick’s approximation form, so I reuse it here. Following graph show the result approximation  (Blue):

R=\frac{(n_t - 1)^2+4n_t(1-\cos \theta)^5+k^2}{(n+1)^2+k^2} -2 n_t \cos \theta (1 - \cos \theta)^\alpha for n_t=1.5, k_t=5 and alpha=7.7

Fresnelschlickcompensation

Result is good, equation is cheap but only work for a particular case and need to be tune for each value (Remember that’s a metal IOR is chromatic and the process need to be done 3 times), so not really useful in practice. Anyway the visual impact of this “dip” is negligible and we can avoid it [8].

When n_i > n_t mean we go to a less dense medium there is  TIR occuring at critical angle. \theta_c=ArcSin(\frac{n_t}{n_i})=ArcSin(\eta). FresnelSchlick don’t handle TIR. This can be fix by using \cos \theta_t instead of \cos \theta_i in the FresnelSchlick when n_i > n_t. Here is a graph showing original Fresnel, FresnelSchlick (both curve an almost merged in brown) and FresnelSchlick with \cos \theta_t (Blue) that I will call FresnelSchlickTIR for n_i=1.33 and n_t=1  (water-air interface):

FresnelSchlickTIR

In hlsl this gives:

float FresnelSchlickTIR(float nt, float ni, float3 n, float3 i)
{
   float R0 = (nt - ni) / (nt + ni);
   R0 *= R0;
   float CosX = dot(n, i);
   if (ni > nt)
   {
      float inv_eta = ni / nt;
      float SinT2 = inv_eta * inv_eta * (1.0f - CosX * CosX);
      if (SinT2 > 1.0f)
      {
         return 1.0f; // TIR
      }
      CosX = sqrt(1.0f - SinT2);
   }

   return R0 + (1.0f - R0) * pow(1.0 - CosX, 5.0);
}

With so much instruction the interest for FresnelSchlick start to decreases.

For an optimized FresnelSchlick formulation see Spherical Gaussian approximation for Blinn-Phong, Phong and Fresnel.

Reflect/Refract vector

Reflect and refract vector are derive from Snell ‘s law [4]:

Reflect:

float3 reflect( float3 i, float3 n )
{
  return i - 2.0 * n * dot(n, i);
}

Refract vector is calculated as:

t= \frac{n_i}{n_t}i+(\frac{n_i}{n_t}\cos \theta_i-\sqrt{1-\sin^2\theta_t}) n

\sin^2\theta_t= \frac{n_i}{n_t}^2(1-\cos^2 \theta_i)

\frac{n_i}{n_t}=\frac{1}{\eta}

The derivation for refraction vector can be found in [4].

float3 refract( float3 i, float3 n, float inv_eta )
{
  float cosi = dot(-i, n);
  float cost2 = 1.0f - inv_eta * inv_eta * (1.0f - cosi*cosi);
  float3 t = inv_eta * i + ((inv_eta * cosi - sqrt(abs(cost2))) * n);
  return t * (float3)(cost2 > 0);
}

Optimized version of refract can be done for specific IOR, see Water drop 3a – Physically based wet surfaces.

Mathematica

All the content of this post has been check in Mathematica. Here are the Mathematica files, one pdf: Fresnel and the notebook : Fresnel_nb (Right click, save, rename “.pdf” to “.nb” as wordpress don’t support zip files).

Reference

[1] Shirley, “Physically based lighting calculations for computer graphics”,  http://www.cs.virginia.edu/~jdl/bib/globillum/shirley_thesis.pdf
[2] Wikipedia, “Fresnel equation”, http://en.wikipedia.org/wiki/Fresnel_equations
[3] Pharr, Humphreys, “Physically based rendering book”, http://www.pbrt.org/
[4]Greve, “Reflections and Refractions in Ray Tracing”,  http://graphics.stanford.edu/courses/cs148-10-summer/docs/2006–degreve–reflection_refraction.pdf
[5] Möller, “Optics” book, http://www.springer.com/physics/optics+%26+lasers/book/978-0-387-26168-3
[6] Siegel, Howell, “Thermal radiation hit transfer, 3rd edition” book, http://www.amazon.com/Thermal-Radiation-Heat-Transfer-5th/dp/1439805334
[7] Larsen,  “Shading a Surface using BRDFs”, http://www2.imm.dtu.dk/~jab/shaderday/ShadingASurfaceUsingBRDFs.ppt
[8] Lazanyi, Szirmay-Kalos, “Fresnel term approximations for Metals” http://sirkan.iit.bme.hu/~szirmay/fresnel.pdf
[9] Schlick, “An inexpensive BRDF model for physically based rendering”, http://www.cs.virginia.edu/~jdl/bib/appearance/analytic%20models/schlick94b.pdf

About these ads

8 Responses to Memo on Fresnel equations

  1. seblagarde says:

    V1.1 : Minor change in Reflect/refract vector section

  2. SebH says:

    Nice! Good to see all these data gathered in one place!
    Tiny detail concerning the Schilck approximation: the pow(x,5) can be translated to x2=x*x; x5=x2*x2*x; This was always faster on gen3 in the use cases I have encountered. :)

  3. seblagarde says:

    v1.2 : Add 3 hlsl code for Dielectric-Conductor interface and Dieletric-Dieltricinterface

  4. seblagarde says:

    v1.3: Add a statement for the last form of exact fresnel equation:
    If g is imaginary, this indicates a total internal reflection and R=1

    And update the schlick approximation with the other formulation without R(0):

    Example for dielectric-dielectric:

    R=\frac{(n_t - n_i)^2+4 n_t n_i(1-\cos \theta)^5}{(n_t + n_i)^2}

    Example for dielectric-conductor which is the formulation provide in Lazanyi et al [8].

    R=\frac{(n_t - n_i)^2+4 n_t n_i(1-\cos \theta)^5+k_t^2}{(n_t+n_i)^2+k_t^2}

  5. Mika Ishikawa says:

    Thank you for the great memo, it is really really helpful! I just would like to confirm is it ok to reference your code in our own projects or is there a copyright restriction or other similar restrictions?

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: