Skip to content

Commit

Permalink
fixed a bug related to light ray - planet intersections not being con…
Browse files Browse the repository at this point in the history 
…sidered
  • Loading branch information
Daniel Shervheim committed Jul 16, 2020
1 parent 63ef956 commit f61ec8d
Show file tree
Hide file tree
Showing 2 changed files with 94 additions and 45 deletions.
133 changes: 88 additions & 45 deletions src/atmosphere.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@

#include <cmath> // acos
#include <iostream> // cout
#include <stdexcept> // runtime_error

Atmosphere::Atmosphere()
{
Expand Down Expand Up @@ -54,12 +55,12 @@ void Atmosphere::PrecomputeTable()
// For both, we use sin() for x and cos() for y, since we are
// measuring the angle from the zenith, rather than the horizon.

double u = x / (double)(table_dimension_ - 1);
double theta_view_dir = acos(-pow(2.0 * u - 1.0, 3.0));
double u = (x+0.5) / (double)(table_dimension_);
double theta_view_dir = acos(TextureCoordinateToCosTheta(u));
vec2 view_dir = vec2(sin(theta_view_dir), cos(theta_view_dir));

double v = y / (double)(table_dimension_ - 1);
double theta_light_dir = acos(-pow(2.0 * v - 1.0, 3.0));
double v = (y+0.5) / (double)(table_dimension_);
double theta_light_dir = acos(TextureCoordinateToCosTheta(v));
vec2 light_dir = vec2(sin(theta_light_dir), cos(theta_light_dir));

double rayleigh = 0.0;
Expand Down Expand Up @@ -142,85 +143,106 @@ void Atmosphere::PrecomputeTableCell(vec2 view_dir, vec2 light_dir, double lambd
rayleigh = 0.0;
mie = 0.0;

// Compute `pa`, the position of an observer in the atmosphere. We use
// the height of an average human as the offset from the surface of the
// planet.
// `pa`, the viewing position. We offset the planet radius by the height
// of an average human to get a typical viewing position.
vec2 pa = vec2(0.0, radius_planet_ + 1.8);

double t_min;
double t_max;
double dist_to_pb;
double t_min, t_max;

// Intersect the viewing ray with the atmosphere.
if (!Utilities::RayCircleIntersection(pa, view_dir, vec2(0.0, 0.0), radius_atmosphere_, t_min, t_max))
{
// TODO: error.
std::cerr << "pa outside the atmosphere. ensure that the atmosphere radius is greater than the planet radius." << std::endl;
}

dist_to_pb = t_max;

// Intersect the viewing ray with the ground. This ensures that we do not
// integrate through the planet to the atmosphere on the opposite side.
if (Utilities::RayCircleIntersection(pa, view_dir, vec2(0.0, 0.0), radius_planet_, t_min, t_max))
{
dist_to_pb = fmin(dist_to_pb, t_min);
}

// Compute `pb`, the intersection of the `view_dir` ray with the atmosphere
// from `pa`.
vec2 pb = pa + (view_dir * t_max);
vec2 pb = pa + (view_dir * dist_to_pb);

double kDeltaView = vec2::magnitude(pb - pa) / (double)view_ray_integration_steps_;
const double kViewDelta = vec2::magnitude(pb - pa) / (double)view_ray_integration_steps_;

double optical_depth_rayleigh = 0.0;
double optical_depth_mie = 0.0;
double optical_depth_ozone = 0.0;

// Integrate from `pa` to `pb`.
for (int i = 0; i < view_ray_integration_steps_; i++)
{
double t = i / (double)(view_ray_integration_steps_ - 1);
double t = (double)i / (double)(view_ray_integration_steps_ - 1.0);
vec2 p = vec2::lerp(pa, pb, t);

// Compute the height of the point `p` from the surface of the planet.
double height_p = vec2::magnitude(p) - radius_planet_;

// TODO: is this ozone formula correct? Ozone concentration does not
// decrease with height...
optical_depth_rayleigh += exp(-height_p / scale_height_rayleigh_) * kDeltaView;
optical_depth_mie += exp(-height_p / scale_height_mie_) * kDeltaView;
optical_depth_ozone += exp(-height_p / scale_height_rayleigh_) * 6e-7 * kDeltaView;
// NOTE: the ozone density function is an approximation. In reality,
// ozone does not decrease exponentially with height.
optical_depth_rayleigh += exp(-height_p / scale_height_rayleigh_) * kViewDelta;
optical_depth_mie += exp(-height_p / scale_height_mie_) * kViewDelta;
optical_depth_ozone += exp(-height_p / scale_height_rayleigh_) * 6e-7 * kViewDelta;

// Compute the optical depth and transmittance from the `pa` to `p`.
double optical_depth_pa_to_p =
optical_depth_rayleigh * Coefficients::GetRayleighExtinctionCoefficient(lambda) +
optical_depth_mie * Coefficients::GetMieExtinctionCoefficient(lambda) +
optical_depth_ozone * Coefficients::GetOzoneExtinctionCoefficient(lambda);
double transmittance_pa_to_p = exp(-optical_depth_pa_to_p);

// Compute `pc`, the intersection of the `light_dir` ray with the atmosphere
// from `p`.
if (!Utilities::RayCircleIntersection(p, light_dir, vec2(0.0, 0.0), radius_atmosphere_, t_min, t_max))
// Intersect the light ray with the atmosphere. We add 1 meter to the
// atmosphere radius, otherwise the intersection test fails when p = pb).
if (!Utilities::RayCircleIntersection(p, light_dir, vec2(0.0, 0.0), radius_atmosphere_+1.0, t_min, t_max))
{
// TODO: error?
std::cerr << "p is outside the atmosphere." << std::endl;
}

vec2 pc = p + (light_dir * t_max);
double dist_to_pc = t_max;
double transmittance_p_to_pc = 0.0;

// Compute the transmittance between `p` and `pc`.
const double kDeltaLight = vec2::magnitude(pc - p) / (double)light_ray_integration_steps_;
double odr_p_to_pc = 0.0;
double odm_p_to_pc = 0.0;
double odo_p_to_pc = 0.0;

for (int j = 0; j < light_ray_integration_steps_; j++)
// Intersect the light ray with the planet. We only want to integrate
// from `p` to `pc` if the light ray is not blocked by the planet.
if (!Utilities::RayCircleIntersection(p, light_dir, vec2(0.0, 0.0), radius_planet_, t_min, t_max))
{
double t_2 = j / (double)(light_ray_integration_steps_ - 1);
vec2 p_2 = vec2::lerp(p, pc, t_2);
double height_p_2 = vec2::magnitude(p_2) - radius_planet_;

odr_p_to_pc += exp(-height_p_2 / scale_height_rayleigh_) * kDeltaLight;
odm_p_to_pc += exp(-height_p_2 / scale_height_mie_) * kDeltaLight;
odo_p_to_pc += exp(-height_p_2 / scale_height_rayleigh_) * 6e-7 * kDeltaLight;
// Compute `pc`, the intersection of the `light_dir` ray with the atmosphere
// from `p`.
vec2 pc = p + (light_dir * dist_to_pc);

// Compute the transmittance between `p` and `pc`.
const double kDeltaLight = vec2::magnitude(pc - p) / (double)light_ray_integration_steps_;
double odr_p_to_pc = 0.0;
double odm_p_to_pc = 0.0;
double odo_p_to_pc = 0.0;

for (int j = 0; j < light_ray_integration_steps_; j++)
{
double t_2 = (double)j / (double)(light_ray_integration_steps_ - 1.0);
vec2 p_2 = vec2::lerp(p, pc, t_2);
double height_p_2 = vec2::magnitude(p_2) - radius_planet_;

odr_p_to_pc += exp(-height_p_2 / scale_height_rayleigh_) * kDeltaLight;
odm_p_to_pc += exp(-height_p_2 / scale_height_mie_) * kDeltaLight;
odo_p_to_pc += exp(-height_p_2 / scale_height_rayleigh_) * 6e-7 * kDeltaLight;
}

double optical_depth_p_to_pc =
odr_p_to_pc * Coefficients::GetRayleighExtinctionCoefficient(lambda) +
odm_p_to_pc * Coefficients::GetMieExtinctionCoefficient(lambda) +
odo_p_to_pc * Coefficients::GetOzoneExtinctionCoefficient(lambda);
transmittance_p_to_pc = exp(-optical_depth_p_to_pc);
}

double optical_depth_p_to_pc =
odr_p_to_pc * Coefficients::GetRayleighExtinctionCoefficient(lambda) +
odm_p_to_pc * Coefficients::GetMieExtinctionCoefficient(lambda) +
odo_p_to_pc * Coefficients::GetOzoneExtinctionCoefficient(lambda);
double transmittance_p_to_pc = exp(-optical_depth_p_to_pc);

double transmittance_pa_to_p_to_pc = transmittance_pa_to_p * transmittance_p_to_pc;

rayleigh += exp(-height_p / scale_height_rayleigh_) * transmittance_pa_to_p_to_pc * kDeltaView;
mie += exp(-height_p / scale_height_mie_) * transmittance_pa_to_p_to_pc * kDeltaView;
rayleigh += exp(-height_p / scale_height_rayleigh_) * transmittance_pa_to_p_to_pc * kViewDelta;
mie += exp(-height_p / scale_height_mie_) * transmittance_pa_to_p_to_pc * kViewDelta;;
}

rayleigh *= Coefficients::GetRayleighScatteringCoefficient(lambda) / (4.0 * M_PI);
Expand Down Expand Up @@ -263,3 +285,24 @@ void Atmosphere::NormalizeTable()
normalization_min_mie_) / range_mie;
}
}

