diff --git a/src/atmosphere.cc b/src/atmosphere.cc index 1851383..b5587f0 100644 --- a/src/atmosphere.cc +++ b/src/atmosphere.cc @@ -83,6 +83,32 @@ void Atmosphere::PrecomputeTable() { NormalizeTable(); } + + // Compute transmittance + if (precomputed_transmittance_table_ != nullptr) + { + delete [] precomputed_transmittance_table_; + } + precomputed_transmittance_table_ = new double[table_dimension_*3]; + + for (int i = 0; i < table_dimension_*3; i += 3) + { + int j = i/3; + double u = (j+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 transmittance = 0.0; + + PrecomputeTransmittance(view_dir, lambda_r_, transmittance); + precomputed_transmittance_table_[i+0] = transmittance; + + PrecomputeTransmittance(view_dir, lambda_g_, transmittance); + precomputed_transmittance_table_[i+1] = transmittance; + + PrecomputeTransmittance(view_dir, lambda_b_, transmittance); + precomputed_transmittance_table_[i+2] = transmittance; + } } int Atmosphere::GetTableLength() @@ -100,6 +126,11 @@ double* Atmosphere::GetPrecomputedMieTable() return precomputed_mie_table_; } +double* Atmosphere::GetPrecomputedTransmittanceTable() +{ + return precomputed_transmittance_table_; +} + void Atmosphere::GetNormalizationFactorsRayleigh(double& min, double& max) { min = normalization_min_rayleigh_; @@ -252,6 +283,66 @@ void Atmosphere::PrecomputeTableCell(vec2 view_dir, vec2 light_dir, double lambd // the pixel shader, because they are constant anyways. } +void Atmosphere::PrecomputeTransmittance(vec2 view_dir, double lambda, double& transmittance) +{ + transmittance = 0.0; + + // `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 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)) + { + 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)) + { + transmittance = 0.0; + return; + } + + // Compute `pb`, the intersection of the `view_dir` ray with the atmosphere + // from `pa`. + vec2 pb = pa + (view_dir * dist_to_pb); + + 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 = (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_; + + 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; + } + + optical_depth_rayleigh *= Coefficients::GetRayleighExtinctionCoefficient(lambda); + optical_depth_mie *= Coefficients::GetMieExtinctionCoefficient(lambda); + optical_depth_ozone *= Coefficients::GetOzoneExtinctionCoefficient(lambda); + + double optical_depth = optical_depth_rayleigh + optical_depth_mie + optical_depth_ozone; + + transmittance = exp(-optical_depth); +} + void Atmosphere::NormalizeTable() { const int kTableLength = GetTableLength(); diff --git a/src/atmosphere.h b/src/atmosphere.h index 022bc87..eb58b04 100644 --- a/src/atmosphere.h +++ b/src/atmosphere.h @@ -43,6 +43,8 @@ class Atmosphere double* GetPrecomputedRayleighTable(); double* GetPrecomputedMieTable(); + double* GetPrecomputedTransmittanceTable(); + void GetNormalizationFactorsRayleigh(double& min, double& max); void GetNormalizationFactorsMie(double& min, double& max); @@ -51,6 +53,8 @@ class Atmosphere double* precomputed_rayleigh_table_ = nullptr; double* precomputed_mie_table_ = nullptr; + double* precomputed_transmittance_table_ = nullptr; + double normalization_min_rayleigh_; double normalization_max_rayleigh_; @@ -58,6 +62,7 @@ class Atmosphere double normalization_max_mie_; void PrecomputeTableCell(vec2 view_dir, vec2 light_dir, double lambda, double& rayleigh, double& mie); + void PrecomputeTransmittance(vec2 view_dir, double lambda, double& transmittance); void NormalizeTable(); diff --git a/src/main.cc b/src/main.cc index 521a67e..563af70 100644 --- a/src/main.cc +++ b/src/main.cc @@ -111,6 +111,29 @@ int main(int argc, char* argv[]) } delete [] pixels; + + // Write transmittance to file. + pixels = new Imf::Rgba[DIM]; + double* transmittance = atmosphere.GetPrecomputedTransmittanceTable(); + char transmittance_path[output_path.size() + 17 + 1]; + strcpy(transmittance_path, (output_path + "transmittance.exr").c_str()); + for (int i = 0; i < DIM*3; i += 3) + { + int j = i / 3; + pixels[j] = Imf::Rgba(transmittance[i+0], transmittance[i+1], transmittance[i+2], 1.0); + } + try + { + Imf::RgbaOutputFile transmittance_file(transmittance_path, DIM, 1, Imf::WRITE_RGBA); + transmittance_file.setFrameBuffer(pixels, 1, DIM); + transmittance_file.writePixels(1); + } + catch (const std::exception &exc) + { + (void)exc; + std::cout << "ERROR: failed to write transmittance table to .exr file." << std::endl; + } + delete [] pixels; } else {