Reading list on ReSTIR

Recently a short video from dark magic programmer Tomasz Stachowiak made the rounds in the graphics programming community, at the sound of jaws hitting the floor in its wake. It shows his recent progress on in his renderer pet project: beautiful real-time global illumination with fast convergence and barely any noise, in a static environment with dynamic lighting.

In a Twitter thread where he discussed some details, one keyword in particular caught my attention: ReSTIR.

ReSTIR stands for “Reservoir-based Spatio-Temporal Importance Resampling” and is a sampling technique published at SIGGRAPH 2020 and getting refined since.

The original publication

Spatiotemporal reservoir resampling for real-time ray tracing with dynamic direct lighting
The publication page includes the recording of the SIGGRAPH presentation, with a well articulated explanation of the technique by main author Benedikt Bitterli.
(same publication hosted on the NVidia website).

Explanations of ReSTIR

Improvements over the original publication

After the initial publication, NVidia published a refined version producing images with less noise at a lower cost, which they call “RTXDI” (for RTX Direct Illumination).

Other limitations

When discussing on Twitter some of the limitations of ReSTIR, Chris Wyman made the following remarks:

To be clear, right now, ReSTIR is a box of razor blades without handles (or a box of unlabeled knobs). It’s extremely powerful, but you have to know what you’re doing. It is not intuitive, if your existing perspective is traditional Monte Carlo (or real-time) sampling techniques.

People sometimes think SIGGRAPH paper = solved. Nope. We’ve learned a lot since the first paper, and our direct lighting is a lot more stable with that knowledge. We’re still learning how to do it well on full-length paths.

And there’s a bunch of edge cases, even in direct lighting, that we know how to solve but haven’t had time to write them up, polish, and demo.

We haven’t actually tried to solve the extra noise at disocclusions in (what I think of as) a very principled way. Right now a world-space structure is probably the best way. I’m pretty sure it can be done without a (formal) world-space structure, just “more ReSTIR.”

Intersection of a ray and a plane

I previously showed the derivation of how to determine the intersection of a plane and a cone. At the time I had to solve that equation, so after doing so I decided to publish it for anyone to use. Given the positive feedback, it seems this was useful, so I might as well continue with a few more.

Here is probably the most basic intersection: a ray and a plane. Solving it is straightforward, which I hope can be seen below. Like last time, I am using vector notation.

  1. We define a ray with its origin $O$ and its direction as a unit vector $\hat{D}$.
    Any point $X$ on the ray at a signed distance $t$ from the origin of the ray verifies: $\vec{X} = \vec{O} + t\vec{D}$.
    When $t$ is positive $X$ is in the direction of the ray, and when $t$ is negative $X$ is in the opposite direction.
  2. We define a plane with a point $S$ on that plane and the normal unit vector $\hat{N}$, perpendicular to the plane.
    The distance between any point $X$ and the plane is $d = \lvert (\vec{X} – \vec{S}) \cdot \vec{N} \rvert$. If this equality is not obvious for you, you can think of it as the distance between $X$ and $S$ along the $\vec{N}$ direction. When $d=0$, it means $X$ is on the plane itself.
  3. We define $P$ the intersection or the ray and the plane, and which we are interested in finding.

