Real-time photo-magnetic imaging

: We previously introduced a new high resolution diffuse optical imaging modality termed, photo-magnetic imaging (PMI). PMI irradiates the object under investigation with near-infrared light and monitors the variations of temperature using magnetic resonance thermometry (MRT). In this paper, we present a real-time PMI image reconstruction algorithm that uses analytic methods to solve the forward problem and assemble the Jacobian matrix much faster. The new algorithm is validated using real MRT measured temperature maps. In fact, it accelerates the reconstruction process by more than 250 times compared to a single iteration of the FEM-based algorithm, which opens the possibility for the real-time PMI.


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][2][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.

Forward: analytical solution to the combined diffusion and heat equation
The PMI forward problem is defined by a system of two equations [6,[13][14][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 bioheat thermal equation is used [10]: (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]: 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].

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]: 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.

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 −1 and 0.82 mm −1 , respectively. Meanwhile, the thermal properties are set to ρ = 1000 kg m −3 , c = 4200 J (kg C o ) −1 , and κ = 0.8.10 −3 W (mm C o ) −1 . 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 2 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. (1) Based on these fittings, we established that the distribution of the kernel J n at any node n(x 0 ,y 0 ) can be expressed as: where J s is the shape of the fitted kernel, obtained after normalizing the Jacobian computed on the node n(x 0 ,y 0 ). A is the amplitude of the kernel at the node position (x 0 ,y 0 ).

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 (6) 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.

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.   (7) 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 : 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).

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 0 ,y 0 ) is obtained using Eq. (5). Finally, the Jacobian, J, is used to update Δμ a using Eq. (4).

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 15meter 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. 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 (φ 0 ) 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 −1 ) 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. Figure 3 shows the first three successive MRT frames (200x200 pixels 2 ) 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.

Results
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].

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 2 ) 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.