double Atmosphere::CosThetaToTextureCoordinate(double cosTheta)
{
// Favors precision near horizon
double sign = cosTheta < 0.0 ? -1.0 : 1.0;
return (sign*pow(abs(cosTheta), 1.0/texture_exponent_) + 1.0)/2.0;

// Linear remapping
// return cosTheta*0.5 + 0.5;
}

double Atmosphere::TextureCoordinateToCosTheta(double textureCoordinate)
{
// Favors precision near horizon
double remapped = textureCoordinate*2.0 - 1.0;
double sign = remapped < 0.0 ? -1.0 : 1.0;
return sign * pow(abs(remapped), texture_exponent_);

// Linear remapping
// return textureCoordinate*2.0 - 1.0;
}
6 changes: 6 additions & 0 deletions src/atmosphere.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,9 @@ class Atmosphere
double scale_height_rayleigh_ = 8000.0; // TODO: explain what these are.
double scale_height_mie_ = 1200.0;

// Exponent for lending more precision to horizon vs zeniths
double texture_exponent_ = 3.0;

Atmosphere();
~Atmosphere();

Expand All @@ -57,6 +60,9 @@ class Atmosphere
void PrecomputeTableCell(vec2 view_dir, vec2 light_dir, double lambda, double& rayleigh, double& mie);

void NormalizeTable();

double TextureCoordinateToCosTheta(double cosTheta);
double CosThetaToTextureCoordinate(double cosTheta);
};

#endif // ATMOSPHERE_H_

0 comments on commit f61ec8d

Please sign in to comment.