Added bruneton original source code for atmospheric scattering

This commit is contained in:
Michaël Lemaire 2013-12-25 16:42:56 +01:00
parent 7d2d022a98
commit 19d4c3444f
18 changed files with 2948 additions and 0 deletions

View file

@ -0,0 +1,723 @@
/**
* 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
*/
#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <cmath>
#include <cstring>
#include <GL/glew.h>
#include <GL/glut.h>
#include "tiffio.h"
#include "vec3.h"
#include "mat4.h"
#include "Main.h"
using namespace std;
// ----------------------------------------------------------------------------
// TOOLS
// ----------------------------------------------------------------------------
void loadTIFF(const char *name, unsigned char *tex)
{
tstrip_t strip = 0;
tsize_t off = 0;
tsize_t n = 0;
TIFF* tf = TIFFOpen(name, "r");
while ((n = TIFFReadEncodedStrip(tf, strip, tex + off, (tsize_t) -1)) > 0) {
strip += 1;
off += n;
};
TIFFClose(tf);
}
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();
}
string line;
while (getline(file, line)) {
*result += line;
*result += '\n';
}
file.close();
return result;
}
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;
}
}
unsigned int loadProgram(const vector<string> &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 + 1];
cout << "loading program " << files[n - 1] << "..." << endl;
bool geo = false;
for (int i = 0; i < n; ++i) {
string* s = loadFile(files[i]);
strs[i] = s;
lines[i + 1] = s->c_str();
if (strstr(lines[i + 1], "_GEOMETRY_") != NULL) {
geo = true;
}
}
lines[0] = "#define _VERTEX_\n";
glShaderSource(vertexShaderId, n + 1, lines, NULL);
glCompileShader(vertexShaderId);
printShaderLog(vertexShaderId);
if (geo) {
unsigned geometryShaderId = glCreateShader(GL_GEOMETRY_SHADER_EXT);
glAttachShader(programId, geometryShaderId);
lines[0] = "#define _GEOMETRY_\n";
glShaderSource(geometryShaderId, n + 1, lines, NULL);
glCompileShader(geometryShaderId);
printShaderLog(geometryShaderId);
glProgramParameteriEXT(programId, GL_GEOMETRY_INPUT_TYPE_EXT, GL_TRIANGLES);
glProgramParameteriEXT(programId, GL_GEOMETRY_OUTPUT_TYPE_EXT, GL_TRIANGLE_STRIP);
glProgramParameteriEXT(programId, GL_GEOMETRY_VERTICES_OUT_EXT, 3);
}
lines[0] = "#define _FRAGMENT_\n";
glShaderSource(fragmentShaderId, n + 1, lines, NULL);
glCompileShader(fragmentShaderId);
printShaderLog(fragmentShaderId);
glLinkProgram(programId);
for (int i = 0; i < n; ++i) {
delete strs[i];
}
delete[] strs;
delete[] lines;
return programId;
}
void drawQuad()
{
glBegin(GL_TRIANGLE_STRIP);
glVertex2f(-1.0, -1.0);
glVertex2f(+1.0, -1.0);
glVertex2f(-1.0, +1.0);
glVertex2f(+1.0, +1.0);
glEnd();
}
// ----------------------------------------------------------------------------
// PRECOMPUTATIONS
// ----------------------------------------------------------------------------
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;
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);
}
void loadData()
{
//return;
cout << "loading Earth texture..." << endl;
glActiveTexture(GL_TEXTURE0 + reflectanceUnit);
glGenTextures(1, &reflectanceTexture);
glBindTexture(GL_TEXTURE_2D, reflectanceTexture);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
unsigned char *tex = new unsigned char[2500 * 1250 * 4];
loadTIFF("earth.tiff", tex);
glBindBuffer(GL_PIXEL_UNPACK_BUFFER_ARB, 0);
glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, 2500, 1250, 0, GL_RGBA, GL_UNSIGNED_BYTE, tex);
glGenerateMipmapEXT(GL_TEXTURE_2D);
delete[] tex;
}
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_ARB, 0);
glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB16F_ARB, 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_ARB, 0);
glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB16F_ARB, 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_ARB, 0);
glTexImage3D(GL_TEXTURE_3D, 0, GL_RGBA16F_ARB, 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_ARB, 0);
glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB16F_ARB, 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_ARB, 0);
glTexImage3D(GL_TEXTURE_3D, 0, GL_RGB16F_ARB, 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_ARB, 0);
glTexImage3D(GL_TEXTURE_3D, 0, GL_RGB16F_ARB, 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_ARB, 0);
glTexImage3D(GL_TEXTURE_3D, 0, GL_RGB16F_ARB, RES_MU_S * RES_NU, RES_MU, RES_R, 0, GL_RGB, GL_FLOAT, NULL);
vector<string> files;
files.push_back("Main.h");
files.push_back("common.glsl");
files.push_back("transmittance.glsl");
transmittanceProg = loadProgram(files);
files.clear();
files.push_back("Main.h");
files.push_back("common.glsl");
files.push_back("irradiance1.glsl");
irradiance1Prog = loadProgram(files);
files.clear();
files.push_back("Main.h");
files.push_back("common.glsl");
files.push_back("inscatter1.glsl");
inscatter1Prog = loadProgram(files);
files.clear();
files.push_back("Main.h");
files.push_back("common.glsl");
files.push_back("copyIrradiance.glsl");
copyIrradianceProg = loadProgram(files);
files.clear();
files.push_back("Main.h");
files.push_back("common.glsl");
files.push_back("copyInscatter1.glsl");
copyInscatter1Prog = loadProgram(files);
files.clear();
files.push_back("Main.h");
files.push_back("common.glsl");
files.push_back("inscatterS.glsl");
jProg = loadProgram(files);
files.clear();
files.push_back("Main.h");
files.push_back("common.glsl");
files.push_back("irradianceN.glsl");
irradianceNProg = loadProgram(files);
files.clear();
files.push_back("Main.h");
files.push_back("common.glsl");
files.push_back("inscatterN.glsl");
inscatterNProg = loadProgram(files);
files.clear();
files.push_back("Main.h");
files.push_back("common.glsl");
files.push_back("copyInscatterN.glsl");
copyInscatterNProg = loadProgram(files);
files.clear();
files.push_back("Main.h");
files.push_back("common.glsl");
files.push_back("earth.glsl");
drawProg = loadProgram(files);
glUseProgram(drawProg);
glUniform1i(glGetUniformLocation(drawProg, "reflectanceSampler"), reflectanceUnit);
glUniform1i(glGetUniformLocation(drawProg, "transmittanceSampler"), transmittanceUnit);
glUniform1i(glGetUniformLocation(drawProg, "irradianceSampler"), irradianceUnit);
glUniform1i(glGetUniformLocation(drawProg, "inscatterSampler"), inscatterUnit);
glUniform1i(glGetUniformLocation(drawProg, "debug2DSampler"), deltaEUnit);
glUniform1i(glGetUniformLocation(drawProg, "debug3DSampler"), deltaJUnit);
cout << "precomputations..." << endl;
glGenFramebuffersEXT(1, &fbo);
glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, fbo);
glReadBuffer(GL_COLOR_ATTACHMENT0_EXT);
glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
// computes transmittance texture T (line 1 in algorithm 4.1)
glFramebufferTextureEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, transmittanceTexture, 0);
glViewport(0, 0, TRANSMITTANCE_W, TRANSMITTANCE_H);
glUseProgram(transmittanceProg);
drawQuad();
// computes irradiance texture deltaE (line 2 in algorithm 4.1)
glFramebufferTextureEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, deltaETexture, 0);
glViewport(0, 0, SKY_W, SKY_H);
glUseProgram(irradiance1Prog);
glUniform1i(glGetUniformLocation(irradiance1Prog, "transmittanceSampler"), transmittanceUnit);
drawQuad();
// computes single scattering texture deltaS (line 3 in algorithm 4.1)
// Rayleigh and Mie separated in deltaSR + deltaSM
glFramebufferTextureEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, deltaSRTexture, 0);
glFramebufferTextureEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT1_EXT, deltaSMTexture, 0);
unsigned int bufs[2] = { GL_COLOR_ATTACHMENT0_EXT, GL_COLOR_ATTACHMENT1_EXT };
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();
}
glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT1_EXT, GL_TEXTURE_2D, 0, 0);
glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
// copies deltaE into irradiance texture E (line 4 in algorithm 4.1)
glFramebufferTextureEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, irradianceTexture, 0);
glViewport(0, 0, SKY_W, SKY_H);
glUseProgram(copyIrradianceProg);
glUniform1f(glGetUniformLocation(copyIrradianceProg, "k"), 0.0);
glUniform1i(glGetUniformLocation(copyIrradianceProg, "deltaESampler"), deltaEUnit);
drawQuad();
// copies deltaS into inscatter texture S (line 5 in algorithm 4.1)
glFramebufferTextureEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, 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();
}
// loop for each scattering order (line 6 in algorithm 4.1)
for (int order = 2; order <= 4; ++order) {
// computes deltaJ (line 7 in algorithm 4.1)
glFramebufferTextureEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, deltaJTexture, 0);
glViewport(0, 0, RES_MU_S * RES_NU, RES_MU);
glUseProgram(jProg);
glUniform1f(glGetUniformLocation(jProg, "first"), order == 2 ? 1.0 : 0.0);
glUniform1i(glGetUniformLocation(jProg, "transmittanceSampler"), transmittanceUnit);
glUniform1i(glGetUniformLocation(jProg, "deltaESampler"), deltaEUnit);
glUniform1i(glGetUniformLocation(jProg, "deltaSRSampler"), deltaSRUnit);
glUniform1i(glGetUniformLocation(jProg, "deltaSMSampler"), deltaSMUnit);
for (int layer = 0; layer < RES_R; ++layer) {
setLayer(jProg, layer);
drawQuad();
}
// computes deltaE (line 8 in algorithm 4.1)
glFramebufferTextureEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, deltaETexture, 0);
glViewport(0, 0, SKY_W, SKY_H);
glUseProgram(irradianceNProg);
glUniform1f(glGetUniformLocation(irradianceNProg, "first"), order == 2 ? 1.0 : 0.0);
glUniform1i(glGetUniformLocation(irradianceNProg, "transmittanceSampler"), transmittanceUnit);
glUniform1i(glGetUniformLocation(irradianceNProg, "deltaSRSampler"), deltaSRUnit);
glUniform1i(glGetUniformLocation(irradianceNProg, "deltaSMSampler"), deltaSMUnit);
drawQuad();
// computes deltaS (line 9 in algorithm 4.1)
glFramebufferTextureEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, deltaSRTexture, 0);
glViewport(0, 0, RES_MU_S * RES_NU, RES_MU);
glUseProgram(inscatterNProg);
glUniform1f(glGetUniformLocation(inscatterNProg, "first"), order == 2 ? 1.0 : 0.0);
glUniform1i(glGetUniformLocation(inscatterNProg, "transmittanceSampler"), transmittanceUnit);
glUniform1i(glGetUniformLocation(inscatterNProg, "deltaJSampler"), deltaJUnit);
for (int layer = 0; layer < RES_R; ++layer) {
setLayer(inscatterNProg, layer);
drawQuad();
}
glEnable(GL_BLEND);
glBlendEquationSeparate(GL_FUNC_ADD, GL_FUNC_ADD);
glBlendFuncSeparate(GL_ONE, GL_ONE, GL_ONE, GL_ONE);
// adds deltaE into irradiance texture E (line 10 in algorithm 4.1)
glFramebufferTextureEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, irradianceTexture, 0);
glViewport(0, 0, SKY_W, SKY_H);
glUseProgram(copyIrradianceProg);
glUniform1f(glGetUniformLocation(copyIrradianceProg, "k"), 1.0);
glUniform1i(glGetUniformLocation(copyIrradianceProg, "deltaESampler"), deltaEUnit);
drawQuad();
// adds deltaS into inscatter texture S (line 11 in algorithm 4.1)
glFramebufferTextureEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, inscatterTexture, 0);
glViewport(0, 0, RES_MU_S * RES_NU, RES_MU);
glUseProgram(copyInscatterNProg);
glUniform1i(glGetUniformLocation(copyInscatterNProg, "deltaSSampler"), deltaSRUnit);
for (int layer = 0; layer < RES_R; ++layer) {
setLayer(copyInscatterNProg, layer);
drawQuad();
}
glDisable(GL_BLEND);
// DEBUG Stop
/*glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0);
glFinish();
cout << "ready." << endl;
glUseProgram(drawProg);
return;*/
}
glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0);
glFinish();
cout << "ready." << endl;
glUseProgram(drawProg);
}
void recompute()
{
glDeleteTextures(1, &transmittanceTexture);
glDeleteTextures(1, &irradianceTexture);
glDeleteTextures(1, &inscatterTexture);
glDeleteTextures(1, &deltaETexture);
glDeleteTextures(1, &deltaSRTexture);
glDeleteTextures(1, &deltaSMTexture);
glDeleteTextures(1, &deltaJTexture);
glDeleteProgram(transmittanceProg);
glDeleteProgram(irradiance1Prog);
glDeleteProgram(inscatter1Prog);
glDeleteProgram(copyIrradianceProg);
glDeleteProgram(copyInscatter1Prog);
glDeleteProgram(jProg);
glDeleteProgram(irradianceNProg);
glDeleteProgram(inscatterNProg);
glDeleteProgram(copyInscatterNProg);
glDeleteProgram(drawProg);
glDeleteFramebuffersEXT(1, &fbo);
precompute();
}
// ----------------------------------------------------------------------------
// RENDERING
// ----------------------------------------------------------------------------
int width, height;
int oldx, oldy;
int move;
vec3f s(0.0, -1.0, 0.0);
double lon = 0.0;
double lat = 0.0;
double theta = 0.0;
double phi = 0.0;
double d = Rg;
vec3d position;
mat4d view;
double exposure = 0.4;
void updateView()
{
double co = cos(lon);
double so = sin(lon);
double ca = cos(lat);
double sa = sin(lat);
vec3d po = vec3d(co*ca, so*ca, sa) * Rg;
vec3d px = vec3d(-so, co, 0);
vec3d py = vec3d(-co*sa, -so*sa, ca);
vec3d pz = vec3d(co*ca, so*ca, sa);
double ct = cos(theta);
double st = sin(theta);
double cp = cos(phi);
double sp = sin(phi);
vec3d cx = px * cp + py * sp;
vec3d cy = -px * sp*ct + py * cp*ct + pz * st;
vec3d cz = px * sp*st - py * cp*st + pz * ct;
position = po + cz * d;
if (position.length() < Rg + 0.01) {
position.normalize(Rg + 0.01);
}
view = mat4d(cx.x, cx.y, cx.z, 0,
cy.x, cy.y, cy.z, 0,
cz.x, cz.y, cz.z, 0,
0, 0, 0, 1);
view = view * mat4d::translate(-position);
}
void redisplayFunc()
{
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
float h = position.length() - Rg;
float vfov = 2.0 * atan(float(height) / float(width) * tan(80.0 / 180 * M_PI / 2.0)) / M_PI * 180;
mat4f proj = mat4f::perspectiveProjection(vfov, float(width) / float(height), 0.1 * h, 1e5 * h);
mat4f iproj = proj.inverse();
mat4d iview = view.inverse();
vec3d c = iview * vec3d(0.0, 0.0, 0.0);
mat4f iviewf = mat4f(iview[0][0], iview[0][1], iview[0][2], iview[0][3],
iview[1][0], iview[1][1], iview[1][2], iview[1][3],
iview[2][0], iview[2][1], iview[2][2], iview[2][3],
iview[3][0], iview[3][1], iview[3][2], iview[3][3]);
glUniform3f(glGetUniformLocation(drawProg, "c"), c.x, c.y, c.z);
glUniform3f(glGetUniformLocation(drawProg, "s"), s.x, s.y, s.z);
glUniformMatrix4fv(glGetUniformLocation(drawProg, "projInverse"), 1, true, iproj.coefficients());
glUniformMatrix4fv(glGetUniformLocation(drawProg, "viewInverse"), 1, true, iviewf.coefficients());
glUniform1f(glGetUniformLocation(drawProg, "exposure"), exposure);
drawQuad();
glutSwapBuffers();
}
// ----------------------------------------------------------------------------
// USER INTERFACE
// ----------------------------------------------------------------------------
void reshapeFunc(int x, int y)
{
width = x;
height = y;
glViewport(0, 0, x, y);
glutPostRedisplay();
}
void mouseClickFunc(int button, int state, int x, int y)
{
oldx = x;
oldy = y;
int modifiers = glutGetModifiers();
bool ctrl = (modifiers & GLUT_ACTIVE_CTRL) != 0;
bool shift = (modifiers & GLUT_ACTIVE_SHIFT) != 0;
if (ctrl) {
move = 0;
} else if (shift) {
move = 1;
} else {
move = 2;
}
}
void mouseMotionFunc(int x, int y)
{
if (move == 0) {
phi += (oldx - x) / 500.0;
theta += (oldy - y) / 500.0;
theta = max(0.0, min(M_PI, theta));
updateView();
oldx = x;
oldy = y;
} else if (move == 1) {
double factor = position.length() - Rg;
factor = factor / Rg;
lon += (oldx - x) / 400.0 * factor;
lat -= (oldy - y) / 400.0 * factor;
lat = max(-M_PI / 2.0, min(M_PI / 2.0, lat));
updateView();
oldx = x;
oldy = y;
} else if (move == 2) {
float vangle = asin(s.z);
float hangle = atan2(s.y, s.x);
vangle += (oldy - y) / 180.0 * M_PI / 4;
hangle += (oldx - x) / 180.0 * M_PI / 4;
s.x = cos(vangle) * cos(hangle);
s.y = cos(vangle) * sin(hangle);
s.z = sin(vangle);
oldx = x;
oldy = y;
}
}
void specialKeyFunc(int c, int x, int y)
{
switch (c) {
case GLUT_KEY_PAGE_UP:
d = d * 1.05;
updateView();
break;
case GLUT_KEY_PAGE_DOWN:
d = d / 1.05;
updateView();
break;
case GLUT_KEY_F5:
recompute();
glViewport(0, 0, width, height);
break;
}
}
void keyboardFunc(unsigned char c, int x, int y)
{
if (c == 27) {
::exit(0);
} else if (c == '+') {
exposure *= 1.1;
} else if (c == '-') {
exposure /= 1.1;
}
}
void idleFunc()
{
glutPostRedisplay();
}
int main(int argc, char* argv[])
{
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
glutInitWindowSize(1024, 768);
glutCreateWindow("Precomputed Atmospheric Scattering");
glutCreateMenu(NULL);
glutDisplayFunc(redisplayFunc);
glutReshapeFunc(reshapeFunc);
glutMouseFunc(mouseClickFunc);
glutMotionFunc(mouseMotionFunc);
glutSpecialFunc(specialKeyFunc);
glutKeyboardFunc(keyboardFunc);
glutIdleFunc(idleFunc);
glewInit();
loadData();
precompute();
updateView();
glutMainLoop();
}

