Fast 3D optical reconstruction in turbid media using spatially modulated light.

A method to perform fast 3-D optical reconstruction, based on structured light, in thick samples is demonstrated and experimentally validated. The experimental and reconstruction procedure, based on Finite Elements Method, used to reconstruct absorbing heterogeneities, with arbitrary arrangement in space, is discussed. In particular we demonstrated that a 2D sampling of the source Fourier plane is required to improve the imaging capability.

Diffuse Optical Tomography (DOT) is a promising technique to map molecular contrast, encoded into an optical variation, throughout a biological organism. The optical contrast can be due to anomalies of the absorption, scattering or fluorescence properties of the medium. For the case where a fluorescence signal is measured, DOT is generally referred to as Fluorescence Diffuse Optical Tomography (FDOT) or Fluorescence Mediated Tomography (FMT). DOT has been demonstrated for breast tumor detection [1, 2], brain functional imaging [3, 4] and small animal imaging [5, 6]. In the latter case, FMT is generally adopted even if the measurement of absorption anomalies finds direct application, such as the case of vascularization studies. Moreover, the reconstruction of the absorption and scattering properties of the sample allows a better quantification and localization of the fluorescence signal.

DOT, in all its implementations, has the advantage to be non-invasive, non-ionizing and highly selective with respect to the spectroscopic features of tissue constituents. However, the main drawback is represented by the low spatial resolution due to the strong multiple scattering of optical radiation within biological tissue in the visible/NIR range.

Among different approaches proposed to improve the reconstruction quality, dense sampling of the illumination and detection spaces is usually considered. This method generally consists of raster scanning a point light source over the sample. Then, for each injection point, a serial (e.g. by an optical fiber) or parallel (e.g. by a CCD) detection is carried out. In both cases, this leads to a huge data set, which requires long acquisition and computational times to solve the inverse problem. This represents a strong limitation for in vivo applicability.

In order to reduce the computational time of the inverse problem, it has been recently proposed to compress the measured data set by using a Fourier [7,8] or wavelets [9] encoding. In particular, analytical expressions have been developed which exploit symmetry properties of particular geometries such as a slab [7,10]. This approach has been successfully applied to data sets acquired by a raster scan of the source in the case of imaging of intrinsic absorption [11,12] and fluorescent concentration of contrast agents [13].

Recently, methods that use structured light, to perform fast wide field mapping of the optical properties of biological samples, have been proposed [14,15]. The key concept is that a highly scattering medium behaves as a low pass filter, leading to a limitation of the propagating Fourier components through the sample. The same concept has been also proven with time-resolved spatially modulated light [16, 17] and for fluorescence diffuse optical tomography by Joshi et al. [18].

As demonstrated by Schotland et al., all the relevant spatial information obtained by acquiring with point-like scanning sources, can be captured by recording only a few spatial Fourier components with wide field structured illumination [7]. Therefore, it has been proposed to perform optical tomography by sending structured illumination to generate data only within the low-spatial frequencies that are measurable [19, 20]. This approach allows to reduce both acquisition and reconstruction times, whilst retaining the same reconstruction quality. This novel imaging technique has been experimentally validated for the reconstruction, in reflectance geometry, of sub-surface heterogeneities [21] and in transmittance [22, 23]. In all these cases, spatial modulation in one direction and degenerate cylindrical inclusions, orientated parallel to the modulation direction, have been adopted. To the best of our knowledge, a tomographic approach based on structured light technique has not been applied on thick samples containing inclusions arbitrarily oriented with respect to the source modulation direction.

In this work we present and experimentally validate a fast 3D reconstruction method, based on structured light, in thick samples. In particular, we aim at reconstructing absorbing heterogeneities with arbitrary arrangement in space. To this end, we employ a number of source patterns forming a basis that spans the input plane, in particular a 2D cartesian sampling of the Fourier plane, i.e. a set of 2D sinusoids at different spatial frequencies. The reconstruction procedure is based on the solution of the diffusion equation by means of the Finite Elements Method (FEM) in order to make the technique suitable for arbitrary geometries. It is worth noting that the proposed reconstruction method is particularly suitable for small animal imaging, where the thickness of the animal is low enough to allow the detection of a reasonable number of spatial frequencies in transmittance geometry.

