Quantitative characterization of turbidity by radiative transfer based reflectance imaging

A new and noncontact approach of multispectral reflectance imaging has been developed to inversely determine the absorption coefficient of μa, the scattering coefficient of μs and the anisotropy factor g of a turbid target from one measured reflectance image. The incident beam was profiled with a diffuse reflectance standard for deriving both measured and calculated reflectance images. A GPU implemented Monte Carlo code was developed to determine the parameters with a conjugate gradient descent algorithm and the existence of unique solutions was shown. We noninvasively determined embedded region thickness in heterogeneous targets and estimated in vivo optical parameters of nevi from 4 patients between 500 and 950nm for melanoma diagnosis to demonstrate the potentials of quantitative reflectance imaging. © 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement OCIS codes: (170.0110) Imaging systems; (120.5700) Reflection; (290.4210) Multiple scattering; (170.6510) Spectroscopy, tissue diagnostics. References and links 1. H. C. van de Hulst, Multiple light scattering: tables, formulas, and applications (Academic Press, 1980). 2. P. Lemaillet, C. C. Cooksey, J. Hwang, H. Wabnitz, D. Grosenick, L. Yang, and D. W. Allen, “Correction of an adding-doubling inversion algorithm for the measurement of the optical parameters of turbid media,” Biomed. Opt. Express 9(1), 55–71 (2018). 3. B. Guo, Y. Wang, C. Peng, H. Zhang, G. Luo, H. Le, C. Gmachl, D. Sivco, M. Peabody, and A. Cho, “Laserbased mid-infrared reflectance imaging of biological tissues,” Opt. Express 12(1), 208–219 (2004). 4. T. Tkaczyk, M. Rahman, V. Mack, K. Sokolov, J. Rogers, R. Richards-Kortum, and M. Descour, “High resolution, molecular-specific, reflectance imaging in optically dense tissue phantoms with structuredillumination,” Opt. Express 12(16), 3745–3758 (2004). 5. D. Roblyer, R. Richards-Kortum, K. Sokolov, A. K. El-Naggar, M. D. Williams, C. Kurachi, and A. M. Gillenwater, “Multispectral optical imaging device for in vivo detection of oral neoplasia,” J. Biomed. Opt. 13(2), 024019 (2008). 6. H. M. Leung and A. F. Gmitro, “Fluorescence and reflectance spectral imaging system for a murine mammary window chamber model,” Biomed. Opt. Express 6(8), 2887–2894 (2015). 7. N. Dimitriadis, B. Grychtol, M. Theuring, T. Behr, C. Sippel, and N. C. Deliolanis, “Spectral and temporal multiplexing for multispectral fluorescence and reflectance imaging using two color sensors,” Opt. Express 25(11), 12812–12829 (2017). 8. A. Esteva, B. Kuprel, R. A. Novoa, J. Ko, S. M. Swetter, H. M. Blau, and S. Thrun, “Dermatologist-level classification of skin cancer with deep neural networks,” Nature 542(7639), 115–118 (2017). 9. B. Fei, G. Lu, X. Wang, H. Zhang, J. V. Little, M. R. Patel, C. C. Griffith, M. W. El-Diery, and A. Y. Chen, “Label-free reflectance hyperspectral imaging for tumor margin assessment: a pilot study on surgical specimens of cancer patients,” J. Biomed. Opt. 22(8), 1–7 (2017). 10. M. Pilz, S. Honold, and A. Kienle, “Determination of the optical properties of turbid media by measurements of the spatially resolved reflectance considering the point-spread function of the camera system,” J. Biomed. Opt. 13(5), 054047 (2008). Vol. 9, No. 5 | 1 May 2018 | BIOMEDICAL OPTICS EXPRESS 2081 #322632 https://doi.org/10.1364/BOE.9.002081 Journal © 2018 Received 6 Feb 2018; revised 20 Mar 2018; accepted 26 Mar 2018; published 4 Apr 2018 11. P. B. Jones, H. K. Shin, D. A. Boas, B. T. Hyman, M. A. Moskowitz, C. Ayata, and A. K. Dunn, “Simultaneous multispectral reflectance imaging and laser speckle flowmetry of cerebral blood flow and oxygen metabolism in focal cerebral ischemia,” J. Biomed. Opt. 13(4), 044007 (2008). 12. H. Cen and R. Lu, “Optimization of the hyperspectral imaging-based spatially-resolved system for measuring the optical properties of biological materials,” Opt. Express 18(16), 17412–17432 (2010). 13. E. Vitkin, V. Turzhitsky, L. Qiu, L. Guo, I. Itzkan, E. B. Hanlon, and L. T. Perelman, “Photon diffusion near the point-of-entry in anisotropically scattering turbid media,” Nat. Commun. 2, 587 (2011). 14. S. F. Bish, N. Rajaram, B. Nichols, and J. W. Tunnell, “Development of a noncontact diffuse optical spectroscopy probe for measuring tissue optical properties,” J. Biomed. Opt. 16(12), 120505 (2011). 15. S. Kawauchi, I. Nishidate, Y. Uozumi, H. Nawashiro, H. Ashida, and S. Sato, “Diffuse light reflectance signals as potential indicators of loss of viability in brain tissue due to hypoxia: charge-coupled-device-based imaging and fiber-based measurement,” J. Biomed. Opt. 18(1), 015003 (2013). 16. D. Abookasis, B. Volkov, and M. S. Mathews, “Closed head injury-induced changes in brain pathophysiology assessed with near-infrared structured illumination in a mouse model,” J. Biomed. Opt. 18(11), 116007 (2013). 17. B. S. Nichols, A. Llopis, G. M. Palmer, S. S. McCachren 3rd, O. Senlik, D. Miller, M. A. Brooke, N. M. Jokerst, J. Geradts, R. Greenup, and N. Ramanujam, “Miniature spectral imaging device for wide-field quantitative functional imaging of the morphological landscape of breast tumor margins,” J. Biomed. Opt. 22(2), 026007 (2017). 18. N. Joshi, C. Donner, and H. W. Jensen, “Noninvasive measurement of scattering anisotropy in turbid materials by nonnormal incident illumination,” Opt. Lett. 31(7), 936–938 (2006). 19. P. Ghassemi, L. T. Moffatt, J. W. Shupp, and J. C. Ramella-Roman, “A new approach for optical assessment of directional anisotropy in turbid media,” J. Biophotonics 9(1-2), 100–108 (2016). 20. J. Q. Lu, X. H. Hu, and K. Dong, “Modeling of the rough-interface effect on a converging light beam propagating in a skin tissue phantom,” Appl. Opt. 39(31), 5890–5897 (2000). 21. X. Ma, J. Q. Lu, and X. H. Hu, “Effect of surface roughness on determination of bulk tissue optical parameters,” Opt. Lett. 28(22), 2204–2206 (2003). 22. X. Ma, J. Q. Lu, H. Ding, and X. H. Hu, “Bulk optical parameters of porcine skin dermis at eight wavelengths from 325 to 1557 nm,” Opt. Lett. 30(4), 412–414 (2005). 23. C. Chen, J. Q. Lu, K. Li, S. Zhao, R. S. Brock, and X. H. Hu, “Numerical study of reflectance imaging using a parallel Monte Carlo method,” Med. Phys. 34(7), 2939–2948 (2007). 24. X. Chen, Y. Feng, J. Q. Lu, X. Liang, J. Ding, Y. Du, and X. H. Hu, “Fast method for inverse determination of optical parameters from two measured signals,” Opt. Lett. 38(12), 2095–2097 (2013). 25. X. Liang, M. Li, J. Q. Lu, C. Huang, Y. Feng, Y. Sa, J. Ding, and X. H. Hu, “Spectrophotometric determination of turbid optical parameters without using an integrating sphere,” Appl. Opt. 55(8), 2079–2085 (2016). 26. W. W. Hager and H. Zhang, “A new conjugate gradient method with guaranteed descent and an efficient line search,” SIAM J. Optim. 16(1), 170–192 (2005). 27. Y.-H. Dai and C.-X. Kou, “A nonlinear conjugate gradient algorithm with an optimal property and an improved Wolfe line search,” SIAM J. Optim. 23(1), 296–320 (2013). 28. C. Chen, J. Q. Lu, H. Ding, K. M. Jacobs, Y. Du, and X. H. Hu, “A primary method for determination of optical parameters of turbid samples and application to intralipid between 550 and 1630 nm,” Opt. Express 14(16), 7420–7435 (2006). 29. Z. Song, K. Dong, X. H. Hu, and J. Q. Lu, “Monte Carlo simulation of converging laser beams propagating in biological materials,” Appl. Opt. 38(13), 2944–2949 (1999). 30. H. Ding, J. Q. Lu, K. M. Jacobs, and X. H. Hu, “Determination of refractive indices of porcine skin tissues and intralipid at eight wavelengths between 325 and 1557 nm,” J. Opt. Soc. Am. A 22(6), 1151–1157 (2005). 31. G. Kortum, Reflectance spectroscopy. Principles, methods, applications (Springer, 1969). 32. V. R. Weidner and J. J. Hsia, “Reflection properties of pressed polytetrafluoroethylene powder,” J. Opt. Soc. Am. 71(7), 856–861 (1981). 33. S. L. Jacques, J. C. Ramella-Roman, and K. Lee, “Imaging skin pathology with polarized light,” J. Biomed. Opt. 7(3), 329–340 (2002). 34. R. Novak, “Loop Optimization for Divergence Reduction on GPUs with SIMT Architecture,” IEEE Trans. Parallel Distrib. Syst. 26(6), 1633–1642 (2015). 35. C. M. Balch, J. E. Gershenwald, S. J. Soong, J. F. Thompson, M. B. Atkins, D. R. Byrd, A. C. Buzaid, A. J. Cochran, D. G. Coit, S. Ding, A. M. Eggermont, K. T. Flaherty, P. A. Gimotty, J. M. Kirkwood, K. M. McMasters, M. C. Mihm, Jr., D. L. Morton, M. I. Ross, A. J. Sober, and V. K. Sondak, “Final version of 2009 AJCC melanoma staging and classification,” J. Clin. Oncol. 27(36), 6199–6206 (2009). 36. J. Q. Lu, C. Chen, D. W. Pravica, R. S. Brock, and X. H. Hu, “Validity of a closed-form diffusion solution in P1 approximation for reflectance imaging with an oblique beam of arbitrary profile,” Med. Phys. 35(9), 3979–3987


