mitsuba-renderer
[✨ feature] Emissive Volume Rendering
Description
This PR adds emissive event sampling in the null-scattering volpathmis integrator for sampling the VRE and MIS for emission. It supersedes the previous PR https://github.com/mitsuba-renderer/mitsuba3/pull/40.
Adds emissive volume support by adding a volume emitter. The interface adds a volumelight emitter to a given shape that is then sampled as part of the medium sampling routines in volpath/volpathmis/prbvolpath.
The current implementation modifies none of the interfaces of endpoint and thus does not work with MIS as-of-yet. I am unsure of the best way to implement this in terms of modifying the interface. The cleanest solution seems to be that the endpoint sampling methods take an extra float to sample a 3d point, but whether this should be implemented by changing the existing Point2f sample to Point3f or adding an additional argument is unclear to me.
The current implementation is feature complete except for NEE for volume emitters, I have added and tested the python interfaces by bringing the prbvolpath implementation to parity with volpath.
Checklist
- [x] My code follows the style guidelines of this project
- [x] My changes generate no new warnings ~~- generates warnings about unused parameters in
volumelight, these methods will be implemented and the parameters used.~~ There are no longer any new warnings. - [x] My code also compiles for
cuda_*andllvm_*variants. I have compiled every single variant and run all the tests, and everything passes except for the tutorial notebooks, but these throw an unicode char encoding error so I think it's unrelated. I will check with a fresh master copy if the same error occurs, might be a windows only issue. - [x] I have commented my code
- [x] I have made corresponding changes to the documentation
- [x] I have added tests that prove my fix is effective or that my feature works - ~~this is in-progress, have added emissive volume scenes and tests that I will make a PR for in the mitsuba-renderer/mitsuba-data repo.~~ Created the following PR to merge emissive volume tests: mitsuba-renderer/mitsuba-data#17
- [x] I cleaned the commit history and removed any "Merge" commits
- [x] I give permission that the Mitsuba 3 project may redistribute my contributions under the terms of its license
Notes
- The implementation uses non-analogue probabilities as described in Monte Carlo Methods for Volumetric Light Transport Simulation. This leads to more consistent results when a given radiance channel is zero and better convergence in all the tests I've done.
- Both
volpathandprbvolpathhave implementations of the analogue, maximum and mean event probability calculations for testing. The issue of inconsistency only appears with analogue probabilities in non-scalar modes, and only for volpathmis and prbvolpath. - Sorry for the delay in the implementation! The move to Mitsuba 3 meant I had to rewrite a lot of the code, but I managed to be significantly more thorough this time around so everything compiles and all the tests pass.
- I've attached an example render as well!
I've identified the issue with the failing tutorial notebook, it's specifically related to the prbvolpath integrator and my incorrect implementation of accumulating radiance. I'm not quite sure what the correct implementation is as my understanding of it is that the radiance accumulated from the medium should be treated similarly to intersection with emitters, but this seems to break the differentiability of the rest of the medium parameters.
Hi @Microno95
I quickly skimmed through this PR. It's looking good! I'll need to spend more time to carefully review it, and think about some of the apparently necessary interface change :smile:
No worries! I've made a list of the necessary features for NEE sampling volumes, hopefully it helps with deciding on the interface and I'm very happy to discuss it in more detail as well 😊
Extensions to shape:
- Sample a 3d point inside a shape
- Check that a point is inside a shape
- Calculate (or estimate) the volume of a shape
I think these should belong to the shape class, same as sample_position, surface_area and pdf_position, but generalising to the enclosing volume. The only issue is how to handle shapes that are not watertight, but I guess that also has issues with other aspects so the onus falls on the user.
And then the emitter calls to the appropriate method while forwarding the extra random sample as needed.
Nice work! Sorry I somehow missed this new PR.
Just a high level comment for now: This adds quite a lot of code complexity - would it make sense to have special integrators for emissive media? It's somewhat of a special case and the existing volpath integrators are already hard to parse (e.g., for someone coming from using NeRF...)
Also, would it make sense to move some of those "medium probabilities" functions to the medium class?
I have significantly simplified the sampling of emission and reduced the complexity of the integrators in my latest code that I will push soon.
Generally, I think it does make sense to move these to the medium class as the different medium event sampling modes are better suited for different media. e.g. if a medium has 0 scattering for, say the green channel, the analogue probabilities actually underestimate the number of scatters when using hero wavelength sampling and lead to significant fireflies.
Turns out you can sample radiance at every medium interaction (at least when not using NEE for medium emitters, that's a separate case to be worked on), and this has simplified the inclusion of radiative media effects to a single additional line in both volpath and volpathmis integrators and only 4 additional lines in prbvolpath (a few more for the derivative tracking).
I have held off on pushing any updates as I was waiting on how I should go about updating the endpoint/emitter interface to support a third random sample.
I see, thanks for clarifying.
I don't think adding an extra sample to the endpoint interface is a big deal.
How do you deal with sampling the volume of complex shapes? Do you simply sample the bounding box and fail the sample if it falls outside the shape? I think currently it's still tricky to use an arbitrary number of random samples within a virtual function call (e.g., to do rejection sampling)
Does it make sense to add the 3rd volume sample by changing the Point2f sample to a Point3f sample? It seems the least cluttered way and is the most logical extension. But happy to also make it an additional parameter in the interface.
I have not yet decided on the best way to sample triangle mesh-based shapes. My two ideas were to either sample within the bounding box and fail the sample if it falls outside as you said, or to build a discrete distrubution using a simple tetrahedral decomposition of the mesh.
The upside of using a tetrahedral representation (or really any discretisation of the volume) is that it'll be faster to sample and easier to check if a sample falls inside/outside the shape.
The alternative is iterating over every triangle to calculate a winding number or using ray intersections (calculating the number of crossings). The issue with ray intersections is in regards to choosing a direction for testing and I've found that it can be problematic depending on the mesh and winding number calculation scales linearly with the number of triangles.
I think changing the sample to be Point3f isn't a big deal, better than adding another parameter.
I don't think adding a tetrahedral representation makes sense (from an engineering point of view). As far as I know - correct me on that - computing a tetrahedralization is far from trivial and not something that can easily be implemented in. a standalone helper function.
I would suggest to check if a point is inside the mesh by simply tracing a ray in a fixed (or randomized) direction and checking if the triangle we hit is facing away from the ray direction. This should work for watertight meshes, and that either way is the assumption for any of the volumetric integrators to work correctly.
Hey @dvicini, I am in the process of implementing the interface for Next Event Estimation of volume emitters, but I'm struggling with estimating the pdf of a direction sampled in the emitter, specifically implementing the sample_direction and pdf_direction methods.
The issue I'm having is that while I know the pdf of sampling a given point in the medium, the pdf of a direction is with respect to the solid angle whereas the pdf of sampling a point is with respect to the volume. I'm not sure how to make these two commensurate.
I know the following relation should hold
$$ \int_\Omega p(\vec{\omega}) \mathrm{d\Omega} = \int_V p(x) \mathrm{dV} $$
i.e. the integral of the solid angle pdf over all solid angles should be equal to the integral of the volume pdf over all volumes (and both should be equal to 1 for a normalised pdf.). Knowing that, I can rewrite the equation as follows.
$$ \int_\Omega p(\vec{\omega}) \mathrm{d\Omega} = \int_\Omega\int_r p(x) r^2\mathrm{dr}\mathrm{d\Omega} $$
And (this is the part I'm unsure about), this should imply that
$$ p(\vec{\omega}) \mathrm{d\Omega} =\int_r p(x) r^2\mathrm{dr}\mathrm{d\Omega} $$
If all that is correct, then the pdf of sampling a given direction should be straightforward to compute based on the line segments that pass through the medium - easily found via ray intersections. I'm not sure how this is to be treated in the case of medium vs bsdf interactions given that $\mathrm{d\Omega}$ has two different domains, how to handle the case where a given point is inside the medium and how to handle refractive interfaces (do they matter?).
The Jacobian to change between solid angle and volume sampling is indeed $r^2$. This is exactly the geometry term in a volume.
I am not sure I fully understand you question. If this is primarily about MIS, you will just want to make sure that both PDFs are in the same measure.
The conceptually simplest way is to integrate the volume PDF over the ray segment, but that is costly and potentially inaccurate (e.g. can be done using ray marching). It's good to keep in mind that MIS weights don't have to use exact PDFs, they just need to be computed consistently to ensure unbiasedness.
Another option could be to compute the MIS in volume measure, I think that would forego the need to integrate over segments (potentially at an increase in noise? but with cheaper per-sample cost). In that case, you would just have to convert the solid angle PDF to volume PDF using the geometry term.
Right, I think I understand better now. I was thinking I'd need to convert the volume PDF to a solid angle PDF in order to do MIS with the volume emitter, but you're saying I can do it the other way around (which makes perfect sense).
In terms of converting the solid angle PDF to a volume PDF, would it be correct to divide $p(\vec{\omega})$ by $r^2$ or does it require something else? i.e. the MIS is computed via $\mathrm{MIS}(r^{-2}p(\vec\omega), p(\vec x))$
Yes, this sounds good to me.
You will need to make sure that the MIS is consistent with the evaluation of emission after direction sampling. I.e., the simplest method would again be to evaluate the emission along the entire ray segment.
It seems that simply converting the solid angle PDF to a volume PDF does not work, or perhaps I'm doing something wrong. I've verified that integrating the volume PDF over the ray length works, but it is slower (by a lot, around 20x-40x due to the ray casts required) and raises OptiX warnings in CUDA modes.
At first I thought it would be a simple change of variables using the differential element $$r^{-2}\frac{\mathrm{dV}}{\mathrm{d\Omega}}$$, but this gives very different results compared to integrating over the ray length and is inconsistent with the previous results obtained without NEE.
Is it perhaps to do with the fact that $r^{-2}p(\vec{\omega})dV$ is a non-normalisable PDF? i.e. the following integral does not converge
$$ \int_V r^{-2} p(\vec{\omega}) \mathrm{dV} = \int_V \left(r^{-2}p(\vec{\omega})\right) \mathrm{r^2\sin\theta dr d\theta d\phi} = \int_V p(\vec{\omega})\mathrm{\sin\theta dr d\theta d\phi} $$
The current implementation correctly tracks the ray entry and exit points in a given shape and calculates the exact directional pdf. In llvm_* modes this is quite slow compared to not using NEE by around 5x while in cuda_* modes the discrepancy is much less at around 1.5x slower with a noticeable reduction in the noise for the same sample count.
# ~~~~ Was wrong about this ~~~~ #
There is one big issue in that in llvm_* modes, the ray entry and exit points lead to different results compared to cuda_* modes, most likely due to the different handling of ray intersections by Embree vs OptiX and thus the current implementation only gives correct results for cuda_* modes.
I am working on improving and fixing the implementation, but I would welcome a second set of eyes over the code to see if there's anything I'm missing in how I'm accumulating the pdf. # ~~~~ ~~~~~~~~ ~~~~ #
Okay, nevermind, it seems that both are incorrect for an as-of-yet unknown reason by the same amount. With OptiX I was getting the correct results in that case due to a bug so I'll investigate down that path.
Looking further into, it seems that cuda_* modes gave the correct results when the ray intersection routines were disabled giving a pdf of 0 - this was due to how I initialized the nested volume scene allowing me to do ray intersections. Ultimately this gave the correct result as it effectively disabled NEE.
I've attached two renders, with and without NEE (512x512@4096spp), that show the discrepancy between the NEE and non-NEE results. This is of a sphere shape with a purely absorbing medium with a radiance of 1 and an extinction coefficient of 1 at the origin with a radius of 1.
The discrepancy definitely has to do with how the PDF is computed, but I'm not sure what's wrong as I've verified that the calculations are correct - by computing the analytic solution to the directional PDF and comparing it to the ray traced solution for the sphere shape as per our earlier discussions. From the emission gathered at the floor and back wall, it looks like the NEE version is underestimating the illumination from the medium - perhaps there is a term that has to be accounted for.
Without NEE:
With NEE:
I've also attached the .exr files themselves for better numerical comparison in TEV. emitting_sphere_renders.zip
Hmm, really hard to say what's wrong just from these two images. I would say that for a homogeneous sphere of emissive volume it should really be possible to get a perfect match, as it doesn't have any of the complexity of complex shapes or delta tracking for heterogeneous media. My only suggestion would be to further reduce the complexity: 1) only render direct illumination, no scattering and 2) completely disable absorption. I.e., can you match the total emission from an emissive sphere volume that does not absorb any light?
With a homogeneous sphere where $\sigma_t=\sigma_s=0$ and $\sigma_n=1$ (in order to gather emission from the medium), the renders still diverge for any ray max depth $>1$. For a max ray depth of 1, so just purely illumination from emitter intersections, the two renders agree perfectly. For a max ray depth of 2 or greater, they diverge with the NEE renders looking similarly dimmer as before.
In the literature, the closest I've seen to the implementation of the solid angle PDF is in Equation 15 of Line Integration for Rendering Heterogeneous Emissive Volumes where they estimate the solid angle PDF by integrating the volume PDF along the ray direction with the appropriate jacobian.
Perhaps there is additional consideration in how NEE estimates of illumination are combined with Unidirectional estimates for volume emitters.
I should add that a purely emissive medium is not possible with the current homogeneous medium handling in volpath/volpathmis due to the way distance sampling in the medium is computed and how homogeneous media are optimised in the code paths. i.e. to gather emission from homogeneous emissive media with no absorption, there is no way to sample free-flight distances and gather samples at different points along the ray since the ray immediately exits the medium. Thus I achieve the renders by using a heterogeneous medium with $\sigma_{\left[t,s\right]}=0$ and $\sigma_n=1$.
I see, yes it makes sense that for this setting you need to use some non-zero null density to even sample the emission contributions. I would investigate this case further, this seems like a good minimal test case to debug.
I've looked further into it, and I think there's actually a bug in how volpath handles NEE. I've created a case where a surface emitter with a diffuse reflectance of 0 gives incorrect results with volpath when emitter sampling is enabled, regardless of whether or not there are media in the scene at all.
Here's the reference render using the Path integrator:
And here's the render from VolPath:
I've also included the renders from volpathmis and prbvolpath as well.
volpathmis:
prbvolpath:
I've tested it and the issue exists in the master branch so I am sure it is unrelated to my changes to the code. I've attached the scene here:
minimal_example.zip
Oh, that's a good catch. Not sure what's wrong here, but maybe something related to the ShadowEpsilon or so. I will take a look, would be good to fix
Had a chance to take a look while wrapping my head around the sample_emitter implementation, and as a byproduct, figured out the issue. Essentially, the sample_emitter generates a path with a null vertex on the surface of the emitter before returning; in the case of a surface emitter, this is easily alleviated by avoiding the generation of a null vertex whenever the remaining distance is below the ShadowEpsilon threshold after a surface intersection before https://github.com/mitsuba-renderer/mitsuba3/blob/ffcde97708f3310727eb04422608972c2628e503/src/integrators/volpath.cpp#L419-L425
remaining_dist = ds.dist * (1.f - math::ShadowEpsilon<Float>) - total_dist;
ray.maxt = remaining_dist;
active &= (remaining_dist > math::ShadowEpsilon<Float>);
if (dr::none_or<false>(active))
break;
I think the problem is a bit more subtle. By adding the epsilon to the check for the remaining distance, you are essentially doubling the shadow epsilon. The epsilon is already included in the computation of remaining_dist. The better solution is to replace the spawn_ray function call with a call to spawn_ray_to, which has a better epsilon calculation. Then, we can also replace the use of ds.dist with the maxt set within spawn_ray_to. This is more consistent with what the path integrator does. I will make a PR
Ah that makes sense, I will leave that to you and rebase on top when it's in master.
Regarding NEE with the volume emitter, I am thinking that perhaps the issue lies with how the path to the sampled emission point is constructed. Normally with surface emitters a random set of null vertices are constructed from the interaction point to the emitter, but this is all with the assumption that the emitter interaction occurs on a surface so null paths through media and interfaces up to the last vertex, $p_n$, are all quite sensible.
When a medium emitter is sampled, $p_n$ is inside the medium and so some of the assumptions of the sample_emitter method are broken. For example, mei.t > remaining_dist doesn't necessarily imply that the ray escaped the medium as assumed in https://github.com/mitsuba-renderer/mitsuba3/blob/ffcde97708f3310727eb04422608972c2628e503/src/integrators/volpath.cpp#L391
My idea is that up until the second-to-last vertex, $p_{n-1}$, we construct null-paths and then the connection from $p_{n-1}$ to $p_n$ is constructed deterministically whenever the sampled free flight distance exceeds the distance to $p_n$. The problem is that I am not sure what the PDF of such a path would be.
EDIT: I no longer think that makes sense as a solution. I think the problem may require formulating the constructed NEE path construction in the integral formulation a-la https://cs.dartmouth.edu/wjarosz/publications/miller19null.html and then finding the pdf of a given constructed path. This is hefty 😅
In order to simplify some of testing, I separated out emitter sampling and unidirectional sampling in volpath. Separately, the two methods give near identical results for an analytic shape like a sphere even in the presence of absorption and volumetric scattering.
Here are two (low-res) renders to show the results:
| Unidirectional Sampling | Emitter Sampling |
|---|---|
The convergence for emitter sampling is atrocious, it requires an order of magnitude more samples for this case than unidirectional sampling - most likely due to the fact that I can control how densely unidirectional sampling accumulates emission via $\sigma_n$, but emitter sampling is indifferent to it. Either way, I can at least show the implementation of the sampling procedure is correct.
There are now two issues to resolve:
- Rejection-like sampling of the mesh isn't working and I'm not exactly sure why. It seems that setting the PDF to $0$ when outside the shape and $V^{-1}$ when inside gives incorrect results - it significantly underestimates the emission compared to unidirectional sampling.
- MIS between the two sampling procedures is problematic as the PDF of a given direction is not compatible with the fact that I'm accumulating radiance from a single point rather than the emission along the ray segment passing through the medium in that direction. For this, I think a possible solution is to use the same direction sampling procedure, but instead of accumulating the radiance from a point, accumulate the radiance along the ray in
sample_emitter. I will try that and see how it goes...
I apologise for turning this pull request comment thread into a stream of thought type log, but it has helped me work through the code and I'm hoping that it provides some insight into the implementation when it's finalised.
I have a version of the code that works, it uses emitter sampling to sample the emitting sphere media and estimates the radiance along the ray through the medium, and can MIS the results with unidirectional sampling. This only works in volpath at the moment, but it works:
| Unidirectional Sampling | Emitter Sampling | Multiple Importance Sampling |
|---|---|---|
Nice, looks good! FYI, the PR fixing the epsilon issues has been merged in the meantime.
Hey @dvicini,
So I've got a solution implemented for NEE for both the volpath and prbvolpath integrators and they agree with each other in terms of results and with the no-NEE integrators. I'm still working on the volpathmis solution, but how emitter sampling is done is somewhat more complicated as the radiance along the ray passing through the medium is estimated by the medium stepping routine in the emitter sampling code.
I'll hopefully have that sorted in the next day or two. What I wanted to ask about is how to correctly implement the adjoint for the Path Replay Backpropagation version of the volpath integrator. I've an implementation in the repository, but I'm not sure it is entirely correct and I'm not sure how best to check it - would using the inverse_rendering/volume_optimization.ipynb notebook be a good start by trying to do the optimisation with an emissive medium instead of an absorbing/scattering one or is there something more basic that occurs to you?
EDIT: Another question I had is in regards to delta/ratio tracking. Which ray segment does the transmittance estimator apply to? If there are, say 3 vertices along the ray (x0, x1 and x2) which are all medium interactions, does the estimator computed at x1 apply to the segment x0 to x1 or x1 to x2? I think it's the latter, but I wanted to double-check.
Hi,
That sounds great! Yes, I think the volume_optimization notebook should be a good starting point to work with the PRB version of the integrator. You should also take a look at the test_ad_integrators.py file to see how we run finite-difference validation for our differentiable integrators. For the existing PRB integrators, we implemented both a forward and reverse-mode version. The forward mode version can be used to output "gradient images" which allow to more easily visually compare to a finite difference gradient (for 1 single parameter in this case).
I am not sure i understand the question about delta/ratio tracking. We use delta tracking throughout the medium and ratio tracking to evaluate the transmittance term during NEE.
Regarding my question, it was simply about which ray segment the computed transmittance estimate applies to. More specifically if an interaction happens at vertex x1, is the transmittance relevant to the emission from x1 the estimated transmittance before the null-scattering update (i.e. the $1 - \mu/\bar\mu$ factor) or after (turns out it's before, it makes sense but I did double check using Blender).
I've got the gradients in the PRBVolpath integrator sorta working, when loop recording is disabled the test with finite differencing passes, but when it is enabled, python crashes with Fatal Python Error: Aborted. And it also crashes in Debug builds regardless of whether Loop Recording is enabled or not. Bit baffling, might not be to do with my changes though as I did rebase the repo on top of master recently.
Furthermore, I added multiple AD integrator tests for the volumetric path tracing integrators and I'm finding that every volpath integrator passes all of the tests in forward-mode but not the backward mode when the target variable is $\sigma_t$ or the albedo and the emitter is a constant emitter. It is quite strange...
Also, is forward-mode differentation meant to be enabled for PRBVolpath? It throws an error when run in forward-mode, but your message seemed to imply that it should also work in forward-mode.