Experimental investigation of perturbation Monte-Carlo based derivative estimation for imaging low-scattering tissue

Experimental results for imaging the low-scattering tissue phantoms based on the derivative estimation through perturbation MonteCarlo (pMC) method are presented. It is proven that pMC-based methods give superior reconstructions compared to diffusion-based reconstruction methods. An easy way to estimate the Jacobian using analytical expression obtained from perturbation Monte-Carlo method is employed. Simulation studies on the same objects, considered in the experiment, are performed and corresponding results are found to be in reasonable agreement with the experimental studies. It is shown that inter-parameter cross talk in diffusion based methods lead to false results for the low-scattering tissue, where as the pMC-based method gives accurate results. © 2005 Optical Society of America OCIS codes: (170.0170) Medical optics and biotechnology; (100.3190) Inverse problems; (170.3660) Light propagation in tissues; (170.4580) Optical diagnostics for medicine; (170.5280) Photon migration; (170.6920) Time-resolved imaging; (170.7050) Turbid media. References and links 1. J. C. Hebden, D. A. Boas, J. S. George, and Anthony J. Durkin, “Topics in Biomedical Optics: Introduction,” Appl. Opt. 42, 2869-2870 (2003). 2. B. Chance, R. R. Alfano, B. J. Tromberg, M. Tamura, and E. M. Sevick-Muraca, ed., Optical Tomography and Spectroscopy of Tissue IV, Proc. SPIE 4250 (2001). 3. B. Chance and R. R. Alfano, ed., Optical tomography and spectroscopy of tissue: Theory, instrumentation, model, and human studies, Proc. SPIE 2979 (1997). 4. B. Chance and R. R. Alfano, ed., Optical tomography, photon migration and spectroscopy of tissue and model media: theory, human studies, and instrumentation, Proc. SPIE 2389 (1995). 5. A. Yodh and B. Chance, “Spectroscopy and Imaging with diffusing light,” Phy. Today 48, 34–40 (1995). 6. J. C. Hebden, S. R. Arridge, and D. T. Delpy, “Optical imaging in medicine. I. Experimental techniques,” Phys. Med. Biol. 42, 825–840 (1997). 7. S. R. Arridge and J. C. Hebden, “Optical imaging in medicine. II. Modelling and reconstruction,” Phys. Med. Biol. 42, 841–853 (1997). 8. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41–R93 (1999). 9. R. Chandrasekhar, Radiation Transfer (Oxford, Clarendon, 1950). 10. A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE Press, New York, 1997). 11. M. Firbank, S. R. Arridge, M. Schweiger, and D. T. Delpy, “An investigation of light transport through scattering bodies with nonscattering regions,” Phys. Med. Biol. 41, 767–783 (1996). 12. J. Ripoll, N. Nieto-Vesperinas, S. R. Arridge, and H. Dehghani, “Boundary conditions for light propagation in diffuse media with nonscattering regions,” J. Opt. Soc. Am. A 17, 1671–1681 (2000). 13. H. Dehghani, S. R. Arridge, M. Schweiger, and D. T. Delpy, “Optical tomography in the presence of void regions,” J. Opt. Soc. Am. A 17, 1659–1670 (2000). 14. Y. Xu, Q. Zhang and H. Jiang, “Optical image reconstruction of non-scattering and low scattering heterogeneities in turbid media based on the diffusion approximation model,” J. Opt. A: Pure Appl. Opt. 6, 29–35 (2004). (C) 2005 OSA 7 February 2005 / Vol. 13, No. 3 / OPTICS EXPRESS 985 #5877 $15.00 US Received 30 November 2004; revised 24 January 2005; accepted 31 January 2005 15. E. Okada, M. Firbank, M. Schweiger, S. R. Arridge, M. Cope, and D. T. Delpy, “Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head,” Appl. Opt. 36, 21–31 (1997). 16. M. Wolf, M. Keel, V. Dietz, K. von Siebenthal, H. U. Bucher, and O. Baenziger, “The influence of a clear layer on near-infrared spectrophotometry measurements using a liquid neonatal head phantom,” Phys. Med. Biol. 44, 1743–1753 (1999). 17. H. Dehghani and D. T. Delpy, “Near-infrared spectroscopy of the adult head: effect of scattering and absorbing obstructions in the cerebrospinal fluid layer on light distribution in the tissue,” Appl. Opt. 39, 4721–4729 (2000). 18. E. Okada and D. T. Delpy, “Near-infrared light propagation in an adult head model. I. Modeling of low-level scattering in the cerebrospinal fluid layer,” Appl. Opt. 42, 2906–2914 (2003). 19. E. Okada and D. T. Delpy, “Near-infrared light propagation in an adult head model. II. Effect of superficial tissue thickness on the sensitivity of the near-infrared spectroscopy signal,” Appl. Opt. 42, 2915–2922 (2003). 20. A. D. Klose, V. Prapavat, O. Minet, J. Beuthan, and G. Muller, “RA diagnostics applying optical tomography in frequency-domain,” Proc. SPIE 3196, 194–204 (1997). 21. A. D. Klose, A. H. Hielscher, K. M. Hanson, and J. Beuthan, “Three-dimensional optical tomography of a finger joint model for diagnostic of rheumatoid arthritis,” Proc. SPIE 3566, 151–159(1998). 22. Y. Xu, N. V. Iftimia, H. Jiang, L. L. Key, and M. B. Bolster, "Imaging of in vitro and in vivo bones and joints with continuous-wave diffuse optical tomography," Opt. Express 8, 447-451 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-7-447 23. O. Dorn, “A transport-backtransport method for optical tomography,” Inverse Problems 14, 1107-1130 (1998). 24. A. D. Klose and A. H. Hielscher, “Iterative reconstruction scheme for optical tomography based on the equation of radiative transfer,” Med. Phys. 26, 1698-1707 (1999). 25. 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, 691-713 (2002). 26. A. D. Klose and A. H. Hielscher, “Optical tomography using the time-independent equation of radiative transfer – Part 2:inverse model,” J. Quant. Spectrosc. Radiat. Transf. 72, 715-732 (2002). 27. A. D. Klose and A. H. Hielscher, “Quasi-Newton methods in optical tomographic image reconstruction,” Inverse Problems 14, 387-403 (2003). 28. M. Firbank, E. Okada, and D. T. Delpy, “A theoretical study of the signal contribution of regions of the adult head to near infrared spectroscopy studies of visual evoked responses,” Neuroimage 8, 69–78 (1998). 29. S. R. Arridge, H. Dehghani, M. Schweiger, and E. Okada, “The finite element model for the propagation of light scattering media: a direct method for domain with nonscattering regions,” Med. Phys. 27, 252–264 (2000). 30. H. Jiang, "Optical image reconstruction based on the third-order diffusion equations ," Opt. Express 4, 241246 (1999), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-4-8-241 31. N. Herschkowitz, “Brain development in the fetus, neonate and infant, Biol. Neonate 54, 1-19 (1988) 32. S. R. Arridge and W. R. B. Lionheart, “Nonuniqueness in diffusion–based optical tomography,” Opt. Lett. 23, 882–884 (1998). 33. Y. Pei, H. L. Graber, and R. L. Barbour, "Normalized-constraint algorithm for minimizing inter-parameter crosstalk in DC optical tomography," Opt. Express 9, 97-109 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-9-2-97 34. J. C. Ye, K. J. Webb, R. P. Millane, and T. J. Downar, “Modified distorted Born iterative method with an approximate Frechet derivative for optical diffusion tomography,” J. Opt. Soc. Am. A 16, 1814-1826 (1999). 35. A. Sassaroli, C. Blumetti, F. Martelli, L. Alianelli, D. Contini, A. Ismaelli, and G. Zaccanti, “Monte Carlo procedure for investigating light propagation and imaging of highly scattering media,” Appl. Opt. 37, 7392– 7400 (1998). 36. C. K. Hayakawa, J. Spanier, F. Bevilacqua, A. K. Dunn, J. S. You, B. J. Tromberg, and V. Venugopalan, “Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues,” Opt. Lett. 26, 1335-1337 (2001). 37. Y. Phaneendra Kumar and R. M. Vasu, “Reconstruction of optical properties of low-scattering tissue using derivative estimated through perturbation Monte-Carlo method,” J. Biomed. Opt. 9, 1002-1012 (2004). 38. F. Natterer and F. Wübbeling, “A propagation-backpropagation method for ultrasound tomography,” Inverse Problems 11, 1225-1232 (1995). 39. S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. 20, 299-309 (1993). 40. A. H. Hielscher, A. D. Klose, and K. M. Hanson, “Gradient-based iterative image reconstruction scheme for time-resolved optical tomography,” IEEE Trans. Med. Imag. 18, 262-271 (1999). 41. J. Spanier and E. M. Gelbard, Monte Carlo Principles and Neutron Transport Problems (Addison-Wesley, Reading, Mass., 1969). 42. B. Jain, P. K. Gupta, V. A. Podzyavnikov, and V. K. Chevokin, “Development and characterization of a UVvisible streak camera and its use for time resolved fluorescence studies on human tissues,” Proc. National Laser Symposium, B.A.R.C., Mumbai, India, January 17-19, 1996, National Laser Program, Department of Atomic Energy, Government of India, C3-C4 (1996). (C) 2005 OSA 7 February 2005 / Vol. 13, No. 3 / OPTICS EXPRESS 986 #5877 $15.00 US Received 30 November 2004; revised 24 January 2005; accepted 31 January 2005 43. B. Kanmani and R. M. Vasu, “Diffuse optical tomography using intensity measurements and the a priori acquired regions of interest: theory and simulations,” Phy. Med. Biol. 50, 247-264 (2005). 44. D. A. Boas, J. P. Culver, J. J. Stott, and A. K. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10, 159-170 (2002), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-3-159