Introduction
Quantitative characterization of turbidity attracts intense research interests for tasks ranging from materials analysis to disease diagnosis [1,2].Imaging reflected light at multiple wavelengths provides an approach of big data (> 10 4 signals) for such characterization.
Benefits include noninvasive and noncontact nature of measurement and simplicity of instrumentation given rapid progress in sensor technology.Most approaches of reflectance imaging realized to date, however, rely on pattern analysis to characterize targets for recognition and classification [3][4][5][6][7][8][9].While these methods require no modeling of light-matter interaction, their applications depend sensitively on configurations of illumination and detection.Simplified models based on Beer-Lambert law or diffusion approximations have been developed to inversely determine optical parameters of absorption coefficient μ a and reduced scattering coefficient μ s ' [10][11][12][13][14][15][16][17].Improvement with the radiative transfer (RT) theory has been reported for inverse determination of μ a , scattering coefficient μ s and anisotropy factor g of homogeneous turbid samples [18,19].Despite these advances, optical model based reflectance imaging is far from the point for translation into practical systems to characterize turbidity according to the RT theory [1].First, accurate and fast general-purpose methods need to be established for imaging system calibration and modeling spatial distribution of reflected light that allow robust inverse solutions.It is also critically important to develop effective algorithms for determination of μ a , μ s and g from the measured image.Finally, a new method should have the capacity to determine spatial distribution of optical and size parameters in heterogeneous turbid targets, which is necessary in applications including diagnosis and delineate normal and cancerous tissues.
In this report, we present a new approach of reflectance imaging that combines multispectral image acquisition, RT theory based Monte Carlo (MC) simulations of light transport [20][21][22][23], graphics processing unit (GPU) implementation of MC simulation [24,25] and a conjugate gradient [26,27] based inverse algorithm to retrieve a set of target parameters P. The axial-symmetric Henyey-Greenstein function p HG (cosθ) was used as the scattering phase function that is fully characterized by the anisotropy factor g [1].Thus the optical parameters of a turbid target can be collectively represented by the notation P = {μ a , μ s , g} for a homogeneous phantom or P = {μ a1 , μ s1 , g 1 , μ a2 , μ s2 , g 2 } for a 2-zone heterogeneous phantom.Furthermore, size parameters like thickness D for an embedded zone in the 2-zone phantoms can be included in the inverse determination to yield morphology information.To validate the new method, measurements with homogeneous turbid phantoms have been carried out for comparison with the parameters of thin disk samples made of similar materials by a method of integrating sphere measurement and MC based inverse algorithm [28].Multispectral reflectance images have been acquired from 4 patients suspected of cutaneous melanoma to estimate optical parameters and examine their potentials for cancer diagnosis.

