Attempt to use symbolic computation for renderer verification

Posted on Jul 7, 2017
Experimental

(This article is written for ray-tracing camp 5 Advent Calender)

Introduction

Despite its simpleness in theory, implementing rendering systems is known to be difficult. For instance, although bidirectional path tracing (BDPT) [Veach and Guibas 1994] is relatively simple in concept and formulation, it is deceivingly difficult to implement correctly. The cause of the difficulties lies the fact that the rendering systems have a limited verifiability. Due to the use of Monte Carlo method, the output of the renderers are almost always indeterministic. Many exceptional handling slipped from the formulation could make the implementation of the renderer complex and difficult. Taking an example of BDPT, handling pinhole camera is an exception of the formulation as the position of the camera is spatially degenerated (i.e, a point). The probability of sampling the position on the camera through direction sampling is almost surely zero so we need to care about the weight computation not to include the probability from the strategy when we consider the paths generated with the other strategies.

Unfortunately, as far as I know, there is no silver bullet for the difficulties. One can of course use various standard debugging testing in software development, like print debugging or using debuggers, or adopt the software development techniques like unit testing. Unfortunately, however, any of them can partially solve the difficulties. Print debugging is clearly out of its depth when we realize the fact that the result of the rendering is essentially the accumulation of the millions of iterations. Debuggers are rather useful as we can make use richer runtime information and the flow control, but essentially their capability is similar to that of print debugging. More importantly and crucially, neither of them can solve the systematic errors, that is, error on the specification, algorithms, etc. caused by the human him/herself. Basically, the developers just need to confront there difficulties with the combination of several approaches.

In this article, we focus on one aspect of the difficulties: the difference between the formulation and the actual code. Rendering techniques are typically describes with a mathematical formulation rather than the code. So if we want to implement the technique, we need to interpret the formulation before implementation. The systematic error tends to contaminate the code in this interpretation phase. A typical way to retain correctness is to compare the result with the result with the another implementation. If we want to implement the famous techniques like path tracing, it is relatively easy because we can readily find the existing open source implementations. Even if the technique is unknown, we can at least compare the result with another existing techniques because the rendering eventually solves the same governing equation.

This kind of comparison-with-result approaches is rather top-down, as we keep verifiability based on the final result. The top-down approaches, however, find it hard to implement the class of the techniques that requires to implement several components in order to obtain the final result and that each of components is hard to verify by themselves. At best, in such a case, the top-down approach can only find the existence of the defects. So, here we want to think about the alternative approach, say, the bottom-up approach. The bottom-up approach compares the result per component. The reference for the comparison is, for instance, the analytical solution of the simple cases which has been usually done with the hand-written calculations.

In this article, we will introduce a way to obtain the reference in the bottom-up approach with the symbolic computation. The symbolic computation is a way of computation by manipulating mathematical objects symbolically, rather than numerically. We will demonstrate the usefulness of the symbolic computation in the context of renderer implementations with several examples. The examples focus on the measure conversion. The operation frequently happens in the renderer code and its implementation is highly error-prone because it involves in many derivative calculations.

This article is inspired by the recently published paper [Anderson et al. 2017]. They developed an embedded domain specific language which facilities the symbolic computation to automatically generate the code from the sampler. The approach in more general context is also published recently [Bhat et al. 2017]. Similar to these researches, in this article, we will peek the power of symbolic computation and attempt to apply the technology to the well-known computations in the renderer computation, which have been done with the hand-written calculations. Ultimately, we want to integrate the symbolic computation with the verification process of the renderer development. I hope this article becomes a step toward the goal.

Symbolic computation

The symbolic computation, a.k.a. computer algebra, is a way of computation by manipulating mathematical expressions in the symbolic way. Conversely and usually, the standard way to compute something in the scientific computation is to use the numerical computation. For instance, if we want to compute a mathematical expression $\sin^2{\theta}+\cos^2{\theta}$ using computers. In the numerical computation, the expression can be computed if the value assigned for $\theta$ is given. The result should be near to one irrespective to the value of $\theta$, but not exactly one because of the numerical errors.

On the other hand, the symbolic handles the mathematical expression without binding the actual value into the expressions. The computation proceeds with changing the shape of the expressions. With the example above, the expression can be converted as $\sin^2{\theta}+\cos^2{\theta}=1$ using the trigonometric identity. This is how the symbolic computation goes, and, as you know, this is what we usually do with paper and pen!

I want to note that the symbolic and numerical computation are not the complementary concepts. We can even compute the final value using numerical computation using the partial result of the symbolic computations. This way of thinking is useful when we want to actually apply the symbolic computation in the real usage.

Applications range among the various field of the scientific computation. Various computer algebra systems, such as Mathematica, are used to solve the various problem in the real world. Needless to say, the rendering is one of the applications of the scientific computation. There is no reason not to use symbolic computation in rendering.

In the following section, we will explain the usage with the actual code, using a python code written with the SymPy library. SymPy is a python library for symbolic computation. It has well-defined capability to handle typical usage, such as basic symbolic arithmetic, calculus or algebra. The library has good integration into the Jupiter notebook, with pretty formatting with LaTeX. And what is better, SymPy is free! (licensed under New BSD License)

Application to renderer verification

The purpose of the article is to attempt the symbolic computation to the renderer verification. In this article, we focus on one mathematical aspect that is often observed in the formulation of the rendering techniques: change of probability densities through transformations.

The light transport is written as a integral w.r.t. to some measure. For instance, using project solid angle in the famous spherical formulation, the outgoing radiance can be written as an integral w.r.t. the solid angle measure $\sigma$: $$ \begin{equation} L_o(\mathbf{x},\omega_o) = L_e(\mathbf{x},\omega_o) + \int_{\Omega} f_r(\mathbf{x},\omega,\omega_o) L_i(\mathbf{x},\omega) \cos{\theta}\; d\sigma(\omega), \end{equation} $$ skipping the explanation of the terms. In the following description, we will use the shorthand notations such as $L_o\equiv L_o(\mathbf{x},\omega_o)$, $f_r(\omega)\equiv f_s(\mathbf{x},\omega,\omega_o)$, etc.

The Monte Carlo estimate of the integral should match the measure of the density. In other words, the probability density function must be defined on the same measure as the target integral. In this example, therefore, an estimate can be written as $$ \begin{equation} L_o \approx \langle L_o \rangle \equiv L_e + \frac{1}{N}\sum_{i=1}^N \frac{f_r(\omega_i)L_i(\omega_i)\cos{\theta}}{p_\sigma(\omega_i)}. \end{equation} $$ Here $p_\sigma$ is the pdf w.r.t. the solid angle measure $\sigma$, which is same as the original integral.

However, it is often the case that the density is not directly sampled from the distribution with the target measure. For instance, the direction is often sampled from the uniform distribution (with Lebesgue measure $\lambda$) but the pdf must actually be evaluated in the solid angle measure. $$ \begin{equation} Pr(\omega\in\Omega) = \int_\Omega p_\sigma(\omega) d\sigma(\omega) = \int_{\mathcal{U}} \underbrace{p_\sigma(\omega(\mathbf{u})) \frac{d\sigma}{d\lambda}}_{p_\lambda(\mathbf{u})} d\lambda(\mathbf{u}) = Pr(\mathbf{u}\in\mathcal{U}) \end{equation} $$ We need to compute the derivative $\frac{d\sigma}{d\lambda}$ to fill in the gap between the sampled density and the measure of the target integral. In the following examples, we will introduce several examples of the derivatives like this using the symbolic computation.

Example: Uniform sampling on a unit disk

Let us consider the transformation of densities from numbers taken from the uniform distributions $u_1, u_2 \sim \mathcal{U}(0,1)$ and an unit disk parametrized by $\phi\in [0,2\pi]$ and $r\in [0,1]$. Here we want to consider polar mapping between two spaces. The relationship between $(u_1, u_2)$ and $(\phi, r)$ is written as $$ \begin{equation} \phi = 2\pi u_1,\; r = \sqrt{u_2}. \end{equation} $$ Applying change-of-variable formula, the integral on the space of $(\phi,r)$ reads $$ \begin{equation} \int p(\phi,r) d\phi dr = \int \underbrace{p(\phi(u_1,u_2),r(u_1,u_2)) \left|\frac{\partial(\phi,r)}{\partial(u_1,u_2)}\right|}_{p(u_1, u_2)} du_1du_2. \end{equation} $$ Here the Jacobian can be calculated as $$ \begin{equation} \left|\frac{\partial(\phi,r)}{\partial(u_1,u_2)}\right| = \begin{vmatrix} \frac{\partial \phi}{\partial u_1} & \frac{\partial \phi}{\partial u_2} \\ \frac{\partial r}{\partial u_1} & \frac{\partial r}{\partial u_2} \end{vmatrix} = \begin{vmatrix} 2\pi & 0 \\ 0 & \frac{1}{2}{u_2}^{-1/2} \end{vmatrix} = \pi{u_2}^{-1/2}. \end{equation} $$ Again, we consider the transform $(\phi,r)$ to $(x,y)\in\mathbb{R}^2$: $$ \begin{equation} x = r\cos{\phi},\; y = r\sin{\phi}. \end{equation} $$ Similarly, applysing change-of-variable formula, we obtain $$ \begin{equation} \int p(x,y) dxdy = \int \underbrace{p(x(\phi,r),y(\phi,r)) \left|\frac{\partial(x,y)}{\partial(r,\phi)}\right|}_{p(r,\phi)} d\phi dr. \end{equation} $$ The Jacobian is $$ \begin{equation} \left|\frac{\partial(x,y)}{\partial(\phi,r)}\right| = \begin{vmatrix} \frac{\partial x}{\partial \phi} & \frac{\partial x}{\partial r} \\ \frac{\partial y}{\partial \phi} & \frac{\partial y}{\partial r} \end{vmatrix} = \begin{vmatrix} -r\sin{\phi} & \cos{\phi} \\ r\cos{\phi} & \sin{\phi} \end{vmatrix} = r. \end{equation} $$ Applying chain rule, we can obtain the uniform density $$ \begin{equation} p(x,y) = \left( \left|\frac{\partial(x,y)}{\partial(r,\phi)}\right| \left|\frac{\partial(r,\phi)}{\partial(u_1,u_2)}\right| \right)^{-1}p(u_1,u_2) = \left( r \cdot \frac{\pi}{\sqrt{u_2}} \right)^{-1} \cdot 1 = \frac{1}{\pi} \end{equation} $$ as desired.

Let’s mimic this process with SymPy. In the first place, we will import SymPy and enable pretty printing with LaTeX. If you try this script in Jupyter notebook, the output of the mathematical expressions should be printed with LaTeX typesetting.

1 from sympy import *
2 init_printing()

Next we will declare the variable for $u_1$ and $u_2$ with symbol command. The variables u1 and u2 are just python variables. In the following script, the mathematical expression can be written as expressions in python.

1 u1, u2 = symbols('u_1 u_2')

We can express the transformations $\phi(u_1,u_2)$ and $r(u_1,u_2)$ with python expressions.

1 rho_u1u2 = 2 * pi * u1
2 r_u1u2 = sqrt(u2)

Fortunately, SymPy has a library to do various matrix operations, including computing Jacobian matrix. Using this feature, we can compute the Jacobi matrix in Eq.6 as

1 J_rhophi_u1u2 = Matrix([rho_u1u2,r_u1u2]).jacobian([u1,u2])

Given the Jacobi matrix, the computation of Jacobian is straightforward using predefined function det to compute the determinant.

1 J_rhophi_u1u2_det = abs(J_rhophi_u1u2.det())

You will notice that displaying J_rhophi_u1u2_det indicates $\pi\left|1 / \sqrt{u_2}\right|$ instead of $\pi / \sqrt{u_2}$. We can remove the unnecessary outer $|\cdot|$ using refine function with the assumption of the variable. The assumption can be made using the context with assuming(...):. Here we only impose an assumption: $u_2>0$.

1 with assuming(Q.positive(u2)):
2     J_rhophi_u1u2_det = refine(J_rhophi_u1u2_det)

Similarly, we can compute the Jacobian for the transformation $x(\phi,r), y(\phi,r)$.

1 phi, r = symbols('\\phi r')
2 x_phir = r * cos(phi)
3 y_phir = r * sin(phi)
4 J_xy_phir = Matrix([x_phir,y_phir]).jacobian([phi,r])

Again if we display J_xy_rhir_det we obtain $\left|r\sin^2{\phi}+r\cos^2{\phi}\right|$ instead of $r$. The expression containing trigonometric identities like this can be simplified with trigsimp function. Also we impose the positive assumption for $r$ to remove $|\cdot|$.

1 with assuming(Q.positive(r)):
2     J_xy_phir_det = refine(trigsimp(J_xy_phir_det))

Finally we can apply chain rule with multiplying two Jacobians. The expression J_xy_rhir_det is a function of phi and r so in order for the cancel out to happen properly phi and r must be substituted by rho_u1u2 and r_u1u2 respectively. The substitution can be done with subs function.

1 J_xy_u1u2_det = J_xy_phir_det.subs({rho: rho_u1u2, r: r_u1u2}) * J_rhophi_u1u2_det

Displaying J_xy_u1u2_det shall read $\pi$ as expected.

Example: Surface integral

The integral in the rendering context often reads the surface integral, e.g. the integral w.r.t. solid angle measure, or area measure. In this example we will derive the change of the measure between the solid angle $\omega$ and the spherical coordinates $(\theta,\phi)$ on the unit sphere: $$ \begin{equation} d\omega = \sin{\theta}\; d\theta d\phi \end{equation} $$ We cannot directly apply the change-of-variable formula as in the previous example because the dimension changes with the transformation. Instead, this relationship is usually obtained by observing the geometric relationship between differential elements. This way of derivation, however, cannot be directly handled in the symbolic computation.

Generalized change-of-variable formula with matrix volume

In order to make it possible to handle this kind of the surface integral in the symbolic computation system, in this section, we will introduce the generalized way to change measures between the space of different dimensions [Ben-Israel 1999].

Reviewing the basic calculus, the change-of-variable formula for the transformation $\phi: \mathcal{U}\to\mathcal{V}$, where $\mathcal{U}$ and $\mathcal{V}$ is a set in $\mathbb{R}^n$, is written as $$ \begin{equation} \int_\mathcal{V} f(\mathbf{v}) d\mathbf{v} = \int_\mathcal{U} f(\phi(\mathbf{v}(\mathbf{u}))) \left| J_\phi(\mathbf{u}) \right| d\mathbf{u}, \end{equation} $$ where $J_\phi$ is the Jacobi matrix and $\left| \cdot \right|$ is the determinant. The generalized formula extends it in the case that $\mathcal{U}$ and $\mathcal{V}$ is in the space of different dimensions, specifically, $\mathcal{U}\in\mathcal{R}^n$ and $\mathcal{V}\in\mathcal{R}^m$ with $n>m$, by replacing $\left| J_\phi\right|$ with the matrix volume $\mathrm{vol}\,J_\phi$: $$ \begin{equation} \int_\mathcal{V} f(\mathbf{v}) d\mathbf{v} = \int_\mathcal{U} f(\phi(\mathbf{v}(\mathbf{u})))\; \mathrm{vol}\, J_\phi(\mathbf{u})\; d\mathbf{u}. \end{equation} $$ If the Jacobi matrix $J_\phi$ is full rank, the volume is written as a $m\times n$ matrix: $$ \begin{equation} \mathrm{vol}\,J_\phi = \sqrt{\det J_\phi^T J_\phi }, \end{equation} $$ or generally with an arbitrary-ranked matrix $J_\phi$: $$ \begin{equation} \mathrm{vol}\,J_\phi = (\text{Product of the singular values of $J_\phi$}). \end{equation} $$ This quantity is computationally feasible in the symbolic computation system. Many of the system already have a feature to compute singular values as a library.

Derivation