Introduction
Biomedical optical imaging [1][2][3][4][5][6][7][8] in the recent years has been the topic of much interest to many researchers around the globe.Near Infra-red (NIR) optical imaging is mainly used to image thick tissue such as human brain, breast and joints.Main aim of NIR imaging is to get a spatial distribution of optical properties from the measurements done on the boundary of tissue.Light propagation in tissue can be accurately described by Radiative Transfer Equation (RTE) [9], which is an integro-differential equation.RTE is difficult to solve either analytically or numerically.The diffusion approximation [10] is generally applied to RTE and solved therein to reconstruct the optical properties.Even though diffusion approximation will be accurately able to describe the light propagation in tissue in most of the cases, it fails in the cases of low-scattering tissue [11][12][13][14], such as Cerebral Spinal Fluid (CSF) layer of the brain [15][16][17][18][19] and Synovial fluid layer in the joints [20][21][22].Numerical solutions from RTE are proven to be computationally expensive [23][24][25][26][27]. Hybrid radiosity finite-element-and higherorder diffusion models [28][29][30] require a priori knowledge of the boundaries of the lowscattering regions.This a priori knowledge is difficult to obtain in most of the practical situations, for example neonatal brain grows immensely in both shape and complexity immediately after birth [31].One of the main issues of diffuse optical tomography is the inter parameter cross-talk [32][33][34], which is due to underlying non-uniqueness in the inverse problem.
To image these low scattering tissues, there have been attempts to incorporate Monte Carlo (MC) simulation for doing reconstruction problem [35][36][37] (inverse problem of optical tomography).MC simulation is considered to be equivalent of implementing the RTE, which is accurate to describe the light propagation in these cases.Here again, MC simulation takes large computational time involved for taking many millions of photons through the tissue.Inversion methods that make use of repeated application of the forward model and updating of the derivative needed for guiding the solution to convergence would become too expensive to be a useful option.This problem has been recently addressed [36,37] making use of the perturbation Monte-Carlo (pMC) method to extract Jacobian, which is the derivative of measurements (forward model) with respect to the optical absorption (µ a ) and scattering coefficient (µ s ).
An analytical expression for the perturbation in the detected photon weight consequent to changes in optical properties in a sub-region, borrowed from neutron-transport theory, is used to update both Jacobian and the computed forward data [36,37].After discretization, such a procedure will require only a single Monte-Carlo (MC) simulation with the derivative-as well as forward projection data extraction handled by the pMC, which is rapidly done.In Ref. [36], this is used to solve a simple two region inverse problem of photon migration in a heterogeneous slab of thickness comparable to the transport mean free path, an object which falls under the transport-regime, where diffusion equation (DE) fails to hold.An extension of the pMC approach to construct a Jacobian matrix with basic MC simulation as forward model is presented in Ref. [37] for use in a perturbation-based method to reconstruct transportregime low-scattering objects with discretized absorption-and scattering coefficient distributions.
With numerical simulations the utility of a pMC-based iterative reconstruction strategy is shown [37] to reconstruct low-scattering objects with single-or multiple inhomogeneities.Verification for the utility of the above algorithm by reconstructing µ a -and/or µ s distributions from experimental time-domain data obtained from low-scattering phantoms with single-or multiple inhomogeneities is presented here.The propagation-backpropagation (PBP) strategy [38], where data is handled one source location at a time is used in the iterative reconstruction scheme, which reduces the dimension of the problem and the overall computation time.Reconstructed images of tissue-equivalent phantoms using the experimental data are presented.Even with the total-intensity measurements it is demonstrated that, one can clearly distinguish absorption-and scattering heterogeneities, in which case, diffusion based reconstruction methods require additional constrains [33]. ) with large no. of photons and store n & l for each of the pixel-subregion.

