Open Access
Add to BibSonomy
Abstract
Diffuse optical tomography (DOT) and photoacoustic tomography (PAT) are functional imaging modalities that provide absorption coefficient maps of the tissue. Spatial resolution of DOT is relatively low due to light scattering characteristics of the tissue. On the other hand, although PAT can resolve regions of different absorptions with a high spatial resolution, measuring the absolute value of optical absorptions using PAT is challenging due to unknown light fluence distribution in the tissue. Development of image guidance techniques using a priori information of imaging target structure has been shown to increase the accuracy of DOT. PAT is one such method that can be used as a complementary modality to serve as a guide for DOT image reconstruction. On the other hand, estimated fluence map provided by DOT can be used to quantitatively correct PAT images. In this study we introduce a mutually-guided imaging system for fast and simultaneous optical and photoacoustic measurements of tissue absorption map, where DOT is guided by the PAT image and vice versa. Using the obtained absorption map of the tissue, we then estimate the tissue scattering map. We conducted this study using a series of simulations on digital phantoms and demonstrated the effectiveness of the proposed method.
© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Diffuse optical tomography (DOT) [1–5] and photoacoustic tomography (PAT) [6–13] are used for non-invasive optical functional imaging. DOT provides spatial distribution of intrinsic tissue optical properties through a model-based reconstruction algorithm using near infrared (NIR) measurements along or near tissue boundaries [14–17]. In DOT, accurate quantification of optical properties for image reconstruction is challenging due to intense NIR photon scattering leading to low spatial resolution making the detection of small targets nearly impossible or significantly distorted [16,17]. PAT can resolve regions of different optical contrast with high spatial resolution, by ultrasound detection of thermoelastically generated photoacoustic pressure waves through short laser pulses illumination [18–24]. Although PAT excels at localizing absorbers because of tissue’s significantly lower acoustic scattering than optical scattering, the variability of fluence distribution of light, makes quantitative measurement of absorption in PAT challenging [25].
To overcome the shortcoming of DOT, hybrid imaging approaches that combine DOT with other imaging techniques have been investigated [26–33]. Structural a priori information provided by complementary modalities such as ultrasound [26,27], x-ray [28,29] and MRI [30–33] have been used to guide DOT image reconstruction [16,17]. However, each of these modalities have limitations: ultrasound is not sensitive to absorption, and therefore cannot adequately reconstruct the absorption map of the tissue, and, x-ray and MRI systems are difficult to make into an integrated co-planar system with DOT. PA imaging has potential to be used as a guide for DOT because it measures absorption and it can easily be integrated into a DOT system. Kumavor et al. [34] introduced a frequency-domain PA-guided DOT approach. In their work, locations of imaging targets at different depths were determined using PA imaging. The PA images were then used as a priori information to guide DOT reconstruction of absorption and scattering coefficients. However, the quality of results degraded with depth. Similarly, Wang et al. [35] proposed a dual imaging system consisting of separate PAT and DOT devices, imaging the same object. Each plane was first imaged by the PAT device, PA images were then used as a priori information to guide the DOT device to image the same plane. In this study, PAT and DOT measurements were not simultaneous as the sample moved from the PAT device to DOT device, which potentially introduced plane co-registration errors between the two. This becomes even more untenable when attempting functional imaging of in vivo subjects. On the other hand, there are methods developed to enable quantitative PAT including, deep learning algorithms [36,37], the eigenspectra method [38] and wavelength optimization method [39,40]. Fluence compensation [41–46], is another widely-used methodology to provide the estimation of light fluence distribution in the tissue in order to enable quantitative PAT. One approach to estimate the tissue fluence map is to use Beer-Lambert law to estimate the depth-dependent light attenuation according to the tissue optical properties [47,48] that are not readily accessible. Another approach is to use DOT [45,46,49].
Here, we introduce a fully integrated, co-planar mutually-guided imaging system for fast and simultaneous optical and photoacoustic measurements of tissue absorption map, where DOT is guided by the PAT image and vice versa. The goal is to overcome low resolution of DOT and improve quantification accuracy of PAT. This hybrid design allows for multi-plane imaging with inherent co-registration of PAT and DOT data. The proposed modality is fast, with ∼1Hz frame rate (hemodynamic response in human tissue is in the order of a few seconds [50,51]), which makes it suitable for functional imaging studies [4]. Using the proposed design, different hybrid imaging schemes are feasible; (i) DOT-compensated-PAT (DOT-PAT): uses the reconstructed fluence map from DOT to compensate the PAT image, (ii) PAT-guided-DOT (PAT-DOT): uses the structural information provided by PAT as a priori knowledge for DOT reconstruction, and (iii) PAT-guided-DOT-compensated-PAT (PAT-DOT-PAT): uses the fluence map obtained from PAT-guided-DOT to compensate the PAT image. We conducted these studies using a series of simulations on digital phantoms, and investigated the performance of each scheme. Using the obtained absorption map of the tissue, we then estimated the tissue scattering map.
In the following sections we first describe the concept of the proposed experimental setup in detail including its components and sequences with practical considerations being taken into account. Simulation setup is then defined based on this experimental setup, i.e. the location of light sources and detectors as well as ultrasound transducers. Then a 2D digital phantom which is used for this study is described, and the reconstruction algorithms are explained. Next section presents the results and investigates the performance of each scheme. Then we present the algorithm based on which the scattering map is estimated. Lastly, using the vertical scanning capability of this setup, we study the 3D imaging capability of the proposed modality.
2. Materials and methods
The simulation setup is described in section 2.1. The synthetic phantom properties are discussed in section 2.2. In section 2.3, the procedure of simulating PAT and DOT are explained and finally, image reconstruction method is discussed in section 2.4.
2.1 Simulation setup
The conceptual design for co-planar mutually-guided PAT/DOT is depicted in Fig. 1(a). Several practical considerations have been taken into account for the proposed setup. A tank of 50 mm diameter houses a gelatin-based phantom of the same diameter. Two lasers are used, one for PAT, and one for DOT imaging; with a more sophisticated control system, using a single laser for both modalities is also possible. Illumination wavelength of 800 nm is chosen due to its relatively high penetration depth [4,17,52]. The lasers generate pulses with the repetition rate of 100 Hz. Light from the PAT laser is directly coupled into a bundle of 16 optical fibers of 1mm diameter each. Following the method in [53], light from the DOT laser couples into a 1×16 optical switch with 16, 1mm optical fibers as its output. 48 holes are embedded in the circumference of the wall of the tank: 16, 2mm diameter holes for housing the optical fibers needed for DOT and PAT lasers illumination [shown in red in Fig. 1(a)], 16, 1mm holes for optical fibers needed for DOT detection [shown in gray in Fig. 1(a)], and 16, 1mm holes for ultrasound (US) transducers [shown in black in Fig. 1(a)]. DOT and PAT fibers, as well as DOT photodetectors and US transducers, are distributed equidistantly around the circumference of the phantom, as depicted in Fig. 1(a). These 3 sets alternate in the circumference of the tank wall with 22.5° between each component of its respective set. US transducers have an element size of 1 mm. The US signals are amplified using 16 amplifiers before being sent to the data acquisition (DAQ) system. The tank is rotated using a servo motor and a gear-bearing system, which allows vertical movement of the tank, while the phantom remains motionless. Details of the mechanical design of the rotating system have been given in our prior publication [54].
Fig. 1. Simulation setup and imaging protocol. (a) Schematic of the concept design of the co-planar mutually guided PAT and DOT system. Red and gray lines are illumination and detection fibers, respectively; black line is ultrasound transducer wire, (b) imaging protocol for DOT and PAT, (c) digital phantom of a blood vessel-like imaging target in gelatin.
The servo motor smoothly rotates the tank 22.5° before reversing direction and returning to the start position. One second is required for rotation in each direction. At the start of each rotation, the DOT laser sends one pulse to each of the DOT optical fibers, through the optical switch, in 160 ms. With each pulse, all 16 detection optical fibers collect backscattered light and send it to photodetectors. Then, the DOT laser is turned off, PAT laser is turned on, and servo motor starts moving for 680 ms. The PAT laser sends 68 pulses during this time. US transducers detect PA generated signals after each pulse. Finally, after 840 ms the servo motor stops and the DOT laser illumination and backscattered signal detection is repeated for the remaining160 ms. Figure 1(b) shows the pulse sequence protocol.
2.2 Synthetic phantom design
In this study, we use a 2D digital phantom as depicted in Fig. 1(c). The background 50 mm diameter circle represents a high scattering medium similar to a gelatin phantom containing light scattering material such as TiO2. The absorption coefficient of the background is set to $\mu _{a,true}^{back}$ = 0.01 mm−1. A blood vessel-like imaging target is embedded in the background with an absorption coefficient of $\mu _{a,true}^{vessel}$ = 0.5 mm−1 similar to that of blood [52]. The blood vessel-like imaging model is part of the k-wave software library [55]. The reduced scattering coefficient for background and imaging target are considered $\mu _{s,true}^{\prime}{back}$ = 1.35 mm−1 and $\mu _{s,true}^{\prime}{vessel}$=2.7 mm−1, respectively, following [52,56], to mimic a biological tissue such as breast. Since the reduced scattering coefficient is much greater than the absorption coefficient of the background, we are allowed to use the diffusion equation [17].
2.3 Simulated data
The generation process of the simulated data collected from US transducers and photodetectors within 1 second of working the mutually-guided PAT/DOT system are described here. Such process is detailed in the dashed green box in the flowchart shown in Fig. 2.
Fig. 2. Flowchart of data generation and image reconstruction in mutually-guided PAT/DOT simulations. Green-dashed box shows the procedures for generating PAT signals and DOT signals. DOT signals are only simulated at the beginning and end of servo motor rotation. Orange-dashed box shows reconstruction procedure. ϕ is the current angle of the tank. PAT-DOT, DOT-PA and PAT-DOT-PAT are described in Introduction Section.
2.3.1 DOT data generation at photodetectors
A circle of 5 cm radius is meshed with 96 ring, 57001 nodes and 112860 triangular elements and the phantom [shown in Fig. 1(c)] is mapped onto this mesh. 16, 1mm isotropic sources and 16, 1mm Gaussian detectors are defined according to the configuration described previously. Using TOAST++ [57], which solves the diffusion equation, we calculate the amount of diffused light that reaches each detection optical fiber with each DOT laser pulse. DOT requires separate measurements for each source-detector channel; that means a source fiber delivers a pulse while all detection fibers collect the diffused backscattered light. This procedure continues for all individual sources. It is important that throughout this procedure the source-detector configuration remains unchanged. To achieve this, we conduct DOT measurements at the beginning and end of each servo motor rotation [see the protocol depicted in Fig. 1(b)]. For each set of measurement we use 16 pulses. During DOT data acquisition, the system is fixed. Of the 100 pulses generated in one second, 32 are used for DOT and the remaining for PAT. The measurements are done in continuous-wave approach where detectors only measure the intensity of the light reaching them [4].
2.3.2 PA data generation
To obtain PA data, first we calculate the initial pressure, p0, generated by the laser illumination at each location inside the phantom. ${p_0}$ propagates until reaching the transducers. Initial pressure is calculated based on the equation, ${p_0} = \mathrm{\Gamma }F{\mu _a}$, where, F is light fluence, ${\mu _a}$ is absorption coefficient and Γ is Gruneisen parameter. In order to calculate ${p_0}$, the fluence distribution needs to be computed. For each laser pulse, there are 16 optical fibers around the phantom that illuminate the phantom at the same time. We used TOAST++ software to find fluence distribution inside the phantom. The output energy of each fiber is set to 1 mJ, and Γ set to 0.2. Propagation of PA waves is simulated by k-space method [58] using k-wave software [55]. PA signals recorded by each transducer are used for PAT image reconstruction. In the PA simulations, the grid size is set to 500×500 with 0.1 mm grid unit. A perfectly matched layer (PML) with thickness of 20 grid units is considered outside the main computational grid with the attenuation parameter, PMLAlpha, set to 2, and Courant-Friedrichs-Lewy (CFL) number set to 0.6. The maximum supported spatial frequency of this grid is 7.5 MHz. US transducers are simulated as 1mm planar transducers, each consisting of 10 grid units. Directionality of the transducers is achieved by spatially averaging over all grid points within each transducer. Such consideration will make the results closer to the experiment. Speed of sound and density of the phantom media are considered as that of water as 1500 m/s and 1 g/cm3, respectively. PA signals are recorded with a sampling rate of 25 MS/s.
2.4 Image reconstruction
2.4.1 PAT image reconstruction
The reconstruction procedure is shown in the orange dashed box in Fig. 2. It takes 1 second for the system to rotate 22.5° (1/16 of circle). During this time there are 68 laser pulses for PAT data collection. At each of the 68 time steps, the data from 16 channels of the DAQ are collected, resulting in 68×16 = 1088 view angles. We use the time-reversal reconstruction algorithm [59,60] (k-wave toolbox [55], Matlab) for image reconstruction.
2.4.2 Mask generation
We use MATLAB built-in functions to generate a mask from the reconstructed PA image described in section 2.4.1 using an opening-closing algorithm, described in [61]. The algorithm outlines the following: Initially, an opening operation is applied to the logarithmic PA image using the built-in imopen MATLAB function. The resulting image is subtracted from the original image followed by another opening operation and a closing operation; closing operation is performed using the built-in imclose MATLAB function. Both functions use an octogon structural element created by strel function. The window size for the first opening, second opening and the closing operations is set to 12, 3 and 9, respectively. Finally, a binary mask is generated by thresholding. The window size in the morphological operations and the threshold value to generate the mask are chosen empirically.
2.4.3 DOT, PAT-DOT, DOT-PAT and PAT-DOT-PAT reconstruction
We use TOAST++ for DOT image reconstruction. The inverse solver uses a model-based approach to iteratively update an initial distribution of absorption and scattering coefficients. Here, reduced scattering coefficient of the phantom is fixed and set to $\mu _{s,initial}^{\prime}$=1 mm−1 (a typical value for a biological tissue [52]) and only absorption coefficient is updated. Using the PAT mask generated in section 2.4.2 as an initial estimate for the DOT image reconstruction, two individual PAT-DOT images are obtained at the beginning and end of the servo motor full rotation. The two images are then averaged; resulting in a single PAT-DOT image for each 1 second of imaging. This is a continuous process. The initial values of absorption coefficient for the background and vessel, defined in the binary mask, are set to $\mu _{a,initial}^{back}$ = 0.005 mm−1 and $\mu _{a,initial}^{vessel}$ = 0.025 mm−1, respectively. We use the Gauss-Newton (GN) scheme [62] with Krylov solver for PAT-DOT reconstruction. The maximum iteration number and the convergence criterion are set to 100 and 10−6, respectively. Total variation method is used for regularization, with the regularization parameter β = 0.01 and hyperparameter τ = 10−6. Regularization prior is the binary PAT mask. Diffusivity factor and Perona-Malik threshold are set to 0.1 and 0.25, respectively. Similarly, the DOT image without PAT guide is obtained with initial values for absorption and reduced scattering coefficients, set to ${\mu _{a,initial}}$ = 0.005 mm−1, and $\mu _{s,initial}^{\prime}$ = 1 mm−1, respectively. Conjugate gradient (CG) method [63] and total variation regularization method are used for DOT reconstruction. The solver and its regularization parameters are the same as the ones used in PAT-DOT reconstruction; the regularization prior, diffusivity factor and Perona-Malik threshold are not used.
In addition to the absorption map, two fluence maps are also obtained from DOT as well as PAT-DOT. In DOT-PAT the fluence map obtained from DOT is used for the compensation of PAT image, and in PAT-DOT-PAT the fluence map obtained from PAT-DOT is used for the compensation of PAT image. The procedure of compensation is performed by dividing the PAT image by the given fluence map of a particular scheme. To assess the performance of different schemes, we compare the absorption maps obtained from DOT, PAT-DOT, DOT-PAT and PAT-DOT-PAT.
2.4.4 Scattering coefficient estimation
The obtained absorption map from DOT-PAT scheme described in section 2.4.3 is used in a DOT image reconstruction for the estimation of the scattering map. The scattering coefficient initial value is 1 mm−1, and the reconstruction parameters are same as the ones described in DOT reconstruction in section 2.4.3. In the reconstruction of the scattering map we realized that in situations where the absorption coefficient of a region is much greater than the background (>10), the scattering coefficient is not correctly estimated, probably because the gradient of the objective function in DOT reconstruction approaches zero. Therefore, in order to evaluate the capability of our system in estimating the scattering coefficient, we model five regions with diameters 1 mm, 1.5 mm and 2 mm and a high scattering coefficient of 7 mm−1 and the absorption coefficient same as the background, i.e., 0.01 mm−1. This tissue model could mimic the alteration in the tissue structure due to conditions such as thermal damage; thermal damage increases the tissue scattering coefficient several times [64].
3. Results and discussion
All simulations were done using a computer with 8 cores, an Intel i7-7700k CPU and 16GB of RAM. Figure 3(a) shows the ${p_0}$ distribution inside the phantom in a logarithmic scale calculated using the known absorption coefficient map of the phantom and the fluence map obtained from the diffusion equation. The bright regions along the perimeter denotes where the illumination optical fibers were located. Minimum pressure generated inside the vessels was in the order of ∼100 Pa which is close to the noise equivalent pressure level of commercial transducers; this pressure increased to up to several kPa closer to the boundaries of the phantom [see Fig. 3(a)]. Black arrows in Fig. 3(a) shows the locations of the US transducers for one of the steps.
Fig. 3. Initial pressure and PAT image of the synthetic phantom depicted in Fig.1 (c). (a) Initial pressure generated by a single laser pulse; the black arrows are the location of the transducers, (b) reconstructed PA image presented in logarithmic scale, (c) a binary mask extracted from PA image to guide the DOT image reconstruction.
Using the data from 1088 US transducers simultaneously (see section 2.4.1), a PAT image was reconstructed; the result is demonstrated in Fig. 3(b). Figure 3(c) shows the binary structural mask obtained from the logarithmic PAT image. This mask was used as a guide for the PAT-DOT reconstruction.
A side-by-side comparison between different PAT/DOT schemes is given in Figs. 4(a)−4(f). Figure 4(a) shows the true values of the phantom (true initial pressure image). Figure 4(b) shows the DOT image without PA guidance. The image provides an approximate location and extent of the imaging target. However, it is evident that in the absence of a priori knowledge, DOT could not produce an image with fine structures, and that the values of the absorption coefficients of different locations of the imaging target were underestimated. Figure 4(c) shows PAT-DOT image with a finer resolution than the one reconstructed by DOT alone. Figure 4(d) shows compensated PAT image using the fluence map obtained from a circular homogenous model with ${\mu _a}$=$\mu _{a,true}^{back}$=0.01 mm−1 and $\mu _s^{\prime}$=$\mu _{s,true}^{\prime}{back}$=1.35 mm−1. It can be seen that despite using the exact optical properties of the background, the compensated PAT image underestimated the absorption coefficient values. Figures 4(e) and 4(f) show DOT-PAT and PAT-DOT-PAT images, respectively, that are the compensated PAT images using fluence maps obtained from DOT and PA-DOT, respectively. It can be seen that the absorption maps in PAT-DOT and PAT-DOT-PAT images improved quantitatively compared to the one obtained from the homogenous model. Two horizontal intensity profiles [one shown with dashed yellow lines and another with dash-dotted yellow lines in Figs. 4(a)−4(f)] in the true absorption map, DOT, PAT-DOT, homogenous compensated PAT, DOT-PAT, PAT-DOT-PAT images are shown in Fig. 4(g) and Fig. 4(h), respectively. A comparison between the size of the same vessel [shown in dashed yellow lines in Figs. 4(a)−4(f)] defined as full width half maximum (FWHM) of the intensity signal profile have been performed, the results are given in Fig. 4(g). True size of the vessel was 0.7 mm (dashed black line) while the vessel size in the DOT image measured as 12 mm showing a large spread. The size of the vessel in the PA-DOT, DOT-PA and PAT-DOT-PAT images was measured as 0.8 mm, 0.8 mm and 1 mm, respectively. Our observations for DOT and PAT images are as follow [see Fig. 4(h)]: as oppose to the DOT that could not identify the spatial details of the phantom, and that in the PAT image compensated by the homogenous model, the peaks were visible but their values were not accurate, PAT-DOT, DOT-PAT, and PAT-DOT-PAT resolved four correct peaks at their correct locations.
Fig. 4. Quantitative analysis of different schemes of mutually-guided PAT/DOT system. (a) true absorption coefficients of the phantom, (b) DOT image without PAT guide, (c) PAT-DOT image, (d) PAT image compensated by the homogenous model, (e) DOT-PAT image, (f) PAT-DOT-PAT image, (g) intensity profile along the dashed yellow line, (h) intensity profile along the dash-dotted yellow line. PAT-DOT, DOT-PA and PAT-DOT-PAT were described in Introduction Section.
In order to assess the performance of different PAT/DOT schemes, the mean values and standard deviations of both vessel and background region in Figs. 4(b)−4(f) were calculated. Additionally, we defined a parameter, i.e. $erro{r_{{\mu _a}}}$, as root-mean-square-errors between each scheme’s image and the true absorption image [see Eq. (1)]. These quantitative results are given in Table 1.
Table 1. Mean, standard deviation and ${erro}{{r}_{{{\mu }_{a}}}}$ values of the reconstructed absorption maps obtained with PAT/DOT schemes.
Figure 5(a) shows the true fluence map inside the phantom. Figure 5(b) shows the reconstructed fluence map obtained from DOT. Although the absorption map obtained from DOT was not accurate [see Fig. 4(b)], its fluence map represents a close approximation of the true fluence map. Figure 5(c) shows the fluence map obtained from PAT-DOT. PAT-DOT-PAT performed better than DOT-PAT; that may be because PAT-DOT calculated the fluence map more accurately than DOT, this improvement however achieved at the cost of a longer computational time (∼15 minutes for DOT reconstruction versus ∼3 hours for PAT-DOT reconstruction). Figure 5(d) shows the fluence map obtained from the homogenous model that was significantly different from the true fluence map, despite using the true optical properties of the phantom (i.e., $\mu _{a,true}^{back}$, $\; \mu _{s,true}^{\prime}{back}$). This explains the poor performance of PAT compensated image with the homogenous model.
Fig. 5. Fluence distribution maps. (a) True fluence inside the phantom. Fluence maps obtained from, (b) DOT, (c) PAT-DOT, and (d) homogenous model.
Since DOT-PAT and PAT-DOT-PAT schemes use estimated fluence maps to compensate PAT images, the effect of noise in PAT image will be translated directly into the absorption maps obtained from DOT-PAT and PAT-DOT-PAT schemes. In PAT-DOT-PAT scheme, noise affects the quality of the extracted mask and lowers the image quality of the obtained absorption map. In DOT, noisy photodetector signal degrades the obtained absorption maps, however, the obtained fluence maps remained a close approximation to the true fluence map which suggests that DOT-PAT and PAT-DOT-PAT images are less prone to noise than DOT images. These are preliminary results and for confirmation, it requires a more rigorous study.
Using the absorption maps obtained from DOT-PAT scheme, through a DOT reconstruction, we obtained the phantom’s scattering map. The initial value of scattering coefficient was set to $\mu _{s,initial}^{\prime}$=1 mm−1. Figures 6(a) and 6(b) show the true and reconstructed scattering map of the phantom. The reconstruction was able to estimate the scattering coefficient of the background, i.e., 1.35 mm−1 (true value), as 1.17 ± 0.14 mm−1 (estimated value); the underestimation was probably due to the overestimation in the obtained absorption coefficient of the background as can be seen in Table 1.
Fig. 6. Estimation of differential scattering map. (a) True scattering coefficient map of the phantom without inclusion of the high scattering regions, (b) reconstructed scattering coefficient map of the phantom in (a) using the obtained absorption map from DOT-PAT, (c) true scattering coefficient map of the phantom including high scattering regions (depicted by white circles), (d) reconstructed scattering coefficient map of the phantom in (c) using the obtained absorption map from DOT-PAT, and (e) Differential scattering map obtained from subtraction of reconstructed scattering maps in (d) from (b); white arrows indicate the locations of the high scattering regions.
We studied different scenarios with the combination of different absorption and scattering coefficients and learned that the obtained scattering coefficient inside the vessel was not reliably estimated most likely due to the high absorption coefficient of the vessel; for instance, the vessel seen in Fig. 6(b) is an artifact and not the correct scattering map of the vessel. Figure 6(c) shows the same phantom as Fig. 6(a) with five scattering regions with diameters 1 mm, 1.5 mm, 1.5 mm, 2 mm and 2 mm and scattering and absorption coefficients of 7 mm−1 and 0.01 mm−1, respectively. Figure 6(d) shows the reconstructed scattering map of the true scattering map shown in Fig. 6(c). By subtracting the scattering maps in Fig. 6(d) from Fig. 6(b), the differential scattering map shown in Fig. 6(e) is obtained. It can be seen that the regions with high scattering coefficient are well detected in the differential scattering map.
Finally, we created a 3D synthetic phantom to investigate the vertical scanning capability of the setup. A complex 3D vessel-like phantom is made as shown in Fig. 7(a). The vessel was 33 mm in z direction. The extent of it in x and y directions was about 3 cm. The absorption coefficient of the vessel considered 0.5 mm−1. The vessel was placed inside a cylinder as background tissue with an absorption coefficient of 0.01 mm−1. The diameter of the cylinder was 5 cm similar to our previous 2D phantom and its height was 33 mm, the same as the vessel. The cylindrical background tissue is not shown in Fig. 7(a) for vessel to be more visible. The resolution of the true phantom was 0.1 mm. We scanned 67 planes in 0.5 mm intervals. For each plane we followed the procedure explained previously for DOT-PA reconstruction (see section 2.4.3). Figure 7(b) shows the resulting 3D reconstructed image by combining all 67 DOT-PAT images and after removing the background using the threshold value, 0.1. The mean and standard deviation values of the reconstructed absorption coefficient of the vessel was 0.326 ± 0.166 mm−1. We also reconstructed the 3D phantom with only DOT. The result was almost a uniform distribution with the value close to the background (i.e., vessels were not resolve).
Fig. 7. Reconstruction of 3D digital phantom. (a) 3D volume of true absorption coefficients of a vessel, and (b) DOT-PAT 3D volume of absorption coefficients of a vessel. The vertical resolution in these volumes is 0.5 mm.
DOT draws the attention of clinical researchers as a functional imaging method due to (i) its noninvasive nature, (ii) high temporal resolution and (iii) ability to measure the concentration of different tissue chromophores. Due to the highly scattering nature of the tissue, DOT has a poor spatial resolution especially at deep structures inside the tissue; this leads to inaccurate reconstruction of the absorption coefficient and consequently concentration of chromophores. One way to overcome this limitation, is to utilize the help of a higher spatial resolution imaging modality, to reconstruct the tissue structure within which DOT image reconstruction can be performed. PAT has a great potential for this purpose, mainly because of the fact that it can resolve regions with different absorption coefficients, can image deep into the tissue, and be implemented with a similar system configuration to DOT [65–67]. On the other hand, unlike DOT which directly measures the absorption coefficient, PAT measures the absorption coefficient that has been convolved with the fluence, making quantification of absorption coefficient challenging. Thus, with the aid of the fluence map obtained from DOT, the quantification accuracy of PAT can greatly be improved. Using the mutually-guided system we proposed, DOT is guided by the PAT image and vice versa, via three hybrid schemes: PAT-DOT, DOT-PAT, and PAT-DOT-PAT. Using the PAT-DOT scheme the vessel’s absorption reconstruction error was decreased by 29% and 33% compared to DOT, and PAT compensated image with homogenous model, respectively. The error was reduced by 41% and 44%, respectively when DOT-PAT used, and reduced by 45% and 48% when PAT-DOT-PAT used. DOT-PAT and PAT-DOT-PAT had the best performance among others in reconstructing the absorption map in terms of both quantitative accuracy and spatial resolution (see Fig. 4 and see the values of $erro{r_{{\mu _a}}}$ in Table 1). Using the absorption maps obtained from DOT-PAT, we then estimated the scattering maps of the phantom with and without inclusion of high scattering regions. The high scattering regions were detected in the obtained differential scattering map.
Some of the potential capabilities of the proposed system are as follows (of note that these claims are confirmed only when the system is implemented in practice): the system is able to measure both absorption and scattering maps of the tissue; the concurrent imaging capability of the system eliminates the potential error in co-registration of planes due to the movement; it can provide differential scattering map of the tissue to monitor any potential tissue structural damage; it can provide functional imaging capability due to the fast image acquisition (∼ 1 frame/s); the vertical scanning capability makes it especially suitable for long imaging objects such as breast.
To use the diffusion equation it is required the optical fibers of the system to be in contact with the imaging target or to be in a close proximity of it, meaning that the diameter of the construct has to be object specific. In some cases such condition may not be satisfied. One possible solution is a design which allows radial movement of fibers using a mechanism to push or pull fibers according to the size of the imaging target.
4. Conclusion
In this study, we simulated a co-planar mutually-guided PAT/DOT imaging modality for fast and concurrent DOT and PAT measurements where DOT reconstruction was guided by the structural information provided by the PAT image and in turn, the PAT image was compensated using the fluence map obtained from DOT reconstruction. We conducted a series of simulations on digital phantoms and assessed the performance of PAT-DOT, DOT-PAT and PAT-DOT-PAT imaging schemes on obtaining absorption coefficient maps. DOT-PAT and PAT-DOT-PAT produced the most accurate absorption maps in terms of both quantitative accuracy and spatial resolution. We also estimated the tissue scattering map and demonstrated the functional capability of our system to monitor tissue damage. With the fast growing application of DOT, as well as a great potential application of PAT, in breast and brain imaging, we believe the proposed metrology could greatly improve the diagnostic value of both DOT and PAT and facilitates the clinical translation of these promising modalities. Implementation and performance evaluation of the proposed system is the next step plan.
Funding
National Institutes of Health (R01EB027769-01, R01EB028661-01).
Disclosures
The authors declare no conflicts of interest.
References
1. D. R. Leff, O. J. Warren, L. C. Enfield, A. Gibson, T. Athanasiou, D. K. Patten, J. Hebden, G. Z. Yang, and A. Darzi, “Diffuse optical imaging of the healthy and diseased breast: A systematic review,” Breast Cancer Res. Treat. 108(1), 9–22 (2008). [CrossRef]
2. R. Choe, S. D. Konecky, A. Corlu, K. Lee, T. Durduran, D. R. Busch, S. Pathak, B. J. Czerniecki, J. Tchou, D. L. Fraker, A. DeMichele, B. Chance, S. R. Arridge, M. Schweiger, J. P. Culver, M. D. Schnall, M. E. Putt, M. A. Rosen, and A. G. Yodh, “Differentiation of benign and malignant breast tumors by in-vivo three-dimensional parallel-plate diffuse optical tomography,” J. Biomed. Opt. 14(2), 024020 (2009). [CrossRef]
3. B. J. Tromberg, A. Cerussi, N. Shah, M. Compton, A. Durkin, D. Hsiang, J. Butler, and R. Mehta, “Diffuse optics in breast cancer: Detecting tumors in pre-menopausal women and monitoring neoadjuvant chemotherapy,” Breast Cancer Res. 7(6), 279–285 (2005). [CrossRef]
4. F. Scholkmann, S. Kleiser, A. J. Metz, R. Zimmermann, J. Mata Pavia, U. Wolf, and M. Wolf, “A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology,” NeuroImage 85, 6–27 (2014). [CrossRef]
5. A. P. Gibson, T. Austin, N. L. Everdell, M. Schweiger, S. R. Arridge, J. H. Meek, J. S. Wyatt, D. T. Delpy, and J. C. Hebden, “Three-dimensional whole-head optical tomography of passive motor evoked responses in the neonate,” NeuroImage 30(2), 521–528 (2006). [CrossRef]
6. A. R. Mohammadi-Nejad, M. Mahmoudzadeh, M. S. Hassanpour, F. Wallois, O. Muzik, C. Papadelis, A. Hansen, H. Soltanian-Zadeh, J. Gelovani, and M. Nasiriavanaki, “Neonatal brain resting-state functional connectivity imaging modalities,” J. Photoacoust. 10, 1–19 (2018). [CrossRef]
7. L. De Liao, M. L. Li, H. Y. Lai, Y. Y. I. Shih, Y. C. Lo, S. Tsang, P. C. P. Chao, C. T. Lin, F. S. Jaw, and Y. Y. Chen, “Imaging brain hemodynamic changes during rat forepaw electrical stimulation using functional photoacoustic microscopy,” NeuroImage 52(2), 562–570 (2010). [CrossRef]
8. M. Gerling, Y. Zhao, S. Nania, K. J. Norberg, C. S. Verbeke, B. Englert, R. V. Kuiper, Å. Bergström, M. Hassan, A. Neesse, J. M. Löhr, and R. L. Heuchel, “Real-time assessment of tissue hypoxia In Vivo with combined photoacoustics and high-frequency ultrasound,” Theranostics 4(6), 604–613 (2014). [CrossRef]
9. J. Reber, M. Willershäuser, A. Karlas, K. Paul-Yuan, G. Diot, D. Franz, T. Fromme, S. V. Ovsepian, N. Bézière, E. Dubikovskaya, D. C. Karampinos, C. Holzapfel, H. Hauner, M. Klingenspor, and V. Ntziachristos, “Non-invasive Measurement of Brown Fat Metabolism Based on Optoacoustic Imaging of Hemoglobin Gradients,” Cell Metab. 27(3), 689–701.e4 (2018). [CrossRef]
10. C. Wood, K. Harutyunyan, D. R. T. Sampaio, M. Konopleva, and R. Bouchard, “Photoacoustic-based oxygen saturation assessment of murine femoral bone marrow in a preclinical model of leukemia,” J. Photoacoust. 14(2019), 31–36 (2019). [CrossRef]
11. M. R. Tomaszewski, M. Gehrung, J. Joseph, I. Quiros-Gonzalez, J. A. Disselhorst, and S. E. Bohndiek, “Oxygen-enhanced and dynamic contrast- enhanced optoacoustic tomography provide surrogate biomarkers of tumor vascular function, hypoxia, and necrosis,” Cancer Res. 78(20), 5980–5991 (2018). [CrossRef]
12. V. Ermolayev, X. L. Dean-Ben, S. Mandal, V. Ntziachristos, and D. Razansky, “Simultaneous visualization of tumour oxygenation, neovascularization and contrast agent perfusion by real-time three-dimensional optoacoustic tomography,” Eur. Radiol. 26(6), 1843–1851 (2016). [CrossRef]
13. K. Avanaki and J. Gelovani, “Ultrasound and multispectral photoacoustic systems and methods for brain and spinal cord imaging through acoustic windows,” U.S. patent US20200077974A1 (2019).
14. 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]
15. S. R. Arridge and J. C. Hebden, “Optical imaging in medicine II. Modelling and reconstruction-Phy in Med and Bio.pdf,” Phys. Med. Biol. 42, 841–853 (1997). [CrossRef]
16. S. Arridge, “TOPICAL REVIEW: Optical tomography in medical imaging,, Vol: 15, Issue: 2,” Inverse Probl. 15(2), R41–R93 (1999). [CrossRef]
17. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. Dimarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Process. Mag. 18(6), 57–75 (2001). [CrossRef]
18. L. V. Wang, Photoacoustic Imaging and Spectroscopy (CRC Press, 2009).
19. X. Wang, Y. Pang, G. Ku, X. Xie, G. Stoica, and L. V. Wang, “Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain,” Nat. Biotechnol. 21(7), 803–806 (2003). [CrossRef]
20. M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum. 77(4), 041101 (2006). [CrossRef]
21. L. V. Wang, “Photoacoustic Tomography: In Vivo Imaging from Organelles to Organs,” Science 335(6075), 1458–1462 (2012). [CrossRef]
22. E. Z. Zhang, J. G. Laufer, R. B. Pedley, and P. C. Beard, “In vivo high-resolution 3D photoacoustic imaging of superficial vascular anatomy,” Phys. Med. Biol. 54(4), 1035–1046 (2009). [CrossRef]
23. M. Nasiriavanaki, J. Xia, H. Wan, A. Q. Bauer, J. P. Culver, and L. V. Wang, “High-resolution photoacoustic tomography of resting-state functional connectivity in the mouse brain,” Proc. Natl. Acad. Sci. 111(1), 21–26 (2014). [CrossRef]
24. S. Mahmoodkalayeh, H. Z. Jooya, A. Hariri, Y. Zhou, Q. Xu, M. A. Ansari, and M. R. N. Avanaki, “Low Temperature-Mediated Enhancement of Photoacoustic Imaging Depth,” Sci. Rep. 8(1), 4873–4879 (2018). [CrossRef]
25. B. Cox, J. G. Laufer, S. R. Arridge, and P. C. Beard, “Quantitative spectroscopic photoacoustic imaging: a review,” J. Biomed. Opt. 17(6), 061202 (2012). [CrossRef]
26. N. Lv, N. He, Y. Wu, C. Xie, Y. Wang, Y. Kong, W. Wei, and P. Wu, “Effect of vascular haemoglobin concentrations on ultrasound-guided diffuse optical tomography in differentiating benign from malignant breast lesions,” Eur. Radiol. 24(11), 2848–2856 (2014). [CrossRef]
27. Q. Zhu, N. Chen, and S. H. Kurtzman, “Imaging tumor angiogenesis by use of combined near-infrared diffusive light and ultrasound,” Opt. Lett. 28(5), 337 (2003). [CrossRef]
28. J. Selb, S. A. Carp, G. Boverman, E. L. Miller, D. H. Brooks, R. H. Moore, D. B. Kopans, and D. A. Boas, “Combined Optical and X-ray Tomosynthesis Breast Imaging,” Radiology 258(1), 89–97 (2011). [CrossRef]
29. A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chaves, J. Stott, Q. Zhang, T. Wu, M. Chorlton, R. H. Moore, D. B. Kopans, and D. A. Boas, “Tomographic optical breast imaging guided by three-dimensional mammography,” Appl. Opt. 42(25), 5181 (2003). [CrossRef]
30. B. A. Brooksby, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Near-Infrared (NIR) Tomography Breast Image Reconstruction With A Priori Structural Information From MRI: Algorithm Development for Reconstructing Heterogeneities,” IEEE J. Sel. Top. Quantum Electron. 9(2), 199–209 (2003). [CrossRef]
31. S. Merritt, F. Bevilacqua, A. J. Durkin, D. J. Cuccia, R. Lanning, B. J. Tromberg, G. Gulsen, H. Yu, J. Wang, and O. Nalcioglu, “Coregistration of diffuse optical spectroscopy and magnetic resonance imaging in a rat tumor model,” Appl. Opt. 42(16), 2951 (2003). [CrossRef]
32. R. J. Cooper, M. Caffini, J. Dubb, Q. Fang, A. Custo, D. Tsuzuki, B. Fischl, W. Wells, I. Dan, and D. A. Boas, “Validating atlas-guided DOT: A comparison of diffuse optical tomography informed by atlas and subject-specific anatomies,” NeuroImage 62(3), 1999–2006 (2012). [CrossRef]
33. X. Intes, J. Yu, A. G. Yodh, and B. Chance, “Development and evaluation of a multi-wavelength multi-channel time resolved optical instrument for NIR/MRI mammography co-registration,” Proc. IEEE Annu. Northeast Bioeng. Conf. NEBEC 2002-Janua, 91–92 (2002).
34. P. D. Kumavor, C. Xu, A. Aguirre, J. Gamelin, Y. Ardeshirpour, B. Tavakoli, S. Zanganeh, U. Alqasemi, Y. Yang, and Q. Zhu, “Target detection and quantification using a hybrid hand-held diffuse optical tomography and photoacoustic tomography system,” J. Biomed. Opt. 16(4), 049802 (2011). [CrossRef]
35. Y. Wang, J. Li, T. Lu, L. Zhang, Z. Zhou, H. Zhao, and F. Gao, “Combined diffuse optical tomography and photoacoustic tomography for enhanced functional imaging of small animals: a methodological study on phantoms,” Appl. Opt. 56(2), 303 (2017). [CrossRef]
36. J. Gröhl, T. Kirchner, T. Adler, and L. Maier-Hein, “Estimation of blood oxygenation with learned spectral decoloring for quantitative photoacoustic imaging (LSD-qPAI),” arXiv: 1902.05839, (2019).
37. C. Cai, K. Deng, C. Ma, and J. Luo, “End-to-end deep neural network for optical inversion in quantitative photoacoustic imaging,” Opt. Lett. 43(12), 2752 (2018). [CrossRef]
38. S. Tzoumas, A. Nunes, I. Olefir, S. Stangl, P. Symvoulidis, S. Glasl, C. Bayer, G. Multhoff, and V. Ntziachristos, “Eigenspectra optoacoustic tomography achieves quantitative blood oxygenation imaging deep in tissues,” Nat. Commun. 7(1), 12121 (2016). [CrossRef]
39. R. Hochuli, L. An, P. C. Beard, and B. T. Cox, “Estimating blood oxygenation from photoacoustic images: can a simple linear spectroscopic inversion ever work?” J. Biomed. Opt. 24(12), 1 (2019). [CrossRef]
40. N. Dana, T. Sowers, A. Karpiouk, D. Vanderlaan, and S. Emelianov, “Optimization of dual-wavelength intravascular photoacoustic imaging of atherosclerotic plaques using Monte Carlo optical modeling,” J. Biomed. Opt. 22(10), 1 (2017). [CrossRef]
41. B. T. Cox, S. R. Arridge, and P. C. Beard, “Estimating chromophore distributions from multiwavelength photoacoustic images,” J. Opt. Soc. Am. A 26(2), 443 (2009). [CrossRef]
42. E. Malone, S. Powell, B. T. Cox, and S. Arridge, “Reconstruction-classification method for quantitative photoacoustic tomography,” J. Biomed. Opt. 20(12), 126004 (2015). [CrossRef]
43. R. Hochuli, S. Powell, S. Arridge, and B. Cox, “Quantitative photoacoustic tomography using forward and adjoint Monte Carlo models of radiance,” J. Biomed. Opt. 21(12), 126004 (2016). [CrossRef]
44. M. A. Naser, D. R. T. Sampaio, N. M. Munoz, C. A. Wood, T. M. Mitcham, W. Stefan, K. V. Sokolov, T. Z. Pavan, R. Avritscher, and R. R. Bouchard, “Improved Photoacoustic-Based Oxygen Saturation Estimation With SNR-Regularized Local Fluence Correction,” IEEE Trans. Med. Imaging 38(2), 561–571 (2019). [CrossRef]
45. L. Ulrich, L. Ahnen, H. G. Akarçay, S. S. Majos, M. Jaeger, K. G. Held, M. Wolf, and M. Frenz, “Spectral correction for handheld optoacoustic imaging by means of near-infrared optical tomography in reflection mode,” J. Biophotonics 12(1), e201800112–16 (2019). [CrossRef]
46. R. J. Zemp, “Quantitative photoacoustic tomography with multiple optical sources,” Appl. Opt. 49(18), 3566–3572 (2010). [CrossRef]
47. J. C. Ranasinghesagara and R. J. Zemp, “Combined photoacoustic and oblique-incidence diffuse reflectance system for quantitative photoacoustic imaging in turbid media,” J. Biomed. Opt. 15(4), 046016 (2010). [CrossRef]
48. W. C. Vogt, X. Zhou, R. Andriani, K. A. Wear, T. J. Pfefer, and B. S. Garra, “Photoacoustic oximetry imaging performance evaluation using dynamic blood flow phantoms with tunable oxygen saturation,” Biomed. Opt. Express 10(2), 449 (2019). [CrossRef]
49. L. Yin, Q. Wang, Q. Zhang, and H. Jiang, “Tomographic imaging of absolute optical absorption coefficient in turbid media using combined photoacoustic and diffusing light measurements,” Opt. Lett. 32(17), 2556 (2007). [CrossRef]
50. T. J. Huppert, R. D. Hoge, S. G. Diamond, M. A. Franceschini, and D. A. Boas, “A temporal comparison of BOLD, ASL, and NIRS hemodynamic responses to motor stimuli in adult humans,” NeuroImage 29(2), 368–382 (2006). [CrossRef]
51. C. M. Carpenter, R. Rakow-Penner, S. Jiang, B. W. Pogue, G. H. Glover, and K. D. Paulsen, “Monitoring of hemodynamic changes induced in the healthy breast through inspired gas stimuli with MR-guided diffuse optical imaging,” Med. Phys. 37(4), 1638–1646 (2010). [CrossRef]
52. S. L. Jacques, “Optical Properties of Biological Tissues: A Review,” Phys. Med. Biol. 58(11), R37–R61 (2013). [CrossRef]
53. M. Flexman, “Dynamic Digital Optical Tomography for Cancer Imaging and Therapy Monitoring,” Columbia University (2012).
54. M. Zafar, K. Kratkiewicz, R. Manwar, and M. Avanaki, “Development of Low-Cost Fast Photoacoustic Computed Tomography: System Characterization and Phantom Study,” Appl. Sci. 9(3), 374 (2019). [CrossRef]
55. B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. 15(2), 021314 (2010). [CrossRef]
56. J. L. Sandell and T. C. Zhu, “A review of in-vivo optical properties of human tissues and its impact on PDT,” J. Biophotonics 4(11-12), 787 (2011). [CrossRef]
57. M. Schweiger and S. Arridge, “The Toast++ software suite for forward and inverse modeling in optical tomography,” J. Biomed. Opt. 19(4), 040801 (2014). [CrossRef]
58. B. T. Cox, J. G. Laufer, K. P. Kostli, and P. C. Beard, “Experimental validation of photoacoustic k-Space propagation models,” Photons Plus Ultrasound Imaging Sens. 5320, 238 (2004). [CrossRef]
59. B. E. Treeby, E. Z. Zhang, and B. T. Cox, “Photoacoustic tomography in absorbing acoustic media using time reversal,” Inverse Probl. 26(11), 115003 (2010). [CrossRef]
60. P. Omidi, M. Zafar, M. Mozaffarzadeh, A. Hariri, X. Haung, M. Orooji, and M. Nasiriavanaki, “A novel dictionary-based image reconstruction for photoacoustic computed tomography,” Appl. Sci. 8(9), 1570 (2018). [CrossRef]
61. R. Gonzalez, R. Woods, and S. Eddins, Digital Image Processing Using MATLAB, 2nd Edition (2009).
62. M. Schweiger, S. R. Arridge, and I. Nissilä, “Gauss-Newton method for image reconstruction in diffuse optical tomography,” Phys. Med. Biol. 50(10), 2365–2386 (2005). [CrossRef]
63. S. R. Arridge and M. Schweiger, “A gradient-based optimisation scheme for optical tomography,” Opt. Express 2(6), 213 (1998). [CrossRef]
64. V. K. Nagarajan and B. Yu, “Monitoring of tissue optical properties during thermal coagulation of ex vivo tissues,” Lasers Surg. Med. 48(7), 686–694 (2016). [CrossRef]
65. S. Mahmoodkalayeh, M. A. Ansari, and M. Nasiriavanaki, “A new illumination scheme for photoacoustic computed tomography,” Proc. SPIE 10494, 104946T (2018).
66. M. Zarei, S. Mahmoodkalayeh, M. A. Ansari, K. Kratkiewicz, R. Manwar, and M. R. Nasiri Avanaki, “Simultaneous photoacoustic tomography guided diffuse optical tomography: a numerical study,” Proc. SPIE 10878, 213 (2019). [CrossRef]
67. M. Zarei, M. A. Ansari, S. Mahmoodkalayeh, and M. R. Nasiriavanaki, “Hybrid co-planar diffuse optical tomography and photoacoustic imaging,” Proc. SPIE 10494, 245 (2018). [CrossRef]
