Mindgames/Library/PackageCache/com.unity.render-pipelines.core@6.9.0/ShaderLibrary/VolumeRendering.hlsl

   1 #ifndef UNITY_VOLUME_RENDERING_INCLUDED
   2 #define UNITY_VOLUME_RENDERING_INCLUDED
   3
   4 // Reminder:
   5 // OpticalDepth(x, y) = Integral{x, y}{Extinction(t) dt}
   6 // Transmittance(x, y) = Exp(-OpticalDepth(x, y))
   7 // Transmittance(x, z) = Transmittance(x, y) * Transmittance(y, z)
   8 // Integral{a, b}{Transmittance(0, t) dt} = Transmittance(0, a) * Integral{a, b}{Transmittance(0, t - a) dt}
   9
  10 real TransmittanceFromOpticalDepth(real opticalDepth)
  11 {
  12     return exp(-opticalDepth);
  13 }
  14
  15 real OpacityFromOpticalDepth(real opticalDepth)
  16 {
  17     return 1 - TransmittanceFromOpticalDepth(opticalDepth);
  18 }
  19
  20 //
  21 // ---------------------------------- Deep Pixel Compositing ---------------------------------------
  22 //
  23
  24 // TODO: it would be good to improve the perf and numerical stability
  25 // of approximations below by finding a polynomial approximation.
  26
  27 // input = {radiance, opacity}
  28 // Note that opacity must be less than 1 (not fully opaque).
  29 real4 LinearizeRGBA(real4 value)
  30 {
  31     // See "Deep Compositing Using Lie Algebras".
  32     // log(A) = {OpticalDepthFromOpacity(A.a) / A.a * A.rgb, -OpticalDepthFromOpacity(A.a)}.
  33     // We drop redundant negations.
  34     real a = value.a;
  35     real d = -log(1 - a);
  36     real r = (a >= REAL_EPS) ? (d * rcp(a)) : 1; // Prevent numerical explosion
  37     return real4(r * value.rgb, d);
  38 }
  39
  40 // input = {radiance, optical_depth}
  41 // Note that opacity must be less than 1 (not fully opaque).
  42 real4 LinearizeRGBD(real4 value)
  43 {
  44     // See "Deep Compositing Using Lie Algebras".
  45     // log(A) = {A.a / OpacityFromOpticalDepth(A.a) * A.rgb, -A.a}.
  46     // We drop redundant negations.
  47     real d = value.a;
  48     real a = 1 - exp(-d);
  49     real r = (a >= REAL_EPS) ? (d * rcp(a)) : 1; // Prevent numerical explosion
  50     return real4(r * value.rgb, d);
  51 }
  52
  53 // output = {radiance, opacity}
  54 // Note that opacity must be less than 1 (not fully opaque).
  55 real4 DelinearizeRGBA(real4 value)
  56 {
  57     // See "Deep Compositing Using Lie Algebras".
  58     // exp(B) = {OpacityFromOpticalDepth(-B.a) / -B.a * B.rgb, OpacityFromOpticalDepth(-B.a)}.
  59     // We drop redundant negations.
  60     real d = value.a;
  61     real a = 1 - exp(-d);
  62     real i = (a >= REAL_EPS) ? (a * rcp(d)) : 1; // Prevent numerical explosion
  63     return real4(i * value.rgb, a);
  64 }
  65
  66 // input = {radiance, optical_depth}
  67 // Note that opacity must be less than 1 (not fully opaque).
  68 real4 DelinearizeRGBD(real4 value)
  69 {
  70     // See "Deep Compositing Using Lie Algebras".
  71     // exp(B) = {OpacityFromOpticalDepth(-B.a) / -B.a * B.rgb, -B.a}.
  72     // We drop redundant negations.
  73     real d = value.a;
  74     real a = 1 - exp(-d);
  75     real i = (a >= REAL_EPS) ? (a * rcp(d)) : 1; // Prevent numerical explosion
  76     return real4(i * value.rgb, d);
  77 }
  78
  79 //
  80 // ----------------------------- Homogeneous Participating Media -----------------------------------
  81 //
  82
  83 real OpticalDepthHomogeneousMedium(real extinction, real intervalLength)
  84 {
  85     return extinction * intervalLength;
  86 }
  87
  88 real TransmittanceHomogeneousMedium(real extinction, real intervalLength)
  89 {
  90     return TransmittanceFromOpticalDepth(OpticalDepthHomogeneousMedium(extinction, intervalLength));
  91 }
  92
  93 // Integral{a, b}{TransmittanceFromOpticalDepth(0, t - a) dt}.
  94 real TransmittanceIntegralHomogeneousMedium(real extinction, real intervalLength)
  95 {
  96     // Note: when multiplied by the extinction coefficient, it becomes
  97     // Albedo * (1 - TransmittanceFromOpticalDepth(d)) = Albedo * Opacity(d).
  98     return rcp(extinction) - rcp(extinction) * exp(-extinction * intervalLength);
  99 }
 100
 101 //
 102 // ----------------------------------- Height Fog --------------------------------------------------
 103 //
 104
 105 // Can be used to scale base extinction and scattering coefficients.
 106 real ComputeHeightFogMultiplier(real height, real baseHeight, real2 heightExponents)
 107 {
 108     real h    = max(height - baseHeight, 0);
 109     real rcpH = heightExponents.x;
 110
 111     return exp(-h * rcpH);
 112 }
 113
 114 // Optical depth between two endpoints.
 115 real OpticalDepthHeightFog(real baseExtinction, real baseHeight, real2 heightExponents,
 116                            real cosZenith, real startHeight, real intervalLength)
 117 {
 118     // Height fog is composed of two slices of optical depth:
 119     // - homogeneous fog below 'baseHeight': d = k * t
 120     // - exponential fog above 'baseHeight': d = Integrate[k * e^(-(h + z * x) / H) dx, {x, 0, t}]
 121
 122     real H          = heightExponents.y;
 123     real rcpH       = heightExponents.x;
 124     real Z          = cosZenith;
 125     real absZ       = max(abs(cosZenith), REAL_EPS);
 126     real rcpAbsZ    = rcp(absZ);
 127
 128     real endHeight  = startHeight + intervalLength * Z;
 129     real minHeight  = min(startHeight, endHeight);
 130     real h          = max(minHeight - baseHeight, 0);
 131
 132     real homFogDist = clamp((baseHeight - minHeight) * rcpAbsZ, 0, intervalLength);
 133     real expFogDist = intervalLength - homFogDist;
 134     real expFogMult = exp(-h * rcpH) * (1 - exp(-expFogDist * absZ * rcpH)) * (rcpAbsZ * H);
 135
 136     return baseExtinction * (homFogDist + expFogMult);
 137 }
 138
 139 // This version of the function assumes the interval of infinite length.
 140 real OpticalDepthHeightFog(real baseExtinction, real baseHeight, real2 heightExponents,
 141                            real cosZenith, real startHeight)
 142 {
 143     real H          = heightExponents.y;
 144     real rcpH       = heightExponents.x;
 145     real Z          = cosZenith;
 146     real absZ       = max(abs(cosZenith), REAL_EPS);
 147     real rcpAbsZ    = rcp(absZ);
 148
 149     real minHeight  = (Z >= 0) ? startHeight : -rcp(REAL_EPS);
 150     real h          = max(minHeight - baseHeight, 0);
 151
 152     real homFogDist = max((baseHeight - minHeight) * rcpAbsZ, 0);
 153     real expFogMult = exp(-h * rcpH) * (rcpAbsZ * H);
 154
 155     return baseExtinction * (homFogDist + expFogMult);
 156 }
 157
 158 real TransmittanceHeightFog(real baseExtinction, real baseHeight, real2 heightExponents,
 159                             real cosZenith, real startHeight, real intervalLength)
 160 {
 161     real od = OpticalDepthHeightFog(baseExtinction, baseHeight, heightExponents,
 162                                     cosZenith, startHeight, intervalLength);
 163     return TransmittanceFromOpticalDepth(od);
 164 }
 165
 166 real TransmittanceHeightFog(real baseExtinction, real baseHeight, real2 heightExponents,
 167                             real cosZenith, real startHeight)
 168 {
 169     real od = OpticalDepthHeightFog(baseExtinction, baseHeight, heightExponents,
 170                                     cosZenith, startHeight);
 171     return TransmittanceFromOpticalDepth(od);
 172 }
 173
 174 //
 175 // ----------------------------------- Phase Functions ---------------------------------------------
 176 //
 177
 178 real IsotropicPhaseFunction()
 179 {
 180     return INV_FOUR_PI;
 181 }
 182
 183 real HenyeyGreensteinPhasePartConstant(real anisotropy)
 184 {
 185     real g = anisotropy;
 186
 187     return INV_FOUR_PI * (1 - g * g);
 188 }
 189
 190 real HenyeyGreensteinPhasePartVarying(real anisotropy, real cosTheta)
 191 {
 192     real g = anisotropy;
 193     real x = 1 + g * g - 2 * g * cosTheta;
 194     real f = rsqrt(max(x, REAL_EPS)); // x^(-1/2)
 195
 196     return f * f * f; // x^(-3/2)
 197 }
 198
 199 real HenyeyGreensteinPhaseFunction(real anisotropy, real cosTheta)
 200 {
 201     return HenyeyGreensteinPhasePartConstant(anisotropy) *
 202            HenyeyGreensteinPhasePartVarying(anisotropy, cosTheta);
 203 }
 204
 205 real CornetteShanksPhasePartConstant(real anisotropy)
 206 {
 207     real g = anisotropy;
 208
 209     return INV_FOUR_PI * 1.5 * (1 - g * g) / (2 + g * g);
 210 }
 211
 212 real CornetteShanksPhasePartVarying(real anisotropy, real cosTheta)
 213 {
 214     real g = anisotropy;
 215     real x = 1 + g * g - 2 * g * cosTheta;
 216     real f = rsqrt(max(x, REAL_EPS)); // x^(-1/2)
 217     real h = (1 + cosTheta * cosTheta);
 218
 219     // Note that this function is not perfectly isotropic for (g = 0).
 220     return h * (f * f * f); // h * x^(-3/2)
 221 }
 222
 223 // A better approximation of the Mie phase function.
 224 // Ref: Henyey–Greenstein and Mie phase functions in Monte Carlo radiative transfer computations
 225 real CornetteShanksPhaseFunction(real anisotropy, real cosTheta)
 226 {
 227     return CornetteShanksPhasePartConstant(anisotropy) *
 228            CornetteShanksPhasePartVarying(anisotropy, cosTheta);
 229 }
 230
 231 //
 232 // --------------------------------- Importance Sampling -------------------------------------------
 233 //
 234
 235 // Samples the interval of homogeneous participating medium using the closed-form tracking approach
 236 // (proportionally to the transmittance).
 237 // Returns the offset from the start of the interval and the weight = (transmittance / pdf).
 238 // Ref: Monte Carlo Methods for Volumetric Light Transport Simulation, p. 5.
 239 void ImportanceSampleHomogeneousMedium(real rndVal, real extinction, real intervalLength,
 240                                        out real offset, out real weight)
 241 {
 242     // pdf    = extinction * exp(extinction * (intervalLength - t)) / (exp(intervalLength * extinction) - 1)
 243     // pdf    = extinction * exp(-extinction * t) / (1 - exp(-extinction * intervalLength))
 244     // weight = TransmittanceFromOpticalDepth(t) / pdf
 245     // weight = exp(-extinction * t) / pdf
 246     // weight = (1 - exp(-extinction * intervalLength)) / extinction
 247     // weight = OpacityFromOpticalDepth(extinction * intervalLength) / extinction
 248
 249     real x = 1 - exp(-extinction * intervalLength);
 250     real c = rcp(extinction);
 251
 252     weight = x * c;
 253     offset = -log(1 - rndVal * x) * c;
 254 }
 255
 256 // Implements equiangular light sampling.
 257 // Returns the distance from the origin of the ray, the squared distance from the light,
 258 // and the reciprocal of the PDF.
 259 // Ref: Importance Sampling of Area Lights in Participating Medium.
 260 void ImportanceSamplePunctualLight(real rndVal, real3 lightPosition, real lightSqRadius,
 261                                    real3 rayOrigin, real3 rayDirection,
 262                                    real tMin, real tMax,
 263                                    out real t, out real sqDist, out real rcpPdf)
 264 {
 265     real3 originToLight         = lightPosition - rayOrigin;
 266     real  originToLightProjDist = dot(originToLight, rayDirection);
 267     real  originToLightSqDist   = dot(originToLight, originToLight);
 268     real  rayToLightSqDist      = originToLightSqDist - originToLightProjDist * originToLightProjDist;
 269
 270     // Virtually offset the light to modify the PDF distribution.
 271     real sqD  = max(rayToLightSqDist + lightSqRadius, REAL_EPS);
 272     real rcpD = rsqrt(sqD);
 273     real d    = sqD * rcpD;
 274     real a    = tMin - originToLightProjDist;
 275     real b    = tMax - originToLightProjDist;
 276     real x    = a * rcpD;
 277     real y    = b * rcpD;
 278
 279 #if 0
 280     real theta0   = FastATan(x);
 281     real theta1   = FastATan(y);
 282     real gamma    = theta1 - theta0;
 283     real tanTheta = tan(theta0 + rndVal * gamma);
 284 #else
 285     // Same but faster:
 286     // atan(y) - atan(x) = atan((y - x) / (1 + x * y))
 287     // tan(atan(x) + z)  = (x * cos(z) + sin(z)) / (cos(z) - x * sin(z))
 288     // Both the tangent and the angle  cannot be negative.
 289     real tanGamma = abs((y - x) * rcp(max(0, 1 + x * y)));
 290     real gamma    = FastATanPos(tanGamma);
 291     real z        = rndVal * gamma;
 292     real numer    = x * cos(z) + sin(z);
 293     real denom    = cos(z) - x * sin(z);
 294     real tanTheta = numer * rcp(denom);
 295 #endif
 296
 297     real tRelative = d * tanTheta;
 298
 299     sqDist = sqD + tRelative * tRelative;
 300     rcpPdf = gamma * rcpD * sqDist;
 301     t      = originToLightProjDist + tRelative;
 302
 303     // Remove the virtual light offset to obtain the real geometric distance.
 304     sqDist = max(sqDist - lightSqRadius, REAL_EPS);
 305 }
 306
 307 //
 308 // ------------------------------------ Miscellaneous ----------------------------------------------
 309 //
 310
 311 // Absorption coefficient from Disney: http://blog.selfshadow.com/publications/s2015-shading-course/burley/s2015_pbs_disney_bsdf_notes.pdf
 312 real3 TransmittanceColorAtDistanceToAbsorption(real3 transmittanceColor, real atDistance)
 313 {
 314     return -log(transmittanceColor + REAL_EPS) / max(atDistance, REAL_EPS);
 315 }
 316
 317 #endif // UNITY_VOLUME_RENDERING_INCLUDED