In this section a description of the set-up, shown in Fig. 1, is presented. The light is generated by a He:Ne laser (Melles Griot model 05LHP991) emitting at 633 nm with a power of about 10 mW. The light is first coupled by a lens (F220SMA-B, Thorlabs, USA; f=11.0 mm, NA=0.25) into an optical fiber (F2) (M41L01, Thorlabs, USA; Ø 600 μm, NA= 0.48) and then into a glass rod to create a flat top beam. By means of a lens (L1) and a total internal reflection prism the light is imaged on a Digital Micromirror Device (DMD) (Vialux GmbH, Germany). The DMD consists of an array of tiltable mirrors which can be independently controlled. Acting on the duty cycle of each mirror it is possible to spatially modulate the light profile with an 8 bit dynamic range. The pattern, to be projected to the sample, is uploaded and consists of bi-dimensional sinusoidal images. By an optical system (L2), the pattern generated by the DMD is magnified and projected on a parallelepiped sample. The area of illumination is 64 × 64 mm2 to match the dimensions of the phantom. By means of an objective lens (L3) (f = 50 mm, f# = 2.8, Nikon Co., Japan) the diffused light, exiting on the other side of the sample, is imaged on a low noise high sensitivity 16 bit CCD camera (Versarray 512, Princeton Instruments, USA) cooled to -70 C. The whole system is housed in a light-proof cabinet.

Fig. 1.

(a) Experimental set-up; (b) Phantom; (c) Source patterns along x and y from frequency k = 0 rad mm-1 to k = 0.40 rad mm-1.

The applied bi-dimensional sinusoidal modulated light has spatial frequency components k and k along x and y, respectively. For each spatial frequency, the measurement is repeated for 3 different spatial shifts. Figure 1(c) shows the first 4 patterns along x, y and their combinations. Experimental measurements were carried out on the homogeneous and inhomogeneous samples. The whole set-up is fully computer controlled in order to automatically change the frequency and shift of the sinusoidal input patterns and to acquire the corresponding image of the diffused light exiting the sample.

The sample consists of a parallelepiped (64×64×15 mm3) made of epoxy resin. In order to simulate mouse tissue-like optical parameters, TiO2 powder (Sigma-Aldrich, USA) as the scatterer and black toner (Infotec, Italy) as the absorber were added. The phantom has been characterized by a time-resolved spectroscopy system for diffusing media based on Time Correlated Single Photon Counting technique [24], giving the following optical parameters: μ = 0.012 mm-1 and μ′ = 0.81 mm-1.

In order to simulate a 3D geometrical structure, two non perpendicular cylindrical holes were drilled and filled with different inhomogeneities. The two cylinders have the same diameter of 2 mm and their axes are located at different depths. In particular the horizontal inclusion (A) is parallel to the x axis at the depth of 5 mm, while the oblique inclusion (B) is inclined by an angle of 62° with respect to the horizontal plane and located at depth of 11 mm. We will refer to these two inclusions as parallel (A) and oblique (B), respectively.

Two different approaches were exploited in order to change the optical parameters of the inhomogeneities. In the first case the holes were filled with a distilled water solution of Intralipid (IL) (Pharmacia, Italy) and ink (Rotring, Germany) in order to simulate scattering and absorption, respectively. Measurements were carried out by using the same scattering coefficient of the background and an absorption coefficient of μ = 0.24 mm-1 (contrast of 20:1 with respect to the background). In the following, we will refer to this as the first phantom. In the second case, two pieces of a totally absorbing tube was inserted into the oblique hole. The length of the two pieces is 36 mm and 11 mm, respectively, while the distance between them is 14 mm. Between the two pieces a liquid solution with the same optical parameters of the background was inserted. The parallel inclusion was also filled with a totally absorbing tube. In the following, we will refer to this as the second phantom. Finally, in order to obtain a homogeneous phantom the holes were filled with the liquid solution having the same optical properties as the background.

In this section a description of the data analysis and reconstruction algorithms is provided. Figure 2 schematically shows the steps leading to the optical reconstruction. The first step is to extract the modulated components of the light exiting the phantom and to remove the DC component due to the non-zero average of the light source intensity in the input plane. This demodulation was performed by the Phasor method. The basic idea is to sequentially inject at least three input signals at the same spatial frequencies but with different spatial shifts. By exploiting the images acquired at different spatial shifts, the modulated component of the output signal can be extracted [22].

Fig. 2.

Scheme of data analysis and reconstruction procedure.

Once the modulated component is extracted, a 2D Fourier Transform of the image is performed and a few spatial frequencies are selected. Although this is not a requirement of the method, we decided to use the same basis to encode the source and detector plane. Therefore, at this step, the dimension of the experimental data set is given by the product of the number of input frequency components along x and y (n = nin·nin) and the number of frequencies in the output plane along x and y (n = nout·nout). It is worth noting that the data set dimension is independent from the size of the original images.

The same procedure is repeated with the homogeneous and inhomogeneous phantoms. The difference of the Fourier amplitudes between homogeneous and inhomogenous cases, represent the input experimental data for the inverse problem.

The forward model corresponding to the experiment uses solutions to the diffusion equation in the slab domain Ω with inward photon flux J- on surface Γ1 and outward photon flux J+ on surface Γ2 that is detected by the camera

where is the diffusion coefficient, Φ the photon density, ζ a constant depending on refractive index, and ν the outward normal of the surface. Equations (1)–(3) can be summarised as a linear mapping Λ from the incoming to outgoing fluxes

where 𝓖 represents the Green’s operator that solves Eq. (1) with boundary condition [Eq. (2)] and 𝓜 is the measurement operator that restricts Φ to the detector surface Γ2.

We implement the above model using a finite element method (FEM) based on the public 25]. We use a regular grid of n × n × n voxel elements and n, nodes, with trilinear shape functions {u(x,y,z);k= 1…n}, resulting in a discrete model

