src/Shaders/PathtraceBase.fs

   1 #ifdef PATH_TRACING
   2
   3 ///////////////////////////////////////////////////////////////////////////////////////
   4 // Specific data types
   5
   6 //! Describes local space at the hit point (visualization space).
   7 struct SLocalSpace
   8 {
   9   //! Local X axis.
  10   vec3 AxisX;
  11
  12   //! Local Y axis.
  13   vec3 AxisY;
  14
  15   //! Local Z axis.
  16   vec3 AxisZ;
  17 };
  18
  19 //! Describes material properties (BSDF).
  20 struct SMaterial
  21 {
  22   //! Weight of the Lambertian BRDF.
  23   vec4 Kd;
  24
  25   //! Weight of the reflection BRDF.
  26   vec3 Kr;
  27
  28   //! Weight of the transmission BTDF.
  29   vec3 Kt;
  30
  31   //! Weight of the Blinn BRDF (and roughness).
  32   vec4 Ks;
  33
  34   //! Fresnel coefficients.
  35   vec3 Fresnel;
  36
  37   //! Absorption color and intensity of the media.
  38   vec4 Absorption;
  39 };
  40
  41 ///////////////////////////////////////////////////////////////////////////////////////
  42 // Support subroutines
  43
  44 //=======================================================================
  45 // function : LocalSpace
  46 // purpose  : Generates local space for the given normal
  47 //=======================================================================
  48 SLocalSpace LocalSpace (in vec3 theNormal)
  49 {
  50   vec3 anAxisX = cross (vec3 (0.f, 1.f, 0.f), theNormal);
  51   vec3 anAxisY = cross (vec3 (1.f, 0.f, 0.f), theNormal);
  52
  53   float aSqrLenX = dot (anAxisX, anAxisX);
  54   float aSqrLenY = dot (anAxisY, anAxisY);
  55
  56   if (aSqrLenX > aSqrLenY)
  57   {
  58     anAxisX *= inversesqrt (aSqrLenX);
  59     anAxisY = cross (anAxisX, theNormal);
  60   }
  61   else
  62   {
  63     anAxisY *= inversesqrt (aSqrLenY);
  64     anAxisX = cross (anAxisY, theNormal);
  65   }
  66
  67   return SLocalSpace (anAxisX, anAxisY, theNormal);
  68 }
  69
  70 //=======================================================================
  71 // function : toLocalSpace
  72 // purpose  : Transforms the vector to local space from world space
  73 //=======================================================================
  74 vec3 toLocalSpace (in vec3 theVector, in SLocalSpace theSpace)
  75 {
  76   return vec3 (dot (theVector, theSpace.AxisX),
  77                dot (theVector, theSpace.AxisY),
  78                dot (theVector, theSpace.AxisZ));
  79 }
  80
  81 //=======================================================================
  82 // function : fromLocalSpace
  83 // purpose  : Transforms the vector from local space to world space
  84 //=======================================================================
  85 vec3 fromLocalSpace (in vec3 theVector, in SLocalSpace theSpace)
  86 {
  87   return theVector.x * theSpace.AxisX +
  88          theVector.y * theSpace.AxisY +
  89          theVector.z * theSpace.AxisZ;
  90 }
  91
  92 //=======================================================================
  93 // function : convolve
  94 // purpose  : Performs a linear convolution of the vector components
  95 //=======================================================================
  96 float convolve (in vec3 theVector, in vec3 theFactor)
  97 {
  98   return dot (theVector, theFactor) * (1.f / max (theFactor.x + theFactor.y + theFactor.z, 1e-15f));
  99 }
 100
 101 //=======================================================================
 102 // function : sphericalDirection
 103 // purpose  : Constructs vector from spherical coordinates
 104 //=======================================================================
 105 vec3 sphericalDirection (in float theCosTheta, in float thePhi)
 106 {
 107   float aSinTheta = sqrt (1.f - theCosTheta * theCosTheta);
 108
 109   return vec3 (aSinTheta * cos (thePhi),
 110                aSinTheta * sin (thePhi),
 111                theCosTheta);
 112 }
 113
 114 //=======================================================================
 115 // function : fresnelSchlick
 116 // purpose  : Computes the Fresnel reflection formula using
 117 //            Schlick's approximation.
 118 //=======================================================================
 119 vec3 fresnelSchlick (in float theCosI, in vec3 theSpecularColor)
 120 {
 121   return theSpecularColor + (UNIT - theSpecularColor) * pow (1.f - theCosI, 5.f);
 122 }
 123
 124 //=======================================================================
 125 // function : fresnelDielectric
 126 // purpose  : Computes the Fresnel reflection formula for dielectric in
 127 //            case of circularly polarized light (Based on PBRT code).
 128 //=======================================================================
 129 float fresnelDielectric (in float theCosI,
 130                          in float theCosT,
 131                          in float theEtaI,
 132                          in float theEtaT)
 133 {
 134   float aParl = (theEtaT * theCosI - theEtaI * theCosT) /
 135                 (theEtaT * theCosI + theEtaI * theCosT);
 136
 137   float aPerp = (theEtaI * theCosI - theEtaT * theCosT) /
 138                 (theEtaI * theCosI + theEtaT * theCosT);
 139
 140   return (aParl * aParl + aPerp * aPerp) * 0.5f;
 141 }
 142
 143 #define ENVIRONMENT_IOR 1.f
 144
 145 //=======================================================================
 146 // function : fresnelDielectric
 147 // purpose  : Computes the Fresnel reflection formula for dielectric in
 148 //            case of circularly polarized light (based on PBRT code)
 149 //=======================================================================
 150 float fresnelDielectric (in float theCosI, in float theIndex)
 151 {
 152   float anEtaI = theCosI > 0.f ? 1.f : theIndex;
 153   float anEtaT = theCosI > 0.f ? theIndex : 1.f;
 154
 155   float aSinT = (anEtaI / anEtaT) * sqrt (1.f - theCosI * theCosI);
 156
 157   if (aSinT >= 1.f)
 158   {
 159     return 1.f;
 160   }
 161
 162   float aCosT = sqrt (1.f - aSinT * aSinT);
 163
 164   return fresnelDielectric (abs (theCosI), aCosT, anEtaI, anEtaT);
 165 }
 166
 167 //=======================================================================
 168 // function : fresnelConductor
 169 // purpose  : Computes the Fresnel reflection formula for conductor in case
 170 //            of circularly polarized light (based on PBRT source code)
 171 //=======================================================================
 172 float fresnelConductor (in float theCosI, in float theEta, in float theK)
 173 {
 174   float aTmp = 2.f * theEta * theCosI;
 175
 176   float aTmp1 = theEta * theEta + theK * theK;
 177
 178   float aSPerp = (aTmp1 - aTmp + theCosI * theCosI) /
 179                  (aTmp1 + aTmp + theCosI * theCosI);
 180
 181   float aTmp2 = aTmp1 * theCosI * theCosI;
 182
 183   float aSParl = (aTmp2 - aTmp + 1.f) /
 184                  (aTmp2 + aTmp + 1.f);
 185
 186   return (aSPerp + aSParl) * 0.5f;
 187 }
 188
 189 #define FRESNEL_SCHLICK    -0.5f
 190 #define FRESNEL_CONSTANT   -1.5f
 191 #define FRESNEL_CONDUCTOR  -2.5f
 192 #define FRESNEL_DIELECTRIC -3.5f
 193
 194 //=======================================================================
 195 // function : fresnelMedia
 196 // purpose  : Computes the Fresnel reflection formula for general medium
 197 //            in case of circularly polarized light.
 198 //=======================================================================
 199 vec3 fresnelMedia (in float theCosI, in vec3 theFresnelCoeffs)
 200 {
 201   if (theFresnelCoeffs.x > FRESNEL_SCHLICK)
 202   {
 203     return fresnelSchlick (abs (theCosI), theFresnelCoeffs);
 204   }
 205
 206   if (theFresnelCoeffs.x > FRESNEL_CONSTANT)
 207   {
 208     return vec3 (theFresnelCoeffs.z);
 209   }
 210
 211   if (theFresnelCoeffs.x > FRESNEL_CONDUCTOR)
 212   {
 213     return vec3 (fresnelConductor (abs (theCosI), theFresnelCoeffs.y, theFresnelCoeffs.z));
 214   }
 215
 216   return vec3 (fresnelDielectric (theCosI, theFresnelCoeffs.y));
 217 }
 218
 219 //=======================================================================
 220 // function : transmitted
 221 // purpose  : Computes transmitted direction in tangent space
 222 //            (in case of TIR returned result is undefined!)
 223 //=======================================================================
 224 void transmitted (in float theIndex, in vec3 theIncident, out vec3 theTransmit)
 225 {
 226   // Compute relative index of refraction
 227   float anEta = (theIncident.z > 0.f) ? 1.f / theIndex : theIndex;
 228
 229   // Handle total internal reflection for transmission
 230   float aSinT2 = anEta * anEta * (1.f - theIncident.z * theIncident.z);
 231
 232   // Compute transmitted ray direction
 233   float aCosT = sqrt (1.f - min (aSinT2, 1.f)) * (theIncident.z > 0.f ? -1.f : 1.f);
 234
 235   theTransmit = normalize (vec3 (-anEta * theIncident.x,
 236                                  -anEta * theIncident.y,
 237                                   aCosT));
 238 }
 239
 240 //////////////////////////////////////////////////////////////////////////////////////////////
 241 // Handlers and samplers for materials
 242 //////////////////////////////////////////////////////////////////////////////////////////////
 243
 244 //=======================================================================
 245 // function : handleLambertianReflection
 246 // purpose  : Handles Lambertian BRDF, with cos(N, PSI)
 247 //=======================================================================
 248 float handleLambertianReflection (in vec3 theInput, in vec3 theOutput)
 249 {
 250   return max (0.f, theInput.z) * (1.f / M_PI);
 251 }
 252
 253 //=======================================================================
 254 // function : handleBlinnReflection
 255 // purpose  : Handles Blinn glossy BRDF, with cos(N, PSI)
 256 //=======================================================================
 257 vec3 handleBlinnReflection (in vec3 theInput, in vec3 theOutput, in vec3 theFresnelCoeffs, in float theExponent)
 258 {
 259   vec3 aWeight = ZERO;
 260
 261   // Compute half-angle vector
 262   vec3 aHalf = theInput + theOutput;
 263
 264   if (aHalf.z < 0.f)
 265     aHalf = -aHalf;
 266
 267   float aLength = dot (aHalf, aHalf);
 268
 269   if (aLength <= 0.f)
 270     return ZERO;
 271
 272   aHalf *= inversesqrt (aLength);
 273
 274   // Compute Fresnel reflectance
 275   float aCosDelta = dot (theOutput, aHalf);
 276
 277   vec3 aFresnel = fresnelMedia (aCosDelta, theFresnelCoeffs);
 278
 279   // Compute fraction of microfacets that reflect light
 280   float aCosThetaH = max (0.f, aHalf.z);
 281
 282   float aFraction = (theExponent + 2.f) * (M_PI / 2.f) * pow (aCosThetaH, theExponent);
 283
 284   // Compute geometry attenuation term (already includes cos)
 285   float aCosThetaI = max (0.f, theInput.z);
 286   float aCosThetaO = max (0.f, theOutput.z);
 287
 288   float aGeom = min (1.f, 2.f * aCosThetaH / max (0.f, aCosDelta) * min (aCosThetaO, aCosThetaI));
 289
 290   return aCosThetaO < 1.0e-3f ? ZERO :
 291     aFraction * aGeom / (4.f * aCosThetaO) * aFresnel;
 292 }
 293
 294 //=======================================================================
 295 // function : handleMaterial
 296 // purpose  : Returns BSDF value for specified material, with cos(N, PSI)
 297 //=======================================================================
 298 vec3 handleMaterial (in SMaterial theMaterial, in vec3 theInput, in vec3 theOutput)
 299 {
 300   return theMaterial.Kd.rgb * handleLambertianReflection (theInput, theOutput) +
 301     theMaterial.Ks.rgb * handleBlinnReflection (theInput, theOutput, theMaterial.Fresnel, theMaterial.Ks.w);
 302 }
 303
 304 //=======================================================================
 305 // function : sampleLambertianReflection
 306 // purpose  : Samples Lambertian BRDF, W = BRDF * cos(N, PSI) / PDF(PSI)
 307 //=======================================================================
 308 void sampleLambertianReflection (in vec3 theOutput, out vec3 theInput)
 309 {
 310   float aKsi1 = RandFloat();
 311   float aKsi2 = RandFloat();
 312
 313   float aTemp = sqrt (aKsi2);
 314
 315   theInput = vec3 (aTemp * cos (2.f * M_PI * aKsi1),
 316                    aTemp * sin (2.f * M_PI * aKsi1),
 317                    sqrt (1.f - aKsi2));
 318
 319   theInput.z = mix (-theInput.z, theInput.z, step (0.f, theOutput.z));
 320 }
 321
 322 // Types of bounces
 323 #define NON_SPECULAR_BOUNCE 0
 324 #define SPEC_REFLECT_BOUNCE 1
 325 #define SPEC_REFRACT_BOUNCE 2
 326
 327 #define IS_NON_SPEC_BOUNCE(theBounce) (theBounce == 0)
 328 #define IS_ANY_SPEC_BOUNCE(theBounce) (theBounce != 0)
 329 #define IS_REFL_SPEC_BOUNCE(theBounce) (theBounce == 1)
 330 #define IS_REFR_SPEC_BOUNCE(theBounce) (theBounce == 2)
 331
 332 //=======================================================================
 333 // function : sampleSpecularTransmission
 334 // purpose  : Samples specular BTDF, W = BRDF * cos(N, PSI) / PDF(PSI)
 335 //=======================================================================
 336 vec3 sampleSpecularTransmission (in vec3 theOutput, out vec3 theInput,
 337   out int theBounce, in vec3 theWeight, in vec3 theFresnelCoeffs)
 338 {
 339   vec3 aFresnel = fresnelMedia (theOutput.z, theFresnelCoeffs);
 340
 341   float aProbability = convolve (aFresnel, theWeight);
 342
 343   // Check if transmission takes place
 344   theBounce = RandFloat() <= aProbability ?
 345     SPEC_REFLECT_BOUNCE : SPEC_REFRACT_BOUNCE;
 346
 347   // Sample input direction
 348   if (theBounce == SPEC_REFLECT_BOUNCE)
 349   {
 350     theInput = vec3 (-theOutput.x,
 351                      -theOutput.y,
 352                       theOutput.z);
 353
 354     theWeight = aFresnel * (1.f / aProbability);
 355   }
 356   else
 357   {
 358     transmitted (theFresnelCoeffs.y, theOutput, theInput);
 359
 360     theWeight = (UNIT - aFresnel) * (1.f / (1.f - aProbability));
 361   }
 362
 363   return theWeight;
 364 }
 365
 366 //=======================================================================
 367 // function : sampleSpecularReflection
 368 // purpose  : Samples specular BRDF, W = BRDF * cos(N, PSI) / PDF(PSI)
 369 //=======================================================================
 370 vec3 sampleSpecularReflection (in vec3 theOutput, out vec3 theInput, in vec3 theFresnelCoeffs)
 371 {
 372   // Sample input direction
 373   theInput = vec3 (-theOutput.x,
 374                    -theOutput.y,
 375                     theOutput.z);
 376
 377   return fresnelMedia (theOutput.z, theFresnelCoeffs);
 378 }
 379
 380 #define MIN_COS 1.0e-20f
 381
 382 //=======================================================================
 383 // function : sampleBlinnReflection
 384 // purpose  : Samples Blinn BRDF, W = BRDF * cos(N, PSI) / PDF(PSI)
 385 //            The BRDF is a product of three main terms, D, G, and F,
 386 //            which is then divided by two cosine terms. Here we perform
 387 //            importance sample the D part of the Blinn model; trying to
 388 //            develop a sampling procedure that accounted for all of the
 389 //            terms would be complex, and it is the D term that accounts
 390 //            for most of the variation.
 391 //=======================================================================
 392 vec3 sampleBlinnReflection (in vec3 theOutput, out vec3 theInput, in vec3 theFresnelCoeffs, in float theExponent)
 393 {
 394   vec3 aWeight = ZERO;
 395
 396   // Generate two random variables
 397   float aKsi1 = RandFloat();
 398   float aKsi2 = RandFloat();
 399
 400   // Compute sampled half-angle vector for Blinn distribution
 401   float aCosThetaH = pow (aKsi1, 1.f / (theExponent + 1.f));
 402
 403   vec3 aHalf = sphericalDirection (aCosThetaH, aKsi2 * 2.f * M_PI);
 404
 405   if (aHalf.z < 0)
 406   {
 407     aHalf = -aHalf;
 408   }
 409
 410   // Compute incident direction by reflecting about half-vector
 411   float aCosDelta = dot (theOutput, aHalf);
 412
 413   vec3 anInput = 2.f * aCosDelta * aHalf - theOutput;
 414
 415   if (theOutput.z * anInput.z <= 0.f)
 416   {
 417     return ZERO;
 418   }
 419
 420   theInput = anInput;
 421
 422   // Compute Fresnel reflectance
 423   vec3 aFresnel = fresnelMedia (aCosDelta, theFresnelCoeffs);
 424
 425   // Compute geometry attenuation term
 426   float aCosThetaI = max (MIN_COS, theInput.z);
 427   float aCosThetaO = max (MIN_COS, theOutput.z);
 428
 429   float aGeom = min (max (MIN_COS, aCosDelta), 2.f * aCosThetaH * min (aCosThetaO, aCosThetaI));
 430
 431   // Compute weight of the ray sample
 432   return aFresnel * ((theExponent + 2.f) / (theExponent + 1.f) * aGeom / aCosThetaO);
 433 }
 434
 435 //=======================================================================
 436 // function : sampleMaterial
 437 // purpose  : Samples specified composite material (BSDF)
 438 //=======================================================================
 439 void sampleMaterial (in SMaterial theMaterial,
 440                      in vec3      theOutput,
 441                      out vec3     theInput,
 442                      inout vec3   theWeight,
 443                      inout int    theBounce)
 444 {
 445   // Compute the probability of ray reflection
 446   float aPd = convolve (theMaterial.Kd.rgb, theWeight);
 447   float aPs = convolve (theMaterial.Ks.rgb, theWeight);
 448   float aPr = convolve (theMaterial.Kr.rgb, theWeight);
 449   float aPt = convolve (theMaterial.Kt.rgb, theWeight);
 450
 451   float aReflection = aPd + aPs + aPr + aPt;
 452
 453   // Choose BSDF component to sample
 454   float aKsi = aReflection * RandFloat();
 455
 456   theBounce = NON_SPECULAR_BOUNCE;
 457
 458   if (aKsi < aPd) // diffuse reflection
 459   {
 460     sampleLambertianReflection (theOutput, theInput);
 461
 462     theWeight *= theMaterial.Kd.rgb * (aReflection / aPd);
 463   }
 464   else if (aKsi < aPd + aPs) //  glossy reflection
 465   {
 466     theWeight *= theMaterial.Ks.rgb * (aReflection / aPs) *
 467       sampleBlinnReflection (theOutput, theInput, theMaterial.Fresnel, theMaterial.Ks.w);
 468   }
 469   else if (aKsi < aPd + aPs + aPr) //  specular reflection
 470   {
 471     theWeight *= theMaterial.Kr.rgb * (aReflection / aPr) *
 472       sampleSpecularReflection (theOutput, theInput, theMaterial.Fresnel);
 473
 474     theBounce = SPEC_REFLECT_BOUNCE; // specular bounce
 475   }
 476   else //  specular transmission
 477   {
 478     theWeight *= theMaterial.Kt.rgb * (aReflection / aPt) *
 479       sampleSpecularTransmission (theOutput, theInput, theBounce, theWeight, theMaterial.Fresnel);
 480   }
 481
 482   // path termination for extra small weights
 483   theWeight = mix (theWeight, ZERO, float (aReflection < 1e-3f));
 484 }
 485
 486 //////////////////////////////////////////////////////////////////////////////////////////////
 487 // Handlers and samplers for light sources
 488 //////////////////////////////////////////////////////////////////////////////////////////////
 489
 490 //=======================================================================
 491 // function : handlePointLight
 492 // purpose  :
 493 //=======================================================================
 494 float handlePointLight (in vec3 theInput, in vec3 theToLight, in float theRadius, in float theDistance)
 495 {
 496   float aDistance = dot (theToLight, theToLight);
 497
 498   float aCosMax = inversesqrt (1.f + theRadius * theRadius / aDistance);
 499
 500   return float (aDistance < theDistance * theDistance) *
 501     step (aCosMax, dot (theToLight, theInput) * inversesqrt (aDistance));
 502 }
 503
 504 //=======================================================================
 505 // function : handleDirectLight
 506 // purpose  :
 507 //=======================================================================
 508 float handleDirectLight (in vec3 theInput, in vec3 theToLight, in float theCosMax)
 509 {
 510   return step (theCosMax, dot (theInput, theToLight));
 511 }
 512
 513 //=======================================================================
 514 // function : sampleLight
 515 // purpose  : General sampling function for directional and point lights
 516 //=======================================================================
 517 vec3 sampleLight (in vec3 theToLight, inout float theDistance, in bool isInfinite, in float theSmoothness, inout float thePDF)
 518 {
 519   SLocalSpace aSpace = LocalSpace (theToLight * (1.f / theDistance));
 520
 521   // for point lights smoothness defines radius
 522   float aCosMax = isInfinite ? theSmoothness :
 523     inversesqrt (1.f + theSmoothness * theSmoothness / (theDistance * theDistance));
 524
 525   float aKsi1 = RandFloat();
 526   float aKsi2 = RandFloat();
 527
 528   float aTmp = 1.f - aKsi2 * (1.f - aCosMax);
 529
 530   vec3 anInput = vec3 (cos (2.f * M_PI * aKsi1),
 531                        sin (2.f * M_PI * aKsi1),
 532                        aTmp);
 533
 534   anInput.xy *= sqrt (1.f - aTmp * aTmp);
 535
 536   thePDF *= (aCosMax < 1.f) ? 1.f / (2.f * M_PI) / (1.f - aCosMax) : 1.f;
 537
 538   return normalize (fromLocalSpace (anInput, aSpace));
 539 }
 540
 541 // =======================================================================
 542 // function : Latlong
 543 // purpose  : Converts world direction to environment texture coordinates
 544 // =======================================================================
 545 vec2 Latlong (in vec3 thePoint)
 546 {
 547   float aPsi = acos (-thePoint.z);
 548
 549   float aPhi = atan (thePoint.y, thePoint.x) + M_PI;
 550
 551   return vec2 (aPhi * 0.1591549f,
 552                aPsi * 0.3183098f);
 553 }
 554
 555 // =======================================================================
 556 // function: intersectLight
 557 // purpose : Checks intersections with light sources
 558 // =======================================================================
 559 vec3 intersectLight (in SRay theRay, in bool isViewRay, in int theBounce, in float theDistance)
 560 {
 561   vec3 aRadiance = ZERO;
 562
 563   if ((isViewRay || IS_REFR_SPEC_BOUNCE(theBounce)) && uSphereMapForBack == 0)
 564   {
 565     aRadiance = BackgroundColor().xyz;
 566   }
 567   else
 568   {
 569     aRadiance = FetchEnvironment (Latlong (theRay.Direct)).xyz;
 570   }
 571
 572   // Apply gamma correction (gamma is 2)
 573   aRadiance = aRadiance * aRadiance * float (theDistance == MAXFLOAT);
 574
 575   for (int aLightIdx = 0; aLightIdx < uLightCount && (isViewRay || IS_ANY_SPEC_BOUNCE(theBounce)); ++aLightIdx)
 576   {
 577     vec4 aLight = texelFetch (
 578       uRaytraceLightSrcTexture, LIGHT_POS (aLightIdx));
 579     vec4 aParam = texelFetch (
 580       uRaytraceLightSrcTexture, LIGHT_PWR (aLightIdx));
 581
 582     if (aLight.w != 0.f) // point light source
 583     {
 584       aRadiance += aParam.rgb * handlePointLight (
 585         theRay.Direct, aLight.xyz - theRay.Origin, aParam.w /* radius */, theDistance);
 586     }
 587     else if (theDistance == MAXFLOAT) // directional light source
 588     {
 589       aRadiance += aParam.rgb * handleDirectLight (theRay.Direct, aLight.xyz, aParam.w /* angle cosine */);
 590     }
 591   }
 592
 593   return aRadiance;
 594 }
 595
 596 #define MIN_THROUGHPUT   vec3 (0.02f)
 597 #define MIN_CONTRIBUTION vec3 (0.01f)
 598
 599 #define MATERIAL_KD(index)      (18 * index + 11)
 600 #define MATERIAL_KR(index)      (18 * index + 12)
 601 #define MATERIAL_KT(index)      (18 * index + 13)
 602 #define MATERIAL_KS(index)      (18 * index + 14)
 603 #define MATERIAL_LE(index)      (18 * index + 15)
 604 #define MATERIAL_FRESNEL(index) (18 * index + 16)
 605 #define MATERIAL_ABSORPT(index) (18 * index + 17)
 606
 607 // Enables expiremental russian roulette sampling
 608 #define RUSSIAN_ROULETTE
 609
 610 //=======================================================================
 611 // function : PathTrace
 612 // purpose  : Calculates radiance along the given ray
 613 //=======================================================================
 614 vec4 PathTrace (in SRay theRay, in vec3 theInverse)
 615 {
 616   float aRaytraceDepth = MAXFLOAT;
 617
 618   vec3 aRadiance   = ZERO;
 619   vec3 aThroughput = UNIT;
 620
 621   int aBounce = 0; // type of previous hit point
 622   int aTrsfId = 0; // offset of object transform
 623
 624   bool isInMedium = false;
 625
 626   for (int aDepth = 0; aDepth < NB_BOUNCES; ++aDepth)
 627   {
 628     SIntersect aHit = SIntersect (MAXFLOAT, vec2 (ZERO), ZERO);
 629
 630     ivec4 aTriIndex = SceneNearestHit (theRay, theInverse, aHit, aTrsfId);
 631
 632     // check implicit path
 633     vec3 aLe = intersectLight (theRay,
 634       aDepth == 0 /* is view ray */, aBounce, aHit.Time);
 635
 636     if (any (greaterThan (aLe, ZERO)) || aTriIndex.x == -1)
 637     {
 638       aRadiance += aThroughput * aLe; break; // terminate path
 639     }
 640
 641     vec3 aInvTransf0 = texelFetch (uSceneTransformTexture, aTrsfId + 0).xyz;
 642     vec3 aInvTransf1 = texelFetch (uSceneTransformTexture, aTrsfId + 1).xyz;
 643     vec3 aInvTransf2 = texelFetch (uSceneTransformTexture, aTrsfId + 2).xyz;
 644
 645     // compute geometrical normal
 646     aHit.Normal = normalize (vec3 (dot (aInvTransf0, aHit.Normal),
 647                                    dot (aInvTransf1, aHit.Normal),
 648                                    dot (aInvTransf2, aHit.Normal)));
 649
 650     theRay.Origin += theRay.Direct * aHit.Time; // get new intersection point
 651
 652     // Evaluate depth on first hit
 653     if (aDepth == 0)
 654     {
 655       // For polygons that are parallel to the screen plane, the depth slope
 656       // is equal to 1, resulting in small polygon offset. For polygons that
 657       // that are at a large angle to the screen, the depth slope tends to 1,
 658       // resulting in a larger polygon offset
 659       float aPolygonOffset = uSceneEpsilon * EPS_SCALE /
 660         max (abs (dot (theRay.Direct, aHit.Normal)), MIN_SLOPE);
 661
 662       // Hit point in NDC-space [-1,1] (the polygon offset is applied in the world space)
 663       vec4 aNDCPoint = uViewMat * vec4 (theRay.Origin + theRay.Direct * aPolygonOffset, 1.f);
 664
 665       aRaytraceDepth = (aNDCPoint.z / aNDCPoint.w) * 0.5f + 0.5f;
 666     }
 667
 668     // fetch material (BSDF)
 669     SMaterial aMaterial = SMaterial (
 670       vec4 (texelFetch (uRaytraceMaterialTexture, MATERIAL_KD      (aTriIndex.w))),
 671       vec3 (texelFetch (uRaytraceMaterialTexture, MATERIAL_KR      (aTriIndex.w))),
 672       vec3 (texelFetch (uRaytraceMaterialTexture, MATERIAL_KT      (aTriIndex.w))),
 673       vec4 (texelFetch (uRaytraceMaterialTexture, MATERIAL_KS      (aTriIndex.w))),
 674       vec3 (texelFetch (uRaytraceMaterialTexture, MATERIAL_FRESNEL (aTriIndex.w))),
 675       vec4 (texelFetch (uRaytraceMaterialTexture, MATERIAL_ABSORPT (aTriIndex.w))));
 676
 677 #ifdef USE_TEXTURES
 678     if (aMaterial.Kd.w >= 0.f)
 679     {
 680       vec4 aTexCoord = vec4 (SmoothUV (aHit.UV, aTriIndex), 0.f, 1.f);
 681
 682       vec4 aTrsfRow1 = texelFetch (
 683         uRaytraceMaterialTexture, MATERIAL_TRS1 (aTriIndex.w));
 684       vec4 aTrsfRow2 = texelFetch (
 685         uRaytraceMaterialTexture, MATERIAL_TRS2 (aTriIndex.w));
 686
 687       aTexCoord.st = vec2 (dot (aTrsfRow1, aTexCoord),
 688                            dot (aTrsfRow2, aTexCoord));
 689
 690       vec3 aTexColor = textureLod (
 691         sampler2D (uTextureSamplers[int (aMaterial.Kd.w)]), aTexCoord.st, 0.f).rgb;
 692
 693       aMaterial.Kd.rgb *= aTexColor * aTexColor; // de-gamma correction (for gamma = 2)
 694     }
 695 #endif
 696
 697     // compute smooth normal
 698     vec3 aNormal = SmoothNormal (aHit.UV, aTriIndex);
 699
 700     aNormal = normalize (vec3 (dot (aInvTransf0, aNormal),
 701                                dot (aInvTransf1, aNormal),
 702                                dot (aInvTransf2, aNormal)));
 703
 704     SLocalSpace aSpace = LocalSpace (aNormal);
 705
 706     // account for self-emission (not stored in the material)
 707     aRadiance += aThroughput * texelFetch (
 708       uRaytraceMaterialTexture, MATERIAL_LE (aTriIndex.w)).rgb;
 709
 710     if (uLightCount > 0 && convolve (aMaterial.Kd.rgb + aMaterial.Ks.rgb, aThroughput) > 0.f)
 711     {
 712       int aLightIdx = min (int (floor (RandFloat() * uLightCount)), uLightCount - 1);
 713
 714       vec4 aLight = texelFetch (
 715         uRaytraceLightSrcTexture, LIGHT_POS (aLightIdx));
 716       vec4 aParam = texelFetch (
 717         uRaytraceLightSrcTexture, LIGHT_PWR (aLightIdx));
 718
 719       // 'w' component is 0 for infinite light and 1 for point light
 720       aLight.xyz -= mix (ZERO, theRay.Origin, aLight.w);
 721
 722       float aPDF = 1.f / uLightCount, aDistance = length (aLight.xyz);
 723
 724       aLight.xyz = sampleLight (aLight.xyz, aDistance,
 725         aLight.w == 0.f /* is infinite */, aParam.w /* max cos or radius */, aPDF);
 726
 727       vec3 aContrib = (1.f / aPDF) * aParam.rgb /* Le */ * handleMaterial (
 728           aMaterial, toLocalSpace (aLight.xyz, aSpace), toLocalSpace (-theRay.Direct, aSpace));
 729
 730       if (any (greaterThan (aContrib, MIN_CONTRIBUTION))) // first check if light source is important
 731       {
 732         SRay aShadow = SRay (theRay.Origin + aLight.xyz * uSceneEpsilon, aLight.xyz);
 733
 734         aShadow.Origin += aHit.Normal * mix (
 735           -uSceneEpsilon, uSceneEpsilon, step (0.f, dot (aHit.Normal, aLight.xyz)));
 736
 737         float aVisibility = SceneAnyHit (aShadow,
 738           InverseDirection (aLight.xyz), aLight.w == 0.f ? MAXFLOAT : aDistance);
 739
 740         aRadiance += aVisibility * aThroughput * aContrib;
 741       }
 742     }
 743
 744     vec3 anInput;
 745
 746     sampleMaterial (aMaterial,
 747       toLocalSpace (-theRay.Direct, aSpace), anInput, aThroughput, aBounce);
 748
 749     if (isInMedium)
 750     {
 751       aThroughput *= exp (-aHit.Time *
 752         aMaterial.Absorption.w * (UNIT - aMaterial.Absorption.rgb));
 753     }
 754
 755     isInMedium = IS_REFR_SPEC_BOUNCE(aBounce) ? !isInMedium : isInMedium;
 756
 757 #ifndef RUSSIAN_ROULETTE
 758     if (all (lessThan (aThroughput, MIN_THROUGHPUT)))
 759     {
 760       aDepth = INVALID_BOUNCES; // terminate path
 761     }
 762 #else
 763     float aSurvive = aDepth < 3 ? 1.f : min (dot (LUMA, aThroughput), 0.95f);
 764
 765     if (RandFloat() > aSurvive)
 766     {
 767       aDepth = INVALID_BOUNCES; // terminate path
 768     }
 769
 770     aThroughput /= aSurvive;
 771 #endif
 772
 773     anInput = normalize (fromLocalSpace (anInput, aSpace));
 774
 775     theRay = SRay (theRay.Origin + anInput * uSceneEpsilon +
 776       aHit.Normal * mix (-uSceneEpsilon, uSceneEpsilon, step (0.f, dot (aHit.Normal, anInput))), anInput);
 777
 778     theInverse = InverseDirection (anInput);
 779   }
 780
 781   gl_FragDepth = aRaytraceDepth;
 782
 783   return vec4 (aRadiance, aRaytraceDepth);
 784 }
 785
 786 #endif