From e9146f00a320568bdd922988387bd25e7fd732e6 Mon Sep 17 00:00:00 2001 From: Dario Seyb Date: Tue, 23 Feb 2016 22:35:22 +0100 Subject: [PATCH] moved atmosphere precomputation to the GPU --- README.txt | 13 + data/shader/Atmosphere.fsh | 9 +- data/shader/atmosphere/common.glsl | 288 +++++++++ data/shader/atmosphere/copyInscatter1.glsl | 88 +++ data/shader/atmosphere/copyInscatterN.glsl | 87 +++ data/shader/atmosphere/copyIrradiance.glsl | 64 ++ data/shader/atmosphere/inscatter1.glsl | 124 ++++ data/shader/atmosphere/inscatterN.glsl | 109 ++++ data/shader/atmosphere/inscatterS.glsl | 167 +++++ data/shader/atmosphere/irradiance1.glsl | 62 ++ data/shader/atmosphere/irradianceN.glsl | 93 +++ data/shader/atmosphere/transmittance.glsl | 75 +++ data/shader/planets/earth/earth_texture.glsl | 16 +- .../terrain/dispMappingWithWater_fsh.glsl | 4 +- .../include/ACGL/OpenGL/Objects/Texture.hh | 23 +- .../engine/scene/OrbitalSimulationSystem.cpp | 594 +++++++++--------- 16 files changed, 1499 insertions(+), 317 deletions(-) create mode 100644 README.txt create mode 100644 data/shader/atmosphere/common.glsl create mode 100644 data/shader/atmosphere/copyInscatter1.glsl create mode 100644 data/shader/atmosphere/copyInscatterN.glsl create mode 100644 data/shader/atmosphere/copyIrradiance.glsl create mode 100644 data/shader/atmosphere/inscatter1.glsl create mode 100644 data/shader/atmosphere/inscatterN.glsl create mode 100644 data/shader/atmosphere/inscatterS.glsl create mode 100644 data/shader/atmosphere/irradiance1.glsl create mode 100644 data/shader/atmosphere/irradianceN.glsl create mode 100644 data/shader/atmosphere/transmittance.glsl diff --git a/README.txt b/README.txt new file mode 100644 index 0000000..134a72c --- /dev/null +++ b/README.txt @@ -0,0 +1,13 @@ +This is a simplified version of the implementation of our "Precomputed Atmospheric Scattering" article (EGSR 2008). This version uses a perfectly spherical terrain for rendering, so there are no light shafts and the code for them is not included. Also the aerial perspective rendering is simplified, since only a single lookup in the precomputed S table is necessary (instead of two for an irregular terrain). For these reasons the rendering performance is higher than what is described in the paper. The code to precompute our 2D and 4D tables is fully included, without any simplification. + +User Interface +- mouse: + drag = move the Sun + SHIFT + drag = move the camera (longitude, latitude) + CTRL + drag = move the camera (azimut, site) +- page up / page down : increase / decrease distance to ground +- + / - : increase / decrease exposure parameter for tone mapping +- F5 : reload all the shaders and recompute the precomputed tables +- ESC : exit + +you can use F5 to see the effect of changing the atmospheric parameters defined in common.glsl diff --git a/data/shader/Atmosphere.fsh b/data/shader/Atmosphere.fsh index 2f6f1f8..a46f2ef 100644 --- a/data/shader/Atmosphere.fsh +++ b/data/shader/Atmosphere.fsh @@ -9,7 +9,7 @@ uniform sampler2D uSamplerColor; float Rg; float Rt; -const float ISun = 3.0; +const float ISun = 4.0; const vec3 betaR = vec3(5.8e-3, 1.35e-2, 3.31e-2); const float mieG = 0.8; @@ -192,7 +192,7 @@ vec4 color() { vec3 viewDir = normalize(vPosition - uCameraPosition); float planetRadius = uEmissiveColor.r; //6360; - float atmosphereRadius = uEmissiveColor.g * 1.0115; //6778; + float atmosphereRadius = uEmissiveColor.g; //6778; vec3 camSpaceWorldPos = depth != 1 ? unpackWorldPosition(depth) - uCameraPosition : @@ -248,9 +248,6 @@ vec4 color() { //return vec4(t/planetRadius); - - //t = -r * mu - sqrt(r * r * (mu * mu - 1.0) + planetRadius * planetRadius); - vec3 g = x - vec3(0.0, 0.0, planetRadius + 10.0); float a = v.x * v.x + v.y * v.y - v.z * v.z; float b = 2.0 * (g.x * v.x + g.y * v.y - g.z * v.z); @@ -271,7 +268,7 @@ vec4 color() { fogColor = max(vec3(0.0), inscatter(x, worldPos, t, v, s, r, mu, attenuation)); vec3 backgroundColor = texture(uSamplerColor, screenTexCoord).rgb; - fogColor += max(vec3(0.0), attenuation * backgroundColor); + fogColor += max(vec3(0.0), backgroundColor); fogAmount = 1.0; diff --git a/data/shader/atmosphere/common.glsl b/data/shader/atmosphere/common.glsl new file mode 100644 index 0000000..eb07331 --- /dev/null +++ b/data/shader/atmosphere/common.glsl @@ -0,0 +1,288 @@ +/** + * 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 + */ + +// ---------------------------------------------------------------------------- +// PHYSICAL MODEL PARAMETERS +// ---------------------------------------------------------------------------- +const float Rg = 6360.0; +const float Rt = 6420.0; +const float RL = 6421.0; + +const int TRANSMITTANCE_W = 256; +const int TRANSMITTANCE_H = 64; + +const int SKY_W = 64; +const int SKY_H = 16; + +const int RES_R = 32; +const int RES_MU = 128; +const int RES_MU_S = 32; +const int RES_NU = 8; + +const float AVERAGE_GROUND_REFLECTANCE = 0.1; + +// Rayleigh +const float HR = 8.0; +const vec3 betaR = vec3(5.8e-3, 1.35e-2, 3.31e-2); + +// Mie +// DEFAULT +const float HM = 1.2; +const vec3 betaMSca = vec3(4e-3); +const vec3 betaMEx = betaMSca / 0.9; +const float mieG = 0.8; +// CLEAR SKY +/*const float HM = 1.2; +const vec3 betaMSca = vec3(20e-3); +const vec3 betaMEx = betaMSca / 0.9; +const float mieG = 0.76;*/ +// PARTLY CLOUDY +/*const float HM = 3.0; +const vec3 betaMSca = vec3(3e-3); +const vec3 betaMEx = betaMSca / 0.9; +const float mieG = 0.65;*/ + +// ---------------------------------------------------------------------------- +// NUMERICAL INTEGRATION PARAMETERS +// ---------------------------------------------------------------------------- + +const int TRANSMITTANCE_INTEGRAL_SAMPLES = 500; +const int INSCATTER_INTEGRAL_SAMPLES = 50; +const int IRRADIANCE_INTEGRAL_SAMPLES = 32; +const int INSCATTER_SPHERICAL_INTEGRAL_SAMPLES = 16; + +const float M_PI = 3.141592657; + +// ---------------------------------------------------------------------------- +// PARAMETERIZATION OPTIONS +// ---------------------------------------------------------------------------- + +#define TRANSMITTANCE_NON_LINEAR +#define INSCATTER_NON_LINEAR + +// ---------------------------------------------------------------------------- +// PARAMETERIZATION FUNCTIONS +// ---------------------------------------------------------------------------- + +uniform sampler2D transmittanceSampler; + +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 = gl_FragCoord.y / float(TRANSMITTANCE_H); + muS = gl_FragCoord.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) { + r = Rg + (gl_FragCoord.y - 0.5) / (float(SKY_H) - 1.0) * (Rt - Rg); + muS = -0.2 + (gl_FragCoord.x - 0.5) / (float(SKY_W) - 1.0) * (1.0 + 0.2); +} + +vec4 texture4D(sampler3D table, float r, float mu, float muS, float nu) +{ + float H = sqrt(Rt * Rt - Rg * Rg); + float rho = sqrt(r * r - Rg * Rg); +#ifdef INSCATTER_NON_LINEAR + float rmu = r * mu; + float delta = rmu * rmu - r * r + Rg * Rg; + 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 texture3D(table, vec3((uNu + uMuS) / float(RES_NU), uMu, uR)) * (1.0 - lerp) + + texture3D(table, vec3((uNu + uMuS + 1.0) / float(RES_NU), uMu, uR)) * lerp; +} + +void getMuMuSNu(float r, vec4 dhdH, out float mu, out float muS, out float nu) { + float x = gl_FragCoord.x - 0.5; + float y = gl_FragCoord.y - 0.5; +#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 texture2D(transmittanceSampler, uv).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 texture2D(sampler, uv).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); +} diff --git a/data/shader/atmosphere/copyInscatter1.glsl b/data/shader/atmosphere/copyInscatter1.glsl new file mode 100644 index 0000000..b376598 --- /dev/null +++ b/data/shader/atmosphere/copyInscatter1.glsl @@ -0,0 +1,88 @@ +/** + * 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 + */ +// copies deltaS into S (line 5 in algorithm 4.1) + +uniform float r; +uniform vec4 dhdH; +uniform int layer; + +uniform sampler3D deltaSRSampler; +uniform sampler3D deltaSMSampler; + +#ifdef _VERTEX_ + + +const vec2 quad[4] = vec2[] ( + vec2(-1.0, 1.0), + vec2(-1.0,-1.0), + vec2( 1.0, 1.0), + vec2( 1.0,-1.0) +); + +void main() +{ + vec2 p = quad[ gl_VertexID ]; + gl_Position = vec4(p, 0.0, 1.0); +} + +#endif + +#ifdef _GEOMETRY_ +#extension GL_EXT_geometry_shader4 : enable + +void main() { + gl_Position = gl_PositionIn[0]; + gl_Layer = layer; + EmitVertex(); + gl_Position = gl_PositionIn[1]; + gl_Layer = layer; + EmitVertex(); + gl_Position = gl_PositionIn[2]; + gl_Layer = layer; + EmitVertex(); + EndPrimitive(); +} + +#endif + +#ifdef _FRAGMENT_ + +void main() { + vec3 uvw = vec3(gl_FragCoord.xy, float(layer) + 0.5) / vec3(ivec3(RES_MU_S * RES_NU, RES_MU, RES_R)); + vec4 ray = texture3D(deltaSRSampler, uvw); + vec4 mie = texture3D(deltaSMSampler, uvw); + gl_FragColor = vec4(ray.rgb, mie.r); // store only red component of single Mie scattering (cf. "Angular precision") +} + +#endif diff --git a/data/shader/atmosphere/copyInscatterN.glsl b/data/shader/atmosphere/copyInscatterN.glsl new file mode 100644 index 0000000..62cec53 --- /dev/null +++ b/data/shader/atmosphere/copyInscatterN.glsl @@ -0,0 +1,87 @@ +/** + * 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 + */ + +// adds deltaS into S (line 11 in algorithm 4.1) +uniform float r; +uniform vec4 dhdH; +uniform int layer; + +uniform sampler3D deltaSSampler; + +#ifdef _VERTEX_ + + +const vec2 quad[4] = vec2[] ( + vec2(-1.0, 1.0), + vec2(-1.0,-1.0), + vec2( 1.0, 1.0), + vec2( 1.0,-1.0) +); + +void main() +{ + vec2 p = quad[ gl_VertexID ]; + gl_Position = vec4(p, 0.0, 1.0); +} + +#endif + +#ifdef _GEOMETRY_ +#extension GL_EXT_geometry_shader4 : enable + +void main() { + gl_Position = gl_PositionIn[0]; + gl_Layer = layer; + EmitVertex(); + gl_Position = gl_PositionIn[1]; + gl_Layer = layer; + EmitVertex(); + gl_Position = gl_PositionIn[2]; + gl_Layer = layer; + EmitVertex(); + EndPrimitive(); +} + +#endif + +#ifdef _FRAGMENT_ + +void main() { + float mu, muS, nu; + getMuMuSNu(r, dhdH, mu, muS, nu); + vec3 uvw = vec3(gl_FragCoord.xy, float(layer) + 0.5) / vec3(ivec3(RES_MU_S * RES_NU, RES_MU, RES_R)); + gl_FragColor = vec4(texture3D(deltaSSampler, uvw).rgb / phaseFunctionR(nu), 0.0); +} + +#endif diff --git a/data/shader/atmosphere/copyIrradiance.glsl b/data/shader/atmosphere/copyIrradiance.glsl new file mode 100644 index 0000000..eb8bfd9 --- /dev/null +++ b/data/shader/atmosphere/copyIrradiance.glsl @@ -0,0 +1,64 @@ +/** + * 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 + */ + +// copies deltaE into E (lines 4 and 10 in algorithm 4.1) +uniform float k; // k=0 for line 4, k=1 for line 10 +uniform sampler2D deltaESampler; + +#ifdef _VERTEX_ + + +const vec2 quad[4] = vec2[] ( + vec2(-1.0, 1.0), + vec2(-1.0,-1.0), + vec2( 1.0, 1.0), + vec2( 1.0,-1.0) +); + +void main() +{ + vec2 p = quad[ gl_VertexID ]; + gl_Position = vec4(p, 0.0, 1.0); +} + +#endif + +#ifdef _FRAGMENT_ + +void main() { + vec2 uv = gl_FragCoord.xy / vec2(SKY_W, SKY_H); + gl_FragColor = k * texture2D(deltaESampler, uv); // k=0 for line 4, k=1 for line 10 +} + +#endif diff --git a/data/shader/atmosphere/inscatter1.glsl b/data/shader/atmosphere/inscatter1.glsl new file mode 100644 index 0000000..4990789 --- /dev/null +++ b/data/shader/atmosphere/inscatter1.glsl @@ -0,0 +1,124 @@ +/** + * 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 + */ + +// computes single scattering (line 3 in algorithm 4.1) +uniform float r; +uniform vec4 dhdH; +uniform int layer; + +#ifdef _VERTEX_ +const vec2 quad[4] = vec2[] ( + vec2(-1.0, 1.0), + vec2(-1.0,-1.0), + vec2( 1.0, 1.0), + vec2( 1.0,-1.0) +); + +void main() +{ + vec2 p = quad[ gl_VertexID ]; + gl_Position = vec4(p, 0.0, 1.0); +} + +#endif + +#ifdef _GEOMETRY_ +#extension GL_EXT_geometry_shader4 : enable + +void main() { + gl_Position = gl_PositionIn[0]; + gl_Layer = layer; + EmitVertex(); + gl_Position = gl_PositionIn[1]; + gl_Layer = layer; + EmitVertex(); + gl_Position = gl_PositionIn[2]; + gl_Layer = layer; + EmitVertex(); + EndPrimitive(); +} + +#endif + +#ifdef _FRAGMENT_ + +void integrand(float r, float mu, float muS, float nu, float t, out vec3 ray, out vec3 mie) { + ray = vec3(0.0); + mie = vec3(0.0); + float ri = sqrt(r * r + t * t + 2.0 * r * mu * t); + float muSi = (nu * t + muS * r) / ri; + ri = max(Rg, ri); + if (muSi >= -sqrt(1.0 - Rg * Rg / (ri * ri))) { + vec3 ti = transmittance(r, mu, t) * transmittance(ri, muSi); + ray = exp(-(ri - Rg) / HR) * ti; + mie = exp(-(ri - Rg) / HM) * ti; + } +} + +void inscatter(float r, float mu, float muS, float nu, out vec3 ray, out vec3 mie) { + ray = vec3(0.0); + mie = vec3(0.0); + float dx = limit(r, mu) / float(INSCATTER_INTEGRAL_SAMPLES); + float xi = 0.0; + vec3 rayi; + vec3 miei; + integrand(r, mu, muS, nu, 0.0, rayi, miei); + for (int i = 1; i <= INSCATTER_INTEGRAL_SAMPLES; ++i) { + float xj = float(i) * dx; + vec3 rayj; + vec3 miej; + integrand(r, mu, muS, nu, xj, rayj, miej); + ray += (rayi + rayj) / 2.0 * dx; + mie += (miei + miej) / 2.0 * dx; + xi = xj; + rayi = rayj; + miei = miej; + } + ray *= betaR; + mie *= betaMSca; +} + +void main() { + vec3 ray; + vec3 mie; + float mu, muS, nu; + getMuMuSNu(r, dhdH, mu, muS, nu); + inscatter(r, mu, muS, nu, ray, mie); + // store separately Rayleigh and Mie contributions, WITHOUT the phase function factor + // (cf "Angular precision") + gl_FragData[0].rgb = ray; + gl_FragData[1].rgb = mie; +} + +#endif diff --git a/data/shader/atmosphere/inscatterN.glsl b/data/shader/atmosphere/inscatterN.glsl new file mode 100644 index 0000000..69c90f2 --- /dev/null +++ b/data/shader/atmosphere/inscatterN.glsl @@ -0,0 +1,109 @@ +/** + * 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 + */ + +// computes higher order scattering (line 9 in algorithm 4.1) + +uniform float r; +uniform vec4 dhdH; +uniform int layer; + +uniform sampler3D deltaJSampler; + +#ifdef _VERTEX_ + +const vec2 quad[4] = vec2[] ( + vec2(-1.0, 1.0), + vec2(-1.0,-1.0), + vec2( 1.0, 1.0), + vec2( 1.0,-1.0) +); + +void main() +{ + vec2 p = quad[ gl_VertexID ]; + gl_Position = vec4(p, 0.0, 1.0); +} + +#endif + +#ifdef _GEOMETRY_ +#extension GL_EXT_geometry_shader4 : enable + +void main() { + gl_Position = gl_PositionIn[0]; + gl_Layer = layer; + EmitVertex(); + gl_Position = gl_PositionIn[1]; + gl_Layer = layer; + EmitVertex(); + gl_Position = gl_PositionIn[2]; + gl_Layer = layer; + EmitVertex(); + EndPrimitive(); +} + +#endif + +#ifdef _FRAGMENT_ + +vec3 integrand(float r, float mu, float muS, float nu, float t) { + float ri = sqrt(r * r + t * t + 2.0 * r * mu * t); + float mui = (r * mu + t) / ri; + float muSi = (nu * t + muS * r) / ri; + return texture4D(deltaJSampler, ri, mui, muSi, nu).rgb * transmittance(r, mu, t); +} + +vec3 inscatter(float r, float mu, float muS, float nu) { + vec3 raymie = vec3(0.0); + float dx = limit(r, mu) / float(INSCATTER_INTEGRAL_SAMPLES); + float xi = 0.0; + vec3 raymiei = integrand(r, mu, muS, nu, 0.0); + for (int i = 1; i <= INSCATTER_INTEGRAL_SAMPLES; ++i) { + float xj = float(i) * dx; + vec3 raymiej = integrand(r, mu, muS, nu, xj); + raymie += (raymiei + raymiej) / 2.0 * dx; + xi = xj; + raymiei = raymiej; + } + return raymie; +} + +void main() { + float mu, muS, nu; + getMuMuSNu(r, dhdH, mu, muS, nu); + gl_FragColor.rgb = inscatter(r, mu, muS, nu); +} + +#endif + diff --git a/data/shader/atmosphere/inscatterS.glsl b/data/shader/atmosphere/inscatterS.glsl new file mode 100644 index 0000000..87c994b --- /dev/null +++ b/data/shader/atmosphere/inscatterS.glsl @@ -0,0 +1,167 @@ +/** + * 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 + */ + +// computes deltaJ (line 7 in algorithm 4.1) +uniform float r; +uniform vec4 dhdH; +uniform int layer; + +uniform sampler2D deltaESampler; +uniform sampler3D deltaSRSampler; +uniform sampler3D deltaSMSampler; +uniform float first; + +#ifdef _VERTEX_ + + +const vec2 quad[4] = vec2[] ( + vec2(-1.0, 1.0), + vec2(-1.0,-1.0), + vec2( 1.0, 1.0), + vec2( 1.0,-1.0) +); + +void main() +{ + vec2 p = quad[ gl_VertexID ]; + gl_Position = vec4(p, 0.0, 1.0); +} + +#endif + +#ifdef _GEOMETRY_ +#extension GL_EXT_geometry_shader4 : enable + +void main() { + gl_Position = gl_PositionIn[0]; + gl_Layer = layer; + EmitVertex(); + gl_Position = gl_PositionIn[1]; + gl_Layer = layer; + EmitVertex(); + gl_Position = gl_PositionIn[2]; + gl_Layer = layer; + EmitVertex(); + EndPrimitive(); +} + +#endif + +#ifdef _FRAGMENT_ + +const float dphi = M_PI / float(INSCATTER_SPHERICAL_INTEGRAL_SAMPLES); +const float dtheta = M_PI / float(INSCATTER_SPHERICAL_INTEGRAL_SAMPLES); + +void inscatter(float r, float mu, float muS, float nu, out vec3 raymie) { + r = clamp(r, Rg, Rt); + mu = clamp(mu, -1.0, 1.0); + muS = clamp(muS, -1.0, 1.0); + float var = sqrt(1.0 - mu * mu) * sqrt(1.0 - muS * muS); + nu = clamp(nu, muS * mu - var, muS * mu + var); + + float cthetamin = -sqrt(1.0 - (Rg / r) * (Rg / r)); + + vec3 v = vec3(sqrt(1.0 - mu * mu), 0.0, mu); + float sx = v.x == 0.0 ? 0.0 : (nu - muS * mu) / v.x; + vec3 s = vec3(sx, sqrt(max(0.0, 1.0 - sx * sx - muS * muS)), muS); + + raymie = vec3(0.0); + + // integral over 4.PI around x with two nested loops over w directions (theta,phi) -- Eq (7) + for (int itheta = 0; itheta < INSCATTER_SPHERICAL_INTEGRAL_SAMPLES; ++itheta) { + float theta = (float(itheta) + 0.5) * dtheta; + float ctheta = cos(theta); + + float greflectance = 0.0; + float dground = 0.0; + vec3 gtransp = vec3(0.0); + if (ctheta < cthetamin) { // if ground visible in direction w + // compute transparency gtransp between x and ground + greflectance = AVERAGE_GROUND_REFLECTANCE / M_PI; + dground = -r * ctheta - sqrt(r * r * (ctheta * ctheta - 1.0) + Rg * Rg); + gtransp = transmittance(Rg, -(r * ctheta + dground) / Rg, dground); + } + + for (int iphi = 0; iphi < 2 * INSCATTER_SPHERICAL_INTEGRAL_SAMPLES; ++iphi) { + float phi = (float(iphi) + 0.5) * dphi; + float dw = dtheta * dphi * sin(theta); + vec3 w = vec3(cos(phi) * sin(theta), sin(phi) * sin(theta), ctheta); + + float nu1 = dot(s, w); + float nu2 = dot(v, w); + float pr2 = phaseFunctionR(nu2); + float pm2 = phaseFunctionM(nu2); + + // compute irradiance received at ground in direction w (if ground visible) =deltaE + vec3 gnormal = (vec3(0.0, 0.0, r) + dground * w) / Rg; + vec3 girradiance = irradiance(deltaESampler, Rg, dot(gnormal, s)); + + vec3 raymie1; // light arriving at x from direction w + + // first term = light reflected from the ground and attenuated before reaching x, =T.alpha/PI.deltaE + raymie1 = greflectance * girradiance * gtransp; + + // second term = inscattered light, =deltaS + if (first == 1.0) { + // first iteration is special because Rayleigh and Mie were stored separately, + // without the phase functions factors; they must be reintroduced here + float pr1 = phaseFunctionR(nu1); + float pm1 = phaseFunctionM(nu1); + vec3 ray1 = texture4D(deltaSRSampler, r, w.z, muS, nu1).rgb; + vec3 mie1 = texture4D(deltaSMSampler, r, w.z, muS, nu1).rgb; + raymie1 += ray1 * pr1 + mie1 * pm1; + } else { + raymie1 += texture4D(deltaSRSampler, r, w.z, muS, nu1).rgb; + } + + // light coming from direction w and scattered in direction v + // = light arriving at x from direction w (raymie1) * SUM(scattering coefficient * phaseFunction) + // see Eq (7) + raymie += raymie1 * (betaR * exp(-(r - Rg) / HR) * pr2 + betaMSca * exp(-(r - Rg) / HM) * pm2) * dw; + } + } + + // output raymie = J[T.alpha/PI.deltaE + deltaS] (line 7 in algorithm 4.1) +} + +void main() { + vec3 raymie; + float mu, muS, nu; + getMuMuSNu(r, dhdH, mu, muS, nu); + inscatter(r, mu, muS, nu, raymie); + gl_FragColor.rgb = raymie; +} + +#endif + diff --git a/data/shader/atmosphere/irradiance1.glsl b/data/shader/atmosphere/irradiance1.glsl new file mode 100644 index 0000000..55ee5f8 --- /dev/null +++ b/data/shader/atmosphere/irradiance1.glsl @@ -0,0 +1,62 @@ +/** + * 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 + */ + +// computes ground irradiance due to direct sunlight E[L0] (line 2 in algorithm 4.1) +#ifdef _VERTEX_ + + +const vec2 quad[4] = vec2[] ( + vec2(-1.0, 1.0), + vec2(-1.0,-1.0), + vec2( 1.0, 1.0), + vec2( 1.0,-1.0) +); + +void main() +{ + vec2 p = quad[ gl_VertexID ]; + gl_Position = vec4(p, 0.0, 1.0); +} + +#endif + +#ifdef _FRAGMENT_ + +void main() { + float r, muS; + getIrradianceRMuS(r, muS); + gl_FragColor = vec4(transmittance(r, muS) * max(muS, 0.0), 0.0); +} + +#endif diff --git a/data/shader/atmosphere/irradianceN.glsl b/data/shader/atmosphere/irradianceN.glsl new file mode 100644 index 0000000..f092b98 --- /dev/null +++ b/data/shader/atmosphere/irradianceN.glsl @@ -0,0 +1,93 @@ +/** + * 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 + */ +// computes ground irradiance due to skylight E[deltaS] (line 8 in algorithm 4.1) + +uniform sampler3D deltaSRSampler; +uniform sampler3D deltaSMSampler; +uniform float first; + +#ifdef _VERTEX_ + +const vec2 quad[4] = vec2[] ( + vec2(-1.0, 1.0), + vec2(-1.0,-1.0), + vec2( 1.0, 1.0), + vec2( 1.0,-1.0) +); + +void main() +{ + vec2 p = quad[ gl_VertexID ]; + gl_Position = vec4(p, 0.0, 1.0); +} + +#endif + +#ifdef _FRAGMENT_ + +const float dphi = M_PI / float(IRRADIANCE_INTEGRAL_SAMPLES); +const float dtheta = M_PI / float(IRRADIANCE_INTEGRAL_SAMPLES); + +void main() { + float r, muS; + getIrradianceRMuS(r, muS); + vec3 s = vec3(max(sqrt(1.0 - muS * muS), 0.0), 0.0, muS); + + vec3 result = vec3(0.0); + // integral over 2.PI around x with two nested loops over w directions (theta,phi) -- Eq (15) + for (int iphi = 0; iphi < 2 * IRRADIANCE_INTEGRAL_SAMPLES; ++iphi) { + float phi = (float(iphi) + 0.5) * dphi; + for (int itheta = 0; itheta < IRRADIANCE_INTEGRAL_SAMPLES / 2; ++itheta) { + float theta = (float(itheta) + 0.5) * dtheta; + float dw = dtheta * dphi * sin(theta); + vec3 w = vec3(cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta)); + float nu = dot(s, w); + if (first == 1.0) { + // first iteration is special because Rayleigh and Mie were stored separately, + // without the phase functions factors; they must be reintroduced here + float pr1 = phaseFunctionR(nu); + float pm1 = phaseFunctionM(nu); + vec3 ray1 = texture4D(deltaSRSampler, r, w.z, muS, nu).rgb; + vec3 mie1 = texture4D(deltaSMSampler, r, w.z, muS, nu).rgb; + result += (ray1 * pr1 + mie1 * pm1) * w.z * dw; + } else { + result += texture4D(deltaSRSampler, r, w.z, muS, nu).rgb * w.z * dw; + } + } + } + + gl_FragColor = vec4(result, 0.0); +} + +#endif diff --git a/data/shader/atmosphere/transmittance.glsl b/data/shader/atmosphere/transmittance.glsl new file mode 100644 index 0000000..a6f46c3 --- /dev/null +++ b/data/shader/atmosphere/transmittance.glsl @@ -0,0 +1,75 @@ +/** + * 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 + */ + +// computes transmittance table T using Eq (5) + +#ifdef _VERTEX_ + +const vec2 quad[4] = vec2[] ( + vec2(-1.0, 1.0), + vec2(-1.0,-1.0), + vec2( 1.0, 1.0), + vec2( 1.0,-1.0) +); + +void main() +{ + vec2 p = quad[ gl_VertexID ]; + gl_Position = vec4(p, 0.0, 1.0); +} +#else + +float opticalDepth(float H, float r, float mu) { + float result = 0.0; + float dx = limit(r, mu) / float(TRANSMITTANCE_INTEGRAL_SAMPLES); + float xi = 0.0; + float yi = exp(-(r - Rg) / H); + for (int i = 1; i <= TRANSMITTANCE_INTEGRAL_SAMPLES; ++i) { + float xj = float(i) * dx; + float yj = exp(-(sqrt(r * r + xj * xj + 2.0 * xj * r * mu) - Rg) / H); + result += (yi + yj) / 2.0 * dx; + xi = xj; + yi = yj; + } + return mu < -sqrt(1.0 - (Rg / r) * (Rg / r)) ? 1e9 : result; +} + +void main() { + float r, muS; + getTransmittanceRMu(r, muS); + vec3 depth = betaR * opticalDepth(HR, r, muS) + betaMEx * opticalDepth(HM, r, muS); + gl_FragColor = vec4(exp(-depth), 0.0); // Eq (5) +} + +#endif diff --git a/data/shader/planets/earth/earth_texture.glsl b/data/shader/planets/earth/earth_texture.glsl index 546dc20..6377235 100644 --- a/data/shader/planets/earth/earth_texture.glsl +++ b/data/shader/planets/earth/earth_texture.glsl @@ -38,13 +38,11 @@ vec4 texture_getColor(vec3 texCoord3d, vec3 upVectorNormalized, vec3 normalNorma float noiseOct4 = 0.; float noiseOct5 = 0.; - // if (depth == 1.) { - noiseOct1 = snoise_3d(32.* texCoord3d); - noiseOct2 = snoise_3d(64.* texCoord3d); - noiseOct3 = snoise_3d(128. * texCoord3d); - //noiseOct4 = snoise_3d(2*128.*128. * texCoord3d); - //noiseOct5 = snoise_3d(4*128.*128. * texCoord3d); - // } + noiseOct1 = snoise_3d(32.*128. * texCoord3d); + noiseOct2 = snoise_3d(64.*128. * texCoord3d); + noiseOct3 = snoise_3d(128.*128. * texCoord3d); + + float heightPertubation = 1./2.*noiseOct1 + 1./4.*noiseOct2 + 1./8.*noiseOct3; // just some test: // if (noiseOct1 < 0. || noiseOct2 < 0. || noiseOct3 < 0. || noiseOct4 < 0. || noiseOct5 < 0. ) { @@ -63,7 +61,7 @@ vec4 texture_getColor(vec3 texCoord3d, vec3 upVectorNormalized, vec3 normalNorma // + getDetailHeightNoise(texCoord3d) // )*.9 + heightPertubation * .1; - float height =getHeightNoise(texCoord3d)/2.*.9; + float height =getHeightNoise(texCoord3d)/2.*.9 + heightPertubation * 0.05; //testing stuff: //return vec4(normalNormalized.x,normalNormalized.y, normalNormalized.z, 1.); @@ -79,7 +77,7 @@ vec4 texture_getColor(vec3 texCoord3d, vec3 upVectorNormalized, vec3 normalNorma if (height> SNOW_LEVEL && phiPert <= ROCK_SLOPE) { - color = vec4(1., 1., 1., 1.); + color = vec4(1., 1., 1., 1.) * 0.9; } else if (height> ROCK_LEVEL || phiPert > ROCK_SLOPE) { color = vec4(.1, .1, .1, 1.); diff --git a/data/shader/terrain/dispMappingWithWater_fsh.glsl b/data/shader/terrain/dispMappingWithWater_fsh.glsl index dfee3af..861ad4b 100644 --- a/data/shader/terrain/dispMappingWithWater_fsh.glsl +++ b/data/shader/terrain/dispMappingWithWater_fsh.glsl @@ -39,7 +39,7 @@ vec4 color() { vec4 color = vec4(1.); if (isWater) { - color = vec4(0., 0.2, 1., 1.); + color = vec4(0.1, 0.2, 1., 1.); } else { color = texture_getColor(pos, gNormal, approxNormal); } @@ -51,7 +51,7 @@ vec4 color() { vec4 emissive() { //return uEmissiveColor; if (isWater) { - return vec4(0., 0., 0.01, 1.); + return vec4(0., 0., 0.001, 1.); } else { return vec4(0., 0., 0., 1.); } diff --git a/extern/acgl/include/ACGL/OpenGL/Objects/Texture.hh b/extern/acgl/include/ACGL/OpenGL/Objects/Texture.hh index 022bc61..14e35d9 100644 --- a/extern/acgl/include/ACGL/OpenGL/Objects/Texture.hh +++ b/extern/acgl/include/ACGL/OpenGL/Objects/Texture.hh @@ -76,15 +76,16 @@ public: glGenTextures(1, &mObjectName); } - TextureBase(GLenum _target, GLenum _internalFormat) - : mObjectName(0), - mTarget(_target), - mWidth(0), - mHeight(1), - mDepth(1), - mInternalFormat(_internalFormat) - { + TextureBase(GLenum _target, GLenum _internalFormat, GLuint _objectName = 0) + : mObjectName(_objectName), + mTarget(_target), + mWidth(0), + mHeight(1), + mDepth(1), + mInternalFormat(_internalFormat) { + if (mObjectName == 0) { glGenTextures(1, &mObjectName); + } } virtual ~TextureBase(void) @@ -348,7 +349,10 @@ ACGL_SMARTPOINTER_TYPEDEFS(Texture1D) class Texture2D : public TextureBase { public: + Texture2D(GLenum _internalFormat, GLuint _objectName) : TextureBase(GL_TEXTURE_2D, _internalFormat, _objectName) {} + Texture2D( GLenum _internalFormat = GL_RGBA ) : TextureBase( GL_TEXTURE_2D, _internalFormat ) {} + Texture2D( const glm::uvec2 &_size, GLenum _internalFormat = GL_RGBA ) : TextureBase( GL_TEXTURE_2D, _internalFormat ) { resize( _size ); @@ -372,7 +376,10 @@ ACGL_SMARTPOINTER_TYPEDEFS(Texture2D) class Texture3D : public TextureBase { public: + Texture3D(GLenum _internalFormat, GLuint _objectName) : TextureBase(GL_TEXTURE_3D, _internalFormat, _objectName) {} + Texture3D( GLenum _internalFormat = GL_RGBA ) : TextureBase( GL_TEXTURE_3D, _internalFormat ) {} + Texture3D( const glm::uvec3 &_size, GLenum _internalFormat = GL_RGBA ) : TextureBase( GL_TEXTURE_3D, _internalFormat ) { resize( _size ); diff --git a/src/game/src/engine/scene/OrbitalSimulationSystem.cpp b/src/game/src/engine/scene/OrbitalSimulationSystem.cpp index ce77d24..9ac628d 100644 --- a/src/game/src/engine/scene/OrbitalSimulationSystem.cpp +++ b/src/game/src/engine/scene/OrbitalSimulationSystem.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include @@ -23,6 +24,56 @@ #include +// ---------------------------------------------------------------------------- +// PRECOMPUTATIONS +// ---------------------------------------------------------------------------- +const float Rg = 6360.0; +const float Rt = 6420.0; +const float RL = 6421.0; + +const int TRANSMITTANCE_W = 256; +const int TRANSMITTANCE_H = 64; + +const int SKY_W = 64; +const int SKY_H = 16; + +const int RES_R = 32; +const int RES_MU = 128; +const int RES_MU_S = 32; +const int RES_NU = 8; + +const int reflectanceUnit = 0; +const int transmittanceUnit = 1; +const int irradianceUnit = 2; +const int inscatterUnit = 3; +const int deltaEUnit = 4; +const int deltaSRUnit = 5; +const int deltaSMUnit = 6; +const int deltaJUnit = 7; + +unsigned int reflectanceTexture;//unit 0, ground reflectance texture +unsigned int transmittanceTexture;//unit 1, T table +unsigned int irradianceTexture;//unit 2, E table +unsigned int inscatterTexture;//unit 3, S table +unsigned int deltaETexture;//unit 4, deltaE table +unsigned int deltaSRTexture;//unit 5, deltaS table (Rayleigh part) +unsigned int deltaSMTexture;//unit 6, deltaS table (Mie part) +unsigned int deltaJTexture;//unit 7, deltaJ table + +unsigned int transmittanceProg; +unsigned int irradiance1Prog; +unsigned int inscatter1Prog; +unsigned int copyIrradianceProg; +unsigned int copyInscatter1Prog; +unsigned int jProg; +unsigned int irradianceNProg; +unsigned int inscatterNProg; +unsigned int copyInscatterNProg; + +unsigned int fbo; + +unsigned int drawProg; + bool OrbitalSimulationSystem::startup() { RESOLVE_DEPENDENCY(m_settings); RESOLVE_DEPENDENCY(m_events); @@ -31,6 +82,13 @@ bool OrbitalSimulationSystem::startup() { m_events->subscribe([this](const SimulateEvent &e) { simulateOrbitals(e.dt, e.totalTime); }); + m_events->subscribe<"DrawUI"_sh>([this]() { + ImGui::Begin("Atmosphere", 0, ImGuiWindowFlags_AlwaysAutoResize); + ImGui::Image((void*)transmittanceTexture, glm::vec2{ TRANSMITTANCE_W, TRANSMITTANCE_H }, { 0, 1 }, { 1, 0 }); + ImGui::End(); + + }, 1); + // Make a default geometry (TODO) auto vaoSphere = VertexArrayObjectCreator("icosphere2.obj").create(); auto texture = Texture2DFileManager::the()->get(Texture2DCreator("checkerboard.png")); @@ -242,292 +300,271 @@ void OrbitalSimulationSystem::simulateOrbitals(float dt, float totalTime) { } } -glm::dvec2 getDistance(glm::dvec3 o, glm::dvec3 l, double r) { - double a = -glm::dot(l, o); - double d = a * a - glm::dot(o, o) + r * r; - if (d < 0) { - return {-1, -1}; +using namespace std; +string* loadFile(const string &fileName) { + string* result = new string(); + ifstream file(fileName.c_str()); + if (!file) { + std::cerr << "Cannot open file " << fileName << endl; + throw exception(); } - - double dsqrt = glm::sqrt(d); - - // We are inside the atmosphere, return the distance to the far intersection - return {a + dsqrt, a - dsqrt}; -} - -const int TRANSMITTANCE_INTEGRAL_SAMPLES = 25; - -glm::dvec3 computeTransmittance(glm::dvec3 constBetaRayleigh, - double constBetaMie, double planetRadius, - double atmosphereRadius, - double averageDensityHeight, double viewAngle, - double height) { - viewAngle = glm::acos(viewAngle); - - glm::dvec3 pos = {0, height, 0}; - glm::dvec3 dir = glm::normalize(glm::rotateX(glm::dvec3(0, 1, 0), viewAngle)); - - glm::dvec3 result = {0, 0, 0}; - int stepCount = TRANSMITTANCE_INTEGRAL_SAMPLES; - double stepSize = getDistance(pos, dir, atmosphereRadius).x / stepCount; - - double opticalDepthRayleigh = 0; - double opticalDepthMie = 0; - - // Integrate - for (int i = 0; i < stepCount; i++) { - height = glm::length(pos); - - double altitude = (height - planetRadius); - - opticalDepthRayleigh += glm::exp(-altitude / 8000.0) * stepSize; - opticalDepthMie += glm::exp(-altitude / 1200.0) * stepSize; - - pos += dir * stepSize; + string line; + while (getline(file, line)) { + *result += line; + *result += '\n'; } - result = constBetaRayleigh * opticalDepthRayleigh + - glm::dvec3(constBetaMie / 0.9) * opticalDepthMie; - result = glm::exp(-result); + file.close(); return result; } -SharedTextureData computeTransmittanceTexture(double atmosphereRadius, - double planetRadius, - glm::dvec3 constBetaRayleigh, - double constBetaMie, - double averageDensityHeight) { - const int WIDTH = 256; - const int HEIGHT = 64; - - auto transmittance = new GLubyte[WIDTH * HEIGHT * 4]; - auto pixel = transmittance; - for (int y = 0; y < HEIGHT; y++) { - double r = (double(y + 1) / HEIGHT); - double height = r * r * (atmosphereRadius - planetRadius) + planetRadius; - for (int x = 0; x < WIDTH; x++) { - double viewAngle = -0.15 + - glm::tan(1.5 * (double(x + 1) / WIDTH)) / - glm::tan(1.5) * (1.0 + 0.15); - - auto val = - glm::vec4(computeTransmittance(constBetaRayleigh, constBetaMie, - planetRadius, atmosphereRadius, - averageDensityHeight - planetRadius, - viewAngle, height), - 1.0); - pixel[0] = (GLubyte)(val.r * 255); - pixel[1] = (GLubyte)(val.g * 255); - pixel[2] = (GLubyte)(val.b * 255); - pixel[3] = (GLubyte)(val.a * 255); - pixel += 4; - } +void printShaderLog(int shaderId) { + int logLength; + glGetShaderiv(shaderId, GL_INFO_LOG_LENGTH, &logLength); + if (logLength > 0) { + char *log = new char[logLength]; + glGetShaderInfoLog(shaderId, logLength, &logLength, log); + cout << string(log); + delete[] log; } - - SharedTextureData result = std::make_shared(); - result->setWidth(WIDTH); - result->setHeight(HEIGHT); - result->setData(transmittance); - result->setFormat(GL_RGBA); - result->setType(GL_UNSIGNED_BYTE); - - //stbi_write_png("transmittance.png", WIDTH, HEIGHT, 4, transmittance, 0); - - return result; } -const int INSCATTERING_INTEGRAL_SAMPLES = 50; - -glm::dvec3 sampleTransmittance(GLubyte* transmittance, double r, - double mu, double Rg, double Rt) { - int h = 56; - int w = 256; - - - double uR = sqrt((r - Rg) / (Rt - Rg)); - double uMu = atan((mu + 0.15) / (1.0 + 0.15) * tan(1.5)) / 1.5; - - int x = int((w - 1) * uMu); - int y = int((h - 1) * uR); +unsigned int loadProgram(const vector &files) { + unsigned int programId = glCreateProgram(); + unsigned int vertexShaderId = glCreateShader(GL_VERTEX_SHADER); + unsigned int fragmentShaderId = glCreateShader(GL_FRAGMENT_SHADER); + glAttachShader(programId, vertexShaderId); + glAttachShader(programId, fragmentShaderId); + + int n = files.size(); + string **strs = new string*[n]; + const char** lines = new const char*[n + 2]; + bool geo = false; + for (int i = 0; i < n; ++i) { + string* s = loadFile("data/shader/atmosphere/" + files[i]); + strs[i] = s; + lines[i + 2] = s->c_str(); + if (strstr(lines[i + 2], "_GEOMETRY_") != NULL) { + geo = true; + } + } - x = glm::clamp(x, 0, w - 1); - y = glm::clamp(y, 0, h - 1); + lines[0] = "#version 130\n"; + lines[1] = "#define _VERTEX_\n"; + glShaderSource(vertexShaderId, n + 2, lines, NULL); + glCompileShader(vertexShaderId); + printShaderLog(vertexShaderId); + + if (geo) { + unsigned geometryShaderId = glCreateShader(GL_GEOMETRY_SHADER); + glAttachShader(programId, geometryShaderId); + lines[1] = "#define _GEOMETRY_\n"; + glShaderSource(geometryShaderId, n + 2, lines, NULL); + glCompileShader(geometryShaderId); + printShaderLog(geometryShaderId); + glProgramParameteri(programId, GL_GEOMETRY_INPUT_TYPE_EXT, GL_TRIANGLES); + glProgramParameteri(programId, GL_GEOMETRY_OUTPUT_TYPE_EXT, GL_TRIANGLE_STRIP); + glProgramParameteri(programId, GL_GEOMETRY_VERTICES_OUT_EXT, 3); + } - auto pixel = transmittance + (x + y * w) * 4; + lines[1] = "#define _FRAGMENT_\n"; + glShaderSource(fragmentShaderId, n + 2, lines, NULL); + glCompileShader(fragmentShaderId); + printShaderLog(fragmentShaderId); - return{ double(pixel[0]) / 255, double(pixel[1]) / 255, double(pixel[2]) / 255 }; -} + glLinkProgram(programId); -// 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 -glm::vec3 sampleTransmittance(GLubyte* transmittance, double r, - double mu, glm::vec3 v, glm::vec3 x0, double Rg, double Rt) { - glm::vec3 result; - float r1 = glm::length(x0); - float mu1 = glm::dot(x0, v) / r; - if (mu > 0.0) { - result = glm::min(sampleTransmittance(transmittance, r, mu, Rg, Rt) / sampleTransmittance(transmittance, r1, mu1, Rg, Rt), 1.0); - } else { - result = glm::min(sampleTransmittance(transmittance, r1, -mu1, Rg, Rt) / sampleTransmittance(transmittance, r, -mu, Rg, Rt), 1.0); + for (int i = 0; i < n; ++i) { + delete strs[i]; } - return result; -} + delete[] strs; + delete[] lines; -// nearest intersection of ray r,mu with ground or top atmosphere boundary -// mu=cos(ray zenith angle at ray origin) -double limit(double r, double mu, double Rg, double RL) { - double dout = -r * mu + glm::sqrt(r * r * (mu * mu - 1.0) + RL * RL); - double delta2 = r * r * (mu * mu - 1.0) + Rg * Rg; - if (delta2 >= 0.0) { - double din = -r * mu - glm::sqrt(delta2); - if (din >= 0.0) { - dout = glm::min(dout, din); - } - } - return dout; + return programId; } -// 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) -glm::dvec3 sampleTransmittance(GLubyte* transmittance, double r, double mu, double d, double Rg, double Rt) { - glm::dvec3 result; - double r1 = glm::sqrt(r * r + d * d + 2.0 * r * mu * d); - double mu1 = (r * mu + d) / r1; - if (mu > 0.0) { - result = glm::min(sampleTransmittance(transmittance, r, mu, Rg, Rt) / sampleTransmittance(transmittance, r1, mu1, Rg, Rt), 1.0); - } else { - result = glm::min(sampleTransmittance(transmittance, r1, -mu1, Rg, Rt) / sampleTransmittance(transmittance, r, -mu, Rg, Rt), 1.0); - } - return result; +void drawQuad() { + VertexArrayObject vao; + vao.bind(); + glDrawArrays(GL_TRIANGLE_STRIP, 0, 4); } -void integrand(GLubyte* transmittance, double r, double mu, double muS, double nu, double t, double Rg, double Rt, double HR, double HM, glm::dvec3& ray, glm::dvec3& mie) { - ray = glm::dvec3(0.0); - mie = glm::dvec3(0.0); - - double ri = glm::sqrt(r * r + t * t + 2.0 * r * mu * t); - double muSi = (nu * t + muS * r) / ri; - ri = glm::max(Rg, ri); - if (muSi >= -glm::sqrt(1.0 - Rg * Rg / (ri * ri))) { - glm::vec3 ti = sampleTransmittance(transmittance, r, mu, t, Rg, Rt) * sampleTransmittance(transmittance, ri, muSi, Rg, Rt); - ray = exp(-(ri - Rg) / HR) * ti; - mie = exp(-(ri - Rg) / HM) * ti; - } -} -glm::dvec4 inscatter(GLubyte* transmittance, float r, float mu, float muS, float nu, double Rg, double Rt, double RL, double HR, double HM, glm::dvec3 betaR, double betaM) { - glm::dvec3 ray = glm::dvec3(0.0); - glm::dvec3 mie = glm::dvec3(0.0); - double dx = limit(r, mu, Rg, RL) / double(INSCATTERING_INTEGRAL_SAMPLES); - glm::dvec3 rayi; - glm::dvec3 miei; - integrand(transmittance, r, mu, muS, nu, 0.0, Rg, Rt, HR, HM, rayi, miei); - for (int i = 1; i <= INSCATTERING_INTEGRAL_SAMPLES; ++i) { - double xj = double(i) * dx; - glm::dvec3 rayj; - glm::dvec3 miej; - integrand(transmittance, r, mu, muS, nu, xj, Rg, Rt, HR, HM, rayj, miej); - ray += (rayi + rayj) / 2.0 * dx; - mie += (miei + miej) / 2.0 * dx; - rayi = rayj; - miei = miej; - } - ray *= betaR; - mie *= betaM; - return glm::dvec4{ ray, mie.r}; + +void setLayer(unsigned int prog, int layer) { + double r = layer / (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)); + double dmin = Rt - r; + double dmax = sqrt(r * r - Rg * Rg) + sqrt(Rt * Rt - Rg * Rg); + double dminp = r - Rg; + double dmaxp = sqrt(r * r - Rg * Rg); + glUniform1f(glGetUniformLocation(prog, "r"), float(r)); + glUniform4f(glGetUniformLocation(prog, "dhdH"), float(dmin), float(dmax), float(dminp), float(dmaxp)); + glUniform1i(glGetUniformLocation(prog, "layer"), layer); } -SharedTextureData -computeInscatteringTexture(GLubyte* transmittance, - double atmosphereRadius, double planetRadius, - glm::dvec3 constBetaRayleigh, double constBetaMie, - double g) { - const int WIDTH = 256; - const int HEIGHT = 128; - const int DEPTH = 32; - const int VIEW_SUN_LAYERS = 8; - const int PIXELS_PER_SUN_LAYER = WIDTH / VIEW_SUN_LAYERS; - - auto inscattering = new float[WIDTH * HEIGHT * DEPTH * 4]; - SharedTextureData inscatteringData = std::make_shared(); - inscatteringData->setWidth(WIDTH); - inscatteringData->setHeight(HEIGHT); - inscatteringData->setDepth(DEPTH); - inscatteringData->setData((GLubyte*)inscattering); - inscatteringData->setFormat(GL_RGBA); - inscatteringData->setType(GL_FLOAT); - - double Rg = planetRadius; - double Rt = atmosphereRadius; - - std::thread computeThreads[DEPTH]; - - for (int z = 0; z < DEPTH; z++) { - - computeThreads[z] = std::thread([&](int z) { - - double r = double(z + 1) / (DEPTH); - r = r * r; - r = sqrt(Rg * Rg + r * (Rt * Rt - Rg * Rg)) + - (z == 0 ? 0.01 : (z == DEPTH - 1 ? -0.001 : 0.0)); - - double dmin = Rt - r; - double dmax = sqrt(r * r - Rg * Rg) + sqrt(Rt * Rt - Rg * Rg); - double dminp = r - Rg; - double dmaxp = sqrt(r * r - Rg * Rg); - - glm::dvec4 dhdH = {dmin, dmax, dminp, dmaxp}; - - auto pixel = inscattering + z * WIDTH * HEIGHT * 4; - for (int y = 0; y < HEIGHT; y++) { - - for (int x = 0; x < WIDTH; x++) { - - double mu = 0; - double muS = 0; - double nu = 0; - - if (double(y) < double(HEIGHT) / 2.0) { - double d = 1.0 - double(y) / (double(HEIGHT) / 2.0 - 1.0); - d = glm::min(glm::max(dhdH.z, d * dhdH.w), dhdH.w * 0.999); - mu = (Rg * Rg - r * r - d * d) / (2.0 * r * d); - mu = glm::min(mu, -sqrt(1.0 - (Rg / r) * (Rg / r)) - 0.001); - } else { - double d = (double(y) - double(HEIGHT) / 2.0) / - (double(HEIGHT) / 2.0 - 1.0); - d = glm::min(glm::max(dhdH.x, d * dhdH.y), dhdH.y * 0.999); - mu = (Rt * Rt - r * r - d * d) / (2.0 * r * d); - } - muS = glm::mod(double(x), double(PIXELS_PER_SUN_LAYER)) / - (double(PIXELS_PER_SUN_LAYER) - 1.0); - muS = glm::tan((2.0 * muS - 1.0 + 0.26) * 1.1) / glm::tan(1.26 * 1.1); - nu = -1.0 + - glm::floor(double(x) / double(PIXELS_PER_SUN_LAYER)) / - (double(VIEW_SUN_LAYERS) - 1.0) * 2.0; - - auto val = - glm::vec4(inscatter(transmittance, r, mu, muS, nu, planetRadius, - atmosphereRadius, atmosphereRadius + 100, - 8000, 1200, constBetaRayleigh, constBetaMie)); - pixel[0] = val.r; - pixel[1] = val.g; - pixel[2] = val.b; - pixel[3] = val.a; - pixel += 4; - } - } - }, z); - - //stbi_write_png(("inscattering_" + std::to_string(z) + ".png").c_str(), WIDTH, HEIGHT, 4, inscattering + z * WIDTH * HEIGHT * 4, 0); +void precompute() { + glActiveTexture(GL_TEXTURE0 + transmittanceUnit); + glGenTextures(1, &transmittanceTexture); + glBindTexture(GL_TEXTURE_2D, transmittanceTexture); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE); + glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0); + glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB16F, TRANSMITTANCE_W, TRANSMITTANCE_H, 0, GL_RGB, GL_FLOAT, NULL); + + glActiveTexture(GL_TEXTURE0 + irradianceUnit); + glGenTextures(1, &irradianceTexture); + glBindTexture(GL_TEXTURE_2D, irradianceTexture); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE); + glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0); + glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB16F, SKY_W, SKY_H, 0, GL_RGB, GL_FLOAT, NULL); + + glActiveTexture(GL_TEXTURE0 + inscatterUnit); + glGenTextures(1, &inscatterTexture); + glBindTexture(GL_TEXTURE_3D, inscatterTexture); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE); + glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0); + glTexImage3D(GL_TEXTURE_3D, 0, GL_RGBA16F, RES_MU_S * RES_NU, RES_MU, RES_R, 0, GL_RGB, GL_FLOAT, NULL); + + glActiveTexture(GL_TEXTURE0 + deltaEUnit); + glGenTextures(1, &deltaETexture); + glBindTexture(GL_TEXTURE_2D, deltaETexture); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE); + glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0); + glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB16F, SKY_W, SKY_H, 0, GL_RGB, GL_FLOAT, NULL); + + glActiveTexture(GL_TEXTURE0 + deltaSRUnit); + glGenTextures(1, &deltaSRTexture); + glBindTexture(GL_TEXTURE_3D, deltaSRTexture); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE); + glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0); + glTexImage3D(GL_TEXTURE_3D, 0, GL_RGB16F, RES_MU_S * RES_NU, RES_MU, RES_R, 0, GL_RGB, GL_FLOAT, NULL); + + glActiveTexture(GL_TEXTURE0 + deltaSMUnit); + glGenTextures(1, &deltaSMTexture); + glBindTexture(GL_TEXTURE_3D, deltaSMTexture); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE); + glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0); + glTexImage3D(GL_TEXTURE_3D, 0, GL_RGB16F, RES_MU_S * RES_NU, RES_MU, RES_R, 0, GL_RGB, GL_FLOAT, NULL); + + glActiveTexture(GL_TEXTURE0 + deltaJUnit); + glGenTextures(1, &deltaJTexture); + glBindTexture(GL_TEXTURE_3D, deltaJTexture); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE); + glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0); + glTexImage3D(GL_TEXTURE_3D, 0, GL_RGB16F, RES_MU_S * RES_NU, RES_MU, RES_R, 0, GL_RGB, GL_FLOAT, NULL); + + vector files; + files.push_back("common.glsl"); + files.push_back("transmittance.glsl"); + transmittanceProg = loadProgram(files); + + files.clear(); + files.push_back("common.glsl"); + files.push_back("irradiance1.glsl"); + irradiance1Prog = loadProgram(files); + + files.clear(); + files.push_back("common.glsl"); + files.push_back("inscatter1.glsl"); + inscatter1Prog = loadProgram(files); + + files.clear(); + files.push_back("common.glsl"); + files.push_back("copyIrradiance.glsl"); + copyIrradianceProg = loadProgram(files); + + files.clear(); + files.push_back("common.glsl"); + files.push_back("copyInscatter1.glsl"); + copyInscatter1Prog = loadProgram(files); + + files.clear(); + files.push_back("common.glsl"); + files.push_back("inscatterS.glsl"); + jProg = loadProgram(files); + + files.clear(); + files.push_back("common.glsl"); + files.push_back("irradianceN.glsl"); + irradianceNProg = loadProgram(files); + + files.clear(); + files.push_back("common.glsl"); + files.push_back("inscatterN.glsl"); + inscatterNProg = loadProgram(files); + + files.clear(); + files.push_back("common.glsl"); + files.push_back("copyInscatterN.glsl"); + copyInscatterNProg = loadProgram(files); + + glGenFramebuffers(1, &fbo); + glBindFramebuffer(GL_FRAMEBUFFER, fbo); + glReadBuffer(GL_COLOR_ATTACHMENT0); + glDrawBuffer(GL_COLOR_ATTACHMENT0); + + // computes transmittance texture T (line 1 in algorithm 4.1) + glFramebufferTexture(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, transmittanceTexture, 0); + glViewport(0, 0, TRANSMITTANCE_W, TRANSMITTANCE_H); + glUseProgram(transmittanceProg); + drawQuad(); + + // computes single scattering texture deltaS (line 3 in algorithm 4.1) + // Rayleigh and Mie separated in deltaSR + deltaSM + glFramebufferTexture(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, deltaSRTexture, 0); + glFramebufferTexture(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT1, deltaSMTexture, 0); + unsigned int bufs[2] = { GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1 }; + glDrawBuffers(2, bufs); + glViewport(0, 0, RES_MU_S * RES_NU, RES_MU); + glUseProgram(inscatter1Prog); + glUniform1i(glGetUniformLocation(inscatter1Prog, "transmittanceSampler"), transmittanceUnit); + for (int layer = 0; layer < RES_R; ++layer) { + setLayer(inscatter1Prog, layer); + drawQuad(); } - - for (auto& thread : computeThreads) { - thread.join(); + glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT1, GL_TEXTURE_2D, 0, 0); + glDrawBuffer(GL_COLOR_ATTACHMENT0); + + // copies deltaS into inscatter texture S (line 5 in algorithm 4.1) + glFramebufferTexture(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, inscatterTexture, 0); + glViewport(0, 0, RES_MU_S * RES_NU, RES_MU); + glUseProgram(copyInscatter1Prog); + glUniform1i(glGetUniformLocation(copyInscatter1Prog, "deltaSRSampler"), deltaSRUnit); + glUniform1i(glGetUniformLocation(copyInscatter1Prog, "deltaSMSampler"), deltaSMUnit); + for (int layer = 0; layer < RES_R; ++layer) { + setLayer(copyInscatter1Prog, layer); + drawQuad(); } - return inscatteringData; + glBindFramebuffer(GL_FRAMEBUFFER, 0); + + glFinish(); } void precomputeAtmosphereLookupTextures(double molecularNumberDensity, @@ -536,43 +573,16 @@ void precomputeAtmosphereLookupTextures(double molecularNumberDensity, double atmosphereRadius, double averageDensityHeight, double g, SharedTexture2D& transmittanceText, SharedTexture3D& inscatteringText) { - const double wavelengthIndependentFactor = - (8.0 * glm::pow3(M_PI) * glm::pow2(glm::pow3(airIOR) - 1.0)); - const glm::dvec3 constBetaRayleigh = - wavelengthIndependentFactor / - (3 * molecularNumberDensity * glm::pow4(rgbWavelengths)); - const double constBetaMie = wavelengthIndependentFactor / - (3 * molecularNumberDensity * glm::pow(10.0, 36.0)); - - SharedTextureData transmittance; - auto transmittanceCompute = std::thread([&]() { - transmittance = computeTransmittanceTexture(atmosphereRadius, planetRadius, - constBetaRayleigh, constBetaMie, - averageDensityHeight); - }); - transmittanceCompute.join(); - - - SharedTextureData inscattering; - - auto inscatteringCompute = std::thread([&]() { - inscattering = - computeInscatteringTexture(transmittance->getData(), atmosphereRadius, planetRadius, - constBetaRayleigh, constBetaMie, g); - }); - - inscatteringCompute.join(); - - transmittanceText = std::make_shared(GL_RGBA); + + precompute(); + transmittanceText = std::make_shared(GL_RGBA, transmittanceTexture); transmittanceText->setMinFilter(GL_LINEAR); transmittanceText->setMagFilter(GL_LINEAR); - transmittanceText->setImageData(transmittance); - inscatteringText = std::make_shared(GL_RGBA); + inscatteringText = std::make_shared(GL_RGBA, inscatterTexture); inscatteringText->setMinFilter(GL_LINEAR); inscatteringText->setMagFilter(GL_LINEAR); - inscatteringText->setImageData(inscattering); } // Mass is in solar masses! -- GitLab