# A Review of Image Reconstruction Algorithms for Diffuse Optical Tomography
## 3 Forward Process
### Monte Carlo Simulation  
Generally, the MC method is a numerical technique used to solve deterministic problems such as partial differential equations using probabilistic methods [88; 89]. In the forward process of computing light propagation, the movements of photons with scattering and absorption by biological tissues were stochastically simulated. Monte Carlo modeling of light transport in multilayered tissues (MCML) developed by [PERSON] et al. and other software available online are widely used [90; 91; 92; 93]. In the MC method, a photon packet with unit energy \(w\) is launched into a medium comprising layers with optical properties. The length of the free path in which the photon packet is not scattered is determined by an exponentially distributed random variable. The ratio of photons maintaining the direction and surviving the absorption at distance \(l\) is given by \(p(l)=\mu_{t}\exp(-\mu_{l}l)\), where \(\mu_{t}=\mu_{a}+\mu_{s}\). The energy of the photon packet decreases with absorption at the position reached by the photon packet. The absorbed energy \((\mu_{a}/\mu_{t})w\) is recorded at this position. The photon packet with the remaining energy \((\mu_{s}/\mu_{t})w\) continues to travel. The change in the direction of the photon packet by a scattering event is determined randomly using the Henry-Greenstein function \(p(\textbf{s},\textbf{s}^{\prime})\). The photon packet is reflected randomly by obeying Fresnel reflectance at the boundaries of the layers. A large number of photon packets are traced such that the recorded absorbed light distribution converges to a good approximation of the light propagation. Comparisons between MC and DA have been reported in many studies [94; 95; 96; 97]. [PERSON] et al. proposed the perturbation MC (pMC) method for the computation of the Jacobian matrix \(l\) to solve the inverse problem using a gradient-based nonlinear optimization process [98]. In the pMC method, the standard MC method is first used to record the energies of the detected photon packets moving from the light source to detectors with scattering events. To compute the changes in the measurements, the recorded energy \(w\) of the photon traveling through the region with the perturbations of \(\mu_{s}\rightarrow\overline{\mu}_{s}\), \(\mu_{a}\rightarrow\overline{\mu}_{a}\), and \(\mu_{t}\rightarrow\overline{\mu}_{t}\) is modified as  
\[\overline{w}=w\bigg{(}\frac{\overline{\mu}_{s}/\overline{\mu}_{t}}{\mu_{s}/ \mu_{t}}\bigg{)}^{\zeta}\bigg{(}\frac{\overline{\mu}_{t}}{\mu_{t}}\bigg{)}^{ \zeta}\exp\{-(\overline{\mu}_{t}-\mu_{t})L\}, \tag{17}\]  
where \(L\) is the total length of the photon packet traveling through the region perturbed by \(\zeta\) scattering events. The matrix \(J\) can be obtained using Equation (17). The performance of the pMC method was numerically examined to determine the optical properties of an object with two layers. [PERSON] and [PERSON] demonstrated image reconstruction for a heterogeneous medium with the background and inclusions having a low \(\mu_{s}\) of approximately 0.5 mm\({}^{-1}\), because of which the DA did not hold [99]. A nonlinear iterative optimization method, called the conjugate gradient squared method, was employed. Matrix \(J\) was computed using the pMC method. \(\mu_{a}\) and \(\mu_{s}\) values of the inclusions were simultaneously reconstructed with/without prior information regarding the locations of the inclusions in the 2D numerical simulations. Image reconstruction using the pMC method succeeded in imaging inclusions that could not be reconstructed using DA. [PERSON] et al. simulated light propagation in a 3D realistic adult head model using the MC method [100]. Then, using the Rytov approximation, which is a series expansion of the logarithm of \(\Phi\), the linearized forward process was formulated, as shown in Equation (7).