Imaging the optical properties of turbid media with single-pixel detection based on the Kubelka-Munk model

We present a diffuse optical imaging system with structured illumination and integrated detection based on the Kubelka-Munk light propagation model for the spatial characterization of scattering and absorption properties of turbid media. The proposed system is based on the application of single-pixel imaging techniques. Our strategy allows to retrieve images of the absorption and scattering properties of a turbid media slab by using integrating spheres with photodiodes as bucket de-tectors. We validate our idea by imaging the absorption and scattering coefﬁcients of a spatially heterogeneous phantom. optical properties, absorption coefﬁcient µ , scattering µ and anisotropy g , turbid ﬁeld optical imaging spectroscopy, optical used tissues, metabolic in laser

In many cases it is not sufficient to know the bulk absorption and scattering properties and, instead, the spatial distribution is required. In those cases, where 2D images of the optical properties are recovered, we speak of diffuse optical imaging (DOI). Unlike ballistic imaging techniques, that use non-scattered or slightly scattered light (like, for example, confocal microscopy, optical coherence tomography, time-gated imaging or multiphoton microscopy [2][3][4]), DOI techniques work in the diffuse regime and images are recovered using information provided by highly scattered photons. Due to multiple scattering events, the initial spatial information carried by photons is lost, hence it is usually necessary to resort to a light propagation model in order to image the spatial distribution of the absorption and scattering properties of the non-homogeneous turbid media [5].
Although light propagation in turbid media may be simulated by Monte Carlo techniques or mathematically modeled with a theory based on the radiative transfer equation (RTE), these models have a high computational cost or lack of analytic or numerical solutions without the use of approximations [1, 2,5]. A more simple one-dimensional model is provided by the two-flux Kubelka-Munk (KM) model [6][7][8]. This model is still widely used because of its relative simplicity, as it allows to express the scattering and absorption coefficients directly in terms of the measured reflection and transmission, and its acceptable prediction accuracy. However, conventional measurements in combination with the KM model are usually performed with integrating spheres, and it is not possible to resolve the spatial characteristics of the optical properties [9].
In this letter, we present a DOI system based on the KM model that uses structured illumination and integrated detection for the spatial characterization of scattering and absorption properties of non-homogeneous turbid media. The proposed system is based on the concept of single-pixel camera [10] and, in particular, on it's tolerance to the presence of diffusers between the object to be characterized and the bucket detector [11]. This allows to retrieve images of the optical properties of the turbid media slab by using integrating spheres with photodiodes as bucket detectors. The spatial resolution is achieved by projecting a set of microstructured light patterns onto the sample. Those patterns are codified by a digital micromirror device (DMD). Maps of total reflectance and transmittance, to apply the KM method, are generated from the photodiode signals by using single-pixel imaging (SPI) algorithms. This system is used to characterize the optical properties of a solid epoxy resin phantom with several scattering and absorption inhomogeneities and a slice of a cherry to distinguish two different types of tissue.
The KM model assumes that light propagation in turbid media can be modeled by two fluxes, which counterpropagate in the turbid media. The first flux (i) travels in the direction of the incident light and is decreased by absorption and scattering, while it is increased by backscattering of the counterpropagating flux (j). An analogous balance is made for the second or counterpropagating flux (j). The spatial evolution of those fluxes are described by the following system of differential equations: where S and K are the KM scattering and absorption coefficients. The KM model is based on several assumptions: the turbid media is isotropic and homogeneous (S and K are constant throughout the medium), light loss at the edges is neglected as infinite lateral extension is assumed, there is no reflection at the boundaries and illumination incident on the surface is completely diffuse [7,12,13].
For a finite turbid media slab of thickness d, integrating the KM equations yields the following expressions for the diffuse reflectance and transmittance: Where we defined the secondary constants a and b as Measurement of R d and T d yields the KM parameters S and K via the inversion of Eqs. (2-3) Likewise, S and K have been related to the light transport theory coefficients (absorption coefficient µ a and the reduced scattering coefficient µ s ≡ µ s (1 − g)) by the following empirical non-linear relations [14]: where S, K, µ a and µ s are expressed in units of cm −1 . The relations in Eq. (6) can be properly inverted to obtain the coefficients µ s and µ a from the calculated KM coefficients. Besides, we would like to remark that these relations are valid not only in the a diffusive (µ s µ a ) but also in the nondiffusive regime (µ s µ a ), and for nondiffuse incident light [14]. Thus, by choosing these relations we can overcome the model's limitations of diffuse irradiance and isotropic scattering.
Additionally, by applying Beer's law and measuring the collimated transmittance (T c ), the extinction coefficient (µ t = µ s + µ a ) is obtained via: where R F is the Fresnel coefficient of reflection for normal incidence, R F = [(n − 1) / (n + 1)] 2 ; and n is the relative index of the sample and surrounding media. This additional measurement of µ t yields in a straightforward way to obtain the values of the scattering coefficient µ s and the scattering anisotropy g: Thus, all three radiative transfer parameters (µ a , µ s and g) can be determined from experimental measurements of diffuse transmittance, T d , diffuse reflectance, R d , and collimated transmittance, T c , of the sample. The steps of this general process are represented in Fig. 1. The experimental setup is shown schematically in Fig. 2. A broad beam from a collimated deep red (λ = 660 nm) LED source (M660L4, Thorlabs) impinges onto a digital micromirror devide (DMD) (V-7001, ViALUX). The DMD implements a set of illumination patterns, which are functions of the 2-D Walsh-Hadamard (WH) basis. The generated light patterns are projected by a 4−f optical system, constituted by lenses L 1 and L 2 , onto the turbid media slab. The detection stage of this system consist of two integrating spheres (IS236A, Thorlabs), equipped with their corresponding photodiode (SM05PD1B, Thorlabs), collecting and measuring the diffuse reflected and transmitted light for every illumination pattern. By using a double integrating sphere arrangement we are capable of measuring simultaneously reflectance and transmittance [15]. The use of integrating spheres together with the SPI setup allows to obtain images of the distributions of the forward and backward scattered light in all angles, which are computationally reconstructed from the measured electric signal. This is feasible because SPI techniques allow using non-pixelated sensors such as the whole integrating sphere with a photodiode. Reference measurements are taken into account to obtain maps of absolute reflectance and transmittance. The transport coefficients are calculated from the previous reflectance and transmittance maps by performing KM calculations pixel-by-pixel. That is, spatial resolution in the reflectance and transmittance measurements is obtained by using structured light and single-pixel detection. But spatial maps of the transport coefficients are obtained by processing independently each pixel with the equations from the KM model. Double-integrating sphere corrections for multiple light exchange between the spheres are performed [16]. If only diffuse reflectance and diffuse transmittance measurements are performed, the spatial distribution of µ a and µ s are determined. However, the additional measurement of the collimated transmittance permits to determine the spatial distribution of the three transport parameters (µ a , µ s and g). Nevertheless, the acquisition of the collimated transmittance may limit the sample thickness to a few mean free paths (l ≡ µ −1 t ). This collimated transmittance is measured simultaneously with a third photodiode (DET36A, Thorlabs). The signal of all three photodiodes is amplified (PDA200C, Thorlabs) and digitalized by a data acquisition system (DAQ) (NI6251, National Instruments). The proyection of the patterns and the data acquisition is controlled by a homemade LabView software.
Patterns are modulated by 768 ×768 micromirrors of the DMD, each with a pixel pitch of 13.7 µm. Both lenses L 1 and L 2 (see Fig. 2) have the same focal length ( f 1 = f 2 = 100 mm), thus the optical 4−f system has unity magnification. A circular pupil is superimposed to the patterns to adjust the sampled area of the phantom to a size slightly smaller than the circular ports of the integrating spheres. In this way, we avoid artifact errors in the periphery of the measured images produced by light diffused too close to the port edges. This results in a circular imaged area of about 3.5 cm 2 . Walsh-Hadamard functions of order N = 128 are used, and therefore a set of 2×16384 patterns (each one consisting of 128×128 pixels) are projected on the sample. The factor 2 on the number of patterns arises because of the binary modulation nature of the DMD, as each Walsh-Hadamard function is codified into two patterns, one for the positive and one for the negative component. For the DMD working at 20 kHz, the acquisition time for all three images is about 1.6 s. It is worth mentioning that in SPI higher resolution could be obtained at the expense of a longer acquisition time.
In a first set of experiments, and in order to validate the system, the sample to be characterized (Fig. 3 (top)) is a slab of epoxy resin (∅ = 52.45 mm, thickness d = 7.55 mm, index of refraction n = 1.56) with TiO 2 nanoparticles (Titanium(IV) oxide, rutile nanopowder, < 100 nm, Sigma-Aldrich) as scattering agent. This phantom is made with a recipe similar to that in Ref. [17]. The bulk optical properties were measured by oblique-incidence reflectometry (OIR) [18] providing values of µ a of about 0.6 cm −1 and µ s of about 2.6 cm −1 . As depicted in Fig. 3 (bottom), two inclusions were inserted in the otherwise homogeneous slab. The first object is an absorbing and scattering heterogeneity, consisting in a hole of ∅ = 5.4 mm filled with an epoxy resin mixture with a higher nanoparticle concentration (µ a ∼ 0.1 cm −1 , µ s ∼ 10.5 cm −1 ). The second object is a 1.1 mm thick fragment of an absorptive neutral density filter which was placed inside the sample (serving as an absorption object).
Diffuse reflectance, diffuse transmittance and collimated transmittance images are shown in Fig. 4 (a)-(c), respectively. While calculated maps of the corresponding reduced scattering coefficient, absorption coefficient and scattering anisotropy are illustrated in Fig. 4 (d)-(f). For this sample, the reduced scattering coefficient map is shown instead of the scattering coefficient map in order to compare with the values provided by OIR. Each image is plotted within the minimum and maximum pixel value to highlight the spatial fluctuations of the measured parameters. In order to compare with the known bulk properties we will discuss about the properties of the different regions in the sample. A qualitative distinction of the absorption object can be made in the absorption coefficient map (Fig. 4 (d)). The absorption object presents a higher absorption coefficient value (µ a ≈ 0.65 cm −1 ) in contrast to the background. As expected because of the thinness of the absorbing fragment, no significant variation in the scattering properties is introduced by the absorbing object. For the scattering object, the average value of the reduced scattering coefficient is about 10.1 cm −1 . This result, together with an estimated value of the absorption coefficient of about 0.09 cm −1 , shows good agreement with the measured OIR values. The anisotropy parameter of about 0.31 also shows that the higher scattering nanoparticle concentration outcomes not only in a higher scattering coefficient value, but also in a more isotropic scattering. In the background turbid media, the absorption and scattering parameters are µ a ≈ 0.42 cm −1 and µ s ≈ 2.2 cm −1 . These results reproduce the qualitative behavior but don't match exactly the OIR values. However, as OIR relies on the diffusion theory (which requires µ s µ a ) no complete accuracy was expected for OIR values for the background turbid media.
In a second set of experiments, we imaged a 0.85 cm thick slice of a cherry. The area imaged in Fig. 5 (a)-(c) corresponds to the mesocarp layer (the central part) of the fruit. In this layer, and due to the different optical properties, we can distinguish to different tissues or parts. The first one is the parenchyma, which is the bulk or ground tissue of the fruit. This tissue is made of larger cells with a relative thin cell wall. The calculated optical properties maps (Fig. 5 (d)-(f)) show the values µ a ≈ 0.5 cm −1 , µ s ≈ 7.2 cm −1 and g ≈ 0.92 for this tissue. The second one is the vascular tissue responsible for the fluid and nutrient transport for the cells of the parenchyma. This tissue is characterized by a higher scattering coefficient (of about 9.1 cm −1 ), a lower absorption coefficient (about 0.26 cm −1 ) and also a lower anisotropy parameter (about 0.76).
In conclusion, we have presented a fast and inexpensive diffuse optical imaging system with a single-pixel detection approach that allows to measure maps of the spatially varying scattering coefficient, absorption coefficient and scattering anisotropy. Moreover, the system was tested with a heterogeneous phantom and latter used to characterize different tissues in a slice of an organic sample (a cherry). By using the relatively simple Kubelka-Munk model, it was not our aim to provide completely accurate characterization, rather than to separate absorption from scattering optical properties. In this sense, future work will be devoted to finding a more accurate model to retrieve scattering and absorption maps by only measuring reflectance. Besides, the KM model provides no depth information and therefore the determined properties are averaged over the whole depth of the sample. Depth information could be obtained via a multiple view approach. Additionally, the data acquisition process can be speed up by the application of compressive sensing techniques. The simplicity of the detection system makes possible to choose more complex sensors attached to the integrating spheres, such as fiber spectrometers for performing hyperspectral imaging.