Skip to content

Commit

Permalink
Updated all volume integrators to correctly handle medium emitters an…
Browse files Browse the repository at this point in the history 
…d added appropriate python interfaces/documentation.
  • Loading branch information
@Microno95
Microno95 committed Apr 25, 2023
1 parent 2e0def0 commit e3a5a8c
Show file tree
Hide file tree
Showing 12 changed files with 185 additions and 61 deletions.
6 changes: 6 additions & 0 deletions include/mitsuba/python/docstr.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2339,6 +2339,8 @@ static const char *__doc_mitsuba_EmitterFlags_SpatiallyVarying = R"doc(The emiss

static const char *__doc_mitsuba_EmitterFlags_Surface = R"doc(The emitter is attached to a surface (e.g. area emitters))doc";

static const char *__doc_mitsuba_EmitterFlags_Medium = R"doc(The emitter is attached to a medium (e.g. medium emitters))doc";

static const char *__doc_mitsuba_Emitter_Emitter = R"doc()doc";

static const char *__doc_mitsuba_Emitter_class = R"doc()doc";
Expand Down Expand Up @@ -4105,6 +4107,8 @@ static const char *__doc_mitsuba_Medium_intersect_aabb = R"doc(Intersects a ray

static const char *__doc_mitsuba_Medium_is_homogeneous = R"doc(Returns whether this medium is homogeneous)doc";

static const char *__doc_mitsuba_Medium_is_emitter = R"doc(Returns whether this medium is an emitter)doc";

static const char *__doc_mitsuba_Medium_m_has_spectral_extinction = R"doc()doc";

static const char *__doc_mitsuba_Medium_m_id = R"doc(Identifier (if available))doc";
Expand All @@ -4125,6 +4129,8 @@ static const char *__doc_mitsuba_Medium_operator_new_2 = R"doc()doc";

static const char *__doc_mitsuba_Medium_phase_function = R"doc(Return the phase function of this medium)doc";

static const char *__doc_mitsuba_Medium_emitter = R"doc(Return the emitter associated with this medium)doc";