W c
For the PBP approach data from the next source is used after each update, otherwise use full data The work presented here is as follows:-In section 2, a brief overview of the iterative reconstruction scheme is presented which uses the Jacobian to calculate the update vectors.A brief description of the calculation of the Jacobian matrix using the pMC is also given.Then, in section 3, the details of the time-domain experiments conducted are described.The tissueequivalent phantoms and the geometry of data collection are also given in section 3.In section 4, reconstruction results from the experimental data along with simulation results are presented.Finally, to show the effectiveness of this approach, one set of reconstruction results from the simulated data using the diffusion-equation (DE) based reconstruction method [39] are also presented.
Here only PBP method is used where data from a single source (and many detectors) are reconstructed to arrive at an updated object which is used as the initial estimate for reconstructing data from the next source.In Ref. [37], the effectiveness of the PBP approach is shown for doing inverse problem compared to conventional approach.Since the number of equations in the PBP approach is smaller, the reconstruction problem will be too ill-posed to converge to a solution when both scattering-and absorption inhomogeneities are considered simultaneously without a priori knowledge on location [37].For reduction in ill-posedness of the inverse problem and for better convergence, an a priori information is used in all the cases, except one, presented here.Simulation studies for these objects are also presented in this work (details are given in Section 4).As a part of the simulation studies, it is shown that with or without a priori information, a pMC-based approach will be able to give same quality of reconstructions.

