Open Access
Share with Facebook
Tweet This
Post on reddit
Share with LinkedIn
Add to CiteULike
Add to Mendeley
Add to BibSonomyAbstract
Photoacoustic imaging (PAI) is emerging as a potentially powerful imaging tool with multiple applications. Image reconstruction for PAI has been relatively limited because of limited or no modeling of light delivery to deep tissues. This work demonstrates a numerical approach to quantitative photoacoustic image reconstruction that minimizes depth and spectrally derived artifacts. We present the first time-domain quantitative photoacoustic image reconstruction algorithm that models optical sources through acoustic data to create quantitative images of absorption coefficients. We demonstrate quantitative accuracy of less than 5% error in large 3 cm diameter 2D geometries with multiple targets and within 22% error in the largest size quantitative photoacoustic studies to date (6cm diameter). We extend the algorithm to spectral data, reconstructing 6 varying chromophores to within 17% of the true values. This quantitiative PA tomography method was able to improve considerably on filtered-back projection from the standpoint of image quality, absolute, and relative quantification in all our simulation geometries. We characterize the effects of time step size, initial guess, and source configuration on final accuracy. This work could help to generate accurate quantitative images from both endogenous absorbers and exogenous photoacoustic dyes in both preclinical and clinical work, thereby increasing the information content obtained especially from deep-tissue photoacoustic imaging studies.
© 2016 Optical Society of America
1. Introduction
Photoacoustic imaging (PAI) is an exciting, emerging imaging modality that has the ability to combine the high intrinsic contrast of optical imaging with the excellent spatial resolution of ultrasound imaging [1–4]. The technique is non-invasive and can be performed in real time when used in pulse-echo mode like traditional ultrasound. It is an excellent avenue for interrogating angiogenesis and blood content in highly aggressive tumors, and shows great potential for molecular imaging [5]. PAI and Photoacoustic Tomography (PAT) have been developing for many years in both clinical and preclinical regimes with applications in skin [6] and breast cancer [7,8], lymph node mapping [9–11], and endoscopy [12,13]. These technique relies on absorption of light to create contrast, a phenomenon that is governed by both tissue optical properties and light fluence. Light is typically attenuated approximately an order of magnitude per cm in tissue. However, conventional PA techniques ignore this, causing both spatial and spectral distortion artifacts during image reconstruction [14]. These artifacts can cause deep objects to appear dimly or fail to be seen at all and cause incorrect quantification of tissue chromophores such as hemoglobin or oxygen saturation. If PA is to be realized as a deep tissue imaging technique and move from the preclinical imaging domain, these distortions must be addressed. To this end, we propose a model-based tomographic method to map tissue absorption coefficients based on acoustic data. We demonstrate accurate, quantitative images of absorption and spectrally derived chromophores at unprecedented depths.
Various algorithms have been developed for overcoming depth-related artifacts, but these techniques have limitations. Most algorithms use Filtered Back Projection (FBP), which generates qualitative approximations of absorbed energy. While they can produce high contrast images, they tend to break down when realistic absorbers (<10x contrast) are located greater than 3cm from a light source due to exponential depth attenuation. In an effort to improve on these by recovering quantitative images of tissue absorption, many experimental methods have been developed that either update the FBP reconstruction or scale it based on fluence distributions. Previously, Cox et al. demonstrated an algorithm that corrects a photoacoustic image based on light fluence by iteratively fitting for the absorption coefficient, μa, using a FBP-generated initial pressure map as a starting point [15]. Banerjee et al. were able to recover maps of μa from the initial pressure map and boundary pressure measurements [16]. Razansky and colleagues modify their FBP algorithm to include an analytic calculation of fluence distribution where μa is approximated based on a calibration technique [17,18]. Bauer et al. showed improved reconstruction using hybrid diffuse optical tomography and PAT in large phantoms [19]. These methods are memory efficient and reasonably effective in some geometries, but they rely on the FBP reconstruction as part of the algorithm. Therefore, while these methods improve quantitative accuracy, they have difficulty addressing artifacts resulting from light attenuation in deep tissue.
In an effort to reconstruct a quantitative domain based solely on boundary data, other groups have developed model-based imaging techniques. Bu et al. developed a method to boost contrast to noise ratio in low fluence regions using an iterative model to recover accurate but only qualitative images [20]. Yuan and Jiang were able to reconstruct optical energy density in a 1cm3 volume with a known background concentration within 10% of the true value using a frequency domain algorithm [21]. Time domain reconstruction methods in Quantitiative Photoacoustic Tomography (qPAT) were introduced by Yao and Jiang, showing accurate reconstructions of absorbed optical energy in both simulations and phantoms with multiple targets in 3 cm radius circles [22]. While these studies have shown potential to quantify photoacoustic images, none of them have been able to generate quantitative images of tissue absorption coefficient based on optical source excitation. Finally, Cox et al. were able to demonstrate model-based recovery of μa, μs’, or spectrally derived chromophores from multi-wavelength qPAT images in 4x8mm rectangular regions [23]. Here, optical sources are used to reconstruct optical properties, but algorithm performance in regions deeper than 1 cm is unknown. Deep imaging is a requirement for clinical translation. We aim to address all of these issues in an effort to make quantitative photoacoustic imaging clinically possible.
In this work, we present a model-based method for recovering absolute absorption coefficient maps from deep tissues (>4cm) with no prior knowledge of optical properties, initial pressure distributions, or structure within the domain. A diffusion and k-space acoustic model is used to reconstruct absolute absorption coefficients based purely on acoustic boundary data. We demonstrate the ability to recover optical properties and spectral chromophore concentrations simultaneously in 2 circular domains (3 and 6cm diameter) and characterize the algorithm with respect to initial guess and time resolution. The method is able to mitigate the problems of both spectral and spatial distortion that affect FBP and could be used to accurately recover absolute absorber concentrations, making it especially suitable to deep tissue clinical PA imaging.
2. Methods
Photoacoustic imaging requires that a short, nanosecond pulse of laser light illuminates a sample. As the light is absorbed, local heating causes thermoelastic expansion and generates pressure waves. These waves propagate through the tissue and can be measured non-invasively by ultrasound detectors at the tissue boundary. The forward model used in this work consists of two separate stages to solve for the acoustic pressure data given the optical properties, optical excitation and acoustic propagation. The inverse problem uses multiple forward models to minimize the difference between measured data and simulation to recover the optical properties.
2.1 Simulation geometry
We show improvement over conventional reconstruction methods in five different scenarios that show the effects of depth attenuation, error in initial guess, time resolution, source configuration, and spectral decomposition. The simulation geometry is a 30mm diameter circle with four distributed optical sources arranged at 90 degree increments. Simulated acoustic point detectors are arranged every 3 degrees surrounding the entire medium as shown in Fig. 1(a). The mesh consisted of 17,869 nodes with a node spacing of 0.2mm. Optical fluence was generated using a background μa of 0.01mm−1 and μs’ of 1mm−1 (Fig. 1(b)). The targets were 2mm in diameter and placed at depths of 1cm and 1.5cm into the medium. Their optical properties were μa = 0.1 mm−1 and μs’ = 1 mm−1. We assumed a constant Gruneisen parameter of 0.1 for the whole domain to set the initial pressure after light modeling (Fig. 1(c)). The initial pressure distribution seeded acoustic wave propagation and pressure data was recorded at each transducer location with 20ns between data points. For optical reconstruction, data was down sampled ten times for computational savings. Speed of sound was constant at 1500m/s. For reconstruction, the mesh was set to uniform values equal to the background concentration for both μa and μs’. The reconstruction was then allowed to run until the error change between iterations was less than 1% or 40 iterations were reached. For comparison, we also used the reconstructed this same data using an FBP reconstruction algorithm with the same time resolution that the data was generated [24]. Before back projection, data was passed through a frequency domain Ram-Lak filter. Results were evaluated based on quantitative accuracy of the reconstruction. Unless otherwise stated, these properties were used for all simulations.
Fig. 1 Sources and detectors surround model domain with three inclusions (a). A diffusion model calculates light fluence (b) and initial pressure maps (c) based on optical sources (blue dots). Sound waves propagate to transducer locations (green dots) to generate forward data (d). A newton method minimizes model-data mismatch to converge on a quantitative solution (e). Normalized FBP (f) is shown shown for comparison.
2.2 Forward model
Optical excitation is roughly five orders of magnitude faster than sound propagation. As a result, it is acceptable to use a steady state solution to the light diffusion equation,
to solve for optical fluence distribution Φ. Here, a source S with frequency ω describes light fluence through the turbid media. The diffusion coefficient,is a function of the absorption coefficient, μa, and the reduced scattering coefficient, μs’. The optical diffusion equation is an acceptable simplification of the radiative transport equation under the assumptions that optical radiance is only linearly isotropic, and that the rate of change of the photon flux is much lower than the photon collision frequency. The first assumption holds for regions far (>1 mm) from light sources and tissue boundaries, and is acceptable in small animal imaging where the tissue volume is sufficiently large. The second assumption is true because scattering occurs much more frequently than absorption in tissue in the Near-Infrared (NIR) regime [25]. Since the goal of this work is to extend the utility of quantitative photoacoustic imaging in deep tissues, the diffusion equation is a valid choice as it simplifies the model with negligible penalty to accuracy.Optical excitation sources are considered isotropic and placed one scattering distance, approximately 1mm, into the medium to satisfy assumptions. Light propagation is modeled using the open source software package, Nirfast, over a discretized triangular mesh using the finite element method [26]. The conversion of light absorption to sound energy is represented by,
which describes the initial pressure, p0, as a function of the Gruneisen parameter, Γ, the fluence, Φ, and the absorption coefficient, μa. The Gruneisen parameter is a measure of the conversion efficiency from light absorption to pressure and was assumed to be a constant 0.1 throughout the medium. After the optical fluence distribution is obtained for every node in the mesh, this equation is used to solve for the initial pressure at each location. From that point, the acoustic wave can be described by,where H(r,t) is the amount of thermal energy converted to pressure at a position r and time t. Cp is the isobaric specific heat, B is the isobaric volume expansion coefficient, and vs is the acoustic speed in the medium [27,28]. This equation is solved for all space over sufficient time using the Kwave Toolbox, an open source software platform for modeling photoacoustic signals [24,29]. All simulations assume a lossless medium with uniform speed of sound. Simulated ultrasound transducers are placed outside the medium and pressure data is measured as the pressure at that location versus time. Uniform frequency response was assumed in the transducers.2.3 Image reconstruction
The final image is reconstructed using the time-dependent Levenburg-Marquardt formulation to iteratively minimize the difference between modeled data and true data with 1% noise added [30]. The Jacobian matrix, J, is calculated for each time step using the adjoint method [26] and is a measure of the sensitivity of each transducer to each location in the mesh. The full, time-dependent Jacobian was structured as in Tichauer et al [30]. Here, the model-data mismatch is broken into time bins, and a separate Jacobian is built for each of those time bins. Time referencing of the signal is done by convolving each line of the Jacobian with a window function specific to the time step and distance from each detector. Once J has been calculated for every time step, the optical properties of the medium are updated using the equation,
Where c is a vector of the parameters being reconstructed on, I is the identity matrix, and δ is the model-data-mismatch. Lambda is a regularization term included to limit the influence of noise on the data by reducing the ill-posedness of the image reconstruction problem. Optical properties are updated after every iteration and if the change between iteration is less than 1%, the reconstruction is stopped.2.4 Deep object simulation
One of the largest potential areas for fluence-compensated and quantitative reconstruction methods to improve on FBP is in very deep targets. This experiment demonstrates the ability of the qPAT algorithm to resolve identical targets at increasing distances from the light source. A larger geometry we required, using a 60mm diameter mesh with 7,987 nodes, and 0.6mm node spacing, shown in Fig. 2(a). One distributed light source is arranged at the bottom of the circle and 120 acoustic detectors are equally spaced around the object. In four different cases, objects of 5mm diameter were placed at depths of 10mm, 20mm, 30mm, and 40mm. Data was reconstructed using a time resolution of 200ns. True values of inclusion and background were μa = 0.1mm−1 and 0.01mm−1, respectively, and μs’ = 1mm−1. Data was reconstructed with the initial guess set uniformly to the background value. For comparison, the full resolution data sets were also reconstructed using an FBP method.
Fig. 2 A large geometry of 60mm diameter (a) is used to reconstruct a single 5mm diameter target of 10x contrast. One optical source is located at the bottom of the circle and US transducers surround the entire medium. Target depth is varied from 10mm to 40mm and reconstructed by FBP (b) and qPAT (c). Median recovered absorption coefficient from qPAT (d) is within 22% of the true value as compared with FBP quantification (e).
2.5 Source configuration
Next, we investigated the improvement in quantitative accuracy when using multiple source distributions to reconstruct the image. The three-target simulation was reconstructed with 2 distributed sources positioned 90 degrees apart, and 4 sources distributed 90 degrees apart. The images were reconstructed with sources illuminated separately and simultaneously. During separate illumination reconstructions, multi-source Jacobians were constructed using multiple single-source illuminations. Data and the corresponding Jacobian are inverted at once, processing all source-detector information simultaneously. Quantitative accuracy was assessed versus the true distribution in all regions.
2.6 Spectral reconstruction
The qPAT algorithm can also be extended to spectral reconstructions for imaging of both intrinsic contrast and exogenous dyes. We assess the quantitative accuracy in spectral image recovery in a 30mm diameter circle. The three-target geometry is reconstructed at a 10nm increment from 650nm to 950nm as that is the wavelength range typical of many photo acoustic lasers. The algorithm used four light sources, 120 acoustic detectors, and a time resolution of 200ns. After all reconstructions, the results are fit to known absorption spectra of oxy-hemoglobin (HbO), deoxy-hemoglobin (Hb), water, and AlexaFluor750 using Matlab’s lsqnonneg algorithm. True background concentrations of HbO = 1µM, Hb = 1µM, Water = 40%, and AF750 = 0.1% were used as the initial guess across the entire domain. The targets each had concentrations of HbO = 5, 8, 15µM, Hb = 15, 12, 5µM, water = 60, 70, 80%, and AF750 = 6, 8, 10%. A separate fit is done for each node of the mesh and median values of each chromophore in each region are compared with the true values. After fitting, concentrations were constrained to be greater than 0.001. Each wavelength is also reconstructed with FBP, and the results are passed through the same fitting algorithm.
2.7 Algorithm characterization
To characterize the error propagation from the incorrect initial assignment of both μa and μs’ through the reconstruction, we used the three-target simulation with a time resolution of 200ns and varied the initial μa 60% to 140% from the true value of 0.01mm−1. In a second set of reconstructions, the initial μs’ was varied 60% to 140% of the true value of 1mm−1. All reconstructions were allowed to run for 40 iterations or until the projection error change was less than 1% between iterations. After completion, the quantitative accuracy of the median of the background and each target was assessed versus the true value.
The next set of simulations was to try to understand the relationship between time resolution, computational cost, and quantitative accuracy. The multiple target simulation described in the previous section was reconstructed with time resolutions of 40ns, 60ns, 80ns, 100ns, 200ns, 400ns, 600ns, 800ns, and 1000ns. Assuming a speed of sound of 1500m/s, these time steps lead to a theoretical spatial resolution of 0.06mm, 0.12mm, 0.15mm, 0.3mm, 0.6mm, 0.9mm, 1.2mm, and 1.5mm. For consistency, each reconstruction was allowed to run for up to 40 iterations or until the projection error change was less than 1% between iterations. The performance metrics used were quantitative accuracy and computation time for a single iteration.
3. Results
Experiments were completed to demonstrate quantitative accuracy and improvement over filtered back projection in a 3cm circle. Additional studies show maximal depth recovery, source configuration, and spectral image reconstruction. Finally, three experiments were completed to characterize the effect of initial guess and time resolution on the completed solution.
3.1 Three target reconstruction
The reconstruction shown in Fig. 1 demonstrates the qPAT algorithm’s ability to quantitatively reconstruct absorbing targets at multiple depths using pressure data generated from optical sources. Three targets shadow each other and the fluence map in Fig. 1(b) shows a log order decrease in fluence from one target to the next, despite being only 1cm apart. The initial guess was a uniform absorption of 0.01mm−1. A time step of 200ns corresponded to an iteration time of 2.1 minutes. The problem was allowed to run for up to 40 iterations, decreasing the projection error from 5.4 to 0.036. The qPAT reconstruction was able to recover the background and targets 1, 2 and 3 with absolute error of 7.00%, 0.42%, 1.27%, and 0.96% respectively. The target to background ratio of the target regions was 9.4x, 9.5x, and 9.3x. FBP was unable to provide accurate absolute estimates. The recovered contrast to background ratio of the 3 targets was 5.3x, 4.1x, and 5.1x. The correct value was 10x for all targets. This image is normalized and displayed along with their error quantitation in Fig. 1(f).
3.2 Deep object simulation
Using a 60mm diameter domain, data was reconstructed using both FBP and the qPAT algorithm in four different cases. A single 5mm, 10x contrast target was placed 10mm, 20mm, 30mm, and 40mm away from the single light source. Figure 2(a) shows the geometry, and images from both methods. The qPAT method was able to recover the target concentration with errors of 1%, 9%, 21%, and 22% with increasing depth. The target to background contrast decreased from 10.1x to 7.8x with increasing depth. In the FBP reconstructions, the target-to-background contrast decreased from 2.3x to 1x with increasing depth. Qualitatively, it is difficult to see the 30mm and 40mm targets within the images. All quantifications, and images are shown in Fig. 2.
3.3 Source configuration
The three-target environment was reconstructed using 2 and 4 source distributions both simultaneously and sequentially, as shown in Fig. 3. The reconstructions had a uniform initial guess of μa = 0.01mm−1 and μs’ = 1mm−1, and a time resolution of 200ns. Adding additional sources improved the reconstruction accuracies in the deeper targets as can be seen in Fig. 3. The two-view image had absolute errors of 8.1%, 11.7%, 27.2%, and 22.3% for the background and targets 1, 2, and 3 respectively. The same regions in the two-source image had absolute errors of 3.5%, 4.2%, 12.9%, and 9.8%. The four-view and four source images absolute errors of 5.7%, 12.3%, 22.1%, 14.1%, and 4.4%, 2.5%, 2.4%, 1.6% respectively. Relative contrast of tumor to background ratios for all inclusions was between 7.9x and 9.9x for all source-detector configurations. FBP was again unable to provide accurate absolute estimates and images were normalized for viewing. The recovered contrast to background ratio in the two-view and two-source images was 3.9x, 2.6x, and 2.4x in targets 1, 2, and 3 respectively. In the four-view and four-source images, contrast improved to 5.3x, 4.1x, and 5.1x in the same targets. All quantifications, images, and errors are shown in Fig. 3.
Fig. 3 Data from multiple source configurations was combined into reconstructions of the same 3 targets. Fluence maps to illustrate each illumination pattern (a) and reconstructed image qPAT (b) and FBP (c) are shown. Median recovered mua from qPAT within each region is quantified in (d) and absolute error is shown in (e). Normalized FBP concentrations should be 0 and 1 for background and targets (f).
3.4 Spectral reconstruction
The qPAT algorithm also performed well in spectral image recovery. The imaging domain contained three targets of varying chromophore concentration. 31 single-wavelength reconstructions were combined to quantify oxy and deoxy hemoglobin, total hemoglobin, StO2, water, and exogenous AlexaFluor750 dye. Images, quantification, and errors are shown in Fig. 4. Using the qPAT algorithm, targets were reconstructed with errors spanning 3 to 17% from the true value, with the median error being 8.3%. Chromophore quantification error did not correlate with target depth (Fig. 4(e)-4(g)).
Fig. 4 Spectral reconstruction using 31 wavelengths. Images of true distribution, qPAT image, and normalized FBP image are shown for oxy-hemoglobin (a), deoxy-hemoglobin (b), water (c), and Alexafluor750 (d). True concentrations (e), median qPAT quantifications (f), and error in qPAT reconstructions (g) are within 17% of the truth for all targets. Normalized true concentration (h), FBP medians (i), and FBP error (j) are affected by target depth.
The FBP algorithm reconstructed the same single-wavelength data sets and was reduced to chromophore concentrations. This method was able to generate qualitative images but failed to show any quantitative value. FBP reconstructed targets with quantification error ranging from 1.2% to 117%. The median error was 36%. Chromophore concentration error was strongly correlated with signal-to-background ratio (Fig. 4(h)-4(j)).
3.5 Algorithm characterization
Next, we investigated the effect of incorrect initial guess of both absorption and scatter on the algorithm. The three-target environment was used with a time resolution of 200ns. Data was simulated using a background concentration of μa = 0.01mm−1, and μs’ = 1mm−1. Initial guess was uniform and varied from 60% to 140% of the true value for both μa and μs’. The results for absorption are shown in Fig. 5, where an incorrect initial μa did not strongly affect quantification accuracy. Error concentrations in all regions range from 0.4% to 14%. Relative tumor to background ratios were between 9 and 10x. Images, all quantifications and errors are shown in Fig. 5.
Fig. 5 The effect of the initial guess of μa on image quality is shown (a) and quantified in (b). Median values from within each region based on the true distribution are calculated. Absolute error for background and all anomalies (c) show that the method is robust to incorrect initial guess of μa.
An incorrect initial μs’ also affects the quantitative accuracy of the qPAT algorithm quantification. Initial μa was set to the background value. Error in the background region ranged from 2.2% to 22% but was more pronounced in the target regions, ranging from 0.3% to 108%. The deepest target was the most affected, with errors being generally twice that of the shallower targets. The relative tumor to background ratios were between 5.7x to 16x as scatter was increased from 0.6 to 1.4. All of these results are summarized in Fig. 6.
Fig. 6 The effect of the initial guess of μs’ on image quality is shown (a) and quantified in (b). Median values from within each region based on the true distribution are calculated. Absolute error for background and all anomalies (c) show that incorrect assignment of μs’ can have large effects on overall accuracy.
The three-target domain was reconstructed with time resolutions of 20ns, 40ns, 60ns, 80ns, 100ns, 200ns, 400ns, and 600ns, 800ns, and 1000ns. The images for all of these time intervals are shown in Fig. 7. Each resolution was allowed to run for 40 iterations. Reconstructions using time steps of 200ns and below were consistent, with quantification errors ranging from 0.11% to 8% for the three targets. Time steps above 200ns were not able to accurately reconstruct the 2mm diameter targets, with the 1000ns time point being unable to reconstruct the targets at all. The iteration time of these reconstructions increased from 45 seconds to 9 minutes.
Fig. 7 The effect of time resolution on image quality is shown (a) for time resolutions from 40ns to 1000ns. Reconstruction error increases exponentially when time resolution is not sufficient (b). Median values from within each region based on the true distribution are calculated (c) and absolute error for background and all anomalies (d) show that target quantification accuracy is stable with sufficient time resolution.
4. Discussion
Fluence compensation is an important aspect of photoacoustic imaging that will be critical as the technique moves towards clinical use in large organs [7,31]. Since light can be attenuated an order of magnitude for every centimeter that it travels though tissue, quantitation of optical properties and spectrally derived chromophore concentrations are subject to large errors beyond 10 mm when using FBP and similar reconstruction techniques. This work seeks to address errors caused by structural and spectral distortion by iteratively solving for the tissue absorption coefficient and incorporating fluence fields to markedly improve the absolute quantitation accuracy. The proposed algorithm uses a diffusion model to simulate light propagation through the tissue volume. Other studies have used the Radiation Transport Equation (RTE) and Monte-Carlo models for light distribution but given the large tissue volume and the iterative solution, the diffusion equation is both accurate and appropriate [32]. In Nirfast, point sources are interpolated to the nearest node, meaning that the simulated surface fluence is not technically homogeneous. However, since the domain is large, the light will be diffuse within 1mm of the surface and is acceptable for our set up. Light is converted to sound using a homogeneous Gruneisen coefficient of 0.1, lower than measured values, to err on the side of a less efficient light to sound conversion [33]. Then, sound waves are propagated to the detectors using time steps in the 100s of nanoseconds. Like all time-domain iterative Newton methods, this algorithm’s accuracy and computational complexity are subject to choices of initial guess and time resolution.
4.1 Algorithm performance
As demonstrated in Fig. 1, the algorithm is able to reconstruct an image with high quantitative accuracy in a circle with a 3cm diameter, one of the largest domains used in quantitative photoacoustic reconstruction [22]. When compared with FBP, absolute quantification is within 1.5% of the true value, whereas FBP is not quantitative at all. Target to background ratios had a median error of 6% for qPAT vs. 49% for FBP. The qPAT algorithm produced a quantitative image with nearly the correct contrast enhancement, regardless of target depth (Fig. 1).
In a similar case, the qPAT algorithm was also able to reconstruct realistic targets with only 10x contrast from up to 4cm depth with high quantitative accuracy (Fig. 2). These simulations were done using a 6cm diameter circle with only one light source, the largest qPAT domain in the literature to our knowledge. The geometry allowed for targets to be placed 10, 20, 30, and 40mm from the light source. We were able to clearly distinguish each target, with quantitative accuracy better than 23%, even in the deepest target. The large domain greatly affected the FBP method and was unable to produce a nice image for any of the data sets. The 30 and 40mm targets were especially dim or not recovered at all, with a target to background ratio of nearly 1:1 in both, versus approximately 7.5:1 for qPAT. Other works have demonstrated photoacoustic sensitivity at 40mm and 50mm using wires or photoacoustic dye in phantom studies [31,34,35]. However, wires and concentrated dyes have an extremely high contrast to background ratio that is likely not seen in vivo. This study is able to quantify an absorber with only 10% dye concentration (versus 1% in the background) with high accuracy.
Finally, though all of the targets in this work were given contrast of 10x, the algorithm is also capable or reconstructing smaller contrasts of 2x. This algorithm performed well in two additional simulations (data not shown). The first is a 3cm diameter circle with a 2mm diameter target, 2x contrast, 1cm depth. It was reconstructed with less than 1% absolute error in both regions. The second case was a 6cm diameter circle with a 4mm diameter target, 2x contrast, 3cm depth. It was reconstructed with 7.5% absolute error, consistent with the other simulations in this paper.
4.2 Source configuration
One of the common ways to improve photoacoustic imaging has been to use the highest possible laser power to generate the biggest pressures and SNR. However, light fluence is capped at 30mJ/cm2 due to the ANSI exposure limits and this method is still limited by high light attenuation in tissue. Including additional light sources and detectors can greatly improve the accuracy of all reconstruction methods by increasing the number of measurements, and by giving each measurement a more unique light distribution [36–38]. We evaluated the difference between using two and four illuminations in two different configurations: multiple views and multiple sources. Multiple views simulate the object being rotated 90 degrees, and data sets are collected sequentially then combined. In the other case, multiple sources are illuminated simultaneously as if using a branched laser fiber. Adding additional views into the reconstruction improved accuracy in the targets from median 27% to 14% for two and four views, respectively. Additional sources gave a greater improvement, likely due to a higher overall fluence. Error was reduced from a median 9.8% to 2.4% with two and four sources respectively. The additional projections also improved the FBP reconstruction, but accuracy is very dependent on depth and target-to-background ratios are only half of what the qPAT algorithm recovers. Regardless of reconstruction methods, future system designs should strive for uniform light distribution and simultaneous illumination as some commercial products have done [39].
4.3 Spectral reconstructions
While spectral photo acoustics are certainly not new, it is important that qPAT is able to improve on FBP and mitigate spectral distortion through fluence compensation. Spectral endpoints of HbT, StO2, and water are more clinically relevant, and could provide insight into the functional status of a suspicious lesion [40]. This example of a spectral reconstruction demonstrates that it is possible to recover spectral images with similar accuracy to the single wavelength case. These images are generated from four sources at 31 wavelengths in a 3cm diameter circle. After all reconstructions, oxy and deoxy-hemoglobin, water, and AF750 are unmixed and evaluated. Figure 4 shows all chromophores are recovered within 17% of the correct values across all targets. StO2 and AF750 were accurate within 5%. The FBP images that were analyzed were not able to give the same image quality or quantification accuracy. Targets were identifiable in all images except for water. Background enhancement near the sources could be seen in the Hb and water images, where the contrast to background ratios were the lowest. Water was especially difficult, likely because the contrast was subtle, only 1.5x-2x the background, as opposed to HbO, which was 10x the background. Though targets can be seen on the FBP images, the relative concentrations of the targets are neither accurate to the background nor each other. In this example, accuracy is lost with depth, as seen in the HbO and AF750 images, and spectral results are suspect across all images. With the inclusion of fluence compensation, the qPAT method is able to unmix all targets accurately with a median error of 8.3%. Fluence-compensated algorithms are important for both endogenous and exogenous quantification, and will certainly need to be incorporated into future preclinical and clinical systems [17].
4.3 Algorithm characterization
Though fluence compensation allows for absolute quantification and is a drastic improvement over FBP, there are drawbacks that stem from use of an iterative algorithm. The qPAT algorithm is robust to the initial choice of μa, as incorrect initial μa does not greatly affect the final image. When μa was overestimated by 40%, error in absolute recovered μa concentration had a median error of 9% in the three targets. 40% underestimation led to a 6% median error, depending on the depth of the target. Background concentrations were more accurate with lower initial μa, though the error was never more than 13%. As photoacoustic imaging is only sensitive to light absorption, it is expected that the algorithm will perform well in the case of an incorrect guess to absorption. The range that was tested is applicable to many types of tissue such as breast and prostate [41,42] and is encouraging towards work in these areas.
This algorithm does not reconstruct μs’ since PAI has no way to decouple absorption from scatter. Others have used the spectral changes of scattering amplitude and power to reconstruct them in small domains, but the signals are very subtle [23]. We chose to use a continuous wave approach due to the large size of our domains. Therefore, the initial guess of μs’ is the final answer and its choice will influence the recovered solutions. As shown in Fig. 6, a large error in μs’ causes incorrect fluence fields and leads to large errors. An incorrectly chosen μs’ leads to approximately 1:1 error when it is underestimated and 1.5:1 error when it is overestimated. As such, the choice of μs’ is important to the performance of the algorithm. Again, it may be possible to adequately guess the scattering coefficient based on the literature as done in many optical studies. Alternatively, it could be possible to measure scatter from a time-correlated photo diode to ensure an accurate guess.
Lastly, the time resolution must be sufficiently fine since the spatial resolution is correlated to the length of the time step. For targets of 2mm diameter, a 200ns (0.3mm travel at 1500m/s) time step was the largest step that could be taken and still recover accurate values. Ideally, a much smaller time step could be used, but the reconstruction time increases exponentially with time resolution and is not more accurate for these targets. If smaller targets were used, a shorter time resolution would be required. The Jacobian matrix must be calculated at each time step since the spatial sensitivity of the detectors varies in time. This calculation is a large portion of the computational cost and size scales with time resolution. Reconstruction time is minimized with acceptable loss in accuracy, through the selection of appropriate time resolution. Larger domains will require additional time steps to maintain very high resolution and may approach storage limits. This method could be extendable to 3D domains but would likely require compression or a conjugate gradient approach [20,43]. However, as 2D tomographic geometries are used in many PAI systems, the qPAT algorithm is still applicable here.
Though there is considerable improvement in quantitative accuracy from the qPAT algorithm, the reconstruction time is a drawback in comparison to FBP. In our studies, parallelized spectral reconstructions reduced reconstruction time by 8x using a distributed computing cluster. A larger system would make it possible to reduce this time to at least 30x reduction. Further improvements could be made by implementing code onto a GPU or dedicated reconstruction hardware and we will consider this approach in the future.
This work could also be improved in the future through more efficient jacobian calculation, as was implemented in the diffuse optical tomography community, where Wu et al. demonstrate up to tenfold reduction in processing time [44]. Another alternative would be to implement structural information from an alternate image modality for a priori image guidance. Since photoacoustic imaging is easily paired with ultrasound imaging, it would be logical to incorporate prior information for regularization or weighting during the update step [45]. Future studies on this topic will continue to explore these ideas and we expect considerable time reduction through them.
5. Conclusion
This work has demonstrated a numerical approach to quantitative photoacoustic image reconstruction that minimizes depth and spectrally derived artifacts when compared with filtered back projection. A diffusion-based light model and time-domain acoustic wave propagation generate data and reconstruct quantitative images of absorption coefficient. We demonstrate quantitative accuracy of less than 5% error in large 2D geometries with multiple targets and within 23% error in the largest quantitative photoacoustic studies to date. We extend the algorithm to multi-wavelength data, reconstructing six spectral varying chromophores to within 17% of the true values. The qPAT method was able to improve considerably on filtered back projection from the standpoint of image quality, absolute, and relative quantification. We characterize the effects of time step size, initial guess, and source configuration on final accuracy. This work could help to generate accurate quantitative images from both endogenous absorbers and exogenous dyes in both clinical and preclinical work, thereby increasing the information content obtained from photoacoustic imaging studies.
Acknowledgments
Thanks to grant support from NCI CCNE-TR U54 CA119367 (SSG), CCNE-T U54 U54CA151459, NIBIB BRP 5-RO1-EBB000312 (SSG), NCI T32, NCI U01 EDRN CA152737, Canary Foundation, Ben & Catherine Ivy Foundation, Sir Peter Michael Foundation.
References and links
1. S. Zackrisson, S. M. W. Y. van de Ven, and S. S. Gambhir, “Light in and Sound Out: Emerging Translational Strategies for Photoacoustic Imaging,” Cancer Res. 74(4), 979–1004 (2014). [CrossRef] [PubMed]
2. L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science 335(6075), 1458–1462 (2012). [CrossRef] [PubMed]
3. P. Beard, “Biomedical photoacoustic imaging,” Interface Focus 1(4), 602–631 (2011). [CrossRef] [PubMed]
4. C. Li and L. V. Wang, “Photoacoustic tomography and sensing in biomedicine,” Phys. Med. Biol. 54(19), R59–R97 (2009). [CrossRef] [PubMed]
5. A. De la Zerda, C. Zavaleta, S. Keren, S. Vaithilingam, S. Bodapati, Z. Liu, J. Levi, B. R. Smith, T.-J. Ma, O. Oralkan, Z. Cheng, X. Chen, H. Dai, B. T. Khuri-Yakub, and S. S. Gambhir, “Carbon nanotubes as photoacoustic molecular imaging agents in living mice,” Nat. Nanotechnol. 3(9), 557–562 (2008). [CrossRef] [PubMed]
6. Y. Zhou, W. Xing, K. I. Maslov, L. A. Cornelius, and L. V. Wang, “Handheld photoacoustic microscopy to detect melanoma depth in vivo,” Opt. Lett. 39(16), 4731–4734 (2014). [CrossRef] [PubMed]
7. M. Heijblom, D. Piras, W. Xia, J. C. G. van Hespen, J. M. Klaase, F. M. van den Engh, T. G. van Leeuwen, W. Steenbergen, and S. Manohar, “Visualizing breast cancer using the Twente photoacoustic mammoscope: what do we learn from twelve new patient measurements?” Opt. Express 20(11), 11582–11597 (2012). [CrossRef] [PubMed]
8. R. Li, P. Wang, L. Lan, F. P. Lloyd Jr, C. J. Goergen, S. Chen, and J.-X. Cheng, “Assessing breast tumor margin by multispectral photoacoustic tomography,” Biomed. Opt. Express 6(4), 1273–1281 (2015). [CrossRef] [PubMed]
9. C. Kim, K. H. Song, F. Gao, and L. V. Wang, “Sentinel Lymph Nodes and Lymphatic Vessels: Noninvasive Dual-Modality in Vivo Mapping by Using Indocyanine Green in Rats--Volumetric Spectroscopic Photoacoustic Imaging and Planar Fluorescence Imaging,” Radiology 255(2), 442–450 (2010). [CrossRef] [PubMed]
10. K. H. Song, E. W. Stein, J. A. Margenthaler, and L. V. Wang, “Noninvasive photoacoustic identification of sentinel lymph nodes containing methylene blue in vivo in a rat model,” J. Biomed. Opt. 13(5), 054033 (2008). [CrossRef] [PubMed]
11. A. Garcia-Uribe, T. N. Erpelding, A. Krumholz, H. Ke, K. Maslov, C. Appleton, J. A. Margenthaler, and L. V. Wang, “Dual-Modality Photoacoustic and Ultrasound Imaging System for Noninvasive Sentinel Lymph Node Detection in Patients with Breast Cancer,” Sci. Rep. 5, 15748 (2015). [CrossRef] [PubMed]
12. J.-M. Yang, C. Li, R. Chen, B. Rao, J. Yao, C.-H. Yeh, A. Danielli, K. Maslov, Q. Zhou, K. K. Shung, and L. V. Wang, “Optical-resolution photoacoustic endomicroscopy in vivo,” Biomed. Opt. Express 6(3), 918–932 (2015). [CrossRef] [PubMed]
13. J.-M. Yang, C. Favazza, R. Chen, J. Yao, X. Cai, K. Maslov, Q. Zhou, K. K. Shung, and L. V. Wang, “Simultaneous functional photoacoustic and ultrasonic endoscopy of internal organs in vivo,” Nat. Med. 18(8), 1297–1302 (2012). [CrossRef] [PubMed]
14. 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] [PubMed]
15. B. T. Cox, S. R. Arridge, K. P. Köstli, and P. C. Beard, “Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method,” Appl. Opt. 45(8), 1866–1875 (2006). [CrossRef] [PubMed]
16. B. Banerjee, S. Bagchi, R. M. Vasu, and D. Roy, “Quantitative photoacoustic tomography from boundary pressure measurements: noniterative recovery of optical absorption coefficient from the reconstructed absorbed energy map,” J. Opt. Soc. Am. A 25(9), 2347–2356 (2008). [CrossRef] [PubMed]
17. D. Razansky and V. Ntziachristos, “Hybrid photoacoustic fluorescence molecular tomography using finite-element-based inversion,” Med. Phys. 34(11), 4293–4301 (2007). [CrossRef] [PubMed]
18. D. Razansky, A. Buehler, and V. Ntziachristos, “Volumetric real-time multispectral optoacoustic tomography of biomarkers,” Nat. Protoc. 6(8), 1121–1129 (2011). [CrossRef] [PubMed]
19. A. Q. Bauer, R. E. Nothdurft, T. N. Erpelding, L. V. Wang, and J. P. Culver, “Quantitative photoacoustic imaging: correcting for heterogeneous light fluence distributions using diffuse optical tomography,” J. Biomed. Opt. 16(9), 096016 (2011). [CrossRef] [PubMed]
20. S. Bu, Z. Liu, T. Shiina, K. Kondo, M. Yamakawa, K. Fukutani, Y. Someda, and Y. Asao, “Model-Based Reconstruction Integrated with Fluence Compensation for Photoacoustic Tomography,” IEEE Trans. Biomed. Eng. 59(5), 1354–1363 (2012). [CrossRef] [PubMed]
21. Z. Yuan and H. Jiang, “Three-dimensional finite-element-based photoacoustic tomography: reconstruction algorithm and simulations,” Med. Phys. 34(2), 538–546 (2007). [CrossRef] [PubMed]
22. L. Yao and H. Jiang, “Finite-element-based photoacoustic tomography in time domain,” J. Opt. A, Pure Appl. Opt. 11(8), 085301 (2009). [CrossRef]
23. 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–455 (2009). [CrossRef] [PubMed]
24. 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] [PubMed]
25. S. L. Jacques and B. W. Pogue, “Tutorial on diffuse light transport,” J. Biomed. Opt. 13, 041302 (2008).
26. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng. 25(6), 711–732 (2009). [CrossRef] [PubMed]
27. M. B. Moffett, “Theoretical Acoustics. Philip M. Morse and K. Uno Ingard. McGraw-Hill, New York, 1968. xxii, 938 pp., illus. $23. International Series in Pure and Applied Physics,” Science 170(3954), 156–157 (1970). [CrossRef]
28. G. J. Diebold, T. Sun, and M. I. Khan, “Photoacoustic monopole radiation in one, two, and three dimensions,” Phys. Rev. Lett. 67(24), 3384–3387 (1991). [CrossRef] [PubMed]
29. B. T. Cox, S. Kara, S. R. Arridge, and P. C. Beard, “k-space propagation models for acoustically heterogeneous media: application to biomedical photoacoustics,” J. Acoust. Soc. Am. 121(6), 3453–3464 (2007). [CrossRef] [PubMed]
30. K. M. Tichauer, R. W. Holt, F. El-Ghussein, Q. Zhu, H. Dehghani, F. Leblond, and B. W. Pogue, “Imaging workflow and calibration for CT-guided time-domain fluorescence tomography,” Biomed. Opt. Express 2(11), 3021–3036 (2011). [CrossRef] [PubMed]
31. S. R. Kothapalli, T.-J. Ma, S. Vaithilingam, O. Oralkan, B. T. Khuri-Yakub, and S. S. Gambhir, “Deep tissue photoacoustic imaging using a miniaturized 2-D capacitive micromachined ultrasonic transducer array,” IEEE Trans. Biomed. Eng. 59(5), 1199–1204 (2012). [CrossRef] [PubMed]
32. B. W. Pogue, S. Geimer, T. O. McBride, S. Jiang, U. L. Osterberg, and K. D. Paulsen, “Three-dimensional simulation of near-infrared diffusion in tissue: boundary condition and geometry analysis for finite-element image reconstruction,” Appl. Opt. 40(4), 588–600 (2001). [CrossRef] [PubMed]
33. D.-K. Yao, C. Zhang, K. Maslov, and L. V. Wang, “Photoacoustic measurement of the Grüneisen parameter of tissue,” J. Biomed. Opt. 19(1), 017007 (2014). [CrossRef] [PubMed]
34. C. Kim, T. N. Erpelding, L. Jankovic, M. D. Pashley, and L. V. Wang, “Deeply penetrating in vivo photoacoustic imaging using a clinical ultrasound array system,” Biomed. Opt. Express 1(1), 278–284 (2010). [CrossRef] [PubMed]
35. G. Ku and L. V. Wang, “Deeply penetrating photoacoustic tomography in biological tissues enhanced with an optical contrast agent,” Opt. Lett. 30(5), 507–509 (2005). [CrossRef] [PubMed]
36. A. Pulkkinen, B. T. Cox, S. R. Arridge, J. P. Kaipio, and T. Tarvainen, “Quantitative photoacoustic tomography using illuminations from a single direction,” J. Biomed. Opt. 20(3), 036015 (2015). [CrossRef] [PubMed]
37. H. Gao, J. Feng, and L. Song, “Limited-view multi-source quantitative photoacoustic tomography,” Inverse Probl. 31(6), 065004 (2015). [CrossRef]
38. N. Song, C. Deumié, and A. Da Silva, “Considering sources and detectors distributions for quantitative photoacoustic tomography,” Biomed. Opt. Express 5(11), 3960–3974 (2014). [CrossRef] [PubMed]
39. D. Razansky, “Multispectral Optoacoustic Tomography -Volumetric Color Hearing in Real Time,” IEEE J. Sel. Top. Quantum Electron. 18(3), 1234–1243 (2012). [CrossRef]
40. M. A. Mastanduno, J. Xu, F. El-Ghussein, S. Jiang, H. Yin, Y. Zhao, K. Wang, F. Ren, J. Gui, B. W. Pogue, and K. D. Paulsen, “MR-guided near infrared spectral tomography increases diagnostic performance of breast MRI,” Clin. Cancer Res. 21, 3906–3912 (2015).
41. M. Solonenko, R. Cheung, T. M. Busch, A. Kachur, G. M. Griffin, T. Vulcan, T. C. Zhu, H. W. Wang, S. M. Hahn, and A. G. Yodh, “In vivo reflectance measurement of optical properties, blood oxygenation and motexafin lutetium uptake in canine large bowels, kidneys and prostates,” Phys. Med. Biol. 47(6), 857–873 (2002). [PubMed]
42. A. Cerussi, N. Shah, D. Hsiang, A. Durkin, J. Butler, and B. J. Tromberg, “In vivo absorption, scattering, and physiologic properties of 58 malignant breast tumors determined by broadband diffuse optical spectroscopy,” J. Biomed. Opt. 11(4), 044005 (2006). [CrossRef] [PubMed]
43. A. Rosenthal, D. Razansky, and V. Ntziachristos, “Quantitative Optoacoustic Signal Extraction Using Sparse Signal Representation,” IEEE Trans. Med. Imaging 28(12), 1997–2006 (2009). [CrossRef] [PubMed]
44. X. Wu, A. T. Eggebrecht, S. L. Ferradal, J. P. Culver, and H. Dehghani, “Fast and efficient image reconstruction for high density diffuse optical imaging of the human brain,” Biomed. Opt. Express 6(11), 4567–4584 (2015). [CrossRef] [PubMed]
45. J. D. Gruber, A. Paliwal, V. Krishnaswamy, H. Ghadyani, M. Jermyn, J. A. O’Hara, S. C. Davis, J. S. Kerley-Hamilton, N. W. Shworak, E. V. Maytin, T. Hasan, and B. W. Pogue, “System development for high frequency ultrasound-guided fluorescence quantification of skin layers,” J. Biomed. Opt. 15(2), 026028 (2010). [CrossRef] [PubMed]