2013-12-08 17:36:39 +00:00
|
|
|
#include "AtmosphereModelBruneton.h"
|
|
|
|
|
2012-12-02 11:08:56 +00:00
|
|
|
/*
|
|
|
|
* Atmospheric scattering, based on E. Bruneton and F.Neyret work.
|
|
|
|
* http://evasion.inrialpes.fr/~Eric.Bruneton/
|
|
|
|
*/
|
|
|
|
|
2013-11-07 08:37:11 +00:00
|
|
|
#include <cassert>
|
|
|
|
#include <cmath>
|
|
|
|
#include <cstdio>
|
|
|
|
#include <cstdlib>
|
2015-12-30 23:36:22 +00:00
|
|
|
#include "Maths.h"
|
2013-10-26 21:19:49 +00:00
|
|
|
#include "System.h"
|
2013-12-09 10:59:57 +00:00
|
|
|
#include "ParallelWork.h"
|
2013-11-07 08:37:11 +00:00
|
|
|
#include "PackStream.h"
|
2013-12-08 19:54:34 +00:00
|
|
|
#include "Scenery.h"
|
2013-11-15 23:27:40 +00:00
|
|
|
#include "AtmosphereDefinition.h"
|
2013-12-08 19:54:34 +00:00
|
|
|
#include "AtmosphereRenderer.h"
|
2015-09-24 22:12:31 +00:00
|
|
|
#include "AtmosphereResult.h"
|
2013-12-08 19:54:34 +00:00
|
|
|
#include "SoftwareRenderer.h"
|
|
|
|
#include "WaterRenderer.h"
|
2015-12-29 23:20:20 +00:00
|
|
|
#include "CelestialBodyDefinition.h"
|
2013-12-08 19:54:34 +00:00
|
|
|
#include "LightComponent.h"
|
|
|
|
#include "LightStatus.h"
|
2013-12-09 21:16:00 +00:00
|
|
|
#include "Texture2D.h"
|
|
|
|
#include "Texture4D.h"
|
|
|
|
#include "CacheFile.h"
|
2015-09-21 21:12:43 +00:00
|
|
|
#include "FloatNode.h"
|
2013-12-08 19:54:34 +00:00
|
|
|
|
2015-12-21 16:48:30 +00:00
|
|
|
#define WORKAROUND_OFFSET 0.1
|
|
|
|
|
2012-12-13 21:33:42 +00:00
|
|
|
/*********************** Constants ***********************/
|
2012-12-02 11:08:56 +00:00
|
|
|
|
2016-01-14 23:39:33 +00:00
|
|
|
static constexpr double Rg = Scenery::EARTH_RADIUS;
|
|
|
|
static constexpr double Rt = Rg + 60.0;
|
|
|
|
static constexpr double RL = Rt + 1.0;
|
|
|
|
static constexpr double ISun = 100.0;
|
|
|
|
static constexpr double AVERAGE_GROUND_REFLECTANCE = 0.1;
|
2012-12-02 11:08:56 +00:00
|
|
|
|
2013-03-10 19:55:24 +00:00
|
|
|
#if 1
|
2012-12-12 13:21:46 +00:00
|
|
|
#define RES_MU 128
|
|
|
|
#define RES_MU_S 32
|
|
|
|
#define RES_R 32
|
|
|
|
#define RES_NU 8
|
2015-12-21 16:48:30 +00:00
|
|
|
#define SKY_W 256
|
|
|
|
#define SKY_H 64
|
|
|
|
#define TRANSMITTANCE_W 512
|
|
|
|
#define TRANSMITTANCE_H 128
|
2012-12-15 10:14:57 +00:00
|
|
|
#define TRANSMITTANCE_INTEGRAL_SAMPLES 500
|
2015-12-21 16:48:30 +00:00
|
|
|
#define INSCATTER_INTEGRAL_SAMPLES 100
|
|
|
|
#define IRRADIANCE_INTEGRAL_SAMPLES 64
|
|
|
|
#define INSCATTER_SPHERICAL_INTEGRAL_SAMPLES 32
|
2012-12-15 19:45:19 +00:00
|
|
|
#else
|
|
|
|
#define RES_MU 64
|
|
|
|
#define RES_MU_S 16
|
|
|
|
#define RES_R 16
|
|
|
|
#define RES_NU 8
|
|
|
|
#define SKY_W 64
|
|
|
|
#define SKY_H 16
|
|
|
|
#define TRANSMITTANCE_W 256
|
|
|
|
#define TRANSMITTANCE_H 64
|
|
|
|
#define TRANSMITTANCE_INTEGRAL_SAMPLES 100
|
|
|
|
#define INSCATTER_INTEGRAL_SAMPLES 10
|
|
|
|
#define IRRADIANCE_INTEGRAL_SAMPLES 16
|
|
|
|
#define INSCATTER_SPHERICAL_INTEGRAL_SAMPLES 8
|
|
|
|
#endif
|
2012-12-13 21:33:42 +00:00
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
Texture2D *_transmittanceTexture = NULL;
|
|
|
|
Texture2D *_irradianceTexture = NULL;
|
|
|
|
Texture4D *_inscatterTexture = NULL;
|
2012-12-12 13:21:46 +00:00
|
|
|
|
2012-12-24 13:59:17 +00:00
|
|
|
/* Rayleigh */
|
2012-12-02 11:08:56 +00:00
|
|
|
static const double HR = 8.0;
|
|
|
|
static const Color betaR = {5.8e-3, 1.35e-2, 3.31e-2, 1.0};
|
|
|
|
|
2012-12-24 13:59:17 +00:00
|
|
|
/* Mie */
|
2015-12-21 16:48:30 +00:00
|
|
|
// TODO Use the good ones, determined by weather definition
|
2012-12-24 13:59:17 +00:00
|
|
|
/* DEFAULT */
|
2015-12-21 16:48:30 +00:00
|
|
|
/*static const double HM = 1.2;
|
2012-12-12 13:21:46 +00:00
|
|
|
static const Vector3 betaMSca = {4e-3, 4e-3, 4e-3};
|
|
|
|
static const Vector3 betaMEx = {4e-3 / 0.9, 4e-3 / 0.9, 4e-3 / 0.9};
|
2015-12-21 16:48:30 +00:00
|
|
|
static const double mieG = 0.8;*/
|
2012-12-24 13:59:17 +00:00
|
|
|
/* CLEAR SKY */
|
2015-12-21 16:48:30 +00:00
|
|
|
static const double HM = 1.2;
|
2012-12-24 13:59:17 +00:00
|
|
|
static const Vector3 betaMSca = {20e-3, 20e-3, 20e-3};
|
|
|
|
static const Vector3 betaMEx = {20e-3 / 0.9, 20e-3 / 0.9, 20e-3 / 0.9};
|
2015-12-21 16:48:30 +00:00
|
|
|
static const double mieG = 0.76;
|
2012-12-24 13:59:17 +00:00
|
|
|
/* PARTLY CLOUDY */
|
|
|
|
/*static const double HM = 3.0;
|
|
|
|
static const Vector3 betaMSca = {3e-3, 3e-3, 3e-3};
|
|
|
|
static const Vector3 betaMEx = {3e-3 / 0.9, 3e-3 / 0.9, 3e-3 / 0.9};
|
|
|
|
static const double mieG = 0.65;*/
|
2012-12-02 11:08:56 +00:00
|
|
|
|
2012-12-13 21:33:42 +00:00
|
|
|
/*********************** Shader helpers ***********************/
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
#define step(_a_, _b_) ((_b_) < (_a_) ? 0.0 : 1.0)
|
2012-12-02 11:08:56 +00:00
|
|
|
#define sign(_a_) ((_a_) < 0.0 ? -1.0 : ((_a_) > 0.0 ? 1.0 : 0.0))
|
2015-11-09 21:30:46 +00:00
|
|
|
#define mix(_x_, _y_, _a_) ((_x_) * (1.0 - (_a_)) + (_y_) * (_a_))
|
|
|
|
static inline double min(double a, double b) {
|
2012-12-13 21:33:42 +00:00
|
|
|
return a < b ? a : b;
|
|
|
|
}
|
2015-11-09 21:30:46 +00:00
|
|
|
static inline double max(double a, double b) {
|
2012-12-13 21:33:42 +00:00
|
|
|
return a > b ? a : b;
|
|
|
|
}
|
2015-11-09 21:30:46 +00:00
|
|
|
static inline Color vec4mix(Color v1, Color v2, double a) {
|
2012-12-12 13:21:46 +00:00
|
|
|
v1.r = mix(v1.r, v2.r, a);
|
|
|
|
v1.g = mix(v1.g, v2.g, a);
|
|
|
|
v1.b = mix(v1.b, v2.b, a);
|
|
|
|
v1.a = mix(v1.a, v2.a, a);
|
|
|
|
return v1;
|
|
|
|
}
|
2015-11-09 21:30:46 +00:00
|
|
|
static inline double clamp(double x, double minVal, double maxVal) {
|
|
|
|
if (x < minVal) {
|
2012-12-12 13:21:46 +00:00
|
|
|
x = minVal;
|
|
|
|
}
|
|
|
|
return (x > maxVal) ? maxVal : x;
|
|
|
|
}
|
2015-11-09 21:30:46 +00:00
|
|
|
static inline double smoothstep(double edge0, double edge1, double x) {
|
2012-12-12 13:21:46 +00:00
|
|
|
double t = clamp((x - edge0) / (edge1 - edge0), 0.0, 1.0);
|
|
|
|
return t * t * (3.0 - 2.0 * t);
|
|
|
|
}
|
2015-11-09 21:30:46 +00:00
|
|
|
static inline void _fixVec4Min(Color *vec, double minVal) {
|
|
|
|
if (vec->r < minVal) {
|
|
|
|
vec->r = minVal;
|
|
|
|
}
|
|
|
|
if (vec->g < minVal) {
|
|
|
|
vec->g = minVal;
|
|
|
|
}
|
|
|
|
if (vec->b < minVal) {
|
|
|
|
vec->b = minVal;
|
|
|
|
}
|
|
|
|
if (vec->a < minVal) {
|
|
|
|
vec->a = minVal;
|
|
|
|
}
|
2012-12-12 13:21:46 +00:00
|
|
|
}
|
2015-11-09 21:30:46 +00:00
|
|
|
static inline Color vec4max(Color vec, double minVal) {
|
|
|
|
if (vec.r < minVal) {
|
|
|
|
vec.r = minVal;
|
|
|
|
}
|
|
|
|
if (vec.g < minVal) {
|
|
|
|
vec.g = minVal;
|
|
|
|
}
|
|
|
|
if (vec.b < minVal) {
|
|
|
|
vec.b = minVal;
|
|
|
|
}
|
|
|
|
if (vec.a < minVal) {
|
|
|
|
vec.a = minVal;
|
|
|
|
}
|
2012-12-12 13:21:46 +00:00
|
|
|
return vec;
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static inline Vector3 vec3(double x, double y, double z) {
|
2012-12-12 13:21:46 +00:00
|
|
|
Vector3 result;
|
|
|
|
result.x = x;
|
|
|
|
result.y = y;
|
|
|
|
result.z = z;
|
|
|
|
return result;
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static inline Color vec4(double r, double g, double b, double a) {
|
2012-12-12 13:21:46 +00:00
|
|
|
Color result;
|
|
|
|
result.r = r;
|
|
|
|
result.g = g;
|
|
|
|
result.b = b;
|
|
|
|
result.a = a;
|
|
|
|
return result;
|
|
|
|
}
|
|
|
|
|
2012-12-13 21:33:42 +00:00
|
|
|
/*********************** Texture manipulation ***********************/
|
2012-12-12 13:21:46 +00:00
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static Color _texture4D(Texture4D *tex, double r, double mu, double muS, double nu) {
|
|
|
|
if (r < Rg + 0.00000001)
|
|
|
|
r = Rg + 0.00000001;
|
2012-12-12 20:11:57 +00:00
|
|
|
double H = sqrt(Rt * Rt - Rg * Rg);
|
|
|
|
double rho = sqrt(r * r - Rg * Rg);
|
|
|
|
double rmu = r * mu;
|
|
|
|
double delta = rmu * rmu - r * r + Rg * Rg;
|
2015-12-17 00:13:20 +00:00
|
|
|
Color cst = (rmu < 0.0 && delta > 0.0) ? vec4(1.0, 0.0, 0.0, 0.5 - 0.5 / to_double(RES_MU))
|
|
|
|
: vec4(-1.0, H * H, H, 0.5 + 0.5 / to_double(RES_MU));
|
|
|
|
double uR = 0.5 / to_double(RES_R) + rho / H * (1.0 - 1.0 / to_double(RES_R));
|
|
|
|
double uMu = cst.a + (rmu * cst.r + sqrt(delta + cst.g)) / (rho + cst.b) * (0.5 - 1.0 / to_double(RES_MU));
|
2015-11-09 21:30:46 +00:00
|
|
|
double uMuS =
|
2015-12-17 00:13:20 +00:00
|
|
|
0.5 / to_double(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 / to_double(RES_MU_S));
|
2013-03-09 09:54:31 +00:00
|
|
|
|
2013-12-09 21:16:00 +00:00
|
|
|
return tex->getLinear(uMu, uMuS, nu, uR);
|
2012-12-12 20:11:57 +00:00
|
|
|
}
|
2012-12-02 11:08:56 +00:00
|
|
|
|
2012-12-13 21:33:42 +00:00
|
|
|
/*********************** Physics functions ***********************/
|
|
|
|
|
2012-12-02 11:08:56 +00:00
|
|
|
/* Rayleigh phase function */
|
2015-11-09 21:30:46 +00:00
|
|
|
static double _phaseFunctionR(double mu) {
|
2015-12-30 23:36:22 +00:00
|
|
|
return (3.0 / (16.0 * Maths::PI)) * (1.0 + mu * mu);
|
2012-12-02 11:08:56 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
/* Mie phase function */
|
2015-11-09 21:30:46 +00:00
|
|
|
static double _phaseFunctionM(double mu) {
|
2015-12-30 23:36:22 +00:00
|
|
|
return 1.5 * 1.0 / (4.0 * Maths::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);
|
2012-12-02 11:08:56 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
/* approximated single Mie scattering (cf. approximate Cm in paragraph "Angular precision") */
|
2015-11-09 21:30:46 +00:00
|
|
|
static Color _getMie(Color rayMie) {
|
2012-12-12 13:21:46 +00:00
|
|
|
Color result;
|
|
|
|
|
|
|
|
result.r = rayMie.r * rayMie.a / max(rayMie.r, 1e-4) * (betaR.r / betaR.r);
|
|
|
|
result.g = rayMie.g * rayMie.a / max(rayMie.r, 1e-4) * (betaR.r / betaR.g);
|
|
|
|
result.b = rayMie.b * rayMie.a / max(rayMie.r, 1e-4) * (betaR.r / betaR.b);
|
|
|
|
result.a = 1.0;
|
|
|
|
|
|
|
|
return result;
|
2012-12-02 11:08:56 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
/* 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 */
|
2015-11-09 21:30:46 +00:00
|
|
|
static double _opticalDepth(double H, double r, double mu, double d) {
|
2012-12-02 11:08:56 +00:00
|
|
|
double a = sqrt((0.5 / H) * r);
|
|
|
|
double ax = a * (mu);
|
|
|
|
double ay = a * (mu + d / r);
|
|
|
|
double axs = sign(ax);
|
|
|
|
double ays = sign(ay);
|
|
|
|
double axq = ax * ax;
|
|
|
|
double ayq = ay * ay;
|
|
|
|
double x = ays > axs ? exp(axq) : 0.0;
|
|
|
|
double yx = axs / (2.3193 * fabs(ax) + sqrt(1.52 * axq + 4.0));
|
2013-01-16 14:26:46 +00:00
|
|
|
double yy = ays / (2.3193 * fabs(ay) + sqrt(1.52 * ayq + 4.0)) * exp(-d / H * (d / (2.0 * r) + mu));
|
2012-12-02 11:08:56 +00:00
|
|
|
return sqrt((6.2831 * H) * r) * exp((Rg - r) / H) * (x + yx - yy);
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static inline void _getTransmittanceUV(double r, double mu, double *u, double *v) {
|
|
|
|
if (r < Rg + 0.00000001)
|
|
|
|
r = Rg + 0.00000001;
|
2012-12-13 21:33:42 +00:00
|
|
|
double dr = (r - Rg) / (Rt - Rg);
|
2012-12-15 19:45:19 +00:00
|
|
|
*v = sqrt(dr);
|
2012-12-12 20:11:57 +00:00
|
|
|
*u = atan((mu + 0.15) / (1.0 + 0.15) * tan(1.5)) / 1.5;
|
|
|
|
}
|
|
|
|
|
|
|
|
/* transmittance(=transparency) of atmosphere for infinite ray (r,mu)
|
|
|
|
(mu=cos(view zenith angle)), intersections with ground ignored */
|
2015-11-09 21:30:46 +00:00
|
|
|
static Color _transmittance(double r, double mu) {
|
2012-12-12 20:11:57 +00:00
|
|
|
double u, v;
|
|
|
|
_getTransmittanceUV(r, mu, &u, &v);
|
2013-12-09 21:16:00 +00:00
|
|
|
return _transmittanceTexture->getLinear(u, v);
|
2012-12-12 20:11:57 +00:00
|
|
|
}
|
|
|
|
|
2012-12-13 21:33:42 +00:00
|
|
|
/* 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) */
|
2015-11-09 21:30:46 +00:00
|
|
|
static Color _transmittance3(double r, double mu, double d) {
|
2012-12-13 21:33:42 +00:00
|
|
|
Color result, t1, t2;
|
|
|
|
double r1 = sqrt(r * r + d * d + 2.0 * r * mu * d);
|
|
|
|
double mu1 = (r * mu + d) / r1;
|
2015-11-09 21:30:46 +00:00
|
|
|
if (mu > 0.0) {
|
2012-12-13 21:33:42 +00:00
|
|
|
t1 = _transmittance(r, mu);
|
|
|
|
t2 = _transmittance(r1, mu1);
|
2015-11-09 21:30:46 +00:00
|
|
|
} else {
|
2012-12-13 21:33:42 +00:00
|
|
|
t1 = _transmittance(r1, -mu1);
|
|
|
|
t2 = _transmittance(r, -mu);
|
|
|
|
}
|
|
|
|
result.r = min(t1.r / t2.r, 1.0);
|
|
|
|
result.g = min(t1.g / t2.g, 1.0);
|
|
|
|
result.b = min(t1.b / t2.b, 1.0);
|
|
|
|
result.a = 1.0;
|
2012-12-12 13:21:46 +00:00
|
|
|
return result;
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static void _getIrradianceRMuS(double x, double y, double *r, double *muS) {
|
|
|
|
*r = Rg + y *(Rt - Rg);
|
|
|
|
*muS = -0.2 + x *(1.0 + 0.2);
|
2012-12-13 21:33:42 +00:00
|
|
|
}
|
|
|
|
|
2012-12-12 13:21:46 +00:00
|
|
|
/* nearest intersection of ray r,mu with ground or top atmosphere boundary
|
|
|
|
* mu=cos(ray zenith angle at ray origin) */
|
2015-11-09 21:30:46 +00:00
|
|
|
static double _limit(double r, double mu) {
|
2012-12-12 13:21:46 +00:00
|
|
|
double dout = -r * mu + sqrt(r * r * (mu * mu - 1.0) + RL * RL);
|
|
|
|
double delta2 = r * r * (mu * mu - 1.0) + Rg * Rg;
|
2015-11-09 21:30:46 +00:00
|
|
|
if (delta2 >= 0.0) {
|
2012-12-12 20:11:57 +00:00
|
|
|
double din = -r * mu - sqrt(delta2);
|
2012-12-12 13:21:46 +00:00
|
|
|
if (din >= 0.0) {
|
|
|
|
dout = min(dout, din);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return dout;
|
|
|
|
}
|
|
|
|
|
2012-12-13 21:33:42 +00:00
|
|
|
/* 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 */
|
2015-11-09 21:30:46 +00:00
|
|
|
static Vector3 _analyticTransmittance(double r, double mu, double d) {
|
2012-12-13 21:33:42 +00:00
|
|
|
Vector3 result;
|
|
|
|
|
2013-01-16 14:26:46 +00:00
|
|
|
result.x = exp(-betaR.r * _opticalDepth(HR, r, mu, d) - betaMEx.x * _opticalDepth(HM, r, mu, d));
|
|
|
|
result.y = exp(-betaR.g * _opticalDepth(HR, r, mu, d) - betaMEx.y * _opticalDepth(HM, r, mu, d));
|
|
|
|
result.z = exp(-betaR.b * _opticalDepth(HR, r, mu, d) - betaMEx.z * _opticalDepth(HM, r, mu, d));
|
2012-12-13 21:33:42 +00:00
|
|
|
|
|
|
|
return result;
|
|
|
|
}
|
|
|
|
|
|
|
|
/* transmittance(=transparency) of atmosphere for infinite ray (r,mu)
|
|
|
|
(mu=cos(view zenith angle)), or zero if ray intersects ground */
|
2015-11-09 21:30:46 +00:00
|
|
|
static Color _transmittanceWithShadow(double r, double mu) {
|
2012-12-13 21:33:42 +00:00
|
|
|
return mu < -sqrt(1.0 - (Rg / r) * (Rg / r)) ? COLOR_BLACK : _transmittance(r, mu);
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static void _texCoordToMuMuSNu(double x, double y, double z, double r, Color dhdH, double *mu, double *muS,
|
|
|
|
double *nu) {
|
2013-03-09 11:06:39 +00:00
|
|
|
double d;
|
2012-12-30 16:06:11 +00:00
|
|
|
|
2015-12-17 00:13:20 +00:00
|
|
|
x /= to_double(RES_MU);
|
|
|
|
y /= to_double(RES_MU_S);
|
|
|
|
z /= to_double(RES_NU);
|
2013-03-09 11:06:39 +00:00
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
if (x < 0.5) {
|
2013-03-09 11:06:39 +00:00
|
|
|
d = 1.0 - x / 0.5;
|
2012-12-13 21:33:42 +00:00
|
|
|
d = min(max(dhdH.b, d * dhdH.a), dhdH.a * 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);
|
2015-11-09 21:30:46 +00:00
|
|
|
} else {
|
2013-03-09 11:06:39 +00:00
|
|
|
d = (x - 0.5) / 0.5;
|
2012-12-13 21:33:42 +00:00
|
|
|
d = min(max(dhdH.r, d * dhdH.g), dhdH.g * 0.999);
|
|
|
|
*mu = (Rt * Rt - r * r - d * d) / (2.0 * r * d);
|
|
|
|
}
|
2013-03-09 11:06:39 +00:00
|
|
|
*muS = tan((2.0 * y - 1.0 + 0.26) * 1.1) / tan(1.26 * 1.1);
|
|
|
|
*nu = -1.0 + z / 2.0;
|
2012-12-12 13:21:46 +00:00
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static void _getIrradianceUV(double r, double muS, double *uMuS, double *uR) {
|
2012-12-15 10:14:57 +00:00
|
|
|
*uR = (r - Rg) / (Rt - Rg);
|
|
|
|
*uMuS = (muS + 0.2) / (1.0 + 0.2);
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static Color _irradiance(Texture2D *sampler, double r, double muS) {
|
2012-12-15 10:14:57 +00:00
|
|
|
double u, v;
|
|
|
|
_getIrradianceUV(r, muS, &u, &v);
|
2013-12-09 21:16:00 +00:00
|
|
|
return sampler->getLinear(u, v);
|
2012-12-15 10:14:57 +00:00
|
|
|
}
|
|
|
|
|
2012-12-30 16:06:11 +00:00
|
|
|
/*********************** transmittance.glsl ***********************/
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static void _getTransmittanceRMu(double x, double y, double *r, double *muS) {
|
2012-12-30 16:06:11 +00:00
|
|
|
*r = Rg + (y * y) * (Rt - Rg);
|
|
|
|
*muS = -0.15 + tan(1.5 * x) / tan(1.5) * (1.0 + 0.15);
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static double _opticalDepthTransmittance(double H, double r, double mu) {
|
2012-12-30 16:06:11 +00:00
|
|
|
double result = 0.0;
|
2015-12-17 00:13:20 +00:00
|
|
|
double dx = _limit(r, mu) / to_double(TRANSMITTANCE_INTEGRAL_SAMPLES);
|
2012-12-30 16:06:11 +00:00
|
|
|
double yi = exp(-(r - Rg) / H);
|
|
|
|
int i;
|
|
|
|
for (i = 1; i <= TRANSMITTANCE_INTEGRAL_SAMPLES; ++i) {
|
2015-12-17 00:13:20 +00:00
|
|
|
double xj = to_double(i) * dx;
|
2012-12-30 16:06:11 +00:00
|
|
|
double yj = exp(-(sqrt(r * r + xj * xj + 2.0 * xj * r * mu) - Rg) / H);
|
|
|
|
result += (yi + yj) / 2.0 * dx;
|
|
|
|
yi = yj;
|
|
|
|
}
|
|
|
|
return mu < -sqrt(1.0 - (Rg / r) * (Rg / r)) ? 1e9 : result;
|
|
|
|
}
|
2012-12-13 21:33:42 +00:00
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static void _precomputeTransmittanceTexture() {
|
2012-12-12 13:21:46 +00:00
|
|
|
int x, y;
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
for (x = 0; x < TRANSMITTANCE_W; x++) {
|
|
|
|
for (y = 0; y < TRANSMITTANCE_H; y++) {
|
2012-12-12 13:21:46 +00:00
|
|
|
double r, muS;
|
2015-12-17 00:13:20 +00:00
|
|
|
_getTransmittanceRMu(to_double(x + 0.5) / TRANSMITTANCE_W, to_double(y + 0.5) / TRANSMITTANCE_H, &r, &muS);
|
2012-12-12 13:21:46 +00:00
|
|
|
double depth1 = _opticalDepthTransmittance(HR, r, muS);
|
|
|
|
double depth2 = _opticalDepthTransmittance(HM, r, muS);
|
|
|
|
Color trans;
|
|
|
|
trans.r = exp(-(betaR.r * depth1 + betaMEx.x * depth2));
|
|
|
|
trans.g = exp(-(betaR.g * depth1 + betaMEx.y * depth2));
|
|
|
|
trans.b = exp(-(betaR.b * depth1 + betaMEx.z * depth2));
|
|
|
|
trans.a = 1.0;
|
2013-12-09 21:16:00 +00:00
|
|
|
_transmittanceTexture->setPixel(x, y, trans); /* Eq (5) */
|
2012-12-12 13:21:46 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2012-12-30 16:06:11 +00:00
|
|
|
/*********************** irradiance1.glsl ***********************/
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static void _precomputeIrrDeltaETexture(Texture2D *destination) {
|
2012-12-13 21:33:42 +00:00
|
|
|
int x, y;
|
|
|
|
|
|
|
|
/* Irradiance program */
|
2015-11-09 21:30:46 +00:00
|
|
|
for (x = 0; x < SKY_W; x++) {
|
|
|
|
for (y = 0; y < SKY_H; y++) {
|
2012-12-30 16:06:11 +00:00
|
|
|
double r, muS;
|
2012-12-13 21:33:42 +00:00
|
|
|
Color trans, irr;
|
2015-12-17 00:13:20 +00:00
|
|
|
_getIrradianceRMuS(to_double(x) / SKY_W, to_double(y) / SKY_H, &r, &muS);
|
2012-12-13 21:33:42 +00:00
|
|
|
trans = _transmittance(r, muS);
|
|
|
|
|
|
|
|
irr.r = trans.r * max(muS, 0.0);
|
|
|
|
irr.g = trans.g * max(muS, 0.0);
|
|
|
|
irr.b = trans.b * max(muS, 0.0);
|
|
|
|
irr.a = 1.0;
|
|
|
|
|
2013-12-09 21:16:00 +00:00
|
|
|
destination->setPixel(x, y, irr);
|
2012-12-13 21:33:42 +00:00
|
|
|
}
|
2012-12-12 20:11:57 +00:00
|
|
|
}
|
2012-12-12 13:21:46 +00:00
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static void _getLayerParams(int layer, double *_r, Color *_dhdH) {
|
2012-12-13 21:33:42 +00:00
|
|
|
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);
|
2012-12-12 13:21:46 +00:00
|
|
|
|
2012-12-18 16:20:38 +00:00
|
|
|
*_r = r;
|
|
|
|
_dhdH->r = dmin;
|
|
|
|
_dhdH->g = dmax;
|
|
|
|
_dhdH->b = dminp;
|
|
|
|
_dhdH->a = dmaxp;
|
2012-12-13 21:33:42 +00:00
|
|
|
}
|
2012-12-12 13:21:46 +00:00
|
|
|
|
2012-12-13 21:33:42 +00:00
|
|
|
/*********************** inscatter1.glsl ***********************/
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static void _integrand1(double r, double mu, double muS, double nu, double t, Color *ray, Color *mie) {
|
2012-12-13 21:33:42 +00:00
|
|
|
double ri = sqrt(r * r + t * t + 2.0 * r * mu * t);
|
|
|
|
double muSi = (nu * t + muS * r) / ri;
|
|
|
|
ri = max(Rg, ri);
|
2015-11-09 21:30:46 +00:00
|
|
|
if (muSi >= -sqrt(1.0 - Rg * Rg / (ri * ri))) {
|
2012-12-13 21:33:42 +00:00
|
|
|
Color t1, t2;
|
|
|
|
t1 = _transmittance3(r, mu, t);
|
|
|
|
t2 = _transmittance(ri, muSi);
|
|
|
|
double fR = exp(-(ri - Rg) / HR);
|
|
|
|
double fM = exp(-(ri - Rg) / HM);
|
|
|
|
ray->r = fR * t1.r * t2.r;
|
|
|
|
ray->g = fR * t1.g * t2.g;
|
|
|
|
ray->b = fR * t1.b * t2.b;
|
|
|
|
mie->r = fM * t1.r * t2.r;
|
|
|
|
mie->g = fM * t1.g * t2.g;
|
|
|
|
mie->b = fM * t1.b * t2.b;
|
2015-11-09 21:30:46 +00:00
|
|
|
} else {
|
2012-12-13 21:33:42 +00:00
|
|
|
ray->r = ray->g = ray->b = 0.0;
|
|
|
|
mie->r = mie->g = mie->b = 0.0;
|
|
|
|
}
|
2012-12-12 13:21:46 +00:00
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static void _inscatter1(double r, double mu, double muS, double nu, Color *ray, Color *mie) {
|
2012-12-13 21:33:42 +00:00
|
|
|
ray->r = ray->g = ray->b = 0.0;
|
|
|
|
mie->r = mie->g = mie->b = 0.0;
|
2015-12-17 00:13:20 +00:00
|
|
|
double dx = _limit(r, mu) / to_double(INSCATTER_INTEGRAL_SAMPLES);
|
2012-12-13 21:33:42 +00:00
|
|
|
Color rayi;
|
|
|
|
Color miei;
|
|
|
|
_integrand1(r, mu, muS, nu, 0.0, &rayi, &miei);
|
|
|
|
int i;
|
2015-11-09 21:30:46 +00:00
|
|
|
for (i = 1; i <= INSCATTER_INTEGRAL_SAMPLES; ++i) {
|
2015-12-23 23:16:50 +00:00
|
|
|
double xj = to_double(i) * dx;
|
2012-12-13 21:33:42 +00:00
|
|
|
Color rayj;
|
|
|
|
Color miej;
|
|
|
|
_integrand1(r, mu, muS, nu, xj, &rayj, &miej);
|
|
|
|
ray->r += (rayi.r + rayj.r) / 2.0 * dx;
|
|
|
|
ray->g += (rayi.g + rayj.g) / 2.0 * dx;
|
|
|
|
ray->b += (rayi.b + rayj.b) / 2.0 * dx;
|
|
|
|
mie->r += (miei.r + miej.r) / 2.0 * dx;
|
|
|
|
mie->g += (miei.g + miej.g) / 2.0 * dx;
|
|
|
|
mie->b += (miei.b + miej.b) / 2.0 * dx;
|
|
|
|
rayi = rayj;
|
|
|
|
miei = miej;
|
|
|
|
}
|
|
|
|
ray->r *= betaR.r;
|
|
|
|
ray->g *= betaR.g;
|
|
|
|
ray->b *= betaR.b;
|
|
|
|
mie->r *= betaMSca.x;
|
|
|
|
mie->g *= betaMSca.y;
|
|
|
|
mie->b *= betaMSca.z;
|
2012-12-02 11:08:56 +00:00
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
typedef struct {
|
|
|
|
Texture4D *ray;
|
|
|
|
Texture4D *mie;
|
2012-12-18 16:20:38 +00:00
|
|
|
} Inscatter1Params;
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static int _inscatter1Worker(ParallelWork *, int layer, void *data) {
|
|
|
|
Inscatter1Params *params = (Inscatter1Params *)data;
|
2012-12-18 16:20:38 +00:00
|
|
|
|
|
|
|
double r;
|
|
|
|
Color dhdH;
|
|
|
|
_getLayerParams(layer, &r, &dhdH);
|
|
|
|
|
2013-03-09 11:06:39 +00:00
|
|
|
int x, y, z;
|
2015-11-09 21:30:46 +00:00
|
|
|
for (x = 0; x < RES_MU; x++) {
|
|
|
|
for (y = 0; y < RES_MU_S; y++) {
|
|
|
|
for (z = 0; z < RES_NU; z++) {
|
2013-03-09 11:06:39 +00:00
|
|
|
Color ray = COLOR_BLACK;
|
|
|
|
Color mie = COLOR_BLACK;
|
|
|
|
double mu, muS, nu;
|
2015-12-17 00:13:20 +00:00
|
|
|
_texCoordToMuMuSNu(to_double(x), to_double(y), to_double(z), r, dhdH, &mu, &muS, &nu);
|
2013-03-09 11:06:39 +00:00
|
|
|
_inscatter1(r, mu, muS, nu, &ray, &mie);
|
|
|
|
/* store separately Rayleigh and Mie contributions, WITHOUT the phase function factor
|
|
|
|
* (cf "Angular precision") */
|
2013-12-09 21:16:00 +00:00
|
|
|
params->ray->setPixel(x, y, z, layer, ray);
|
|
|
|
params->mie->setPixel(x, y, z, layer, mie);
|
2013-03-09 11:06:39 +00:00
|
|
|
}
|
2012-12-13 21:33:42 +00:00
|
|
|
}
|
|
|
|
}
|
2012-12-18 16:20:38 +00:00
|
|
|
return 1;
|
2012-12-15 10:14:57 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
/*********************** inscatterS.glsl ***********************/
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static Color _inscatterS(double r, double mu, double muS, double nu, int first, Texture2D *deltaE, Texture4D *deltaSR,
|
|
|
|
Texture4D *deltaSM) {
|
2012-12-15 10:14:57 +00:00
|
|
|
Color raymie = COLOR_BLACK;
|
|
|
|
|
2015-12-30 23:36:22 +00:00
|
|
|
double dphi = Maths::PI / to_double(INSCATTER_SPHERICAL_INTEGRAL_SAMPLES);
|
|
|
|
double dtheta = Maths::PI / to_double(INSCATTER_SPHERICAL_INTEGRAL_SAMPLES);
|
2012-12-15 10:14:57 +00:00
|
|
|
|
|
|
|
r = clamp(r, Rg, Rt);
|
|
|
|
mu = clamp(mu, -1.0, 1.0);
|
|
|
|
muS = clamp(muS, -1.0, 1.0);
|
|
|
|
double var = sqrt(1.0 - mu * mu) * sqrt(1.0 - muS * muS);
|
|
|
|
nu = clamp(nu, muS * mu - var, muS * mu + var);
|
|
|
|
|
|
|
|
double cthetamin = -sqrt(1.0 - (Rg / r) * (Rg / r));
|
|
|
|
|
|
|
|
Vector3 v = vec3(sqrt(1.0 - mu * mu), 0.0, mu);
|
|
|
|
double sx = v.x == 0.0 ? 0.0 : (nu - muS * mu) / v.x;
|
|
|
|
Vector3 s = vec3(sx, sqrt(max(0.0, 1.0 - sx * sx - muS * muS)), muS);
|
|
|
|
|
|
|
|
/* integral over 4.PI around x with two nested loops over w directions (theta,phi) -- Eq (7) */
|
|
|
|
int itheta;
|
2015-11-09 21:30:46 +00:00
|
|
|
for (itheta = 0; itheta < INSCATTER_SPHERICAL_INTEGRAL_SAMPLES; ++itheta) {
|
2015-12-17 00:13:20 +00:00
|
|
|
double theta = (to_double(itheta) + 0.5) * dtheta;
|
2012-12-15 10:14:57 +00:00
|
|
|
double ctheta = cos(theta);
|
|
|
|
|
|
|
|
double greflectance = 0.0;
|
|
|
|
double dground = 0.0;
|
|
|
|
Color gtransp = {0.0, 0.0, 0.0, 0.0};
|
2015-11-09 21:30:46 +00:00
|
|
|
if (ctheta < cthetamin) {
|
2012-12-15 10:14:57 +00:00
|
|
|
/* if ground visible in direction w
|
|
|
|
* compute transparency gtransp between x and ground */
|
2015-12-30 23:36:22 +00:00
|
|
|
greflectance = AVERAGE_GROUND_REFLECTANCE / Maths::PI;
|
2012-12-15 10:14:57 +00:00
|
|
|
dground = -r * ctheta - sqrt(r * r * (ctheta * ctheta - 1.0) + Rg * Rg);
|
|
|
|
gtransp = _transmittance3(Rg, -(r * ctheta + dground) / Rg, dground);
|
|
|
|
}
|
|
|
|
|
|
|
|
int iphi;
|
2015-11-09 21:30:46 +00:00
|
|
|
for (iphi = 0; iphi < 2 * INSCATTER_SPHERICAL_INTEGRAL_SAMPLES; ++iphi) {
|
2015-12-17 00:13:20 +00:00
|
|
|
double phi = (to_double(iphi) + 0.5) * dphi;
|
2012-12-15 10:14:57 +00:00
|
|
|
double dw = dtheta * dphi * sin(theta);
|
|
|
|
Vector3 w = vec3(cos(phi) * sin(theta), sin(phi) * sin(theta), ctheta);
|
|
|
|
|
2013-12-11 10:32:10 +00:00
|
|
|
double nu1 = s.dotProduct(w);
|
|
|
|
double nu2 = v.dotProduct(w);
|
2012-12-15 10:14:57 +00:00
|
|
|
double pr2 = _phaseFunctionR(nu2);
|
|
|
|
double pm2 = _phaseFunctionM(nu2);
|
|
|
|
|
|
|
|
/* compute irradiance received at ground in direction w (if ground visible) =deltaE */
|
|
|
|
Vector3 gnormal;
|
|
|
|
gnormal.x = dground * w.x / Rg;
|
|
|
|
gnormal.y = dground * w.y / Rg;
|
2012-12-30 16:06:11 +00:00
|
|
|
gnormal.z = (r + dground * w.z) / Rg;
|
2013-12-11 10:32:10 +00:00
|
|
|
Color girradiance = _irradiance(deltaE, Rg, gnormal.dotProduct(s));
|
2012-12-15 10:14:57 +00:00
|
|
|
|
|
|
|
Color 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.r = greflectance * girradiance.r * gtransp.r;
|
|
|
|
raymie1.g = greflectance * girradiance.g * gtransp.g;
|
|
|
|
raymie1.b = greflectance * girradiance.b * gtransp.b;
|
|
|
|
|
|
|
|
/* second term = inscattered light, =deltaS */
|
2015-11-09 21:30:46 +00:00
|
|
|
if (first) {
|
2012-12-15 10:14:57 +00:00
|
|
|
/* first iteration is special because Rayleigh and Mie were stored separately,
|
|
|
|
* without the phase functions factors; they must be reintroduced here */
|
|
|
|
double pr1 = _phaseFunctionR(nu1);
|
|
|
|
double pm1 = _phaseFunctionM(nu1);
|
2012-12-24 11:24:27 +00:00
|
|
|
Color ray1 = _texture4D(deltaSR, r, w.z, muS, nu1);
|
|
|
|
Color mie1 = _texture4D(deltaSM, r, w.z, muS, nu1);
|
2013-01-10 21:21:56 +00:00
|
|
|
raymie1.r += ray1.r * pr1 + mie1.r * pm1;
|
|
|
|
raymie1.g += ray1.g * pr1 + mie1.g * pm1;
|
|
|
|
raymie1.b += ray1.b * pr1 + mie1.b * pm1;
|
2015-11-09 21:30:46 +00:00
|
|
|
} else {
|
2012-12-24 11:24:27 +00:00
|
|
|
Color col = _texture4D(deltaSR, r, w.z, muS, nu1);
|
2012-12-15 10:14:57 +00:00
|
|
|
raymie1.r += col.r;
|
|
|
|
raymie1.g += col.g;
|
|
|
|
raymie1.b += col.b;
|
|
|
|
}
|
|
|
|
|
|
|
|
/* 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.r += raymie1.r * (betaR.r * exp(-(r - Rg) / HR) * pr2 + betaMSca.x * exp(-(r - Rg) / HM) * pm2) * dw;
|
|
|
|
raymie.g += raymie1.g * (betaR.g * exp(-(r - Rg) / HR) * pr2 + betaMSca.y * exp(-(r - Rg) / HM) * pm2) * dw;
|
|
|
|
raymie.b += raymie1.b * (betaR.b * exp(-(r - Rg) / HR) * pr2 + betaMSca.z * exp(-(r - Rg) / HM) * pm2) * dw;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
/* output raymie = J[T.alpha/PI.deltaE + deltaS] (line 7 in algorithm 4.1) */
|
|
|
|
return raymie;
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
typedef struct {
|
|
|
|
Texture4D *result;
|
|
|
|
Texture2D *deltaE;
|
|
|
|
Texture4D *deltaSR;
|
|
|
|
Texture4D *deltaSM;
|
2012-12-18 16:20:38 +00:00
|
|
|
int first;
|
|
|
|
} jParams;
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static int _jWorker(ParallelWork *, int layer, void *data) {
|
|
|
|
jParams *params = (jParams *)data;
|
2012-12-18 16:20:38 +00:00
|
|
|
|
|
|
|
double r;
|
|
|
|
Color dhdH;
|
|
|
|
_getLayerParams(layer, &r, &dhdH);
|
|
|
|
|
2013-03-09 11:06:39 +00:00
|
|
|
int x, y, z;
|
2015-11-09 21:30:46 +00:00
|
|
|
for (x = 0; x < RES_MU; x++) {
|
|
|
|
for (y = 0; y < RES_MU_S; y++) {
|
|
|
|
for (z = 0; z < RES_NU; z++) {
|
2013-03-09 11:06:39 +00:00
|
|
|
Color raymie;
|
|
|
|
double mu, muS, nu;
|
2015-12-17 00:13:20 +00:00
|
|
|
_texCoordToMuMuSNu(to_double(x), to_double(y), to_double(z), r, dhdH, &mu, &muS, &nu);
|
2013-03-09 11:06:39 +00:00
|
|
|
raymie = _inscatterS(r, mu, muS, nu, params->first, params->deltaE, params->deltaSR, params->deltaSM);
|
2013-12-09 21:16:00 +00:00
|
|
|
params->result->setPixel(x, y, z, layer, raymie);
|
2013-03-09 11:06:39 +00:00
|
|
|
}
|
2012-12-15 10:14:57 +00:00
|
|
|
}
|
|
|
|
}
|
2012-12-18 16:20:38 +00:00
|
|
|
return 1;
|
2012-12-15 10:14:57 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
/*********************** irradianceN.glsl ***********************/
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
void _irradianceNProg(Texture2D *destination, Texture4D *deltaSR, Texture4D *deltaSM, int first) {
|
2012-12-15 10:14:57 +00:00
|
|
|
int x, y;
|
2015-12-30 23:36:22 +00:00
|
|
|
double dphi = Maths::PI / to_double(IRRADIANCE_INTEGRAL_SAMPLES);
|
|
|
|
double dtheta = Maths::PI / to_double(IRRADIANCE_INTEGRAL_SAMPLES);
|
2015-11-09 21:30:46 +00:00
|
|
|
for (x = 0; x < SKY_W; x++) {
|
|
|
|
for (y = 0; y < SKY_H; y++) {
|
2012-12-15 10:14:57 +00:00
|
|
|
double r, muS;
|
|
|
|
int iphi;
|
2015-12-17 00:13:20 +00:00
|
|
|
_getIrradianceRMuS(to_double(x) / SKY_W, to_double(y) / SKY_H, &r, &muS);
|
2012-12-15 10:14:57 +00:00
|
|
|
Vector3 s = vec3(max(sqrt(1.0 - muS * muS), 0.0), 0.0, muS);
|
|
|
|
|
|
|
|
Color result = COLOR_BLACK;
|
|
|
|
/* integral over 2.PI around x with two nested loops over w directions (theta,phi) -- Eq (15) */
|
2015-11-09 21:30:46 +00:00
|
|
|
for (iphi = 0; iphi < 2 * IRRADIANCE_INTEGRAL_SAMPLES; ++iphi) {
|
2015-12-17 00:13:20 +00:00
|
|
|
double phi = (to_double(iphi) + 0.5) * dphi;
|
2012-12-15 10:14:57 +00:00
|
|
|
int itheta;
|
2015-11-09 21:30:46 +00:00
|
|
|
for (itheta = 0; itheta < IRRADIANCE_INTEGRAL_SAMPLES / 2; ++itheta) {
|
2015-12-17 00:13:20 +00:00
|
|
|
double theta = (to_double(itheta) + 0.5) * dtheta;
|
2012-12-15 10:14:57 +00:00
|
|
|
double dw = dtheta * dphi * sin(theta);
|
|
|
|
Vector3 w = vec3(cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta));
|
2013-12-11 10:32:10 +00:00
|
|
|
double nu = s.dotProduct(w);
|
2015-11-09 21:30:46 +00:00
|
|
|
if (first) {
|
2012-12-15 10:14:57 +00:00
|
|
|
/* first iteration is special because Rayleigh and Mie were stored separately,
|
|
|
|
without the phase functions factors; they must be reintroduced here */
|
|
|
|
double pr1 = _phaseFunctionR(nu);
|
|
|
|
double pm1 = _phaseFunctionM(nu);
|
2012-12-24 11:24:27 +00:00
|
|
|
Color ray1 = _texture4D(deltaSR, r, w.z, muS, nu);
|
|
|
|
Color mie1 = _texture4D(deltaSM, r, w.z, muS, nu);
|
2012-12-15 10:14:57 +00:00
|
|
|
result.r += (ray1.r * pr1 + mie1.r * pm1) * w.z * dw;
|
|
|
|
result.g += (ray1.g * pr1 + mie1.g * pm1) * w.z * dw;
|
|
|
|
result.b += (ray1.b * pr1 + mie1.b * pm1) * w.z * dw;
|
2015-11-09 21:30:46 +00:00
|
|
|
} else {
|
2012-12-24 11:24:27 +00:00
|
|
|
Color col = _texture4D(deltaSR, r, w.z, muS, nu);
|
2012-12-15 10:14:57 +00:00
|
|
|
result.r += col.r * w.z * dw;
|
|
|
|
result.g += col.g * w.z * dw;
|
|
|
|
result.b += col.b * w.z * dw;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2013-12-09 21:16:00 +00:00
|
|
|
destination->setPixel(x, y, result);
|
2012-12-15 10:14:57 +00:00
|
|
|
}
|
|
|
|
}
|
2012-12-13 21:33:42 +00:00
|
|
|
}
|
|
|
|
|
2012-12-15 19:45:19 +00:00
|
|
|
/*********************** inscatterN.glsl ***********************/
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
typedef struct {
|
|
|
|
Texture4D *destination;
|
|
|
|
Texture4D *deltaJ;
|
2012-12-24 11:24:27 +00:00
|
|
|
} InscatterNParams;
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static Color _integrand2(Texture4D *deltaJ, double r, double mu, double muS, double nu, double t) {
|
2012-12-15 19:45:19 +00:00
|
|
|
double ri = sqrt(r * r + t * t + 2.0 * r * mu * t);
|
|
|
|
double mui = (r * mu + t) / ri;
|
|
|
|
double muSi = (nu * t + muS * r) / ri;
|
|
|
|
Color c1, c2;
|
2012-12-24 11:24:27 +00:00
|
|
|
c1 = _texture4D(deltaJ, ri, mui, muSi, nu);
|
2012-12-15 19:45:19 +00:00
|
|
|
c2 = _transmittance3(r, mu, t);
|
|
|
|
c1.r *= c2.r;
|
|
|
|
c1.g *= c2.g;
|
|
|
|
c1.b *= c2.b;
|
2012-12-30 16:06:11 +00:00
|
|
|
c1.a = 1.0;
|
2012-12-15 19:45:19 +00:00
|
|
|
return c1;
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static Color _inscatterN(Texture4D *deltaJ, double r, double mu, double muS, double nu) {
|
2012-12-15 19:45:19 +00:00
|
|
|
Color raymie = COLOR_BLACK;
|
2015-12-17 00:13:20 +00:00
|
|
|
double dx = _limit(r, mu) / to_double(INSCATTER_INTEGRAL_SAMPLES);
|
2012-12-24 11:24:27 +00:00
|
|
|
Color raymiei = _integrand2(deltaJ, r, mu, muS, nu, 0.0);
|
2012-12-15 19:45:19 +00:00
|
|
|
int i;
|
2015-11-09 21:30:46 +00:00
|
|
|
for (i = 1; i <= INSCATTER_INTEGRAL_SAMPLES; ++i) {
|
2015-12-23 23:16:50 +00:00
|
|
|
double xj = to_double(i) * dx;
|
2012-12-24 11:24:27 +00:00
|
|
|
Color raymiej = _integrand2(deltaJ, r, mu, muS, nu, xj);
|
2012-12-15 19:45:19 +00:00
|
|
|
raymie.r += (raymiei.r + raymiej.r) / 2.0 * dx;
|
2012-12-30 16:06:11 +00:00
|
|
|
raymie.g += (raymiei.g + raymiej.g) / 2.0 * dx;
|
|
|
|
raymie.b += (raymiei.b + raymiej.b) / 2.0 * dx;
|
2012-12-15 19:45:19 +00:00
|
|
|
raymiei = raymiej;
|
|
|
|
}
|
|
|
|
return raymie;
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static int _inscatterNWorker(ParallelWork *, int layer, void *data) {
|
|
|
|
InscatterNParams *params = (InscatterNParams *)data;
|
2012-12-24 11:24:27 +00:00
|
|
|
|
2012-12-18 16:20:38 +00:00
|
|
|
double r;
|
|
|
|
Color dhdH;
|
|
|
|
_getLayerParams(layer, &r, &dhdH);
|
|
|
|
|
2013-03-09 11:06:39 +00:00
|
|
|
int x, y, z;
|
2015-11-09 21:30:46 +00:00
|
|
|
for (x = 0; x < RES_MU; x++) {
|
|
|
|
for (y = 0; y < RES_MU_S; y++) {
|
|
|
|
for (z = 0; z < RES_NU; z++) {
|
2013-03-09 11:06:39 +00:00
|
|
|
double mu, muS, nu;
|
2015-12-17 00:13:20 +00:00
|
|
|
_texCoordToMuMuSNu(to_double(x), to_double(y), to_double(z), r, dhdH, &mu, &muS, &nu);
|
2013-12-09 21:16:00 +00:00
|
|
|
params->destination->setPixel(x, y, z, layer, _inscatterN(params->deltaJ, r, mu, muS, nu));
|
2013-03-09 11:06:39 +00:00
|
|
|
}
|
2012-12-15 19:45:19 +00:00
|
|
|
}
|
|
|
|
}
|
2012-12-18 16:20:38 +00:00
|
|
|
return 1;
|
2012-12-15 19:45:19 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
/*********************** copyInscatterN.glsl ***********************/
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
typedef struct {
|
|
|
|
Texture4D *source;
|
|
|
|
Texture4D *destination;
|
2012-12-18 16:20:38 +00:00
|
|
|
} CopyInscatterNParams;
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static int _copyInscatterNWorker(ParallelWork *, int layer, void *data) {
|
|
|
|
CopyInscatterNParams *params = (CopyInscatterNParams *)data;
|
2012-12-18 16:20:38 +00:00
|
|
|
|
|
|
|
double r;
|
|
|
|
Color dhdH;
|
|
|
|
_getLayerParams(layer, &r, &dhdH);
|
|
|
|
|
2013-03-09 11:06:39 +00:00
|
|
|
int x, y, z;
|
2015-11-09 21:30:46 +00:00
|
|
|
for (x = 0; x < RES_MU; x++) {
|
|
|
|
for (y = 0; y < RES_MU_S; y++) {
|
|
|
|
for (z = 0; z < RES_NU; z++) {
|
2013-03-09 11:06:39 +00:00
|
|
|
double mu, muS, nu;
|
2015-12-17 00:13:20 +00:00
|
|
|
_texCoordToMuMuSNu(to_double(x), to_double(y), to_double(z), r, dhdH, &mu, &muS, &nu);
|
2013-12-09 21:16:00 +00:00
|
|
|
Color col1 = params->source->getPixel(x, y, z, layer);
|
|
|
|
Color col2 = params->destination->getPixel(x, y, z, layer);
|
2013-03-09 11:06:39 +00:00
|
|
|
col2.r += col1.r / _phaseFunctionR(nu);
|
|
|
|
col2.g += col1.g / _phaseFunctionR(nu);
|
|
|
|
col2.b += col1.b / _phaseFunctionR(nu);
|
2013-12-09 21:16:00 +00:00
|
|
|
params->destination->setPixel(x, y, z, layer, col2);
|
2013-03-09 11:06:39 +00:00
|
|
|
}
|
2012-12-15 19:45:19 +00:00
|
|
|
}
|
|
|
|
}
|
2012-12-18 16:20:38 +00:00
|
|
|
return 1;
|
2012-12-15 19:45:19 +00:00
|
|
|
}
|
|
|
|
|
2012-12-13 21:33:42 +00:00
|
|
|
/*********************** Final getters ***********************/
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static inline Color _applyInscatter(Color inscatter, Color attmod, Color samp) {
|
2012-12-23 22:05:12 +00:00
|
|
|
inscatter.r = inscatter.r - attmod.r * samp.r;
|
|
|
|
inscatter.g = inscatter.g - attmod.g * samp.g;
|
|
|
|
inscatter.b = inscatter.b - attmod.b * samp.b;
|
|
|
|
inscatter.a = inscatter.a - attmod.a * samp.a;
|
|
|
|
return vec4max(inscatter, 0.0);
|
|
|
|
}
|
|
|
|
|
2012-12-02 11:08:56 +00:00
|
|
|
/* inscattered light along ray x+tv, when sun in direction s (=S[L]-T(x,x0)S[L]|x0) */
|
2015-11-09 21:30:46 +00:00
|
|
|
static Color _getInscatterColor(Vector3 *_x, double *_t, Vector3 v, Vector3 s, double *_r, double *_mu,
|
|
|
|
Vector3 *attenuation) {
|
2012-12-02 11:08:56 +00:00
|
|
|
Color result;
|
2013-12-11 10:32:10 +00:00
|
|
|
double r = _x->getNorm();
|
|
|
|
double mu = _x->dotProduct(v) / r;
|
2012-12-02 11:08:56 +00:00
|
|
|
double d = -r * mu - sqrt(r * r * (mu * mu - 1.0) + Rt * Rt);
|
2013-01-22 20:50:37 +00:00
|
|
|
attenuation->x = attenuation->y = attenuation->z = 0.0;
|
2015-11-09 21:30:46 +00:00
|
|
|
if (d > 0.0) {
|
2012-12-12 13:21:46 +00:00
|
|
|
/* if x in space and ray intersects atmosphere
|
|
|
|
move x to nearest intersection of ray with top atmosphere boundary */
|
|
|
|
_x->x += d * v.x;
|
|
|
|
_x->y += d * v.y;
|
|
|
|
_x->z += d * v.z;
|
|
|
|
*_t -= d;
|
2012-12-02 11:08:56 +00:00
|
|
|
mu = (r * mu + d) / Rt;
|
|
|
|
r = Rt;
|
|
|
|
}
|
2012-12-12 13:21:46 +00:00
|
|
|
double t = *_t;
|
|
|
|
Vector3 x = *_x;
|
2015-11-09 21:30:46 +00:00
|
|
|
if (r <= Rt) {
|
2012-12-12 13:21:46 +00:00
|
|
|
/* if ray intersects atmosphere */
|
2013-12-11 10:32:10 +00:00
|
|
|
double nu = v.dotProduct(s);
|
|
|
|
double muS = x.dotProduct(s) / r;
|
2012-12-02 11:08:56 +00:00
|
|
|
double phaseR = _phaseFunctionR(nu);
|
|
|
|
double phaseM = _phaseFunctionM(nu);
|
2012-12-13 21:33:42 +00:00
|
|
|
Color inscatter = vec4max(_texture4D(_inscatterTexture, r, mu, muS, nu), 0.0);
|
2015-11-09 21:30:46 +00:00
|
|
|
if (t > 0.0) {
|
2013-12-11 10:32:10 +00:00
|
|
|
Vector3 x0 = x.add(v.scale(t));
|
|
|
|
double r0 = x0.getNorm();
|
|
|
|
double rMu0 = x0.dotProduct(v);
|
2012-12-02 11:08:56 +00:00
|
|
|
double mu0 = rMu0 / r0;
|
2013-12-11 10:32:10 +00:00
|
|
|
double muS0 = x0.dotProduct(s) / r0;
|
2012-12-12 13:21:46 +00:00
|
|
|
/* avoids imprecision problems in transmittance computations based on textures */
|
2012-12-02 11:08:56 +00:00
|
|
|
*attenuation = _analyticTransmittance(r, mu, t);
|
2015-12-21 16:48:30 +00:00
|
|
|
if (r0 > Rg + 0.01) {
|
2012-12-12 13:21:46 +00:00
|
|
|
/* computes S[L]-T(x,x0)S[L]|x0 */
|
|
|
|
Color attmod = {attenuation->x, attenuation->y, attenuation->z, attenuation->x};
|
2012-12-13 21:33:42 +00:00
|
|
|
Color samp = _texture4D(_inscatterTexture, r0, mu0, muS0, nu);
|
2012-12-12 13:21:46 +00:00
|
|
|
inscatter = _applyInscatter(inscatter, attmod, samp);
|
2015-11-09 21:30:46 +00:00
|
|
|
/* avoids imprecision problems near horizon by interpolating between two points above and below horizon
|
|
|
|
*/
|
2015-12-21 16:48:30 +00:00
|
|
|
const double EPS = 0.004;
|
2012-12-12 13:21:46 +00:00
|
|
|
double muHoriz = -sqrt(1.0 - (Rg / r) * (Rg / r));
|
2015-11-09 21:30:46 +00:00
|
|
|
if (fabs(mu - muHoriz) < EPS) {
|
2012-12-12 13:21:46 +00:00
|
|
|
double a = ((mu - muHoriz) + EPS) / (2.0 * EPS);
|
2012-12-02 11:08:56 +00:00
|
|
|
|
|
|
|
mu = muHoriz - EPS;
|
|
|
|
r0 = sqrt(r * r + t * t + 2.0 * r * t * mu);
|
|
|
|
mu0 = (r * mu + t) / r0;
|
2012-12-13 21:33:42 +00:00
|
|
|
Color inScatter0 = _texture4D(_inscatterTexture, r, mu, muS, nu);
|
|
|
|
Color inScatter1 = _texture4D(_inscatterTexture, r0, mu0, muS0, nu);
|
2012-12-12 13:21:46 +00:00
|
|
|
Color inScatterA = _applyInscatter(inScatter0, attmod, inScatter1);
|
2012-12-02 11:08:56 +00:00
|
|
|
|
|
|
|
mu = muHoriz + EPS;
|
|
|
|
r0 = sqrt(r * r + t * t + 2.0 * r * t * mu);
|
|
|
|
mu0 = (r * mu + t) / r0;
|
2012-12-13 21:33:42 +00:00
|
|
|
inScatter0 = _texture4D(_inscatterTexture, r, mu, muS, nu);
|
|
|
|
inScatter1 = _texture4D(_inscatterTexture, r0, mu0, muS0, nu);
|
2012-12-12 13:21:46 +00:00
|
|
|
Color inScatterB = _applyInscatter(inScatter0, attmod, inScatter1);
|
2012-12-02 11:08:56 +00:00
|
|
|
|
2012-12-12 13:21:46 +00:00
|
|
|
inscatter = vec4mix(inScatterA, inScatterB, a);
|
2012-12-02 11:08:56 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2012-12-12 13:21:46 +00:00
|
|
|
/* avoids imprecision problems in Mie scattering when sun is below horizon */
|
|
|
|
inscatter.a *= smoothstep(0.00, 0.02, muS);
|
2012-12-13 21:33:42 +00:00
|
|
|
Color mie = _getMie(inscatter);
|
2012-12-12 13:21:46 +00:00
|
|
|
result.r = inscatter.r * phaseR + mie.r * phaseM;
|
|
|
|
result.g = inscatter.g * phaseR + mie.g * phaseM;
|
|
|
|
result.b = inscatter.b * phaseR + mie.b * phaseM;
|
2012-12-23 22:05:12 +00:00
|
|
|
result.a = 1.0;
|
2012-12-12 13:21:46 +00:00
|
|
|
_fixVec4Min(&result, 0.0);
|
2015-11-09 21:30:46 +00:00
|
|
|
} else {
|
2012-12-12 13:21:46 +00:00
|
|
|
/* x in space and ray looking in space */
|
2012-12-02 11:08:56 +00:00
|
|
|
result = COLOR_BLACK;
|
|
|
|
}
|
|
|
|
|
|
|
|
*_r = r;
|
|
|
|
*_mu = mu;
|
|
|
|
result.r *= ISun;
|
|
|
|
result.g *= ISun;
|
|
|
|
result.b *= ISun;
|
|
|
|
result.a = 1.0;
|
|
|
|
return result;
|
|
|
|
}
|
|
|
|
|
|
|
|
/* direct sun light for ray x+tv, when sun in direction s (=L0) */
|
2015-11-09 21:30:46 +00:00
|
|
|
static Color _sunColor(Vector3 v, Vector3 s, double r, double mu, double radius) {
|
2013-01-16 15:26:33 +00:00
|
|
|
Color transmittance = r <= Rt ? _transmittanceWithShadow(r, mu) : COLOR_WHITE; /* T(x,xo) */
|
2016-01-15 17:01:01 +00:00
|
|
|
double isun = step(cos(radius), v.dotProduct(s)) * ISun; /* Lsun */
|
2013-01-16 15:26:33 +00:00
|
|
|
transmittance.r *= isun;
|
|
|
|
transmittance.g *= isun;
|
|
|
|
transmittance.b *= isun;
|
|
|
|
transmittance.a = 1.0;
|
|
|
|
return transmittance; /* Eq (9) */
|
2012-12-02 11:08:56 +00:00
|
|
|
}
|
|
|
|
|
2012-12-24 11:24:27 +00:00
|
|
|
/*********************** Cache/debug methods ***********************/
|
2012-12-14 16:16:09 +00:00
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static int _tryLoadCache2D(Texture2D *tex, const char *tag, int order) {
|
2012-12-14 16:16:09 +00:00
|
|
|
int xsize, ysize;
|
2013-12-09 21:16:00 +00:00
|
|
|
tex->getSize(&xsize, &ysize);
|
2012-12-14 16:16:09 +00:00
|
|
|
|
2013-12-09 21:16:00 +00:00
|
|
|
CacheFile cache("atmo-br", "cache", tag, xsize, ysize, 0, 0, order);
|
2015-11-09 21:30:46 +00:00
|
|
|
if (cache.isReadable()) {
|
2013-11-03 12:00:31 +00:00
|
|
|
PackStream stream;
|
2013-12-09 21:16:00 +00:00
|
|
|
stream.bindToFile(cache.getPath());
|
|
|
|
tex->load(&stream);
|
2012-12-24 11:24:27 +00:00
|
|
|
|
2012-12-14 16:16:09 +00:00
|
|
|
return 1;
|
2015-11-09 21:30:46 +00:00
|
|
|
} else {
|
2012-12-14 16:16:09 +00:00
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static void _saveCache2D(Texture2D *tex, const char *tag, int order) {
|
2012-12-14 16:16:09 +00:00
|
|
|
int xsize, ysize;
|
2013-12-09 21:16:00 +00:00
|
|
|
tex->getSize(&xsize, &ysize);
|
2012-12-14 16:16:09 +00:00
|
|
|
|
2013-12-09 21:16:00 +00:00
|
|
|
CacheFile cache("atmo-br", "cache", tag, xsize, ysize, 0, 0, order);
|
2015-11-09 21:30:46 +00:00
|
|
|
if (cache.isWritable()) {
|
2013-11-03 12:00:31 +00:00
|
|
|
PackStream stream;
|
2015-12-21 16:48:30 +00:00
|
|
|
stream.bindToFile(cache.getPath(), true);
|
2013-12-09 21:16:00 +00:00
|
|
|
tex->save(&stream);
|
2012-12-24 11:24:27 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static void _saveDebug2D(Texture2D *tex, const char *tag, int order) {
|
2012-12-24 11:24:27 +00:00
|
|
|
int xsize, ysize;
|
2013-12-09 21:16:00 +00:00
|
|
|
tex->getSize(&xsize, &ysize);
|
2012-12-24 11:24:27 +00:00
|
|
|
|
2013-12-09 21:16:00 +00:00
|
|
|
CacheFile cache("atmo-br", "png", tag, xsize, ysize, 0, 0, order);
|
2015-11-09 21:30:46 +00:00
|
|
|
if (cache.isWritable()) {
|
2013-12-09 21:16:00 +00:00
|
|
|
tex->saveToFile(cache.getPath());
|
2012-12-14 16:16:09 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static int _tryLoadCache4D(Texture4D *tex, const char *tag, int order) {
|
2013-03-09 09:54:31 +00:00
|
|
|
int xsize, ysize, zsize, wsize;
|
2013-12-09 21:16:00 +00:00
|
|
|
tex->getSize(&xsize, &ysize, &zsize, &wsize);
|
2012-12-14 16:16:09 +00:00
|
|
|
|
2013-12-09 21:16:00 +00:00
|
|
|
CacheFile cache("atmo-br", "cache", tag, xsize, ysize, zsize, wsize, order);
|
2015-11-09 21:30:46 +00:00
|
|
|
if (cache.isReadable()) {
|
2013-11-03 12:00:31 +00:00
|
|
|
PackStream stream;
|
2013-12-09 21:16:00 +00:00
|
|
|
stream.bindToFile(cache.getPath());
|
|
|
|
tex->load(&stream);
|
2012-12-24 11:24:27 +00:00
|
|
|
|
2012-12-14 16:16:09 +00:00
|
|
|
return 1;
|
2015-11-09 21:30:46 +00:00
|
|
|
} else {
|
2012-12-14 16:16:09 +00:00
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static void _saveCache4D(Texture4D *tex, const char *tag, int order) {
|
2013-03-09 09:54:31 +00:00
|
|
|
int xsize, ysize, zsize, wsize;
|
2013-12-09 21:16:00 +00:00
|
|
|
tex->getSize(&xsize, &ysize, &zsize, &wsize);
|
2012-12-14 16:16:09 +00:00
|
|
|
|
2013-12-09 21:16:00 +00:00
|
|
|
CacheFile cache("atmo-br", "cache", tag, xsize, ysize, zsize, wsize, order);
|
2015-11-09 21:30:46 +00:00
|
|
|
if (cache.isWritable()) {
|
2013-11-03 12:00:31 +00:00
|
|
|
PackStream stream;
|
2015-12-21 16:48:30 +00:00
|
|
|
stream.bindToFile(cache.getPath(), true);
|
2013-12-09 21:16:00 +00:00
|
|
|
tex->save(&stream);
|
2012-12-24 11:24:27 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
static void _saveDebug4D(Texture4D *tex, const char *tag, int order) {
|
2013-03-09 09:54:31 +00:00
|
|
|
int xsize, ysize, zsize, wsize;
|
2013-12-09 21:16:00 +00:00
|
|
|
tex->getSize(&xsize, &ysize, &zsize, &wsize);
|
2012-12-24 11:24:27 +00:00
|
|
|
|
2013-12-09 21:16:00 +00:00
|
|
|
CacheFile cache("atmo-br", "png", tag, xsize, ysize, zsize, wsize, order);
|
2015-11-09 21:30:46 +00:00
|
|
|
if (cache.isWritable()) {
|
2013-12-09 21:16:00 +00:00
|
|
|
tex->saveToFile(cache.getPath());
|
2012-12-14 16:16:09 +00:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2012-12-13 21:33:42 +00:00
|
|
|
/*********************** Public methods ***********************/
|
2015-11-09 21:30:46 +00:00
|
|
|
int brunetonInit() {
|
2013-03-09 09:54:31 +00:00
|
|
|
int x, y, z, w, order;
|
2015-11-09 21:30:46 +00:00
|
|
|
ParallelWork *work;
|
2012-12-02 11:08:56 +00:00
|
|
|
|
2012-12-24 11:24:27 +00:00
|
|
|
assert(_inscatterTexture == NULL);
|
2012-12-15 19:45:19 +00:00
|
|
|
|
2012-12-24 11:24:27 +00:00
|
|
|
/* TODO Deletes */
|
2013-12-09 21:16:00 +00:00
|
|
|
_transmittanceTexture = new Texture2D(TRANSMITTANCE_W, TRANSMITTANCE_H);
|
|
|
|
_irradianceTexture = new Texture2D(SKY_W, SKY_H);
|
|
|
|
_inscatterTexture = new Texture4D(RES_MU, RES_MU_S, RES_NU, RES_R);
|
2012-12-24 11:24:27 +00:00
|
|
|
|
|
|
|
/* try loading from cache */
|
2015-11-09 21:30:46 +00:00
|
|
|
if (_tryLoadCache2D(_transmittanceTexture, "transmittance", 0) &&
|
|
|
|
_tryLoadCache2D(_irradianceTexture, "irradiance", 0) && _tryLoadCache4D(_inscatterTexture, "inscatter", 0)) {
|
2013-11-15 23:27:40 +00:00
|
|
|
return 1;
|
2012-12-15 19:45:19 +00:00
|
|
|
}
|
2012-12-12 13:21:46 +00:00
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
Texture2D *_deltaETexture = new Texture2D(SKY_W, SKY_H);
|
|
|
|
Texture4D *_deltaSMTexture = new Texture4D(RES_MU, RES_MU_S, RES_NU, RES_R);
|
|
|
|
Texture4D *_deltaSRTexture = new Texture4D(RES_MU, RES_MU_S, RES_NU, RES_R);
|
|
|
|
Texture4D *_deltaJTexture = new Texture4D(RES_MU, RES_MU_S, RES_NU, RES_R);
|
2012-12-24 11:24:27 +00:00
|
|
|
|
|
|
|
/* computes transmittance texture T (line 1 in algorithm 4.1) */
|
|
|
|
_precomputeTransmittanceTexture();
|
|
|
|
_saveDebug2D(_transmittanceTexture, "transmittance", 0);
|
|
|
|
|
2012-12-15 19:45:19 +00:00
|
|
|
/* computes irradiance texture deltaE (line 2 in algorithm 4.1) */
|
2012-12-24 11:24:27 +00:00
|
|
|
_precomputeIrrDeltaETexture(_deltaETexture);
|
|
|
|
_saveDebug2D(_deltaETexture, "deltaE", 0);
|
2012-12-15 19:45:19 +00:00
|
|
|
|
|
|
|
/* computes single scattering texture deltaS (line 3 in algorithm 4.1)
|
|
|
|
* Rayleigh and Mie separated in deltaSR + deltaSM */
|
2012-12-24 11:24:27 +00:00
|
|
|
Inscatter1Params params = {_deltaSRTexture, _deltaSMTexture};
|
2013-12-09 10:59:57 +00:00
|
|
|
work = new ParallelWork(_inscatter1Worker, RES_R, ¶ms);
|
|
|
|
work->perform();
|
|
|
|
delete work;
|
2013-03-09 09:54:31 +00:00
|
|
|
_saveDebug4D(_deltaSRTexture, "deltaSR", 0);
|
|
|
|
_saveDebug4D(_deltaSMTexture, "deltaSM", 0);
|
2012-12-02 11:08:56 +00:00
|
|
|
|
2012-12-15 19:45:19 +00:00
|
|
|
/* copies deltaE into irradiance texture E (line 4 in algorithm 4.1) */
|
|
|
|
/* ??? all black texture (k=0.0) ??? */
|
2013-12-09 21:16:00 +00:00
|
|
|
_irradianceTexture->fill(COLOR_BLACK);
|
2012-12-15 19:45:19 +00:00
|
|
|
|
|
|
|
/* copies deltaS into inscatter texture S (line 5 in algorithm 4.1) */
|
2015-11-09 21:30:46 +00:00
|
|
|
for (x = 0; x < RES_MU; x++) {
|
|
|
|
for (y = 0; y < RES_MU_S; y++) {
|
|
|
|
for (z = 0; z < RES_NU; z++) {
|
|
|
|
for (w = 0; w < RES_R; w++) {
|
2013-12-09 21:16:00 +00:00
|
|
|
Color result = _deltaSRTexture->getPixel(x, y, z, w);
|
|
|
|
Color mie = _deltaSMTexture->getPixel(x, y, z, w);
|
2013-03-09 09:54:31 +00:00
|
|
|
result.a = mie.r;
|
2013-12-09 21:16:00 +00:00
|
|
|
_inscatterTexture->setPixel(x, y, z, w, result);
|
2013-03-09 09:54:31 +00:00
|
|
|
}
|
2012-12-15 19:45:19 +00:00
|
|
|
}
|
2012-12-14 16:16:09 +00:00
|
|
|
}
|
2012-12-15 19:45:19 +00:00
|
|
|
}
|
2013-03-09 09:54:31 +00:00
|
|
|
_saveDebug4D(_inscatterTexture, "inscatter", 0);
|
2012-12-13 21:33:42 +00:00
|
|
|
|
2012-12-15 19:45:19 +00:00
|
|
|
/* loop for each scattering order (line 6 in algorithm 4.1) */
|
2015-11-09 21:30:46 +00:00
|
|
|
for (order = 2; order <= 4; ++order) {
|
2012-12-15 19:45:19 +00:00
|
|
|
/* computes deltaJ (line 7 in algorithm 4.1) */
|
2012-12-24 11:24:27 +00:00
|
|
|
jParams jparams = {_deltaJTexture, _deltaETexture, _deltaSRTexture, _deltaSMTexture, order == 2};
|
2013-12-09 10:59:57 +00:00
|
|
|
work = new ParallelWork(_jWorker, RES_R, &jparams);
|
|
|
|
work->perform();
|
|
|
|
delete work;
|
2013-03-09 09:54:31 +00:00
|
|
|
_saveDebug4D(_deltaJTexture, "deltaJ", order);
|
2012-12-13 21:33:42 +00:00
|
|
|
|
2012-12-15 19:45:19 +00:00
|
|
|
/* computes deltaE (line 8 in algorithm 4.1) */
|
2012-12-24 11:24:27 +00:00
|
|
|
_irradianceNProg(_deltaETexture, _deltaSRTexture, _deltaSMTexture, order == 2);
|
|
|
|
_saveDebug2D(_deltaETexture, "deltaE", order);
|
2012-12-13 21:33:42 +00:00
|
|
|
|
2012-12-15 19:45:19 +00:00
|
|
|
/* computes deltaS (line 9 in algorithm 4.1) */
|
2012-12-24 11:24:27 +00:00
|
|
|
InscatterNParams iparams = {_deltaSRTexture, _deltaJTexture};
|
2013-12-09 10:59:57 +00:00
|
|
|
work = new ParallelWork(_inscatterNWorker, RES_R, &iparams);
|
|
|
|
work->perform();
|
|
|
|
delete work;
|
2012-12-13 21:33:42 +00:00
|
|
|
|
2012-12-15 19:45:19 +00:00
|
|
|
/* adds deltaE into irradiance texture E (line 10 in algorithm 4.1) */
|
2013-12-09 21:16:00 +00:00
|
|
|
_irradianceTexture->add(_deltaETexture);
|
2012-12-24 11:24:27 +00:00
|
|
|
_saveDebug2D(_irradianceTexture, "irradiance", order);
|
2012-12-13 21:33:42 +00:00
|
|
|
|
2012-12-15 19:45:19 +00:00
|
|
|
/* adds deltaS into inscatter texture S (line 11 in algorithm 4.1) */
|
2012-12-24 11:24:27 +00:00
|
|
|
CopyInscatterNParams cparams = {_deltaSRTexture, _inscatterTexture};
|
2013-12-09 10:59:57 +00:00
|
|
|
work = new ParallelWork(_copyInscatterNWorker, RES_R, &cparams);
|
|
|
|
work->perform();
|
|
|
|
delete work;
|
2013-03-09 09:54:31 +00:00
|
|
|
_saveDebug4D(_inscatterTexture, "inscatter", order);
|
2012-12-13 21:33:42 +00:00
|
|
|
}
|
2012-12-24 11:24:27 +00:00
|
|
|
|
|
|
|
_saveCache2D(_transmittanceTexture, "transmittance", 0);
|
|
|
|
_saveCache2D(_irradianceTexture, "irradiance", 0);
|
2013-03-09 09:54:31 +00:00
|
|
|
_saveCache4D(_inscatterTexture, "inscatter", 0);
|
2012-12-24 11:24:27 +00:00
|
|
|
|
2013-12-09 21:16:00 +00:00
|
|
|
delete _deltaETexture;
|
|
|
|
delete _deltaSMTexture;
|
|
|
|
delete _deltaSRTexture;
|
|
|
|
delete _deltaJTexture;
|
2013-11-15 23:27:40 +00:00
|
|
|
|
|
|
|
return 1;
|
2012-12-02 11:08:56 +00:00
|
|
|
}
|
|
|
|
|
2015-10-18 23:39:22 +00:00
|
|
|
static int _inited;
|
2013-11-15 23:27:40 +00:00
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
AtmosphereModelBruneton::AtmosphereModelBruneton(SoftwareRenderer *parent) : parent(parent) {
|
|
|
|
if (not _inited) {
|
2015-10-18 23:39:22 +00:00
|
|
|
_inited = brunetonInit();
|
|
|
|
}
|
2013-12-08 17:36:39 +00:00
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
AtmosphereModelBruneton::~AtmosphereModelBruneton() {
|
2015-09-24 22:12:31 +00:00
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
AtmosphereResult AtmosphereModelBruneton::getSkyColor(Vector3 eye, const Vector3 &direction,
|
2016-01-10 16:14:54 +00:00
|
|
|
const Vector3 &sun_position, const Color &base) const {
|
2016-01-15 17:01:01 +00:00
|
|
|
auto definition = parent->getScenery()->getAtmosphere();
|
2016-01-14 23:39:33 +00:00
|
|
|
Vector3 x = {0.0, Rg + eye.y * Scenery::UNIT_TO_KM, 0.0};
|
2013-12-11 10:32:10 +00:00
|
|
|
Vector3 v = direction.normalize();
|
|
|
|
Vector3 s = sun_position.sub(x).normalize();
|
2012-12-02 11:08:56 +00:00
|
|
|
|
2013-12-11 10:32:10 +00:00
|
|
|
double r = x.getNorm();
|
|
|
|
double mu = x.dotProduct(v) / r;
|
2012-12-02 11:08:56 +00:00
|
|
|
double t = -r * mu - sqrt(r * r * (mu * mu - 1.0) + Rg * Rg);
|
|
|
|
|
2013-03-14 17:29:12 +00:00
|
|
|
AtmosphereResult result;
|
2012-12-02 11:08:56 +00:00
|
|
|
Vector3 attenuation;
|
2016-01-15 17:01:01 +00:00
|
|
|
Color sunColor = _sunColor(v, s, r, mu, definition->childSun()->getAngularRadius()); /* L0 */
|
2013-01-14 19:07:56 +00:00
|
|
|
|
2013-03-18 21:07:57 +00:00
|
|
|
/*result.base.r = base.r + sunColor.r;
|
2013-03-14 17:29:12 +00:00
|
|
|
result.base.g = base.g + sunColor.g;
|
2013-03-18 21:07:57 +00:00
|
|
|
result.base.b = base.b + sunColor.b;*/
|
2013-12-26 14:36:15 +00:00
|
|
|
result.base = base.add(sunColor);
|
2013-03-14 17:29:12 +00:00
|
|
|
result.inscattering = _getInscatterColor(&x, &t, v, s, &r, &mu, &attenuation); /* S[L]-T(x,xs)S[l]|xs */
|
|
|
|
/* TODO Use atmosphere attenuation */
|
2016-01-14 23:39:33 +00:00
|
|
|
result.distance = Scenery::FAR_LIMIT_SCALED;
|
2013-01-16 15:26:33 +00:00
|
|
|
|
2013-12-08 19:54:34 +00:00
|
|
|
result.updateFinal();
|
2013-03-14 17:29:12 +00:00
|
|
|
|
|
|
|
return result;
|
2012-12-24 11:24:27 +00:00
|
|
|
}
|
|
|
|
|
2016-01-10 16:14:54 +00:00
|
|
|
AtmosphereResult AtmosphereModelBruneton::applyAerialPerspective(Vector3 location, const Color &base) const {
|
2016-01-15 17:01:01 +00:00
|
|
|
auto definition = parent->getScenery()->getAtmosphere();
|
2016-01-14 23:07:02 +00:00
|
|
|
Vector3 eye = parent->getCameraLocation();
|
2015-12-17 18:51:40 +00:00
|
|
|
eye.y = max(eye.y, 0.0);
|
|
|
|
location.y = max(location.y, 0.0);
|
2016-01-15 17:01:01 +00:00
|
|
|
Vector3 sun_position = definition->childSun()->getLocation(false).scale(Scenery::UNIT_TO_KM);
|
2012-12-24 11:24:27 +00:00
|
|
|
|
2016-01-14 23:39:33 +00:00
|
|
|
Vector3 direction = location.sub(eye).scale(Scenery::UNIT_TO_KM);
|
2015-10-08 22:43:32 +00:00
|
|
|
double t = direction.getNorm();
|
2015-11-09 21:30:46 +00:00
|
|
|
if (t < 0.000001) {
|
2016-01-14 23:39:33 +00:00
|
|
|
direction = parent->getCameraDirection().scale(0.001 * Scenery::UNIT_TO_KM);
|
2015-10-08 22:43:32 +00:00
|
|
|
t = direction.getNorm();
|
|
|
|
}
|
2013-01-30 22:06:29 +00:00
|
|
|
|
2016-01-14 23:39:33 +00:00
|
|
|
Vector3 x = {0.0, Rg + WORKAROUND_OFFSET + eye.y * Scenery::UNIT_TO_KM, 0.0};
|
2013-12-11 10:32:10 +00:00
|
|
|
Vector3 v = direction.normalize();
|
|
|
|
Vector3 s = sun_position.sub(x).normalize();
|
2012-12-24 11:24:27 +00:00
|
|
|
|
2015-12-21 16:48:30 +00:00
|
|
|
if (v.y > s.y) {
|
|
|
|
v.y = s.y;
|
|
|
|
v.normalize();
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
if (v.y == 0.0) {
|
2013-01-16 15:26:33 +00:00
|
|
|
v.y = -0.000001;
|
|
|
|
}
|
|
|
|
|
2013-12-11 10:32:10 +00:00
|
|
|
double r = x.getNorm();
|
|
|
|
double mu = x.dotProduct(v) / r;
|
2012-12-24 11:24:27 +00:00
|
|
|
|
2013-03-14 17:29:12 +00:00
|
|
|
AtmosphereResult result;
|
2012-12-24 11:24:27 +00:00
|
|
|
Vector3 attenuation;
|
2013-01-14 19:07:56 +00:00
|
|
|
|
2013-03-14 17:29:12 +00:00
|
|
|
result.base = base;
|
|
|
|
result.inscattering = _getInscatterColor(&x, &t, v, s, &r, &mu, &attenuation); /* S[L]-T(x,xs)S[l]|xs */
|
|
|
|
result.attenuation.r = attenuation.x;
|
|
|
|
result.attenuation.g = attenuation.y;
|
|
|
|
result.attenuation.b = attenuation.z;
|
2016-01-14 23:39:33 +00:00
|
|
|
result.distance = t * Scenery::KM_TO_UNIT;
|
2013-01-16 15:26:33 +00:00
|
|
|
|
2013-12-08 19:54:34 +00:00
|
|
|
result.updateFinal();
|
2013-03-14 17:29:12 +00:00
|
|
|
|
|
|
|
return result;
|
2012-12-02 11:08:56 +00:00
|
|
|
}
|
2013-01-20 15:07:45 +00:00
|
|
|
|
2015-12-10 23:36:50 +00:00
|
|
|
bool AtmosphereModelBruneton::getLightsAt(vector<LightComponent> &result, const Vector3 &location) const {
|
2016-01-15 17:01:01 +00:00
|
|
|
auto definition = parent->getScenery()->getAtmosphere();
|
2013-12-08 19:54:34 +00:00
|
|
|
LightComponent sun, irradiance;
|
2013-01-20 15:07:45 +00:00
|
|
|
double muS;
|
|
|
|
|
2016-01-14 23:39:33 +00:00
|
|
|
double altitude = max(location.y * Scenery::UNIT_TO_KM, 0.0);
|
2013-01-30 22:06:29 +00:00
|
|
|
|
2015-08-20 23:19:55 +00:00
|
|
|
double r0 = Rg + WORKAROUND_OFFSET + altitude;
|
2013-01-20 15:07:45 +00:00
|
|
|
Vector3 up = {0.0, 1.0, 0.0};
|
2016-01-15 17:01:01 +00:00
|
|
|
Vector3 sun_position = definition->childSun()->getLocation(false).scale(Scenery::UNIT_TO_KM);
|
2013-01-20 15:07:45 +00:00
|
|
|
Vector3 x = {0.0, r0, 0.0};
|
2013-12-11 10:32:10 +00:00
|
|
|
Vector3 s = sun_position.sub(x).normalize();
|
2013-01-20 15:07:45 +00:00
|
|
|
|
2013-12-11 10:32:10 +00:00
|
|
|
muS = up.dotProduct(s);
|
2015-11-09 21:30:46 +00:00
|
|
|
if (altitude > RL) {
|
2016-01-15 17:01:01 +00:00
|
|
|
sun.color = definition->sun_color;
|
2015-11-09 21:30:46 +00:00
|
|
|
} else {
|
2013-12-26 15:35:29 +00:00
|
|
|
sun.color = _transmittanceWithShadow(r0, muS);
|
|
|
|
}
|
2013-12-11 10:32:10 +00:00
|
|
|
sun.direction = s.scale(-1.0);
|
2013-03-03 18:10:14 +00:00
|
|
|
sun.reflection = ISun;
|
2013-01-20 15:07:45 +00:00
|
|
|
sun.altered = 1;
|
|
|
|
|
2015-09-24 22:12:31 +00:00
|
|
|
result.push_back(sun);
|
2013-01-20 15:07:45 +00:00
|
|
|
|
2016-01-17 17:27:18 +00:00
|
|
|
irradiance.color = _irradiance(_irradianceTexture, r0, muS).scaled(1.0 - min(location.y / Scenery::ATMOSPHERE_WIDTH_SCALED, 1.0));
|
2013-01-22 20:50:37 +00:00
|
|
|
irradiance.direction = VECTOR_DOWN;
|
2013-01-20 15:07:45 +00:00
|
|
|
irradiance.reflection = 0.0;
|
|
|
|
irradiance.altered = 0;
|
|
|
|
|
2015-09-24 22:12:31 +00:00
|
|
|
result.push_back(irradiance);
|
|
|
|
|
|
|
|
return true;
|
2013-01-20 15:07:45 +00:00
|
|
|
}
|
2013-12-21 22:48:54 +00:00
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
Texture2D *AtmosphereModelBruneton::getTextureTransmittance() const {
|
2013-12-21 22:48:54 +00:00
|
|
|
return _transmittanceTexture;
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
Texture2D *AtmosphereModelBruneton::getTextureIrradiance() const {
|
2013-12-21 22:48:54 +00:00
|
|
|
return _irradianceTexture;
|
|
|
|
}
|
|
|
|
|
2015-11-09 21:30:46 +00:00
|
|
|
Texture4D *AtmosphereModelBruneton::getTextureInscatter() const {
|
2013-12-21 22:48:54 +00:00
|
|
|
return _inscatterTexture;
|
|
|
|
}
|