# Implementation of VCM/UPS in Lightmetrica: Part 2

Posted on Dec 18, 2016Development Lightmetrica

(This article is written for Ray-tracing Advent Calendar 2016)

## Introduction of Part 2

Continuing the previous article, in this article we will introduce the implementation part of VCM/UPS in Lightmetrica.
Our code resides in `src/liblightmetrica/renderer/renderer_vcm.cpp`

(link) relative to the root directory of the repository.
We can refer this renderer from the configuration file with `renderer::vcm`

.
From now on we will explain the internal detail of the code.
For the sake of explanation we omit some unnecessary code for debugging, so the code snippets in this article are not exactly same as the original code.
We sometimes refer to the numbering of the equation in the previous articles, so it would be helpful to peek the previous article in the same time.

## VCM/UPS in Lightmetrica

### Code Structure

The structure of the code follows basic interface of renderers in Lightmetrica;
`Renderer_VCM`

class inherited from `Renderer`

class implements `Renderer::Initialize`

and `Renderer::Render`

interfaces.

```
1 class Renderer_VCM final : public Renderer
2 {
3 public:
4
5 LM_IMPL_CLASS(Renderer_VCM, Renderer);
6
7 private:
8
9 // Member variables
10
11 private:
12
13 LM_IMPL_F(Initialize) = [this](const PropertyNode* p) -> bool
14 {
15 // 1. Initialization
16 };
17
18 LM_IMPL_F(Render) = [this](const Scene* scene, Random* initRng, const std::string& outputPath) -> void
19 {
20 // 2. Rendering
21 };
22
23 };
```

### 1. Initialization

`Renderer::Initialize`

loads parameters of the renderer.
`max_num_vertices`

/ `min_num_vertices`

controls the number of path vertices (= path length + 1) to render. This parameter is useful for debugging, e.g., focusing on specific lighting effect only seem in the specific path length.
Our implementation based on the progressive photon density estimation by Knaus and Zwicker [2011]. Therefore the quality of the rendering is controlled by the `num_iteration_pass`

parameter, which specifies the number of iteration in rendering.
One iteration contains photon tracing and ray tracing parts. `num_photon_trace_samples`

and `num_eye_trace_samples`

controls the granularity of these process.
`initial_radius`

specifies the initial merge radius and `alpha`

controls its convergence ratio.
As well as the above parameters,
the `mode`

parameter can switch several fall-back implementations for BDPT or BDPM.
This parameter is mainly for debugging purpose.
In this article, we fixed it to `vcm`

.
The loaded parameters are stored in the class variables like `numIterationPass_`

which is referenced in the implementation.

### 2. Rendering

`Renderer::Initialize`

implements the core logic of the renderer.
The iteration loop is composed of four parts.

```
1 for (long long pass = 0; pass < numIterationPass_; pass++)
2 {
3 // 2-1. Update merge radius
4 // 2-2. Sample light subpaths
5 // 2-3. Construct range query structure
6 // 2-4. Estimate contribution
7 }
```

### 2-1. Update merge radius

Starting from the value specified by `initial_radius`

,
a merge radius used for the range query is updated according to the formulation in Knaus and Zwicker [2011]:

```
1 mergeRadius = pass == 0
2 ? initialRadius_
3 : std::sqrt((alpha_ + pass) / (1_f + pass)) * mergeRadius;
```

### 2-2. Sample light subpaths

In this part we sample light subpaths by the number specified by `num_photon_trace_samples`

.
A subpath is defined as a structure `VCMSubpath`

and sampled by `VCMSubpath::SampleSubpath`

function. The structure contains a vector of `VCMPathVertex`

which contains a path vertex type (`VCMPathVertex::type`

) represented by a bitwised flag defined with `SurfaceInteractionType`

, the surface geometry information such as the position or the normal (`VCMPathVertex::geom`

), and the reference to the primitive (`VCMPathVertex::primitive`

).

`VCMSubpath::SampleSubpath`

deligates the process to `PathSamplerUtils::TraceSubpath`

function, which is a general subpath sampler in our framework.
This function calls the given callback function once a new path vertex is generated.
Here we just copy the information and push to our vertex structure (line 16-26).
Note that we ignores `throughput`

argument which indicates a product of contributions of the sampled subpath so far. This is because of our design choice for our implementation; avoid optimizations.

```
1 struct VCMPathVertex
2 {
3 int type;
4 SurfaceGeometry geom;
5 const Primitive* primitive = nullptr;
6 };
7
8 struct VCMSubpath
9 {
10 std::vector<VCMPathVertex> vertices;
11 auto SampleSubpath(const Scene* scene, Random* rng, TransportDirection transDir, int maxNumVertices) -> void
12 {
13 vertices.clear();
14 PathSamplerUtils::TraceSubpath(
15 scene, rng, maxNumVertices, transDir,
16 [&](int numVertices, const Vec2& rasterPos,
17 const PathSamplerUtils::PathVertex& pv,
18 const PathSamplerUtils::PathVertex& v, SPD& throughput) -> bool
19 {
20 VCMPathVertex v_;
21 v_.type = v.type;
22 v_.geom = v.geom;
23 v_.primitive = v.primitive;
24 vertices.emplace_back(v_);
25 return true;
26 });
27 }
28 };
```

### 2-3. Construct range query structure

Based on the light subpaths generated in `2-2`

, we construct the range query structure.
This structure accelerates the query to obtain a light subpath from a position in the scene surfaces. A light subpath is selected if one of its vertices is within the range of given kernel size from the position. We used kd-tree for the structure implemented as `VCMKdTree`

class.

```
1 VCMKdTree pm(subpathLs);
2 pm.Build();
```

### 2-4. Estimate contribution

Next we trace eye subpaths and estimate intensities contributed to pixels. We used a standard parallelization function `Parallel::For`

in Lightmetrica, which here executes `numEyeTraceSamples_`

iterations in parallel. `Context`

defines thread-specific storage. The random number generator (`Random`

) must be assigned for each thread because it maintains its own internal states.
`Context::film`

accumulates the estimated contribution. Although the most of paths contributes to one pixel, for simplicity of code, we choose a design to accumulate all contribution to one film instead of separating two films as said in Veach’s paper [1997].

After the rendering completes, the films for each of threads are accumulated into one. Here we compute the average of films in `numIterationPass_`

iterations in the progressive way (line 19-24).

```
1 struct Context
2 {
3 Random rng;
4 Film::UniquePtr film{nullptr, nullptr};
5 };
6 std::vector<Context> contexts(Parallel::GetNumThreads());
7 for (auto& ctx : contexts)
8 {
9 ctx.rng.SetSeed(initRng->NextUInt());
10 ctx.film = ComponentFactory::Clone<Film>(film);
11 ctx.film->Clear();
12 }
13 Parallel::For(numEyeTraceSamples_, [&](long long index, int threadid, bool init)
14 {
15 // 2-4-1. Sample subpaths
16 // 2-4-2. Estimate contribution by vertex connection
17 // 2-4-3. Estimate contribution by vertex merging
18 });
19 film->Rescale((Float)(pass) / (1_f + pass));
20 for (auto& ctx : contexts)
21 {
22 ctx.film->Rescale(1_f / (1_f + pass));
23 film->Accumulate(ctx.film.get());
24 }
```

### 2-4-1. Sample subpaths

First we sample subpaths with both directions.
These subpaths are later used for the vertex connection and the vertex merging.
Note that specifically for the vertex connection, we can reuse previously sampled light subpaths in `2-2`

instead of sampling a new one. This, however, needs to compensate correlation among the reused paths and it is dependent on the number of light subpaths sampled for each iteration. Following the policy of not doing optimizations, here we did not adopt that way.

```
1 VCMSubpath subpathE;
2 VCMSubpath subpathL;
3 subpathE.SampleSubpath(scene, &ctx.rng, TransportDirection::EL, maxNumVertices_);
4 subpathL.SampleSubpath(scene, &ctx.rng, TransportDirection::LE, maxNumVertices_);
```

### 2-4-2. Estimate contribution by vertex connection

In this part we estimate contribution by vertex connection (`$\langle I \rangle_{VC}$`

in Eq.8).
The basic structure follows double loop for `s`

and `t`

, each corresponds to the number of vertices in light and eye subpaths.
The ranges of these values are clipped so that the number of vertices in the connected path is in the range of max(2, `minNumVertices_`

) and min(`nL`

+`nE`

, `maxNumVertices_`

), where `nE`

and `nL`

is the number of vertices in the sampled subpaths.

```
1 const int nE = (int)(subpathE.vertices.size());
2 for (int t = 1; t <= nE; t++)
3 {
4 const int nL = (int)(subpathL.vertices.size());
5 const int minS = Math::Max(0, Math::Max(2 - t, minNumVertices_ - t));
6 const int maxS = Math::Min(nL, maxNumVertices_ - t);
7 for (int s = minS; s <= maxS; s++)
8 {
9 // 2-4-2-1. Connect vertices and create a full path
10 // 2-4-2-2. Evaluate contribution
11 // 2-4-2-3. Evaluate path PDF
12 // 2-4-2-4. Evaluate MIS weight
13 // 2-4-2-5. Accumulate contribution
14 }
15 }
```