where 𝖪 is the n × n system matrix of the FEM model, and 𝖬 is the nΓ × n measurement matrix with ones at the location of the nodes corresponding to surface Γ2 and zero elsewhere. In this work, the following numerical values were adopted: n = 64, n = 64, n = 16, n = 65536 and nΓ = 3969. The measurement on the CCD camera is modelled by projection from the surface Γ2

Here we take 𝖯 to be the identity operator since the camera is orthographically viewing the surface Γ2, although a more general method using ray-casting is possible for more general domains [9].

The applied inward currents J- are taken to be trigonometric functions J- = e-i where k = (k,k) is the wave number on Γ1 and ρ1 = (x1,y1) is the spatial position on Γ1. The measurements used in the reconstruction are the Fourier components of the projected currents y given by

To derive a linear reconstruction method, we consider the variation in the data due to a perturbation in absorption δμa by differentiating Eq. (7)

where Φ, and Φ+ are the forward and adjoint fields respectively, corresponding to incoming -i and Fourier measurement pattern e-i and is the variation in the FEM system matrix due to the perturbation in μa. We represent the image x = δμa in the same regular voxel grid as the FEM solutions, so that the change due to each voxel is given by the integral of the FEM shape functions and leads to the construction of the (n · n) × n Jacobian matrix 𝖨. The linear reconstruction problem is now represented

where g is the experimental data corresponding to the model data y. We use a Krylov solver method for Eq. (9) to give

where α is a regularisation parameter.