Measurement of reflectance and profile images
A multispectral system has been constructed to image light reflected from a target illuminated with an 175W xenon fiber optic light source (ASB-XE-175EX, Spectral Products).The beam output from the fiber was collimated and incident on the target with 2w = 14mm in diameter and about 2° in diverging angle.A steel stencil S of stripes lines was inserted in front of the collimating lens L to change the incident beam profile from top-hat to a grating pattern, as shown in Fig. 1.Interference filters of 10nm in bandwidth were assembled on two rotating wheels to produce monochromatic illumination between 500 to 950nm for phantom study in steps of around 20nm.A camera lens of 25mm in focal length and 20mm in aperture diameter A was used to relay the light distribution at the target surface to the sensor of a cooled 14-bit cooled CCD camera (KX32ME, Apogee Instruments).The illumination arm was set at an incident angle of θ 0 from the normal of target surface while the imaging arm was oriented along the normal to reduce contribution of specular reflection.
The imaging system was first aligned by adjusting h in Fig. 1(a) to ensure a conjugate relation between the object plane of z = 0 at the target surface and the image plane of camera sensor, which yielded the maximum contrast in acquired images with the grating pattern.Then reflection images I r (x, y; λ) were acquired from a sample target at selected wavelengths where (x, y) are the pixel coordinates of the camera sensor corresponding to the target surface of (x, y, 0) in the field-of-view (FOV).Afterward, the target was replaced with a diffuse reflectance standard to measure the incident beam profile as I s (x, y; λ).Two standards were used with calibrated reflectance R s (λ) centered at 40% or 80% supplied by vendor (SRS-40-020, SRS-80-020, Labsphere).In each multispectral measurement, background images of I rb (x, y) or I sb (x, y) were also acquired with the incident beam blocked.The measured reflectance images R m (x, y; λ) and beam profile images S m (x, y; λ) were obtained as follows max, ( , ; ) ( , ; ) ( , ; ) ( ), ( )

Calculation of reflectance images by Monte Carlo simulation
Modeling of light transport in a turbid target to obtain calculated reflectance image was carried out with an individual photon tracking MC (iMC) code by tracking individual photons according to the RT theory.The code has been developed on the basis of previous ones [20-25, 28, 29] that tracks photons in the finite volume of the target as depicted in Fig. 1(a).The iMC code was validated against the numerical solution of RT equation in the cases of homogeneous turbid phantoms [1] (data not shown).With the individual photon tracking design, the iMC code can be implemented for parallel computing by GPU and perturbation based inverse determination of optical parameters [24,25].We briefly describe the logic flow of the iMC code here for simulation of reflectance imaging with details of the algorithm and validation results given elsewhere [20,24,28,29].It first imports the measured profile S m (x, y; λ) of the incident beam to generate a proportional photon density distribution ρ i (x, y, 0_; λ) at each grid point of the target surface or air-side of x-y plane at z = 0. Photons of total number N 0 are injected into the phantom to represent the incident beam.N 0 needs to be large enough so that the pixel variances in calculated reflectance images are less than 1%.After removing specularly reflected photons from N 0 by the Fresnel equation for unpolarized light, a total pathlength L a is obtained for each of remaining photons to be tracked in the turbid target before its first entry.The probability density of the {L a } set is an exponential function characterized by μ a .Individual photon tracking then starts by updating free pathlength L s,i + 1 by another exponential probability density function by μ s and propagation direction by p HG (cosθ) of single parameter g after i th scattering event [24,29].Tracking ends for a photon when either its accumulated pathlength L s exceeds L a or it exits the target.At each interface between air and the target or different zones of a heterogeneous target, the same Fresnel equation is used to determine the fate of tracked photon for reflection or refraction [23].The refractive index n of the target was set to a constant value of n = 1.40 based on the measured value at λ = 633nm for the silicone-based phantoms used in this study.The n values were determined using a coherent reflectance method [30] between 500 and 950nm and its variation in the range was found to be less than 4%.
The diffusely reflected photon density ρ r (x, y, 0_; P) was determined from photons exiting the target on the air side where P is the optical parameter set as P = {μ a , μ s , g} for homogeneous phantoms and P = (μ a1 , μ s1 , g 1 ; μ a2 , μ s2 , g 2 ) for 2-zone heterogeneous phantoms.Those photons registered for ρ r need to locate inside the FOV defined by the camera lens at z = −h and the solid angle of viewing ΔΩ(x, y, h) subtended by the lens of numerical aperture NA = sinα, as marked in Fig. 1(a).With ρ i and ρ r , we define the calculated reflectance image R c (x, y, -h; λ, P) as max, ( , , ; ) ( , , ; , ) , ( ) where max, ( ) max[ ( , ,0 ; )], In deriving Eq. ( 3), we assumed that the diffuse reflectance standard behaves as an ideal Lambertian reflector which is a very good approximation for small viewing angles [31,32].We verified the Lambertian behavior of the diffuse reflectance standards by measurement of the reflected light intensity recorded at the center region of the CCD sensor as a function of cosθ 0 for 40° ≤ θ 0 ≤ 80°, which was found to be linear with a slope proportional to the calibrated reflectance as expected.

Sample preparation and patient imaging
Two types of silicone based turbid suspensions with light and dark appearance have been prepared for making phantoms with TiO 2 powders (213581000, Acros Organics) and different concentrations of brown pigment powders (Pbr7, Kama Pigments) mixed in silicone polymer (RTV615A, MG Chemicals).The suspensions were stirred for 1 week to ensure homogeneity before casting and curing.Homogeneous turbid phantoms were prepared in cylindrical shaped molds of 40mm in outside diameter and 10mm in thickness.Multiple thin disk copies were made for each homogeneous phantom out of the same suspension with 18mm in diameter and 0.1 to 1mm in thickness for integrating sphere based measurements to determine their optical parameter set P [28].In this method, the collimated transmittance was first measured with four thin disks of different thickness to determine the attenuation coefficient μ t = μ a + μ s .This was followed by the measurement of diffuse reflectance and transmittance from one disk sample.Finally, a MC code developed for the integrating sphere method was used to determined P = (μ a , μ s , g) at different wavelengths [28].These data allow calibration of the reflectance imaging method for P measurement with the homogeneous phantoms.
Heterogeneous phantoms of melanoma have also been made with the same two suspensions.The light suspensions were first cast into a mold to prepare host phantoms, each of which has the same outside diameter and thickness as the homogenous phantoms and a cylindrically shaped cavity of 8mm in diameter and variable depth D m at the center.After curing, the dark suspension was filled into the cavity as the embedded region and pressed with a cover glass to ensure that the completed phantoms have flat surfaces and contain no visible air bubbles.Heterogeneous phantoms of three different D m were selected for this study to investigate the feasibility of reflectance imaging method for noninvasive determination of D m .Figure 1(b) shows one heterogeneous and two homogenous turbid phantoms.
To further examine the clinical potentials, we have acquired multispectral images from 4 patients of dysplastic melanocytic nevi in the Leo W. Jenkins Cancer Center under an IRB protocol approved by Brody school of Medicine, East Carolina University.For patient imaging, the interference filter wheels were replaced by a long-pass filter with the edge at 500nm and one of two liquid-crystal-tunable-filters (LCTF, VIS-07-20-STD-110, NIR-07-20-STD-110, CRI) of 10nm bandwidth to rapidly tune the wavelengths of incident beam during imaging with 15nm stepsizes.Figure 2(a) shows the patient imaging system with the incident beam angle of θ 0 = 45°.To reduce motion of FOV during imaging, a locking device was developed that consists of an optical window of 25mm in diameter at the end of illumination arm and a holding ring taped on the imaged skin site of patient as shown in Fig. 2(b).Before imaging, the optical window was coated with a very thin layer of clear ultrasound gel to ensure no-air contact with the patient skin.Each of illumination and imaging arms included a linear polarizer for acquiring cross-polarized reflection images of I r (x, y; λ) to reduce the light backscattered from superficial skin layers [33].

