ieesqrt

7 files changed, 62 insertions, 35 deletions

diff --git a/common.glsl b/common.glsl
index 3d28814..13dd0c2 100644
--- a/common.glsl
+++ b/common.glsl
@@ -83,10 +83,37 @@ const float M_PI = 3.141592657;
 

 uniform sampler2D transmittanceSampler;

 

+float ieeesqrt(float v)

+{

+  float ret;

+  if (v < 0.0) {

+    ret = 0.0/0.0;

+  } else {

+    ret = sqrt(v);

+  }

+  return ret;

+}

+

+vec2 ieeesqrt(vec2 v)

+{

+  return vec2(ieeesqrt(v.x), ieeesqrt(v.y));

+}

+

+float ieeepow(float b, float e)

+{

+  float ret;

+  if (b < 0.0 && floor(e) != e) {

+    ret = 0.0/0.0;

+  } else {

+    ret = pow(b, e);

+  }

+  return ret;

+}

+

 vec2 getTransmittanceUV(float r, float mu) {

     float uR, uMu;

 #ifdef TRANSMITTANCE_NON_LINEAR

-	uR = sqrt((r - Rg) / (Rt - Rg));

+	uR = ieeesqrt((r - Rg) / (Rt - Rg));

 	uMu = atan((mu + 0.15) / (1.0 + 0.15) * tan(1.5)) / 1.5;

 #else

 	uR = (r - Rg) / (Rt - Rg);

@@ -110,14 +137,14 @@ void getIrradianceRMuS(out float r, out float muS) {
 

 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);

+    float H = ieeesqrt(Rt * Rt - Rg * Rg);

+    float rho = ieeesqrt(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));

