Article Index

Les modèles de matériaux réalistes les plus courants de nos jours sont les modèles basés sur la théorie des micro-facettes. Les matériaux de type Disney font partie de cette catégorie. L'idée est de penser une surface comme un ensemble de particules plus ou moins petites, ayant chacun une normale différente (celles d'un miroir seront homogènes tandis que celles du carton ou du béton ne le sont pas) donnant un aspect plus ou moins rugueux. On peut définir cette répartition des normales par une loi de distribution, et calculer les phénomènes qui interviennent, notamment le masquage et l'ombrage. L'angle de vue du matériau joue aussi un rôle crucial et est pris en compte par le modèle Disney.

Objectif

L'objectif est donc d'implémenter les matériaux physiquement réaliste a base de micro-facette : Plastic, Metal, Substrate et Disney. On peut ainsi voir dans les figures suivantes les effets de la variation de rugosité sur le rendu avec un accroissement de la tache spéculaire différente selon le modèle :

  • Sphères avec un matériau de type Substrate (du plus rugueux au plus lisse).
  • Sphères avec un matériau de type Plastic (du plus rugueux au plus lisse).

 

BDSF

En rendu, la BSDF (Bidirectional Scattering Distribution Function) est la fonction déterminant la probabilité qu’un rayon lumineux soit dispersé dans une direction souhaitée. Celle-ci peut se distinguer en deux types de BxDF : La BRDF qui prend en compte les interactions à l’impact du rayon et la BTDF qui simule les rayons traversant le matériaux. Nous nous concentrerons surtout sur les interactions à la surface. La BSDF peut ainsi se modéliser par la somme de deux types de réaction à la surface : une interaction diffuse qui va répartir la lumière dans toutes les directions et une interaction spéculaire qui va concentrer une grande partie de l’énergie dans une seule direction. On peut modéliser cette fonction ainsi :

$bsdf(x, w_i, w_o) = bsdf_{diff} + bsdf_{spec}$

Pour la suite ont utiliseras les symboles suivant :

$x$ Position de l'intersection sur la scène
$w_i$ La direction du rayon lumineux
$w_o$ La direction du point vers la camera
$n_x$ La normal de la surface en x
$h_x$ Le demi-vecteur à la surface en x
$c_{d,x}$ La couleur diffuse de la surface en x
$c_{s,x}$ La couleur spéculaire de la surface en x
$\alpha_x$ La rugosité de la surface en x
$n_x$ L'indice de réfraction de la surface en x (plastique)
$n_x + ik_x$ L'indice de réfraction complexe de la surface en x (métal)

Modèle de Cook-Torrance

L’interaction spéculaire à la surface étant dépendante de sa géométrie, la mise en place d’un modèle physiquement réaliste va donc passer par le choix d’une BSDF spéculaire représentant ce type de surface. Le modèle de Cook-Torrance [3] nous permet de résoudre ce problème :

$\rho_{ct}(x, w_i, w_o) = \frac{D \ast G \ast F}{4 \ast (w_o \cdot n_x) \ast (w_i \cdot n_x)}$

On voit que ce modèle repose sur trois fonctions : D, G et F. La fonction D est la fonction de distribution des normales à la surface dans une direction donnée. Par exemple si 50% des normales sont dans la direction $\lambda$ alors $D(\lambda)$ vaut $0,5$. De nombreuses versions de cette fonction existent dans la littérature, mais nous utiliserons celle de Trowbridge-Reitz :