Effects of imaging parameters on the calculated reflectance image R c
The measured and calculated reflectance images, R m and R c , defined in Eqs. ( 1) and (3), differ by the planes of light or photon distributions.The former was collected with the camera lens at z = -h and the latter was calculated at z = 0_.Since the imaging process in this study is noncoherent in nature and far from the diffraction limit, we assumed that the point-spread function of camera lens is shift-invariant and distortion-free.Hence the difference between R m and R c is mainly due to the variation of solid angle ΔΩ(x, y, h) across FOV.While R c can be simulated at z = -h by eliminating the diffusely reflected photons exiting at (x, y, 0_) outside ΔΩ(x, y, h), such an approach reduces significantly the registered photon numbers because of the very small solid angle ratio of ΔΩ(x, y, h)/2π, about 10 −4 for h≥400mm, and requires large values of N 0 to reduce the pixel variances of R c .To improve simulation speed, we investigated the h dependence of R c (x, y)| h = R c (x, y, -h; λ, P) of homogeneous phantoms for both beam profiles of top-hat and grating.The effects of incident angle θ 0 from 0° to 45°, beam radius w, lens aperture diameter A and optical parameter P have also been investigated on R c with two sets of examples presented in Fig. 3.
In Fig. 3(a), linear profiles of R c (x, y)| h collected by a camera lens were plotted against the x-axis with an ideal top-hat beam incident on a homogeneous phantom with θ 0 = 0°.For h = 0.2mm, the profile is obviously limited by the lens aperture of A = 10mm.But as h exceeds 16mm, R c (x,0)| h profiles become alike and the peak pixel value ratio of R c (0,0)| h / R c (0, 0)| h = 0.2mm approaches to a constant of S = 1.60, as shown in Fig. 1(c).Similar results were obtained with grating beam profile as shown in Fig. 1(b) but S value increases to 1.86 as θ 0 increases to 45°.Collectively, these data show clearly that the image pattern is nearly independent of h for values larger than 100mm while the peak value of R c increases with h monotonically in Fig. 3(c) with a slope approaching to 0. Consequently, we define a scaling constant S = R c (x, y)| h /R c (x, y)| h = 0 for h = 400 mm for obtaining R c (x, y, -h; λ, P) by S•R c (x, y, 0_; λ, P).Figures 3(a) and 3(b) display scaling examples for both incident beam profiles.We also studied the dependence of the S on different P and found the relative change of S to be less than 5% in the concerned parameter ranges (data not shown).Based on these results, MC simulations of R c (x, y; 0_; λ, P) for this study were carried out at h = 0 and then scaled to h = 400mm by multiplying with S = 1.86 for θ 0 = 45° before comparison to R m .This method drastically reduces the required values of N 0 by 4 orders of magnitude to about 1x10 8 and simulation time.

Rapid simulation of the reflectance image R c by GPU implementation
The iMC code for R c (x, y, 0_; λ, P) has been implemented in a CUDA (8.0, Nvidia) enabled GPU version for parallel executions on a GPU board (Tesla K20, Nvidia), which is denoted as iMC-GPU.In this code, photons injected into the target were tracked with a singleinstruction-multiple-threads (SIMT) design.Compared to the single-instruction multiple-data (SIMD) paradigm, SIMT executes the programmatically equal code as different threads on GPU cores concurrently with data branching [34] using statistically independent sequences of random numbers.To improve performance, we grouped individual photon tracking to the SIMT based streaming architects.Before tracking starts, the total pathlength L a were obtained for all injected photons which were then sorted into groups of similar L a .The threads for tracking photons in the same group were submitted to and processed in the same streaming multiprocessor (SM) of the K20 boards that include 13 SMs of total 2496 cores.Each SM can handle up to 1920 threads or 10 threads per core.With this and other improvements over our earlier version [24], the current iMC-GPU code can have a speedup ranging between 30 to 50 relative to the single-CPU-core execution on a CPU of 3.4GHz (i7-3770, Intel).For N 0 = 1x10 8 , the wall-clock time for iMC-GPU simulation of one reflectance image of R c of a homogeneous phantom were 2.3s and 49s with P = (0.64mm −1 , 0.66mm −1 , 0.96) and (0.0088mm −1 , 1.29mm −1 , 0.36), respectively.

