\begin{algorithm}[t]
	\caption{\label{alg:diffRadiances}
		Unbiased Monte Carlo estimator of $\dLv$
	}
	\begin{algorithmic}[1]
		\Function{$\dLv$}{$\bx$, $\dbx$, $\bom$, $\dbom$}
			\State Calculate $\dist$ and $\ddist$ based on $\bx$, $\dbx$, $\bom$, and $\dbom$
				\Comment{Appendix~\ref{sec:devD}}
			\State $\dLv \gets \dLdir(\bx, \bom)$
				\Comment{\eqref{eqn:rte3d_Ldir_diff}}
			\State $\throughput \gets 1$
				\Comment{path throughput}
			\While{true}
				\State $\tau \gets -\log(\mathrm{rand}())/\sigT$ \label{alg1:line:free_distance}
				\If{$\tau < \dist$} \Comment{volumetric event}
					\State $\bx \gets \bx - \tau\bom$ \label{alg1:line:vol_start}
						\Comment{$\bx \in \domain \setminus \boundary$}
					\State $\dbx \gets \dbx - \tau\dbom$
					\State $\throughput \gets \throughput (\sigS/\sigT)$
					\State Estimate $v_0 = \int_{\Sph} \dphase(-\bom', \bom) \,\Lind(\bx, \bom') \,\D\bom$
						\label{alg1:line:vol_v0}
					\State Estimate $v_1 = \int_{\Sph} \devtext[\phase(-\bom', \bom) \,\Ldir(\bx, \bom')] \,\D\bom$
						\label{alg1:line:vol_v1}
					\State Draw $\by$ from $\boundary(\bx)$ with probability $p(\by)$
						\label{alg1:line:vol_v2_start}
					\State $\bom' \gets (\by \to \bx)$
					\State Calculate $\gnorm$, $\dbom'$ using Eqs.~(\ref{eqn:rte3d_sph_norm}, \ref{eqn:rte3d_sph_vel})
					\State $\Delta\Lv \gets \Delta\Lv(\bx, \dbx, \bom', \dbom')$ \Comment{Alg.~\ref{alg:deltaRadiance}}
						\label{alg1:line:vol_deltaL}
					\State $v_2 \gets | \langle \gnorm, \dbom' \rangle | \, \phase(-\bom', \bom) \, \Delta\Lv \, \frac{\vis(\bx, \by) \,\sin\theta}{\| \by - \bx \| \, p(\by)}$
						\label{alg1:line:vol_v2_end}
					\State $\dLv \gets \dLv + (v_0 + v_1 + v_2) \, \throughput$
					\State Draw $\bom' \sim \phase(-\bom', \bom)$ \label{alg1:line:vol_recurse}
					\State $\bom \gets \bom'$, $\dbom \gets \bm{0}$ \label{alg1:line:vol_end}
				\Else \Comment{interfacial event}
					\State $\bx \gets \bx - \dist\bom$
						\Comment{$\bx \in \boundary$}
					\State $\dbx \gets \dbx - \ddist\bom - \dist\dbom$ \label{alg1:line:surf_start}
					\State Estimate $v_0 = \int_{\Sph} \dbsdf(\bx, -\bom', \bom) \,\Lind(\bx, \bom') \,\D\bom$
						\label{alg1:line:surf_dLr_start}
					\State Estimate $v_1 = \int_{\Sph} \devtext[\bsdf(\bx, -\bom', \bom) \,\Ldir(\bx, \bom')] \,\D\bom$		
						\label{alg1:line:surf_v1}
					\State Draw $\by$ from $\boundary(\bx)$ with probability $p(\by)$
					\State $\bom' \gets (\by \to \bx)$
					\State Calculate $\gnorm$, $\dbom'$ using Eqs.~(\ref{eqn:rte3d_sph_norm}, \ref{eqn:rte3d_sph_vel})
					\State $\Delta\Lv \gets \Delta\Lv(\bx, \dbx, \bom', \dbom')$ \Comment{Alg.~\ref{alg:deltaRadiance}}
						\label{alg1:line:surf_deltaL}
					\State $v_2 \gets | \langle \gnorm, \dbom' \rangle | \, \bsdf(\bx, -\bom', \bom) \, \Delta\Lv \, \frac{\vis(\bx, \by) \,\sin\theta}{\| \by - \bx \| \, p(\by)}$
						\label{alg1:line:surf_dLr_end}
					\State Estimate $\Lins = \Lins(\bx, \bom)$ \Comment{Alg.~\ref{alg:inScatteredRadiance}}
						\label{alg1:line:surf_Lins}
					\State Estimate $\Lr = \Lr(\bx, \bom)$
						\label{alg1:line:surf_Lr}
					\State $\dLv \gets \dLv + [v_0 + v_1 + v_2 + \ddist \, (\sigS \,\Lins - \sigT \,\Lr)] \,\throughput$
					\State Draw $\bom'$ with probability $p(\bom')$
					\State $\throughput \gets \throughput \,(\bsdf(\bx, -\bom', \bom)/p(\bom'))$
					\State $\bom \gets \bom'$, $\dbom \gets \bm{0}$ \label{alg1:line:surf_end}
				\EndIf
				\State Calculate $\dist$ and $\ddist$ based on $\bx$, $\dbx$, $\bom$, and $\dbom$
			\EndWhile
			\State \Return $\dLv$
		\EndFunction
	\end{algorithmic}
\end{algorithm}
