"\n" "// Returns net particle density to the atmosphere top\n" "//\n" "// | Zenith .\'\n" "// | angle .\'\n" "// | .\'\n" "// | .\'\n" "// | .\'\n" "// _ _ _ _ _|.\'\n" "// A |\n" "// fAltitude | |\n" "// _V_ _ _ _|\n" "// . \' \' \' .\n" "// .\' \'.\n" "// \' \'\n" "float2 GetNetParticleDensity(in float fAltitude, // Altitude (height above the sea level) of the starting point\n" " in float fCosZenithAngle, // Cosine of the zenith angle\n" " in float fAtmBottomAltitude, // Altitude of the bottom atmosphere boundary\n" " in float fAtmAltitudeRangeInv // 1 / (fAtmTopAltitude - fAtmBottomAltitude)\n" " )\n" "{\n" " float fNormalizedAltitude = (fAltitude - fAtmBottomAltitude) * fAtmAltitudeRangeInv;\n" " return g_tex2DOccludedNetDensityToAtmTop.SampleLevel(g_tex2DOccludedNetDensityToAtmTop_sampler, float2(fNormalizedAltitude, fCosZenithAngle*0.5+0.5), 0);\n" "}\n" "\n" "void ApplyPhaseFunctions(inout float3 f3RayleighInscattering,\n" " inout float3 f3MieInscattering,\n" " in float cosTheta)\n" "{\n" " f3RayleighInscattering *= g_MediaParams.f4AngularRayleighSctrCoeff.rgb * (1.0 + cosTheta*cosTheta);\n" " \n" " // Apply Cornette-Shanks phase function (see Nishita et al. 93):\n" " // F(theta) = 1/(4*PI) * 3*(1-g^2) / (2*(2+g^2)) * (1+cos^2(theta)) / (1 + g^2 - 2g*cos(theta))^(3/2)\n" " // f4CS_g = ( 3*(1-g^2) / (2*(2+g^2)), 1+g^2, -2g, 1 )\n" " float fDenom = rsqrt( dot(g_MediaParams.f4CS_g.yz, float2(1.0, cosTheta)) ); // 1 / (1 + g^2 - 2g*cos(theta))^(1/2)\n" " float fCornettePhaseFunc = g_MediaParams.f4CS_g.x * (fDenom*fDenom*fDenom) * (1.0 + cosTheta*cosTheta);\n" " f3MieInscattering *= g_MediaParams.f4AngularMieSctrCoeff.rgb * fCornettePhaseFunc;\n" "}\n" "\n" "// This function computes atmospheric properties in the given point\n" "void GetAtmosphereProperties(in float3 f3Pos,\n" " in float3 f3EarthCentre,\n" " in float fEarthRadius,\n" " in float fAtmBottomAltitude,\n" " in float fAtmAltitudeRangeInv,\n" " in float4 f4ParticleScaleHeight,\n" " in float3 f3DirOnLight,\n" " out float2 f2ParticleDensity,\n" " out float2 f2NetParticleDensityToAtmTop)\n" "{\n" " // Calculate the point height above the SPHERICAL Earth surface:\n" " float3 f3EarthCentreToPointDir = f3Pos - f3EarthCentre;\n" " float fDistToEarthCentre = length(f3EarthCentreToPointDir);\n" " f3EarthCentreToPointDir /= fDistToEarthCentre;\n" " float fHeightAboveSurface = fDistToEarthCentre - fEarthRadius;\n" "\n" " f2ParticleDensity = exp( -fHeightAboveSurface * f4ParticleScaleHeight.zw );\n" "\n" " // Get net particle density from the integration point to the top of the atmosphere:\n" " float fCosSunZenithAngleForCurrPoint = dot( f3EarthCentreToPointDir, f3DirOnLight );\n" " f2NetParticleDensityToAtmTop = GetNetParticleDensity(fHeightAboveSurface, fCosSunZenithAngleForCurrPoint, fAtmBottomAltitude, fAtmAltitudeRangeInv);\n" "}\n" "\n" "// This function computes differential inscattering for the given particle densities \n" "// (without applying phase functions)\n" "void ComputePointDiffInsctr(in float2 f2ParticleDensityInCurrPoint,\n" " in float2 f2NetParticleDensityFromCam,\n" " in float2 f2NetParticleDensityToAtmTop,\n" " out float3 f3DRlghInsctr,\n" " out float3 f3DMieInsctr)\n" "{\n" " // Compute total particle density from the top of the atmosphere through the integraion point to camera\n" " float2 f2TotalParticleDensity = f2NetParticleDensityFromCam + f2NetParticleDensityToAtmTop;\n" " \n" " // Get optical depth\n" " float3 f3TotalRlghOpticalDepth = g_MediaParams.f4RayleighExtinctionCoeff.rgb * f2TotalParticleDensity.x;\n" " float3 f3TotalMieOpticalDepth = g_MediaParams.f4MieExtinctionCoeff.rgb * f2TotalParticleDensity.y;\n" " \n" " // And total extinction for the current integration point:\n" " float3 f3TotalExtinction = exp( -(f3TotalRlghOpticalDepth + f3TotalMieOpticalDepth) );\n" "\n" " f3DRlghInsctr = f2ParticleDensityInCurrPoint.x * f3TotalExtinction;\n" " f3DMieInsctr = f2ParticleDensityInCurrPoint.y * f3TotalExtinction; \n" "}\n" "\n" "void ComputeInsctrIntegral(in float3 f3RayStart,\n" " in float3 f3RayEnd,\n" " in float3 f3EarthCentre,\n" " in float fEarthRadius,\n" " in float fAtmBottomAltitude,\n" " in float fAtmAltitudeRangeInv,\n" " in float4 f4ParticleScaleHeight,\n" " in float3 f3DirOnLight,\n" " in uint uiNumSteps,\n" " inout float2 f2NetParticleDensityFromCam,\n" " inout float3 f3RayleighInscattering,\n" " inout float3 f3MieInscattering)\n" "{\n" " // Evaluate the integrand at the starting point\n" " float2 f2PrevParticleDensity = float2(0.0, 0.0);\n" " float2 f2NetParticleDensityToAtmTop;\n" " GetAtmosphereProperties(f3RayStart,\n" " f3EarthCentre,\n" " fEarthRadius,\n" " fAtmBottomAltitude,\n" " fAtmAltitudeRangeInv,\n" " f4ParticleScaleHeight,\n" " f3DirOnLight,\n" " f2PrevParticleDensity,\n" " f2NetParticleDensityToAtmTop);\n" "\n" " float3 f3PrevDiffRInsctr = float3(0.0, 0.0, 0.0);\n" " float3 f3PrevDiffMInsctr = float3(0.0, 0.0, 0.0);\n" " ComputePointDiffInsctr(f2PrevParticleDensity, f2NetParticleDensityFromCam, f2NetParticleDensityToAtmTop, f3PrevDiffRInsctr, f3PrevDiffMInsctr);\n" "\n" " float fRayLen = length(f3RayEnd - f3RayStart);\n" "\n" " // We want to place more samples when the starting point is close to the surface,\n" " // but for high altitudes linear distribution works better.\n" " float fStartAltitude = length(f3RayStart - f3EarthCentre) - fEarthRadius;\n" " float pwr = lerp(2.0, 1.0, saturate((fStartAltitude - fAtmBottomAltitude) * fAtmAltitudeRangeInv));\n" "\n" " float fPrevSampleDist = 0.0;\n" " for (uint uiSampleNum = 1u; uiSampleNum <= uiNumSteps; ++uiSampleNum)\n" " {\n" " // Evaluate the function at the end of each section and compute the area of a trapezoid\n" " \n" " // We need to place more samples closer to the start point and fewer samples farther away.\n" " // I tried to use more scientific approach (see the bottom of the file), however\n" " // it did not work out because uniform media assumption is inapplicable.\n" " // So instead we will use ad-hoc approach (power function) that works quite well.\n" " float r = pow(float(uiSampleNum) / float(uiNumSteps), pwr);\n" " float3 f3CurrPos = lerp(f3RayStart, f3RayEnd, r);\n" " float fCurrSampleDist = fRayLen * r;\n" " float fStepLen = fCurrSampleDist - fPrevSampleDist;\n" " fPrevSampleDist = fCurrSampleDist;\n" "\n" " float2 f2ParticleDensity;\n" " GetAtmosphereProperties(f3CurrPos,\n" " f3EarthCentre,\n" " fEarthRadius,\n" " fAtmBottomAltitude,\n" " fAtmAltitudeRangeInv,\n" " f4ParticleScaleHeight,\n" " f3DirOnLight,\n" " f2ParticleDensity,\n" " f2NetParticleDensityToAtmTop);\n" "\n" " // Accumulate net particle density from the camera to the integration point:\n" " f2NetParticleDensityFromCam += (f2PrevParticleDensity + f2ParticleDensity) * (fStepLen / 2.0);\n" " f2PrevParticleDensity = f2ParticleDensity;\n" "\n" " float3 f3DRlghInsctr, f3DMieInsctr;\n" " ComputePointDiffInsctr(f2ParticleDensity, f2NetParticleDensityFromCam, f2NetParticleDensityToAtmTop, f3DRlghInsctr, f3DMieInsctr);\n" "\n" " f3RayleighInscattering += (f3DRlghInsctr + f3PrevDiffRInsctr) * (fStepLen / 2.0);\n" " f3MieInscattering += (f3DMieInsctr + f3PrevDiffMInsctr) * (fStepLen / 2.0);\n" "\n" " f3PrevDiffRInsctr = f3DRlghInsctr;\n" " f3PrevDiffMInsctr = f3DMieInsctr;\n" " }\n" "}\n" "\n" "\n" "void IntegrateUnshadowedInscattering(in float3 f3RayStart, \n" " in float3 f3RayEnd,\n" " in float3 f3ViewDir,\n" " in float3 f3EarthCentre,\n" " in float fEarthRadius,\n" " in float fAtmBottomAltitude,\n" " in float fAtmAltitudeRangeInv,\n" " in float4 f4ParticleScaleHeight,\n" " in float3 f3DirOnLight,\n" " in uint uiNumSteps,\n" " out float3 f3Inscattering,\n" " out float3 f3Extinction)\n" "{\n" " float2 f2NetParticleDensityFromCam = float2(0.0, 0.0);\n" " float3 f3RayleighInscattering = float3(0.0, 0.0, 0.0);\n" " float3 f3MieInscattering = float3(0.0, 0.0, 0.0);\n" " ComputeInsctrIntegral( f3RayStart,\n" " f3RayEnd,\n" " f3EarthCentre,\n" " fEarthRadius,\n" " fAtmBottomAltitude,\n" " fAtmAltitudeRangeInv,\n" " f4ParticleScaleHeight,\n" " f3DirOnLight,\n" " uiNumSteps,\n" " f2NetParticleDensityFromCam,\n" " f3RayleighInscattering,\n" " f3MieInscattering);\n" "\n" " float3 f3TotalRlghOpticalDepth = g_MediaParams.f4RayleighExtinctionCoeff.rgb * f2NetParticleDensityFromCam.x;\n" " float3 f3TotalMieOpticalDepth = g_MediaParams.f4MieExtinctionCoeff.rgb * f2NetParticleDensityFromCam.y;\n" " f3Extinction = exp( -(f3TotalRlghOpticalDepth + f3TotalMieOpticalDepth) );\n" "\n" " // Apply phase function\n" " // Note that cosTheta = dot(DirOnCamera, LightDir) = dot(ViewDir, DirOnLight) because\n" " // DirOnCamera = -ViewDir and LightDir = -DirOnLight\n" " float cosTheta = dot(f3ViewDir, f3DirOnLight);\n" " ApplyPhaseFunctions(f3RayleighInscattering, f3MieInscattering, cosTheta);\n" "\n" " f3Inscattering = f3RayleighInscattering + f3MieInscattering;\n" "}\n" "\n" "\n" "// BACKUP\n" "//\n" "// To better distribute samples for numerical integration I tried used the following simplifying assumptions:\n" "// - Uniform media density\n" "// - No light extinction\n" "// Under these assumptions the total amount of light inscattered from 0 to distance t will be computed by\n" "//\n" "// I(t) = Integral(exp(-a*x), 0 -> t) = 1/a * (1 - exp(-a*t))\n" "//\n" "// We want to find sampling that produces even distribution for function I(t):\n" "//\n" "// I(tn) = I(D) * n / N, where D is the total distance and N is the total sample count.\n" "//\n" "// From here:\n" "//\n" "// 1/a * (1 - exp(-a*tn)) = I(D) * n / N\n" "// 1 - exp(-a*tn) = a * I(D) * n / N\n" "// exp(-a*tn) = 1 - a * I(D) * n / N\n" "// tn = -1/a * log(1 - a * I(D) * n / N) = -1/a * log(1 - n/N * (1 - exp(-a*D)))\n" "//\n" "// This unfortunately did not work because homogeneous media assumption is inappropriate:\n" "// - when a = 1e-5, the distribution is basically linear\n" "// - when a = 1e-4, the first N-1 samples cover first 50% distance, while the last one covers the remaining 50% distance\n"