\section{Details on Our Monte Carlo Estimator}
\label{sec:mc_ext}
%
We now provide more details on the estimation of radiance change~$\Delta\Lv$ and in-scattered radiance $\Lins$.

\paragraph{Estimating $\Delta\Lv$}
Computing the edge integrals (Lines~\ref{alg1:line:vol_dLins_3} and \ref{alg1:line:surf_dLr_3} of Algorithm~\ref{alg:diffRadiances})  requires estimating the change of radiance $\Delta\Lv$ across discontinuity boundaries (Line~\ref{alg2:line:deltaL} of Algorithm~\ref{alg:edgeIntegral}).

Recall that the discontinuity of $\Lv$, under the assumptions stated at the beginning of \S\ref{sec:diff_rte}, arises due to the geometric discontinuities resulting from boundary, silhouette, and sharp edges (Figure~\ref{fig:edge_types}).
Based on this observation, Eqs.~\eqref{eqn:rte3d_deltaL_silhouette} and \eqref{eqn:rte3d_deltaL_sharp} derived in Appendix~\ref{ssec:deltaL_deriv} evaluate $\Delta\Lv$ under two cases.
These equations can further be estimated using Monte Carlo integration, as shown in Algorithm~\ref{alg:deltaL}.

Discontinuities caused by boundary and silhouette edges~\eqref{eqn:rte3d_deltaL_silhouette} are handled by Lines~\ref{alg3:line:case1_start}--\ref{alg3:line:case1_end} of Algorithm~\ref{alg:deltaL}.
Let $\dist^\mathrm{near} := \min(\dist^-, \dist^+)$ and $\dist^\mathrm{far} := \max(\dist^-, \dist^+)$.
Then, due to Eq.~\eqref{eqn:rte3d_deltaL_silhouette}, it holds that
%
\begin{equation}
	\label{eqn:rte3d_deltaL_silhouette_est}
	\begin{split}	
		& |\Delta\Lv(\bx, \bom)|\\
		\!\!\!\!=\;& \tint_{\dist^\mathrm{near}}^{\dist^\mathrm{far}} \trans(\bx', \bx) \,\sigS(\bx') \, \Lins(\bx', \bom) \,\D\tau
		+ \trans(\bx^\mathrm{far}, \bx) \,\Ls(\bx^\mathrm{far}, \bom)\\
		&  - \trans(\bx^\mathrm{near}, \bx) \,\Ls(\bx^\mathrm{near}, \bom),
	\end{split}
\end{equation}
%
where $\bx^\mathrm{near} := \bx - \dist^\mathrm{near}\bom$ and $\bx^\mathrm{far} := \bx - \dist^\mathrm{far}\bom$.
After obtaining $|\Delta\Lv|$, it is easy to verify that $\Delta\Lv$ equals $|\Delta\Lv|$ if $\dist^- > \dist^+$ and $-|\Delta\Lv|$ otherwise.
Further, Eq.~\eqref{eqn:rte3d_deltaL_silhouette_est} can be estimated by tracing two side paths: one for the first two terms on the RHS side (Lines~\ref{alg3:line:sidepath_2_start}--\ref{alg3:line:sidepath_2_end}) and the other for the last term (Line~\ref{alg3:line:sidepath_1}).

Lastly, discontinuities caused by sharp edges~\eqref{eqn:rte3d_deltaL_sharp} can be handled by tracing only one side path, as shown in Lines~\ref{alg3:line:case1_start}--\ref{alg3:line:case1_end} of Algorithm~\ref{alg:deltaL}.

\paragraph{Estimating $\Lins$}
To estimate $\Lins(\bx, \bom)$ given by Eq.~\eqref{eqn:rte3d_Lins} for interior points $\bx \in \domain\setminus\boundary$ and Eq.~\eqref{eqn:rte3d_Lins_boundary} for boundary ones $\bx \in \boundary$, we randomly draw a direction $\bom'$ according to the phase function at $\bx$ and estimate $\Lv(\bx, \bom')$ or $\Ls(\bx, \bom')$ using standard volumetric path tracing.
Algorithm~\ref{alg:inScatteredRadiance} contains the pseudocode for this process.

\input{alg2.tex}
\input{alg3.tex}