Iterative reconstruction algorithm
Complete details of the iterative reconstruction algorithm are presented in Ref. [37].Here this is briefly explained.The steps involved in the iterative reconstruction procedure are shown in the flow chart of Fig. 1.It has two iterations, one main one (the outer iteration) and another subsidiary one (the inner iteration) [40].The i th solution vector (µ a i , µ s i ) given to the forward model (the Monte-Carlo procedure is repeated with 2.1 million photons to model the forward propagation) predicts the computed data (W c ).The experimental data (W e ) plugged in help us find ∆W (= W e -W c ) that is used to set-up the update-equation (Eq.( 1)).In this work, integrated (or DC) intensity (Eq.(4.8) in Ref. [8]) is used as the data.The inner iteration outputs the update vector (∆µ a , ∆µ s ), which is used to arrive at the new-solution µ a i+1 = µ a i + ∆µ a and µ s i+1 = µ s i + ∆µ s .The conjugate gradients squared (CGS) method is used to solve this minimization problem (inner iteration).
Analytic expressions for changes introduced in detected photon weights owing to optical property perturbation in a sub-region in the object are made use of to calculate both the rate of change of photon weights with respect to optical properties and the new photon weights [36,37,41].Specifically if µ a and µ s are perturbed in certain locations and become where n is the number of collisions the photon undergoes in the perturbed region, and l is the path-length traversed therein and µ t = µ a + µ s .Eq. ( 1) provides a way to estimate the changes in measured photon density (or weight) owing to changes in µ a and µ s in a specified location.
Here the object is discretized into 81x81 pixels with freedom for µ a and µ s values in these pixels to move towards their values at convergence.If data from only M detectors is considered for each of the S source positions the Jacobian matrix connecting the change in measurements at the boundary to perturbations in µ a and µ s , has dimensions either (SxM)x{(81x81)x2} or Mx{(81x81)x2} (when data from each source is handled separately which is the PBP approach ), depending on whether the handling of the data is from all the sources together or one source at a time.For constructing such a Jacobian, using Eq. ( 1), there are 2x(81x81) regions of possibly different µ a and µ s values centered around each pixel.Approximately circular sub-regions containing 12 pixels are introduced surrounding each pixel in the original object domain (for the boundary pixels this is done by extending the boundary), with the object assigned homogeneous background values, µ a b and µ s b for the calculations of the Jacobian.A large number of photons (28.8 million) are taken through the object to determine the average number of collisions, n, and the length traversed, l, in each of the sub-regions entered around each of the pixels for the detected photons.The set of n & l values determined are frozen and used to update the Jacobian matrix during the course of the iterations.To arrive at a particular element of the Jacobian matrix, which is the rate of change of data at a detector with respect to the optical properties at a certain pixel, the stored n & l values corresponding to that pixel sub-region are used to determine (1).More details about the weight matrix (Jacobian) generation can be found in Section 3 of Ref. [37].