### 2-4-2-1. Connect vertices and create a full path

After fixing `t`

and `s`

, we connect two subpaths and construct a *full* path.
The actual logic for path connection is implemented in `VCMPath::ConnectSubpath`

function.
For both-positive `s`

and `t`

, the connection process needs to check visibility between `s`

-th vertex in the light subpath and `t`

-th vertex in the eye subpath. `Scene::Visible`

function checks this (line 26).
We do not need this visibility query if one of `s`

or `t`

is zero, indicating a path sampling strategy that directly hit emitters. For instance, `s=0`

corresponds to path tracing. In this case, we only need to check the last vertex in the corresponding direction is on an emitter.
If the path connection fails, in both cases the evaluation process is terminated because the resulting paths do not carry any contribution.

```
1 VCMPath fullpath;
2 if (!fullpath.ConnectSubpaths(scene, subpathL, subpathE, s, t)) { continue; }
```

```
1 struct VCMPath
2 {
3 ...
4 auto ConnectSubpaths(const Scene* scene, const VCMSubpath& subpathL, const VCMSubpath& subpathE, int s, int t) -> bool
5 {
6 assert(s >= 0);
7 assert(t >= 0);
8 vertices.clear();
9 if (s == 0 && t > 0)
10 {
11 vertices.insert(vertices.end(), subpathE.vertices.rend() - t, subpathE.vertices.rend());
12 if ((vertices.front().primitive->Type() & SurfaceInteractionType::L) == 0) { return false; }
13 vertices.front().type = SurfaceInteractionType::L;
14 }
15 else if (s > 0 && t == 0)
16 {
17 vertices.insert(vertices.end(), subpathL.vertices.begin(), subpathL.vertices.begin() + s);
18 if ((vertices.back().primitive->Type() & SurfaceInteractionType::E) == 0) { return false; }
19 vertices.back().type = SurfaceInteractionType::E;
20 }
21 else
22 {
23 const auto& vL = subpathL.vertices[s - 1];
24 const auto& vE = subpathE.vertices[t - 1];
25 if (vL.geom.infinite || vE.geom.infinite) { return false; }
26 if (!scene->Visible(vL.geom.p, vE.geom.p)) { return false; }
27 vertices.insert(vertices.end(), subpathL.vertices.begin(), subpathL.vertices.begin() + s);
28 vertices.insert(vertices.end(), subpathE.vertices.rend() - t, subpathE.vertices.rend());
29 }
30 return true;
31 }
32 ...
33 };
```

### 2-4-2-2. Evaluate contribution

Here we evaluate the measurement contribution ($f$ in Eq.1) for the connected (full) path. This logic is implemented in `VCMPath::EvaluateF`

function.
This function takes arguments `s`

and `merge`

.
Some might be temped to evaluate this function always with some fixed strategy (e.g., `s=0`

) irrespective to the actual sampling strategies, but we cannot because of the following reason.

First, our code of the BSDF evaluation (`BSDF::EvaluateDirection`

) compensates for non-symmetries for shading normals and specular refractions [Veach 1997]. This makes the measurement contribution function to be dependent on the sampling strategy.
This means we must match the strategy of evaluating the same strategy with which the path is sampled.

Second, in the case that a path contains vertices on specular surfaces, we implicitly handles delta functions in the BSDF evaluations. In the sequence of the evaluation of contribution, the delta functions will be canceled out except for the BSDF evaluations related to connecting vertices. In order to compensate for this oddity, our implementation introduces special flag to ignore delta functions when evaluating BSDFs, which again makes it dependent on the sampling strategy.

`2-4-2-2-1`

and `2-4-2-2-2`

computes sequence of `$f_s\cdot G$`

from the edges of E and L respectably. `2-4-2-2-3`

takes care of the remaining terms.
Note that we omits the visibility check if the `merge`

flag is enabled.

```
1 const auto f = fullpath.EvaluateF(s, false);
2 if (f.Black()) { continue; }
```

