# Bringing an Accurate Fresnel to Real-Time Rendering¶

#### Laurent Belcour (Unity), Megane Bati (IOGS), and Pascal Barla (Inria)¶

This notebook is the implementation & details guide for Bringing an Accurate Fresnel to Real-Time Rendering: A Preintegrable Decomposition. It provides all the details to reproduce our code and results for off-line rendering. For real-time rendering, we also provide a reference code to compute the directional albedo. Please refer to that code to reproduce the new split-sum integral form.

You can freely copy/paste the python scripts here to reproduce our work. To work, our code requires that you have the python package below installed. You can install those packages using pip for example.

Note that this notebook requires that you download our Fresnel python package fresnel_tools.py.

In [21]:
# %matplotlib inline
%matplotlib notebook
import numpy as np
import scipy.optimize
from scipy.interpolate import interp2d
import matplotlib as mpl
from matplotlib.ticker import MaxNLocator
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from matplotlib.widgets import Slider
from mpl_toolkits.mplot3d import Axes3D

# Custom package
import fresnel_tools as ft

mpl.rcParams['font.family']=['serif']
mpl.rcParams['font.serif']=['Linux Libertine'] + mpl.rcParams['font.serif']


## A Basis Decomposition for Fresnel¶

In this work, we which to reproduce Fresnel reflectance_reflection_and_transmission_coefficients) using a linear decomposition of data driven basis elements. That is, if the unpolarized reflected power of a planar surface (denoted Fresnel reflectance) is given by $F(\theta)$, we will express it as: $$F(\theta) = \sum_i c_i b_i(\theta),$$ where $b_i(\theta)$ are basis functions and $c_i$ are coefficients.

Note that Fresnel reflectance is a wavelength depent function $F(\theta; \lambda)$. In such a case, we will incorporate wavelength dependence in the coefficients $c_i$, but not in the basis elements $b_i(\theta)$: $$F(\theta; \lambda) = \sum_i c_i(\lambda) b_i(\theta).$$ In such form, the basis elements can be factored out of spectral preintegration (when one whish to extract RGB reflectance out of spectral data): $$\int_\lambda F(\theta; \lambda) s(\lambda) \mbox{d}\lambda = \sum_i \left[ \int_\lambda c_i(\lambda) s(\lambda) \mbox{d}\lambda \right] b_i(\theta).$$

## How Many Basis Elements?¶

In this Section, we analyze how many basis elements are required to reproduce Fresnel curves with high acuracy. To do so, we extract the eigenvalues of a set of Fresnel curves and plot their relative value to the sum of all eigenvalues.

In [22]:
'''
Compute the SVD for Fresnel curves and extract the eigenvalues

'''
nb_points = 1000
nb_IOR    = 100

'''
We define the curves with respect to the elevation.
For real-time implementation, it makes more sense to define
since it avoids an acos call.
'''
cT        = np.linspace(0.0, 1.0, nb_points, endpoint=True)
T         = np.arccos(cT)

subrange  = np.linspace(0, nb_points-1, 20, dtype=int)
IOR       = np.linspace(1.01, 3.0, nb_IOR)
KAPPA     = np.zeros((nb_IOR, ))
curves    = np.ndarray( (nb_IOR, nb_points) )

for k,ior in enumerate(IOR):
F  = ft.fresnel(cT, ior, KAPPA[k])
R0 = ft.fresnel(1.0, ior, KAPPA[k])
curves[k,:] = (F - R0) / (1.0-R0)

u, s, vh = np.linalg.svd(curves)
eigv = s[0:20] / np.sum(s)

'''
Display the progression of eigenvalues
'''
fig, axs = plt.subplots(1, 2, figsize=(9,4))
handles = []
labels  = []

axs[0].grid()
axs[0].axis([0, 0.5*np.pi, 0, 1.0])
for k in np.linspace(0, 10, endpoint=False, dtype=int):
handle, = axs[0].plot(T, curves[k, :])
axs[0].set_ylabel('Reflectance')
axs[0].set_xlabel(r'$\theta$ in radians')

axs[1].xaxis.set_major_locator(MaxNLocator(integer=True))
axs[1].grid()
axs[1].axis([0, 19, eigv[0]-0.01, 1.01])
handle2, = axs[1].plot([0, 19], [0.99, 0.99])
handle1, = axs[1].plot(np.cumsum(eigv), marker='o')
axs[1].set_ylabel('Cumulative eigenvalues')
axs[1].set_xlabel(r'Number of basis elements')

plt.tight_layout()
plt.show()


As we can see, for dieletrics materials with $\eta \in [1.01, 3.0]$, only two basis elements different from the constant are required to correctly capture $99\%$ of the Fresnel reflectance curves' eigenvalues.

## Compute the Decomposition¶

In this Section, we compute our decomposition using a specific set of Index of Refractions (IORs). We used a thousand of different curves of both dielectric materials with $\eta_2 > \eta_1$ and conductors (spanned using the reflectivity / edge-tint parameterization) to build our decomposition.

In [23]:
'''
Compute the basis for Fresnel and define the projection operator.

'''
nb_points = 1000
nb_IOR    = 1000

'''
Gulbrandsen conversion routines

Gulbrandsen. Artist Friendly Metallic Fresnel. JCGT 2014.
http://jcgt.org/published/0003/04/03/paper.pdf

More details on it later.
'''
def RGtoNK2(r, g):
if r == 0.0:
return 1.0, 0.0
n = g*(1.0-r)/(1.0+r) + (1.0-g)*(1.0+np.sqrt(r))/(1.0-np.sqrt(r))
k2 = 0.0 if g<=0.0 else (r*(1.0+n)*(1.0+n) - (n-1.0)*(n-1.0)) / (1.0-r)
return n, k2

def RGtoNK(r, g):
n, k2 = RGtoNK2(r, g)
return n, np.sqrt(k2)

'''
We define the curves with respect to the elevation.
For real-time implementation, it makes more sense to define
since it avoids an acos call.
'''
cT        = np.linspace(0.0, 1.0, nb_points, endpoint=True)
T         = np.arccos(cT)

subrange  = np.linspace(0, nb_points-1, 20, dtype=int)
IOR       = np.linspace(1.1, 3.0, nb_IOR)
KAPPA     = np.zeros((nb_IOR, ))
curves    = np.ndarray( (IOR.shape[0], nb_points) )

'''
We use a small percentage (1%) of conducting materials (g=1)
and a large percentage (99%) of dieletrics to build the basis.
This helps to better reconstruct both of them.
'''
alpha = 0.01
for k in range(1, nb_IOR):
r = k / nb_IOR

if r <= alpha:
r = r/alpha
r = (r) / (1.0-1.0E-5) + 1.0E-5
g = 1.0
IOR[k], KAPPA[k] = RGtoNK(r, g)
else:
r = (r - alpha)/(1.0-alpha)
IOR[k], KAPPA[k] = 1.01+(3.0-1.01)*r,0.0

'''
Build the Fresnel reflectance curves from the set of IOR
set above. We remove F(0deg) from the curves and rescale
it to map between 0 to 1.
'''
for k,ior in enumerate(IOR):
F  = ft.fresnel(cT, ior, KAPPA[k])
R0 = ft.fresnel(1.0, ior, KAPPA[k])
curves[k,:] = (F - R0) / (1.0-R0)

u, s, basis = np.linalg.svd(curves)

# Define the zeroth and first componnents
basis0 = np.ones_like(basis[0, :])
basis1 = np.abs(basis[0, :])
basis1 = basis1 / max(basis1)

for k,ior in enumerate(IOR):
F  = ft.fresnel(cT, ior, KAPPA[k])
R0 = ft.fresnel(1.0, ior, KAPPA[k])
curves[k,:] = (F - R0)
curves[k,:]  = curves[k,:]/(1.0-R0) - basis1

u, s, basis = np.linalg.svd(curves)

# Define the second and third componnents
basis2 = basis[0] / np.abs(basis[0]).max()
basis3 = basis[1] / np.abs(basis[1]).max()

