# A Review of Image Reconstruction Algorithms for Diffuse Optical Tomography
## 3 Forward Process  
### Radiative Transfer Equation  
Light propagation through biological tissues is mathematically described by the radiative transfer equation (RTE), as follows [23]:  
\[\left\{\frac{1}{c}\frac{\partial}{\partial t}+\mathbf{s}\cdot\
abla+\mu_{a}( \mathbf{r})+\mu_{s}(\mathbf{r})\right\}I(\mathbf{r},\mathbf{s},t)=\mu_{s}( \mathbf{r})\int_{4\pi}p(\mathbf{s}\cdot\mathbf{s}^{\prime})I(\mathbf{r}, \mathbf{s}^{\prime},t)d\Omega^{\prime}+q(\mathbf{r},\mathbf{s},t), \tag{10}\]  
where \(c\) is the speed of light, \(I(\mathbf{r}\)**,\(\mathbf{s}\),\(t)\)** the radiance (light intensity), and \(q(\mathbf{r}\)**,\(\mathbf{s}\),\(t)\)** the light source term. \(\mathbf{r}\) and \(\mathbf{s}\) represent the position and direction of light, respectively. \(p(\mathbf{s}\cdot\mathbf{s}^{\prime})\) is the scattering phase function, which represents the probability that the direction of light \(\mathbf{s}^{\prime}\) is changed to \(\mathbf{s}\) by the scattering event. The Henyey-Greenstein function is commonly used as \(p(\mathbf{s}\)**,\(\mathbf{s}^{\prime})\)** in DOT image reconstruction. Equation (10) represents the time-dependent case. The time-derivative term on the left side of Equation (10) is removed in the case of time-independent CW measurement. Fourier-transformed versions of Equation (10) are used for FD measurements. The computation of the RTE does not impose any limitations on the optical properties. Therefore, RTE can compute light propagation more accurately than other methods, such as DA, and can be the best choice for image reconstruction for objects that include void-like regions, such as the trachea in the neck, which can affect DOT imaging of thyroid cancer [34; 35]. [PERSON] and [PERSON] employed time-independent RTE for image reconstruction [36; 37]. For the computation of the RTE, they employed the upwind-difference discrete-ordinate method by discretizing the angular and spatial variables in the RTE. After experimental validation of the forward process with a rectangular phantom, including the void-like region [36], image reconstruction using the gradient method, referred to as model-based iterative image reconstruction (MOBIR), was performed. Images of \(\mu_{a}\) and \(\mu_{s}\) were obtained in a phantom containing void-like regions [37]. [PERSON] also reported MOBIR with a forward process using FEM [38]. [PERSON] et al. used the RTE in FD measurements during the forward process [39]. The RTE was computed using the FEM, in which stream diffusion modification was utilized. In the optimization process, the cost function with total variation regularization (see Section 4.1) was minimized using the Gaussian-Newton method. Image reconstruction was successful in the 2D numerical simulations of various cases, including low-scattering blobs and low-scattering and low-absorption gaps. [PERSON] and [PERSON] proposed RTE-based image reconstruction for the FD measurement of an object, such as an embryo that consists of weakly and highly scattering regions [40]. Assuming that \(\mu_{a}+\mu_{s}\) is significantly different between the weak and high scattering regions, an approximation of the RTE solution was derived for the forward process. The optimization process with Tikhonov regularization was demonstrated using 3D numerical simulations. [PERSON] et al. proposed image reconstruction with linearization employing the Rytov approximation accompanied by [PERSON]'s function calculated using the \(F_{N}\) method [41]. The \(F_{N}\) method was originally used to compute the 1D RTE using singular-eigenfunction expansion, which was extended to 3D computation by employing the rotated reference frame method [42]. The image reconstruction method was examined in phantom experiments using CCD measurements. [PERSON] also reported numerical simulations of image reconstruction based on the Born approximation and 3D FN computation of the RTE for a slab geometry with spatially oscillating structured light [43].