This GLSL fragment shader is used to render our demo scene, which consists of a sphere S on a purely spherical planet P. It is rendered by "ray tracing", i.e. the vertex shader outputs the view ray direction, and the fragment shader computes the intersection of this ray with the spheres S and P to produce the final pixels. The fragment shader also computes the intersection of the light rays with the sphere S, to compute shadows, as well as the intersections of the view ray with the shadow volume of S, in order to compute light shafts.
Our fragment shader has the following inputs and outputs:
uniform vec3 camera; uniform float exposure; uniform vec3 white_point; uniform vec3 earth_center; uniform vec3 sun_direction; uniform vec2 sun_size; in vec3 view_ray; layout(location = 0) out vec4 color;
It uses the following constants, as well as the following atmosphere
rendering functions, defined externally (by the Model's
GetShader() shader). The USE_LUMINANCE option is used
to select either the functions returning radiance values, or those returning
luminance values (see model.h).
const float PI = 3.14159265;
const vec3 kSphereCenter = vec3(0.0, 0.0, 1000.0) / kLengthUnitInMeters;
const float kSphereRadius = 1000.0 / kLengthUnitInMeters;
const vec3 kSphereAlbedo = vec3(0.8);
const vec3 kGroundAlbedo = vec3(0.0, 0.0, 0.04);
#ifdef USE_LUMINANCE
#define GetSolarRadiance GetSolarLuminance
#define GetSkyRadiance GetSkyLuminance
#define GetSkyRadianceToPoint GetSkyLuminanceToPoint
#define GetSunAndSkyIrradiance GetSunAndSkyIlluminance
#endif
vec3 GetSolarRadiance();
vec3 GetSkyRadiance(vec3 camera, vec3 view_ray, float shadow_length,
vec3 sun_direction, out vec3 transmittance);
vec3 GetSkyRadianceToPoint(vec3 camera, vec3 point, float shadow_length,
vec3 sun_direction, out vec3 transmittance);
vec3 GetSunAndSkyIrradiance(
vec3 p, vec3 normal, vec3 sun_direction, out vec3 sky_irradiance);
The functions to compute shadows and light shafts must be defined before we can use them in the main shader function, so we define them first. Testing if a point is in the shadow of the sphere S is equivalent to test if the corresponding light ray intersects the sphere, which is very simple to do. However, this is only valid for a punctual light source, which is not the case of the Sun. In the following function we compute an approximate (and biased) soft shadow by taking the angular size of the Sun into account:
float GetSunVisibility(vec3 point, vec3 sun_direction) {
vec3 p = point - kSphereCenter;
float p_dot_v = dot(p, sun_direction);
float p_dot_p = dot(p, p);
float ray_sphere_center_squared_distance = p_dot_p - p_dot_v * p_dot_v;
float distance_to_intersection = -p_dot_v - sqrt(
kSphereRadius * kSphereRadius - ray_sphere_center_squared_distance);
if (distance_to_intersection > 0.0) {
// Compute the distance between the view ray and the sphere, and the
// corresponding (tangent of the) subtended angle. Finally, use this to
// compute an approximate sun visibility.
float ray_sphere_distance =
kSphereRadius - sqrt(ray_sphere_center_squared_distance);
float ray_sphere_angular_distance = -ray_sphere_distance / p_dot_v;
return smoothstep(1.0, 0.0, ray_sphere_angular_distance / sun_size.x);
}
return 1.0;
}
The sphere also partially occludes the sky light, and we approximate this effect with an ambient occlusion factor. The ambient occlusion factor due to a sphere is given in Radiation View Factors (Isidoro Martinez, 1995). In the simple case where the sphere is fully visible, it is given by the following function:
float GetSkyVisibility(vec3 point) {
vec3 p = point - kSphereCenter;
float p_dot_p = dot(p, p);
return
1.0 + p.z / sqrt(p_dot_p) * kSphereRadius * kSphereRadius / p_dot_p;
}
To compute light shafts we need the intersections of the view ray with the shadow volume of the sphere S. Since the Sun is not a punctual light source this shadow volume is not a cylinder but a cone (for the umbra, plus another cone for the penumbra, but we ignore it here):
Noting, as in the above figure, $\bp$ the camera position, $\bv$ and $\bs$ the unit view ray and sun direction vectors and $R$ the sphere radius (supposed to be centered on the origin), the point at distance $d$ from the camera is $\bq=\bp+d\bv$. This point is at a distance $\delta=-\bq\cdot\bs$ from the sphere center along the umbra cone axis, and at a distance $r$ from this axis given by $r^2=\bq\cdot\bq-\delta^2$. Finally, at distance $\delta$ along the axis the umbra cone has radius $\rho=R-\delta\tan\alpha$, where $\alpha$ is the Sun's angular radius. The point at distance $d$ from the camera is on the shadow cone only if $r^2=\rho^2$, i.e. only if \begin{equation} (\bp+d\bv)\cdot(\bp+d\bv)-((\bp+d\bv)\cdot\bs)^2= (R+((\bp+d\bv)\cdot\bs)\tan\alpha)^2 \end{equation} Developping this gives a quadratic equation for $d$: \begin{equation} ad^2+2bd+c=0 \end{equation} where
void GetSphereShadowInOut(vec3 view_direction, vec3 sun_direction,
out float d_in, out float d_out) {
vec3 pos = camera - kSphereCenter;
float pos_dot_sun = dot(pos, sun_direction);
float view_dot_sun = dot(view_direction, sun_direction);
float k = sun_size.x;
float l = 1.0 + k * k;
float a = 1.0 - l * view_dot_sun * view_dot_sun;
float b = dot(pos, view_direction) - l * pos_dot_sun * view_dot_sun -
k * kSphereRadius * view_dot_sun;
float c = dot(pos, pos) - l * pos_dot_sun * pos_dot_sun -
2.0 * k * kSphereRadius * pos_dot_sun - kSphereRadius * kSphereRadius;
float discriminant = b * b - a * c;
if (discriminant > 0.0) {
d_in = max(0.0, (-b - sqrt(discriminant)) / a);
d_out = (-b + sqrt(discriminant)) / a;
// The values of d for which delta is equal to 0 and kSphereRadius / k.
float d_base = -pos_dot_sun / view_dot_sun;
float d_apex = -(pos_dot_sun + kSphereRadius / k) / view_dot_sun;
if (view_dot_sun > 0.0) {
d_in = max(d_in, d_apex);
d_out = a > 0.0 ? min(d_out, d_base) : d_base;
} else {
d_in = a > 0.0 ? max(d_in, d_base) : d_base;
d_out = min(d_out, d_apex);
}
} else {
d_in = 0.0;
d_out = 0.0;
}
}
Using these functions we can now implement the main shader function, which computes the radiance from the scene for a given view ray. This function first tests if the view ray intersects the sphere S. If so it computes the sun and sky light received by the sphere at the intersection point, combines this with the sphere BRDF and the aerial perspective between the camera and the sphere. It then does the same with the ground, i.e. with the planet sphere P, and then computes the sky radiance and transmittance. Finally, all these terms are composited together (an opacity is also computed for each object, using an approximate view cone - sphere intersection factor) to get the final radiance.
We start with the computation of the intersections of the view ray with the shadow volume of the sphere, because they are needed to get the aerial perspective for the sphere and the planet:
void main() {
// Normalized view direction vector.
vec3 view_direction = normalize(view_ray);
// Tangent of the angle subtended by this fragment.
float fragment_angular_size =
length(dFdx(view_ray) + dFdy(view_ray)) / length(view_ray);
float shadow_in;
float shadow_out;
GetSphereShadowInOut(view_direction, sun_direction, shadow_in, shadow_out);
// Hack to fade out light shafts when the Sun is very close to the horizon.
float lightshaft_fadein_hack = smoothstep(
0.02, 0.04, dot(normalize(camera - earth_center), sun_direction));
We then test whether the view ray intersects the sphere S or not. If it does,
we compute an approximate (and biased) opacity value, using the same
approximation as in GetSunVisibility:
// Compute the distance between the view ray line and the sphere center,
// and the distance between the camera and the intersection of the view
// ray with the sphere (or NaN if there is no intersection).
vec3 p = camera - kSphereCenter;
float p_dot_v = dot(p, view_direction);
float p_dot_p = dot(p, p);
float ray_sphere_center_squared_distance = p_dot_p - p_dot_v * p_dot_v;
float distance_to_intersection = -p_dot_v - sqrt(
kSphereRadius * kSphereRadius - ray_sphere_center_squared_distance);
// Compute the radiance reflected by the sphere, if the ray intersects it.
float sphere_alpha = 0.0;
vec3 sphere_radiance = vec3(0.0);
if (distance_to_intersection > 0.0) {
// Compute the distance between the view ray and the sphere, and the
// corresponding (tangent of the) subtended angle. Finally, use this to
// compute the approximate analytic antialiasing factor sphere_alpha.
float ray_sphere_distance =
kSphereRadius - sqrt(ray_sphere_center_squared_distance);
float ray_sphere_angular_distance = -ray_sphere_distance / p_dot_v;
sphere_alpha =
min(ray_sphere_angular_distance / fragment_angular_size, 1.0);
We can then compute the intersection point and its normal, and use them to get the sun and sky irradiance received at this point. The reflected radiance follows, by multiplying the irradiance with the sphere BRDF:
vec3 point = camera + view_direction * distance_to_intersection;
vec3 normal = normalize(point - kSphereCenter);
// Compute the radiance reflected by the sphere.
vec3 sky_irradiance;
vec3 sun_irradiance = GetSunAndSkyIrradiance(
point - earth_center, normal, sun_direction, sky_irradiance);
sphere_radiance =
kSphereAlbedo * (1.0 / PI) * (sun_irradiance + sky_irradiance);
Finally, we take into account the aerial perspective between the camera and the sphere, which depends on the length of this segment which is in shadow:
float shadow_length =
max(0.0, min(shadow_out, distance_to_intersection) - shadow_in) *
lightshaft_fadein_hack;
vec3 transmittance;
vec3 in_scatter = GetSkyRadianceToPoint(camera - earth_center,
point - earth_center, shadow_length, sun_direction, transmittance);
sphere_radiance = sphere_radiance * transmittance + in_scatter;
}
In the following we repeat the same steps as above, but for the planet sphere P instead of the sphere S (a smooth opacity is not really needed here, so we don't compute it. Note also how we modulate the sun and sky irradiance received on the ground by the sun and sky visibility factors):
// Compute the distance between the view ray line and the Earth center,
// and the distance between the camera and the intersection of the view
// ray with the ground (or NaN if there is no intersection).
p = camera - earth_center;
p_dot_v = dot(p, view_direction);
p_dot_p = dot(p, p);
float ray_earth_center_squared_distance = p_dot_p - p_dot_v * p_dot_v;
distance_to_intersection = -p_dot_v - sqrt(
earth_center.z * earth_center.z - ray_earth_center_squared_distance);
// Compute the radiance reflected by the ground, if the ray intersects it.
float ground_alpha = 0.0;
vec3 ground_radiance = vec3(0.0);
if (distance_to_intersection > 0.0) {
vec3 point = camera + view_direction * distance_to_intersection;
vec3 normal = normalize(point - earth_center);
// Compute the radiance reflected by the ground.
vec3 sky_irradiance;
vec3 sun_irradiance = GetSunAndSkyIrradiance(
point - earth_center, normal, sun_direction, sky_irradiance);
ground_radiance = kGroundAlbedo * (1.0 / PI) * (
sun_irradiance * GetSunVisibility(point, sun_direction) +
sky_irradiance * GetSkyVisibility(point));
float shadow_length =
max(0.0, min(shadow_out, distance_to_intersection) - shadow_in) *
lightshaft_fadein_hack;
vec3 transmittance;
vec3 in_scatter = GetSkyRadianceToPoint(camera - earth_center,
point - earth_center, shadow_length, sun_direction, transmittance);
ground_radiance = ground_radiance * transmittance + in_scatter;
ground_alpha = 1.0;
}
Finally, we compute the radiance and transmittance of the sky, and composite together, from back to front, the radiance and opacities of all the objects of the scene:
// Compute the radiance of the sky.
float shadow_length = max(0.0, shadow_out - shadow_in) *
lightshaft_fadein_hack;
vec3 transmittance;
vec3 radiance = GetSkyRadiance(
camera - earth_center, view_direction, shadow_length, sun_direction,
transmittance);
// If the view ray intersects the Sun, add the Sun radiance.
if (dot(view_direction, sun_direction) > sun_size.y) {
radiance = radiance + transmittance * GetSolarRadiance();
}
radiance = mix(radiance, ground_radiance, ground_alpha);
radiance = mix(radiance, sphere_radiance, sphere_alpha);
color.rgb =
pow(vec3(1.0) - exp(-radiance / white_point * exposure), vec3(1.0 / 2.2));
color.a = 1.0;
}