Phantom fabrication
Tissue-equivalent phantoms are fabricated as described below: The background material is obtained by mixing 100 parts of Lapox A-53 epoxy resin with 10 parts of K-6 hardener (Cibatul Limited, Gujarat, India).Scattering is introduced by adding TiO 2 powder (Dupont, Ti-Pure R-902) and absorption through India ink of appropriate quantities as required.
The phantoms with inhomogeneities are fabricated in two steps: In the first a cylinder is cast which is designed to have background optical properties, µ a b = 0.008 mm -1 and µ s b = 0.05 mm -1 .For this, while mixing the epoxy resin to hardener 469 mg of TiO 2 powder and 8 mg of 2% India ink are also added and thoroughly mixed.The resin is allowed to set for 24 hours at room temperature (approx.26 o C).Inhomogeneity is introduced as smaller diameter rods of different optical properties.In the experiments conducted here, the inhomogeneity values at its centre in absorption (µ a in ) and scattering (µ s in ) are varied from 0.008-0.021mm -1 and 0.005-0.14mm -1 respectively.For this in a 110ml mixture of epoxy resin and hardener, 262mg TiO 2 and 1.6mg 2% Indian ink are added.Both the background cylinder and the rods are thoroughly degassed in a vacuum chamber.The rod, which is to serve as the inhomogeneity, is machined and fitted exactly into a hole drilled in the background cylinder.A composite object with multiple inhomogeneities is also prepared following the same procedure.
The details of tissue-equivalent phantoms fabricated using the above procedure are as follows: Overall diameter = 80 mm.Inclusion diameter = 10 mm.The background optical properties, µ a b = 0.008 mm -1 and µ s b = 0.05 mm -1 .The optical properties of the absorptionand scattering inclusion are µ a in = 0.021 mm -1 and µ s in = 0.14 mm -1 respectively.The refractive index of the material, which is considered uniform, is 1.53.Henyey-Greenstien phase function is used to describe the scattering.The anisotropy factor g is kept at 0.75.

