Image reconstruction in diffuse optical tomography based on simplified spherical harmonics approximation

The use of higher order approximations to the Radiative transport equation, through simplified spherical harmonics expansion (SPN) in optical tomography are presented. It is shown that, although the anisotropy factor can be modeled in the forward problem, its sensitivity to the measured boundary data is limited to superficial regions and more importantly, due to uniqueness of the inverse problem it cannot be determined using frequency domain data. Image reconstruction through the use of higher ordered models is presented. It is demonstrated that at higher orders (for example SP7) the image reconstruction becomes highly under-determined due to the large increase in the number of unknowns which cannot be adequately recovered. However, reconstruction of diffuse parameters, namely optical absorption and reduced scatter have shown to be more accurate where only the sensitivity matrix used in the inverse problem is based on SPN method and image reconstruction is limited to these two diffuse parameters. ©2009 Optical Society of America OCIS codes: (100.3190) Inverse problems; (170.3660) Light propagation in tissues References and links 1. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50(4), R1–R43 (2005). 2. D. R. Leff, O. J. Warren, L. C. Enfield, A. Gibson, T. Athanasiou, D. K. Patten, J. Hebden, G. Z. Yang, and A. Darzi, “Diffuse optical imaging of the healthy and diseased breast: a systematic review,” Breast Cancer Res. Treat. 108(1), 9–22 (2008). 3. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Process. Mag. 18(6), 57–75 (2001). 4. R. Choe, A. Corlu, K. Lee, T. Durduran, S. D. Konecky, M. Grosicka-Koptyra, S. R. Arridge, B. J. Czerniecki, D. L. Fraker, A. DeMichele, B. Chance, M. A. Rosen, and A. G. Yodh, “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med. Phys. 32(4), 1128–1139 (2005). 5. H. Dehghani, B. W. Pogue, S. P. Poplack, and K. D. Paulsen, “Multiwavelength three-dimensional near-infrared tomography of the breast: initial simulation, phantom, and clinical results,” Appl. Opt. 42(1), 135–145 (2003). 6. H. Jiang, Y. Xu, N. Iftimia, J. Eggert, K. Klove, L. Baron, and L. Fajardo, “Three-dimensional optical tomographic imaging of breast in a human subject,” IEEE Trans. Med. Img. 20(12), 1334–1340 (2001). 7. L. C. Enfield, A. P. Gibson, N. L. Everdell, D. T. Delpy, M. Schweiger, S. R. Arridge, C. Richardson, M. Keshtgar, M. Douek, and J. C. Hebden, “Three-dimensional time-resolved optical mammography of the uncompressed breast,” Appl. Opt. 46(17), 3628–3638 (2007). 8. B. W. Zeff, B. R. White, H. Dehghani, B. L. Schlaggar, and J. P. Culver, “Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography,” Proc. Natl. Acad. Sci. U.S.A. 104(29), 12169–12174 (2007). 9. D. A. Boas, K. Chen, D. Grebert, and M. A. Franceschini, “Improving the diffuse optical imaging spatial resolution of the cerebral hemodynamic response to brain activation in humans,” Opt. Lett. 29(13), 1506–1508 (2004). 10. J. P. Culver, A. M. Siegel, J. J. Stott, and D. A. Boas, “Volumetric diffuse optical tomography of brain activity,” Opt. Lett. 28(21), 2061–2063 (2003). 11. J. C. Hebden, A. Gibson, T. Austin, R. M. Yusof, N. Everdell, D. T. Delpy, S. R. Arridge, J. H. Meek, and J. S. Wyatt, “Imaging changes in blood volume and oxygenation in the newborn infant brain using three-dimensional optical tomography,” Phys. Med. Biol. 49(7), 1117–1130 (2004). #117752 $15.00 USD Received 25 Sep 2009; revised 13 Nov 2009; accepted 12 Dec 2009; published 18 Dec 2009 (C) 2009 OSA 21 December 2009 / Vol. 17, No. 26 / OPTICS EXPRESS 24208 12. A. Y. Bluestone, G. Abdoulaev, C. Schmitz, R. L. Barbour, and A. H. Hielscher, “Three-dimensional optical tomography of hemodynamics in the human head,” Opt. Express 9(6), 272–286 (2001). 13. A. D. Klose, and A. H. Hielscher, “Fluorescence tomography with simulated data based on the equation of radiative transfer,” Opt. Lett. 28(12), 1019–1021 (2003). 14. V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed. Eng. 8(1), 1–33 (2006). 15. E. M. Sevick-Muraca, G. Lopez, J. S. Reynolds, T. L. Troy, and C. L. Hutchinson, “Fluorescence and absorption contrast mechanisms for biomedical optical imaging using frequency-domain techniques,” Photochem. Photobiol. 66(1), 55–64 (1997). 16. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50(4225), 4241 (2005). 17. C. H. Contag, and M. H. Bachmann, “Advances in in vivo bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4(1), 235–260 (2002). 18. C. Kuo, O. Coquoz, T. L. Troy, H. Xu, and B. W. Rice, “Three-dimensional reconstruction of in vivo bioluminescent sources based on multispectral imaging,” JBO 12, 24007 (2007). 19. G. Wang, W. Cong, K. Durairaj, X. Qian, H. Shen, P. Sinn, E. Hoffman, G. McLennan, and M. Henry, “In vivo mouse studies with bioluminescence tomography,” Opt. Express 14(17), 7801–7809 (2006). 20. H. Dehghani, S. C. Davis, and B. W. Pogue, “Spectrally resolved bioluminescence tomography using the reciprocity approach,” Med. Phys. 35(11), 4863–4871 (2008). 21. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15(2), R41–R93 (1999). 22. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near Infrared Optical Tomography using NIRFAST: Algorithms for Numerical Model and Image Reconstruction Algorithms,” Communications in Numerical Methods in Engineering (2008). 23. A. H. Hielscher, A. D. Klose, and K. M. Hanson, “Gradient-based iterative image reconstruction scheme for time-resolved optical tomography,” IEEE Trans. Med. Imaging 18(3), 262–271 (1999). 24. J. Riley, H. Dehghani, M. Schweiger, S. Arridge, J. Ripoll, and M. Nieto-Vesperinas, “3D optical tomography in the presence of void regions,” Opt. Express 7(13), 462–467 (2000). 25. S. R. Arridge, H. Dehghani, M. Schweiger, and E. Okada, “The finite element model for the propagation of light in scattering media: a direct method for domains with nonscattering regions,” Med. Phys. 27(1), 252–264 (2000). 26. H. Dehghani, S. R. Arridge, M. Schweiger, and D. T. Delpy, “Optical tomography in the presence of void regions,” J. Opt. Soc. Am. 17(9), 1659–1670 (2000). 27. A. H. Hielscher, R. E. Alcouffe, and R. L. Barbour, “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol. 43(5), 1285– 1302 (1998). 28. J. Ripoll, M. Nieto-Vesperinas, S. R. Arridge, and H. Dehghani, “Boundary conditions for light propagation in diffusive media with nonscattering regions,” J. Opt. Soc. Am. A 17(9), 1671–1681 (2000). 29. K. M. Case, and P. F. Zweifel, Linear Transport Theory (Addidon-Wesley, Reading, MA, 1967). 30. H. B. Jiang, “Optical image reconstruction based on the third-order diffusion equations,” Opt. Express 4(8), 241– 246 (1999). 31. S. Wright, M. Schweiger, and S. R. Arridge, “Reconstruction in optical tomography using the P-N approximations,” Meas. Sci. Technol. 18(1), 79–86 (2007). 32. E. D. Aydin, C. R. E. de Oliveira, and A. J. H. Goddard, “A comparison between transport and diffusion calculations using a finite element-spherical harmonics radiation transport method,” Med. Phys. 29(9), 2013– 2023 (2002). 33. A. D. Klose, U. Netz, J. Beuthan, and A. H. Hielscher, “Optical tomography using the time-independent equation of radiative transfer — Part 1: forward model,” J. Quant. Spectrosc. Radiat. Transf. 72(5), 691–713 (2002). 34. S. Cahandrasekhar, Radiative Transfer (Clarendon Press, London, 1950). 35. A. D. Klose, and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comput. Phys. 220(1), 441–470 (2006). 36. M. Chu, K. Vishwanath, A. D. Klose, and H. Dehghani, “Light transport in biological tissue using threedimensional frequency-domain simplified spherical harmonics equations,” Phys. Med. Biol. 54(8), 2493–2509 (2009). 37. G. Marquez, L. V. Wang, S. P. Lin, J. A. Schwartz, and S. L. Thomsen, “Anisotropy in the absorption and scattering spectra of chicken breast tissue,” Appl. Opt. 37(4), 798–804 (1998). 38. S. Nickell, M. Hermann, M. Essenpreis, T. J. Farrell, U. Krämer, and M. S. Patterson, “Anisotropy of light propagation in human skin,” Phys. Med. Biol. 45(10), 2873–2886 (2000). 39. S. R. Arridge, and M. Schweiger, “Photon-measurement density functions. Part2: Finite-element-method calculations,” Appl. Opt. 34(34), 8026–8037 (1995). 40. S. R. Arridge, and W. R. B. Lionheart, “Nonuniqueness in diffusion-based optical tomography,” Opt. Lett. 23(11), 882–884 (1998). 41. H. Dehghani, B. W. Pogue, J. Shudong, B. Brooksby, and K. D. Paulsen, “Three-dimensional optical tomography: resolution in small-object imaging,” Appl. Opt. 42(16), 3117–3128 (2003). 42. M. E. Eames, and H. Dehghani, “Wavelength dependence of sensitivity in spectral diffuse optical imaging: effect of normalization on image reconstruction,” Opt. Express 16(22), 17780–17791 (2008). #117752 $15.00 USD Received 25 Sep 2009; revised 13 Nov 2009; accepted 12 Dec 2009; published 18 Dec 2009 (C) 2009 OSA 21 December 2009 / Vol. 17, No. 26 / OPTICS EXPRESS 24209