View file

@ -0,0 +1,48 @@
/**
* 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
*/
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;

View file

@ -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

View file

@ -0,0 +1,25 @@
TEMPLATE = app
CONFIG += console
CONFIG -= app_bundle
CONFIG -= qt
TARGET = paysages-experiment-bruneton
include(../../common.pri)
QMAKE_CXXFLAGS -= -std=c++11
unix {
CONFIG += link_pkgconfig
PKGCONFIG += gl glu glew libtiff-4
LIBS += -lglut
}
SOURCES += \
Main.cpp
HEADERS += \
Main.h \
mat4.h \
vec3.h

View file

@ -0,0 +1,280 @@
/**
* 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 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);
}
#ifdef _FRAGMENT_
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
}
#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);
}
#ifdef _FRAGMENT_
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);
}
#endif
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;
}
#ifdef _FRAGMENT_
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
}
#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);
}

View file

@ -0,0 +1,79 @@
/**
* 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_
void main() {
gl_Position = gl_Vertex;
}
#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

View file

@ -0,0 +1,78 @@
/**
* 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_
void main() {
gl_Position = gl_Vertex;
}
#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

View file

@ -0,0 +1,55 @@
/**
* 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_
void main() {
gl_Position = gl_Vertex;
}
#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

View file

@ -0,0 +1,232 @@
/**
* 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
*/
#define FIX
const float ISun = 100.0;
uniform vec3 c;
uniform vec3 s;
uniform mat4 projInverse;
uniform mat4 viewInverse;
uniform float exposure;
uniform sampler2D reflectanceSampler;//ground reflectance texture
uniform sampler2D irradianceSampler;//precomputed skylight irradiance (E table)
uniform sampler3D inscatterSampler;//precomputed inscattered light (S table)
uniform sampler2D debug2DSampler;
uniform sampler3D debug3DSampler;
varying vec2 coords;
varying vec3 ray;
#ifdef _VERTEX_
void main() {
coords = gl_Vertex.xy * 0.5 + 0.5;
ray = (viewInverse * vec4((projInverse * gl_Vertex).xyz, 0.0)).xyz;
gl_Position = gl_Vertex;
}
#else
//inscattered light along ray x+tv, when sun in direction s (=S[L]-T(x,x0)S[L]|x0)
vec3 inscatter(inout vec3 x, inout float t, vec3 v, vec3 s, out float r, out float mu, out vec3 attenuation) {
vec3 result;
r = length(x);
mu = dot(x, v) / r;
float d = -r * mu - sqrt(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;
t -= d;
mu = (r * mu + d) / Rt;
r = Rt;
}
if (r <= Rt) { // if ray intersects atmosphere
float nu = dot(v, s);
float muS = dot(x, s) / r;
float phaseR = phaseFunctionR(nu);
float phaseM = phaseFunctionM(nu);
vec4 inscatter = max(texture4D(inscatterSampler, r, mu, muS, nu), 0.0);
if (t > 0.0) {
vec3 x0 = x + t * v;
float r0 = length(x0);
float rMu0 = dot(x0, v);
float mu0 = rMu0 / r0;
float muS0 = dot(x0, s) / r0;
#ifdef FIX
// avoids imprecision problems in transmittance computations based on textures
attenuation = analyticTransmittance(r, mu, t);
#else
attenuation = transmittance(r, mu, v, x0);
#endif
if (r0 > Rg + 0.01) {
// computes S[L]-T(x,x0)S[L]|x0
inscatter = max(inscatter - attenuation.rgbr * texture4D(inscatterSampler, r0, mu0, muS0, nu), 0.0);
#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));
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);
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);
mu0 = (r * mu + t) / r0;
inScatter0 = texture4D(inscatterSampler, r, mu, muS, nu);
inScatter1 = texture4D(inscatterSampler, r0, mu0, muS0, nu);
vec4 inScatterB = max(inScatter0 - attenuation.rgbr * inScatter1, 0.0);
inscatter = mix(inScatterA, inScatterB, a);
}
#endif
}
}
#ifdef FIX
// avoids imprecision problems in Mie scattering when sun is below horizon
inscatter.w *= smoothstep(0.00, 0.02, muS);
#endif
result = max(inscatter.rgb * phaseR + getMie(inscatter) * phaseM, 0.0);
} else { // x in space and ray looking in space
result = vec3(0.0);
}
return result * ISun;
}
//ground radiance at end of ray x+tv, when sun in direction s
//attenuated bewteen ground and viewer (=R[L0]+R[L*])
vec3 groundColor(vec3 x, float t, vec3 v, vec3 s, float r, float mu, vec3 attenuation)
{
vec3 result;
if (t > 0.0) { // if ray hits ground surface
// ground reflectance at end of ray, x0
vec3 x0 = x + t * v;
float r0 = length(x0);
vec3 n = x0 / r0;
vec2 coords = vec2(atan(n.y, n.x), acos(n.z)) * vec2(0.5, 1.0) / M_PI + vec2(0.5, 0.0);
vec4 reflectance = texture2D(reflectanceSampler, coords) * vec4(0.2, 0.2, 0.2, 1.0);
if (r0 > Rg + 0.01) {
reflectance = vec4(0.4, 0.4, 0.4, 0.0);
}
// direct sun light (radiance) reaching x0
float muS = dot(n, s);
vec3 sunLight = transmittanceWithShadow(r0, muS);
// precomputed sky light (irradiance) (=E[L*]) at x0
vec3 groundSkyLight = irradiance(irradianceSampler, r0, muS);
// light reflected at x0 (=(R[L0]+R[L*])/T(x,x0))
vec3 groundColor = reflectance.rgb * (max(muS, 0.0) * sunLight + groundSkyLight) * ISun / M_PI;
// 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);
groundColor += reflectance.w * max(waterBrdf, 0.0) * sunLight * ISun;
}
result = attenuation * groundColor; //=R[L0]+R[L*]
} else { // ray looking at the sky
result = vec3(0.0);
}
return result;
}
// direct sun light for ray x+tv, when sun in direction s (=L0)
vec3 sunColor(vec3 x, float t, vec3 v, vec3 s, float r, float mu) {
if (t > 0.0) {
return vec3(0.0);
} else {
vec3 transmittance = r <= Rt ? transmittanceWithShadow(r, mu) : vec3(1.0); // T(x,xo)
float isun = step(cos(M_PI / 180.0), dot(v, s)) * ISun; // Lsun
return transmittance * isun; // Eq (9)
}
}
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);
return L;
}
void main() {
vec3 x = c;
vec3 v = normalize(ray);
float r = length(x);
float mu = dot(x, v) / r;
float t = -r * mu - sqrt(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);
bool cone = d > 0.0 && abs(x.z + d * v.z - Rg) <= 10.0;
if (t > 0.0) {
if (cone && d < t) {
t = d;
}
} else if (cone) {
t = d;
}
vec3 attenuation;
vec3 inscatterColor = inscatter(x, t, v, s, r, mu, attenuation); //S[L]-T(x,xs)S[l]|xs
vec3 groundColor = groundColor(x, t, v, s, r, mu, attenuation); //R[L0]+R[L*]
vec3 sunColor = sunColor(x, t, v, s, r, mu); //L0
gl_FragColor = vec4(HDR(sunColor + groundColor + inscatterColor), 1.0); // Eq (16)
//gl_FragColor = texture3D(inscatterSampler,vec3(coords,(s.x+1.0)/2.0));
//gl_FragColor = vec4(texture2D(irradianceSampler,coords).rgb*5.0, 1.0);
//gl_FragColor = texture2D(transmittanceSampler,coords);
//gl_FragColor = vec4(texture2D(debug2DSampler,coords).rgb * 1.0, 1.0);
//gl_FragColor = texture3D(debug3DSampler,vec3(coords,0.5))*1000.0;
}
#endif