i1  = np.argmax( basis2 )
t1  = T[ i1 ]
ct1 = cT[ i1 ]
print('Critical angle for basis2 = ', ct1, t1, 180.0/np.pi*t1, i1)

i2  = np.argmax( +basis3 )
i3  = np.argmax( -basis3 )
i2  = i2 if abs(i2-i1) > abs(i3-i1) else i3
t2  = T[ i2 ]
ct2 = cT[ i2 ]
print('Critical angle for basis3 = ', ct2, t2, 180.0/np.pi*t2, i2)

'''
Define the projection operator that maps a complex IOR (eta+i*kappa)
to the coefficients of the decomposition.
'''
def fresnel2coeffs(eta12, kappa=0.0):
coeffs    = np.zeros((4))
coeffs[0] = ft.fresnel(1.0, eta12, kappa)
coeffs[1] = 1.0 - coeffs[0] if coeffs[0] else 0.0

dF  = np.zeros((1,2))
dF[0,0]  = ft.fresnel(ct1, eta12, kappa) - (coeffs[0] + coeffs[1]*basis1[i1])
dF[0,1]  = ft.fresnel(ct2, eta12, kappa) - (coeffs[0] + coeffs[1]*basis1[i2])

F12 = np.zeros((2,2))
F12[0,0] = basis2[i1]
F12[0,1] = basis2[i2]
F12[1,0] = basis3[i1]
F12[1,1] = basis3[i2]
A = np.linalg.inv(F12)

coeffs[2:4] = np.dot(dF, A)
return coeffs

'''
Display the basis elements and range of value for the coefficients for various
values of eta and a dieletric material.
'''
plt.subplots(figsize=(9,4))
handles = []
labels  = []

plt.grid()
handle0, = plt.plot(T, basis0)
handle1, = plt.plot(T, basis1)
handle2, = plt.plot(T, basis2)
handle3, = plt.plot(T, basis3)
plt.legend([handle0, handle1, handle2, handle3], [r'$b_0(\theta)$', r'$b_1(\theta)$', r'$b_2(\theta)$', r'$b_3(\theta)$'])
plt.xlabel(r'angle of incidence $\theta$')
plt.axis([0, 0.5*np.pi, -0.9, 1.1])
plt.tight_layout()
plt.show()
# plt.savefig('../figures/basis.pdf')

Critical angle for basis2 =  0.12812812812812813 1.4423150071582151 82.63856263854684 128
Critical angle for basis3 =  0.43243243243243246 1.123607581568183 64.37797225275827 432

In [24]:
'''
Compare the basis elements with the one of Shlick and Lazanyi.
'''
plt.subplots(ncols=1, figsize=(9,4))
handles   = []
labels    = []

plt.grid()
plt.title('basis for $F(\\theta)$')
handle, = plt.plot(T, basis1, color='C1')
handles += [handle]
labels  += ['Ours vs  Shlick']
plt.plot(T, np.power(1.0-cT, 5.0),  dashes=[1, 3], color='C1')

handle, = plt.plot(T, basis2)
handles += [handle]
labels  += ['Ours vs  Lazanyi']
lazanyi = cT * np.power(1.0-cT, 6.0)
plt.plot(T, lazanyi/lazanyi.max(), color=handle.get_color(), dashes=[1, 3])
plt.legend(handles, labels, loc='upper left')
plt.axis([0, 0.5*np.pi, 0.0, 1.0])
plt.tight_layout()
plt.show()


## Reconstruction of Fresnel Curves¶

In this Section, we use the decomposition defined above to fit different Fresnel curves. We start with dielectric materials, and then decompose conductors. Interestingly, and even though we did not use conductors to compute the decomposition, conductors are also well approximated by the decomposition.

In [25]:
'''
Display the fitting of dieletric reflectances
'''
plt.subplots(figsize=(9,4))
plt.grid()
handles = []
labels = []

for k,ior in enumerate([1.0, 1.01, 1.1, 1.5, 2.0]):
F = ft.fresnel(cT, ior)
handle, = plt.plot(T, F)
handles.append(handle)