Introduction
Near infrared (NIR) optical tomography is a non-invasive imaging modality in which the optical properties within a volume of interest can be reconstructed using measured transmission and reflectance NIR data. These calculated optical maps can then be used to derive functional and structural information about the tissue being imaged [1][2][3].
In recent years, several NIR imaging systems have been developed with specific applications such as the detection and characterization of breast cancer [1,[4][5][6][7] and functional imaging of the brain [8][9][10][11][12]. More recently the use of biological markers for function specific optical imaging has led to the development of novel imaging systems and algorithms that aim to recover images of function specific activity within small animals using, for example, fluorescence [13][14][15] or bioluminescence markers [16][17][18][19][20].
Most NIR imaging systems rely on the use of a model based image reconstruction algorithm in which the light distribution within the imaging domain is approximated to provide a match to measured boundary data [21][22][23] and therefore providing an estimate to the optical parameters within the medium being imaged. The most commonly applied model for the estimation of NIR light propagation in tissue is achieved through the Diffusion Approximation (DA) to the Radiative Transport Equation (RTE) through its robustness in implementation and computational speed and flexibility.
In order to accurately apply the DA, scattering effects within the medium must be dominant over the absorption, i.e. µ s '>>µ a , where µ s ' is the reduced scattering coefficient and µ a is the absorption coefficient. Whilst these conditions are generally applicable in most soft tissue, the presence of non-scattering regions, such as the cerebrospinal fluid found within the brain, or small source-detector separation, as encountered in small animal imaging, can mean that the DA is not valid for all cases [24][25][26][27][28].
To accurately model NIR light propagation in small geometries, or models with high absorption (such as for bioluminescence imaging), higher ordered forward models are needed whereby the problem is no longer limited by the 'diffuse' approximation. Whilst Monte Carlo simulations which are based on probabilistic models can produce accurate results, they tend to be computationally slow to be of clinical use. Alternatively, numerical approximations to the RTE can be used whereby light propagation is estimated through well established particle transport models based on complex integral-differential equations, which can account for both spatial and directional component of light propagation. The spherical harmonics (P N ) approximation of the RTE [29], for example, expands the angular components of radiance into a series of spherical harmonics. This method has been successfully applied in optical imaging but although more appropriate in conditions where the DA is not valid, it is still less than ideal due to the computational cost incurred [29][30][31][32]. The discrete ordinates method (S N ) is also a common implementation of the RTE often encountered in nuclear physics which discretises the angular components into a number of discrete directions. This method has also been used successfully in optical imaging application but as in the latter case results in a much higher computational cost as compared to the DA [27,33,34].
More recently, the application of Simplified Spherical Harmonics (SP N ) approximation to the RTE has been applied and studied [35,36]. The use of SP N methods have been demonstrated to give accurate solutions for small geometries and in cases where the source/detector separation is small. The SP N method has the advantage of requiring (N + 1)/2 equations where N is the number of Legendre Polynomials. This is in contrast to (N + 1) 2 equations for the P N approximation, and N(N + 2) equations for the S N method, where N is the number of allowed directions.
Due to the computational complexities outlined above, the vast majority of existing model based reconstruction algorithms make use of the DA. The image reconstructions therefore aim to form images of the absorption coefficient, µ a and the diffusion coefficient κ = 1/3(µ a + µ s ') where µ s ' = (1-g)µ s and µ s and g are the scattering coefficient and anisotropic factor respectively. One drawback of the application of the DA is that through the introduction of the reduced scattering coefficient, information about the scattering phase function is lost [32] and assumed to be constant at g = 0.9 which is accepted as valid for soft tissue. Little work has been done to investigate the effect of varying anisotropy factor within image reconstruction and whether it is at all possible to approximate its value within the inverse problem. The effect of anisotropic structures that affect the propagation of near-infrared photons according to their direction is of current interest and structures have been experimentally shown to cause direction dependence of light propagation in cases of chicken breast tissue and human skin [37,38]. The theoretical considerations have so far concentrated on developing ways to model the anisotropic propagation of NIR light, and less emphasis has been put on exploring the implications of anisotropic effects in image reconstruction.
In this study, SP N based models are used to study the effects of anisotropic scattering factor on image reconstruction. More specifically the use of SP N methods have been utilized to evaluate the sensitivity of measured boundary data from a frequency domain (FD) measurement system (where amplitude and phase of the fluence are used at typically 100 MHz modulation frequency) to not only optical absorption and scatter but also anisotropy factor. Through calculation of error-norm (residual) from a theoretical model, the uniqueness of this problem is demonstrated. Furthermore, a framework for the use of SP N methods in image reconstruction is presented.

Forward model
The forward problem is modeled using the SP N approximation to the RTE. The SP N equations can be derived by taking the planar spherical harmonics approximation and replacing the 1D derivatives with their 3D counterparts [35,36]. The SP 7 approximation results in a set of four coupled equations: where Q is the isotropic source term, n φ are composite moments of the fluence and µ an is the n th order absorption coefficient given by: where the total attenuation coefficient µ t is given by the sum of the absorption and scattering coefficients (µ t = µ a + µ s ). In practice, however, it is only possible to experimentally record the total fluence which is calculated from the composite moments as If we assume that the composite moments are in fact measureable, the SP N approximation introduces the possibility of reconstructing a wider range of optical parameters. Whereas the DA allows the recovery of just µ a and µ s ', the SP N approximation allows the recovery of µ a , µ a1 , µ a2, µ a3 , … µ an and, from these values µ s and g can be calculated through Eq. (2). It is important to note that although it is possible to measure the angular dependence of the fluence at the boundary through the use of optical fibers whose numerical aperture can be adjusted, it is not clear how the composite moments as shown in Eq. (1) and 3 may be measured experimentally. Nonetheless, for the purpose of the theoretical study presented here, it will be assumed that these moments can be measured. The SP N equations have been implemented using the NIRFAST package and previously validated using Monte Carlo models [36].

Image reconstruction
The goal of the inverse problem is the recovery of the optical properties, µ, by minimizing the difference between the measured data Φ M and data calculated by the forward model Φ C using a modified Tinkhonov minimization approach given by [22] ( ) ( ) where NM is the number of measurements, NN is the number of nodes in the FEM, µ 0 is an initial estimate of the optical properties and λ is the Tikhonov regularization parameter. By considering only the first order terms and applying the Levenberg-Marquardt procedure, this leads to an update vector ( ) where 2 λ λ = and J is the Jacobian matrix ( ) δ δµ Φ , often referred to as the weight matrix or sensitivity matrix. The Jacobian defines the relationship between changes in calculated boundary data due to small changes in the optical properties. If both phase and amplitude data is available, the structure of the Jacobian becomes  I  I  I   I  I  I   J   I  I  I   δ  δ  δ  δµ  δµ  δµ  δθ  δθ  δθ  δµ  δµ  δµ  δ  δ  δ  δµ  δµ  δµ  δθ  δθ  δθ  δµ δµ δµ where δlnI i /δµ j are sub-matrices that define changes in the ith measurement of log amplitude due to a change in optical properties at node j, and δθ i /δµ j define changes in the ith measurement of phase due to changes in optical properties at node j. In this study, the Jacobian as presented in sections 3.1 -3.3 is constructed using the perturbation method and for section 3.4 these are calculated using the Reciprocity approach which was validated using the perturbation method [39]. As the optical properties of the FEM mesh were defined by absorption coefficient, scattering coefficient and anisotropy factor, the final Jacobian was made up of three separate kernels with the form

Sensitivity mapping
The sensitivity of the SP 7 models to changes in optical properties were studied using a circular geometry of 43 mm radius with 16 equally spaced and collocated sources and detectors as shown in Fig. 1. The Jacobian was constructed using the perturbation method before being mapped onto the mesh coordinates to produce a map of sensitivity to changes in optical properties. The reference data was calculated using homogeneous optical properties of µ a = 0.01 mm −1 , µ s = 10 mm −1 , g = 0.9 and refractive index (n) = 1.33 using a modulation frequency of 100MHz. The sources were modeled as point sources located at 1 mm inside the boundary, to correspond with one reduced scattering distance as in DA. The total sensitivity of log amplitude for all sources and detectors pairs to changes in absorption coefficient is mapped in Fig. 2(a). Areas of high sensitivity exist in the regions surrounding the sources and detectors with a much lower background sensitivity decreasing, as expected, towards the centre of the domain. The negative magnitudes of sensitivity confirm that increases in absorption lead to reductions in measurements of intensity. Figure 2(b) shows the sensitivity of same data-type to scattering coefficient and of similar trend to the previous case. The sensitivity scale for scattering is again negative but has a much lower magnitude than that of absorption. Figure 2(c) shows that the sensitivity of intensity measurements to changes in anisotropy factor and is again highly localized to the sources and detectors. The magnitude of sensitivity of g is several orders lower than that of both absorption and scattering.
Figures 2(a)-(b) show that increases in both the absorption and scattering coefficients lead to decreases in intensity measurements. This suggests that it is not possible to distinguish between the two effects. A decrease in intensity, for example, could be caused by either a small increase in absorption or a large increase in scattering.   It is evident that the sensitivity is greater in the regions surrounding the sources and detectors as compared to the regions in between them. Additionally, the sensitivity although increasing to a maximum at a few scattering depths inside, it then decreases towards the centre of the domain, which is consistent with previous findings [39]. The scattering sensitivity map, Fig. 3(b), shows regions of very high sensitivity near the sources and detectors with much lower sensitivity elsewhere. The magnitude of scattering sensitivity is again much lower than that of absorption but with a positive scale. Sensitivity of phase measurements to changes in anisotropy factor is again limited to the regions surrounding the sources and detectors with a much lower magnitude as compared to absorption and scatter.
The results show that using log amplitude data alone, it is not possible to separate changes due to perturbations in absorption and scatter properties since both result in the same effect (i.e. increase in either absorption or scatter results in a decrease in log amplitude measurements). With the introduction of phase data, however, it can be seen that changes in amplitude and scatter have opposite effects on such data and so it is possible to distinguish between perturbations in absorption or scatter. However, since scattering and anisotropy factor have opposite sensitivities for both data-types, an increase in one parameter can be compensated by a decrease in the other and vice-versa. Therefore the problem appears nonunique for reconstruction of scatter and anisotropy factor.
3.2 Uniqueness of boundary data Next the uniqueness of the imaging problem is analyzed. The problem of non-uniqueness arises when more than one set of optical property distributions lead to an identical set of boundary data, in which case, the simultaneous recovery of multiple optical properties, such as the absorption and scattering coefficients, is not possible [40].
The previously used circular mesh, Fig. 1, was used to test the uniqueness of the SP 7 equations. A circular anomaly, with radius 10 mm, was inserted into the mesh as shown in Fig. 4. The optical properties within the anomaly were perturbed with either the absorption coefficient ranging between 0.0125mm −1 to 0.0375mm −1 , scattering coefficient ranging between 10mm −1 to 30mm −1 and anisotropy factor ranged between 0 and 1 and corresponding boundary data for all possible combinations were calculated. These wide range of values are chosen in order to allow a comprehensive study of their effect on problem uniqueness. Reference data was also calculated using the optical properties of µ a = 0.025 mm −1 , µ s = 20 mm −1 and g = 0.9. The absolute error (L1 norm) for each data-type between the reference and each set of perturbed data was the defined as where ref ϕ is the reference homogeneous unperturbed boundary data (either phase or log amplitude) and µ ϕ is the perturbed data. Figure 5(a) shows the map of error (Eq. (9) in amplitude data for each combination of µ a and µ s . A large band of optical property combinations exists for which the error between the reference data and perturbed data falls to zero. These equivalent combinations lead to identical amplitude data at the boundary and, as such, amplitude data alone is insufficient to recover both absorption and scattering properties simultaneously. By introducing phase measurements, however, the number of non-unique solutions can be minimized. Figure 5(b) shows the error map of phase measurements for the same range of optical properties. In this case, the optical property combinations leading to identical boundary data lie within a much more localized region. By combining the two data types, however, Fig. 5(c), the solution can converge on a smaller range of optical properties that lead to the perturbed boundary data.
The study was then repeated with varying absorption coefficients and anisotropy factors. Figure 6(a) shows that a large band of non-unique optical property combinations exists for amplitude measurements. The errors in phase measurements, Fig. 6(b), show an even larger range of equivalent combinations. When combining the two data types, however, it can be seen that the number of a µ and g combinations leading to the reference data actually falls within a smaller narrow region, Fig. 6(c).
The map of error in amplitude measurements with varying µ s and g, are also shown in Fig.  7. Figure 7(a), shows a large band of equivalent combinations of the two parameters for log amplitude data. The error in phase data actually shows two large regions of non-unique optical property combinations. The band of equivalent combinations in the amplitude error map is parallel to one of the bands in the phase error map. As such, even when combining the two data types, a large band of equivalent combinations still exists, Fig. 7(c), and therefore, even with both phase and amplitude data, it is not possible to distinguish between changes due to s µ and g.
By direct visual comparison between Figs. 6(c) and 7(c), it is evident that there exists a large range of absorption, scatter and anisotropy factor that would give rise to the same boundary measurements, even when considering both data-types and hence indicating the non-uniqueness of the problem to reconstruct all three parameters.

Uniqueness with the introduction of higher order terms
The use of higher ordered SP N approximations introduces a number of new variables, µ an , which are defined as in Eq. (2). These higher ordered absorption coefficients are dependent on both µ s and g and as such, may make it possible to separate the effects due to changes in the two parameters, if it is assumed that the composite moments of the boundary data can be measured. To test the uniqueness of the higher absorption orders, the values of µ an for n = 1:4 (as available for SP 5 ) were calculated using optical properties of µ a = 0.01 mm −1 , µ s = 20 mm −1 and g = 0.5 as a reference. The values of µ an were then re-calculated with values of the scattering coefficient ranging from 1mm −1 to 40mm −1 and the anisotropy factor ranging from 0.1 to 0.9. Figure 8(a) shows the difference between the reference value of µ a1 and those calculated for each combination of scattering coefficient and anisotropy factor. There is a large band of optical property combinations for which the residual falls to zero and therefore result in the same value of µ a1 . Figure 8(b) shows the same differences between values of µ a2 . There is still a large number of scattering coefficient / anisotropy factor combinations that lead to the same value in this case, although slightly smaller than for µ a1 . Figures 8(c)-(d) show similar plots for µ a3 and µ a4 respectively. It can be seen that as higher order moments of absorption are considered, the number of scattering coefficient / anisotropy factor combinations leading to the same value reduces, but does not converge to a particular combination. As such, the scattering coefficient and anisotropy factors cannot be simultaneously identified.

Multi-parameter image reconstruction using the SP N approximation
SP N based image reconstruction algorithms were developed for N = 1, 3, 5 and 7 using the modified Tikhonov minimization method. Test boundary data was generated using the same circular geometry, Fig. 1, with µ a = 0.01 mm −1 , µ s = 10 mm −1 , g = 0.9 and n = 1.33. A highly absorbing anomaly, with µ a = 0.02 mm −1 and a highly scattering anomaly with µ s = 20 mm −1 , were inserted into the mesh as shown in Fig. 9. Both anomalies have the same anisotropic factor and refractive index as the background.
Each of the SP N reconstruction algorithms (for N = 1, 3, 5 and 7) was capable of reconstructing different optical properties, depending on the order N. In order to test the capability of each of the SP N models, the optical properties listed in Table 1 were reconstructed, and from these, values of µ a and µ s were extracted. For the inverse model, the same mesh as the forward model was used without added noise to ensure that any errors in the reconstructed optical maps were only due to the complexity of the higher ordered approximations. Each of the algorithms were allowed to continue until there was less than a 0.1% change in successive iterations. The initial regularization parameter λ (Eq. (5) for all algorithms was set to 10 times the maximum diagonal of the matrix J T J and was reduced by a factor of 10 0.25 at each iteration [22]. Figure 10(a) shows the reconstructed images generated using an SP 1 based reconstruction algorithm and SP 1 forward data. The location of the two anomalies has been accurately reconstructed. The background optical properties have been accurately recovered although the values of the two anomalies have been over-estimated. The reconstructed image shows a maximum absorption coefficient of 0.0249 mm −1 and a maximum scattering coefficient of 23.9 mm −1 . The SP 1 reconstruction required 32 iterations at 3.1 second per iteration to recover the two relevant unknowns.  The SP3 reconstruction (using SP3 data) Fig. 10(b) has also accurately recovered the location and shape of the two anomalies. The recovered absorption and scattering coefficients have underestimated the target values, returning 0.0189mm −1 and 17.3mm −1 respectively. Unlike the SP1 image, the SP3 reconstruction shows signs of cross talk between the absorption and scattering images. The SP3 reconstruction required 13 iterations at 7.3 second per iteration to recover the four relevant unknowns.
Reconstruction using the SP 5 model (using SP 5 data) is shown in Fig. 10(c). In this case, the optical properties have again been overestimated with a maximum absorption coefficient of 0.0224mm −1 and a maximum scattering coefficient of 21.4mm −1 . The cross talk between the absorption and scattering images is still present. The SP 5 reconstruction required 34 iterations at 54 second per iteration to recover the six relevant unknowns.
The optical maps generated by the SP 7 model (using SP 7 data), Fig. 10(d) contains artifact throughout the domain. The location of the highly absorbing target has been recovered with poor size and contrast accuracy. The reconstruction failed to recover the highly scattering target. It is likely that the failure of the SP 7 reconstruction is due to the highly underdetermined nature of the problem. The SP 7 model requires 8 unique variables to be calculated at each node, i.e. 14280 unknowns, based on just 1920 boundary measurements. The SP 7 reconstruction required 20 iterations at 52 second per iteration to recover the eight relevant unknowns.

Multi-parameter image reconstruction and hard priori information
The under-determined nature of the image reconstruction problems can be minimized with the use of prior information. MRI data can be used, for example, to determine regions of similar tissue types and therefore, regions of homogeneous optical properties [41]. This then reduces the problem of calculating optical properties at every node to just a small number of known regions.
The previous study was repeated using region-based reconstruction for the SP N methods with N = 1, 3 and 5. The SP 7 model has been omitted from further studies as the increase in accuracy over the SP 5 model has been shown to be very small and neglectable [35].
It can be seen (Fig. 11) that the use of prior information, whereby instead of reconstruction the unknown parameters at each spatial variable, homogeneous values are estimated for each unknown region (in this case 3) has enabled the SP 5 reconstruction (when using SP 5 forward data) to accurately recover both targets. This is due to the vast reduction in the number of unknowns from 6 times the number of nodes (10710) to 6 times the number of regions (18). The same findings (not shown) were found for SP 1 and SP 3 cases. Fig. 11. Recovered optical map generated using SP5 reconstruction with prior information.

Diffusion based image reconstruction using the SP N approximation
As discussed earlier, the SP N approximations, Eqs. (1)(a)-(d), are based on composite moments of fluence with the total flux calculated using Eq. (2). The SP N reconstruction algorithms introduced so far were therefore based on various moments of the phase and amplitude boundary data which attempted to recover all moments of absorption and scattering coefficients. In practice, however, there are no experimental methods to easily measure the angular components of amplitude and phase and so image reconstructions must be based on total values of fluence. A new reconstruction algorithm was therefore developed that used the SP 5 model but only required measurements of phase and amplitude as available from the total fluence and not composite moments.
The new reconstruction algorithm uses the SP 5 forward model to calculate the composite moments of fluence. The total fluence is then calculated using Eq. (3). The Jacobian is then constructed using Eqs. (8a) and (8b) where log of amplitude and absolute phase are extracted from the calculated total fluence. By eliminating the composite moments of fluence, however, it is only possible to reconstruct for µ a . and µ s (assuming g = 0.9). This new reconstruction algorithm is tested using SP 5 forward data and the resulting optical map is shown in Fig.  12(a). The reconstruction has performed well recovering the optical parameters within 10% of the target values. For comparison, the same forward data is also used to reconstruct optical parameters using the DA based algorithm, Fig. 12(b). In this case, the target values have been over-estimated and boundary artifact has been introduced, indicating the maximum errors seen from high order model mismatch, is as expected near the source / detector positions, which lead to image inaccuracy and artifacts.

Discussions and conclusions
Forward models for image reconstruction based on the simplified spherical harmonic approximation have been presented. Sensitivity maps for amplitude and phase data are calculated, Figs. 2 and 3, using SP 7 and it is shown that the sensitivity of measured boundary data for changes in the anisotropy factor are many orders of magnitude lower than those of absorption and scattering coefficients. It is also shown that the sensitivity to both scattering coefficient anisotropy factor is highly localized to the regions adjacent to the source and detector (Figs. 2 and 3). This suggests that beyond a few scattering lengths, the light distribution is fully diffuse and as such, distant changes in scattering and anisotropy factors have little effect on measured data. Although it is possible to normalize the sensitivity maps prior to use in image reconstruction [42], it is also evident that an increase in scattering parameter will have an effect on boundary data which is opposite of that due to anisotropy factor. This in turn indicates that an increase in one parameter can be compensated by a decrease on the other and vice-versa. This effect is unlike that seen for absorption and scattering, whereby the two parameters can be separated using log amplitude and phase data, but it is unclear whether all 3 parameters can be resolved using these 2 data-types. Through the use of error maps, whereby the difference in some reference data with respect to data in the presence of an anomaly is calculated, it is shown that using frequency domain data, it is possible to distinguish, with some accuracy, the absorption and scattering properties of the medium, Fig. 5, which has also been presented previously [21]. The same principles have been applied to higher order models and additionally the effect of the anisotropy factor has been studied as well as the optical absorption and scatter. It is demonstrated, Fig. 6 and 7, that using frequency domain data, a large and non-unique range of all three optical properties can lead to the same measured data, indicating that it is not possible to reconstruct all three parameters using only 2 data-types. This finding is in line with earlier results, Fig. 3, and demonstrate that it is difficult to separate for both scattering and anisotropic factors simultaneously using frequency domain data. Further work is needed to investigate whether other data-types, for example using time resolved data, and / or spectral data, would allow separation of these variables.
Residual maps for different orders of absorption coefficient (µ an ) as available using SP N models have been calculated to investigate the possibility of resolving scattering and anisotropy factors, if µ an coefficients can be measured and calculated. It is shown, Fig. 8, that even if all higher moments of the absorption coefficients can be calculated, there still exists a large range of scattering and anisotropy factors that would lead to the same µ an values, thus indicating that these cannot be easily resolved.
Reconstruction algorithms based on the SP N approximation have been developed and tested. For image reconstruction, the same FEM model as for the forward data was used and no noise was added. The inverse crime has been strictly committed, since we are concerned with the accuracy of each model in determining the unknown parameters associated with each model. Additionally, noise free data have only used to only highlight the error seen due to model mismatch, rather than the effects of noise within the data, which can show similar trends. It was shown, Fig. 10(a)-(c), that for N = 1, 3 and 5 the reconstructions performed well, with the reconstructed values being within 24% of expected values with the worse results being obtained from SP 1 and absorption coefficient. This accuracy could be further improved by the optimization of the regularization parameter and stopping criteria. The N = 7 model, however, failed to accurately recover the target values, Fig. 10(d). The absorption coefficient, although located, is extremely over estimated, with some cross-talk from the scattering value. The reconstructed scatter image contains high valued artifacts and can be considered inaccurate. By limiting the presented results to these larger domain, it is also worth noting that the majority of image artifacts are seen at the boundary, near the source locations, where SP 1 solution is known the be less accurate. It is therefore expected that the errors seen due to the lack of higher order approximations to be substantially more significant in small geometry imaging experiments. The poor performance of the SP 7 model was most likely due to the under-determined nature of the problem. As stated earlier, the SP 7 model contains 8 unique unknown variables which need to be calculated at each spatial location (FEM node), i.e. 14280 unknowns, based on just 1920 boundary measurements (240 log amplitude and 240 Phase, based on 16 co-located source and detectors and 4 composite moments of fluence). In order to eliminate the large degree of freedoms, the accuracy of the reconstruction algorithms was further improved by creating a better determined problem. This was achieved for the SP 1 , SP 3 and SP 5 approximations with the use of prior information. The prior information was used to identify a number of homogeneous regions, reducing the number of properties to be recovered. This simplification of the problem led to improvements in all three of the reconstructions (for N = 1, 3 and 5). Both the absorbing and scattering targets were accurately recovered by all models to within 2%, Fig. 11.
The SP N equations are based on composite moments of fluence, Eq. (2). In reality, however, experimental systems can only measure the total fluence at the boundary of the domain and so the full SP N reconstruction algorithms are of limited use in image reconstruction where all composite moments are needed. An alternative method has been proposed in which the forward models are based on absolute fluence but the Jacobian matrix was calculated using the SP 5 model to allow the calculation of absorption (µ a ) and scattering (µ s ) coefficient only. The diffusion parameter based images reconstructed from simulated SP 5 data whereby the Jacobian is based on either SP 1 or SP 5 , Fig. 12, indicate that although both models can be used, the higher order model is more accurate both in terms of quantitative and qualitative analysis. The target values of the test problem calculated using SP 5 are recovered with more accuracy as compared to the SP 1 based model with much less artifact. Additionally, eliminating the use of complex moments also results in decreased computation time and memory requirements.
The results presented in this work indicate that for image reconstruction whereby the DA is less valid, in for example, small animal imaging and / or where the absorption coefficient is more dominant, the higher order models based on simplified spherical harmonics can be used to generate the sensitivity matrix for diffusion based image reconstruction, without the additional computation complexity in terms of the number of unknown parameters. The incorporation of these more accurate models can however allow for a better accuracy in terms of light propagation models.