Binary file not shown.

View file

@ -0,0 +1,117 @@
/**
* 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_
void main() {
gl_Position = gl_Vertex;
}
#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

View file

@ -0,0 +1,100 @@
/**
* 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_
void main() {
gl_Position = gl_Vertex;
}
#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

View file

@ -0,0 +1,158 @@
/**
* 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_
void main() {
gl_Position = gl_Vertex;
}
#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

View file

@ -0,0 +1,53 @@
/**
* 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_
void main() {
gl_Position = gl_Vertex;
}
#endif
#ifdef _FRAGMENT_
void main() {
float r, muS;
getIrradianceRMuS(r, muS);
gl_FragColor = vec4(transmittance(r, muS) * max(muS, 0.0), 0.0);
}
#endif

View file

@ -0,0 +1,85 @@
/**
* 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_
void main() {
gl_Position = gl_Vertex;
}
#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

View file

@ -0,0 +1,429 @@
/**
* 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.
*/
#ifndef MAT4_H_
#define MAT4_H_
#include "vec3.h"
#ifndef M_PI
#define M_PI 3.141592657
#endif
/**
* A 4x4 matrix.
*/
template <typename type> class mat4
{
protected:
union {
type m[4][4];
type _m[16];
};
public:
/**
* Creates a new, uninitialized matrix.
*/
mat4();
/**
* Creates a new matrix with the given components. The first index is the
* row index, the second one is the column index.
*/
mat4(type m00, type m01, type m02, type m03,
type m10, type m11, type m12, type m13,
type m20, type m21, type m22, type m23,
type m30, type m31, type m32, type m33);
/**
* Returns the coefficients of this matrix.
*/
const type* coefficients() const;
/**
* Returns the row of this matrix whose index is given.
*/
const type* operator[](int iRow) const;
/**
* Returns true is this matrix is equal to the given matrix.
*/
bool operator==(const mat4& m2) const;
/**
* Returns true is this matrix is different from the given matrix.
*/
bool operator!=(const mat4& m2) const;
/**
* Returns the sum of this matrix and of the given matrix.
*/
mat4 operator+(const mat4& m2) const;
/**
* Returns the difference of this matrix and of the given matrix.
*/
mat4 operator-(const mat4& m2) const;
/**
* Returns the product of this matrix and of the given matrix.
*/
mat4 operator*(const mat4& m2) const;
/**
* Returns the product of this matrix and of the given vector. The given
* vector w coordinate is set to 1, and the 4 vector result is converted
* to a 3 vector by dividing its xyz components by its w component.
*/
vec3<type> operator*(const vec3<type>& v) const;
/**
* Returns the product of this matrix and of the given scalar.
*/
mat4 operator*(type f) const;
/**
* Returns the transpose of this matrix.
*/
mat4 transpose(void) const;
/**
* Returns the adjoint of this matrix.
*/
mat4 adjoint() const;
/**
* Returns the inverse of this matrix.
*/
mat4 inverse() const;
/**
* Returns the determinant of this matrix.
*/
float determinant() const;
/**
* Returns the translation matrix corresponding to the given translation
* vector.
*/
static mat4 translate(const vec3<type>& v);
/**
* Returns the perspective projection matrix corresponding to the given
* projection parameters.
*
* @param fovy vertical field of view in degrees.
* @param aspect aspect ratio of the projection window.
* @param zNear near clipping plane.
* @param zFar far clipping plane.
*/
static mat4 perspectiveProjection(type fovy, type aspect, type zNear, type zFar);
};
/**
* A 4x4 matrix with float components.
*/
typedef mat4<float> mat4f;
/**
* A 4x4 matrix with double components.
*/
typedef mat4<double> mat4d;
template <typename type>
inline mat4<type>::mat4()
{
}
template <typename type>
inline mat4<type>::mat4(type m00, type m01, type m02, type m03,
type m10, type m11, type m12, type m13,
type m20, type m21, type m22, type m23,
type m30, type m31, type m32, type m33)
{
m[0][0] = m00;
m[0][1] = m01;
m[0][2] = m02;
m[0][3] = m03;
m[1][0] = m10;
m[1][1] = m11;
m[1][2] = m12;
m[1][3] = m13;
m[2][0] = m20;
m[2][1] = m21;
m[2][2] = m22;
m[2][3] = m23;
m[3][0] = m30;
m[3][1] = m31;
m[3][2] = m32;
m[3][3] = m33;
}
template <typename type>
inline const type* mat4<type>::coefficients() const
{
return _m;
}
template <typename type>
inline const type* mat4<type>::operator[](int iRow) const
{
//assert(iRow < 4);
return m[iRow];
}
template <typename type>
inline bool mat4<type>::operator==(const mat4<type>& m2) const
{
if (
m[0][0] != m2.m[0][0] || m[0][1] != m2.m[0][1] || m[0][2] != m2.m[0][2] || m[0][3] != m2.m[0][3] ||
m[1][0] != m2.m[1][0] || m[1][1] != m2.m[1][1] || m[1][2] != m2.m[1][2] || m[1][3] != m2.m[1][3] ||
m[2][0] != m2.m[2][0] || m[2][1] != m2.m[2][1] || m[2][2] != m2.m[2][2] || m[2][3] != m2.m[2][3] ||
m[3][0] != m2.m[3][0] || m[3][1] != m2.m[3][1] || m[3][2] != m2.m[3][2] || m[3][3] != m2.m[3][3])
return false;
return true;
}
template <typename type>
inline bool mat4<type>::operator!=(const mat4<type>& m2) const
{
if (
m[0][0] != m2.m[0][0] || m[0][1] != m2.m[0][1] || m[0][2] != m2.m[0][2] || m[0][3] != m2.m[0][3] ||
m[1][0] != m2.m[1][0] || m[1][1] != m2.m[1][1] || m[1][2] != m2.m[1][2] || m[1][3] != m2.m[1][3] ||
m[2][0] != m2.m[2][0] || m[2][1] != m2.m[2][1] || m[2][2] != m2.m[2][2] || m[2][3] != m2.m[2][3] ||
m[3][0] != m2.m[3][0] || m[3][1] != m2.m[3][1] || m[3][2] != m2.m[3][2] || m[3][3] != m2.m[3][3])
return true;
return false;
}
template <typename type>
inline mat4<type> mat4<type>::operator+(const mat4<type>& m2) const
{
mat4<type> r;
r.m[0][0] = m[0][0] + m2.m[0][0];
r.m[0][1] = m[0][1] + m2.m[0][1];
r.m[0][2] = m[0][2] + m2.m[0][2];
r.m[0][3] = m[0][3] + m2.m[0][3];
r.m[1][0] = m[1][0] + m2.m[1][0];
r.m[1][1] = m[1][1] + m2.m[1][1];
r.m[1][2] = m[1][2] + m2.m[1][2];
r.m[1][3] = m[1][3] + m2.m[1][3];
r.m[2][0] = m[2][0] + m2.m[2][0];
r.m[2][1] = m[2][1] + m2.m[2][1];
r.m[2][2] = m[2][2] + m2.m[2][2];
r.m[2][3] = m[2][3] + m2.m[2][3];
r.m[3][0] = m[3][0] + m2.m[3][0];
r.m[3][1] = m[3][1] + m2.m[3][1];
r.m[3][2] = m[3][2] + m2.m[3][2];
r.m[3][3] = m[3][3] + m2.m[3][3];
return r;
}
template <typename type>
inline mat4<type> mat4<type>::operator-(const mat4<type>& m2) const
{
mat4 r;
r.m[0][0] = m[0][0] - m2.m[0][0];
r.m[0][1] = m[0][1] - m2.m[0][1];
r.m[0][2] = m[0][2] - m2.m[0][2];
r.m[0][3] = m[0][3] - m2.m[0][3];
r.m[1][0] = m[1][0] - m2.m[1][0];
r.m[1][1] = m[1][1] - m2.m[1][1];
r.m[1][2] = m[1][2] - m2.m[1][2];
r.m[1][3] = m[1][3] - m2.m[1][3];
r.m[2][0] = m[2][0] - m2.m[2][0];
r.m[2][1] = m[2][1] - m2.m[2][1];
r.m[2][2] = m[2][2] - m2.m[2][2];
r.m[2][3] = m[2][3] - m2.m[2][3];
r.m[3][0] = m[3][0] - m2.m[3][0];
r.m[3][1] = m[3][1] - m2.m[3][1];
r.m[3][2] = m[3][2] - m2.m[3][2];
r.m[3][3] = m[3][3] - m2.m[3][3];
return r;
}
template <typename type>
inline mat4<type> mat4<type>::operator*(const mat4<type>& m2) const
{
mat4 r;
r.m[0][0] = m[0][0] * m2.m[0][0] + m[0][1] * m2.m[1][0] + m[0][2] * m2.m[2][0] + m[0][3] * m2.m[3][0];
r.m[0][1] = m[0][0] * m2.m[0][1] + m[0][1] * m2.m[1][1] + m[0][2] * m2.m[2][1] + m[0][3] * m2.m[3][1];
r.m[0][2] = m[0][0] * m2.m[0][2] + m[0][1] * m2.m[1][2] + m[0][2] * m2.m[2][2] + m[0][3] * m2.m[3][2];
r.m[0][3] = m[0][0] * m2.m[0][3] + m[0][1] * m2.m[1][3] + m[0][2] * m2.m[2][3] + m[0][3] * m2.m[3][3];
r.m[1][0] = m[1][0] * m2.m[0][0] + m[1][1] * m2.m[1][0] + m[1][2] * m2.m[2][0] + m[1][3] * m2.m[3][0];
r.m[1][1] = m[1][0] * m2.m[0][1] + m[1][1] * m2.m[1][1] + m[1][2] * m2.m[2][1] + m[1][3] * m2.m[3][1];
r.m[1][2] = m[1][0] * m2.m[0][2] + m[1][1] * m2.m[1][2] + m[1][2] * m2.m[2][2] + m[1][3] * m2.m[3][2];
r.m[1][3] = m[1][0] * m2.m[0][3] + m[1][1] * m2.m[1][3] + m[1][2] * m2.m[2][3] + m[1][3] * m2.m[3][3];
r.m[2][0] = m[2][0] * m2.m[0][0] + m[2][1] * m2.m[1][0] + m[2][2] * m2.m[2][0] + m[2][3] * m2.m[3][0];
r.m[2][1] = m[2][0] * m2.m[0][1] + m[2][1] * m2.m[1][1] + m[2][2] * m2.m[2][1] + m[2][3] * m2.m[3][1];
r.m[2][2] = m[2][0] * m2.m[0][2] + m[2][1] * m2.m[1][2] + m[2][2] * m2.m[2][2] + m[2][3] * m2.m[3][2];
r.m[2][3] = m[2][0] * m2.m[0][3] + m[2][1] * m2.m[1][3] + m[2][2] * m2.m[2][3] + m[2][3] * m2.m[3][3];
r.m[3][0] = m[3][0] * m2.m[0][0] + m[3][1] * m2.m[1][0] + m[3][2] * m2.m[2][0] + m[3][3] * m2.m[3][0];
r.m[3][1] = m[3][0] * m2.m[0][1] + m[3][1] * m2.m[1][1] + m[3][2] * m2.m[2][1] + m[3][3] * m2.m[3][1];
r.m[3][2] = m[3][0] * m2.m[0][2] + m[3][1] * m2.m[1][2] + m[3][2] * m2.m[2][2] + m[3][3] * m2.m[3][2];
r.m[3][3] = m[3][0] * m2.m[0][3] + m[3][1] * m2.m[1][3] + m[3][2] * m2.m[2][3] + m[3][3] * m2.m[3][3];
return r;
}
template <typename type>
inline vec3<type> mat4<type>::operator*(const vec3<type>& v) const
{
vec3<type> r;
float fInvW = 1.0 / (m[3][0] * v.x + m[3][1] * v.y + m[3][2] * v.z + m[3][3]);
r.x = (m[0][0] * v.x + m[0][1] * v.y + m[0][2] * v.z + m[0][3]) * fInvW;
r.y = (m[1][0] * v.x + m[1][1] * v.y + m[1][2] * v.z + m[1][3]) * fInvW;
r.z = (m[2][0] * v.x + m[2][1] * v.y + m[2][2] * v.z + m[2][3]) * fInvW;
return r;
}
template <typename type>
inline mat4<type> mat4<type>::operator*(type f) const
{
mat4<type> r;
r.m[0][0] = m[0][0] * f;
r.m[0][1] = m[0][1] * f;
r.m[0][2] = m[0][2] * f;
r.m[0][3] = m[0][3] * f;
r.m[1][0] = m[1][0] * f;
r.m[1][1] = m[1][1] * f;
r.m[1][2] = m[1][2] * f;
r.m[1][3] = m[1][3] * f;
r.m[2][0] = m[2][0] * f;
r.m[2][1] = m[2][1] * f;
r.m[2][2] = m[2][2] * f;
r.m[2][3] = m[2][3] * f;
r.m[3][0] = m[3][0] * f;
r.m[3][1] = m[3][1] * f;
r.m[3][2] = m[3][2] * f;
r.m[3][3] = m[3][3] * f;
return r;
}
template <typename type>
inline mat4<type> mat4<type>::transpose(void) const
{
return mat4(m[0][0], m[1][0], m[2][0], m[3][0],
m[0][1], m[1][1], m[2][1], m[3][1],
m[0][2], m[1][2], m[2][2], m[3][2],
m[0][3], m[1][3], m[2][3], m[3][3]);
}
template <typename type>
inline static float
MINOR(const mat4<type>& m, int r0, int r1, int r2, int c0, int c1, int c2)
{
return m[r0][c0] *(m[r1][c1] * m[r2][c2] - m[r2][c1] * m[r1][c2]) -
m[r0][c1] *(m[r1][c0] * m[r2][c2] - m[r2][c0] * m[r1][c2]) +
m[r0][c2] *(m[r1][c0] * m[r2][c1] - m[r2][c0] * m[r1][c1]);
}
template <typename type>
mat4<type> mat4<type>::adjoint() const
{
return mat4(MINOR(*this, 1, 2, 3, 1, 2, 3),
-MINOR(*this, 0, 2, 3, 1, 2, 3),
MINOR(*this, 0, 1, 3, 1, 2, 3),
-MINOR(*this, 0, 1, 2, 1, 2, 3),
-MINOR(*this, 1, 2, 3, 0, 2, 3),
MINOR(*this, 0, 2, 3, 0, 2, 3),
-MINOR(*this, 0, 1, 3, 0, 2, 3),
MINOR(*this, 0, 1, 2, 0, 2, 3),
MINOR(*this, 1, 2, 3, 0, 1, 3),
-MINOR(*this, 0, 2, 3, 0, 1, 3),
MINOR(*this, 0, 1, 3, 0, 1, 3),
-MINOR(*this, 0, 1, 2, 0, 1, 3),
-MINOR(*this, 1, 2, 3, 0, 1, 2),
MINOR(*this, 0, 2, 3, 0, 1, 2),
-MINOR(*this, 0, 1, 3, 0, 1, 2),
MINOR(*this, 0, 1, 2, 0, 1, 2));
}
template <typename type>
mat4<type> mat4<type>::inverse() const
{
return adjoint() * (1.0f / determinant());
}
template <typename type>
float mat4<type>::determinant() const
{
return m[0][0] * MINOR(*this, 1, 2, 3, 1, 2, 3) -
m[0][1] * MINOR(*this, 1, 2, 3, 0, 2, 3) +
m[0][2] * MINOR(*this, 1, 2, 3, 0, 1, 3) -
m[0][3] * MINOR(*this, 1, 2, 3, 0, 1, 2);
}
template <typename type>
inline mat4<type> mat4<type>::translate(const vec3<type> &v)
{
return mat4<type>(1, 0, 0, v.x,
0, 1, 0, v.y,
0, 0, 1, v.z,
0, 0, 0, 1);
}
template <typename type>
inline mat4<type> mat4<type>::perspectiveProjection(type fovy, type aspect, type zNear, type zFar)
{
type f = (type) 1 / tan(fovy * M_PI / 180.0 / 2);
return mat4<type>(f / aspect, 0, 0, 0,
0, f, 0, 0,
0, 0, (zFar + zNear) / (zNear - zFar), (2*zFar*zNear) / (zNear - zFar),
0, 0, -1, 0);
}
#endif /*MAT4_H_*/