Data gathering
A schematic of the setup used for experiments is shown in Fig. 2. The second harmonic output of a ps Nd:YAG laser (Continuum PY61C-10, 532 nm, 27 ps, 10 Hz) is used to pump a 100 cm long Raman cell filled with H 2 at 30 Atm. pressure to generate 683 nm Stokes output by Stimulated Raman Scattering (SRS).Typical conversion efficiency of the set up is ~20% with pulse-to-pulse fluctuation of ~ 8% in first stokes output at 683nm.The Stokes output is separated from pump and the higher order Stokes and anti-Stokes output by proper high and low pass filters.The spatially filtered Stokes beam is used to illuminate the object mounted on a jig with provision for rotating to any angle.A multimode fiber (Φ = 1000 µm) of 1 meter length is used for collecting the scattered light exiting from the phantom.The detector fiber position can be independently adjusted to any position around the phantom.The Stokes light emerging from the scattering medium gets temporally elongated due to multiple scattering.A single shot streak camera [42] is used to record the photon arrival histogram (the TPSF).Data is collected for 13 detector positions, with an angular spacing of 5 o covering an angle of 60 o diametrically opposite to the source position for each of 12 source positions at spacing of 30 o covering the phantom fully.Averaging of 50 data sets for each source-detector location is performed, after making suitable correction for temporal jitter, before the data is subjected to further analysis.

Experiments
An initial Monte-Carlo simulation with 28.8 million photons having an extended boundary and homogeneous background values of µ a b and µ s b is performed to store n & l, the average number of collisions and length of traverse of photons, in the 12-pixels circular sub-regions surrounding each of the pixels.This data set is needed for calculation of the Jacobian as given in section 2 (see Fig. 1).For forward computed data (W c ), integrated intensity, an MC simulation with 2.1 million photons is done on the object with same source-detector positions with the same detector diameter.PBP strategy is used in combination with a priori information to reduce the computation effort.Here (in PBP) data from one view is considered sufficient to solve the iteration, which resulted in an updated object for use with data from the next view.All the views are considered cyclically.With µ a -and µ s inhomogeneities simultaneously present, without assumption of locations, the PBP-based problem is too illposed to result in a solution [37].The convergence criterion used is that the norm of the difference between computed data and the experimental data should be below a, preset, small value.Since Monte-Carlo simulation is used, a stochastic forward model, the reconstructed images required post-processing.A local median filter of dimension 5x5 pixels is employed to smoothen the images.
As stated earlier, experimental data (integrated intensity) is collected for three phantoms: first one with only absorption inhomogeneity, the second only with scattering inhomogeneity and third one with both at different locations.In each of the cases, the TPSF's measured are integrated to obtain the total photon weight (integrated intensity) measured for that sourcedetector pair (to serve as W e ).Data from a particular view is input to the iteration of Fig. 1, which outputs the update vector [∆µ a , ∆µ s ].The updated µ a and µ s are used along with the data from the next 'view' to continue the iterations (PBP approach).In all the experimental cases considered, it is assumed that the location of inhomogeneity is a priori known to be within a region of 30x30 pixels.The initial estimates of (µ a , µ s ) for the reconstruction algorithm Figure 3 gives the reconstruction of µ a -distribution from the experimental data obtained from the tissue phantom with only µ a -inhomogeneity.The unknowns are the value of µ a at these spatial a priori known locations.What is assumed a priori is that the inhomogeneity is somewhere inside a region of 30x30 pixels centered at (50, 41) making number of unknowns as 900.As PBP approach is used Jacobian for this problem has a dimension of 13x900.The original object shown in Fig. 3(a) has a µ a -inhomogeneity.The reconstructed µ a -distribution is