Since $P$ is both on the ray and on the plane, we can write: $$ \left\{ \begin{array}{l} \vec{P}=\vec{O} + t\vec{D} \\ \lvert (\vec{P} – \vec{S}) \cdot \vec{N} \rvert = 0 \end{array} \right. $$ Because the distance $d$ from the plane is $0$, the absolute value is irrelevant here. We can just write: $$ \left\{ \begin{array}{l} \vec{P}=\vec{O} + t\vec{D} \\
(\vec{P} – \vec{S}) \cdot \vec{N} = 0 \end{array} \right. $$ All we have to do is replace $P$ with $\vec{O} + t\vec{D}$ in the second equation, and reorder the terms to get $t$ on one side.
$$ (\vec{O} + t\vec{D} – \vec{S}) \cdot \vec{N} = 0 $$ $$ \vec{O} \cdot \vec{N} + t\vec{D} \cdot \vec{N} – \vec{S} \cdot \vec{N} = 0 $$ $$ t\vec{D} \cdot \vec{N} = \vec{S} \cdot \vec{N} – \vec{O} \cdot \vec{N} $$ $$ t = \frac{(\vec{S} – \vec{O}) \cdot \vec{N}}{ \vec{D} \cdot \vec{N} } $$

A question to ask ourselves is: what about the division by $0$? Looking at the diagram, we can see that $\vec{D} \cdot \vec{N} = 0$ means the ray is parallel to the plane, and there is no solution unless $O$ is already on the plane. Otherwise, the ray intersects the plane for the value of $t$ written above. That’s it, we’re done.

Note: There are several, equivalent, ways of representing a plane. If your plane is not defined by a point $S$ and a normal vector $\hat{N}$, but rather with a distance to the origin $s$ and a normal vector $\hat{N}$, you can notice that $s = \vec{S} \cdot \vec{N}$ and simplify the result above, which becomes: $$ t = \frac{s – \vec{O} \cdot \vec{N}}{ \vec{D} \cdot \vec{N} } $$


Signed distance to a plane

For the sake of simplicity, in the above we defined the distance to the plane as an absolute value. It is possible however to define it as a signed value: $d = (\vec{X} – \vec{S}) \cdot \vec{N}$. In this case $d>0$ means $X$ is somewhere on the side of the plane pointed by $\vec{N}$, while $d<0$ means $X$ is on the opposite side of the plane.

Distances that can be negative are called signed distances, and they are a foundation of Signed Distance Fields (SDF).

Long hiatus

Last week I was lucky enough to attend SIGGRAPH 2018, in Vancouver. My colleagues and I were presenting on a booth the work we had done, a VR story with a distinctive comic book look. I was also invited to participate to a panel session on demoscene, where I shared some lessons learned while making the 64k intro H – Immersion. The event brought a certain sense of conclusion to this work, aside from filling me with inspiration and motivation to try new things.

It has been a long time since I last posted anything here. For the last two years the majority of my spare time went into making that 64k intro. In fact the last post, “Intersection of a ray and a cone”, was related to it. I was implementing volumetric lighting for the underwater scenes, and wanted to resolve cones of light with ray tracing, before marching inside those cones. LLB and I have talked about the creation process in two making-of articles: “A dive into the making of Immersion”, and “Texturing in a 64kB intro”.

During that time, a lot of new things have happened in the computer graphics community. It has been difficult to keep track of everything. The last topic I started experimenting with is point cloud and mesh capture from photos; I might expend on it here in the future. I also want to experiment with DIY motion capture. Anyway, it’s time to resume posting here.

Intersection of a ray and a cone

Some time ago I needed to solve analytically the intersection of a ray and a cone. I was surprised to see that there are not that many resources available; there are some, but not nearly as many as on the intersection of a ray and a sphere for example. Add to it that they all use their own notation and that I lack math exercise, after a bit of browsing I decided I needed to write a proof by myself to get a good grasp of the result.

So here goes, the solution to the intersection of a ray and a cone, in vector notation.

  1. We define a ray with its origin $O$ and its direction as a unit vector $\hat{D}$.
    Any point $X$ on the ray at a signed distance $t$ from the origin of the ray verifies: $\vec{X} = \vec{O} + t\vec{D}$.
    When $t$ is positive $X$ is in the direction of the ray, and when $t$ is negative $X$ is in the opposite direction.
  2. We define a cone with its tip $C$, its axis as a unit vector $\hat{V}$ in the direction of increasing radius, and $\theta$ the half angle between the axis and the surface.
    Any point $X$ on the cone verifies: $(\vec{X} – \vec{C}) \cdot \vec{V} = \lVert \vec{X} – \vec{C} \rVert \cos\theta$
  3. Finally we define $P$ the intersection or the ray and the cone, and which we are interested in finding.

$P$ verifies both equations, so we can write:

$$
\left\{
\begin{array}{l}
\vec{P}=\vec{O} + t\vec{D} \\
\frac{ \vec{P} – \vec{C} }{\lVert \vec{P} – \vec{C} \rVert} \vec{V} = \cos\theta
\end{array}
\right.
$$

We can multiply the second equation by itself to work with it, then reorder things a bit.

$$
\left\{
\begin{array}{l}
\vec{P}=\vec{O} + t\vec{D} \\
\frac{ ((\vec{P} – \vec{C}) \cdot \vec{V})^2 }{ (\vec{P} – \vec{C}) \cdot (\vec{P} – \vec{C}) } = \cos^2\theta
\end{array}
\right.
$$

$$
\left\{
\begin{array}{l}
\vec{P}=\vec{O} + t\vec{D} \\
((\vec{P} – \vec{C}) \cdot \vec{V})^2 – (\vec{P} – \vec{C}) \cdot (\vec{P} – \vec{C}) \cos^2\theta = 0
\end{array}
\right.
$$

Remember the mouthful earlier about $\hat{V}$ being in the direction of increasing radius? By elevating $\cos\theta$ to square, we’re making negative values of $\cos$ positive: values of $\theta$ beyond 90° become indistinguishable from values below 90°. This has the side effect of turning it into the equation of not one, but two cones sharing the same axis, tip and angle, but in opposite directions. We’ll fix that later.

We replace $\vec{P}$ with $\vec{O} + t\vec{D}$ and work the equation until we get a good old quadratic function that we can solve.

$$
\require{cancel}
((\vec{O} + t\vec{D} – \vec{C})\cdot\vec{V})^2 – (\vec{O} + t\vec{D} – \vec{C}) \cdot (\vec{O} + t\vec{D} – \vec{C}) \cos^2\theta = 0
$$

$$
((t\vec{D} + \vec{CO})\cdot\vec{V})^2 – (t\vec{D} + \vec{CO}) \cdot (t\vec{D} + \vec{CO}) \cos^2\theta = 0
$$

$$
(t\vec{D}\cdot\vec{V} + \vec{CO}\cdot\vec{V})^2 – (t^2\cancel{\vec{D}\cdot\vec{D}} + 2t\vec{D}\cdot\vec{CO} + \vec{CO}\cdot\vec{CO}) \cos^2\theta = 0
$$

$$
(t^2(\vec{D}\cdot\vec{V})^2 + 2t(\vec{D}\cdot\vec{V})(\vec{CO}\cdot\vec{V}) + (\vec{CO}\cdot\vec{V})^2) – (t^2 + 2t\vec{D}\cdot\vec{CO} + \vec{CO}\cdot\vec{CO}) \cos^2\theta = 0
$$

$$
t^2(\vec{D}\cdot\vec{V})^2
+ 2t(\vec{D}\cdot\vec{V})(\vec{CO}\cdot\vec{V})
+ (\vec{CO}\cdot\vec{V})^2
– t^2\cos^2\theta
– 2t\vec{D}\cdot\vec{CO}\cos^2\theta
– \vec{CO}\cdot\vec{CO}\cos^2\theta
= 0
$$

Reorder a bit:

$$
t^2((\vec{D}\cdot\vec{V})^2 – \cos^2\theta)
+ 2t((\vec{D}\cdot\vec{V})(\vec{CO}\cdot\vec{V}) – \vec{D}\cdot\vec{CO}\cos^2\theta)
+ (\vec{CO}\cdot\vec{V})^2 – \vec{CO}\cdot\vec{CO}\cos^2\theta
= 0
$$

There we go, we have our $at^2 + bt + c = 0$ equation, with:

$$
\left\{
\begin{array}{l}
a = (\vec{D}\cdot\vec{V})^2 – \cos^2\theta \\
b = 2\Big((\vec{D}\cdot\vec{V})(\vec{CO}\cdot\vec{V}) – \vec{D}\cdot\vec{CO}\cos^2\theta\Big) \\
c = (\vec{CO}\cdot\vec{V})^2 – \vec{CO}\cdot\vec{CO}\cos^2\theta
\end{array}
\right.
$$

From there, you know the drill: calculate the determinant $\Delta = b^2 – 4ac$ then depending on its value:

  • If $\Delta < 0$, the ray is not intersecting the cone.
  • If $\Delta = 0$, the ray is intersecting the cone once at $t = \frac{-b}{2a}$.
  • If $\Delta > 0$, the ray is intersecting the cone twice, at $t_1 = \frac{-b – \sqrt{\Delta}}{2a}$ and $t_2 = \frac{-b + \sqrt{\Delta}}{2a}$.

But wait! We don’t have one cone but two, so we have to reject solutions that intersect with the shadow cone. $P$ must still verify $\frac{ \vec{P} – \vec{C} }{\lVert \vec{P} – \vec{C} \rVert} \vec{V} = \cos\theta$, or simply, if $\theta < 90°$: $(\vec{P} – \vec{C})\cdot\vec{V} > 0$.

Note that there is also the corner case of the ray tangent to the cone and having an infinity of solutions to consider. I’ve completely swept it under the rug since it doesn’t matter in the context I was, but if it does to you, you’ve been warned about it. Also remember to check the sign of $t$ to know whether $P$ is in the direction of the ray. You may need to determine which of $t_1$ or $t_2$ you want to use, which depends on your use case. For example is your ray origin inside or outside of the cone?


Now for a little sanity test, let’s consider the corner case $C=O$, where the ray origin is the tip of the cone (thanks Rubix for the suggestion!). We have $b=0$ and $c=0$ thus $\Delta=0$ and $t=\frac{-b}{2a}=0$ which is the expected result.

I also tried the cases $\theta=0$ and $\theta=\pi/2$, but expanding $\Delta$ proved too tedious to proceed to the end. So this is left as an exercise, as they say. :)


Finally, to demonstrate that the result is indeed correct, here is a glorious ray traced cone scene on ShaderToy:

I hope this can prove useful to others too.

Oh, and Happy New Year by the way!

Various links on ray tracing

Here are some links related to ray tracing, and more specifically, path tracing.

Some ray tracing related projects or blogs:

Some major publications:

  • The rendering equation, SIGGRAPH 1986, James T. Kajiya. From the paper:

    We present an integral equation which generalizes a variety of known rendering algorithms.
    […]
    We mention that the idea behind the rendering equation is hardly new.
    […]
    However, the form in which we present this equation is well suited for computer graphics, and we believe that this form has not appeared before.

  • Bi-directional path tracing, Compugraphics 1993, Eric P. Lafortune and Yves D. Willems. From the paper:

    The basic idea is that particles are shot at the same time from a selected light source and from the viewing point, in much the same way. All hit points on respective particle paths are then connected using shadow rays and the appropriate contributions are added to the flux of pixel  in question.

  • Optimally Combining Sampling Techniques for Monte Carlo Rendering, SIGGRAPH 1995, Eric Veach and Leonidas J. Guibas. From the abstract:

    We present a powerful alternative for constructing robust Monte Carlo estimators, by combining samples from several distributions in a way that is provably good.

  • Metropolis Light Transport, SIGGRAPH 1997, Eric Veach and Leonidas J. Guibas. From the abstract:

    To render an image, we generate a sequence of light transport paths by randomly mutating a single current path (e.g. adding a new vertex to the path).

  • Robust Monte Carlo methods for light transport simulation, 1998, Erich Veach PhD thesis (432 pages pdf): it presents bidirectional path tracing, and introduces Metropolis Light Transport and Multiple Importance Sampling. From the abstract:

    Our statistical contributions include a new technique called multiple importance sampling, which can greatly increase the robustness of Monte Carlo integration. It uses more than one sampling technique to evaluate an integral, and then combines these samples in a way that is provably close to optimal. This leads to estimators that have low variance for a broad class of integrands. We also describe a new variance reduction technique called efficiency-optimized Russian roulette.

    […]

    The second algorithm we describe is Metropolis light transport, inspired by the Metropolis sampling method from computational physics. Paths are generated by following a random walk through path space, such that the probability density of visiting each path is proportional to the contribution it makes to the ideal image.

Other:

On a slightly different topic, fxguide had a great series of articles on the state of rendering in the film industry, which I previously mentioned.

John Carmack on Oculus at GDC 2015

John Carmack, the CTO of Oculus VR, gave a talk at the Game Developers Conference that just ended this week. Various topics are addressed, including the story behind Samsung’s Gear VR and what’s coming next, the democratization of virtual reality, the work on the API, the unsolved problem of controllers in VR, or the use of real-time ray tracing in VR.

John Carmack’s GDC 2015 talk.

It is a fairly long video (1h30), and as often with him, there are no pictures to see, just hear his personal views and insights on the work he is currently taking care of.

Relativistic and non euclidean space rendering

The Portal series built a full game concept out of non euclidean spaces. Besides being great games, I think it is fascinating how true the tagline “Now you’re thinking with portals” is.

Here are two interesting experiments putting the person in different spaces than we are used to due to real world conditions. This video by Varun Ramesh demonstrates a non-euclidean ray tracer:

This other video by the MIT Game Lab demonstrates OpenRelativity, a Unity toolkit allowing simulation of navigation at relativistic speeds, used for the prototype game A Slower Speed of Light:

Update: Sylvain mentioned in comments that Carl Sagan explains those effects in the following video: