1
0
forked from Rowland/EG
EG/RenderPipelineFile/rpplugins/scattering/shader/eric_bruneton/scattering_common.glsl

344 lines
12 KiB
GLSL

#pragma once
// This file includes defines and functions used for the scattering
// #define SCATTERING_USE_32_BIT 1
#pragma include "render_pipeline_base.inc.glsl"
#extension GL_ARB_shading_language_420pack : enable
// These can be overridden depending on if the shader is a compute shader
// or a regular fragment shader
#ifndef NO_COMPUTE_SHADER
#define PIXEL_X gl_GlobalInvocationID.x
#define PIXEL_Y gl_GlobalInvocationID.y
#else
#define PIXEL_X gl_FragCoord.x
#define PIXEL_Y gl_FragCoord.x
#endif
// This define can be set to 1.0 when writing out the scattering textures, so
// the alpha channel is not 0.0 and the texture cannot get displayed
#define SCAT_DEBUG_ALPHA 0.0
// Actual include:
/**
* Precomputed Atmospheric Scattering
* Copyright (c) 2008 INRIA
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 3. Neither the name of the copyright holders nor the names of its
* contributors may be used to endorse or promote products derived from
* this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
* THE POSSIBILITY OF SUCH DAMAGE.
* Author: Eric Bruneton
*/
// Constants
const float Rg = 6360.0;
const float Rt = 6420.0;
const float RL = 6421.0;
const int TRANSMITTANCE_W = 256 * 4;
const int TRANSMITTANCE_H = 64 * 4;
const int SKY_W = 64 * 4;
const int SKY_H = 16 * 4;
const int RES_R = 32;
const int RES_MU = 128;
const int RES_MU_S = 32;
const int RES_NU = 8;
// ----------------------------------------------------------------------------
// PHYSICAL MODEL PARAMETERS
// ----------------------------------------------------------------------------
const float AVERAGE_GROUND_REFLECTANCE = GET_SETTING(scattering, ground_reflectance);
// Rayleigh
const float HR = GET_SETTING(scattering, rayleigh_height_scale);
const vec3 betaR = vec3(5.8e-3, 1.35e-2, 3.31e-2);
// Mie
const float HM = GET_SETTING(scattering, mie_height_scale);
const vec3 betaMSca = vec3(GET_SETTING(scattering, beta_mie_scattering) * 1e-3);
const float mieG = GET_SETTING(scattering, mie_phase_factor);
// CLEAR SKY
// const float HM = 1.2;
// const vec3 betaMSca = vec3(20e-3);
// const float mieG = 0.76;
// PARTLY CLOUDY
/*const float HM = 3.0;
const vec3 betaMSca = vec3(3e-3);
const float mieG = 0.65;*/
// Beta Mie Extinction
const vec3 betaMEx = betaMSca / 0.9;
// ----------------------------------------------------------------------------
// NUMERICAL INTEGRATION PARAMETERS
// ----------------------------------------------------------------------------
const int TRANSMITTANCE_INTEGRAL_SAMPLES = 500;
const int INSCATTER_INTEGRAL_SAMPLES = 50;
const int IRRADIANCE_INTEGRAL_SAMPLES = 32;
// Needs to be a low value, leads to driver crashes otherwise
const int INSCATTER_SPHERICAL_INTEGRAL_SAMPLES = 8;
// ----------------------------------------------------------------------------
// PARAMETERIZATION OPTIONS
// ----------------------------------------------------------------------------
#define TRANSMITTANCE_NON_LINEAR
#define INSCATTER_NON_LINEAR
// ----------------------------------------------------------------------------
// PARAMETERIZATION FUNCTIONS
// ----------------------------------------------------------------------------
uniform sampler2D transmittanceSampler;
void get_r_dhdh(int layer, out float r, out vec4 dhdh) {
r = float(layer) / float(RES_R - 1.0);
r = r * r;
r = sqrt(Rg * Rg + r * (Rt * Rt - Rg * Rg)) +
(layer == 0 ? 0.01 : (layer == RES_R - 1 ? -0.001 : 0.0));
float dmin = Rt - r;
float dmax = sqrt(r * r - Rg * Rg) + sqrt(Rt * Rt - Rg * Rg);
float dminp = r - Rg;
float dmaxp = sqrt(r * r - Rg * Rg);
dhdh = vec4(dmin, dmax, dminp, dmaxp);
}
vec2 getTransmittanceUV(float r, float mu) {
float uR, uMu;
#ifdef TRANSMITTANCE_NON_LINEAR
uR = sqrt((r - Rg) / (Rt - Rg));
uMu = atan((mu + 0.15) / (1.0 + 0.15) * tan(1.5)) / 1.5;
#else
uR = (r - Rg) / (Rt - Rg);
uMu = (mu + 0.15) / (1.0 + 0.15);
#endif
return vec2(uMu, uR);
}
void getTransmittanceRMu(out float r, out float muS) {
r = PIXEL_Y / float(TRANSMITTANCE_H);
muS = PIXEL_X / float(TRANSMITTANCE_W);
#ifdef TRANSMITTANCE_NON_LINEAR
r = Rg + (r * r) * (Rt - Rg);
muS = -0.15 + tan(1.5 * muS) / tan(1.5) * (1.0 + 0.15);
#else
r = Rg + r * (Rt - Rg);
muS = -0.15 + muS * (1.0 + 0.15);
#endif
}
vec2 getIrradianceUV(float r, float muS) {
float uR = (r - Rg) / (Rt - Rg);
float uMuS = (muS + 0.2) / (1.0 + 0.2);
return vec2(uMuS, uR);
}
void getIrradianceRMuS(out float r, out float muS) {
// POSSIBLE PROBLEMS int(coord) does not match gl_FragCoord ?
r = Rg + int(PIXEL_Y) / (float(SKY_H) - 1.0) * (Rt - Rg);
muS = -0.2 + int(PIXEL_X) / (float(SKY_W) - 1.0) * (1.0 + 0.2);
}
vec4 texture4D(sampler3D table, float r, float mu, float muS, float nu)
{
float Rg_sq = Rg * Rg;
float H = sqrt(Rt * Rt - Rg_sq);
float rho = sqrt(r * r - Rg_sq);
#ifdef INSCATTER_NON_LINEAR
float rmu = r * mu;
float delta = rmu * rmu - r * r + Rg_sq;
vec4 cst = rmu < 0.0 && delta > 0.0 ?
vec4(1.0, 0.0, 0.0, 0.5 - 0.5 / float(RES_MU)) :
vec4(-1.0, H * H, H, 0.5 + 0.5 / float(RES_MU));
float uR = 0.5 / float(RES_R) + rho / H * (1.0 - 1.0 / float(RES_R));
float uMu = cst.w + (rmu * cst.x + sqrt(delta + cst.y)) / (rho + cst.z) *
(0.5 - 1.0 / float(RES_MU));
// paper formula
// float uMuS = 0.5 / float(RES_MU_S) + max((1.0 - exp(-3.0 * muS - 0.6)) / (1.0 - exp(-3.6)), 0.0) * (1.0 - 1.0 / float(RES_MU_S));
// better formula
float uMuS = 0.5 / float(RES_MU_S) + (atan(max(muS, -0.1975) * tan(1.26 * 1.1)) / 1.1 +
(1.0 - 0.26)) * 0.5 * (1.0 - 1.0 / float(RES_MU_S));
#else
float uR = 0.5 / float(RES_R) + rho / H * (1.0 - 1.0 / float(RES_R));
float uMu = 0.5 / float(RES_MU) + (mu + 1.0) / 2.0 * (1.0 - 1.0 / float(RES_MU));
float uMuS = 0.5 / float(RES_MU_S) + max(muS + 0.2, 0.0) / 1.2 * (1.0 - 1.0 / float(RES_MU_S));
#endif
float lerp = (nu + 1.0) / 2.0 * (float(RES_NU) - 1.0);
float uNu = floor(lerp);
lerp = lerp - uNu;
return textureLod(table, vec3((uNu + uMuS) / float(RES_NU), uMu, uR), 0) * (1.0 - lerp) +
textureLod(table, vec3((uNu + uMuS + 1.0) / float(RES_NU), uMu, uR), 0) * lerp;
}
void getMuMuSNu(float r, vec4 dhdH, out float mu, out float muS, out float nu) {
float x = int(PIXEL_X);
float y = int(PIXEL_Y);
#ifdef INSCATTER_NON_LINEAR
if (y < float(RES_MU) / 2.0) {
float d = 1.0 - y / (float(RES_MU) / 2.0 - 1.0);
d = min(max(dhdH.z, d * dhdH.w), dhdH.w * 0.999);
mu = (Rg * Rg - r * r - d * d) / (2.0 * r * d);
mu = min(mu, -sqrt(1.0 - (Rg / r) * (Rg / r)) - 0.001);
} else {
float d = (y - float(RES_MU) / 2.0) / (float(RES_MU) / 2.0 - 1.0);
d = min(max(dhdH.x, d * dhdH.y), dhdH.y * 0.999);
mu = (Rt * Rt - r * r - d * d) / (2.0 * r * d);
}
muS = mod(x, float(RES_MU_S)) / (float(RES_MU_S) - 1.0);
// paper formula
//muS = -(0.6 + log(1.0 - muS * (1.0 - exp(-3.6)))) / 3.0;
// better formula
muS = tan((2.0 * muS - 1.0 + 0.26) * 1.1) / tan(1.26 * 1.1);
nu = -1.0 + floor(x / float(RES_MU_S)) / (float(RES_NU) - 1.0) * 2.0;
#else
mu = -1.0 + 2.0 * y / (float(RES_MU) - 1.0);
muS = mod(x, float(RES_MU_S)) / (float(RES_MU_S) - 1.0);
muS = -0.2 + muS * 1.2;
nu = -1.0 + floor(x / float(RES_MU_S)) / (float(RES_NU) - 1.0) * 2.0;
#endif
}
// ----------------------------------------------------------------------------
// UTILITY FUNCTIONS
// ----------------------------------------------------------------------------
// nearest intersection of ray r,mu with ground or top atmosphere boundary
// mu=cos(ray zenith angle at ray origin)
float limit(float r, float mu) {
float dout = -r * mu + sqrt(r * r * (mu * mu - 1.0) + RL * RL);
float delta2 = r * r * (mu * mu - 1.0) + Rg * Rg;
if (delta2 >= 0.0) {
float din = -r * mu - sqrt(delta2);
if (din >= 0.0) {
dout = min(dout, din);
}
}
return dout;
}
// transmittance(=transparency) of atmosphere for infinite ray (r,mu)
// (mu=cos(view zenith angle)), intersections with ground ignored
vec3 transmittance(float r, float mu) {
vec2 uv = getTransmittanceUV(r, mu);
return textureLod(transmittanceSampler, uv, 0).rgb;
}
// transmittance(=transparency) of atmosphere for infinite ray (r,mu)
// (mu=cos(view zenith angle)), or zero if ray intersects ground
vec3 transmittanceWithShadow(float r, float mu) {
return mu < -sqrt(1.0 - (Rg / r) * (Rg / r)) ? vec3(0.0) : transmittance(r, mu);
}
// transmittance(=transparency) of atmosphere between x and x0
// assume segment x,x0 not intersecting ground
// r=||x||, mu=cos(zenith angle of [x,x0) ray at x), v=unit direction vector of [x,x0) ray
vec3 transmittance(float r, float mu, vec3 v, vec3 x0) {
vec3 result;
float r1 = length(x0);
float mu1 = dot(x0, v) / r;
if (mu > 0.0) {
result = min(transmittance(r, mu) / transmittance(r1, mu1), 1.0);
} else {
result = min(transmittance(r1, -mu1) / transmittance(r, -mu), 1.0);
}
return result;
}
// optical depth for ray (r,mu) of length d, using analytic formula
// (mu=cos(view zenith angle)), intersections with ground ignored
// H=height scale of exponential density function
float opticalDepth(float H, float r, float mu, float d) {
float a = sqrt((0.5 / H) * r);
vec2 a01 = a * vec2(mu, mu + d / r);
vec2 a01s = sign(a01);
vec2 a01sq = a01 * a01;
float x = a01s.y > a01s.x ? exp(a01sq.x) : 0.0;
vec2 y = a01s / (2.3193 * abs(a01) + sqrt(1.52 * a01sq + 4.0)) *
vec2(1.0, exp(-d / H * (d / (2.0 * r) + mu)));
return sqrt((6.2831 * H) * r) * exp((Rg - r) / H) * (x + dot(y, vec2(1.0, -1.0)));
}
// transmittance(=transparency) of atmosphere for ray (r,mu) of length d
// (mu=cos(view zenith angle)), intersections with ground ignored
// uses analytic formula instead of transmittance texture
vec3 analyticTransmittance(float r, float mu, float d) {
return exp(- betaR * opticalDepth(HR, r, mu, d) - betaMEx * opticalDepth(HM, r, mu, d));
}
// transmittance(=transparency) of atmosphere between x and x0
// assume segment x,x0 not intersecting ground
// d = distance between x and x0, mu=cos(zenith angle of [x,x0) ray at x)
vec3 transmittance(float r, float mu, float d) {
vec3 result;
float r1 = sqrt(r * r + d * d + 2.0 * r * mu * d);
float mu1 = (r * mu + d) / r1;
if (mu > 0.0) {
result = min(transmittance(r, mu) / transmittance(r1, mu1), 1.0);
} else {
result = min(transmittance(r1, -mu1) / transmittance(r, -mu), 1.0);
}
return result;
}
vec3 irradiance(sampler2D sampler, float r, float muS) {
vec2 uv = getIrradianceUV(r, muS);
return textureLod(sampler, uv, 0).rgb;
}
// Rayleigh phase function
float phaseFunctionR(float mu) {
return (3.0 / (16.0 * M_PI)) * (1.0 + mu * mu);
}
// Mie phase function
float phaseFunctionM(float mu) {
return 1.5 * 1.0 / (4.0 * M_PI) * (1.0 - mieG * mieG) *
pow(1.0 + (mieG * mieG) - 2.0 * mieG * mu, -3.0 / 2.0) *
(1.0 + mu * mu) / (2.0 + mieG * mieG);
}
// approximated single Mie scattering (cf. approximate Cm in paragraph "Angular precision")
vec3 getMie(vec4 rayMie) { // rayMie.rgb=C*, rayMie.w=Cm,r
return rayMie.rgb * rayMie.w / max(rayMie.r, 1e-4) * (betaR.r / betaR);
}