coeffs = fresnel2coeffs(ior)
P = coeffs[0] + coeffs[1]*basis1 + coeffs[2]*basis2 + coeffs[3]*basis3
plt.scatter(T[subrange], P[subrange], color=handle.get_color(), marker='o')
plt.plot(T, P, dashes=[1, 3], color=handle.get_color())

L2 = np.linalg.norm(F-P)
labels.append('IOR = ' + str(ior) + ', L2 = ' + f'{L2:.2f}')

plt.legend(handles, labels, loc='upper left')
plt.title('Reconstruct dielectrics $F(\\theta)$ with our decomposition')
plt.axis([0, 0.5*np.pi, -0.01, 1.0])
plt.tight_layout()
plt.show()

In [26]:
'''
Display the fitting of conducting reflectances
'''
plt.subplots(figsize=(9,4))
plt.grid()

Etas   = [1.90, 0.46, 0.15, 2.9, 0.29]
Kappas = [6.79, 2.97, 3.47, 3.0, 2.86]
Titles = ['Al', 'Cu', 'Ag', 'Fe', 'Au']

handles = []
labels = []
for k,eta12 in enumerate(Etas):
kappa = Kappas[k]
F = ft.fresnel(cT, eta12, kappa)
handle, = plt.plot(T, F)
handles.append(handle)

coeffs = fresnel2coeffs(eta12, kappa)
P = coeffs[0] + coeffs[1]*basis1 + coeffs[2]*basis2 + coeffs[3]*basis3
plt.scatter(T[subrange], P[subrange], color=handle.get_color(), marker='o')
plt.plot(T, P, dashes=[1, 3], color=handle.get_color())

L2 = np.linalg.norm(F-P)
labels.append(Titles[k] + r', $\eta = ' + str(eta12) + r' + i ' + str(kappa) + '$, L2 = ' + f'{L2:.2f}')

plt.legend(handles, labels, loc=(0.02, 0.2))
plt.axis([0, 0.5*np.pi, 0.45, 1.0])
plt.title('Reconstruct conductors $F(\\theta)$ with our decomposition')
plt.tight_layout()
plt.show()
# plt.savefig('../figures/reconstruct_conductors.pdf')

In [28]:
'''
Interactive editing of the reflectance curve using the coefficients
'''
# Params for sliders
HSLIDERS = 2; VSLIDERS = 1
INIT_VAL = [1.1, 0.0]; MIN_VAL = [1.01, 0.0]; MAX_VAL = [4.0, 4.0]
COLOR = ["lightgrey", "mediumaquamarine"]
LABELS = ["$\\eta$", "$\\kappa$"]
LENGTH, X, Y, BOTTOM = ft.init(HSLIDERS, VSLIDERS)

# Start plot
START = 0.0; #START = np.pi/4.0-0.001

fig, ax = plt.subplots(figsize=(9,4))
plt.grid()

plt.axis([START, np.pi/2.0, 0, 1])
# plt.title("Fresnel along the isoline")
plt.xlabel("incidence angle $\\theta$")
plt.ylabel("Reflectance $F(\\theta)$")
plt.tight_layout()

# Some interesting points
P = ft.fresnel(cT, INIT_VAL[0], INIT_VAL[1])
fresnel, = plt.plot(T, P, 'tab:blue', linestyle='--')
coeffs = fresnel2coeffs(INIT_VAL[0], INIT_VAL[1])
P = coeffs[0] + coeffs[1]*basis1 + coeffs[2]*basis2 + coeffs[3]*basis3
points,  = plt.plot(T, P, 'tab:blue')
P = coeffs[0] + (1.0-coeffs[0])*np.power(1.0-cT, 5.0)
schlick, = plt.plot(T, P, 'orange')
plt.legend(['Reference Fresnel', 'Our Decomposition', 'Schlick\'s Fresnel'])

sliders = []
sliders.append(ft.add_slider(X[0], Y[0], INIT_VAL[0], LABELS[0], COLOR[0], LENGTH, MIN_VAL[0], MAX_VAL[0]))
sliders.append(ft.add_slider(X[0], Y[1], INIT_VAL[1], LABELS[1], COLOR[1], LENGTH, MIN_VAL[1], MAX_VAL[1]))

