c2e5622d50
These are causing quite a big difference in existing files, which is not
easy to address in versioning. Since the goal of removing this was to
simplify things for us and that's not the case, just revert this change.
This reverts commit ab21755aaf.
Ref #139923
Pull Request: https://projects.blender.org/blender/blender/pulls/146336
402 lines
18 KiB
C++
402 lines
18 KiB
C++
/* SPDX-FileCopyrightText: 2022 Fernando García Liñán
|
|
* SPDX-FileCopyrightText: 2011-2025 Blender Authors
|
|
*
|
|
* SPDX-License-Identifier: MIT */
|
|
|
|
/** \file
|
|
* \ingroup intern_sky_modal
|
|
*/
|
|
|
|
/*
|
|
* This code is a converted version of the ShaderToy written by Fernando García Liñán.
|
|
*
|
|
* This shader is the final result of my Master's Thesis.
|
|
* The main contributions are:
|
|
*
|
|
* 1. A spectral rendering technique that only requires 4 wavelength samples to
|
|
* get accurate results.
|
|
* 2. A multiple scattering approximation.
|
|
*
|
|
* Both of these approximations rely on an analytical fit, so they only work for
|
|
* Earth's atmosphere. We make up for it by using a very flexible atmosphere
|
|
* model that is able to represent a wide variety of atmospheric conditions.
|
|
*
|
|
* A brief description of this spectral rendering technique can be found in the
|
|
* following article:
|
|
* https://fgarlin.com/posts/2024-12-06-spectral_sky/
|
|
*
|
|
* The path tracer that has been used as a ground truth can be found at:
|
|
* https://github.com/fgarlin/skytracer
|
|
*/
|
|
|
|
#include <algorithm>
|
|
|
|
#include "sky_math.h"
|
|
#include "sky_nishita.h"
|
|
|
|
using std::min;
|
|
|
|
/* Earth's atmosphere parameters. */
|
|
/* Ground reflectance. */
|
|
static const float4 GROUND_ALBEDO = make_float4(0.3f, 0.3f, 0.3f, 0.3f);
|
|
static const float PHASE_ISOTROPIC = M_1_4PI_F;
|
|
static const float RAYLEIGH_PHASE_SCALE = (3.0f / 16.0f) * M_1_PI_F;
|
|
/* Aerosols anisotropy. */
|
|
static const float G = 0.8f;
|
|
static const float SQR_G = G * G;
|
|
/* Earth radius (km). */
|
|
static const float EARTH_RADIUS = 6371.0f;
|
|
/* Atmosphere thickness (km). */
|
|
static const float ATMOSPHERE_THICKNESS = 100.0f;
|
|
static const float ATMOSPHERE_RADIUS = EARTH_RADIUS + ATMOSPHERE_THICKNESS;
|
|
/* Ray marching steps. Higher steps means increased accuracy but worse performance. */
|
|
static const int TRANSMITTANCE_STEPS = 64;
|
|
static const int IN_SCATTERING_STEPS = 64;
|
|
|
|
/* LUTs. */
|
|
static const int TRANSMITTANCE_RES_X = 256;
|
|
static const int TRANSMITTANCE_RES_Y = 64;
|
|
|
|
/* Spectral data sampled at 630, 560, 490, 430 nm for urban area. */
|
|
static const float4 SUN_SPECTRAL_IRRADIANCE = make_float4(1.679f, 1.828f, 1.986f, 1.307f);
|
|
static const float4 MOLECULAR_SCATTERING_COEFFICIENT_BASE = make_float4(
|
|
6.605e-3f, 1.067e-2f, 1.842e-2f, 3.156e-2f);
|
|
static const float4 OZONE_ABSORPTION_CROSS_SECTION = make_float4(
|
|
3.472e-25f, 3.914e-25f, 1.349e-25f, 11.03e-27f);
|
|
/* Average ozone dobson of monthly mean values. */
|
|
static const float OZONE_MEAN_DOBSON = 334.5f;
|
|
static const float4 AEROSOL_ABSORPTION_CROSS_SECTION = make_float4(
|
|
2.8722e-24f, 4.6168e-24f, 7.9706e-24f, 1.3578e-23f);
|
|
static const float4 AEROSOL_SCATTERING_CROSS_SECTION = make_float4(
|
|
1.5908e-22f, 1.7711e-22f, 2.0942e-22f, 2.4033e-22f);
|
|
static const float AEROSOL_BASE_DENSITY = 1.3681e20f;
|
|
static const float AEROSOL_BACKGROUND_DENSITY = 2e6f;
|
|
static const float AEROSOL_HEIGHT_SCALE = 0.73f;
|
|
/* Spectral to XYZ space conversion matrix. */
|
|
static const float3 SPECTRAL_XYZ[4] = {
|
|
make_float3(53.386917738564668023f, 22.981337506691024754f, 0.0f),
|
|
make_float3(43.904844466369358263f, 71.347795700053393866f, 0.102506867965741307f),
|
|
make_float3(1.6137278251608962005f, 18.422960591455485011f, 31.742921188390805758f),
|
|
make_float3(20.762668673810577145f, 2.3614213523314368527f, 110.48009643252140334f),
|
|
};
|
|
|
|
inline float molecular_phase_function(const float cos_theta)
|
|
{
|
|
return RAYLEIGH_PHASE_SCALE * (1.0f + sqr(cos_theta));
|
|
}
|
|
|
|
inline float aerosol_phase_function(const float cos_theta)
|
|
{
|
|
float den = 1.0f + SQR_G + 2.0f * G * cos_theta;
|
|
return M_1_4PI_F * (1.0f - SQR_G) / (den * sqrtf(den));
|
|
}
|
|
|
|
inline float4 get_molecular_scattering_coefficient(const float h)
|
|
{
|
|
return MOLECULAR_SCATTERING_COEFFICIENT_BASE * expf(-0.07771971f * powf(h, 1.16364243f));
|
|
}
|
|
|
|
inline float4 get_molecular_absorption_coefficient(const float h)
|
|
{
|
|
const float log_h = logf(fmaxf(h, 1e-4f));
|
|
float density = 3.78547397e20f * expf(-sqr(log_h - 3.22261f) * 5.55555555f - log_h);
|
|
return OZONE_ABSORPTION_CROSS_SECTION * OZONE_MEAN_DOBSON * density;
|
|
}
|
|
|
|
inline float get_aerosol_density(const float h)
|
|
{
|
|
float division = AEROSOL_BACKGROUND_DENSITY / AEROSOL_BASE_DENSITY;
|
|
return AEROSOL_BASE_DENSITY * (expf(-h / AEROSOL_HEIGHT_SCALE) + division);
|
|
}
|
|
|
|
inline float3 spectral_to_xyz(const float4 L)
|
|
{
|
|
float3 xyz = make_float3(0.0f, 0.0f, 0.0f);
|
|
for (int i = 0; i < 4; i++) {
|
|
xyz += SPECTRAL_XYZ[i] * L[i];
|
|
}
|
|
return xyz;
|
|
}
|
|
|
|
/* Precomputed data. */
|
|
class SkyMultipleScattering {
|
|
public:
|
|
SkyMultipleScattering(const float air_density,
|
|
const float aerosol_density,
|
|
const float ozone_density)
|
|
: air_density(air_density), aerosol_density(aerosol_density), ozone_density(ozone_density)
|
|
{
|
|
}
|
|
|
|
/* Compute atmosphere's transmittance from the given altitude to the sun. */
|
|
inline float4 get_transmittance(const float cos_theta, const float normalized_altitude) const
|
|
{
|
|
const float3 sun_dir = sun_direction(cos_theta);
|
|
const float distance_to_earth_center = mix(
|
|
EARTH_RADIUS, ATMOSPHERE_RADIUS, normalized_altitude);
|
|
const float3 ray_origin = make_float3(0.0f, 0.0f, distance_to_earth_center);
|
|
const float t_d = ray_sphere_intersection(ray_origin, sun_dir, ATMOSPHERE_RADIUS);
|
|
const float t_step = t_d / TRANSMITTANCE_STEPS;
|
|
|
|
float4 result = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
|
|
for (int step = 0; step < TRANSMITTANCE_STEPS; step++) {
|
|
const float t = (step + 0.5f) * t_step;
|
|
const float3 x_t = ray_origin + sun_dir * t;
|
|
const float altitude = fmaxf(x_t.length() - EARTH_RADIUS, 0.0f);
|
|
float4 aerosol_absorption, aerosol_scattering, molecular_absorption, molecular_scattering;
|
|
get_atmosphere_collision_coefficients(altitude,
|
|
aerosol_absorption,
|
|
aerosol_scattering,
|
|
molecular_absorption,
|
|
molecular_scattering);
|
|
const float4 extinction = aerosol_absorption + aerosol_scattering + molecular_absorption +
|
|
molecular_scattering;
|
|
result += extinction * t_step;
|
|
}
|
|
|
|
return exp(-result);
|
|
}
|
|
|
|
/* Compute in-scattered radiance for the given ray. */
|
|
float4 get_inscattering(const float3 sun_dir,
|
|
const float3 ray_origin,
|
|
const float3 ray_dir,
|
|
const float t_d) const
|
|
{
|
|
const float cos_theta = dot(-ray_dir, sun_dir);
|
|
const float molecular_phase = molecular_phase_function(cos_theta);
|
|
const float aerosol_phase = aerosol_phase_function(cos_theta);
|
|
const float dt = t_d / IN_SCATTERING_STEPS;
|
|
float4 L_inscattering = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
|
|
float4 transmittance = make_float4(1.0f, 1.0f, 1.0f, 1.0f);
|
|
for (int i = 0; i < IN_SCATTERING_STEPS; i++) {
|
|
const float t = (i + 0.5f) * dt;
|
|
const float3 x_t = ray_origin + ray_dir * t;
|
|
const float distance_to_earth_center = x_t.length();
|
|
const float3 zenith_dir = x_t / distance_to_earth_center;
|
|
const float altitude = fmaxf(distance_to_earth_center - EARTH_RADIUS, 0.0f);
|
|
const float normalized_altitude = altitude / ATMOSPHERE_THICKNESS;
|
|
const float sample_cos_theta = dot(zenith_dir, sun_dir);
|
|
float4 aerosol_absorption, aerosol_scattering, molecular_absorption, molecular_scattering;
|
|
get_atmosphere_collision_coefficients(altitude,
|
|
aerosol_absorption,
|
|
aerosol_scattering,
|
|
molecular_absorption,
|
|
molecular_scattering);
|
|
const float4 extinction = aerosol_absorption + aerosol_scattering + molecular_absorption +
|
|
molecular_scattering;
|
|
const float4 transmittance_to_sun = lookup_transmittance(sample_cos_theta,
|
|
normalized_altitude);
|
|
const float4 ms = lookup_multiscattering(
|
|
sample_cos_theta, normalized_altitude, distance_to_earth_center);
|
|
const float4 S = SUN_SPECTRAL_IRRADIANCE *
|
|
(molecular_scattering * (molecular_phase * transmittance_to_sun + ms) +
|
|
aerosol_scattering * (aerosol_phase * transmittance_to_sun + ms));
|
|
const float4 step_transmittance = exp(-dt * extinction);
|
|
/* Energy-conserving analytical integration "Physically Based Sky, Atmosphere and Cloud
|
|
* Rendering in Frostbite" by Sébastien Hillaire. */
|
|
const float4 cut_ext = max(extinction, 1e-7f);
|
|
const float4 S_int = (S - S * step_transmittance) / cut_ext;
|
|
L_inscattering = L_inscattering + transmittance * S_int;
|
|
transmittance *= step_transmittance;
|
|
}
|
|
|
|
return L_inscattering;
|
|
}
|
|
|
|
/* Precompute the transmittance LUT. Must be called before get_inscattering(). */
|
|
void precompute_lut()
|
|
{
|
|
SKY_parallel_for(0, TRANSMITTANCE_RES_Y, 4, [&](const size_t begin, const size_t end) {
|
|
for (int y = begin; y < end; y++) {
|
|
for (int x = 0; x < TRANSMITTANCE_RES_X; x++) {
|
|
const float2 uv = make_float2(x / float(TRANSMITTANCE_RES_X - 1),
|
|
y / float(TRANSMITTANCE_RES_Y - 1));
|
|
transmittance_lut[y][x] = get_transmittance(uv.x * 2.0f - 1.0f, uv.y);
|
|
}
|
|
}
|
|
});
|
|
}
|
|
|
|
protected:
|
|
float4 transmittance_lut[TRANSMITTANCE_RES_Y][TRANSMITTANCE_RES_X];
|
|
float air_density;
|
|
float aerosol_density;
|
|
float ozone_density;
|
|
|
|
/* Compute absorption/scattering coeffients at the given altitude. */
|
|
inline void get_atmosphere_collision_coefficients(const float altitude,
|
|
float4 &aerosol_absorption,
|
|
float4 &aerosol_scattering,
|
|
float4 &molecular_absorption,
|
|
float4 &molecular_scattering) const
|
|
{
|
|
const float local_aerosol_density = get_aerosol_density(altitude) * aerosol_density;
|
|
aerosol_absorption = AEROSOL_ABSORPTION_CROSS_SECTION * local_aerosol_density;
|
|
aerosol_scattering = AEROSOL_SCATTERING_CROSS_SECTION * local_aerosol_density;
|
|
molecular_absorption = get_molecular_absorption_coefficient(altitude) * ozone_density;
|
|
molecular_scattering = get_molecular_scattering_coefficient(altitude) * air_density;
|
|
}
|
|
|
|
inline float4 lookup_multiscattering(float cos_theta, float normalized_height, float d) const
|
|
{
|
|
/* Solid angle subtended by the planet from a point at d distance from the planet center. */
|
|
const float omega = M_2PI_F * (1.0f - sqrtf(1.0f - sqr(EARTH_RADIUS / d)));
|
|
const float4 T_to_ground = lookup_transmittance_at_ground(cos_theta);
|
|
/* We can split the path into Ground <-> Sample <-> Sun.
|
|
* The LUT gives us both T(Sample,Sun) and T(Ground,Sun) = T(Ground,Sample)*T(Sample,Sun),
|
|
* so we can easily compute T(Ground,Sample) from those two. */
|
|
const float4 T_ground_to_sample = lookup_transmittance_to_sun(0.0f) /
|
|
lookup_transmittance_to_sun(normalized_height);
|
|
/* 2nd order scattering from the ground. */
|
|
const float4 L_ground = PHASE_ISOTROPIC * omega * (GROUND_ALBEDO * M_1_PI_F) * T_to_ground *
|
|
T_ground_to_sample * cos_theta;
|
|
/* Fit of Earth's multiple scattering coming from other points in the atmosphere. */
|
|
const float4 L_ms = 0.02f * make_float4(0.217f, 0.347f, 0.594f, 1.0f) *
|
|
(1.0f / (1.0f + 5.0f * expf(-17.92f * cos_theta)));
|
|
return L_ms + L_ground;
|
|
}
|
|
|
|
/* Look up a transmittance from the precomputed LUT. */
|
|
inline float4 lookup_transmittance(const float cos_theta, const float normalized_altitude) const
|
|
{
|
|
const float u = saturate(cos_theta * 0.5f + 0.5f);
|
|
const float v = saturate(normalized_altitude);
|
|
const float x = float(TRANSMITTANCE_RES_X - 1) * u;
|
|
const float y = float(TRANSMITTANCE_RES_Y - 1) * v;
|
|
const int x1 = int(x);
|
|
const int y1 = int(y);
|
|
const int x2 = min(x1 + 1, TRANSMITTANCE_RES_X - 1);
|
|
const int y2 = min(y1 + 1, TRANSMITTANCE_RES_Y - 1);
|
|
const float fx = x - x1;
|
|
const float fy = y - y1;
|
|
const float4 bottom = mix(transmittance_lut[y1][x1], transmittance_lut[y1][x2], fx);
|
|
const float4 top = mix(transmittance_lut[y2][x1], transmittance_lut[y2][x2], fx);
|
|
return mix(bottom, top, fy);
|
|
}
|
|
|
|
/* Specialized versions of lookup_transmittance that skip one interpolation. */
|
|
inline float4 lookup_transmittance_at_ground(const float cos_theta) const
|
|
{
|
|
const float u = saturate(cos_theta * 0.5f + 0.5f);
|
|
const float x = float(TRANSMITTANCE_RES_X - 1) * u;
|
|
const int x1 = int(x);
|
|
const int x2 = min(x1 + 1, TRANSMITTANCE_RES_X - 1);
|
|
const int y = 0;
|
|
const float fx = x - x1;
|
|
return mix(transmittance_lut[y][x1], transmittance_lut[y][x2], fx);
|
|
}
|
|
|
|
inline float4 lookup_transmittance_to_sun(const float normalized_altitude) const
|
|
{
|
|
const float v = saturate(normalized_altitude);
|
|
const float y = float(TRANSMITTANCE_RES_Y - 1) * v;
|
|
const int x = TRANSMITTANCE_RES_X - 1;
|
|
const int y1 = int(y);
|
|
const int y2 = min(y1 + 1, TRANSMITTANCE_RES_Y - 1);
|
|
const float fy = y - y1;
|
|
return mix(transmittance_lut[y1][x], transmittance_lut[y2][x], fy);
|
|
}
|
|
};
|
|
|
|
void SKY_multiple_scattering_precompute_texture(float *pixels,
|
|
int stride,
|
|
int width,
|
|
int height,
|
|
float sun_elevation,
|
|
float altitude,
|
|
float air_density,
|
|
float aerosol_density,
|
|
float ozone_density)
|
|
{
|
|
SkyMultipleScattering sms(air_density, aerosol_density, ozone_density);
|
|
sms.precompute_lut();
|
|
|
|
/* Clamp altitude to avoid numerical issues. */
|
|
altitude = clamp(altitude, 1.0f, 99999.0f) / 1000.0f;
|
|
const int half_width = width / 2;
|
|
const float sun_zenith_cos_angle = cosf(M_PI_2_F - sun_elevation);
|
|
const float3 sun_dir = sun_direction(sun_zenith_cos_angle);
|
|
const int rows_per_task = std::max(1024 / width, 1);
|
|
|
|
SKY_parallel_for(0, height, rows_per_task, [=](const size_t begin, const size_t end) {
|
|
for (int y = begin; y < end; y++) {
|
|
float *pixel_row = pixels + (y * width * stride);
|
|
for (int x = 0; x < half_width; x++) {
|
|
float2 uv = make_float2((x + 0.5f) / width, (y + 0.5f) / height);
|
|
|
|
const float azimuth = M_2PI_F * uv.x;
|
|
/* Apply a non-linear transformation to the elevation to dedicate more texels to the
|
|
* horizon, where having more detail matters. */
|
|
const float l = uv.y * 2.0f - 1.0f;
|
|
/* [-pi/2, pi/2]. */
|
|
const float elev = copysignf(sqr(l), l) * M_PI_2_F;
|
|
const float3 ray_dir = make_float3(
|
|
cosf(elev) * cosf(azimuth), cosf(elev) * sinf(azimuth), sinf(elev));
|
|
const float3 ray_origin = make_float3(0.0f, 0.0f, EARTH_RADIUS + altitude);
|
|
const float atmos_dist = ray_sphere_intersection(ray_origin, ray_dir, ATMOSPHERE_RADIUS);
|
|
const float ground_dist = ray_sphere_intersection(ray_origin, ray_dir, EARTH_RADIUS);
|
|
/* If no ground collision then use the distance to the outer atmosphere, else we have a
|
|
* collision with the ground so we use the distance to it. */
|
|
const float t_d = (ground_dist < 0.0f) ? atmos_dist : ground_dist;
|
|
const float4 L = sms.get_inscattering(sun_dir, ray_origin, ray_dir, t_d);
|
|
const float3 sky = spectral_to_xyz(L);
|
|
|
|
/* Store pixels. */
|
|
const int pos_x = x * stride;
|
|
pixel_row[pos_x] = sky.x;
|
|
pixel_row[pos_x + 1] = sky.y;
|
|
pixel_row[pos_x + 2] = sky.z;
|
|
/* Mirror pixels. */
|
|
const int mirror_x = (width - x - 1) * stride;
|
|
pixel_row[mirror_x] = sky.x;
|
|
pixel_row[mirror_x + 1] = sky.y;
|
|
pixel_row[mirror_x + 2] = sky.z;
|
|
}
|
|
}
|
|
});
|
|
}
|
|
|
|
void SKY_multiple_scattering_precompute_sun(float sun_elevation,
|
|
float angular_diameter,
|
|
float altitude,
|
|
float air_density,
|
|
float aerosol_density,
|
|
float ozone_density,
|
|
float r_pixel_bottom[3],
|
|
float r_pixel_top[3])
|
|
{
|
|
const SkyMultipleScattering sms(air_density, aerosol_density, ozone_density);
|
|
|
|
/* Clamp altitude to avoid numerical issues. */
|
|
altitude = clamp(altitude, 1.0f, 99999.0f) / 1000.0f;
|
|
const float half_angular = angular_diameter / 2.0f;
|
|
const float solid_angle = M_2PI_F * (1.0f - cosf(half_angular));
|
|
const float normalized_altitude = altitude / ATMOSPHERE_THICKNESS;
|
|
|
|
/* Compute 2 pixels for Sun disc: one is the lowest point of the disc, one is the highest. */
|
|
auto get_sun_xyz = [&](const float elevation) {
|
|
const float sun_zenith_cos_angle = cosf(M_PI_2_F - elevation);
|
|
const float4 transmittance_to_sun = sms.get_transmittance(sun_zenith_cos_angle,
|
|
normalized_altitude);
|
|
const float4 spectrum = SUN_SPECTRAL_IRRADIANCE * transmittance_to_sun / solid_angle;
|
|
return spectral_to_xyz(spectrum);
|
|
};
|
|
const float3 bottom = get_sun_xyz(sun_elevation - half_angular);
|
|
const float3 top = get_sun_xyz(sun_elevation + half_angular);
|
|
|
|
/* Store pixels */
|
|
r_pixel_bottom[0] = bottom.x;
|
|
r_pixel_bottom[1] = bottom.y;
|
|
r_pixel_bottom[2] = bottom.z;
|
|
r_pixel_top[0] = top.x;
|
|
r_pixel_top[1] = top.y;
|
|
r_pixel_top[2] = top.z;
|
|
}
|
|
|
|
float SKY_earth_intersection_angle(float altitude)
|
|
{
|
|
/* Calculate intersection angle between line passing through viewpoint and Earth surface. */
|
|
return M_PI_2_F - asinf(EARTH_RADIUS / (EARTH_RADIUS + altitude / 1000.0f));
|
|
}
|