Real-time photo-magnetic imaging

Farouk Nouizi; Hakan Erkol; Alex Luk; Mehmet B. Unlu; Gultekin Gulsen

doi:10.1364/BOE.7.003899

1. Introduction

The ability to probe several centimeters deep biological tissues using near-infrared (NIR) light paved the way for the development of diffuse optical tomography (DOT), which recovers the absorption and scattering properties of biological tissue [1–3]. However, the highly scattering nature of biological tissue drastically degrades the spatial resolution of DOT. In addition, performing measurements only at the boundary of the subjects makes DOT reconstruction a highly ill-posed inverse problem yielding non-unique solutions [4, 5]. To overcome the aforementioned DOT limitations, we previously developed a high resolution diffuse optical imaging modality, termed Photo-Magnetic Imaging (PMI) [6]. Technically, PMI is based on an alternative scheme compared to the conventional light detection at the boundaries. In fact, during the illumination of the subject with a continuous wave (CW) NIR light, PMI monitors the internal temperature variations induced by the optical absorption of the medium. Performing internal measurements reduces the under-determination of the inverse problem making its resolution more stable [7]. The PMI image reconstruction algorithm consists then in the minimization of the difference between the spatiotemporal simulated and measured temperature maps yielding high resolution absorption maps. During the data acquisition, the spatiotemporal distribution of temperature inside the medium is recorded using Magnetic Resonance Thermometry (MRT) [8]. The PMI simulated spatiotemporal temperature maps are obtained by solving the combined diffusion and Penne's bio-heat equations [9, 10]. Although, our previous FEM-based algorithm could provide high resolution absorption images, its computation time needed to be considerably accelerated for real time imaging.

First, using FEM requires meshing the studied medium and generating its basis functions, which is relatively a time consuming process. In addition, assembling the Jacobian matrix is the longest process in the resolution of the inverse problem. Unlike DOT, the use of the adjoint-method to assemble the Jacobian method is inadequate. This is because each MRT pixel is used as a measurement in PMI, which makes the number of measurements equal to the number of unknowns [2, 11]. Therefore, the Jacobian matrix is computed using the perturbation method, by sequentially solving the forward problem for a small variation in the absorption at each FEM node. Obviously, this solution is very time consuming and needed to be accelerated.

Here, we propose a new real-time image reconstruction algorithm that can provide a high resolution absorption map faster than the time required to acquire a set of MRT data [12]. This will be a key enabling feature for some certain PMI applications that require continuous feedback to enable orientation and fast localization of areas of interest such as tumor margin detection. To accelerate the resolution of the forward problem, we recently introduced a new analytical solution that directly provides the continuous spatiotemporal distribution of laser induced temperature [13]. This analytical solution is first obtained for the diffusion equation using an integral approach using the Robin boundary condition. Secondly, using the separation of variables, the heterogeneous heat equation is solved considering the separability of its temporal and spatial parts. For the inverse problem, a new analytical technique is used to assemble the Jacobian matrix. This empirical approach is based on a learning process from a series of extensive fittings performed on Jacobians computed using our FEM-based algorithm. Finally, the performance of our new real-time algorithm is evaluated using experimental data.

2. Methods

2.1 Forward: analytical solution to the combined diffusion and heat equation

The PMI forward problem is defined by a system of two equations [6, 13–15]. First, the CW form diffusion equation models the propagation of light in the medium [16]. Secondly, to model the propagation and dynamics of the temperature induced by the laser, the Penne's bio-heat thermal equation is used [10]:

\begin{array}{l} {- \nabla D (r) \nabla Φ (r) + μ_{a} (r) Φ (r) = S (r) \end{array}

{ρ c \frac{\partial T (r, t)}{\partial t} - \nabla k \nabla T (r, t) = Φ (r) μ_{a} (r)

where μ_a [mm⁻¹],

D (r) = 1 / 3 (μ_{a} + μ_{s}^{'})

[mm] and

μ_{s}^{'}

[mm⁻¹] are the absorption, diffusion and reduced scattering coefficients, respectively.

Φ (r)

[W mm⁻²] is the photon density at position r [mm] and S(r) is the source of light. The term ρ is the density [g mm⁻³], c, specific heat [J (g °C)⁻¹] and k the thermal conductivity [W (mm °C)⁻¹] of the medium. The right hand side of Eq. (2) represents the source of thermal energy induced by the laser [14]. The diffusion equation is solved by introducing a special Green’s function and the Robin boundary condition. Technically, utilizing the properties of the Dirac delta function and this boundary condition leads to an extensive photon density expression that will be used in the source term of Eq. (2) [17]. After implementing the heat source term, the separation of variables method is used to obtain the final temperature expression at any position (r,θ) and at any time t [13]:

T (r, θ, t) = T_{s} + \sum_{m = - \infty}^{\infty} \sum_{l = - \infty}^{\infty} \frac{ρ c}{κ λ_{l}^{2}} ω_{m, l} B_{m} (λ_{l} r) \cos (m θ) [1 - \exp (\frac{- κ}{ρ c} λ_{l}^{2} t)]

where T_s is the ambient temperature and B_m is the first kind Bessel function. The factors ω_m,l are derived from the solution of Eq. (2) and λ_l are obtained from the heat convection boundary condition [13].

2.2 Inverse problem

The PMI high resolution absorption images are obtained by minimizing the quadratic difference between the measured MRT temperatures and the simulated temperatures [6]:

Δ μ_{a} = {(J^{T} J + λ I)}^{- 1} J^{T} {(T_{p}^{M} - T_{p}^{s})}^{2}

where Δμ_a is the update to the absorption unknowns vector. I is the identity matrix, λ is the inversion regularization parameter and J is the Jacobian matrix.

2.3 Analytical sensitivity matrix

Considering that analytical solutions perform only on homogeneous media, we cannot assemble the Jacobian matrix using the perturbation method analytically. Therefore, we introduce a new method based on a learning process performed on Jacobians computed using our FEM-based algorithm. During this empirical approach, we vary the optical and thermal properties of the medium over a wide range and compute the individual Jacobians at each FEM node (J_n, n = 1,2,...N) in order to study the dependence of these Jacobians on the optical and thermal properties as well as their location. The simulations are performed on a 25 mm diameter synthetic circular phantom whose μ_a and $μ_{s}^{'}$ are set to 0.011 mm⁻¹ and 0.82 mm⁻¹, respectively. Meanwhile, the thermal properties are set to ρ = 1000 kg m⁻³, c = 4200 J (kg C^o)⁻¹, and κ = 0.8.10⁻³ W (mm C^o)⁻¹. The light source is positioned at the right side of the phantom. The Jacobian J_n at each node is mapped onto a 200x200 pixel² Cartesian grid. Unlike the banana shape Jacobian matrix obtained in DOT [11], the PMI Jacobian J_n consists in an exponential kernel centered at this particular node n. Afterwards, extensive series of fittings are performed on each Jacobian, obtained with a set of optical and thermal properties. Based on these fittings, we established that the distribution of the kernel J_n at any node n(x₀,y₀) can be expressed as:

J_{n} (x, y) = A (x_{0}, y_{0}) J_{s} (x - x_{0}, y - y_{0})

where J_s is the shape of the fitted kernel, obtained after normalizing the Jacobian computed on the node n(x₀,y₀). A is the amplitude of the kernel at the node position (x₀,y₀).

2.3.1. Shape of the kernel (J_s)

Based on the fitting series performed on the normalized Jacobians, we observed that the shape of the kernel is governed only by the optical and thermal properties and is totally independent from the spatial position within the medium or the strength of the laser source. The kernel fitted on the normalized Jacobian J_n /max(J_n), is given by

J_{s} (x, y) = e^{(- 3.12 μ_{a}^{0.58} - 2.41) {(\frac{κ}{ρ c})}^{0.27} \sqrt{{(x - x_{0})}^{2} + {(y - y_{0})}^{2}}} .

As mentioned above, this expression is validated for a large range of optical and thermal properties. This formula allows describing the Jacobian at any point (x, y) in the Cartesian grid and beats the spatial limitation imposed by the quality of the FEM mesh.

2.3.2. Amplitude of the kernel (A_n)

In addition to the dependences on the optical and thermal properties, the amplitude of the kernel, A, shows a strong dependence on both the spatial position within the medium and the power of the laser. Considering the high complexity of obtaining A analytically, we introduce an alternative time effective technique. This method allows obtaining A for any set of optical and thermal properties and is performed in two steps. First, the so called Total Jacobian, J_T, is computed analytically using the perturbation method. In fact, this is possible by varying the absorption coefficient of the whole medium rather than varying it at each pixel individually. Also, it is important to note that the Total Jacobian, J_T, can be obtained by the sum of all the individual Jacobians J_n. Based on this definition and Eq. (5), J_T can be expressed as the convolution of the normalized shape, J_s(x,y), and an amplitude function, A(x,y), as follows:

J_{T} (x, y) = \sum_{x = 0}^{N_{x}} \sum_{y = 0}^{N_{y}} J_{n} (x, y) = \sum_{x = 0}^{N_{x}} \sum_{y = 0}^{N_{y}} A (x_{0}, y_{0}) J_{s} (x - x_{0}, y - y_{0}) = A \otimes J_{s}

where N_x and N_y are the size of the Cartesian grid describing the domain containing the phantom, in the x- and y- directions, respectively. Therefore, Eq. (7) shows that A(x,y), the amplitude of the kernels at any position within the phantom is simply obtained by deconvolving the normalized kernel J_s from the Total Jacobian, J_T:

A = J_{T} \otimes^{- 1} J_{s} .

It is important to note that both J_s and J_T are noise free, therefore their deconvolution is performed using any deconvolution method. In all the following, the deconvolution is performed using the Lucy-Richardson method. After obtaining the amplitudes, A, the analytic Jacobians are simply assembled using Eq. (5).

2.4 Real-time PMI image reconstruction algorithm

The new analytic Jacobian assembly method is performed following these steps, Fig. 1. First, the shape of the kernel, J_s, is obtained using Eq. (6). Secondly, the Total Jacobian, J_T, is calculated analytically using the perturbation method by varying the absorption coefficient of the whole medium. The amplitude of the Jacobian at each pixel is obtained by simply deconvolving J_s from J_T. The Jacobian J at any point n(x₀,y₀) is obtained using Eq. (5). Finally, the Jacobian, J, is used to update Δμ_a using Eq. (4).

Fig. 1 Diagram and timeline of the real-time image reconstruction algorithm.

Download Full Size | PDF

2.5. Experimental studies

The experiments are performed using a homemade dedicated small animal PMI interface, Fig. 2. It consists of a customized small animal coil, with a side window for illumination. The coil is positioned at the center of a 3 Tesla MR system. The illumination of the phantom is performed using an 808 nm Focuslight FL-FCSE01 high power fiber-coupled laser. The laser, the signal generator for MRI synchronization and all the electronics are positioned in the control room due to the high magnetic fringe field. Light is delivered to the phantom via a 15-meter long optical fiber. To perform a uniform illumination over a 13.5 mm diameter spot, a 32.5 mm diameter Newport optics aspherical lens is used to collimate the laser beam. The laser power is kept under the maximum permissible exposure of skin imposed by the American National Standards Institute.

Fig. 2 a) Schematic of the PMI setup showing the optical instrumentation inside the MRI bore. b) Picture of the PMI setup.

Download Full Size | PDF

During the MRT acquisition, a high resolution phase map is acquired every 8 seconds using a gradient echo sequence with a repetition time (TR) and echo time (TE) set to 100 ms and 16 ms, respectively. First, in order to localize the laser position in the axial direction, a T1 weighted image is acquired. Afterwards, a dynamic imaging set consisting of multiple frames (8 seconds each) is acquired. Prior to turning the laser on, a baseline phase map (ϕ₀) is recorded. Then, the laser is turned on and the rest of the dynamic frames are acquired (ϕ_i, i = 1,2,3). The laser induced temperature is obtained during post processing and used in the PMI image reconstruction algorithm, Eq. (4) [14]. The experiment is performed on a 25 mm diameter cylindrical agarose phantom. The optical and thermal properties of the phantom are set to match the ones used for the numerical phantom used in section 2.3. A 4 mm diameter and approximately four times more absorbent (μ_a = 0.04 mm⁻¹) cylindrical inclusion is embedded 5 mm below the surface. The laser power and irradiated area are matched to those used in the previous experimental studies [14, 15]. To mimic a dynamic change in the absorption, the phantom is rotated at the beginning of each MRT frame by 40 degrees around its axis, Fig. 3.

Fig. 3 The columns show the results at the three successive MRT frames. The size and position of the inclusion are indicated with the blue dashed line circles. The inclusion is four times more absorbent than the background. The first row presents the high resolution MRT temperature maps (a-c) while the second row depicts the reconstructed absorption maps (d-f) obtained using our real-time algorithm at the three different MRT frames.

Download Full Size | PDF

3. Results

Figure 3 shows the first three successive MRT frames (200x200 pixels²) in the top row. It is observed that, temperature rises more within the inclusion because of its higher absorption [14]. The reconstructed absorption maps for each time point obtained with our real-time image reconstruction algorithm are presented in the bottom row of Fig. 3. Despite the noise present in the measured MRT maps, the inclusion is recovered successfully at all three frames and its location shows a very good agreement with the real position indicated by the blue dashed circles.

Our new algorithm allows the reconstruction of a high resolution absorption map every 7 seconds, which is approximately 250 times faster than one single-iteration of our FEM-based algorithm. This considerable acceleration of the reconstruction process allows us to obtain an image at the end of each MRI acquisition frame, making this method a real-time imaging technique. Moreover, the obtained results show that our new algorithm is robust to noise due to the higher amount of data used in the resolution of the inverse problem. In fact, this analytical-based reconstruction method uses the whole MRT data and does not require any mapping to mesh nodes as is the case for our previous FEM-based algorithm. However, slight artifacts at the boundary near the illumination site are visible in the images. Our analytical method models the laser as a point source while the laser beam is collimated and has a diameter of 13.5 mm. This mismatch in the source size yields these slight artifacts particularly near the illuminated surface. Nevertheless, these slight artifacts do not degrade the quantitative accuracy of our method which provides similar results compared to our previous studies using the FEM-based algorithm [14].

4. Conclusion

Previously, we introduced a new imaging technique that provides high resolution optical absorption maps by leveraging the high spatial resolution of MRI and the high sensitivity of optical imaging. The feasibility of PMI was demonstrated both with simulation and experimental data [6, 14]. However, the main limitation of this technique was its long computation time. In fact, the inverse problem of our FEM-based algorithm requires the assembly of the Jacobian matrix using the perturbation method, which requires heavy computational resources and a long computation time.

In this contribution, a real-time analytical reconstruction algorithm is presented for PMI. This algorithm is accelerated using analytical methods which allow the resolution of the forward problem without any medium meshing or system matrix assembly. Also, these analytical solutions provide a continuous solution over the medium, which allows a minimization between the simulated and the high resolution MRT temperature maps avoiding the measurements to FEM mesh mapping. This permits the utilization of the entire high resolution MRT measurements avoiding any data discardment. Essentially, we use these analytical solutions in the fast implementation of the Jacobian matrix which considerably accelerated the computation time. In fact, compared to our previous FEM-based algorithm, an acceleration of over 250 times is obtained. Technically, a PMI high resolution absorption image (200x200 pixel²) is reconstructed in approximately 7 seconds which is even shorter than the 8 seconds necessary to acquire an MRT frame. Considering the high performance of our new real-time image reconstruction algorithm, our future work would be its implementation for multi-wavelength small animal and breast imaging.

Acknowledgments

This research is supported in part by Fulbright awarded to Dr. Nouizi, NIH grants, F31CA171915-01A1, 1R21CA191389, R21EB013387, R01EB008716, P30CA062203, Bogazici University Research funding Grant No. BAP 7126, 11461 and TUBITAK Grant No. 112T253. Portions of this work were presented at the Biomed, OSA 2016, #OW4D.7.

References and links

1. J. C. Hebden, S. R. Arridge, and D. T. Delpy, “Optical imaging in medicine: I. Experimental techniques,” Phys. Med. Biol. 42(5), 825–840 (1997). [CrossRef] [PubMed]

2. S. R. Arridge and J. C. Hebden, “Optical imaging in medicine: II. Modelling and reconstruction,” Phys. Med. Biol. 42(5), 841–853 (1997). [CrossRef] [PubMed]

3. H. Dehghani, S. Srinivasan, B. W. Pogue, and A. Gibson, “Numerical modelling and image reconstruction in diffuse optical tomography,” Philos. Trans. R. Soc. A. 367, 3073–3093 (2009).

4. S. R. Arridge and W. R. B. Lionheart, “Nonuniqueness in diffusion-based optical tomography,” Opt. Lett. 23(11), 882–884 (1998). [CrossRef] [PubMed]

5. F. Nouizi, M. Torregrossa, R. Chabrier, and P. Poulet, “Improvement of absorption and scattering discrimination by selection of sensitive points on temporal profile in diffuse optical tomography,” Opt. Express 19(13), 12843–12854 (2011). [CrossRef] [PubMed]

6. Y. Lin, H. Gao, D. Thayer, A. L. Luk, and G. Gulsen, “Photo-magnetic imaging: resolving optical contrast at MRI resolution,” Phys. Med. Biol. 58(11), 3551–3562 (2013). [CrossRef] [PubMed]

7. E. J. Woo and J. K. Seo, “Magnetic resonance electrical impedance tomography (MREIT) for high-resolution conductivity imaging,” Physiol. Meas. 29(10), R1–R26 (2008). [CrossRef] [PubMed]

8. V. Rieke and K. Butts Pauly, “MR thermometry,” J. Magn. Reson. Imaging 27(2), 376–390 (2008). [CrossRef] [PubMed]

9. S. R. Arridge, “Photon-measurement density functions. Part I: Analytical forms,” Appl. Opt. 34(31), 7395–7409 (1995). [CrossRef] [PubMed]

10. E. H. Wissler, “Pennes’ 1948 paper revisited,” J. Appl. Physiol. 85(1), 35–41 (1998). [PubMed]

11. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15(2), R41–R93 (1999). [CrossRef]

12. F. Nouizi, H. Erkol, A. Luk, Y. Lin, and G. Gulsen, “Analytical Photo Magnetic Imaging,” Optical Tomography and Spectroscopy. Optical Society of America 7, OW4D (2016).

13. H. Erkol, F. Nouizi, A. Luk, M. B. Unlu, and G. Gulsen, “Comprehensive analytical model for CW laser induced heat in turbid media,” Opt. Express 23(24), 31069–31084 (2015). [CrossRef] [PubMed]

14. F. Nouizi, A. Luk, D. Thayer, Y. Lin, S. Ha, and G. Gulsen, “Experimental validation of a high-resolution diffuse optical imaging modality: photomagnetic imaging,” J. Biomed. Opt. 21(1), 016009 (2016). [CrossRef] [PubMed]

15. D. A. Thayer, Y. Lin, A. Luk, and G. Gulsen, “Laser-induced photo-thermal magnetic imaging,” Appl. Phys. Lett. 101(8), 083703 (2012). [CrossRef] [PubMed]

16. S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol. 37(7), 1531–1560 (1992). [CrossRef] [PubMed]

17. H. Erkol, F. Nouizi, M. B. Unlu, and G. Gulsen, “An extended analytical approach for diffuse optical imaging,” Phys. Med. Biol. 60(13), 5103–5121 (2015). [CrossRef] [PubMed]

Real-time photo-magnetic imaging

Abstract

1. Introduction