def update_coeffs(val):
n = sliders[0].val
k = sliders[1].val

P = ft.fresnel(cT, n, k)
fresnel.set_ydata(P)

coeffs = fresnel2coeffs(n, k)
P = coeffs[0] + coeffs[1]*basis1 + coeffs[2]*basis2 + coeffs[3]*basis3
points.set_ydata(P)

P = coeffs[0] + (1.0-coeffs[0])*np.power(1.0-cT, 5.0)
schlick.set_ydata(P)

fig.canvas.draw_idle()

sliders[0].on_changed(update_coeffs)
sliders[1].on_changed(update_coeffs)

plt.show()


## Gulbransen Parametrization¶

The Gulbransen parametrization map a couple $(r, g) \in [0,1]^2$ denoted reflectivity ($r$) and edge-tint ($g$) to the complex Index of Refraction $\eta + i\,\kappa$. This parametrization tries to provide an artistic mapping in a closed range that behave more linearly than the complex IOR. The parametrization is the following:

\begin{align} \eta &= g \dfrac{1-r}{1+r} + \left(1-g\right)\dfrac{1+\sqrt{r}}{1-\sqrt{r}} \\ \kappa &= \sqrt{\dfrac{1}{1-r}\left( r \, \left( 1+n \right)^2 - \left( n-1 \right)^2 \right)} \end{align}
In [29]:
'''
Interactive plot of Gulbransen's parametrization.
'''
R = np.linspace(0.0, 1.0, nb_IOR)
G = np.linspace(0.0, 1.0, nb_IOR)

# Params for sliders
HSLIDERS = 2; VSLIDERS = 1
INIT_VAL = [0.5, 0.0, 0.0, 0.0]; MIN_VAL = [0.0, 0.0]; MAX_VAL = [1.0, 1.0]
COLOR = ["lightgrey", "mediumaquamarine"]
LABELS = ["$r$", "$g$"]
LENGTH, X, Y, BOTTOM = ft.init(HSLIDERS, VSLIDERS)

# Start plot
START = 0.0; #START = np.pi/4.0-0.001

fig, ax = plt.subplots(figsize=(9,4))
plt.grid()
plt.axis([START, np.pi/2.0, 0, 1])
plt.xlabel(r"incidence angle $\theta$")
plt.ylabel(r"Reflectance $F(\theta)$")
plt.tight_layout()

# place a text box in upper left in axes coords
textstr = r'$\eta = 0$, $\kappa = 0$'
textbox = ax.text(0.01, 0.06, textstr, transform=ax.transAxes, fontsize=10, verticalalignment='top')
textstr2 = r'$c_0 = 0$, $c_1 = 0$, $c_2 = 0$, $c_3 = 0$'
textbox2 = ax.text(0.59, 0.06, textstr, transform=ax.transAxes, fontsize=10, verticalalignment='top')

# Some interesting points
coeffs = INIT_VAL
P = 0.0*basis3
points,  = plt.plot(T, P, linewidth=2)
refs,    = plt.plot(T, P, dashes=[1, 3], color='k')
schlick, = plt.plot(T, P)

labels = ['Ours', 'Reference', 'Schlick']
plt.legend([points, refs, schlick], labels, loc='upper left')

def update_rg(val):
r, g = sliders[0].val, sliders[1].val
n, k = RGtoNK(r, g)
str_n = f"{n:.2f}"
str_k = f"{k:.2f}"
textstr = r'$\eta = ' + str_n + r'$, $\kappa = ' + str_k + '$'
textbox.set_text(textstr)

coeffs = fresnel2coeffs(n, k)
resOurs = coeffs[0] + coeffs[1]*basis1 + coeffs[2]*basis2 + coeffs[3]*basis3
points.set_ydata(resOurs)
str_c0 = f"{coeffs[0]:.2f}"
str_c1 = f"{coeffs[1]:.2f}"
str_c2 = f"{coeffs[2]:.2f}"
str_c3 = f"{coeffs[3]:.2f}"
textstr = r'$c_0 = ' + str_c0 + r',$c_1 = ' + str_c1 + r', $c_2 = ' + str_c2 + r',$c_3 = ' + str_c3
textbox2.set_text(textstr)

