# A Review of Image Reconstruction Algorithms for Diffuse Optical Tomography
## 3 Forward Process
### \(P_{n}\) Approximation  
To obtain a good approximation of the RTE, \(P_{N}\) approximation [23] was attempted in the forward process of image reconstruction. The \(P_{N}\) approximation is derived from an orthogonal function expansion using spherical harmonics in Equation (10), in which \(I(\textbf{r,s,}t)\), \(q(\textbf{r,s,}t)\), and \(p(\textbf{s}\cdot\textbf{s}^{\prime})\)[23] are formulated as  
\[I(\textbf{r,s,}t)=\sum_{l=1}^{\infty}\sum_{m=-l}^{l}\left(\frac{2l+1}{4\pi} \right)^{\frac{1}{2}}\psi_{l,m}(\textbf{r,}t)Y_{l,m}(\textbf{s}), \tag{11}\]  
\[q(\textbf{r,s,}t)=\sum_{l=1}^{\infty}\sum_{m=-l}^{l}\left(\frac{2l+1}{4\pi} \right)^{\frac{1}{2}}q_{l,m}(\textbf{r,}t)Y_{l,m}(\textbf{s}), \tag{12}\]  
\[P(\textbf{s}\cdot\textbf{s}^{\prime})=\sum_{l=1}^{\infty}\sum_{m=-l}^{l} \Theta_{l,m}Y_{l,m}^{*}(\textbf{s}^{\prime})Y_{l,m}(\textbf{s}), \tag{13}\]  
where \(\psi_{l,m}\), \(q_{l,m}\), and \(\Theta_{l,m}\) are the coefficients, and \(Y_{l,m}\) the spherical harmonics of order \(l\) and degree \(m\), which are associated with the Legendre polynomial and orthonormal function. \(Y_{l,m}^{*}\) is the complex conjugate of \(Y_{l,m}\). Owing to the orthonormality of the spherical harmonics, the RTE, in which Equations (11)-(13) are substituted, is decoupled into \((N+1)!\) simultaneous equations by multiplying \(Y_{l,m}^{*}\) when the spherical harmonic expansion is terminated at \(l=N\). The \(P_{N}\) approximation is computed with simultaneous equations. [PERSON] et al. compared \(P_{3}\) approximation with DA for the determination of the optical properties from MC simulated measurements and from \(P_{3}\)-approximated measurements in the FD measurement [56]. \(\mu_{a}\) was more accurately determined using \(P_{3}\) approximation, although the DA estimated \(\mu_{s}\)' better than \(P_{3}\) approximation in the highly forward-scattering case. [PERSON] and [PERSON] derived a stable and computationally convenient higher-order DA (\(P_{3}\) approximation) with second-order spatial derivative terms and compared the measurements with a cylindrical phantom and the computational results with the higher-order DA and DA using FEM [57]. [PERSON] and [PERSON] developed a computer program based on their original code, called EVENT, to compute the \(P_{N}\) approximation with FEM and compared the computational results with the \(P_{N}\) approximation and with DA. They succeeded in reconstructingimages from MC-simulated measurements by employing \(P_{7}\) approximation in TD measurements [58]. [PERSON] developed an image reconstruction algorithm based on the third-order diffusion approximation (\(P_{3}\) approximation) with Marquardt and Tikhonov regularizations and succeeded in image reconstruction using 2D numerical simulations [59]. [PERSON] et al. implemented a 3D FEM computation of \(P_{3}\) approximation and compared \(P_{3}\) approximation with \(P_{1}\) approximation and MC simulation. They reconstructed an image of a 3D numerical phantom mimicking a finger joint with cartilage between the bones [60]. [PERSON] et al. performed image reconstruction implementing the FEM calculation of the \(P_{N}\) approximation in the forward process in FD measurement [61]. They employed the even-parity formulation of RTE, which was derived by introducing the light intensity, \(I^{\pm}(\mathbf{s})=\{I(\mathbf{s})\pm I(-\mathbf{s})\}/2\), and the light source, \(q^{\pm}(\mathbf{s})=\{q(\mathbf{s})\pm q(-\mathbf{s})\}/2\). Consequently, the RTE was transformed into the formula of \(I^{+}\) including the second-order spatial derivative. The simultaneous equations of \(P_{N}\) approximation, which are expressed in the vector matrix formula, are discretized using the FEM. Numerical simulations of image reconstruction employing the \(P_{1}\), \(P_{3}\), and \(P_{5}\) approximations were performed, with the measurements computed using the \(P_{7}\) approximation.