$D(h_x, n_x, \alpha_x) = \frac{\alpha^2 X(h_x \cdot n_x)}{\pi((h_x \cdot n_x)^2(\alpha_x^2 +tan^2\frac{1-(w \cdot h_x)^2}{(w \cdot h_x)^2}}$

De plus ont pose la fonction $X$ défini par $X(x) = max(0,min(x,1))$ La fonction G quant à elle définit l’atténuation lumineuse à la surface due à l’auto-occultation des microfacettes. Cette fonction est donc fortement liée à la fonction de distribution des normales à la surface. Nous allons ici aussi utiliser celle de Trowbridge-Reitz :

$G(w, h_x, n_x, \alpha_x) = \frac{2 X( \frac{w \cdot h_x}{w \cdot n_x}) } {1+ \sqrt{1+\alpha_x^2 \frac{1-(w \cdot h_x)^2}{(w \cdot h_x)^2} } }$

Ont a donc $G(w_i, w_o, h_x, n_x, \alpha_x) = G_p(w_o, h_x, n_x, \alpha_x) G_p(w_i, h_x, n_x, \alpha_x)$ La dernière fonction F décrit les phénomènes de réfraction-réflexion à la surface selon le coefficient de Fresnel. Ce coefficient s’exprimant comme la moyenne des réflectances des deux plans de polarisation (parallèle et perpendiculaire) du rayon incident :

$R = \frac{R_s+R_p}{2}$

Avec $R_s$ la réflectance selon le plan perpendiculaire et $R_p$ celle selon le plan parallèle.


BSDF de type plastique

Comme son nom l’indique, cette BSDF tend à modéliser des interactions lumineuses similaires à celle d’un matériau plastique. La BSDF de type Plastique se base sur un modèle de matériau diélectrique standard. Il sera donc modélisé avec un modèle diffusion Lambertien et un spéculaire basé sur un Fresnel diélectrique. Pour faire cela nous allons donc user de la fonction de diffusion Lambertienne et de la fonction spéculaire de Cook-Torrance.

$diff_{Lambert}(x, w_i, w_o) = \frac{c_{d,x}}{\pi}$

Les matériaux plastiques étant de type Diélectrique, le fresnel utilisé par le modèle de Cook-Torrance seradonc en adéquation avec ce type de matériau :

$R_s(n_x, \theta) = \frac{\cos \theta - n_x \sqrt{1- (\frac{1}{n_x}^2(1-\cos^2 \theta))}} {\cos \theta - n_x \sqrt{1+ (\frac{1}{n_x}^2(1-\cos^2 \theta))} }$

$R_p(n_x, \theta) = \frac{\sqrt{1- (\frac{1}{n_x}^2(1-\cos^2 \theta))} - n_x \cos \theta } {\sqrt{1- (\frac{1}{n_x}^2(1-\cos^2 \theta))} + n_x \cos \theta }$

 avec $\theta$ l’angle entre le rayon et la micro-normale à la surface. Dans le cas du calcul du terme spéculaire, on utilisera le demi-vecteur $h_x$.

BSDF de type métal

La BSDF Métallique quant à elle représente des objets conducteurs, elle ne possède donc pas de terme diffus et son spéculaire suit un Fresnel conducteur. Dans le cas d’une surface métallique, l’interaction lumineuse est complètement spéculaire. Nous allons donc seulement utiliser le modèle de Cook-Torrance pour représenter ce type de surface. Contrairement aux matériaux plastiques, le métal est conducteur. Nous choisissons donc le Fresnel adéquat :

$a^2 = 0.25 (\sqrt{(n_x^2 - k_x^2 - \sin^2 \theta)^2 + 4 n_x ^2} - k_x^2 - \sin^2 \theta$

$b^2 = 0.25 (\sqrt{(n_x^2 - k_x^2 - \sin^2 \theta)^2 + 4 n_x ^2} + k_x^2 + \sin^2 \theta$

 $R_s(n_x, k_x, \theta) = \frac {a^2 + b^2 - 2a\cos \theta + \cos^2 \theta} {a^2 + b^2 + 2a\cos \theta + \cos^2 \theta}$

$R_p(n_x, k_x, \theta) = R_s(n_x, k_x, \theta) \ast \frac {a^2+b^2 - 2a \sin \theta \tan \theta + \sin ^2 \theta \tan^2 \theta} {a^2+b^2 + 2a \sin \theta \tan \theta + \sin ^2 \theta \tan^2 \theta}$

avec $\theta$ l’angle du rayon par rapport à la micro-normale de la surface. Dans le cas du calcul du terme spéculaire,on préférera utiliser le demi-vecteur $h_x$

BDSF de type Matte

La BSDF Matte représente un objet complètement diffus, elle ne possède donc aucun terme spéculaire et ne se modélise qu’avec un modèle Oren-Nayar. Afin de représenter une surface complètement matte, le modèle ne doit posséder aucun terme spéculaire. On va donc annuler ce terme dans l’équation de BSDF et mettre un terme diffus tel que décrit par Oren-Nayar :

$s = (w_i \cdot w_o) - (n_x \cdot w_i) \ast (n_x \cdot w_o)$

$t = 1$ si $s \le 0$ sinon $t = max(n_x \cdot w_i, n_x \cdot w_o))$

$A= \frac{1}{\pi + \frac{\pi}{2} - \frac{2}{3}) \ast \alpha_x}$ $B= \frac{\sigma'}{\pi + \frac{\pi}{2} - \frac{2}{3}) \ast \alpha_x}$ 

$diff_{matte}(x, w_i, w_o) = c_{d,x}(n_x \cdot w_i) \ast (A+B \ast \frac{s}{t})$

BSDF de type substrate

Contrairement au modèle Plastique, le modèle Substrate permet aussi de représenter des interactions lumineuses sur plusieurs couches de matériaux. Pour ce type de BSDF nous allons nous servir d’un type de spéculaire de Cook-Torrance modifié. Afinde mieux représenter le facteur géométrique de ce type de surface, le terme $G$ de l’équation à été développé directement dans la formule de Cook-Torrance, ce qui nous donne :

$diff_{substrate}(x, w_i, w_o) = \frac{28 \ast c_{d,x}} {23 \ast \pi} \ast (1-k_s) \ast (1- \frac{1}{2} \ast |w_i \cdot n_x|) ^5$

$spec_{substrate}(x, w_i, w_o) = \frac {D(\alpha_x, w_h) \ast F(c_{s, x}, w_i \cdot w_h)} {4 \ast (w_i \cdot w_h) \ast \max (|w_i \cdot n_x|, |w_o \cdot n_x|)}$


Principled BSDF (Disney)

 Le Principled Shader de Disney fut créé sur l’envie de modéliser un grand nombre de matériaux physiquement réalistes à l’aide d’une seule fonction fortement paramétrisée. Pour faire cela, ils se sont donc basé sur les nombreux types de BSDF déjà disponibles et ont simplement composé ces dernières avec des paramètres définissant l’influence de chaque part sur le rendu final.La fonction ainsi obtenue prend la forme suivante :

$diff_{disney} = (1-\alpha_{met}) ((1-\alpha_{ss}) diffuse_{disney} + \alpha_{ss} \ast SS_{disney} + \alpha_{sheen} \ast sheen_{disney})$

$spec_{disney}(x, w_i, w_o) = \frac { D_{disney} \ast G \ast F_{shlick}} { 4 \cos(w_i \cdot n_x) \cos (w_o \cdot n_x)} + \alpha_{cc} \ast CC_{disney}$

Dans la suite nous utiliseront les symboles suivants :

$c_x$ La couleur de la surface en x
$\alpha_{specTint}$ Paramètre controlant la teinte du speculaire en x
$\alpha_{met}$ Paramètre controlant l'aspet métalique de la surface en x
$\alpha_{rou}$ La rugosité de la surface en x
$\alpha_{anisotropique}$ Paramètre controlant le taux d'anisotropie en x
$\alpha_{ss}$ Paramètre controlant le taux SS en x
$\alpha_{sheen}$ Paramètre controlant le taux de sheen en x
$\alpha_{sheenTint}$ Paramètre controlant la couleur du speculaire en x
$\alpha_{cc}$ Paramètre controlant le taux de ClearCoat en x
$\alpha_{ccGlow}$ Paramètre controlant la brillance du ClearCoat en x
$n_x$ L'indice de réfraction de la surface en x

On voit ainsi l’apparition de différentes fonctions qui vont venir définir les différentes parties du principled BSDF.

Terme diffus : Pour représenter la diffusion de ses matériaux, la BSDF utilise un modèle Lambertien auquelelle multiplie un Fresnel tel que présenté dans le cours SIGGRAPH de Burley en 2015

$diffuse_{disney} = \frac {c_x}{\pi} (1-\frac{1}{2} (1-(w_o \cdot n_x)^5)) (1-\frac{1}{2} (w_i \cdot n_x) ^5)$

Approximation du SubSurface Scattering : Afin d’obtenir un résultat toujours plus réaliste, le termediffus est mélangé avec une approximation d’un rendu modélisant la dispersion sous la surface de l’objet. Pource faire, l’équipe Disney s’est inspiré du modèle présenté par Hanrahan-Kruger

$F_{ss} = \text{lerp}((w_o \cdot n_x)^5, 1, (w_i \cdot h_x)^2 \alpha_x) \ast \text{lerp}((w_i \cdot n_x)^5, 1, (w_o \cdot h_x)^2 \alpha_x)$

$SS_{disney} = 1.25(F_{ss} \ast (\frac{1}{|w_o \cdot n_x| + |w_i \cdot n_x|} - \frac{1}{2}) + \frac{1}{2})$

Modélisation du Sheen : Le Sheen (ou brillance en français) sert à représenter l’effet de brillance qui peutêtre vu é angle rasant sur certains tissus. Cet ajout est donc modélisé comme une BSDF proche d’un Fresnel

$sheen_{disney} = \text{lerp}( \alpha_{sheenTint}, 1, kd_x) (1 - \cos(w_i \cdot h_x))^5$

Distribution des Normales : Le modèle de BSDF de Disney peut aussi représenter des modèles anisotro-piques. De ce fait un recalcul de D a été nécessaire afin de représenter ces surfaces

$aspect = \sqrt{1 - 0.9 \ast \alpha_{anisotropique}}$

$\alpha_{rou, x} = {\alpha_{rou}^2} \div {aspect}$ $\alpha_{rou, y} = {\alpha_{rou}^2} \cdot {aspect}$

$D_{disney} = \frac{1}{\pi} \frac{1}{\alpha_{rou, x} \alpha_{rou, y}} \frac{1}{ ((h_x \cdot \hat{x})^2 \div \alpha_{rou, x} + (h_x \cdot \hat{y})^2 \div \alpha_{rou, y} + (h_x \cdot n_x)^2)^2}$

Fresnel : Selon l’équipe Disney, l’approximation de Schlick est une méthode suffisamment précise pour obtenirles coefficients de Fresnel

$F_0 = (1-\alpha_{met}) \frac{\sqrt{n-1}}{\sqrt{n+1}} ((1-\alpha_{specTint}) + \alpha_{specTint}c_x) + \alpha_{met}c_x$

$F_{shlick}(F_0, \theta) = F_0 + (1-F_0)(1 - \cos \theta)^5$

Avec $F_0$ la réflexion spéculaire à incidence normale et $\theta$ l’angle entre le rayon et la micro-normale ($h_x$ dans le cas présent). Mais l’on remarquera que le modèle peut présenter des erreurs dans le cas de matériaux métalliques.

Addition d’un effet ClearCoat : L’effet ClearCoat définit par l’équipe Disney permet de modéliser unecouche de vernis par dessus la surface de l’objet. Cet effet est obtenu en ajoutant un lobe spéculaire a la BSDF avec une rugosité fixée $\alpha = 0.25$ et un IOR de $1.5$

$CC_{disney} = \frac{1}{4} G(w_i, w_o, h_x, n_x, \frac{1}{4}) \ast F_{shlick}(0.04, w_i \cdot h_x) \ast D_{GTR1}(h_x, n_x, \text{lerp}(\alpha_{ccgloss}, 0.1, 0.001)$

Contrairement aux autres BSDFs, nous choisissons une distribution GTR 1.0 car celle-ci permet des traînées proche de la réalité dans le cas d’une couche translucide :

$D_{GTR1}(h_x, n_x, \alpha) = \frac{\alpha^2-1} { \pi \log(\alpha^2)(1+(\alpha^2-1)(h_x \cdot n_x)^2}$


Example Substrate
vec3 computeMaterialInternal(Material material, vec2 texC, vec3 L, vec3 V, vec3 N, vec3 X, vec3 Y) {
    vec3 W =            normalize(L + V);

    vec2 alpha =        getRoughness(material, texC);
    vec3 Kd =			getKd(material, texC);
    vec3 Ks =			getKs(material, texC);

    float WdN =         abs(dot(W, N));
    float WdX =         dot(W, X);
    float WdY =         dot(W, Y);
    float WdV =         max(dot(W, V), 0.f);
    float WdL =         max(dot(W, L), 0.f);
    float NdV =         abs(dot(V, N)) + 0.005f;
    float NdL =         max(dot(L, N), 0.f);

    float D =   		GTR2_AnisotropicDistribution(WdN, WdX, WdY, alpha.x, alpha.y);

    vec3 diff = 		(28.f / (23.f * Pi)) * Kd * (vec3(1.f) - Ks) *
    						(1.f - pow5(1.f - 0.5f*NdV)) * (1.f - pow5(1.f - 0.5f*NdL));
    
    vec3 spec = 		D / (4 * abs(WdL) * max(NdV, NdL)) * FresnelShlick(WdL, Ks, vec3(1.0f));

    return (diff + spec) * NdL;
}
Example Plastique
vec3 computeMaterialInternal(Material material, vec2 texC, vec3 L, vec3 V, vec3 N, vec3 X, vec3 Y) {
    vec3 W =            normalize(L + V);

    float alpha =       getRoughness(material, texC).r;

    float WdN =         abs(dot(W, N));
    float WdX =         dot(W, X);
    float WdY =         dot(W, Y);
    float WdV =         max(dot(W, V), 0.f);
    float WdL =         max(dot(W, L), 0.f);
    float NdV =         abs(dot(V, N)) + 0.005f;
    float NdL =         max(dot(L, N), 0.f);

    float D =   GTR2_AnisotropicDistribution(WdN, WdX, WdY, alpha, alpha);
    float G =   GGX_GeometryTerm(WdV, NdV, alpha) * GGX_GeometryTerm(WdL, NdL, alpha);
    float F =   FresnelDielectric(WdL, 1.5f, 1.0f);

    vec3 diff = getKd(material, texC) / Pi;
    vec3 spec = getKs(material, texC) * D * G * F  / (4 * NdL * NdV);

    return (spec + diff) * NdL;
}
Example OrenNayar
float A(float sigma2) {
	return 1.0 - 0.5 * (sigma2/(sigma2 + 0.33));
}

float B(float sigma2) {
	return 0.45 * (sigma2/(sigma2 + 0.09));
}


float geomAtt(vec3 l, vec3 n, vec3 e) {
	float NdotE = dot(n, e);
	float NdotL = dot(n, l);
	
	float maxnlne = max(NdotL, NdotE);
	float minnlne = min(NdotL, NdotE);
	float factor = sqrt( (1 - maxnlne*maxnlne) * (1 - minnlne*minnlne) );
	factor = factor / maxnlne;
	
	vec3 emn = normalize(e - n*NdotE);
	vec3 lmn = normalize(l - n*NdotL);
	
	return max( 0, dot(emn,lmn) )*factor;
}


vec3 computeMaterialInternal(Material material, vec2 texC, vec3 L, vec3 V, vec3 N, vec3 X, vec3 Y) {
	float costi = max(dot(N, L), 0);
	float sigma2 = getSigma2(material, texC);

	// is there an "intensity" parameter for lights ?
	return (getKd(material, texC) * InvPi) * costi * ( A(sigma2) + B(sigma2) * geomAtt(L, N, V));
}
Example Metal
vec3 computeMaterialInternal(Material material, vec2 texC, vec3 L, vec3 V, vec3 N, vec3 X, vec3 Y) {
    vec3 W =            normalize(L + V);

    vec2 alpha =        getRoughness(material, texC);

    float WdN =         abs(dot(W, N));
    float WdX =         dot(W, X);
    float WdY =         dot(W, Y);
    float WdV =         max(dot(W, V), 0.f);
    float WdL =         max(dot(W, L), 0.f);
    float NdV =         abs(dot(V, N)) + 0.005f;
    float NdL =         max(dot(L, N), 0.f);

    float D =   GTR2_AnisotropicDistribution(WdN, WdX, WdY, alpha.x, alpha.y);
    float G =   GGX_GeometryTerm(WdV, NdV, alpha.x) * GGX_GeometryTerm(WdL, NdL, alpha.x);
    vec3 F =    FresnelConductor(WdL, vec3(1.0f), getEta(material, texC), getK(material, texC));

    vec3 spec = D * G * F  / (4 * NdL * NdV);

    return spec * NdL;
}
Example Disney
vec3 computeMaterialInternal(Material material, vec2 texC, vec3 L, vec3 V, vec3 N, vec3 X, vec3 Y) {
    vec3 W =            normalize(L + V);

    float WdN =         max(dot(W, N), 0.f);
    float WdX =         dot(W, X);
    float WdY =         dot(W, Y);
    float WdV =         max(dot(W, V), 0.f);
    float WdL =         max(dot(W, L), 0.f);
    float NdV =         abs(dot(V, N)) + 0.005f;
    float NdL =         max(dot(L, N), 0.f);

    vec3 color = 		getColor(material, texC);
    int thin = 			material.thin;

    float metallic =	getMetallic(material, texC);
    float flatness = 	getFlatness(material, texC);
    float roughness =	getRoughness(material, texC).x;
    float sheen =       getSheen(material, texC);
    float sheenTint =   getSheenTint(material, texC);
    float ior =         getIor(material, texC);
    float clearcoat =   getClearcoat(material, texC);
    float gloss =       getClearcoatGloss(material, texC);
    float specTrans = 	getSpecTrans(material, texC);
    float diffTrans =	getDiffTrans(material, texC) / 2.f;
    float specTint =    getSpecularTint(material, texC);

    //Diffuse Term-----------------------------------------------------------------------------------
    
    vec3 diff = vec3(disneyDiffuse(NdV, NdL, WdL, roughness));

    float diffWeight =  (1.f - metallic) * (1.f - specTrans);
    
    if (thin == 1) {
        diff *= diffWeight * (1.f - diffTrans) * (1.f - flatness) * color;
    	diff += disneyFakeSS(diffWeight * (1.f - diffTrans) * flatness * color, roughness, WdL, NdV, NdL);
    }
    else {
    	diff *= diffWeight * color;
    }

    diff += disneyRetro(diffWeight * color, roughness, WdL, NdV, NdL);

    //Calculatio of the diffuse sheen
    diff += diffWeight * disneySheen(color, sheen, sheenTint, WdL);

    //Specular Term---------------------------------------------------------------------------------
    
    vec3 spec;

    float aspect =      sqrt( 1.f - getAnisotropic(material, texC) * 0.9f);
    float ax =          max(0.001f, sqr(roughness) / aspect);
    float ay =          max(0.001f, sqr(roughness) * aspect);
    float nIor =        sqr(ior - 1) / sqr(ior + 1);
    vec3 Cspec0 =       mix(nIor * mix(vec3(1.f), color, specTint), color, metallic);

    float D =   GTR2_AnisotropicDistribution(WdN, WdX, WdY, ax, ay);
    float G =   GGX_GeometryTerm(WdV, NdV, roughness) * GGX_GeometryTerm(WdL, NdL, roughness);
    vec3 F =    mix(vec3(FresnelDielectric(WdL, 1.f, ior)),
                    FresnelShlick(WdL, Cspec0, vec3(1.0f)),
                    metallic
                    );

    spec = D * G * F / (4 * dot(N, L) * dot(N, V));
    spec += disneyClearCoat(clearcoat, gloss, WdN, WdL, WdV, NdV, NdL);

    return (diff + spec) * NdL;
}
code commun

float GTR2_AnisotropicDistribution(float NdH, float HdX, float HdY, float ax, float ay) {
    return 1.0f / (Pi * ax*ay * sqr(sqr(HdX/ax) + sqr(HdY/ay) + NdH*NdH));
}

float chiGGX(float v) {
    return v > 0.0f ? 1.0f : 0.0f;
}

float GGX_DistributionTerm( float NdotH , float m ) {
    float m2 = m * m ;
    float NdotH2 = NdotH * NdotH;
    float f = NdotH2 * m2 + (1 - NdotH2);
    return (chiGGX(NdotH) * m2) / (f * f * Pi);
}

float GGX_GeometryTerm(float VdH, float VdN, float alpha) {
    float chi = chiGGX(VdH / VdN);
    VdH = VdH * VdH;

    float tan2 = (1-VdH)/VdH;
    return (chi * 2) / (1 + sqrt(1+alpha*alpha*tan2));
}

vec3 FresnelShlick(float cosThetaI, vec3 F0, vec3 F90) {
    return F0 + pow5(1 - cosThetaI) * (F90 - F0);
}


float FresnelDielectric(float cosThetaI, float etaI, float etaT) {
    float sinThetaI = sqrt(max(0.f, 1 - cosThetaI * cosThetaI));
    float sinThetaT = etaI / etaT * sinThetaI;

    float cosThetaT = sqrt(max(0.f, 1.f - sinThetaT * sinThetaT));
    float Rparl = ((etaT * cosThetaI) - (etaI * cosThetaT)) /
                  ((etaT * cosThetaI) + (etaI * cosThetaT));
    float Rperp = ((etaI * cosThetaI) - (etaT * cosThetaT)) /
                  ((etaI * cosThetaI) + (etaT * cosThetaT));

    return (Rparl * Rparl + Rperp * Rperp) / 2;
}

vec3 FresnelConductor(float cosThetaI, vec3 etaI, vec3 etaT, vec3 k) {
    vec3 eta = etaT / etaI;
    vec3 etak = k / etaI;

    float cosThetaI2 = cosThetaI * cosThetaI;
    float sinThetaI2 = 1. - cosThetaI2;
    vec3 eta2 = eta * eta;
    vec3 etak2 = etak * etak;

    vec3 t0 = eta2 - etak2 - sinThetaI2;
    vec3 a2plusb2 = sqrt(t0 * t0 + 4 * eta2 * etak2);
    vec3 t1 = a2plusb2 + cosThetaI2;
    vec3 a = sqrt(0.5f * (a2plusb2 + t0));
    vec3 t2 = 2.f * cosThetaI * a;
    vec3 Rs = (t1 - t2) / (t1 + t2);

    vec3 t3 = cosThetaI2 * a2plusb2 + sinThetaI2 * sinThetaI2;
    vec3 t4 = t2 * sinThetaI2;
    vec3 Rp = Rs * (t3 - t4) / (t3 + t4);

    return 0.5f * (Rp + Rs);
}

A vous de faire le reste ;)