Improving vascular imaging with co-planar mutually guided photoacoustic and diffuse optical tomography: a simulation study

: Diﬀuse optical tomography (DOT) and photoacoustic tomography (PAT) are functional imaging modalities that provide absorption coeﬃcient 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 diﬀerent absorptions with a high spatial resolution, measuring the absolute value of optical absorptions using PAT is challenging due to unknown light ﬂuence 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 ﬂuence 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 eﬀectiveness of the proposed method.

To overcome the shortcoming of DOT, hybrid imaging approaches that combine DOT with other imaging techniques have been investigated [26][27][28][29][30][31][32][33]. Structural a priori information provided by complementary modalities such as ultrasound [26,27], x-ray [28,29] and MRI [30][31][32][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][42][43][44][45][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.

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.

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

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 TiO 2 . The absorption coefficient of the background is set to µ back a,true = 0.01 mm −1 . A blood vessel-like imaging target is embedded in the background with an absorption coefficient of µ vessel a,true = 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 µ s,true back = 1.35 mm −1 and µ s,true 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].

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.

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

PA data generation
To obtain PA data, first we calculate the initial pressure, p 0 , 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 = ΓF µ a , where, F is light fluence, µ 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/cm 3 , respectively. PA signals are recorded with a sampling rate of 25 MS/s.

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.

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.

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 µ s,initial =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 µ back a,initial = 0.005 mm −1 and µ vessel a,initial = 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 µ a,initial = 0.005 mm −1 , and µ s,initial = 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.

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

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. 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 Table 1.
where, µ a,recon is the absorption coefficient map in Figs. 4(b)−4(f), and µ a,true is the true absorption coefficient map in Fig. 4(a). n pixel , is the number of pixels in each region. According to the results, DOT-PAT and PAT-DOT-PAT had the best performance among others.  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., µ back a,true , µ s,true back). This explains the poor performance of PAT compensated image with the 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 µ s,initial =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.
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). 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][66][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 error µ 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.

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.

Disclosures
The authors declare no conflicts of interest.