View file

@ -0,0 +1,67 @@
/**
* 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_
void main() {
gl_Position = gl_Vertex;
}
#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

View file

@ -0,0 +1,406 @@
/**
* 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.
*/
#ifndef VEC3_H_
#define VEC3_H_
#include <cassert>
/**
* A 3D vector.
*/
template <typename type> class vec3
{
public:
type x, y, z;
/**
* Creates a new, uninitialized vector.
*/
vec3();
/**
* Creates a new vector with the given coordinates.
*/
vec3(type xi, type yi, type zi);
/**
* Creates a new vector with the given coordinates.
*/
vec3(const type v[3]);
/**
* Creates a new vector as a copy of the given vector.
*/
vec3(const vec3& v);
/**
* Returns the coordinate of this vector whose index is given.
*/
type operator[](const int i) const;
/**
* Returns the coordinate of this vector whose index is given.
*/
type& operator[](const int i);
/**
* Assigns the given vector to this vector.
*/
void operator=(const vec3& v);
/**
* Returns true if this vector is equal to the given vector.
*/
bool operator==(const vec3& v) const;
/**
* Returns true if this vector is different from the given vector.
*/
bool operator!=(const vec3& v) const;
/**
* Returns the sum of this vector and of the given vector.
*/
vec3 operator+(const vec3& v) const;
/**
* Returns the difference of this vector and of the given vector.
*/
vec3 operator-(const vec3& v) const;
/**
* Returns the product of this vector and of the given vector. The
* product is done component by component.
*/
vec3 operator*(const vec3& v) const;
/**
* Returns the product of this vector and of the given scalar.
*/
vec3 operator*(const type scalar) const;
/**
* Returns the division of this vector and of the given vector. The
* division is done component by component.
*/
vec3 operator/(const vec3& v) const;
/**
* Returns the division of this vector and of the given scalar.
*/
vec3 operator/(const type scalar) const;
/**
* Returns the opposite of this vector.
*/
vec3 operator-() const;
/**
* Adds the given vector to this vector.
*/
vec3& operator+=(const vec3& v);
/**
* Substracts the given vector from this vector.
*/
vec3& operator-=(const vec3& v);
/**
* Multiplies this vector by the given scalar.
*/
vec3& operator*=(const type& scalar);
/**
* Divides this vector by the given scalar.
*/
vec3& operator/=(const type& scalar);
/**
* Returns the length of this vector.
*/
type length() const;
/**
* Returns the squared length of this vector.
*/
type squaredlength() const;
/**
* Returns the dot product of this vector and of the given vector.
*/
type dotproduct(const vec3& v) const;
/**
* Normalizes this vector and returns its initial length.
*/
type normalize();
/**
* Normalizes this vector to the given length and returns its initial length.
*/
type normalize(type l);
/**
* Returns he cross product of this vector and of the given vector.
*/
vec3 crossProduct(const vec3& v) const;
/**
* The null vector (0,0,0).
*/
static const vec3 ZERO;
/**
* The unit x vector (1,0,0).
*/
static const vec3 UNIT_X;
/**
* The unit y vector (0,1,0).
*/
static const vec3 UNIT_Y;
/**
* The unit z vector (0,0,1).
*/
static const vec3 UNIT_Z;
};
/**
* A 3D vector with float coordinates.
*/
typedef vec3<float> vec3f;
/**
* A 3D vector with double coordinates.
*/
typedef vec3<double> vec3d;
/**
* A 3D vector with int coordinates.
*/
typedef vec3<int> vec3i;
template <typename type>
inline vec3<type>::vec3()
{
}
template <typename type>
inline vec3<type>::vec3(type xi, type yi, type zi) : x(xi), y(yi), z(zi)
{
}
template <typename type>
inline vec3<type>::vec3(const type v[3]) : x(v[0]), y(v[1]), z(v[2])
{
}
template <typename type>
inline vec3<type>::vec3(const vec3& v) : x(v.x), y(v.y), z(v.z)
{
}
template <typename type>
inline type vec3<type>::operator[](const int i) const
{
//assert(i<3);
return *(&x + i);
}
template <typename type>
inline type& vec3<type>::operator[](const int i)
{
//assert(i<3);
return *(&x + i);
}
template <typename type>
inline void vec3<type>::operator=(const vec3<type>& v)
{
x = v.x;
y = v.y;
z = v.z;
}
template <typename type>
inline bool vec3<type>::operator==(const vec3<type>& v) const
{
return (x == v.x && y == v.y && z == v.z);
}
template <typename type>
inline bool vec3<type>::operator!=(const vec3<type>& v) const
{
return (x != v.x || y != v.y || z != v.z);
}
template <typename type>
inline vec3<type> vec3<type>::operator+(const vec3<type>& v) const
{
return vec3(x + v.x, y + v.y, z + v.z);
}
template <typename type>
inline vec3<type> vec3<type>::operator-(const vec3<type>& v) const
{
return vec3(x - v.x, y - v.y, z - v.z);
}
template <typename type>
inline vec3<type> vec3<type>::operator*(const vec3<type>& v) const
{
return vec3(x * v.x, y * v.y, z * v.z);
}
template <typename type>
inline vec3<type> vec3<type>::operator*(const type scalar) const
{
return vec3(x * scalar, y * scalar, z * scalar);
}
template <typename type>
inline vec3<type> vec3<type>::operator/(const vec3<type>& v) const
{
return vec3(x / v.x, y / v.y, z / v.z);
}
template <typename type>
inline vec3<type> vec3<type>::operator/(const type scalar) const
{
assert(scalar != 0);
type inv = 1 / scalar;
return vec3(x * inv, y * inv, z * inv);
}
template <typename type>
inline vec3<type> vec3<type>::operator-() const
{
return vec3(-x, -y, -z);
}
template <typename type>
inline vec3<type>& vec3<type>::operator+=(const vec3<type>& v)
{
x += v.x;
y += v.y;
z += v.z;
return *this;
}
template <typename type>
inline vec3<type>& vec3<type>::operator-=(const vec3<type>& v)
{
x -= v.x;
y -= v.y;
z -= v.z;
return *this;
}
template <typename type>
inline vec3<type>& vec3<type>::operator*=(const type& scalar)
{
x *= scalar;
y *= scalar;
z *= scalar;
return *this;
}
template <typename type>
inline vec3<type>& vec3<type>::operator/=(const type& scalar)
{
assert(scalar != 0);
type inv = 1 / scalar;
x *= inv;
y *= inv;
z *= inv;
return *this;
}
template <typename type>
inline type vec3<type>::length() const
{
return sqrt(x*x + y*y + z*z);
}
template <typename type>
inline type vec3<type>::squaredlength() const
{
return (x*x + y*y + z*z);
}
template <typename type>
inline type vec3<type>::dotproduct(const vec3<type>& v) const
{
return (x*v.x + y*v.y + z*v.z);
}
template <typename type>
inline type vec3<type>::normalize()
{
type length = sqrt(x * x + y * y + z * z);
type invLength = 1.0 / length;
x *= invLength;
y *= invLength;
z *= invLength;
return length;
}
template <typename type>
inline type vec3<type>::normalize(type l)
{
type length = sqrt(x * x + y * y + z * z);
type invLength = l / length;
x *= invLength;
y *= invLength;
z *= invLength;
return length;
}
template <typename type>
inline vec3<type> vec3<type>::crossProduct(const vec3<type>& v) const
{
return vec3(y*v.z - z*v.y, z*v.x - x*v.z, x*v.y - y*v.x);
}
template <typename type>
const vec3<type> vec3<type>::ZERO(0, 0, 0);
template <typename type>
const vec3<type> vec3<type>::UNIT_X(1, 0, 0);
template <typename type>
const vec3<type> vec3<type>::UNIT_Y(0, 1, 0);
template <typename type>
const vec3<type> vec3<type>::UNIT_Z(0, 0, 1);
#endif /*VEC3_H_*/