```
1 struct VCMPath
2 {
3 ...
4 auto EvaluateF(int s, bool merge) const -> SPD
5 {
6 const int n = (int)(vertices.size());
7 const int t = n - s;
8 assert(n >= 2);
9
10 // 2-4-2-2-1. Evaluate product of f_s from L
11 // 2-4-2-2-2. Evaluate product of f_s from E
12 // 2-4-2-2-3. Evaluate c_{s,t}
13
14 return fL * cst * fE;
15 }
16 ...
17 };
```

### 2-4-2-3. Evaluate path PDF

Next we will evaluate PDF for the connected path (`$p_{s,t}^{VC}$`

in Eq.2 or `$p_{s,t}^{VM}$`

in Eq.7), which is implemented in `VCMPath::EvaluateF`

function.
This function first checks if the given path is samplable with the selected strategy (`2-4-2-3-1`

).
If the path is determined not to be samplable, this function returns zero pdf value.
For instance, a path originated from the geometrically-degenerated emitters such as a point light source cannot be sampled with path tracing (the strategy with `s=0`

).
We need to check this because this function can be called with different strategy as it sampled with, later in the evaluation of MIS weights.

Once the path is checked as samplable with given strategy, the function evaluates sequence of pdf defined for each vertex of the path in (`2-4-2-3-2`

). This evaluate process also depends on the strategy as well as the `merge`

flag.
Note that the `Primitive::EvaluateDirectionPDF`

functions returns pdf in projected solid angle measure, so we need to convert it to the area measure with `PDFVal::ConvertToArea`

function because the target density is defined on the product area measure.
We also note that when `merge`

flag is enabled, we need to evaluate `$\pi r^2$`

in Eq.7.

```
1 const auto p = fullpath.EvaluatePathPDF(scene, s, false, 0_f);
2 if (p.v == 0_f) { return; }
```

```
1 struct VCMPath
2 {
3 ...
4 auto EvaluatePathPDF(const Scene* scene, int s, bool merge, Float radius) const -> PDFVal
5 {
6 const int n = (int)(vertices.size());
7 const int t = n - s;
8 assert(n >= 2);
9
10 // 2-4-2-3-1. Check if the path is samplable by the selected strategy
11 // 2-4-2-3-2. Evaluate path PDF
12
13 return pdf;
14 }
15 ...
16 };
```

### 2-4-2-4. Evaluate MIS weight

In this part we will evaluate MIS weights (`$w_{v,s,t}$`

in Eq.9).
Here `numPhotonTraceSamples`

corresponds to `$N_{VM}$`

.
In line 7, the numerator of Eq.9 evaluates to `ps`

.
The double loop in the next part evaluates the reciprocal of the weight, considering numerical goodness. Our implementation adopts the power heuristic with `$\beta=2$`

so the ratio of `pi`

and `ps`

is multiplied twice in line 19.
Some strategies might not be possible to sample the given path, in this case `EvaluatePathPDF`

function evaluates to zero and omits from the evaluation of the weight (line 16).

```
1 const auto w = fullpath.EvaluateMISWeight_VCM(
2 scene, s, false, mergeRadius, numPhotonTraceSamples_);
```

```
1 struct VCMPath
2 {
3 ...
4 auto EvaluateMISWeight_VCM(const Scene* scene, int s_, bool merge, Float radius, long long numPhotonTraceSamples) const -> Float
5 {
6 const int n = static_cast<int>(vertices.size());
7 const auto ps = EvaluatePathPDF(scene, s_, merge, radius);
8 assert(ps > 0_f);
9
10 Float invw = 0_f;
11 for (int s = 0; s <= n; s++)
12 {
13 for (int type = 0; type < 2; type++)
14 {
15 const auto pi = EvaluatePathPDF(scene, s, type > 0, radius);
16 if (pi > 0_f)
17 {
18 const auto r = pi.v / ps.v;
19 invw += r*r*(type > 0 ? (Float)(numPhotonTraceSamples) : 1_f);
20 }
21 }
22 }
23
24 return 1_f / invw;
25 }
26 ...
27 };
```

### 2-4-2-5. Accumulate contribution

Hooray! We finally set up necessary terms (`f`

,`w`

, and `p`

) to compute the contribution. Here we compute it and accumulate to the film (`Film::Splat`

function). This function takes a raster position in `$[0,1]^2$`

and the contribution to be accumulated. `VCMPath::RasterPosition`

computes the raster position of the path.

```
1 const auto C = f * w / p;
2 ctx.film->Splat(
3 fullpath.RasterPosition(),
4 C * (Float)(film->Width() * film->Height()) / (Float)numEyeTraceSamples_);
```

