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