static const char *__doc_mitsuba_Medium_sample_interaction =
R"doc(Sample a free-flight distance in the medium.
Expand Down
2 changes: 1 addition & 1 deletion include/mitsuba/render/medium.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ NAMESPACE_BEGIN(mitsuba)
template <typename Float, typename Spectrum>
class MI_EXPORT_LIB Medium : public Object {
public:
MI_IMPORT_TYPES(PhaseFunction, Sampler, Scene, Texture, Emitter);
MI_IMPORT_TYPES(PhaseFunction, Sampler, Scene, Texture, Emitter, Volume);

/// Intersects a ray with the medium's bounding box
virtual std::tuple<Mask, Float, Float>
Expand Down
2 changes: 1 addition & 1 deletion src/emitters/volumelight.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,7 @@ class VolumeLight final : public Emitter<Float, Spectrum> {
m_flags = +EmitterFlags::Medium;

dr::set_attr(this, "flags", m_flags);
dr::set_attr(this, "radiance", m_radiance);
}

void traverse(TraversalCallback *callback) override {
Expand All @@ -82,7 +83,6 @@ class VolumeLight final : public Emitter<Float, Spectrum> {

Spectrum eval(const SurfaceInteraction3f &si, Mask active) const override {
MI_MASKED_FUNCTION(ProfilerPhase::EndpointEvaluate, active);

auto radiance = m_radiance->eval(si, active);

return depolarizer<Spectrum>(radiance);
Expand Down
2 changes: 1 addition & 1 deletion src/integrators/ptracer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,7 @@ class ParticleTracerIntegrator final : public AdjointIntegrator<Float, Spectrum>

// 3.a. Infinite emitters
Mask is_infinite = has_flag(emitter->flags(), EmitterFlags::Infinite),
active_e = active && is_infinite;
active_e = active && is_infinite && !has_flag(emitter->flags(), EmitterFlags::Medium);
if (dr::any_or<true>(active_e)) {
/* Sample a direction toward an envmap emitter starting
from the center of the scene (the sensor is not part of the
Expand Down
83 changes: 50 additions & 33 deletions src/integrators/volpath.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,16 +79,16 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator<Float, Spectrum> {
enum VolumeEventProbabilityMode { Analogue, Maximum, Mean };

VolumetricPathIntegrator(const Properties &props) : Base(props) {
std::string sampling_type = props.string("probability_type", "analogue");
if (sampling_type == "analogue") {
m_sampling_type = VolumeEventProbabilityMode::Analogue;
} else if (sampling_type == "max") {
m_sampling_type = VolumeEventProbabilityMode::Maximum;
} else if (sampling_type == "mean") {
m_sampling_type = VolumeEventProbabilityMode::Mean;
std::string sampling_mode = props.string("event_sampling_mode", "analogue");
if (sampling_mode == "analogue") {
m_sampling_mode = VolumeEventProbabilityMode::Analogue;
} else if (sampling_mode == "max") {
m_sampling_mode = VolumeEventProbabilityMode::Maximum;
} else if (sampling_mode == "mean") {
m_sampling_mode = VolumeEventProbabilityMode::Mean;
} else {
Log(Warn, "Sampling Probability Type \"%s\" not recognised, defaulting to \"analogue\" sampling", sampling_type);
m_sampling_type = VolumeEventProbabilityMode::Analogue;
Log(Warn, "Event Sampling Type \"%s\" not recognised, defaulting to \"analogue\" sampling", sampling_mode);
m_sampling_mode = VolumeEventProbabilityMode::Analogue;
}
}

Expand All @@ -113,7 +113,7 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator<Float, Spectrum> {
UInt32 &channel,
VolumeEventProbabilityMode probability_type = 0) const {
Float prob_emission, prob_scatter, prob_null, c;
Spectrum weight_emission(0.0f), weight_scatter(0.0f), weight_null(0.0f);
Float weight_emission(0.0f), weight_scatter(0.0f), weight_null(0.0f);

if (probability_type == VolumeEventProbabilityMode::Analogue) {
std::tie(prob_emission, prob_scatter, prob_null) = medium_probabilities_analog(radiance,
Expand All @@ -134,6 +134,8 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator<Float, Spectrum> {
Throw("Invalid probability type:", "%i", probability_type);
}

dr::masked(prob_emission, dr::max(radiance) > 0.f) = 1.0f;

c = prob_emission + prob_scatter + prob_null;

// {
Expand All @@ -154,13 +156,13 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator<Float, Spectrum> {
// Log(Info, "%s", oss.str());
// }

dr::masked(weight_emission, prob_emission > 0.f) = 1.f / prob_emission;
dr::masked(weight_scatter, prob_scatter > 0.f) = 1.f / prob_scatter;
dr::masked(weight_null, prob_null > 0.f) = 1.f / prob_null;
dr::masked(weight_emission, prob_emission > 0.f) = dr::rcp(prob_emission);
dr::masked(weight_scatter, prob_scatter > 0.f) = dr::rcp(prob_scatter);
dr::masked(weight_null, prob_null > 0.f) = dr::rcp(prob_null);

dr::masked(weight_emission, neq(weight_emission, weight_emission) || !(weight_emission < dr::Infinity<Float>)) = 0.f;
dr::masked(weight_scatter, neq(weight_scatter, weight_scatter) || !(weight_scatter < dr::Infinity<Float>)) = 0.f;
dr::masked(weight_null, neq(weight_null, weight_null) || !(weight_null < dr::Infinity<Float>)) = 0.f;
dr::masked(weight_emission, dr::neq(weight_emission, weight_emission) || !(weight_emission < dr::Infinity<Float>)) = 0.f;
dr::masked(weight_scatter, dr::neq(weight_scatter, weight_scatter) || !(weight_scatter < dr::Infinity<Float>)) = 0.f;
dr::masked(weight_null, dr::neq(weight_null, weight_null) || !(weight_null < dr::Infinity<Float>)) = 0.f;

return std::tuple{ std::tuple{prob_emission, prob_scatter, prob_null}, std::tuple{weight_emission, weight_scatter, weight_null} };
}
Expand All @@ -185,7 +187,7 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator<Float, Spectrum> {
UnpolarizedSpectrum /*sigma_t*/,
const UnpolarizedSpectrum &throughput) const {
Float prob_e = 0.f, prob_s = 0.f, prob_n = 0.f;
prob_e = dr::maximum(dr::max(dr::abs(radiance * throughput)), dr::max(dr::abs(radiance)));
prob_e = dr::max(dr::abs(radiance * dr::maximum(1.f, throughput)));
prob_s = dr::max(dr::abs( sigma_s * throughput));
prob_n = dr::max(dr::abs( sigma_n * throughput));
return { prob_e, prob_s, prob_n};
Expand All @@ -198,7 +200,7 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator<Float, Spectrum> {
UnpolarizedSpectrum /*sigma_t*/,
const UnpolarizedSpectrum &throughput) const {
Float prob_e = 0.f, prob_s = 0.f, prob_n = 0.f;
prob_e = (0.5f * dr::mean(dr::abs(radiance * throughput)) + 0.5f * dr::mean(dr::abs(radiance)));
prob_e = dr::mean(dr::abs(radiance * (0.5f + 0.5f * throughput)));
prob_s = dr::mean(dr::abs( sigma_s * throughput));
prob_n = dr::mean(dr::abs( sigma_n * throughput));
return {prob_e, prob_s, prob_n} ;
Expand Down Expand Up @@ -269,45 +271,59 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator<Float, Spectrum> {
Mask act_null_scatter = false, act_medium_scatter = false,
escaped_medium = false;

// If the medium does not have a spectrally varying extinction,
// we can perform a few optimizations to speed up rendering
Mask is_spectral = active_medium;
Mask not_spectral = false;
if (dr::any_or<true>(active_medium)) {
is_spectral &= medium->has_spectral_extinction();
not_spectral = !is_spectral && active_medium;
}

if (dr::any_or<true>(active_medium)) {
mei = medium->sample_interaction(ray, sampler->next_1d(active_medium), channel, active_medium);
dr::masked(ray.maxt, active_medium && medium->is_homogeneous() && mei.is_valid()) = mei.t;
Mask intersect = needs_intersection && active_medium;
if (dr::any_or<true>(intersect))
dr::masked(si, intersect) = scene->ray_intersect(ray, intersect);
needs_intersection &= !active_medium;

dr::masked(mei.t, active_medium && (si.t < mei.t)) = dr::Infinity<Float>;
auto [tr, free_flight_pdf] = medium->eval_tr_and_pdf(mei, si, active_medium);
Float tr_pdf = index_spectrum(free_flight_pdf, channel);

auto prev_throughput = throughput;
dr::masked(throughput, active_medium) *= dr::select(tr_pdf > 0.f, tr / tr_pdf, 0.f);
if (dr::any_or<true>(is_spectral)) {
auto [tr, free_flight_pdf] = medium->eval_tr_and_pdf(mei, si, active_medium);
Float tr_pdf = index_spectrum(free_flight_pdf, channel);
dr::masked(throughput, active_medium) *= dr::select(tr_pdf > 0.f, tr / tr_pdf, 0.f);
}

escaped_medium = active_medium && !mei.is_valid();
active_medium &= mei.is_valid();
is_spectral &= active_medium;
not_spectral &= active_medium;

// Compute emission, scatter and null event probabilities
auto radiance = medium->get_radiance(mei, active_medium);
auto [probabilities, weights] = medium_probabilities(
medium->get_radiance(mei, active_medium),
radiance,
unpolarized_spectrum(mei.sigma_s),
unpolarized_spectrum(mei.sigma_n),
unpolarized_spectrum(mei.sigma_t),
unpolarized_spectrum(prev_throughput),
channel,
m_sampling_type);
auto [prob_emission, prob_scatter, prob_null] = probabilities;
m_sampling_mode);
auto [emission_prob, scatter_prob, null_prob] = probabilities;
auto [weight_emission, weight_scatter, weight_null] = weights;

// Handle null and real scatter events
Mask null_scatter = sampler->next_1d(active_medium) < prob_emission + prob_null;
Mask null_scatter = sampler->next_1d(active_medium) < (emission_prob + null_prob);

act_null_scatter |= null_scatter && active_medium;
act_medium_scatter |= !null_scatter && active_medium;
act_null_scatter = null_scatter && active_medium;
act_medium_scatter = !null_scatter && active_medium;

if (dr::any_or<true>(act_null_scatter))
{
dr::masked(result, act_null_scatter) += weight_emission * throughput * medium->get_radiance(mei, active_medium);
dr::masked(throughput, act_null_scatter) *= mei.sigma_n * weight_null;
dr::masked(result, act_null_scatter && (dr::max(radiance) > 0.f)) += weight_emission * throughput * radiance;
dr::masked(throughput, is_spectral && act_null_scatter) *= weight_null * mei.sigma_n;

// Move the ray along
dr::masked(ray.o, act_null_scatter) = mei.p;
Expand All @@ -322,7 +338,7 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator<Float, Spectrum> {
act_medium_scatter &= active;

if (dr::any_or<true>(act_medium_scatter)) {
masked(throughput, act_medium_scatter) *= mei.sigma_s * weight_scatter;
masked(throughput, act_medium_scatter) *= weight_scatter * mei.sigma_s;

PhaseFunctionContext phase_ctx(sampler);
auto phase = mei.medium->phase_function();
Expand Down Expand Up @@ -367,7 +383,8 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator<Float, Spectrum> {
Mask count_direct = ray_from_camera || specular_chain;
EmitterPtr emitter = si.emitter(scene);
Mask active_e = active_surface && dr::neq(emitter, nullptr)
&& !(dr::eq(depth, 0u) && m_hide_emitters);
&& !(dr::eq(depth, 0u) && m_hide_emitters)
&& !has_flag(emitter->flags(), EmitterFlags::Medium); // Ignore any emitting media as they are not compatible with MIS yet
if (dr::any_or<true>(active_e)) {
Float emitter_pdf = 1.0f;
if (dr::any_or<true>(active_e && !count_direct)) {
Expand Down Expand Up @@ -572,7 +589,7 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator<Float, Spectrum> {

MI_DECLARE_CLASS()
protected:
VolumeEventProbabilityMode m_sampling_type;
VolumeEventProbabilityMode m_sampling_mode;
};

MI_IMPLEMENT_CLASS_VARIANT(VolumetricPathIntegrator, MonteCarloIntegrator);
Expand Down
32 changes: 17 additions & 15 deletions src/integrators/volpathmis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -227,7 +227,14 @@ class VolpathMisIntegratorImpl final : public MonteCarloIntegrator<Float, Spectr
}

if (dr::any_or<true>(active_medium)) {
Mask null_scatter = sampler->next_1d(active_medium) >= index_spectrum(mei.sigma_t, channel) / index_spectrum(mei.combined_extinction, channel);
auto radiance = medium->get_radiance(mei, active_medium);
auto emission_prob = radiance * dr::mean(0.5f + 0.5f * mis_weight(p_over_f)),
scatter_prob = mei.sigma_t,
null_prob = mei.sigma_n;
dr::masked(emission_prob, dr::max(radiance) > 0.f) = 1.0f;
auto prob_sum = emission_prob + scatter_prob + null_prob;

Mask null_scatter = sampler->next_1d(active_medium) < (index_spectrum(emission_prob, channel) + index_spectrum(null_prob, channel)) / index_spectrum(prob_sum, channel);
act_null_scatter |= null_scatter && active_medium;
act_medium_scatter |= !act_null_scatter && active_medium;
last_event_was_null = act_null_scatter;
Expand All @@ -243,31 +250,25 @@ class VolpathMisIntegratorImpl final : public MonteCarloIntegrator<Float, Spectr
specular_chain |= act_medium_scatter && !sample_emitters;

if (dr::any_or<true>(act_null_scatter)) {
if (dr::any_or<true>(is_spectral)) {
update_weights(p_over_f, mei.sigma_n / mei.combined_extinction, mei.sigma_n, channel, is_spectral && act_null_scatter);
update_weights(p_over_f_nee, 1.0f, mei.sigma_n, channel, is_spectral && act_null_scatter);
}
if (dr::any_or<true>(not_spectral)) {
update_weights(p_over_f, mei.sigma_n, mei.sigma_n, channel, not_spectral && act_null_scatter);
update_weights(p_over_f_nee, 1.0f, mei.sigma_n / mei.combined_extinction, channel, not_spectral && act_null_scatter);
}
auto p_over_f_old = p_over_f;
update_weights(p_over_f_old, emission_prob / prob_sum, 1.0f, channel, act_null_scatter);
dr::masked(result, act_null_scatter) += mis_weight(p_over_f_old) * radiance;
update_weights(p_over_f, null_prob / prob_sum, mei.sigma_n, channel, act_null_scatter);
update_weights(p_over_f_nee, 1.0f, mei.sigma_n, channel, act_null_scatter);

dr::masked(ray.o, act_null_scatter) = mei.p;
dr::masked(si.t, act_null_scatter) = si.t - mei.t;
}

if (dr::any_or<true>(act_medium_scatter)) {
if (dr::any_or<true>(is_spectral))
update_weights(p_over_f, mei.sigma_t / mei.combined_extinction, mei.sigma_s, channel, is_spectral && act_medium_scatter);
if (dr::any_or<true>(not_spectral))
update_weights(p_over_f, mei.sigma_t, mei.sigma_s, channel, not_spectral && act_medium_scatter);
update_weights(p_over_f, scatter_prob / prob_sum, mei.sigma_s, channel, act_medium_scatter);

PhaseFunctionContext phase_ctx(sampler);
auto phase = mei.medium->phase_function();

// --------------------- Emitter sampling ---------------------
valid_ray |= act_medium_scatter;
Mask active_e = act_medium_scatter && sample_emitters;
Mask active_e = act_medium_scatter && sample_emitters; // TODO: Enable MIS with volume emitters
if (dr::any_or<true>(active_e)) {
auto [p_over_f_nee_end, p_over_f_end, emitted, ds] =
sample_emitter(mei, scene, sampler, medium, p_over_f,
Expand Down Expand Up @@ -309,7 +310,8 @@ class VolpathMisIntegratorImpl final : public MonteCarloIntegrator<Float, Spectr
Mask ray_from_camera = active_surface && dr::eq(depth, 0u);
Mask count_direct = ray_from_camera || specular_chain;
EmitterPtr emitter = si.emitter(scene);
Mask active_e = active_surface && dr::neq(emitter, nullptr) && !(dr::eq(depth, 0u) && m_hide_emitters);
Mask active_e = active_surface && dr::neq(emitter, nullptr) && !(dr::eq(depth, 0u) && m_hide_emitters)
&& !has_flag(emitter->flags(), EmitterFlags::Medium);
if (dr::any_or<true>(active_e)) {
if (dr::any_or<true>(active_e && !count_direct)) {
// Get the PDF of sampling this emitter using next event estimation
Expand Down
Loading

0 comments on commit e3a5a8c

Please sign in to comment.