res = ft.fresnel(cT, n, k)
refs.set_ydata(res)

resSchlick = sliders[0].val + (1.0-sliders[0].val)*np.power(1.0-cT, 5.0)
schlick.set_ydata(resSchlick)

fig.canvas.draw_idle()

sliders = []
sliders.append(ft.add_slider(X[0], Y[0], INIT_VAL[0], LABELS[0], COLOR[0], LENGTH, MIN_VAL[0], MAX_VAL[0]))
sliders.append(ft.add_slider(X[0], Y[1], INIT_VAL[1], LABELS[1], COLOR[1], LENGTH, MIN_VAL[1], MAX_VAL[1]))
sliders[0].on_changed(update_rg)
sliders[1].on_changed(update_rg)
update_rg(0.0)

plt.show()


In the following, we study the effect of the Gulbrandsen's parametrization on the coefficients. Since the coefficients now behave linearly, we can use them to indicate if Gulbransen's mapping behave linearly for Fresnel reflectances. To do so, we are going to visualize the parameterization is the space of fitted coefficients $(c_0, c_2)$. Since $c_0$ is controlled by the reflectivity only, the spacing between points on the vertical $(c_2)$ axis will depict non-linear behaviours.

But first, we introduce a remapping of the $g$ parameter that is more linear: \begin{align} \bar{g} = 1 - (1-g)^3 \end{align}

While we won't need this parametrization in the end but it will help to avoid aliasing issues when building our improved parametrization.

In [30]:
fig, axs = plt.subplots(1,2, figsize=(8,4))

nb_POINTS = 100
nb_RG = 50
NEAR_ONE = 1.0 - 1.0E-5

'''
Generate two lines: one for the perfect dielectrics (g=0),
and one for the maximum edge tint (g=1).
'''
R = np.linspace(0.00, NEAR_ONE, nb_POINTS)
c_diel = np.zeros((nb_POINTS, 4))
c_cond = np.zeros((nb_POINTS, 4))
for i,r in enumerate(R):
[eta,kappa] = RGtoNK(r,0.0)
c_diel[i,:] = fresnel2coeffs(eta, 0.0)
[eta,kappa] = RGtoNK(r,1.0)
c_cond[i,:] = fresnel2coeffs(eta, kappa)

'''
Generate a grid of point, that provide the meshing of the
space between the two curves using Gulbransen's
parametrization.
'''
R = np.linspace(0.00, NEAR_ONE, nb_RG)
G = np.linspace(0.00, NEAR_ONE, nb_RG)
coeffs = np.zeros((nb_RG, nb_RG, 4))
for i, r in enumerate(R):
for j, g in enumerate(G):
[eta,kappa] = RGtoNK(r,g)
coeffs[i,j,:] = fresnel2coeffs(eta, kappa)

ax = axs[0]
ax.scatter(coeffs[:,:,0], coeffs[:,:,2], c=coeffs[:,:,0], s=0.1)
ax.plot(c_diel[:,0], c_diel[:,2], c='red')
ax.plot(c_cond[:,0], c_cond[:,2], c='k')
# ax.plot(c_zero[:,0], c_zero[:,2], c='gray')
ax.set_xlabel('$c_0$')
ax.set_ylabel('$c_2$')
ax.grid()
ax.axis([0,1, -0.35, 0.35])
ax.legend(['$g=0$', '$g=1$'], loc='upper right')
ax.set_title('Gulbransen parametrization')

'''
Generate a grid of point, that provide the meshing of the
space between the two curves using a remapped Gulbransen's
parametrization.
'''
R = np.linspace(0.00, NEAR_ONE, nb_RG)
g = np.linspace(0.00, NEAR_ONE, nb_RG)
coeffG = np.zeros((nb_RG, nb_RG, 4))
for i,r in enumerate(R):
for j,g in enumerate(G):
# empirical remapping
g = 1.0-pow(1.0-g,3)
[eta,k2] = RGtoNK2(r,g)
coeffG[i,j,:] = fresnel2coeffs(eta, np.sqrt(k2))

