\section{Preliminaries}
\label{sec:background}
%
We start by briefly reviewing the radiative transfer framework~(\S\ref{ssec:rte}).
We then define the radiance derivative with respect to arbitrary scene parameterizations~(\S\ref{ssec:shape_dev}). Finally, we introduce mathematical background that paves the way for developing our theory of differential radiative transfer (\S\ref{ssec:diff_integral}).
Table~\ref{tab:los} summarizes our notation.

\subsection{Radiative Transfer}
\label{ssec:rte}
%
\emph{Radiative transfer}~\cite{chandrasekhar1960radiative} uses energy conservation principles to model light transport in participating media.
%
At its core is the \emph{radiative transfer equation}~(RTE): 
Consider a medium confined in a volume $\domain \subseteq \Real^3$ with boundary $\boundary$. 
The \emph{interior radiance field}~$\Lv(\bx,\bom)$ is defined over positions $\bx \in \domain \setminus \boundary$ and light-transport directions $\bom\in\Sph$.
The \emph{steady-state} RTE is a linear integral equation on $\Lv$ that can be expressed in operator form as~\cite{case1967linear,spanier1969mc}:
%
\begin{equation}
	\label{eqn:rte3d}
	\Lv = \kernelT\,\kernelC\,\Lv + \src.
\end{equation}
%
Concretely, the \emph{transport operator} $\kernelT$ maps any function $g:(\domain \setminus \boundary) \times \Sph \to \Real_+$ to a new function
%
\begin{equation}
	\label{eqn:rte3d_kernelT}
	(\kernelT\,g)(\bx, \bom) = \tint_0^{\dist} \trans(\bx', \bx) \, g(\bx', \bom) \,\D\tau,
\end{equation}
%
where: $\bx' := \bx - \tau\bom$; $\dist$ is the distance from $\bx$ to the medium's boundary in the direction of $-\bom$,
%
\begin{equation}
	\label{eqn:rte3d_dist}
	\dist = \inf\{ \tau \in \Real_+ \;:\; \bx - \tau\bom \in \boundary \};