Study of the uniqueness of inverse solutions and the effect of the beam profile
The goal of reflectance imaging is to determine inversely the optical parameter set P and thickness of an embedded region from R m (x, y; λ).The inverse problem is defined here as an optimization problem for searching the minimum value δ min of an objective function δ for given λ and n = 1.40.The function δ given below is defined as the mean error between the measured and calculated images 2 , ( , ; ) ( , , ; , ) 1 ( , ) ( , ; ) where R c (x, y, -h; λ, P) = SR c (x, y, 0_; λ, P) as discussed in subsection 3.1 and the sum is taken over all pixels located at target surface of (x, y, 0_) with a weighing factor H and N p is the total number of pixels with H = 1.The factor H is set to 0 if the pixel value of R m is less than 5% of the maximum pixel value and 1 otherwise to exclude very low valued pixels that are in the peripheral region of the FOV because of low incident irradiance.For inverse solutions to converge or δ ≤ δ min , we found that it is critical to use the beam profile image S m (x, y; λ) measured with a diffuse reflectance standard for calculating.The use of S m instead of an analytical function and the same normalization scheme for R m and R c was essential for the inverse problems to have unique solutions.We first investigated the effect of incident beam profiles on the functional form of δ in the parameter space of (μ t , a, g) with a = μ s /μ t as the scattering albedo.Two profiles of top-hat and grating were investigated which are shown in Figs.1(c) and 1(d).At each value of μ t about 200 images of R c were calculated with different values of a = μ s /μ t and g to obtain a contour plot of δ against one measured image R m .
With different μ t we found that a unique solution or single minimum of δ exists for both beam profiles and Fig. 4 presents contour plots as examples for each profile.As can be seen from these plots that the grating profile produces a slightly steeper slope for δ to descend to δ min and the parameter values obtained there are closer to the values determined by the integrating sphere based method.We attribute the advantage of grating profile over the tophat to the larger variation in the reflectance values which makes δ more sensitive to the parameter change.For this reason, we chose the grating profile for subsequent study with phantoms while for patient imaging the top-hat profile was used for less sensitivity of R m to patient's motion.To examine the convergence of the inverse algorithm with grating profile of incident beam, we also investigated the landscape of δ in the parameter space for additional homogeneous phantoms at selected wavelengths.About 2000 reflectance images were calculated in about 12 steps along each axis of (μ t , a, g) to generate each two-dimensional contour map plot of δ shown in Fig. 5 against each R m measured with a phantom of light appearance.By these data, we confirmed that a unique minimum exists for δ for the surveyed ranges of the parameters for either light or dark phantoms.These results indicate that the objective function δ depends most sensitively on a, less on g and least on μ t .

Determination of P for homogeneous samples
For homogeneous phantoms, a nonlinear conjugate gradient (CG) algorithm [26,27] was adopted to obtain P = (μ a , μ s , g) from R m (x, y; λ) by iterated iMC-GPU simulations.In this algorithm, the solution P was obtained by iterating a trial solution vector of p k = (μ ak , μ sk , g k ) of the following recurrence until δ ≤ δ min 1 , where k = 0,1,… is the iteration index, α k > 0 is the stepsize of iteration and δ min was set to 0.1.The search direction for k th iteration d k is given by where g k = ∇δ k is the gradient of δ in the parameter space and β k is a scalar update parameter to ensure convergence in a nonlinear inverse problem.Initially introduced in [26], the following form of β k was employed in our study [27] 1 1 2 , ( ) y g d  g y  d y d y (8) where the superscript T indicates transpose and y k = g k + 1 − g k .After the direction d k is obtained, a strong Wolfe line search was performed to find α k that satisfies the conditions below [27], where 0 < ρ < σ < 1.Using the above CG algorithm, we iterated GPU-iMC simulations with σ = 2ρ = 0.5 after tests.The initial solution p 0 for each homogeneous phantom at the first wavelength was set to those determined by the integrating sphere method on the thin disk copies of the same suspension materials as the imaging phantoms.With P determined at the first wavelength, the set was used as p 0 for successive wavelengths.It took about 100 to 200 iMC-GPU simulations to converge at δ ≤ δ min for the first wavelength and about 30 to 60 for each of subsequent wavelengths.Local search was performed at selected wavelengths to confirm that the minimum value of δ is reached at the converged parameter values.Figure 6 presents the mean values and standard deviations of optical parameters of the two homogeneous phantoms at 25 wavelengths selected with the interference filters based on three measurements.It also contains the values of the same parameters determined from the thin disk copies of the same suspensions for each phantom by the integrating sphere method for comparison and validation of the reflectance imaging method.The results show good agreements on values of μ a between the two sets of samples which is mainly determined by the molecular absorption of the phantom materials.The increased discrepancies on μ s and g indicate different structures inside the two sets of phantom and thin disk samples on the scales of wavelength that could be due to the variations of embedded air micro-bubbles and actual distribution of TiO 2 particles.Still the wavelength dependence of P yields similar trends on μ s and reduced scattering coefficient μ s ' = (1-g)μ s .Taken together, data in Fig. 6 suggest that the approach of reflectance imaging and MC based inverse solutions is valid and can be used to determine optical parameter set P of turbid phantoms.Fig. 6.The wavelength dependence of P = (μ a , μ s , g) and μ s ' = μ s (1-g) of homogeneous phantoms and thin disk copies made of dark and light suspensions determined by the methods of reflectance imaging (RI) and integrating sphere (IS).The error bars represent standard deviations of three measurements and lines are for visual guide.