ax = axs[1]
ax.scatter(coeffG[:,:,0], coeffG[:,:,2], c=coeffG[:,:,0], s=0.1)
ax.plot(c_diel[:,0], c_diel[:,2], c='red')
ax.plot(c_cond[:,0], c_cond[:,2], c='k')
ax.set_xlabel('$c_0$')
ax.set_ylabel('$c_2$')
ax.grid()
ax.axis([0,1, -0.35, 0.35])
ax.legend(['$g=0$', '$g=1$'], loc='upper right')
ax.set_title('Remapping using $\\bar{g}$')

plt.tight_layout()
plt.show()


This remapping of the Gulbransen parametrization is still not satisfactory: we can see compressions close to the $g = 0$ and $g= 1$ curves.

Such compressions are also visible when visualizing the manifold in 3D.

In [31]:
fig = plt.figure(figsize=(9,4))

nb_RG = 20

'''
Generate the manifold of the coefficients using Gulbransen's
parametrization.
'''
R = np.linspace(0.0, NEAR_ONE, nb_RG)
G = np.linspace(0.0, 1.0, nb_RG)
coeffs = np.zeros((nb_RG, nb_RG, 4))
for i, r in enumerate(R):
for j, g in enumerate(G):
[eta,kappa] = RGtoNK(r,g)
coeffs[i,j,:] = fresnel2coeffs(eta, kappa)

ax.plot_surface(coeffs[:,:,0], coeffs[:,:,2], coeffs[:,:,3])
ax.plot(c_diel[:,0], c_diel[:,2], c_diel[:,3], 'r')
ax.plot(c_cond[:,0], c_cond[:,2], c_cond[:,3], 'k')
ax.view_init(elev=65., azim=-45)

ax.set_xlabel('$c_0$')
ax.set_ylabel('$c_2$')
ax.set_zlabel('$c_3$')
ax.set_xlim(0,1)
ax.set_ylim(-0.4,0.4)
ax.set_zlim(-1,1)
ax.set_title('Gulbransen parametrization')

'''
Generate two lines: one for the perfect dielectrics (g=0),
and one for the maximum edge tint (g=1).
'''
R = np.linspace(0.0, NEAR_ONE, nb_POINTS)
c_diel = np.zeros((nb_POINTS, 4))
c_cond = np.zeros((nb_POINTS, 4))
for i,r in enumerate(R):
[eta,kappa] = RGtoNK(r,0.0)
c_diel[i,:] = fresnel2coeffs(eta, 0.0)
[eta,kappa] = RGtoNK(r,1.0)
c_cond[i,:] = fresnel2coeffs(eta, kappa)

'''
Generate the manifold of the coefficients using the updated
parametrization.
'''
R = np.linspace(0.0, NEAR_ONE, nb_RG)
g = np.linspace(0.0, 1.0, nb_RG)
coeffG = np.zeros((nb_RG, nb_RG, 4))
for i,r in enumerate(R):
for j,g in enumerate(G):
# empirical remapping
g = 1.0-pow(1.0-g,3)
[eta,k2] = RGtoNK2(r,g)
coeffG[i,j,:] = fresnel2coeffs(eta, np.sqrt(k2))

ax.plot(c_diel[:,0], c_diel[:,2], c_diel[:,3], 'r')
ax.plot(c_cond[:,0], c_cond[:,2], c_cond[:,3], 'k')
ax.plot_surface(coeffG[:,:,0], coeffG[:,:,2], coeffG[:,:,3])
ax.view_init(elev=65., azim=-45)

ax.set_xlabel('$c_0$')
ax.set_ylabel('$c_2$')
ax.set_zlabel('$c_3$')
ax.set_xlim(0,1)
ax.set_ylim(-0.4,0.4)
ax.set_zlim(-1,1)
ax.set_title('Remapping using $\\bar{g}$')

plt.show()