\section{Unbiased Monte Carlo Estimator}
\label{sec:mc}
%
We now show how to numerically estimate the scene derivative $\dev\Lv$ of interior radiance $\Lv$.
To this end, we show in \S\ref{ssec:our_est} how individual components $\dLv := \sdev\Lv$ can be computed based on our derivations from \S\ref{sec:diff_rte}.
Then, we discuss the differentiation of pixel intensities in \S\ref{ssec:pixel}.
Lastly, we discuss handling of the total derivative $\dev\Lv$ in \S\ref{ssec:mc_discussion}.

\input{fig_dependency.tex}
\input{fig_comps.tex}

\subsection{Monte Carlo Estimation of $\dLv$}
\label{ssec:our_est}
%
To estimate $\dLv$ for some scene parameter $\param \in \bparam$, we introduce a Monte Carlo technique analogous to unidirectional volumetric path tracing (VPT) that handles $\Lv$ and $\dLv$ jointly.

\subsubsection{Next-Event Estimation.}
%
A variance-reduction technique essential to the efficiency of VPT is next-event estimation (NEE): instead of only ``tallying'' when a light transport path traced from the virtual sensor hits a light source, NEE-based methods tally at each volumetric and interfacial scattering location.
We now show how NEE can be used for estimating $\dLv$.

Let $\Lind := \Lv - \Ldir$.
Then, the in-scattered radiance $\Lins$~\eqref{eqn:rte3d_Lins} and reflected/refracted radiance $\Lr$~\eqref{eqn:rte3d_re} can be expressed as sums of an \emph{indirect} and a \emph{direct} term of the form
%
\begingroup
\small
\begin{equation}
	\label{eqn:rte3d_Lins_Lr_nee}
	\tcbset{top=0mm,bottom=0mm,left=0mm,right=0mm,fonttitle=\bfseries\scriptsize\color{gray},colbacktitle=white,enhanced,attach boxed title to top center={yshift=-1mm},boxed title style={top=0mm,bottom=0mm,left=0mm,right=0mm}}
	%
	\!\!\!\tcboxmath[colframe=purpleD,colback=white,title=Indirect]{\tint_{\Sph} f(\bx, -\bom', \bom) \,\Lind(\bx, \bom') \,\D\bom'}
	+ \tcboxmath[colframe=greenD,colback=white,title=Direct]{\tint_{\Sph} f(\bx, -\bom', \bom) \,\Ldir(\bx, \bom') \,\D\bom',}
\end{equation}
\endgroup
%
where $f := \phase$ for $\Lins$ and $f := \bsdf$ for $\Lr$.

Since $\dLv = \dLind + \dLdir$, the derivatives of $\Lins$ and $\Lr$ given by Eqs.~\eqref{eqn:rte3d_Lins_diff} and \eqref{eqn:rte3d_Lr_diff} can be rewritten as
%
\begin{equation}
	\label{eqn:rte3d_Lins_diff_nee}
	\tcbset{top=0mm,bottom=0mm,left=0mm,right=0mm,fonttitle=\bfseries\scriptsize\color{gray},colbacktitle=white,enhanced,attach boxed title to top center={yshift=-1.5mm},boxed title style={top=0mm,bottom=0mm,left=0mm,right=0mm}}
	\begin{split}
		\dLins(\bx, \bom)
		=\;& \tcboxmath[colframe=purpleD,colback=white,title=Indirect]{\tint_{\Sph} \sdev\big[\phase(\bx, -\bom', \bom) \,\Lind(\bx, \bom') \big] \,\D\bom'} \;+\\[-2pt]
		& \tcboxmath[colframe=greenD,colback=white,title=Direct]{\tint_{\Sph} \sdev\big[\phase(\bx, -\bom', \bom) \, \Ldir(\bx, \bom')\big] \D\bom'} \;+\\[-2pt]
		& \tcboxmath[colframe=orangeD,colback=white,title=Edge]{\begin{multlined}[c]
				\tint_{\gedge(\bx)} \big\langle \gnormp, \sdev(\by \to \bx) \big\rangle \, \phase(\bx, \bx \to \by, \bom)\\
				\Delta\Lv(\bx, \by \to \bx) \, \vis(\bx, \by) \, \tfrac{\sin\theta}{\| \by - \bx \|} \,\D\ell(\by),
			\end{multlined}}
	\end{split}
\end{equation}
%
and
%
\begin{equation}
	\label{eqn:rte3d_Lr_diff_nee}
	\tcbset{top=0mm,bottom=0mm,left=0mm,right=0mm,fonttitle=\bfseries\scriptsize\color{gray},colbacktitle=white,enhanced,attach boxed title to top center={yshift=-1.5mm},boxed title style={top=0mm,bottom=0mm,left=0mm,right=0mm}}
	\begin{split}
		\dLr(\bx, \bom)
		=\;& \tcboxmath[colframe=purpleD,colback=white,title=Indirect]{\tint_{\Sph} \sdev\big[\bsdf(\bx, -\bom', \bom) \,\Lind(\bx, \bom')\big] \,\D\bom'} \;+\\[-2pt]
		& \tcboxmath[colframe=greenD,colback=white,title=Direct]{\tint_{\Sph} \sdev\big[\bsdf(\bx, -\bom', \bom) \, \Ldir(\bx, \bom')\big] \D\bom'} \;+\\[-2pt]
		& \tcboxmath[colframe=orangeD,colback=white,title=Edge]{\begin{multlined}[c]
				\tint_{\gedge(\bx)} \big\langle \gnormp, \sdev(\by \to \bx) \big\rangle \, f_s(\bx, \bx \to \by, \bom)\\
				\Delta\Lv(\bx, \by \to \bx) \, \vis(\bx, \by) \, \tfrac{\sin\theta}{\| \by - \bx \|} \,\D\ell(\by).
			\end{multlined}}
	\end{split}
\end{equation}
%
In Eqs.~(\ref{eqn:rte3d_Lins_Lr_nee}--\ref{eqn:rte3d_Lr_diff_nee}), the integrals involving $\Ldir$ (enclosed in green boxes) act as the \emph{direct} terms and can be estimated using multiple importance sampling (MIS).
The edge integrals (in orange boxes), on the contrary, depend on the differences $\Delta\Lv$ of $\Lv$ across boundary curves~\eqref{eqn:rte3d_delta_L_def} and need to be estimated by tracing additional \emph{side paths} (see Figure~\ref{fig:side_paths} and discussions later).
Lastly, the integrals involving $\Lind$ (in purple boxes) are analogous to \emph{indirect illumination} and need to be estimated recursively based on
%
\begin{equation}
	\label{eqn:rte3d_Lind}
	\tcbset{top=0mm,bottom=0mm,left=0mm,right=0mm,fonttitle=\bfseries\scriptsize\color{gray},colbacktitle=white,enhanced,attach boxed title to top center={yshift=-1mm},boxed title style={top=0mm,bottom=0mm,left=0mm,right=0mm}}
	\begin{split}
		&\Lind(\bx, \bom)\\
		\!\!\!\!=\;& \tint_0^{\dist} \tcboxmath[colframe=redD,colback=white,title=Volumetric]{\trans(\bx', \bx) \,\sigS(\bx') \,\Lins(\bx', \bom)} \,\D\tau +
		\tcboxmath[colframe=yellowD,colback=white,title=Interfacial]{\trans(\bx_0, \bx)\,\Lr(\bx_0, \bom),}
	\end{split}
\end{equation}
%
and
%
\begin{equation}
	\label{eqn:rte3d_Lind_diff}
	\tcbset{top=0mm,bottom=0mm,left=0mm,right=0mm,fonttitle=\bfseries\scriptsize\color{gray},colbacktitle=white,enhanced,attach boxed title to top center={yshift=-1mm},boxed title style={top=0mm,bottom=0mm,left=0mm,right=0mm}}
	\begin{split}
		& \dLind(\bx, \bom) = \dLv(\bx, \bom) - \dLdir(\bx, \bom)\\[-2pt]
		\!\!\!\!=\;& \int_0^{\dist} \tcboxmath[colframe=redD,colback=white,title=Volumetric]{\begin{multlined}[c]
		 	\trans(\bx', \bx) \big[ \sigS(\bx') \,\dLins(\bx', \bom)\\
			+ \left(\dot{\sigS}(\bx') - \intSigT(\bx, \bom, \tau)\,\sigS(\bx')\right) \Lins(\bx', \bom) \big]
		\end{multlined}} \,\D\tau\\[-2pt]
		& + \tcboxmath[colframe=yellowD,colback=white,title=Interfacial]{\begin{multlined}[t]
			\trans(\bx_0, \bx) \big[ -\left( \intSigT(\bx, \bom, \dist) + \ddist \, \sigT(\bx_0) \right) \Lr(\bx_0, \bom)\\
			+ \dLr(\bx_0, \bom) + \ddist \,\sigS(\bx_0) \,\Lins(\bx_0, \bom) \big].
		\end{multlined}}
	\end{split}
\end{equation}

\input{alg1_v2.tex}

\begin{algorithm}[t]
	\caption{\label{alg:edgeIntegral}
		Monte Carlo estimator of the edge integrals
	}
	\begin{algorithmic}[1]
		\Function{edgeIntegral}{$\bx$, $\dbx$, $\bom$, $f$}
		\State Draw $\by$ from $\boundary(\bx)$ with probability $p(\by)$
		\State $\bom' \gets (\by \to \bx)$
		\State Calculate $\gnormp$, $\dbom'$ using Eqs.~(\ref{eqn:rte3d_sph_norm}, \ref{eqn:rte3d_sph_vel})
		\State $\Delta\Lv \gets \textsc{deltaRadiance}(\bx, \dbx, \bom', \dbom')$
		\label{alg2:line:deltaL} \Comment{Appendix~\ref{sec:mc_ext}}
		\State \Return $\langle \gnormp, \dbom' \rangle \, f(\bx, -\bom', \bom) \, \Delta\Lv \, \frac{\vis(\bx, \by) \,\sin\theta}{\| \by - \bx \| \, p(\by)}$
		\EndFunction
	\end{algorithmic}
\end{algorithm}

\begin{figure}[t]
	\centering
	\includegraphics[height=1.45in]{side_paths_3.pdf}
	\caption{\label{fig:side_paths}
		An \textbf{illustration} of the joint estimation of $\Lind(\bx_0, \bom_0)$ and $\dLind(\bx_0, \bom_0)$ using Algorithm~\ref{alg:diffRadiances} for some $\bx_0 \in \domain \setminus \boundary$ and $\bom_0 \in \Sph$.
		Similar to conventional volumetric path tracing, our algorithm begins with drawing a free-flight distance $\tau$ (Line~\ref{alg1:line:free_distance}), yielding a volumetric event at $\bx_1 = \bx_0 - \tau \bom_0$ in this example.
		Then, a new direction $\bom_1$ is drawn (Line~\ref{alg1:line:vol_newdir}) followed by a recursive estimation of the indirect terms $\Lind(\bx_1, \bom_1)$ and $\dLind(\bx_1, \bom_1)$ (Line~\ref{alg1:line:vol_recurse}).
		Next, the direct terms are estimated at $\bx_1$ (Lines~\ref{alg1:line:vol_Lins_1} and \ref{alg1:line:vol_dLins_2}).
		Lastly, the edge integral term $\Bins(\bx_1, \bom_0)$ is handled by drawing another direction $\bom'_1$ in which $\Lv(\bx_1, \bom'_1)$ is discontinuous (Line~\ref{alg1:line:vol_dLins_3}).
		In this example, $\Delta\Lv(\bx_1, \bom'_1)$ arises from a sharp edge and can be estimated by tracing one side path via Eq.~\eqref{eqn:rte3d_deltaL_sharp}.
		Similarly, the joint estimation of $\Lind(\bx_1, \bom_1)$ and $\dLind(\bx_1, \bom_1)$ boils down to (i)~calculating the indirect terms $\Lind(\bx_2, \bom_2)$ and $\dLind(\bx_2, \bom_2)$ recursively; and (ii)~estimating the direct and edge terms.
		In this illustration, $\Delta\Lv(\bx_2, \bom'_2)$ results from a silhouette edge and is estimated by tracing two side paths via Eq.~\eqref{eqn:rte3d_deltaL_silhouette}.
		When estimating $\Lind(\bx_2, \bom_2)$ and $\dLind(\bx_2, \bom_2)$, $\bx_3 \in \boundary$ corresponds to an interfacial event where the in-scattered radiance~$\Lins$ is estimated using Eq.~\eqref{eqn:rte3d_Lins_boundary}.
	}
\end{figure}

\subsubsection{Our Estimator.}
%
Based on Eqs.~(\ref{eqn:rte3d_Lins_Lr_nee}--\ref{eqn:rte3d_Lind_diff}), we introduce an unbiased approach (Algorithm~\ref{alg:diffRadiances}) that jointly estimates $\Lind$ and $\dLind$ using Monte Carlo integration.

\paragraph{Sampling free-flight distance}
Similar to VPT, our algorithm starts with drawing the free-flight distance~$\tau$ (Line~\ref{alg1:line:free_distance} of Algorithm~\ref{alg:diffRadiances}).
If $\tau < \dist$, the \emph{volumetric} terms (in pink boxes) from Eqs.~\eqref{eqn:rte3d_Lind} and \eqref{eqn:rte3d_Lind_diff} are estimated at $\bx' = \bx - \tau\bom$;
otherwise, the \emph{interfacial} terms (in yellow boxes) are calculated at $\bx_0 = \bx - \dist\bom$.
We detail how these terms are estimated in the following.

\paragraph{Volumetric}
The volumetric terms in Eqs.~\eqref{eqn:rte3d_Lind} and \eqref{eqn:rte3d_Lind_diff} are estimated jointly by Lines~\ref{alg1:line:vol_start}--\ref{alg1:line:vol_end} of Algorithm~\ref{alg:diffRadiances}.
When the volumetric terms are picked, $\tau < \dist$, and the probability density of $\tau$ is $p(\tau) = \sigT(\bx')\,\trans(\bx', \bx)$.
Dividing $p(\tau)$ from the integrand cancels out the transmittance terms, yielding
%
\begin{equation}
	\label{eqn:Lind_est_vol}
	\langle \Lind \rangle = \alpha\,\Lins, \quad
	\langle \dLind \rangle = \alpha\,\dLins + \beta\,\Lins,
\end{equation}
%
with $\alpha$ and $\beta$ given by Lines~\ref{alg1:line:vol_alpha} and \ref{alg1:line:vol_beta}.
To estimate $\Lins(\bx', \bom)$ and $\dLins(\bx', \bom)$, the main ingredients of Eq.~\eqref{eqn:Lind_est_vol}, we apply Monte Carlo integration to Eqs.~\eqref{eqn:rte3d_Lins_Lr_nee} and \eqref{eqn:rte3d_Lins_diff_nee} jointly as follows.
%
\begin{myitemizeP}
	\item We estimate the indirect terms by drawing a random $\bom'$ according to the phase function~$\phase$ and obtain $\Lind(\bx', \bom')$ and $\dLind(\bx', \bom')$ recursively (Line~\ref{alg1:line:vol_recurse} of Algorithm~\ref{alg:diffRadiances}). Techniques like Russian roulette can be used here to avoid infinite recursion.
	%
	\item Then, the indirect terms of Eqs.~\eqref{eqn:rte3d_Lins_Lr_nee} and \eqref{eqn:rte3d_Lins_diff_nee} are respectively given by $\Lind(\bx', \bom')$ and
	%
	\[	\begin{split}
			& \sdev\big[\phase(\bx', -\bom', \bom) \,\Lind(\bx, \bom')\big]/\phase(\bx', -\bom', \bom)\\
			=\;& \dLind(\bx, \bom') + \tfrac{\dphase(\bx', -\bom', \bom)}{\phase(\bx', -\bom', \bom)} \Lind(\bx, \bom'),
	\end{split}	\]
	%
	which corresponds to Lines~\ref{alg1:line:vol_Lins_0} and \ref{alg1:line:vol_dLins_1}.
\end{myitemizeP}
%
\begin{myitemizeG}
	\item The direct term can be estimated directly using MIS (Lines~\ref{alg1:line:vol_Lins_1} and \ref{alg1:line:vol_dLins_2}).
\end{myitemizeG}
%
\begin{myitemizeO}
	\item The edge-integral term is computed via the \textsc{EdgeIntegral} function (Line~\ref{alg1:line:vol_dLins_3}) that uniformly samples a point~$\by$ from an edge and calculates the corresponding $\Delta\Lv$ (Algorithm~\ref{alg:edgeIntegral}).
\end{myitemizeO}

\paragraph{Interfacial}
The interfacial terms in Eqs.~\eqref{eqn:rte3d_Lind} and \eqref{eqn:rte3d_Lind_diff} are handled with Lines~\ref{alg1:line:surf_start}--\ref{alg1:line:surf_end} of Algorithm~\ref{alg:diffRadiances}.
Since $p(\tau) = \trans(\bx_0, \bx)$, dividing it from the direct terms cancels out the transmittances, yielding
%
\begin{equation}
	\label{eqn:Lind_est_surf}
	\langle \Lind \rangle = \Lr, \quad
	\langle \dLind \rangle = -\beta\,\Lr + \dLr + \ddist \,\sigS \,\Lins,
\end{equation}
%
with $\beta$ calculated in Line~\ref{alg1:line:surf_beta}.
We estimate $\Lr(\bx_0, \bom)$ and $\dLr(\bx_0, \bom)$, the main ingredients of Eq.~\eqref{eqn:Lind_est_surf}, based on Eqs.~\eqref{eqn:rte3d_Lins_Lr_nee} and \eqref{eqn:rte3d_Lr_diff_nee} as follows.
%
\begin{myitemizeP}
	\item The indirect terms are computed jointly by sampling $\bom'$ according to the local BSDF and acquiring $\Lv(\bx_0, \bom')$ and $\dLv(\bx_0, \bom')$ recursively (Lines~\ref{alg1:line:surf_recurse}--\ref{alg1:line:surf_dLr_1} of Algorithm~\ref{alg:diffRadiances}).
\end{myitemizeP}
%
\begin{myitemizeG}
	\item The direct terms are estimated directly using MIS (Lines~\ref{alg1:line:surf_Lr_1} and \ref{alg1:line:surf_dLr_2}).
\end{myitemizeG}
%
\begin{myitemizeO}
	\item Similar to the volumetric case, the edge-integral term is estimated using the \textsc{EdgeIntegral} function (Line~\ref{alg1:line:surf_dLr_3}).
	%
	\item Additionally, the in-scattered radiance $\Lins(\bx_0, \bom)$, which arises from the differentiation of $\kernelT\,\kernelS\,\Lv$ shown in Eq.~\eqref{eqn:rte3d_ind_diff}, needs to be calculated (Line~\ref{alg1:line:surf_Lins}) according to Eq.~\eqref{eqn:rte3d_Lins_boundary}.
\end{myitemizeO}

The background colors of individual lines in Algorithm~\ref{alg:diffRadiances} indicate the correspondence to the indirect (purple), direct (green), and edge (orange) terms in Eqs.~(\ref{eqn:rte3d_Lins_Lr_nee}--\ref{eqn:rte3d_Lr_diff_nee}).

Lastly, calculating $\Lins(\bx_0, \bom)$ (Line~\ref{alg1:line:surf_Lins} of Algorithm~\ref{alg:diffRadiances}) and $\Delta\Lv$ (Line~\ref{alg2:line:deltaL} of Algorithm~\ref{alg:edgeIntegral}) requires tracing additional \emph{side paths}.
Since both $\Lins(\bx_0, \bom)$ and $\Delta\Lv$ depend on only the interior radiance $\Lv$ but not its derivative $\dLv$, the side paths traced to estimate these terms can be handled using conventional volumetric path tracing.
For more details, please refer to Appendix~\ref{sec:mc_ext}.

We illustrate the execution of Algorithm~\ref{alg:diffRadiances} in Figure~\ref{fig:side_paths}.

\input{pixel.tex}

\subsection{Discussion}
\label{ssec:mc_discussion}
%
\paragraph{Estimating total derivatives}
To estimate the total scene derivative~$\dev\Lv$, Algorithm~\ref{alg:diffRadiances} only needs to be executed once in a \emph{vectorized} fashion.
Specifically, one can replace all partial derivatives with their (vector-valued) total derivative counterparts and perform component-wise calculations between all total-derivative-related quantities.
In this way, all paths traced by Algorithm~\ref{alg:diffRadiances} are effectively shared to calculate derivatives with respect to all scene parameters.

\paragraph{Homogeneous media}
%
Our derivations in \S\ref{sec:diff_rte} and \S\ref{ssec:our_est} enjoy the generality for handling heterogeneous media with spatially varying radiative transfer parameters $\sigT$, $\sigS$, and $\phase$.

On the other hand, for homogeneous media with $\param$-invariant parameters, $\nabla\sigS = \nabla\sigT \equiv \bm{0}$, causing (i)~$\beta \equiv 0$ for Lines~\ref{alg1:line:vol_beta} and \ref{alg1:line:surf_beta} of Algorithm~\ref{alg:diffRadiances}; and (ii)~$\dLind$ and $\dLdir$ in Eqs.~\eqref{eqn:rte3d_Ldir_diff} and \eqref{eqn:rte3d_Lind_diff} to become:
%
\begin{gather}
	\label{eqn:rte3d_Ldir_diff_homo}
	\dLdir(\bx, \bom) = \E^{-\sigT\dist} \left[ \dLe(\bx_0, \bom) - \ddist \,\sigT \,\Le(\bx_0, \bom) \right],
	%
	\\[5pt]
	%
	\label{eqn:rte3d_rte_diff_homo}
	\begin{split}
	& \dLind(\bx, \bom) = \tint_0^{\dist} \E^{-\sigT \tau} \sigS \,\dLins(\bx', \bom) \,\D\tau \; +\\
	& \quad \E^{-\sigT\dist} \left[ \dLr(\bx_0, \bom) + \ddist \left( \sigS \,\Lins(\bx_0, \bom) - \sigT \,\Lr(\bx_0, \bom) \right) \right].
	\end{split}
\end{gather}

\paragraph{Computational efficiency}
Since our estimator is built upon ordinary Monte Carlo (OMC), it has a theoretical convergence rate of $O(1/\sqrt{N})$ with $N$ being the number of (main) sample paths.
In practice, because of the additional side paths needed to be traced in each iteration for estimating $\Delta\Lv$ and $\Lins$, Algorithm~\ref{alg:diffRadiances} is significantly more expensive than conventional VPT by having a time complexity of $O(K^2)$, as opposed to $O(K)$, for each sample path of length $K$.
Designing new Monte Carlo estimators of scene derivatives with better performance is an exciting future research direction.