Noninvasive determination of embedded region thickness
A key potential application for the reflectance imaging approach described here is the noninvasive diagnosis and staging of superficial lesions such as dysplastic nevi and malignant melanoma.In current melanoma staging system the Breslow's thickness or tumor thickness plays a central role, which can only be determined through examination of biopsied tissues [35].To examine this possibility, we performed imaging measurement on the 2-zone heterogeneous phantoms that were made of the two suspensions identical to those for homogeneous phantoms.
Reflectance images R m were acquired at selected wavelengths from each heterogeneous phantom followed by obtaining calculated image R c with the optical parameters presented in Fig. 6 for the two zones and variable D for the zone 2 or the embedded region.Both zones in a heterogeneous phantom are of cylindrical shapes with the zone 1 or the host region being of 40mm in outside diameter and 10mm in thickness and zone 2 embedded in the center of the host with 8mm in diameter and variable thickness D. The thickness of the embedded zone 1, made of dark suspensions, was controlled to be D m = 0.6, 1.2 and 1.7 ± 0.1 mm for each of the three phantoms.

In vivo estimation of optical parameters of melanocytic nevi
Multispectral reflectance imaging and parameter estimation have been performed on 4 patients with melanocytic nevi to test the feasibility for pathological diagnosis.The patients were selected from those suspected of dysplastic nevi and/or melanoma and were prescribed for biopsy.All 4 patients were male Caucasians of fair skin complexion with ages ranging from 31 to 72 years old.Reflection images I r were acquired with p-polarizer in the illumination arm and s-polarizer in the imaging arm before biopsy and typically lasted about 25 minutes for each patient.After patient imaging, beam profile images I s were acquired with the diffuse reflectance standard of R s = 40%.The two sets of images were then processed to obtain reflectance images R m at 31 wavelengths between 500 and 950nm for each patient.Figure 8 presents the reflectance images at selected wavelengths from two patients.After imaging, each imaged lesion was excised surgically with a margin of about 3mm and the excised tissue was fixed in formalin and embedded in paraffin.Tissue slides were stained with hematoxylin and eosin (H&E) for histopathology grading after vertical sectioning.The severity of each lesion was graded on a standard scale of 0 to 8 based on the Clark's level of invasion [35]: 0 for benign lesion, 1 for mild dysplasia, 2 for moderate dysplasia, 3 for severe dysplasia, 4 for melanoma in situ, 5 for invasive melanoma within 1mm from junction, 6 for invasive melanoma between 1 and 2mm, 7 for invasive melanoma between 2 and 4mm and 8 for invasive melanoma beyond 4mm, all from the epidermisdermis junction.For this part of study, we modified a previous iMC code for two-layer skin modeling [20] by a voxel based photon tracking design to handle the irregular shape of the nevi.The melanocytic nevus was treated as the zone 3 of heterogeneous tissues in FOV, with the zone 1 and 2 designated as semi-infinite epidermis and dermis layers respectively.Consequently, the optical parameter set for an imaged skin target consists of 9 parameters given by P = {μ a1 , μ s1 , g 1 ; μ a2 , μ s2 , g 2 ; μ a3 , μ s3 , g 3 }, where the parameter sets of 1, 2 and 3 refer to epidermis, dermis and nevus, respectively.Furthermore, the thickness of the epidermis layer was set to D ep = 0.1mm and that of the nevus D n to a value consistent with the pathology grade to reduce the complexity of the inverse problems of determining P. Compared to the iMC code for modeling 2-zone phantoms, the modified iMC-voxel code took about 5-fold longer time to complete.
While the above optical model of skin nevus is reasonably accurate, GPU implementation of the iMC-voxel code requires substantial resources to complete reconstruction to achieve memory reduction and efficient thread management.Additional difficulty exists in applying the CG algorithm with iMC based forward calculations in the high dimensional space of 9 optical parameters.Instead, we simplified the inverse determination by minimizing the objective function δ sequentially in the order of 9 parameters arranged in P followed with a limited local search in the subspace of {μ a3 , μ s3 , g 3 } for the nevus.This allows estimation of the optical parameters for each imaged nevus with computational times drastically reduce to less than 5 hours per wavelength.In Fig. 9 we presented the estimated optical parameters of the nevus for 4 patients with a stepsize of 30nm in wavelengths.From these initial results, one can see that both absorption coefficient μ a3 and anisotropy factor g 3 in the near-infrared region appear to have larger differences between the patients of low grades and those of high grades.The increased μ a3 and decreased g 3 values for patients of grade 4 and 5 are consistent with the variation of the contrasts in the measured images of R m as shown in Fig. 8, where the nevi exhibit increased degree of darkness for λ > 700nm.