We will derive Eq.11 using the formula. Let $(x,y,z)$ be the points on the surfaces of the unit sphere. The transformation from the spherical coordinate on the unit sphere $(\theta,\phi)$ to $(x,y,z)$ is written as $$ \begin{equation} x = \sin{\theta}\cos{\phi},\; y = \sin{\theta}\sin{\phi},\; z = \cos{\theta}. \end{equation} $$ The Jacobi matrix is written as $$ \begin{equation} \frac{\partial(x,y,z)}{\partial(\theta,\phi)} = \begin{bmatrix} \frac{\partial x}{\partial \theta} & \frac{\partial x}{\partial \phi} \\ \frac{\partial y}{\partial \theta} & \frac{\partial y}{\partial \phi} \\ \frac{\partial z}{\partial \theta} & \frac{\partial z}{\partial \phi} \end{bmatrix} = \begin{bmatrix} \cos{\theta}\cos{\phi} & -\sin{\theta}\sin{\phi} \\ \cos{\theta}\sin{\phi} & \sin{\theta}\cos{\phi} \\ -\sin{\theta} & 0 \end{bmatrix}. \end{equation} $$ Applying Eq.14, the volume can be calculated as $$ \begin{equation} \mathrm{vol}\,\frac{\partial(x,y,z)}{\partial(\theta,\phi)} = \begin{vmatrix} 1 & 0 \\ 0 & \sin^2{\theta} \end{vmatrix}^{\frac{1}{2}} = \sin{\theta}. \end{equation} $$ Therefore, we obtain $$ \begin{equation} \int f(x,y,z) d\omega = \int f(x(\theta,\phi),y(\theta,\phi),z(\theta,\phi)) \sin{\theta}\, d\theta d\phi, \end{equation} $$ as desired.

Let’s follow this process with SymPy. The step until the computation of the Jacobi matrix is similar to the previous example; just declare the variables $\theta$ and $\phi$ with transformations and compute the Jacobi matrix.

1 theta, phi = symbols('\\theta \\phi')
2 x_thetaphi = sin(theta) * cos(phi)
3 y_thetaphi = sin(theta) * sin(phi)
4 z_thetaphi = cos(theta)
5 J_xyz_thetaphi = Matrix([x_thetaphi,y_thetaphi,z_thetaphi]).jacobian([theta,phi])

Given that the computed Jacobi matrix is full rank, we can obtain the volume using Eq.14.

1 J_xyz_thetaphi_vol = sqrt((J_xyz_thetaphi.T * J_xyz_thetaphi).det())

Or alternatively and generally, we can also use Eq.15. SymPy has a useful function singular_value to compute the values. The function requires the assumption that the variables are real number. However, possibly due to a bug, unfortunately, the assumption imposed by with assuming(...) does not apply to the result of the function singular_values. So we need use instead impose the assumption in the symbols function.

1 theta, phi = symbols('\\theta \\phi', real=True)

With this assumption, we can faithfully compute the volume in the generalized manner.

1 from functools import reduce
2 J_xyz_thetaphi_SVs = J_xyz_thetaphi.singular_values()
3 J_xyz_thetaphi_vol = reduce(lambda x,y: x*y, J_xyz_thetaphi_SVs, 1)

The computed value J_xyz_thetaphi_vol again contains some trigonometric functions. Let’s simplify them.

1 J_xyz_thetaphi_vol = trigsimp(J_xyz_thetaphi_vol)

In so far, J_xyz_thetaphi_vol should reads $|\sin{\theta}|$.We are tempted to remove $|\cdot|$ with the assumption that $\theta\in[0,\pi]$:

1 with assuming(Q.positive(theta), Q.negative(theta-pi)):
2 	J_xyz_thetaphi_vol = refine(J_xyz_thetaphi_vol)

SymPy, however, does not support this way of indirect simplification. So we instead give the assumption that $\sin{\theta}>0$.

1 with assuming(Q.positive(sin(theta))):
2 	J_xyz_thetaphi_vol = refine(J_xyz_thetaphi_vol)

Finally, J_xyz_thetaphi_vol displays $\sin{\theta}$ as expected.

Conclusion

In this article, we introduced the use of symbolic computation on the rendering system. We pointed out the usefulness of the symbolic computation on the renderer verification. As an initial attempt, we gave several examples to verify the conversion the measure probability densities between different measures.

As the complexity of rendering system increases, we believe the verification of the system would be more important. I believe, definitely, the use of symbolic computation can alleviate a certain class of difficulties in the renderer implementations. As a future work, I want to investigate the practical uses with rendering images, and the ways to integrate the symbolic computation into the actual development process to enhance the research.

References