Page MenuHome
Differential D8091 Diff 26670 intern/sky/source/sky_nishita.cpp

Changeset View

Standalone View

View Options

intern/sky/source/sky_nishita.cpp

/* /*
* Copyright 2011-2020 Blender Foundation * Copyright 2011-2020 Blender Foundation
* *
* Licensed under the Apache License, Version 2.0 (the "License"); * Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License. * you may not use this file except in compliance with the License.
* You may obtain a copy of the License at * You may obtain a copy of the License at
* *
* http://www.apache.org/licenses/LICENSE-2.0 * http://www.apache.org/licenses/LICENSE-2.0
* *
* Unless required by applicable law or agreed to in writing, software * Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS, * distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and * See the License for the specific language governing permissions and
* limitations under the License. * limitations under the License.
*/ */
#include "sky_model.h"
#include "sky_float3.h" #include "sky_float3.h"
#include "sky_model.h"
/* Constants */ /* Constants */
static const float rayleigh_scale = 8000.0f; // Rayleigh scale height (m) static const float rayleigh_scale = 8000.0f; // Rayleigh scale height (m)
static const float mie_scale = 1200.0f; // Mie scale height (m) static const float mie_scale = 1200.0f; // Mie scale height (m)
static const float mie_coeff = 2e-5f; // Mie scattering coefficient static const float mie_coeff = 2e-5f; // Mie scattering coefficient
static const float mie_G = 0.76f; // aerosols anisotropy static const float mie_G = 0.76f; // aerosols anisotropy
static const float sqr_G = mie_G * mie_G; // squared aerosols anisotropy
static const float earth_radius = 6360000.0f; // radius of Earth (m) static const float earth_radius = 6360000.0f; // radius of Earth (m)
static const float atmosphere_radius = 6420000.0f; // radius of atmosphere (m) static const float atmosphere_radius = 6420000.0f; // radius of atmosphere (m)
static const int steps = 32; // segments per primary ray static const int steps = 32; // segments per primary ray
static const int steps_light = 16; // segments per sun connection ray static const int steps_light = 16; // segments per sun connection ray
static const int num_wavelengths = 21; // number of wavelengths static const int num_wavelengths = 21; // number of wavelengths
static const int max_luminous_efficacy = 683; // maximum luminous efficacy
static const float step_lambda = (num_wavelengths - 1) *
1e-9f; // step between each sampled wavelength
/* irradiance at top of atmosphere */ /* irradiance at top of atmosphere */
static const float irradiance[] = { static const float irradiance[] = {
1.45756829855592995315f, 1.56596305559738380175f, 1.65148449067670455293f, 1.45756829855592995315f, 1.56596305559738380175f, 1.65148449067670455293f,
1.71496242737209314555f, 1.75797983805020541226f, 1.78256407885924539336f, 1.71496242737209314555f, 1.75797983805020541226f, 1.78256407885924539336f,
1.79095108475838560302f, 1.78541550133410664714f, 1.76815554864306845317f, 1.79095108475838560302f, 1.78541550133410664714f, 1.76815554864306845317f,
1.74122069647250410362f, 1.70647127164943679389f, 1.66556087452739887134f, 1.74122069647250410362f, 1.70647127164943679389f, 1.66556087452739887134f,
1.61993437242451854274f, 1.57083597368892080581f, 1.51932335059305478886f, 1.61993437242451854274f, 1.57083597368892080581f, 1.51932335059305478886f,
1.46628494965214395407f, 1.41245852740172450623f, 1.35844961970384092709f, 1.46628494965214395407f, 1.41245852740172450623f, 1.35844961970384092709f,
▲ Show 20 Lines • Show All 47 Lines • ▼ Show 20 Lines
static float3 spec_to_xyz(float *spectrum) static float3 spec_to_xyz(float *spectrum)
{ {
float3 xyz = make_float3(0.0f, 0.0f, 0.0f); float3 xyz = make_float3(0.0f, 0.0f, 0.0f);
for (int i = 0; i < num_wavelengths; i++) { for (int i = 0; i < num_wavelengths; i++) {
xyz.x += cmf_xyz[i][0] * spectrum[i]; xyz.x += cmf_xyz[i][0] * spectrum[i];
xyz.y += cmf_xyz[i][1] * spectrum[i]; xyz.y += cmf_xyz[i][1] * spectrum[i];
xyz.z += cmf_xyz[i][2] * spectrum[i]; xyz.z += cmf_xyz[i][2] * spectrum[i];
} }
return xyz * (20 * 683 * 1e-9f); return xyz * step_lambda * max_luminous_efficacy;
} }
/* Atmosphere volume models */ /* Atmosphere volume models */
static float density_rayleigh(float height) static float density_rayleigh(float height)
{ {
return expf(-height / rayleigh_scale); return expf(-height / rayleigh_scale);
} }
Show All 15 Lines
static float phase_rayleigh(float mu) static float phase_rayleigh(float mu)
{ {
return 3.0f / (16.0f * M_PI_F) * (1.0f + sqr(mu)); return 3.0f / (16.0f * M_PI_F) * (1.0f + sqr(mu));
} }
static float phase_mie(float mu) static float phase_mie(float mu)
{ {
static const float sqr_G = mie_G * mie_G;
return (3.0f * (1.0f - sqr_G) * (1.0f + sqr(mu))) / return (3.0f * (1.0f - sqr_G) * (1.0f + sqr(mu))) /
(8.0f * M_PI_F * (2.0f + sqr_G) * powf((1.0f + sqr_G - 2.0f * mie_G * mu), 1.5)); (8.0f * M_PI_F * (2.0f + sqr_G) * powf((1.0f + sqr_G - 2.0f * mie_G * mu), 1.5));
} }
/* Intersection helpers */ /* Intersection helpers */
static bool surface_intersection(float3 pos, float3 dir) static bool surface_intersection(float3 pos, float3 dir)
{ {
if (dir.z >= 0) if (dir.z >= 0)
Show All 27 Linesstatic float3 ray_optical_depth(float3 ray_origin, float3 ray_dir)
/* Instead of tracking the transmission spectrum across all wavelengths directly, /* Instead of tracking the transmission spectrum across all wavelengths directly,
* we use the fact that the density always has the same spectrum for each type of * we use the fact that the density always has the same spectrum for each type of
* scattering, so we split the density into a constant spectrum and a factor and * scattering, so we split the density into a constant spectrum and a factor and
* only track the factors. */ * only track the factors. */
float3 optical_depth = make_float3(0.0f, 0.0f, 0.0f); float3 optical_depth = make_float3(0.0f, 0.0f, 0.0f);
/* The density of each segment is evaluated at its middle. */ /* The density of each segment is evaluated at its middle. */
float3 P = ray_origin + 0.5f * segment; float3 P = ray_origin + 0.5f * segment;
for (int i = 0; i < steps_light; i++) { for (int i = 0; i < steps_light; i++) {
/* Compute height above sea level. */ /* Compute height above sea level. */
float height = len(P) - earth_radius; float height = len(P) - earth_radius;
/* Accumulate optical depth of this segment (density is assumed to be constant along it). */ /* Accumulate optical depth of this segment (density is assumed to be constant along it). */
float3 density = make_float3( float3 density = make_float3(
density_rayleigh(height), density_mie(height), density_ozone(height)); density_rayleigh(height), density_mie(height), density_ozone(height));
optical_depth += segment_length * density; optical_depth += density;
/* Advance along ray. */ /* Advance along ray. */
P += segment; P += segment;
} }
return optical_depth; return optical_depth * segment_length;
} }
/* Single Scattering implementation */ /* Single Scattering implementation */
static void single_scattering(float3 ray_dir, static void single_scattering(float3 ray_dir,
float3 sun_dir, float3 sun_dir,
float3 ray_origin, float3 ray_origin,
float air_density, float air_density,
float dust_density, float dust_density,
Show All 22 Linesstatic void single_scattering(float3 ray_dir,
/* Compute phase function for scattering and the density scale factor. */ /* Compute phase function for scattering and the density scale factor. */
float mu = dot(ray_dir, sun_dir); float mu = dot(ray_dir, sun_dir);
float3 phase_function = make_float3(phase_rayleigh(mu), phase_mie(mu), 0.0f); float3 phase_function = make_float3(phase_rayleigh(mu), phase_mie(mu), 0.0f);
float3 density_scale = make_float3(air_density, dust_density, ozone_density); float3 density_scale = make_float3(air_density, dust_density, ozone_density);
/* The density and in-scattering of each segment is evaluated at its middle. */ /* The density and in-scattering of each segment is evaluated at its middle. */
float3 P = ray_origin + 0.5f * segment; float3 P = ray_origin + 0.5f * segment;
for (int i = 0; i < steps; i++) { for (int i = 0; i < steps; i++) {
/* Compute height above sea level. */ /* Compute height above sea level. */
float height = len(P) - earth_radius; float height = len(P) - earth_radius;
/* Evaluate and accumulate optical depth along the ray. */ /* Evaluate and accumulate optical depth along the ray. */
float3 density = density_scale * make_float3(density_rayleigh(height), float3 density = density_scale * make_float3(density_rayleigh(height),
density_mie(height), density_mie(height),
density_ozone(height)); density_ozone(height));
Show All 33 Lines for (int i = 0; i < steps; i++) {
/* Advance along ray. */ /* Advance along ray. */
P += segment; P += segment;
} }
} }
/* calculate texture array */ /* calculate texture array */
void SKY_nishita_skymodel_precompute_texture(float *pixels, void SKY_nishita_skymodel_precompute_texture(float *pixels,
int stride, int stride,
int start_y, int start_y,
int end_y, int end_y,
int width, int width,
int height, int height,
float sun_elevation, float sun_elevation,
float altitude, float altitude,
float air_density, float air_density,
float dust_density, float dust_density,
float ozone_density) float ozone_density)
{ {
/* calculate texture pixels */ /* calculate texture pixels */
float spectrum[num_wavelengths]; float spectrum[num_wavelengths];
int half_width = width / 2; int half_width = width / 2;
float3 cam_pos = make_float3(0, 0, earth_radius + altitude); float3 cam_pos = make_float3(0, 0, earth_radius + altitude);
float3 sun_dir = geographical_to_direction(sun_elevation, 0.0f); float3 sun_dir = geographical_to_direction(sun_elevation, 0.0f);
float latitude_step = M_PI_2_F / height; float latitude_step = M_PI_2_F / height;
float longitude_step = M_2PI_F / width; float longitude_step = M_2PI_F / width;
float half_lat_step = latitude_step / 2.0f;
for (int y = start_y; y < end_y; y++) { for (int y = start_y; y < end_y; y++) {
float latitude = latitude_step * y; /* sample more pixels toward the horizon */
float latitude = (M_PI_2_F + half_lat_step) * sqr((float)y / height);
float *pixel_row = pixels + (y * width) * stride; float *pixel_row = pixels + (y * width * stride);
for (int x = 0; x < half_width; x++) { for (int x = 0; x < half_width; x++) {
float longitude = longitude_step * x - M_PI_F; float longitude = longitude_step * x - M_PI_F;
float3 dir = geographical_to_direction(latitude, longitude); float3 dir = geographical_to_direction(latitude, longitude);
single_scattering(dir, sun_dir, cam_pos, air_density, dust_density, ozone_density, spectrum); single_scattering(dir, sun_dir, cam_pos, air_density, dust_density, ozone_density, spectrum);
float3 xyz = spec_to_xyz(spectrum); float3 xyz = spec_to_xyz(spectrum);
pixel_row[x * stride + 0] = xyz.x; int pos_x = x * stride;
pixel_row[x * stride + 1] = xyz.y; pixel_row[pos_x] = xyz.x;
pixel_row[x * stride + 2] = xyz.z; pixel_row[pos_x + 1] = xyz.y;
int mirror_x = width - x - 1; pixel_row[pos_x + 2] = xyz.z;
pixel_row[mirror_x * stride + 0] = xyz.x; /* mirror sky */
pixel_row[mirror_x * stride + 1] = xyz.y; int mirror_x = (width - x - 1) * stride;
pixel_row[mirror_x * stride + 2] = xyz.z; pixel_row[mirror_x] = xyz.x;
pixel_row[mirror_x + 1] = xyz.y;
pixel_row[mirror_x + 2] = xyz.z;
} }
} }
} }
/* Sun disc */ /* Sun disc */
static void sun_radiation(float3 cam_dir, static void sun_radiation(float3 cam_dir,
float altitude, float altitude,
float air_density, float air_density,
float dust_density, float dust_density,
float solid_angle, float solid_angle,
float *r_spectrum) float *r_spectrum)
{ {
float3 cam_pos = make_float3(0, 0, earth_radius + altitude); float3 cam_pos = make_float3(0, 0, earth_radius + altitude);
float3 optical_depth = ray_optical_depth(cam_pos, cam_dir); float3 optical_depth = ray_optical_depth(cam_pos, cam_dir);
/* Compute final spectrum. */ /* Compute final spectrum. */
for (int i = 0; i < num_wavelengths; i++) { for (int i = 0; i < num_wavelengths; i++) {
/* Combine spectra and the optical depth into transmittance. */ /* Combine spectra and the optical depth into transmittance. */
float transmittance = rayleigh_coeff[i] * optical_depth.x * air_density + float transmittance = rayleigh_coeff[i] * optical_depth.x * air_density +
1.11f * mie_coeff * optical_depth.y * dust_density; 1.11f * mie_coeff * optical_depth.y * dust_density;
r_spectrum[i] = (irradiance[i] / solid_angle) * expf(-transmittance); r_spectrum[i] = irradiance[i] * expf(-transmittance) / solid_angle;
} }
} }
void SKY_nishita_skymodel_precompute_sun(float sun_elevation, void SKY_nishita_skymodel_precompute_sun(float sun_elevation,
float angular_diameter, float angular_diameter,
float altitude, float altitude,
float air_density, float air_density,
float dust_density, float dust_density,
float *pixel_bottom, float *r_pixel_bottom,
float *pixel_top) float *r_pixel_top)
{ {
/* definitions */ /* definitions */
float half_angular = angular_diameter / 2.0f; float half_angular = angular_diameter / 2.0f;
float solid_angle = M_2PI_F * (1.0f - cosf(half_angular)); float solid_angle = M_2PI_F * (1.0f - cosf(half_angular));
float spectrum[num_wavelengths]; float spectrum[num_wavelengths];
float bottom = sun_elevation - half_angular; float bottom = sun_elevation - half_angular;
float top = sun_elevation + half_angular; float top = sun_elevation + half_angular;
float elevation_bottom, elevation_top; float elevation_bottom, elevation_top;
float3 pix_bottom, pix_top, sun_dir; float3 pix_bottom, pix_top, sun_dir;
/* compute 2 pixels for sun disc */ /* compute 2 pixels for sun disc */
elevation_bottom = (bottom > 0.0f) ? bottom : 0.0f; elevation_bottom = (bottom > 0.0f) ? bottom : 0.0f;
elevation_top = (top > 0.0f) ? top : 0.0f; elevation_top = (top > 0.0f) ? top : 0.0f;
sun_dir = geographical_to_direction(elevation_bottom, 0.0f); sun_dir = geographical_to_direction(elevation_bottom, 0.0f);
sun_radiation(sun_dir, altitude, air_density, dust_density, solid_angle, spectrum); sun_radiation(sun_dir, altitude, air_density, dust_density, solid_angle, spectrum);
pix_bottom = spec_to_xyz(spectrum); pix_bottom = spec_to_xyz(spectrum);
sun_dir = geographical_to_direction(elevation_top, 0.0f); sun_dir = geographical_to_direction(elevation_top, 0.0f);
sun_radiation(sun_dir, altitude, air_density, dust_density, solid_angle, spectrum); sun_radiation(sun_dir, altitude, air_density, dust_density, solid_angle, spectrum);
pix_top = spec_to_xyz(spectrum); pix_top = spec_to_xyz(spectrum);
/* store pixels */ /* store pixels */
pixel_bottom[0] = pix_bottom.x; r_pixel_bottom[0] = pix_bottom.x;
pixel_bottom[1] = pix_bottom.y; r_pixel_bottom[1] = pix_bottom.y;
pixel_bottom[2] = pix_bottom.z; r_pixel_bottom[2] = pix_bottom.z;
pixel_top[0] = pix_top.x; r_pixel_top[0] = pix_top.x;
pixel_top[1] = pix_top.y; r_pixel_top[1] = pix_top.y;
pixel_top[2] = pix_top.z; r_pixel_top[2] = pix_top.z;
} }