+    float uMu = cst.w + (rmu * cst.x + ieeesqrt(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

@@ -147,7 +174,7 @@ void getMuMuSNu(float r, vec4 dhdH, out float mu, out float muS, out float nu) {
         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);

+        mu = min(mu, -ieeesqrt(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);

@@ -174,10 +201,10 @@ void getMuMuSNu(float r, vec4 dhdH, out float mu, out float muS, out float nu) {
 // 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 dout = -r * mu + ieeesqrt(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);

+        float din = -r * mu - ieeesqrt(delta2);

         if (din >= 0.0) {

             dout = min(dout, din);

         }

@@ -195,7 +222,7 @@ vec3 transmittance(float r, float mu) {
 // 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);

+    return mu < -ieeesqrt(1.0 - (Rg / r) * (Rg / r)) ? vec3(0.0) : transmittance(r, mu);

 }

 

 // transmittance(=transparency) of atmosphere between x and x0

@@ -217,13 +244,13 @@ vec3 transmittance(float r, float mu, vec3 v, vec3 x0) {
 // (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);

+    float a = ieeesqrt((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)));

+    vec2 y = a01s / (2.3193*abs(a01) + ieeesqrt(1.52*a01sq + 4.0)) * vec2(1.0, exp(-d/H*(d/(2.0*r)+mu)));

+    return ieeesqrt((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

@@ -238,7 +265,7 @@ vec3 analyticTransmittance(float r, float mu, float d) {
 // 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 r1 = ieeesqrt(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);

@@ -260,7 +287,7 @@ float phaseFunctionR(float 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);

+	return 1.5 * 1.0 / (4.0 * M_PI) * (1.0 - mieG*mieG) * ieeepow(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")

diff --git a/earth.glsl b/earth.glsl
index e4a2cac..2b37122 100644
--- a/earth.glsl
+++ b/earth.glsl
@@ -64,7 +64,7 @@ vec3 inscatter(inout vec3 x, inout float t, vec3 v, vec3 s, out float r, out flo
     vec3 result;

     r = length(x);

     mu = dot(x, v) / r;

-    float d = -r * mu - sqrt(r * r * (mu * mu - 1.0) + Rt * Rt);

+    float d = -r * mu - ieeesqrt(r * r * (mu * mu - 1.0) + Rt * Rt);

     if (d > 0.0) { // if x in space and ray intersects atmosphere

         // move x to nearest intersection of ray with top atmosphere boundary

         x += d * v;

@@ -96,19 +96,19 @@ vec3 inscatter(inout vec3 x, inout float t, vec3 v, vec3 s, out float r, out flo
 #ifdef FIX

                 // avoids imprecision problems near horizon by interpolating between two points above and below horizon

                 const float EPS = 0.004;

-                float muHoriz = -sqrt(1.0 - (Rg / r) * (Rg / r));

+                float muHoriz = -ieeesqrt(1.0 - (Rg / r) * (Rg / r));

                 if (abs(mu - muHoriz) < EPS) {

                     float a = ((mu - muHoriz) + EPS) / (2.0 * EPS);

 

                     mu = muHoriz - EPS;

-                    r0 = sqrt(r * r + t * t + 2.0 * r * t * mu);

+                    r0 = ieeesqrt(r * r + t * t + 2.0 * r * t * mu);

                     mu0 = (r * mu + t) / r0;

                     vec4 inScatter0 = texture4D(inscatterSampler, r, mu, muS, nu);

                     vec4 inScatter1 = texture4D(inscatterSampler, r0, mu0, muS0, nu);

                     vec4 inScatterA = max(inScatter0 - attenuation.rgbr * inScatter1, 0.0);

 

                     mu = muHoriz + EPS;

-                    r0 = sqrt(r * r + t * t + 2.0 * r * t * mu);

+                    r0 = ieeesqrt(r * r + t * t + 2.0 * r * t * mu);

                     mu0 = (r * mu + t) / r0;

                     inScatter0 = texture4D(inscatterSampler, r, mu, muS, nu);

                     inScatter1 = texture4D(inscatterSampler, r0, mu0, muS0, nu);

@@ -159,8 +159,8 @@ vec3 groundColor(vec3 x, float t, vec3 v, vec3 s, float r, float mu, vec3 attenu
         // water specular color due to sunLight

         if (reflectance.w > 0.0) {

             vec3 h = normalize(s - v);

-            float fresnel = 0.02 + 0.98 * pow(1.0 - dot(-v, h), 5.0);

-            float waterBrdf = fresnel * pow(max(dot(h, n), 0.0), 150.0);

+            float fresnel = 0.02 + 0.98 * ieeepow(1.0 - dot(-v, h), 5.0);

+            float waterBrdf = fresnel * ieeepow(max(dot(h, n), 0.0), 150.0);

             groundColor += reflectance.w * max(waterBrdf, 0.0) * sunLight * ISun;

         }

 

@@ -184,9 +184,9 @@ vec3 sunColor(vec3 x, float t, vec3 v, vec3 s, float r, float mu) {
 

 vec3 HDR(vec3 L) {

     L = L * exposure;

-    L.r = L.r < 1.413 ? pow(L.r * 0.38317, 1.0 / 2.2) : 1.0 - exp(-L.r);

-    L.g = L.g < 1.413 ? pow(L.g * 0.38317, 1.0 / 2.2) : 1.0 - exp(-L.g);

-    L.b = L.b < 1.413 ? pow(L.b * 0.38317, 1.0 / 2.2) : 1.0 - exp(-L.b);

+    L.r = L.r < 1.413 ? ieeepow(L.r * 0.38317, 1.0 / 2.2) : 1.0 - exp(-L.r);

+    L.g = L.g < 1.413 ? ieeepow(L.g * 0.38317, 1.0 / 2.2) : 1.0 - exp(-L.g);

+    L.b = L.b < 1.413 ? ieeepow(L.b * 0.38317, 1.0 / 2.2) : 1.0 - exp(-L.b);

     return L;

 }

 

@@ -196,13 +196,13 @@ void main() {
 

     float r = length(x);

     float mu = dot(x, v) / r;

-    float t = -r * mu - sqrt(r * r * (mu * mu - 1.0) + Rg * Rg);

+    float t = -r * mu - ieeesqrt(r * r * (mu * mu - 1.0) + Rg * Rg);

 

     vec3 g = x - vec3(0.0, 0.0, Rg + 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);

     float c = g.x * g.x + g.y * g.y - g.z * g.z;

-    float d = -(b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);

+    float d = -(b + ieeesqrt(b * b - 4.0 * a * c)) / (2.0 * a);

     bool cone = d > 0.0 && abs(x.z + d * v.z - Rg) <= 10.0;

 

     if (t > 0.0) {

diff --git a/inscatter1.glsl b/inscatter1.glsl
index 5ec4d1a..6eb2486 100644
--- a/inscatter1.glsl
+++ b/inscatter1.glsl
@@ -69,10 +69,10 @@ void main() {
 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 ri = ieeesqrt(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))) {

+    if (muSi >= -ieeesqrt(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;

diff --git a/inscatterN.glsl b/inscatterN.glsl
index c1d49c3..3ab2323 100644
--- a/inscatterN.glsl
+++ b/inscatterN.glsl
@@ -69,7 +69,7 @@ void main() {
 #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 ri = ieeesqrt(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);

diff --git a/inscatterS.glsl b/inscatterS.glsl
index e71f019..668afa3 100644
--- a/inscatterS.glsl
+++ b/inscatterS.glsl
@@ -78,14 +78,14 @@ 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);

+    float var = ieeesqrt(1.0 - mu * mu) * ieeesqrt(1.0 - muS * muS);

     nu = clamp(nu, muS * mu - var, muS * mu + var);

 

-    float cthetamin = -sqrt(1.0 - (Rg / r) * (Rg / r));

+    float cthetamin = -ieeesqrt(1.0 - (Rg / r) * (Rg / r));

 

-    vec3 v = vec3(sqrt(1.0 - mu * mu), 0.0, mu);

+    vec3 v = vec3(ieeesqrt(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);

+    vec3 s = vec3(sx, ieeesqrt(max(0.0, 1.0 - sx * sx - muS * muS)), muS);

 

     raymie = vec3(0.0);

 

@@ -100,7 +100,7 @@ void inscatter(float r, float mu, float muS, float nu, out vec3 raymie) {
         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);

+            dground = -r * ctheta - ieeesqrt(r * r * (ctheta * ctheta - 1.0) + Rg * Rg);

             gtransp = transmittance(Rg, -(r * ctheta + dground) / Rg, dground);

         }

 

diff --git a/irradianceN.glsl b/irradianceN.glsl
index 5e98182..c3a9289 100644
--- a/irradianceN.glsl
+++ b/irradianceN.glsl
@@ -54,7 +54,7 @@ const float dtheta = M_PI / float(IRRADIANCE_INTEGRAL_SAMPLES);
 void main() {

     float r, muS;

     getIrradianceRMuS(r, muS);

-    vec3 s = vec3(sqrt(1.0 - muS * muS), 0.0, muS);

+    vec3 s = vec3(ieeesqrt(1.0 - muS * muS), 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)

diff --git a/transmittance.glsl b/transmittance.glsl
index 9115add..2ede3b2 100644
--- a/transmittance.glsl
+++ b/transmittance.glsl
@@ -61,12 +61,12 @@ float opticalDepth(float H, float r, float mu) {
     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);

+        float yj = exp(-(ieeesqrt(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;

+    return mu < -ieeesqrt(1.0 - (Rg / r) * (Rg / r)) ? 1e9 : result;

 }

 

 void main() {