In Fig. 3 the reconstructed absorption of the first phantom, is presented. Slices at different depths, at step of 1 mm, are shown. Reconstruction was carried out by using 49 input spatial frequencies, due to the combination of 7 components along x and y (from k = 0 rad mm-1 to k = 0.59 rad mm-1). Conversely, the output plane was recreated by adopting a bi-dimensional encoding in the Fourier space, by using 9 spatial frequencies, from k = -0.4 rad mm-1 to k = 0.4 rad mm-1 at step of k = 0.1 rad mm-1, along x and y giving n = 81 measurements per source.

Fig. 3.

Reconstructed absorption at different depths in the case of bi-dimensional sampling of the input sources. Image scale in arbitrary units.

Concerning the resolution in the x-y plane (transversal section), the FWHM dimension of the inclusions is 6.5 mm (at z=11 mm) and 5.8 mm (at z=5 mm) for the inclusion closer to the detector (A) and source (B), respectively (Table 1). Both values are closer to each other even if the two inclusions are located at different depths. Concerning longitudinal direction (along z), the two inclusions can be clearly resolved, even if we observe a shift of the centers of the inclusions towards the border of the phantom. An analogous effect has been observed with optical reconstruction based on data acquired with a standard scanning based instrumentation and it can be attributed to the strong spatial variation of the sensitivity matrix. Concerning the dimension of the inclusions along z, they are clearly overestimated compared with the actual diameter of 2 mm. We believe this can be mainly attributed to the use of a single geometrical view and, therefore, the result can be improved by applying a multiple view approach, e.g. by rotating the sample in the case of a cylinder or more complex geometries.

Table 1.

Table summarizing the FWHM transversal dimension of the parallel (A) and oblique inclusion (B) by using 1D (horizontal and vertical) and 2D modulation.

	1D vertical	1D horizontal	2D
A	8.5 mm	6.4 mm	6.5 mm
B	7.1 mm	8.5 mm	5.8 mm

In order to test the reconstruction with a different geometrical arrangement of the inclusions with respect to the spatial modulation pattern, the first phantom was flipped. In this configuration, inclusion A is closer to the input and B to the output face, respectively. Figure 4 shows the reconstruction result. The number of input/output frequencies is the same as the previous case. Similar results were obtained both concerning longitudinal and transversal resolution.

Fig. 4.

Reconstructed absorption at different depths in the case of bi-dimensional sampling of the input sources with flipped sample. Image scale in arbitrary units.

It is worth comparing the reconstruction results obtained by bi-dimensional Fourier sampling of the sources with those obtained by modulating the source in one dimension only. Therefore, reconstruction was repeated by using the same input spatial frequencies but only along one direction (x or y). It is clear that if only modulation in one direction along the output plane is used, we can not reconstruct objects which are degenerate along that direction. Therefore, concerning the output plane, the same bi-dimensional Fourier encoding by using 9 spatial frequencies was adopted. In Fig. 5 the absorption optical reconstruction of the phantom is shown, by using a 1D spatial modulation of the light source along y [Fig. 5(a)] and along x [Fig. 5(b)], respectively. Slices at different depths, at step of 1 mm, are reported.

Fig. 5.

Reconstructed absorption at different depths in the case of horizontal (a) and vertical (b) source spatially modulation, respectively. Image scale in arbitrary units.

Concerning transversal direction, the FHWM dimension at z=11 mm of the horizontal inclusion (A), closer to the detector, is 6.4 mm [Fig. 5(a)] and 8.5 mm [Fig. 5(b)] in the case of horizontal and vertical spatial modulation, respectively. In the second case, the reconstructed inclusion results to be non uniform and larger. This can be attributed to the degeneracy of the modulation direction with respect to the orientation of the inclusion. Therefore, the spatial information encoded into the higher spatial frequencies is completely lost. On the contrary, in the case of parallel modulation, this leads to an increase of the spatial transversal resolution which is comparable with bi-dimensional sampling case (Table 1).