Discussion
Light penetrating into and emanating from a turbid target in backward scattering directions carries rich information for turbidity characterization.In this report, we present an approach of multispectral reflectance imaging for determining optical and depth parameters of turbid targets within the framework of the RT theory, where the sample can be "non-diffusive" or of low values of a and g.The new approach utilizes an incoherent light source for illumination of the target and can be used to determine μ a , μ s and g of a target from one reflectance image at a selected wavelength.Calibration of the measured and calculated reflectance images with the same beam profile image is critical for convergence of the inverse solution.Development of an accurate and rapid MC method is also important to form a well-posted inverse problem and solve in real times.The imaging approach for optical parameter determination has been validated against a well-established integrating sphere method on homogeneous phantoms.Currently, we are investigating combination of a diffusion solution for reflectance imaging and perturbation to replace iMC-GPU based calculations of R c for further reduction in computational time during iterations [24,36].In additional to accuracy, significant reduction in R c calculation time should enable the new approach of reflectance imaging for spectrophotometric determination of optical parameters within minutes.
One significant use of the reflectance imaging approach is to realize noncontact and noninvasive optical biopsy for diagnosis and staging of lesions.For this purpose, we performed measurement of thickness of an embedded region in heterogeneous phantoms and in vivo estimation of optical parameters of nevi in patients.It is widely recognized that the thickness of a penetrating pigmented nevus in the skin provides the key indicator for diagnosis of a malignant melanoma and patient's prognosis [35].With limited searches in parameter pace, our results presented in Figs.7 and 9 illustrate the possibility of thickness measurement and the optical parameter set of concerned tissues within the RT theory.The large differences in Fig. 9 between the values of g by the two patient groups of low and high lesion grades are particularly interesting in exhibiting the diagnosis value of anisotropy factor or scattered light distribution.It therefore calls for the necessity to apply RT theory instead of various diffusion approximations for taking the advantages of big data nature of imaging measurement.Needless to say, the patient data here are of very limited to carry statistical significance and additional studies are required to obtain conclusive results.While not proved yet, it is plausible that the reflectance imaging approach described here could be markedly improved from sequential 3D inverse searches for optical parameters of different zones in a heterogeneous target to a simultaneous inverse search in a 7D or 10D space of optical and thickness parameters.
where I max, s = max[I s (x, y; λ)−I sb (x, y)] is used to normalize both I r and I s .The measured images were further cropped and pixel binned to produce a 30mm FOV divided into a 101x101 grid for simulation described below.Figures1(b) to 1(d) present photos of reflectance standard, phantoms and examples of I r (x, y; λ) and I s (x, y; λ) before binning.

Fig. 1 .
Fig. 1.(a) Schematic of reflectance imaging and phantom: I: iris; S: stencil; L: collimating lens; WF: wavelength filter wheels; w: beam radius; CL: camera lens; A: lens aperture diameter; h: lens-target distance; (b) photos of (from left to right): two homogeneous phantoms of light and dark appearances, an 80% diffuse reflectance standard and a melanoma phantom; (c) reflection images of a homogeneous phantom of light appearance (I r : left) and an 80% diffuse reflectance standard (I s : right) at θ 0 = 45° and λ = 620nm without stencil for top-hat beam profile; (d) similar as (c) with stencil for grating beam profile.Bars = 10mm for photos in (b) and (c).

Fig. 2 .
Fig. 2. (a) The reflectance imaging unit for patient imaging; (b) a holding ring taped on a patient's back.

Fig. 4 .
Fig. 4. The contour plots of the objective function δ with μ t = 1.66mm −1 for the plots on the left and a = 0.760 for the plots on the right.The measure images R m (x, y; λ) were obtained from a dark homogeneous phantom with θ 0 = 45° and λ = 620nm of grating or top-hat beam profile.

Fig. 5 .
Fig. 5.The contour plots of δ obtained with the reflectance image measured from a light phantom sample with an incident beam of grating profile at θ 0 = 45° and λ = 500nm.
Figure 7 presents the dependence of δ on D for the three heterogeneous phantoms at different wavelengths.The thickness D of the zone 1 determined from R m yielded values close to the measured values for D m = 0.6 and 1.2mm while the difference increases for D m = 1.7mm.Furthermore, the values of D determined from R m at different wavelengths cluster around D m .These results demonstrate the new approach of reflectance imaging has the potential for noninvasive determination of dark zone thickness.The accuracy of D m determination could be increased by averaging D values over multiple wavelengths.

Fig. 7 .
Fig. 7.The objective function δ versus the thickness D of the embedded region in heterogeneous phantoms with reflectance images acquired at different wavelengths.The vertical lines and blue arrows indicate the measured thickness D m and uncertainty and color lines are for visual guide.

Fig. 8 .
Fig. 8. Examples of cross-polarized reflectance images acquired at wavelengths as marked on each image in the unit of nm.Top row: from patient B with pathology grade 1; bottom row: from patient D with pathology grade 4.

Fig. 9 .
Fig. 9.The wavelength dependence of optical parameters of nevus determined from the measured reflectance image Rm of 4 patients with melanocytic nevi.The pathology grades are noted in parentheses.