```
1 struct VCMPath
2 {
3 ...
4 auto RasterPosition() const -> Vec2
5 {
6 const auto& v = vertices[vertices.size() - 1];
7 const auto& vPrev = vertices[vertices.size() - 2];
8 Vec2 rasterPos;
9 v.primitive->sensor->RasterPosition(Math::Normalize(vPrev.geom.p - v.geom.p), v.geom, rasterPos);
10 return rasterPos;
11 }
12 ...
13 };
```

### 2-4-3. Estimate contribution by vertex merging

In this part, we estimate the contribution by vertex merging (`$\langle I \rangle_{VM}$`

in Eq.8).
Note that this has very similar structure to `2-4-2`

. The first loop for `t`

is same as before, and the second loop is replaced with range query. `VCMKdTree::RangeQuery`

function takes a function that is called when a recorded position within the given radius from the given position is found. `si`

corresponds to the index of the light subpaths in `2-2`

, and `vi`

points to its vertex.
Note that the range query takes responsibility for the operation `$\sum_{l=1}^{N_{VM}} \sum_{s,t\geq 2} \chi_{s,t}(\bar{x}_l)$`

in Eq.8.

Line 4 skip density estimation for the point geometries, although it rarely happens in the most of the scenes.
We omits the explanation for the parts `2-4-3-2`

to `2-4-3-5`

because it uses same function with changing parameters as in `2-4-2`

.

```
1 for (int t = 1; t <= nE; t++)
2 {
3 const auto& vE = subpathE.vertices[t - 1];
4 if (vE.primitive->IsDeltaPosition(vE.type)) { continue; }
5 pm.RangeQuery(vE.geom.p, mergeRadius, [&](int si, int vi) -> void
6 {
7 const int s = vi + 1;
8 const int n = s + t - 1;
9 if (n < minNumVertices_ || maxNumVertices_ < n) { return; }
10 // 2-4-3-1. Merge vertices and create a full path
11 // 2-4-3-2. Evaluate contribution
12 // 2-4-3-3. Evaluate path PDF
13 // 2-4-3-4. Evaluate MIS weight
14 // 2-4-3-5. Accumulate contribution
15 });
16 }
```

### 2-4-3-1. Merge vertices and create a full path

Once again after fixing the strategies (`s`

and `t`

), we construct the full path from the light and eye subpaths. We want to be careful of the index of the strategy, the strategy indexed by `(s,t)`

corresponds to the path with `s+t-1`

vertices, not `s+t`

vertices as in the case of vertex connection. And the number of vertices passed to the following functions must be `s-1`

, not `s`

.
`VCMPath::MergeSubpaths`

takes care of the actual process. Note that we no longer needs to check visibility between vertex so the things becomes easier. Again we skip the merging for point geometries in line 11.

```
1 static thread_local VCMPath fullpath;
2 if (!fullpath.MergeSubpaths(subpathLs[si], subpathE, s - 1, t)) { return; }
```

```
1 struct VCMPath
2 {
3 ...
4 auto MergeSubpaths(const VCMSubpath& subpathL, const VCMSubpath& subpathE, int s, int t) -> bool
5 {
6 assert(s >= 1);
7 assert(t >= 1);
8 vertices.clear();
9 const auto& vL = subpathL.vertices[s - 1];
10 const auto& vE = subpathE.vertices[t - 1];
11 if (vL.primitive->IsDeltaPosition(vL.type) || vE.primitive->IsDeltaPosition(vE.type)) { return false; }
12 if (vL.geom.infinite || vE.geom.infinite) { return false; }
13 vertices.insert(vertices.end(), subpathL.vertices.begin(), subpathL.vertices.begin() + s);
14 vertices.insert(vertices.end(), subpathE.vertices.rend() - t, subpathE.vertices.rend());
15 return true;
16 }
17 ...
18 };
```

## Conclusion

Concluding the article, we did not have time to explain all of the code and omits some explanation of the details. So it would be great to read code further for better understands, and of course, feel free to ask me any questions! As said again and again in this article, our implementation does not focus on the optimizations, even on the utterly apparent ones. Please be careful if you adopt (some part of) the code in your productions. We plan to implement the optimized version of VCM/UPS someday, possibly.

Thank you for reading!

## References

- [Veach 1997] E. Veach, Robust Monte Carlo method for light transport simulation, PhD Thesis, Stanford University, 1997.
- [Knaus and Zwicker 2011] Progressive photon mapping: A probabilistic approach, In
*ACM Transactions on Graphics*, 2011.