The oblique inclusion (B), closer to the source, represents an intermediate case because it is neither aligned along x nor y. For the modulation along x [Fig. 5(b)] the inclusion results to be homogeneously reconstructed along its length and the FHWM transversal dimension at z=5 mm is about 7.1 mm. In the case of modulation along y [Fig. 5(a)] the FHWM transversal dimension at z=5 mm is larger than 8.5 mm, even if a precise estimation is difficult due to the irregular structure. Therefore, the transversal direction results to be larger with the modulation along y. These results are in agreement with the previous considerations, since the oblique inclusion forms a smaller angle with the y-axis than the x-axis.

Finally, the reconstruction capability of the proposed technique has been tested on the second phantom containing two totally absorbing inclusions with a more complex shape. The horizontal inclusion is closer to the source and completely filled by the black rod, while the oblique inlcusion is closer to the detector and it consists of two separate parts, as described in Section 2.2. Figure 6 shows the absorption reconstruction at the different depths.

Fig. 6.

Reconstructed absorption at different depths in the case of bi-dimensional sampling of the input sources and a totally absorbing inclusions. Image scale in arbitrary units.

Concerning the oblique inclusion the reconstructed length of the two pieces is about 37 mm and 13 mm, respectively, and their distance is about 13 mm. This is in good agreement with the actual values reported in Section 2.2. Moreover, concerning the reconstructed transversal dimension, the FHWM dimension at z=5 mm of the horizontal inclusion (A) closer to the source, is 6.2 mm, while the FHWM dimension at z=11 mm of the oblique inclusion (B) closer to the detector, is 6.4 mm. Concerning longitudinal spatial resolution, analogous results to the previous cases (Fig. 3 and Fig. 4) can be found regarding the shift of the inclusion towards the borders and its enlargement along z.

Finally it is worth discussing the computational times for the reconstruction. With a 2.4 GHz and 4 Gbytes RAM personal computer the Jacobian calculation needed 114 s, while reconstruction 88 s. The Jacobian calculation depends on both mesh dimension and number of spatial frequencies while reconstruction only on the latter. Hence, the choice of the input/output spatial frequencies is critical for enhancing performance. In this work we operated empirically by reducing the number of input/output spatial frequencies as far as an appreciable change of the reconstruction results were observed. We believe that the optimal choice of the input and output spatial frequencies represents a fundamental step for the optimization of the technique in term of spatial resolution and acquisition/computational times.

Concerning the acquisition time, this is limited by the number of input sources, while it does not depend on the detection which is performed in parallel by a CCD. In our case, the amount of acquisitions is proportional to the number of spatial frequencies (49) multiplied by the number of spatial shifts (3). Hence, the total acquisition time was about 80 s with an integration time of 500 ms per pattern. It is worth noting that a wide field illumination technique, such as structured light approach, allows one to increase the light power impinging on the phantom while preserving a low power density; e.g. by assuming a light source of 100 mW the power density over 64×64 mm2 is about 2mW cm-2. Hence, by using a sample with optical parameters and thickness similar to the one used in this work, an exposure time, for each frame, of less than 50 ms can be easily achieved.

In this work, a fast 3D optical reconstruction, based on the use of spatially modulated light, of absorbing inclusions embedded into thick diffusing samples (15 mm) has been proposed and experimentally validated. The reconstruction algorithm has been implemented by using a Finite Elements Method, and therefore, it is suitable to reconstruct more complex geometries, such as the case of in vivo small animal imaging. In particular, the experimental and data analysis procedure used to extend the applicability of the technique to reconstruct heterogeneities, with arbitrary arrangement in space, has been discussed. We demonstrated that a bi-dimensional sampling of Fourier space of the input source is necessary to improve imaging capability.

Future developments will be devoted to the application of this technique to more complex geometries by using a multiple views tomographic approach and studying the optimal pattern according to the complex surface geometry. Moreover, time-resolved technique will be combined with structured light tomographic method in order to increase data set information and, e.g. to discriminate scattering and absorption contribution. Finally, the proposed technique will be applied to perform Fluorescence Molecular Tomography.