Simulations
Results of the numerical simulations for the same objects considered to those used in the experiments are presented here.The detectors are considered to be of diameter 1mm and simulated experimental data, (W e ), are generated by adding 2% Gaussian noise to the integrated intensity arrived at using MC simulation.Similar to the experimental reconstruction procedure, the homogeneous background optical properties  In the numerical simulation, for the object with only absorption inhomogeneity (Fig. 3(a)), using the PBP approach and the a priori assumption of inhomogeneity location, convergence is obtained in 16 iterations.The result is shown in Fig. 6(a).Without a priori information, the reconstruction problem took longer time (29 iterations, each taking ≈ 43 minutes) to converge and the result, shown in Fig. 6(b), is comparable in accuracy to the one shown in Fig. 6(a).Reconstructed µ a -inhomogeneity value at the center of Fig. 6(a) and 6(b) are 0.018 mm -1 and 0.019 mm -1 respectively.For this reason, in all the other cases considered here, a priori location information is assumed to be known.The reconstruction result of simulations corresponding to the object shown in Fig. 4 Reconstructed µ s -image without a priori information about the location of the inhomogeneity (0.13 mm -1 at the centre).Finally to bring out the effectiveness of pMC-based algorithm, experimental data is used in a diffusion equation-based algorithm, which did not converge to meaningful results.Even with the 2% noise added simulation data, the DE-based algorithm has failed to converge.Without any noise in the simulated integrated intensity, as a test case, the object with two inhomogeneities (Fig. 5

Discussion and conclusions
With these experimental results a confirmation of the usefulness of a Monte-Carlo-based method to reconstruct low-scattering objects, where the diffusion equation fails to hold, is done.The novelty of this method lies in getting the derivative information using a simple analytical equation (Eq. 1) and using the same in subsequent iterations.Further, experimentally, it is shown that cross talk between absorption-and scattering coefficients is negligible in the pMC-based reconstructions.The DE-based reconstructions are prone to cross-talk as the diffusion-and absorption coefficients are not independent of each other.Therefore µ a-and µ s reconstructions based on sensitivity matrices calculated with respect to absorption and diffusion coefficients will have cross-talk which needs to be minimized by adding correction terms to the sensitivity matrix [33,34].In addition, it is also seen that with the DE a measurement data type like intensity has a preferential weightage towards reconstructing µ a and has poor sensitivity towards µ s changes [34].This is not so with an RTE or MC-based method as demonstrated through simulations and experiments here.In comparison, a large false positive is seen in the µ s image at the µ a -inhomogeneity location in the DE-based reconstruction shown in Fig. 10.Usefulness of dividing the reconstruction into a number of sub-problems based on angle of projection (view) in reducing computational complexity is demonstrated which give accurate reconstructions when location of inhomogeneities is a priori known.
The main applications of optical tomography being brain and breast imaging [1][2][3][4][5][6][7], a great challenge in the medical optical imaging is to develop a model which can reconstruct optical properties with a large range of variation.An example is the low-scattering CSF layer surrounding a highly scattering inhomogeneous region.To generalize, the typical problem would be to reconstruct low-scattering inclusions in an otherwise high scattering background.Whereas the RTE, and its equivalent the MC simulation, would hold good for both high and low-scattering regions, the DE is valid only where scattering predominates.One strategy would be to use the DE, which is computationally efficient to implement and switch to the RTE when one encounter the inclusions.This can be successful only if one can provide the correct boundary and source conditions at the interface between the inclusions and the background, which in itself is complicated.Another alternative approach can be to use MC as the forward model and pMC to provide update to both the Jacobian and the computed forward data which is computationally feasible.Reconstruction is limited only to the regions of interest, namely, the low-scattering inclusions, which are a priori obtained by analyzing the intensity data gathered around the objects following the procedure suggested in Ref. [43].Identification of the regions of interest will drastically reduce the problem dimension, application of pMC to update the sensitivity matrix and the forward data will render computational cost manageable and the MC simulation which takes into account angles of scattering with its larger applicability will reconstruct low-scattering inclusions better, i.e., more accurately and with lesser inter-parameter cross-talk.The future direction for this work will be to demonstrate, through simulations and experiment, imaging of low-scattering inclusions in a high scattering background using the above strategy, which would justify the practical utility of the method presented here.
In Ref. [44] a three-dimensional (3D) Monte-Carlo code has been presented for modeling photon transport in adult head.Availability of such efficient 3D models pave the way for extending the 2D demonstration presented here to 3D imaging and using it for practical brain tissue characterization.This method could be combined with and guided by other imaging modalities such as MRI, X-ray CT and ultrasound to complement its capabilities to realize a superior imaging modality which gives more accurate and useful spatial as well as structural information of the human brain than what each one of them can individually provide.

Fig. 1 .Find
Fig. 1.The flowchart used in the iterative reconstruction.The inner loop uses a gradient search algorithm to output update vectors for the optical properties.

Fig. 3 .Fig. 4 .
Fig. 3. Reconstruction of µ a -distribution obtained from PBP approach with experimental data from tissue-equivalent phantom with only µ a -inhomogeneity (a).Original µ a -distribution: Background µ a b and µ s b are 0.008 mm -1 and 0.05 mm -1 respectively, and the inclusion has µ a in = 0.021 mm -1 and µ s in = 0.05 mm -1 .(b).Reconstruction of (a).The reconstructed µ ainhomogeneity value at its centre is µ a in = 0.028 mm -1 .

Fig. 6 .
Fig. 6.Simulation results for Fig. 3(a).(a).Reconstructed µ a -image with a priori information about the location of the inhomogeneity (0.018 mm -1 at the centre).(b).Reconstructed image without a priori information about the location of the inhomogeneity (0.019 mm -1 at the centre).

Fig. 7 .
Fig. 7. Simulation result for Fig. 4(a).Reconstructed µ s -image with a priori information about the location of the inhomogeneity (0.16 mm -1 at the centre).

Fig. 8 .
Fig. 8. Simulation results for Fig. 5(a) & (b).(a).Reconstructed µ a -image with a priori information about the location of the inhomogeneity (0.019 mm -1 at the centre).(b).Reconstructed µ s -image without a priori information about the location of the inhomogeneity (0.13 mm -1 at the centre).

Fig. 9 .
Fig. 9. Comparison of horizontal cross-sections at y = 41 mm of the reconstructed images.(a).Horizontal cross-sections at y = 41 through the centre of images Fig. 3(b) and 6(a).(b).Horizontal cross-sections at y = 41 mm through the centre of images Fig. 4(b) and 7. (c).Horizontal cross-sections at y = 41 mm through the centre of images Fig. 5(c) and 8(a).(d).Horizontal cross-sections at y = 41 mm through the centre of images Fig. 5(d) and 8(b).
(a) and 5(b)) are reconstructed and the results are presented in Fig. 10(a) and 10(b).Even here the reconstruction was carried out with the same a priori information as in Fig. 5(c)and 5(d).The centre of the reconstructed inhomogeneities values in µ a -(Fig.10(a)) and µ s -distribution (Fig. 10(b)) are 0.011 mm -1 and 0.06 mm -1 respectively.Reconstructed µs-inhomogeneity (Fig. 10(b)) is falsely shown in the reconstruction result (original is as shown in Fig. 5(b)).