\end{equation}
%
and $\trans(\bx', \bx)$ is the \emph{transmittance} between $\bx'$ and 
$\bx$, that is, the fraction of light that transmits straight between $\bx'$ and $\bx$ without being absorbed or scattered away:
%
\begin{equation}
	\label{eqn:rte3d_trans}
	\trans(\bx', \bx) = \exp\left( -\tint_0^{\tau} \sigT(\bx - \tau'\bom) \,\D\tau' \right),
\end{equation}
%
with $\sigT$ denoting the medium's \emph{extinction coefficient}.

The \emph{collision operator} $\kernelC$ maps the interior radiance field  $\Lv$ to%
\footnote{In this paper, we follow the convention that all directions point away from $\bx$ when expressing BSDFs and phase functions, yielding the negative sign before $\bom'$ in Eq.~\eqref{eqn:rte3d_Lins}.}
%
\begin{equation}
	\label{eqn:rte3d_Lins}
	(\kernelC\,\Lv)(\bx, \bom) = \sigS(\bx) \,\underbrace{\tint_{\Sph} \phase(\bx, -\bom', \bom) \, \Lv(\bx, \bom') \, \D\bom'}_{=:\;\Lins(\bx, \bom)},
\end{equation}
%
where $\sigS$ and $\phase$ denote the medium's \emph{scattering coefficient} and \emph{single-scattering phase function}, respectively. $\Lins$ is usually termed as the \emph{in-scattered} radiance.%\igkiou{Doesn't in-scattered include the scattering coefficient?}
\footnote{Precisely, $\Lins$ captures the in-scattered radiance before volumetric absorption.}

The \emph{source term} $\src$ of the RTE~\eqref{eqn:rte3d} is sometimes assumed given for problems with simple geometries (hence the name).
Generally, with the presence of nontrivial interfaces, $\src$ equals%
\footnote{
	In Eq.~\eqref{eqn:rte3d_src}, the absorption coefficient~$\sigA$ and radiant emission~$\Lve$ are sometimes combined into a single volumetric source term.
	Here, we follow the work by Kutz~et~al.~\shortcite{Kutz:2017:SDT} and separate the two terms.
	Additionally, when $\domain$ is unbounded, $\dist$ can be infinite for certain $\bx$ and $\bom$, causing the second term on the RHS to vanish.
}
%
\begin{equation}
	\label{eqn:rte3d_src}
	\!\!\!\!\src(\bx, \bom) = \underbrace{\tint_0^{\dist} T(\bx', \bx)\,\sigA(\bx')\,\Lve(\bx', \bom) \,\D\tau}_{=\;(\kernelT\,\sigA\,\Lve)(\bx, \bom)}
	+\; T(\bx_0, \bx)\,\Ls(\bx_0, \bom),
\end{equation}
%
where $\bx_0 := \bx - \dist\bom$ is a point on the medium's boundary, and $\sigA := \sigT - \sigS$ is the \emph{absorption coefficient}.
On the right-hand side~(RHS) of this equation, $\Lve$ represents the medium's \emph{radiant emission}, and $\Ls$ indicates the \emph{interfacial radiance} governed by the \emph{rendering equation}~(RE) for all $\bx \in \boundary$ and $\bom \in \Sph$:
%
\begin{equation}
	\label{eqn:rte3d_re}
	\Ls(\bx, \bom) = \underbrace{\tint_{\Sph} \bsdf(\bx, -\bom', \bom) \, \Lv(\bx, \bom') \,\D\bom'}_{=:\;\Lr(\bx, \bom)} + \,\Le(\bx, \bom),
\end{equation}
%
where $\bsdf$ is the \emph{cosine-weighted BSDF}, $\Le$ denotes the \emph{interfacially emitted radiance}, and $\Lr$ indicates the \emph{reflected/refracted radiance}.

\paragraph{Combining the RTE and the RE}
The source term $Q$ of the RTE~\eqref{eqn:rte3d} can also be expressed compactly 
using operators. 
Substituting Eq.~\eqref{eqn:rte3d_re} into Eq.~\eqref{eqn:rte3d_src} 
and regrouping the resulting terms, we have
%
\begin{equation}
	\label{eqn:rte3d_src1}
	\src = \kernelS\,\Lv + \Ldir,
\end{equation}
%
where the \emph{interfacial scattering operator} $\kernelS$ is defined as
%
\begin{equation}
	\label{eqn:rte3d_kernelS}
	\begin{split}
		(\kernelS\,\Lv)(\bx, \bom) &= T(\bx_0, \bx)\,\Lr(\bx_0, \bom)\\
		&= T(\bx_0, \bx) \tint_{\Sph} \bsdf(\bx_0, -\bom', \bom) \, \Lv(\bx_0, \bom') \,\D\bom',
	\end{split}
\end{equation}
%
which accounts for contributed radiance due to reflection and transmission at the medium's boundary.
The second term in Eq.~\eqref{eqn:rte3d_src1} is
\begin{equation}
	\label{eqn:rte3d_Ldir}
	\Ldir(\bx, \bom) := (\kernelT\,\sigA\,\Lve)(\bx, \bom) + T(\bx_0, \bx) \,\Le(\bx_0, \bom),
\end{equation}
%
indicating the contribution of radiant emission in the medium and from its boundary.

Based on Eqs.~(\ref{eqn:rte3d_src1}--\ref{eqn:rte3d_Ldir}), the
RTE~\eqref{eqn:rte3d} and the RE~\eqref{eqn:rte3d_re} can be combined into an integral
equation on $\Lv$ with an operator form:
%
\begin{myBox}
	\vspace{-1mm}
	\begin{equation}
		\label{eqn:rte3d_Lonly}
		\Lv = (\kernelT\,\kernelC + \kernelS)\,\Lv + \Ldir.
	\end{equation}
\end{myBox}
%
\noindent
Our core theoretical contribution is to derive expressions for the derivatives of the interior radiance $\Lv$ by differentiating individual terms in this equation.

\subsection{Scene Derivatives}
\label{ssec:shape_dev}
%
We assume that the scene specifications (such as object geometries and material properties) are parameterized by a set of $m$ parameters $\bparam = \{\param_1, \param_2, \ldots, \param_m\}$, and each parameter varies continuously.  
The medium's domain and its boundary may also vary with these parameters.
We therefore denote them as $\domain(\bparam)$ and $\boundary(\bparam)$, respectively.

To give a few examples, the translation of an object in 3D along a direction $\hat{\bom}$ can be parameterized by $\by(\param) = \by + \param\hat{\bom}$ for each point $\by$ of the object.
The radiance field $\Lv$ depends on the object position and subsequently the parameter $\param$.
Further, any linear transformed version of the object can be expressed as $\by(\param) = M(\param)\,\by$, where $M \in \Real^{3 \times 3}$ is a transformation matrix of parameter $\param$.

\paragraph{Our problem}
Provided a specific scene parameterization, we aim to develop a theory for computing the total derivative (or gradient) of the interior radiance $\Lv$ with respect to $\bparam$. We denote the total derivative as $\dev\Lv \in \Real^m$ and refer it as the \emph{scene derivative} of $\Lv$.

Let $\Lv(\bx, \bom; \bparam)$ be the interior radiance $\Lv$ at $\bx \in \domain(\bparam) \setminus \boundary(\bparam)$ with the direction $\bom \in \Sph$.
The scene derivative, $\dev\Lv(\bx(\bparam), \bom(\bparam); \bparam)$, is an $m\times1$ vector whose $i$-th element is the partial derivative of $\Lv$ with respect to the $i$-th parameter $\param_i$.
That is,
%
\begin{multline}
	\label{eqn:rte3d_prob}
	\partial_{\param_i} \Lv(\bx(\bparam), \bom(\bparam); \bparam)\\
	= \lim_{\epsilon \to 0} \frac{ \Lv(\bx(\bparam'_i), \bom(\bparam'_i); \bparam'_i) - \Lv(\bx(\bparam), \bom(\bparam); \bparam)}{\epsilon},
\end{multline}
%
where $\bparam'_i := \{\param_1, \ldots, \param_{i - 1}, \param_i + \epsilon, \param_{i + 1}, \ldots, \param_m\}$.
We note that it is crucial to allow both $\bx$ and $\bom$ to depend on the scene parameters~$\bparam$.
As we will show later, this dependence emerges in a number of interesting applications such as (i)~interfacial radiances $\Ls$ on evolving surfaces and (ii)~radiometric measurements from moving sensors.

Without loss of generality, our derivation focuses on the partial derivative of $\Lv$, denoted as $\dLv := \sdev \Lv$, with respect to one scene parameter $\param \in \bparam$.
For notational convenience, we omit the dependencies on $\bparam$ from $\bx$, $\bom$, and the radiances in the rest of the paper.

\subsection{Differentiating Integrals}
\label{ssec:diff_integral}
%
Deriving analytic expressions for the scene derivatives of radiances requires differentiating the integral operators in the RTE~\eqref{eqn:rte3d_Lonly} and the RE~\eqref{eqn:rte3d_re}.
To this end, we use the Reynolds transport theorem~\cite{leal2007advanced},%
\footnote{We present as \thmref{reynolds} a slightly generalized version of the original theorem.}
which originated in fluid mechanics and is a generalization of Leibniz's integral rule for differentiation~\cite{flanders1973differentiation}.

\begin{thm}{Reynolds transport theorem}{reynolds}
	Let $f$ be a scalar-valued function defined on some $n$-dimensional manifold~$\Omega(\param)$ parameterized with some $\param \in \Real$.
	Additionally, let $\interface(\param) \subset \domain(\param)$ be an $(n - 1)$-dimensional manifold given by the union of the \emph{external} boundary~$\boundary(\param)$ and the \emph{internal} one containing the discontinuous locations of $f$.
	Then, it holds that
	%
	\begin{equation}
		\label{eqn:thm1}
		\sdev \left( \int_{\domain(\param)} f \,\D\domain(\param) \right)
		= \tint_{\domain(\param)} \dot{f} \,\D\domain(\param)
		+ \tint_{\interface(\param)} \langle \gnorm, \dbx \rangle \,\Delta f \,\D\interface(\param),
	\end{equation}
	%
	where $\dot{f} := \sdev f$, $\dbx := \sdev\bx$, $\D\domain$ and $\D\interface$ respectively denote the standard measures associated with $\domain$ and $\interface$; and $\langle \cdot, \cdot \rangle$ indicates the dot (inner) product between two vectors.
	Further, $\gnorm$ is the normal direction at each $\bx \in \interface(\param)$, and $\Delta f$ is given by%
	\footnote{
		When $\bx$ approaches $\interface(\param)$ from the exterior of the integral domain~$\domain(\param)$, the corresponding one-sided limit in Eq.~\eqref{eqn:thm1_limits} is set to zero.
	}
	%
	\begin{equation}
		\label{eqn:thm1_limits}
		\Delta f(\bx) := \lim_{\epsilon \to 0^-} f(\bx + \epsilon\,\gnorm)
		- \lim_{\epsilon \to 0^+} f(\bx + \epsilon\,\gnorm),
	\end{equation}
	%
	for all $\bx \in \interface(\param)$.
\end{thm}

Intuitively, \thmref{reynolds} states that the derivative of the integral of a continuous function $f$ over some domain $\domain(\param)$ depends on (i)~the derivative of $f$ itself over $\domain(\param)$ (the first integral); and (ii)~the ``normal velocity'' of the domain boundary and internal discontinuities~$\interface(\param)$ modulated by the differences between function values across these boundaries (the second integral)%
\footnote{
	Previously, Li~et~al.~\shortcite{Li:2018:DMC} used a different but mathematically equivalent method to demonstrate the need of the boundary integral.
	This equivalence was further discussed in his Ph.D. thesis~\shortcite{li2019thesis}.
}.

In a special case, for classical Riemann integrals where $\domain$ is an interval $(a, b) \subset \Real$ and the integrand~$f$ is differentiable everywhere, it holds that $\interface = \{a, b\}$, and the interface integral in Eq.~\eqref{eqn:thm1} reduces to the sum of the integrand evaluated at $a$ and $b$.
Let $\gnorm = 1$, then $\langle \gnorm, \dbx \rangle$ equals $\dot{a}$ at $x = a$ and $\dot{b}$ at $x = b$.
Further, it is easy to verify that $\Delta f(b(\param), \param) = f(b(\param), \param)$ and $\Delta f(a(\param), \param) = -f(a(\param), \param)$, yielding
%
\begin{equation}
	\label{eqn:thm1_red}
	\begin{split}
		\sdev\tint_{a(\param)}^{b(\param)} f(x, \param) \,\D x =\;& \tint_{a(\param)}^{b(\param)} \dot{f}(x, \param) \,\D x\\
		& + \dot{b}(\param) \, f(b(\param), \param) - \dot{a}(\param) \, f(a(\param), \param),
	\end{split}
\end{equation}
%
which is precisely Leibniz's rule for differentiation.

In the rest of this paper, we detail how \thmref{reynolds} can be used to calculate scene derivatives of radiances (\S\ref{sec:diff_rte}), and introduce an unbiased Monte Carlo method for estimating those derivatives (\S\ref{sec:mc}).
