0031096: Visualization, TKOpenGl - support metallic-roughness texture mapping
[occt.git] / src / Shaders / PathtraceBase.fs
1 #ifdef _MSC_VER
2   #define PATH_TRACING // just for editing in MS VS
3
4   #define in
5   #define out
6   #define inout
7
8   typedef struct { float x; float y; } vec2;
9   typedef struct { float x; float y; float z; } vec3;
10   typedef struct { float x; float y; float z; float w; } vec4;
11 #endif
12
13 #ifdef PATH_TRACING
14
15 ///////////////////////////////////////////////////////////////////////////////////////
16 // Specific data types
17
18 //! Describes local space at the hit point (visualization space).
19 struct SLocalSpace
20 {
21   //! Local X axis.
22   vec3 AxisX;
23
24   //! Local Y axis.
25   vec3 AxisY;
26
27   //! Local Z axis.
28   vec3 AxisZ;
29 };
30
31 //! Describes material properties (BSDF).
32 struct SBSDF
33 {
34   //! Weight of coat specular/glossy BRDF.
35   vec4 Kc;
36
37   //! Weight of base diffuse BRDF + base color texture index in W.
38   vec4 Kd;
39
40   //! Weight of base specular/glossy BRDF.
41   vec4 Ks;
42
43   //! Weight of base specular/glossy BTDF + metallic-roughness texture index in W.
44   vec4 Kt;
45
46   //! Fresnel coefficients of coat layer.
47   vec3 FresnelCoat;
48
49   //! Fresnel coefficients of base layer.
50   vec3 FresnelBase;
51 };
52
53 ///////////////////////////////////////////////////////////////////////////////////////
54 // Support subroutines
55
56 //=======================================================================
57 // function : buildLocalSpace
58 // purpose  : Generates local space for the given normal
59 //=======================================================================
60 SLocalSpace buildLocalSpace (in vec3 theNormal)
61 {
62   vec3 anAxisX = vec3 (theNormal.z, 0.f, -theNormal.x);
63   vec3 anAxisY = vec3 (0.f, -theNormal.z, theNormal.y);
64
65   float aSqrLenX = dot (anAxisX, anAxisX);
66   float aSqrLenY = dot (anAxisY, anAxisY);
67
68   if (aSqrLenX > aSqrLenY)
69   {
70     anAxisX *= inversesqrt (aSqrLenX);
71     anAxisY = cross (anAxisX, theNormal);
72   }
73   else
74   {
75     anAxisY *= inversesqrt (aSqrLenY);
76     anAxisX = cross (anAxisY, theNormal);
77   }
78
79   return SLocalSpace (anAxisX, anAxisY, theNormal);
80 }
81
82 //=======================================================================
83 // function : toLocalSpace
84 // purpose  : Transforms the vector to local space from world space
85 //=======================================================================
86 vec3 toLocalSpace (in vec3 theVector, in SLocalSpace theSpace)
87 {
88   return vec3 (dot (theVector, theSpace.AxisX),
89                dot (theVector, theSpace.AxisY),
90                dot (theVector, theSpace.AxisZ));
91 }
92
93 //=======================================================================
94 // function : fromLocalSpace
95 // purpose  : Transforms the vector from local space to world space
96 //=======================================================================
97 vec3 fromLocalSpace (in vec3 theVector, in SLocalSpace theSpace)
98 {
99   return theVector.x * theSpace.AxisX +
100          theVector.y * theSpace.AxisY +
101          theVector.z * theSpace.AxisZ;
102 }
103
104 //=======================================================================
105 // function : convolve
106 // purpose  : Performs a linear convolution of the vector components
107 //=======================================================================
108 float convolve (in vec3 theVector, in vec3 theFactor)
109 {
110   return dot (theVector, theFactor) * (1.f / max (theFactor.x + theFactor.y + theFactor.z, 1e-15f));
111 }
112
113 //=======================================================================
114 // function : fresnelSchlick
115 // purpose  : Computes the Fresnel reflection formula using
116 //            Schlick's approximation.
117 //=======================================================================
118 vec3 fresnelSchlick (in float theCosI, in vec3 theSpecularColor)
119 {
120   return theSpecularColor + (UNIT - theSpecularColor) * pow (1.f - theCosI, 5.f);
121 }
122
123 //=======================================================================
124 // function : fresnelDielectric
125 // purpose  : Computes the Fresnel reflection formula for dielectric in
126 //            case of circularly polarized light (Based on PBRT code).
127 //=======================================================================
128 float fresnelDielectric (in float theCosI,
129                          in float theCosT,
130                          in float theEtaI,
131                          in float theEtaT)
132 {
133   float aParl = (theEtaT * theCosI - theEtaI * theCosT) /
134                 (theEtaT * theCosI + theEtaI * theCosT);
135
136   float aPerp = (theEtaI * theCosI - theEtaT * theCosT) /
137                 (theEtaI * theCosI + theEtaT * theCosT);
138
139   return (aParl * aParl + aPerp * aPerp) * 0.5f;
140 }
141
142 #define ENVIRONMENT_IOR 1.f
143
144 //=======================================================================
145 // function : fresnelDielectric
146 // purpose  : Computes the Fresnel reflection formula for dielectric in
147 //            case of circularly polarized light (based on PBRT code)
148 //=======================================================================
149 float fresnelDielectric (in float theCosI, in float theIndex)
150 {
151   float aFresnel = 1.f;
152
153   float anEtaI = theCosI > 0.f ? 1.f : theIndex;
154   float anEtaT = theCosI > 0.f ? theIndex : 1.f;
155
156   float aSinT2 = (anEtaI * anEtaI) / (anEtaT * anEtaT) * (1.f - theCosI * theCosI);
157
158   if (aSinT2 < 1.f)
159   {
160     aFresnel = fresnelDielectric (abs (theCosI), sqrt (1.f - aSinT2), anEtaI, anEtaT);
161   }
162
163   return aFresnel;
164 }
165
166 //=======================================================================
167 // function : fresnelConductor
168 // purpose  : Computes the Fresnel reflection formula for conductor in case
169 //            of circularly polarized light (based on PBRT source code)
170 //=======================================================================
171 float fresnelConductor (in float theCosI, in float theEta, in float theK)
172 {
173   float aTmp = 2.f * theEta * theCosI;
174
175   float aTmp1 = theEta * theEta + theK * theK;
176
177   float aSPerp = (aTmp1 - aTmp + theCosI * theCosI) /
178                  (aTmp1 + aTmp + theCosI * theCosI);
179
180   float aTmp2 = aTmp1 * theCosI * theCosI;
181
182   float aSParl = (aTmp2 - aTmp + 1.f) /
183                  (aTmp2 + aTmp + 1.f);
184
185   return (aSPerp + aSParl) * 0.5f;
186 }
187
188 #define FRESNEL_SCHLICK    -0.5f
189 #define FRESNEL_CONSTANT   -1.5f
190 #define FRESNEL_CONDUCTOR  -2.5f
191 #define FRESNEL_DIELECTRIC -3.5f
192
193 //=======================================================================
194 // function : fresnelMedia
195 // purpose  : Computes the Fresnel reflection formula for general medium
196 //            in case of circularly polarized light.
197 //=======================================================================
198 vec3 fresnelMedia (in float theCosI, in vec3 theFresnel)
199 {
200   vec3 aFresnel;
201
202   if (theFresnel.x > FRESNEL_SCHLICK)
203   {
204     aFresnel = fresnelSchlick (abs (theCosI), theFresnel);
205   }
206   else if (theFresnel.x > FRESNEL_CONSTANT)
207   {
208     aFresnel = vec3 (theFresnel.z);
209   }
210   else if (theFresnel.x > FRESNEL_CONDUCTOR)
211   {
212     aFresnel = vec3 (fresnelConductor (abs (theCosI), theFresnel.y, theFresnel.z));
213   }
214   else
215   {
216     aFresnel = vec3 (fresnelDielectric (theCosI, theFresnel.y));
217   }
218
219   return aFresnel;
220 }
221
222 //=======================================================================
223 // function : transmitted
224 // purpose  : Computes transmitted direction in tangent space
225 //            (in case of TIR returned result is undefined!)
226 //=======================================================================
227 void transmitted (in float theIndex, in vec3 theIncident, out vec3 theTransmit)
228 {
229   // Compute relative index of refraction
230   float anEta = (theIncident.z > 0.f) ? 1.f / theIndex : theIndex;
231
232   // Handle total internal reflection (TIR)
233   float aSinT2 = anEta * anEta * (1.f - theIncident.z * theIncident.z);
234
235   // Compute direction of transmitted ray
236   float aCosT = sqrt (1.f - min (aSinT2, 1.f)) * sign (-theIncident.z);
237
238   theTransmit = normalize (vec3 (-anEta * theIncident.x,
239                                  -anEta * theIncident.y,
240                                   aCosT));
241 }
242
243 //////////////////////////////////////////////////////////////////////////////////////////////
244 // Handlers and samplers for materials
245 //////////////////////////////////////////////////////////////////////////////////////////////
246
247 //=======================================================================
248 // function : EvalLambertianReflection
249 // purpose  : Evaluates Lambertian BRDF, with cos(N, PSI)
250 //=======================================================================
251 float EvalLambertianReflection (in vec3 theWi, in vec3 theWo)
252 {
253   return (theWi.z <= 0.f || theWo.z <= 0.f) ? 0.f : theWi.z * (1.f / M_PI);
254 }
255
256 #define FLT_EPSILON 1.0e-5f
257
258 //=======================================================================
259 // function : SmithG1
260 // purpose  :
261 //=======================================================================
262 float SmithG1 (in vec3 theDirection, in vec3 theM, in float theRoughness)
263 {
264   float aResult = 0.f;
265
266   if (dot (theDirection, theM) * theDirection.z > 0.f)
267   {
268     float aTanThetaM = sqrt (1.f - theDirection.z * theDirection.z) / theDirection.z;
269
270     if (aTanThetaM == 0.f)
271     {
272       aResult = 1.f;
273     }
274     else
275     {
276       float aVal = 1.f / (theRoughness * aTanThetaM);
277
278       // Use rational approximation to shadowing-masking function (from Mitsuba)
279       aResult = (3.535f + 2.181f * aVal) / (1.f / aVal + 2.276f + 2.577f * aVal);
280     }
281   }
282
283   return min (aResult, 1.f);
284 }
285
286 //=======================================================================
287 // function : EvalBlinnReflection
288 // purpose  : Evaluates Blinn glossy BRDF, with cos(N, PSI)
289 //=======================================================================
290 vec3 EvalBlinnReflection (in vec3 theWi, in vec3 theWo, in vec3 theFresnel, in float theRoughness)
291 {
292   // calculate the reflection half-vec
293   vec3 aH = normalize (theWi + theWo);
294
295   // roughness value -> Blinn exponent
296   float aPower = max (2.f / (theRoughness * theRoughness) - 2.f, 0.f);
297
298   // calculate microfacet distribution
299   float aD = (aPower + 2.f) * (1.f / M_2_PI) * pow (aH.z, aPower);
300
301   // calculate shadow-masking function
302   float aG = SmithG1 (theWo, aH, theRoughness) *
303              SmithG1 (theWi, aH, theRoughness);
304
305   // return total amount of reflection
306   return (theWi.z <= 0.f || theWo.z <= 0.f) ? ZERO :
307     aD * aG / (4.f * theWo.z) * fresnelMedia (dot (theWo, aH), theFresnel);
308 }
309
310 //=======================================================================
311 // function : EvalBsdfLayered
312 // purpose  : Evaluates BSDF for specified material, with cos(N, PSI)
313 //=======================================================================
314 vec3 EvalBsdfLayered (in SBSDF theBSDF, in vec3 theWi, in vec3 theWo)
315 {
316 #ifdef TWO_SIDED_BXDF
317   theWi.z *= sign (theWi.z);
318   theWo.z *= sign (theWo.z);
319 #endif
320
321   vec3 aBxDF = theBSDF.Kd.rgb * EvalLambertianReflection (theWi, theWo);
322
323   if (theBSDF.Ks.w > FLT_EPSILON)
324   {
325     aBxDF += theBSDF.Ks.rgb * EvalBlinnReflection (theWi, theWo, theBSDF.FresnelBase, theBSDF.Ks.w);
326   }
327
328   aBxDF *= UNIT - fresnelMedia (theWo.z, theBSDF.FresnelCoat);
329
330   if (theBSDF.Kc.w > FLT_EPSILON)
331   {
332     aBxDF += theBSDF.Kc.rgb * EvalBlinnReflection (theWi, theWo, theBSDF.FresnelCoat, theBSDF.Kc.w);
333   }
334
335   return aBxDF;
336 }
337
338 //=======================================================================
339 // function : SampleLambertianReflection
340 // purpose  : Samples Lambertian BRDF, W = BRDF * cos(N, PSI) / PDF(PSI)
341 //=======================================================================
342 vec3 SampleLambertianReflection (in vec3 theWo, out vec3 theWi, inout float thePDF)
343 {
344   float aKsi1 = RandFloat();
345   float aKsi2 = RandFloat();
346
347   theWi = vec3 (cos (M_2_PI * aKsi1),
348                 sin (M_2_PI * aKsi1),
349                 sqrt (1.f - aKsi2));
350
351   theWi.xy *= sqrt (aKsi2);
352
353 #ifdef TWO_SIDED_BXDF
354   theWi.z *= sign (theWo.z);
355 #endif
356
357   thePDF *= theWi.z * (1.f / M_PI);
358
359 #ifdef TWO_SIDED_BXDF
360   return UNIT;
361 #else
362   return UNIT * step (0.f, theWo.z);
363 #endif
364 }
365
366 //=======================================================================
367 // function : SampleGlossyBlinnReflection
368 // purpose  : Samples Blinn BRDF, W = BRDF * cos(N, PSI) / PDF(PSI)
369 //            The BRDF is a product of three main terms, D, G, and F,
370 //            which is then divided by two cosine terms. Here we perform
371 //            importance sample the D part of the Blinn model; trying to
372 //            develop a sampling procedure that accounted for all of the
373 //            terms would be complex, and it is the D term that accounts
374 //            for most of the variation.
375 //=======================================================================
376 vec3 SampleGlossyBlinnReflection (in vec3 theWo, out vec3 theWi, in vec3 theFresnel, in float theRoughness, inout float thePDF)
377 {
378   float aKsi1 = RandFloat();
379   float aKsi2 = RandFloat();
380
381   // roughness value --> Blinn exponent
382   float aPower = max (2.f / (theRoughness * theRoughness) - 2.f, 0.f);
383
384   // normal from microface distribution
385   float aCosThetaM = pow (aKsi1, 1.f / (aPower + 2.f));
386
387   vec3 aM = vec3 (cos (M_2_PI * aKsi2),
388                   sin (M_2_PI * aKsi2),
389                   aCosThetaM);
390
391   aM.xy *= sqrt (1.f - aCosThetaM * aCosThetaM);
392
393   // calculate PDF of sampled direction
394   thePDF *= (aPower + 2.f) * (1.f / M_2_PI) * pow (aCosThetaM, aPower + 1.f);
395
396 #ifdef TWO_SIDED_BXDF
397   bool toFlip = theWo.z < 0.f;
398
399   if (toFlip)
400     theWo.z = -theWo.z;
401 #endif
402
403   float aCosDelta = dot (theWo, aM);
404
405   // pick input based on half direction
406   theWi = -theWo + 2.f * aCosDelta * aM;
407
408   if (theWi.z <= 0.f || theWo.z <= 0.f)
409   {
410     return ZERO;
411   }
412
413   // Jacobian of half-direction mapping
414   thePDF /= 4.f * aCosDelta;
415
416   // compute shadow-masking coefficient
417   float aG = SmithG1 (theWo, aM, theRoughness) *
418              SmithG1 (theWi, aM, theRoughness);
419
420 #ifdef TWO_SIDED_BXDF
421   if (toFlip)
422     theWi.z = -theWi.z;
423 #endif
424
425   return (aG * aCosDelta) / (theWo.z * aM.z) * fresnelMedia (aCosDelta, theFresnel);
426 }
427
428 //=======================================================================
429 // function : BsdfPdfLayered
430 // purpose  : Calculates BSDF of sampling input knowing output
431 //=======================================================================
432 float BsdfPdfLayered (in SBSDF theBSDF, in vec3 theWo, in vec3 theWi, in vec3 theWeight)
433 {
434   float aPDF = 0.f; // PDF of sampling input direction
435
436   // We choose whether the light is reflected or transmitted
437   // by the coating layer according to the Fresnel equations
438   vec3 aCoatF = fresnelMedia (theWo.z, theBSDF.FresnelCoat);
439
440   // Coat BRDF is scaled by its Fresnel reflectance term. For
441   // reasons of simplicity we scale base BxDFs only by coat's
442   // Fresnel transmittance term
443   vec3 aCoatT = UNIT - aCoatF;
444
445   float aPc = dot (theBSDF.Kc.rgb * aCoatF, theWeight);
446   float aPd = dot (theBSDF.Kd.rgb * aCoatT, theWeight);
447   float aPs = dot (theBSDF.Ks.rgb * aCoatT, theWeight);
448   float aPt = dot (theBSDF.Kt.rgb * aCoatT, theWeight);
449
450   if (theWi.z * theWo.z > 0.f)
451   {
452     vec3 aH = normalize (theWi + theWo);
453
454     aPDF = aPd * abs (theWi.z / M_PI);
455
456     if (theBSDF.Kc.w > FLT_EPSILON)
457     {
458       float aPower = max (2.f / (theBSDF.Kc.w * theBSDF.Kc.w) - 2.f, 0.f); // roughness --> exponent
459
460       aPDF += aPc * (aPower + 2.f) * (0.25f / M_2_PI) * pow (abs (aH.z), aPower + 1.f) / dot (theWi, aH);
461     }
462
463     if (theBSDF.Ks.w > FLT_EPSILON)
464     {
465       float aPower = max (2.f / (theBSDF.Ks.w * theBSDF.Ks.w) - 2.f, 0.f); // roughness --> exponent
466
467       aPDF += aPs * (aPower + 2.f) * (0.25f / M_2_PI) * pow (abs (aH.z), aPower + 1.f) / dot (theWi, aH);
468     }
469   }
470
471   return aPDF / (aPc + aPd + aPs + aPt);
472 }
473
474 //! Tool macro to handle sampling of particular BxDF
475 #define PICK_BXDF_LAYER(p, k) aPDF = p / aTotalR; theWeight *= k / aPDF;
476
477 //=======================================================================
478 // function : SampleBsdfLayered
479 // purpose  : Samples specified composite material (BSDF)
480 //=======================================================================
481 float SampleBsdfLayered (in SBSDF theBSDF, in vec3 theWo, out vec3 theWi, inout vec3 theWeight, inout bool theInside)
482 {
483   // NOTE: OCCT uses two-layer material model. We have base diffuse, glossy, or transmissive
484   // layer, covered by one glossy/specular coat. In the current model, the layers themselves
485   // have no thickness; they can simply reflect light or transmits it to the layer under it.
486   // We use actual BRDF model only for direct reflection by the coat layer. For transmission
487   // through this layer, we approximate it as a flat specular surface.
488
489   float aPDF = 0.f; // PDF of sampled direction
490
491   // We choose whether the light is reflected or transmitted
492   // by the coating layer according to the Fresnel equations
493   vec3 aCoatF = fresnelMedia (theWo.z, theBSDF.FresnelCoat);
494
495   // Coat BRDF is scaled by its Fresnel term. According to
496   // Wilkie-Weidlich layered BSDF model, transmission term
497   // for light passing through the coat at direction I and
498   // leaving it in O is T = ( 1 - F (O) ) x ( 1 - F (I) ).
499   // For reasons of simplicity, we discard the second term
500   // and scale base BxDFs only by the first term.
501   vec3 aCoatT = UNIT - aCoatF;
502
503   float aPc = dot (theBSDF.Kc.rgb * aCoatF, theWeight);
504   float aPd = dot (theBSDF.Kd.rgb * aCoatT, theWeight);
505   float aPs = dot (theBSDF.Ks.rgb * aCoatT, theWeight);
506   float aPt = dot (theBSDF.Kt.rgb * aCoatT, theWeight);
507
508   // Calculate total reflection probability
509   float aTotalR = (aPc + aPd) + (aPs + aPt);
510
511   // Generate random variable to select BxDF
512   float aKsi = aTotalR * RandFloat();
513
514   if (aKsi < aPc) // REFLECTION FROM COAT
515   {
516     PICK_BXDF_LAYER (aPc, theBSDF.Kc.rgb)
517
518     if (theBSDF.Kc.w < FLT_EPSILON)
519     {
520       theWeight *= aCoatF;
521
522       theWi = vec3 (-theWo.x,
523                     -theWo.y,
524                      theWo.z);
525     }
526     else
527     {
528       theWeight *= SampleGlossyBlinnReflection (theWo, theWi, theBSDF.FresnelCoat, theBSDF.Kc.w, aPDF);
529     }
530
531     aPDF = mix (aPDF, MAXFLOAT, theBSDF.Kc.w < FLT_EPSILON);
532   }
533   else if (aKsi < aTotalR) // REFLECTION FROM BASE
534   {
535     theWeight *= aCoatT;
536
537     if (aKsi < aPc + aPd) // diffuse BRDF
538     {
539       PICK_BXDF_LAYER (aPd, theBSDF.Kd.rgb)
540
541       theWeight *= SampleLambertianReflection (theWo, theWi, aPDF);
542     }
543     else if (aKsi < (aPc + aPd) + aPs) // specular/glossy BRDF
544     {
545       PICK_BXDF_LAYER (aPs, theBSDF.Ks.rgb)
546
547       if (theBSDF.Ks.w < FLT_EPSILON)
548       {
549         theWeight *= fresnelMedia (theWo.z, theBSDF.FresnelBase);
550
551         theWi = vec3 (-theWo.x,
552                       -theWo.y,
553                        theWo.z);
554       }
555       else
556       {
557         theWeight *= SampleGlossyBlinnReflection (theWo, theWi, theBSDF.FresnelBase, theBSDF.Ks.w, aPDF);
558       }
559
560       aPDF = mix (aPDF, MAXFLOAT, theBSDF.Ks.w < FLT_EPSILON);
561     }
562     else // specular transmission
563     {
564       PICK_BXDF_LAYER (aPt, theBSDF.Kt.rgb)
565
566       // refracted direction should exist if we are here
567       transmitted (theBSDF.FresnelCoat.y, theWo, theWi);
568
569       theInside = !theInside; aPDF = MAXFLOAT;
570     }
571   }
572
573   // path termination for extra small weights
574   theWeight = mix (ZERO, theWeight, step (FLT_EPSILON, aTotalR));
575
576   return aPDF;
577 }
578
579 //////////////////////////////////////////////////////////////////////////////////////////////
580 // Handlers and samplers for light sources
581 //////////////////////////////////////////////////////////////////////////////////////////////
582
583 // =======================================================================
584 // function : Latlong
585 // purpose  : Converts world direction to environment texture coordinates
586 // =======================================================================
587 vec2 Latlong (in vec3 thePoint)
588 {
589   float aPsi = acos (-thePoint.z);
590
591   float aPhi = atan (thePoint.y, thePoint.x) + M_PI;
592
593   return vec2 (aPhi * 0.1591549f,
594                aPsi * 0.3183098f);
595 }
596
597 //=======================================================================
598 // function : SampleLight
599 // purpose  : General sampling function for directional and point lights
600 //=======================================================================
601 vec3 SampleLight (in vec3 theToLight, inout float theDistance, in bool isInfinite, in float theSmoothness, inout float thePDF)
602 {
603   SLocalSpace aSpace = buildLocalSpace (theToLight * (1.f / theDistance));
604
605   // for point lights smoothness defines radius
606   float aCosMax = isInfinite ? theSmoothness :
607     inversesqrt (1.f + theSmoothness * theSmoothness / (theDistance * theDistance));
608
609   float aKsi1 = RandFloat();
610   float aKsi2 = RandFloat();
611
612   float aTmp = 1.f - aKsi2 * (1.f - aCosMax);
613
614   vec3 anInput = vec3 (cos (M_2_PI * aKsi1),
615                        sin (M_2_PI * aKsi1),
616                        aTmp);
617
618   anInput.xy *= sqrt (1.f - aTmp * aTmp);
619
620   thePDF = (aCosMax < 1.f) ? (thePDF / M_2_PI) / (1.f - aCosMax) : MAXFLOAT;
621
622   return normalize (fromLocalSpace (anInput, aSpace));
623 }
624
625 //=======================================================================
626 // function : HandlePointLight
627 // purpose  :
628 //=======================================================================
629 float HandlePointLight (in vec3 theInput, in vec3 theToLight, in float theRadius, in float theDistance, inout float thePDF)
630 {
631   float aCosMax = inversesqrt (1.f + theRadius * theRadius / (theDistance * theDistance));
632
633   float aVisibility = step (aCosMax, dot (theInput, theToLight));
634
635   thePDF *= step (-1.f, -aCosMax) * aVisibility * (1.f / M_2_PI) / (1.f - aCosMax);
636
637   return aVisibility;
638 }
639
640 //=======================================================================
641 // function : HandleDistantLight
642 // purpose  :
643 //=======================================================================
644 float HandleDistantLight (in vec3 theInput, in vec3 theToLight, in float theCosMax, inout float thePDF)
645 {
646   float aVisibility = step (theCosMax, dot (theInput, theToLight));
647
648   thePDF *= step (-1.f, -theCosMax) * aVisibility * (1.f / M_2_PI) / (1.f - theCosMax);
649
650   return aVisibility;
651 }
652
653 // =======================================================================
654 // function: IntersectLight
655 // purpose : Checks intersections with light sources
656 // =======================================================================
657 vec3 IntersectLight (in SRay theRay, in int theDepth, in float theHitDistance, out float thePDF)
658 {
659   vec3 aTotalRadiance = ZERO;
660
661   thePDF = 0.f; // PDF of sampling light sources
662
663   for (int aLightIdx = 0; aLightIdx < uLightCount; ++aLightIdx)
664   {
665     vec4 aLight = texelFetch (
666       uRaytraceLightSrcTexture, LIGHT_POS (aLightIdx));
667     vec4 aParam = texelFetch (
668       uRaytraceLightSrcTexture, LIGHT_PWR (aLightIdx));
669
670     // W component: 0 for infinite light and 1 for point light
671     aLight.xyz -= mix (ZERO, theRay.Origin, aLight.w);
672
673     float aPDF = 1.f / uLightCount;
674
675     if (aLight.w != 0.f) // point light source
676     {
677       float aCenterDst = length (aLight.xyz);
678
679       if (aCenterDst < theHitDistance)
680       {
681         float aVisibility = HandlePointLight (
682           theRay.Direct, normalize (aLight.xyz), aParam.w /* radius */, aCenterDst, aPDF);
683
684         if (aVisibility > 0.f)
685         {
686           theHitDistance = aCenterDst;
687           aTotalRadiance = aParam.rgb;
688
689           thePDF = aPDF;
690         }
691       }
692     }
693     else if (theHitDistance == MAXFLOAT) // directional light source
694     {
695       aTotalRadiance += aParam.rgb * HandleDistantLight (
696         theRay.Direct, aLight.xyz, aParam.w /* angle cosine */, aPDF);
697
698       thePDF += aPDF;
699     }
700   }
701
702   if (thePDF == 0.f && theHitDistance == MAXFLOAT) // light source not found
703   {
704     if (theDepth + uSphereMapForBack == 0) // view ray and map is hidden
705     {
706       aTotalRadiance = BackgroundColor().rgb;
707     }
708     else
709     {
710       aTotalRadiance = FetchEnvironment (Latlong (theRay.Direct)).rgb;
711     }
712   #ifdef THE_SHIFT_sRGB
713     aTotalRadiance = pow (aTotalRadiance, vec3 (2.f));
714   #endif
715   }
716   
717   return aTotalRadiance;
718 }
719
720 #define MIN_THROUGHPUT   vec3 (1.0e-3f)
721 #define MIN_CONTRIBUTION vec3 (1.0e-2f)
722
723 #define MATERIAL_KC(index)           (19 * index + 11)
724 #define MATERIAL_KD(index)           (19 * index + 12)
725 #define MATERIAL_KS(index)           (19 * index + 13)
726 #define MATERIAL_KT(index)           (19 * index + 14)
727 #define MATERIAL_LE(index)           (19 * index + 15)
728 #define MATERIAL_FRESNEL_COAT(index) (19 * index + 16)
729 #define MATERIAL_FRESNEL_BASE(index) (19 * index + 17)
730 #define MATERIAL_ABSORPT_BASE(index) (19 * index + 18)
731
732 //! Enables experimental Russian roulette sampling path termination.
733 //! In most cases, it provides faster image convergence with minimal
734 //! bias, so it is enabled by default.
735 #define RUSSIAN_ROULETTE
736
737 //! Frame step to increase number of bounces. This mode is used
738 //! for interaction with the model, when path length is limited
739 //! for the first samples, and gradually increasing when camera
740 //! is stabilizing.
741 #ifdef ADAPTIVE_SAMPLING
742   #define FRAME_STEP 4
743 #else
744   #define FRAME_STEP 5
745 #endif
746
747 //=======================================================================
748 // function : IsNotZero
749 // purpose  : Checks whether BSDF reflects direct light
750 //=======================================================================
751 bool IsNotZero (in SBSDF theBSDF, in vec3 theThroughput)
752 {
753   vec3 aGlossy = theBSDF.Kc.rgb * step (FLT_EPSILON, theBSDF.Kc.w) +
754                  theBSDF.Ks.rgb * step (FLT_EPSILON, theBSDF.Ks.w);
755
756   return convolve (theBSDF.Kd.rgb + aGlossy, theThroughput) > FLT_EPSILON;
757 }
758
759 //=======================================================================
760 // function : PathTrace
761 // purpose  : Calculates radiance along the given ray
762 //=======================================================================
763 vec4 PathTrace (in SRay theRay, in vec3 theInverse, in int theNbSamples)
764 {
765   float aRaytraceDepth = MAXFLOAT;
766
767   vec3 aRadiance   = ZERO;
768   vec3 aThroughput = UNIT;
769
770   int  aTransfID = 0;     // ID of object transformation
771   bool aInMedium = false; // is the ray inside an object
772
773   float aExpPDF = 1.f;
774   float aImpPDF = 1.f;
775
776   for (int aDepth = 0; aDepth < NB_BOUNCES; ++aDepth)
777   {
778     SIntersect aHit = SIntersect (MAXFLOAT, vec2 (ZERO), ZERO);
779
780     ivec4 aTriIndex = SceneNearestHit (theRay, theInverse, aHit, aTransfID);
781
782     // check implicit path
783     vec3 aLe = IntersectLight (theRay, aDepth, aHit.Time, aExpPDF);
784
785     if (any (greaterThan (aLe, ZERO)) || aTriIndex.x == -1)
786     {
787       float aMIS = (aDepth == 0 || aImpPDF == MAXFLOAT) ? 1.f :
788         aImpPDF * aImpPDF / (aExpPDF * aExpPDF + aImpPDF * aImpPDF);
789
790       aRadiance += aThroughput * aLe * aMIS; break; // terminate path
791     }
792
793     vec3 aInvTransf0 = texelFetch (uSceneTransformTexture, aTransfID + 0).xyz;
794     vec3 aInvTransf1 = texelFetch (uSceneTransformTexture, aTransfID + 1).xyz;
795     vec3 aInvTransf2 = texelFetch (uSceneTransformTexture, aTransfID + 2).xyz;
796
797     // compute geometrical normal
798     aHit.Normal = normalize (vec3 (dot (aInvTransf0, aHit.Normal),
799                                    dot (aInvTransf1, aHit.Normal),
800                                    dot (aInvTransf2, aHit.Normal)));
801
802     theRay.Origin += theRay.Direct * aHit.Time; // get new intersection point
803
804     // evaluate depth on first hit
805     if (aDepth == 0)
806     {
807       vec4 aNDCPoint = uViewMat * vec4 (theRay.Origin, 1.f);
808
809       float aPolygonOffset = PolygonOffset (aHit.Normal, theRay.Origin);
810       aRaytraceDepth = (aNDCPoint.z / aNDCPoint.w + aPolygonOffset * POLYGON_OFFSET_SCALE) * 0.5f + 0.5f;
811     }
812
813     SBSDF aBSDF;
814
815     // fetch BxDF weights
816     aBSDF.Kc = texelFetch (uRaytraceMaterialTexture, MATERIAL_KC (aTriIndex.w));
817     aBSDF.Kd = texelFetch (uRaytraceMaterialTexture, MATERIAL_KD (aTriIndex.w));
818     aBSDF.Ks = texelFetch (uRaytraceMaterialTexture, MATERIAL_KS (aTriIndex.w));
819     aBSDF.Kt = texelFetch (uRaytraceMaterialTexture, MATERIAL_KT (aTriIndex.w));
820
821     // fetch Fresnel reflectance for both layers
822     aBSDF.FresnelCoat = texelFetch (uRaytraceMaterialTexture, MATERIAL_FRESNEL_COAT (aTriIndex.w)).xyz;
823     aBSDF.FresnelBase = texelFetch (uRaytraceMaterialTexture, MATERIAL_FRESNEL_BASE (aTriIndex.w)).xyz;
824
825     vec4 anLE = texelFetch (uRaytraceMaterialTexture, MATERIAL_LE (aTriIndex.w));
826
827     // compute smooth normal (in parallel with fetch)
828     vec3 aNormal = SmoothNormal (aHit.UV, aTriIndex);
829     aNormal = normalize (vec3 (dot (aInvTransf0, aNormal),
830                                dot (aInvTransf1, aNormal),
831                                dot (aInvTransf2, aNormal)));
832
833     SLocalSpace aSpace = buildLocalSpace (aNormal);
834
835 #ifdef USE_TEXTURES
836     if (aBSDF.Kd.w >= 0.0 || aBSDF.Kt.w >= 0.0 || anLE.w >= 0.0)
837     {
838       vec4 aTexCoord = vec4 (SmoothUV (aHit.UV, aTriIndex), 0.f, 1.f);
839       vec4 aTrsfRow1 = texelFetch (uRaytraceMaterialTexture, MATERIAL_TRS1 (aTriIndex.w));
840       vec4 aTrsfRow2 = texelFetch (uRaytraceMaterialTexture, MATERIAL_TRS2 (aTriIndex.w));
841       aTexCoord.st = vec2 (dot (aTrsfRow1, aTexCoord),
842                            dot (aTrsfRow2, aTexCoord));
843
844       if (anLE.w >= 0.0)
845       {
846         anLE.rgb *= textureLod (sampler2D (uTextureSamplers[int (anLE.w)]), aTexCoord.st, 0.0).rgb;
847       }
848       if (aBSDF.Kt.w >= 0.0)
849       {
850         vec2 aTexMetRough = textureLod (sampler2D (uTextureSamplers[int (aBSDF.Kt.w)]), aTexCoord.st, 0.0).bg;
851         float aPbrMetal = aTexMetRough.x;
852         float aPbrRough2 = aTexMetRough.y * aTexMetRough.y;
853         aBSDF.Ks.a *= aPbrRough2;
854         // when using metal-roughness texture, global metalness of material (encoded in FresnelBase) is expected to be 1.0 so that Kd will be 0.0
855         aBSDF.Kd.rgb = aBSDF.FresnelBase * (1.0 - aPbrMetal);
856         aBSDF.FresnelBase *= aPbrMetal;
857       }
858       if (aBSDF.Kd.w >= 0.0)
859       {
860         vec4 aTexColor = textureLod (sampler2D (uTextureSamplers[int (aBSDF.Kd.w)]), aTexCoord.st, 0.0);
861         vec3 aDiff = aTexColor.rgb * aTexColor.a;
862         aBSDF.Kd.rgb *= aDiff;
863         aBSDF.FresnelBase *= aDiff;
864         if (aTexColor.a != 1.0)
865         {
866           // mix transparency BTDF with texture alpha-channel
867           aBSDF.Ks.rgb *= aTexColor.a;
868           aBSDF.Kt.rgb = (UNIT - aTexColor.aaa) + aTexColor.a * aBSDF.Kt.rgb;
869         }
870       }
871     }
872 #endif
873
874     if (uLightCount > 0 && IsNotZero (aBSDF, aThroughput))
875     {
876       aExpPDF = 1.f / uLightCount;
877
878       int aLightIdx = min (int (floor (RandFloat() * uLightCount)), uLightCount - 1);
879
880       vec4 aLight = texelFetch (
881         uRaytraceLightSrcTexture, LIGHT_POS (aLightIdx));
882       vec4 aParam = texelFetch (
883         uRaytraceLightSrcTexture, LIGHT_PWR (aLightIdx));
884
885       // 'w' component is 0 for infinite light and 1 for point light
886       aLight.xyz -= mix (ZERO, theRay.Origin, aLight.w);
887
888       float aDistance = length (aLight.xyz);
889
890       aLight.xyz = SampleLight (aLight.xyz, aDistance,
891         aLight.w == 0.f /* is infinite */, aParam.w /* max cos or radius */, aExpPDF);
892
893       aImpPDF = BsdfPdfLayered (aBSDF,
894         toLocalSpace (-theRay.Direct, aSpace), toLocalSpace (aLight.xyz, aSpace), aThroughput);
895
896       // MIS weight including division by explicit PDF
897       float aMIS = (aExpPDF == MAXFLOAT) ? 1.f : aExpPDF / (aExpPDF * aExpPDF + aImpPDF * aImpPDF);
898
899       vec3 aContrib = aMIS * aParam.rgb /* Le */ * EvalBsdfLayered (
900           aBSDF, toLocalSpace (aLight.xyz, aSpace), toLocalSpace (-theRay.Direct, aSpace));
901
902       if (any (greaterThan (aContrib, MIN_CONTRIBUTION))) // check if light source is important
903       {
904         SRay aShadow = SRay (theRay.Origin + aLight.xyz * uSceneEpsilon, aLight.xyz);
905
906         aShadow.Origin += aHit.Normal * mix (
907           -uSceneEpsilon, uSceneEpsilon, step (0.f, dot (aHit.Normal, aLight.xyz)));
908
909         float aVisibility = SceneAnyHit (aShadow,
910           InverseDirection (aLight.xyz), aLight.w == 0.f ? MAXFLOAT : aDistance);
911
912         aRadiance += aVisibility * (aThroughput * aContrib);
913       }
914     }
915
916     // account for self-emission
917     aRadiance += aThroughput * anLE.rgb;
918
919     if (aInMedium) // handle attenuation
920     {
921       vec4 aScattering = texelFetch (uRaytraceMaterialTexture, MATERIAL_ABSORPT_BASE (aTriIndex.w));
922
923       aThroughput *= exp (-aHit.Time * aScattering.w * (UNIT - aScattering.rgb));
924     }
925
926     vec3 anInput = UNIT; // sampled input direction
927
928     aImpPDF = SampleBsdfLayered (aBSDF,
929       toLocalSpace (-theRay.Direct, aSpace), anInput, aThroughput, aInMedium);
930
931     float aSurvive = float (any (greaterThan (aThroughput, MIN_THROUGHPUT)));
932
933 #ifdef RUSSIAN_ROULETTE
934     aSurvive = aDepth < 3 ? aSurvive : min (dot (LUMA, aThroughput), 0.95f);
935 #endif
936
937     // here, we additionally increase path length for non-diffuse bounces
938     if (RandFloat() > aSurvive || all (lessThan (aThroughput, MIN_THROUGHPUT)) || aDepth >= theNbSamples / FRAME_STEP + step (1.f / M_PI, aImpPDF))
939     {
940       aDepth = INVALID_BOUNCES; // terminate path
941     }
942
943 #ifdef RUSSIAN_ROULETTE
944     aThroughput /= aSurvive;
945 #endif
946
947     anInput = normalize (fromLocalSpace (anInput, aSpace));
948
949     theRay = SRay (theRay.Origin + anInput * uSceneEpsilon +
950       aHit.Normal * mix (-uSceneEpsilon, uSceneEpsilon, step (0.f, dot (aHit.Normal, anInput))), anInput);
951
952     theInverse = InverseDirection (anInput);
953   }
954
955   gl_FragDepth = aRaytraceDepth;
956
957   return vec4 (aRadiance, aRaytraceDepth);
958 }
959
960 #endif