Adaptive photoacoustic computed tomography.

For many optical imaging modalities, image qualities are inevitably degraded by wavefront distortions caused by varying light speed. In optical microscopy and astronomy, adaptive optics (AO) has long been applied to compensate for such unwanted aberrations. Photoacoustic computed tomography (PACT), despite relying on the ultrasonic wave for image formation, suffers from the acoustic version of the same problem. However, this problem has traditionally been regarded as an inverse problem of jointly reconstructing both the initial pressure and the sound speed distributions. In this work, we proposed a method similar to indirect wavefront sensing in AO. We argued that wavefront distortions can be extracted and corrected by a frequency domain analysis of local images. In addition to an adaptively reconstructed aberration-free image, the speed of sound map can be subsequently estimated. We demonstrated the method by in silico, phantom, and in vivo experiments.

Introduction

The photoacoustic (PA; or optoacoustic, OA) effect refers to the process of converting light absorption into ultrasound emission. It offers a unique way of measuring optical absorption by ultrasonic detection [1]. Photoacoustic computed tomography (PACT), which employs spatially distributed detection of PA signals and digital image reconstruction, allows tomograms of light absorption to be formed in the optically diffusive regime with relatively high resolution [2].

During PA image reconstruction, the commonly assumed constant speed of sound (SOS) is only an approximation, which generates image artifacts in the presence of acoustic heterogeneity. If the SOS distribution is measured to assist image reconstruction, the image quality can be greatly improved [[3], [4], [5], [6]]. However, the measurement of the SOS requires additional hardware with increased system complexity and cost. Alternatively, it is more desirable to perform joint reconstruction (JR) of the SOS and the PA initial pressure (IP), based upon PA signals alone [[7], [8], [9], [10], [11], [12], [13], [14], [15], [16], [17]]. To do so, the frequency domain [[8], [9], [10]] or time domain [[11], [12], [13]] model-based reconstructions incorporating the acoustic wave equation in heterogeneous media can be adopted. The solution is usually found by iterative optimization to minimize the difference between the measured signals and the predicted ones. Back-projection based JR techniques count on the correlation between images reconstructed using partial data (e.g. two half-ring reconstructions [14], first and second half-time reconstructions [15]) to search for the correct SOS distribution. However, mathematically, JR of the SOS and the IP is fundamentally unstable and the uniqueness of the solution is only guaranteed in some special cases [[18], [19], [20]]. Numerical studies also found that the optimization process in the JR problem is not convex [12], which means that the solution can easily fall into local extrema.

Non-uniform SOS causes wavefront distortions, which ultimately give rise to resolution reduction and image artifacts. Here, comparing the wavefront distortions caused by spatially varying wave speed in PACT with that in light microscopy is enlightening. As shown in Fig. 1, central to the image formation of light microscopy and PACT is a mapping (i.e., ideally representing an identity matrix) from the object domain to the image domain. In optical microscopy, such a mapping is simply done by physical lenses, while in PACT, the mapping is accomplished digitally by reversing the recorded PA waveforms. Fig. 1(b) depicts the imaging process of a standard PACT system based on a ring array [21]. In this case, ultrasound detection is focused in the elevational direction, allowing thin slices of the object to be imaged in 2D (i.e., sectioning). To refocus an image, the back-projection image reconstruction algorithm [22] is commonly used, which can be interpreted as triangulation localization [23]. Specifically, the PA signals are back-projected to arcs that are centered at the positions of the transducer elements. In particular, to image a point source, the radii of the arcs are determined by timing the PA waveforms and assuming a constant speed of sound. Ideally, the arcs are expected to intercept and constructively interfere at the source point to form a bright spot.

Fig. 1

Comparison between optical microscopy and PACT systems. Both imaging processes are affected by wavefront distortion. (a) A wide-field fluorescence microscope. The red dashed lines are the ideal distortion-free wavefronts and the green lines are the distorted wavefronts and light rays. (b) A PACT imaging system. The red dashed circle is the ideal distortion-free photoacoustic wavefront. The green curves are distorted wavefronts. In the digital domain, the green arcs are the back-projection curves in the PACT reconstruction algorithm. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

In fluorescence microscopy, the optical wavefront from a fluorescent point source is often distorted by the spatially varying refractive index of the biological sample, as shown in Fig. 1(a) [24]. The light rays, as depicted by the green lines, fail to converge to the same point on the image plane, resulting in degraded resolutions. In PACT, how the unknown SOS distribution produces wavefront aberration can be understood similarly. As shown in Fig. 1(b), if the acoustic refractive index, i.e. the SOS, is non-uniform, the PA wavefront from a point source will be distorted and eventually back-projected into an irregular pattern [25]. As a result, the reconstructed image features will be distorted and/or split.

The wavefront distortions in optical microscopy can be measured, either directly [[26], [27], [28], [29], [30], [31], [32], [33]] by a Shack-Hartmann wavefront sensor or indirectly [[34], [35], [36], [37], [38], [39], [40], [41]] and subsequently compensated for by a spatial light modulator or a deformable mirror. Known as adaptive optics (AO), this technology was first implemented in astronomic telescopes and then successfully translated into microscopy [42]. In terms of wavefront distortion caused by unknown SOS distribution, the similarity between optical microscopy and PACT suggests that the concept of AO can be extended beyond purely optical imaging. In this paper, we implemented a technique that can solve the JR problem analogous to indirect wavefront measurement in conventional AO: The image is reconstructed patch by patch; within each patch, the wavefront distortion is almost identical (i.e., an “isoplanatic patch”) and can be extracted from the local point spread function (PSF). The advantage of indirect AO is that it does not need a point target to be created as a real “guide star”, while seeks to optimize the compensation wavefront from a number of images at the cost of increased total exposure times [43]. In our method, the PACT image series can be generated digitally rather than physically, and our wavefront extraction is similar to the “phase diversity” method [44]. The local PSF, which has long been regarded as an unknown, can be computationally found from a stack of local images reconstructed with different delays. Thereby, the full image can be better focused via piecewise deconvolution. As a bonus, once the wavefronts of all the patches are determined, they can be used collectively to compute the global SOS map. Unlike all JR approaches which tend to be corrupted by incorrect guess of the SOS distribution, our method bypasses the cumbersome global searching of the SOS map, thus improves the stability and reliability of the solution. To acknowledge the key role played by the concept of AO, we named this technique Adaptive PACT (APACT) and demonstrated in silico, phantom, and in vivo experiments.

Method

Wavefront function and point spread function

In the focused ring array PACT system, 512 transducer elements are evenly arranged on a circle with a diameter of 10 cm. A detailed description of the system and its operating process can be found in [14].

In a homogenous medium with a constant SOS of , the universal back-projection formula [22] directly links , the IP at time zero when the optical energy is absorbed, to , the pressure on a surface that encloses the sources:where is the solid angle subtended by a transducer element at from the source point at . In real practice, the second term on the right is often dominant [22]. The above formula assumes the detectors have infinite bandwidth, uniform angular response, and be distributed over solid angle. In reality, the transducers are band limited and their angular responses are not flat. Also, the array is focused in the elevational direction, converting the reconstruction from 3D to 2D. Due to these limitations, we implemented the following delay and sum (DAS) algorithm instead of Eq. (1) in practice:where and are the location (in 2D) and the output of the n’th transducer, respectively. is the total number of transducers. is the reconstructed image, which is an approximation to . In real practice, Eq. (2) is a reasonable approximation to the back-projection formula, because the response of our transducer can be roughly treated as a time derivative. In addition, although the method described in this paper is based on DAS, it can be applied to the back-projection formula (Eq. (1)) in which additional linear filtering is performed on .

Now, we consider a point source located at and analysis its image. The back-projection arcs, mathematically described by , will intercept at only if the real time of flight (TOF) of the PA signal coincides with that of a perfect spherical wavefront propagating at . We define any deviation from the above ideal case as the wavefront function:Here denotes the angle of a vector pointing from to , and is the ultrasound propagation time from to .

By adding up all the back projection arcs, the image of a point source with unit amplitude (the point spread function, PSF) is obtained, mathematically expressed as . In the DAS algorithm, the back projection arcs are broadened into thin bands whose cross-sectional profiles are determined by the electrical impulse response (EIR) of the transducers . It can be shown that the PSF is mathematically determined by the wavefront function (see Appendix A for further details):

Here is a constant and denotes the position in a shifted coordinate system whose origin coincides with . The unit vector , whose polar angle is , represents the backpropagation direction. As in AO microscopy, the wavefront function changes with . However, it is reasonable to treat to be unchanged inside the isoplanatic patch. Under this assumption, the PSF becomes space-invariant, and thus the local image can be expressed as a convolution:

Wavefront extraction and adaptive reconstruction algorithm

Fig. 2 illustrates the working principle of the algorithm using simulation data, where the red box in Fig. 2(a) depicts one representative isoplanatic patch. The algorithm starts by adding a varying delay during the DAS process, resulting inwhere is the extra delay distance and is the corresponding image of the isoplanatic patch reconstructed by DAS. By doing so, the radii of the back-projection arcs vary by the same amount of and the PSF becomes . It should be noted that similar operations are common in conventional DAS when is adjusted [45] to better focus an image. This can be regarded as the 0th order compensation, and the process is illustrated in Fig. 2(b).

Fig. 2

Diagram of the APACT algorithm. All images are obtained from numerical simulation. (a) Image reconstructed with conventional DAS. The red box is one isoplanatic patch and the blue box is the successive one. (b) Left column: a stack of images reconstructed with increasing ; In each plot, the outline of a point spread function is highlighted in green. Right column: in each plot 16 back-projection arcs are shown to explain how the PSF is formed. The curvatures of the arcs are exaggerated. (c) A 3D rendering of the data cube . (d) The Fourier transform of the images stacked together as in (c) with several radial slices. (e) The Fourier transform of the aberration-free image . (f) The transfer function in which we also provide some slices to show the pattern. (g) The recovered aberration-free image . (h) The extracted wavefront function . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

Now we prove that from a series of images (), we are ready to extract the wavefront function as well as the aberration-free image , whose PSF corresponds to vanished wavefront aberration: . This technique is similar to the phase diversity [44] method in optics. Note that when is varied, the local image will change accordingly. This is because is generated by convolving the initial pressure with a PSF whose shape is determined by and . So is encoded in the datacube (the third dimension is ). In our study, we scanned in 161 discrete steps. The interval was . In Eq. (6) we can see that this interval is small enough for the 40 MHz sampling rate. The maximum length of (both negative and positive value) is chosen to be a quarter of the edge length of the isoplanatic patch. This can guarantee that the defocused local image features are mostly confined within the patch as shown in Fig. 2 (b), which means the back-projection curves do not completely move out of the patch window. To extract the information concealed by the convolution operation, we turn to the spatial frequency domain as shown in Fig. 2(c) and (d). By doing so, convolution is converted to multiplication, and the effect of wavefront distortion is readily modeled by a multiplicative transfer function :where is the 2D Fourier spectrum of an image with a specific delay . is the k-space coordinate in the spatial frequency domain. is the Fourier transform of the aberration-free image . From (4), one can express as a function of . In polar coordinates, it reads (see Appendix B for further details):where is the imaginary unit. The Fourier domain relationship in Eq. (7) is the key to our APACT algorithm. As shown in Fig. 2(d), (e), and (f), we can decompose , into the product of two terms: and . The decomposition can be accomplished using the least-square optimization shown below:

A factor is introduced because in the Fourier domain the noise level is inversely proportional to (see Appendix B for further details). This optimization not only extracts the wavefront function but also estimates the Fourier domain aberration-free image . An inverse Fourier transform is performed to recover the spatial domain image as shown in Fig. 2(g).

The optimization of Eq. (9) is done by exhaustive search for and a linear least square projection for . In AO, the wavefront distortion function is often decomposed into Zernike modes where higher-order modes can be discarded [35]. In our implementation, the searching of is confined to its 0th and 2nd Fourier series coefficients, which can greatly reduce the computation time. In Appendix C we prove that the 1st order coefficient cannot be extracted, the situation bears a resemblance to the tip and tilt Zernike modes [46] that can only result in a position shift of the image but has no influence on image quality.

Image stitching and SOS reconstruction

After all the isoplanatic patches are deconvolved (either sequentially or in parallel) following the above procedure, we stitch them together to get the full image. In addition, we can also use the extracted wavefronts to calculate an SOS distribution.

In our simulation, phantom, and in vivo experiments, the size of the isoplanatic patch was chosen to be . To apply Fourier transform, a Gaussian window was multiplied, confining the effective area to a disk with a full width at half maximum of 1.5 mm. The recovered images inherit the Gaussian modulation, as can be seen in Fig. 2(g). These local images are stitched together, with neighboring patches partially overlapped and superimposed, to generate the full image. By design, a 75% overlap between neighboring patches was employed, as shown in Fig. 2(a) by the red and blue dashed boxes. The overall shape, resulting from adding the Gaussian profiles, will become flat. Because the Gaussian window is close to zero on the patch edges, the interpatch transition will be smooth.

To estimate the SOS map, acoustic refraction was neglected, i.e., ultrasound rays do not bend during propagation and the resulting TOF is computed by a straight line integral of the reciprocal of the SOS [25]. We built linear equations on these reciprocals and solved the equations with regularization. A detailed algorithm is provided in Appendix D. We only computed the SOS in the sample/animal body while the SOS of the surrounding water was estimated by measuring the temperature [47].

Result

Numerical simulation

We used k-Wave toolbox [48] to do the simulation in 2D. To approximate the real case with transducer response, we take the derivative of the simulated signals: . The numerical phantom had an SOS of 1600 m/s with a circular inclusion whose SOS was 1650 m/s, as shown in Fig. 3(e). The SOS of the surrounding water was set to 1500 m/s. We used a vascular mimicking pattern for the PA contrast (Fig. 3(b)).

Fig. 3

Numerical simulation results. (a) The gold standard of the IP distribution. Scale bar: 1 mm. (b) Image reconstructed using conventional DAS for an acoustically homogeneous phantom with SOS = 1500 m/s. (c) Image reconstructed using conventional DAS for the acoustically inhomogeneous phantom shown in (e) with SOS = 1520 m/s. (d) Image reconstructed using APACT. (e) The gold standard of the SOS map. (f) Recovered SOS map. (e) and (f) share the same color bar. (g) The extracted wavefronts are drawn on top of the original image. Red ellipses: estimated wavefronts with a relative fitting error of less than 50%; Purple ellipses: estimated wavefronts with a relative fitting error above 50%. Green ellipses: real wavefronts. Blue circle (bottom left): aberration-free wavefront propagating at 1520 m/s. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

The results of image reconstruction are shown in Fig. 3. Compared to the severely distorted image shown in Fig. 3(c), which was reconstructed by DAS assuming a constant speed of sound (the SOS was already tuned for best visual effect), APACT corrected the aberrations to yield better image quality (Fig. 3(d)). The image quality is comparable to Fig. 3(b), the best aberration-free image we can get when the SOS is homogenous.

The extracted wavefronts and the SOS reconstruction result are also shown. The wavefront distortions are displayed in polar plots as in Fig. 2(h). In Fig. 3(g), each little ellipse represents a snapshot of the acoustic wavefront outside the phantom emanating from a virtual point source located at the center of the ellipse, while the distortion is exaggerated. The ellipses tend to align their major axes radially because traversing the phantom advances the wavefront due to a higher average SOS.

In places with no or weak PA features, it is not possible to use the algorithm to extract the wavefront due to a lack of “guide stars”. In Fig. 3(g) we also plotted the correct wavefronts (only 0th and 2nd order distortions were extracted and drawn), which were calculated directly from Snell’s law. As shown in Fig. 3(g), in the featureless regions, wavefront estimation is prone to errors. Although noise was absent in the simulation, the energy leakage from nearby patches dominated and generated errors. Such unfavorable situations can be identified during the fitting process of Eq. (9), where the mismatch between the two terms is noticeably high irrespective of the choice of . A “relative error” can be defined as the value of the mismatch normalized by the total energy of . To calculate the SOS map, in the simulation we discarded all the wavefronts whose relative error exceeded 50%, as shown in purple in Fig. 3(g). The result of the SOS distribution reconstruction is shown in Fig. 3(f). Compared to the real SOS distribution shown in Fig. 3(e), estimation errors do exist; however, a general similarity is observed.

Phantom experiment

To verify the performance of APACT, we did a phantom experiment. We used a colored leaf vein as the PA target. A photo of the phantom is shown in Fig. 4(c).

Fig. 4

Phantom imaging results. (a) Image reconstructed using conventional DAS. (b) Image reconstructed using APACT. (c) Photo of the phantom. (d) Top-down: zoomed-in views of the regions highlighted in red dashed boxes in (a), (c), and (b). (e) Recovered SOS map. (f) Red ovals: extracted wavefronts with relative fitting errors below 70%. Blue circle (bottom left): aberration-free wavefront propagating at 1505 m/s. Scale bar: 1 mm. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

We embedded the leaf in a 2.5 cm diameter agarose cylinder (solution made by mixing 4 g agarose with 100 g water). The dyed leaf vein was embedded into the phantom perpendicular to the height direction. A cylindrical hole running parallel to the phantom was created, which was filled with water during the experiment. Because agarose has a higher SOS than water, this phantom represents a case of very uneven SOS distribution. The water temperature was measured to be 26℃ during the experiment.

An image reconstructed using conventional DAS (SOS = 1505 m/s) is shown in Fig. 4(a). The laser illumination was at 700 nm. Because the phantom was not optically turbid enough, illumination inside the phantom was not uniform, resulting in lower SNR on the upper left corner of the image. In Fig. 4(b) we presented the result of the image reconstructed using APACT. Compared to Fig. 4(a), the image quality is apparently improved. From the zoomed-in views shown in Fig. 4(d), we can compare the reconstruction results with the real structure of the leaf vein to examine the fidelity of different reconstruction algorithms.

The extracted wavefronts are depicted in Fig. 4(f). In the experiment, due to the existence of noise, we increased the threshold for relative error to 70%, and the survived wavefronts are shown in Fig. 4(f). The missing wavefront plots on the upper-left corner correspond to the region with low SNR.

The reconstruction of the SOS distribution is shown in Fig. 4(e). The regions with higher SOS correlate reasonably well with the solid agar shown in Fig. 4(c). It is worth mentioning that the ellipticity of the ring array, due to manufacturing imperfection, can also contribute to the wavefront function defined in Eq. (3). So to compute the SOS map, should be adjusted to remove this geometrical offset. We measured the geometrical error using a PA point source, and subtracted its 0th and 2nd order components from before drawing Fig. 4(f) and computing the SOS distribution. We followed the same procedure in the in vivo experiment.

In vivo experiment

In the in vivo experiment, we imaged a nude mouse (8-week-old Crl: NU-Foxn 1nu mouse) in the liver region at 1064 nm. (The same data has been used in another work [14] of our lab.) The water temperature was maintained at 31℃ during the experiment. Compared to Fig. 5(a) which was reconstructed using conventional DAS, the APACT reconstruction result shown in Fig. 5(b) shows better quality. For example, the vertical vessels close to the surface of the body are shown as bright spots in Fig. 5(b) but appear to be highly distorted in Fig. 5(a).

Fig. 5

In vivo imaging results. (a) Image reconstructed using conventional DAS. (b) Image reconstructed using APACT. (c, d) The zoomed-in area of (a) and (b) in the red and blue dashed line boxes, respectively. (e) Red ovals: extracted wavefronts with relative fitting errors below 70%. Blue circle (bottom left): aberration-free wavefront propagating at 1515 m/s. (f) Recovered SOS map. Scale bar: 1 mm. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

Two zoomed-in areas with abundant features of blood vessels are shown in Fig. 5(c) and (d). Without a gold standard, it is difficult to judge which reconstruction outperforms the other. However, the APACT images appear to be superior visually with smoother texture and sharper features, which are hallmarks of better focusing. Similar improvements are prevalent in the entire image.

We also chose to use the wavefronts whose relative fitting error is less than 70% to calculate the SOS map. The wavefront map shown in Fig. 5(e) is physically reasonable because the ovals near the edge tend to align their major axes towards the center, while those close to the center are slightly dilated yet less elongated. This trend can be explained by the higher SOS of the mouse than the surrounding water. The computed SOS distribution is shown in Fig. 5(f), indicating a higher SOS value in some parts of the body, which can be anatomically identified as the liver.

Discussion

In this work, we propose and demonstrate a PACT image reconstruction algorithm to remove artifacts caused by non-uniform SOS. Traditional JR methods optimize the SOS distribution together with the PA image iteratively. Without prior knowledge about the SOS distribution, the degrees of freedom of the overall optimization is large. In such cases, the iteration may fall into a local extremum, while the resulting reconstruction is suboptimal. However, our method establishes an entirely new strategy to avoid such a problem. We found that the problem of solving the unknown SOS distribution can be decomposed into independent local wavefront extraction problems. This means the overall image quality does not suffer much from an incorrectly resolved local wavefront.

The estimation accuracy of relies on the sufficiency of the spatial frequency components in all directions in the selected patch. With incomplete frequency components, wavefront estimations are likely to be corrupted. This can be seen in both the simulation and the experimental results. Fortunately, image correction is less affected because the spectrum of the image has little proportion of energy in those unfavorable directions. On the contrary, to compute the SOS map, the set of equations can be underdetermined to begin with, whereas incorrect wavefront estimations simply exacerbate the situation. The regularization we employed in this study assumed that SOS variations are spatially slow, which may not be true in reality. Despite the above difficulties, the major advantage of APACT over JR is that the IP can be recovered reliably if PA features are abundant, and for each recovered patch, the accuracy of the extracted wavefront can be inferred from the fitting error of Eq. (9).

As a newly developed reconstruction algorithm, APACT’s capacity is yet to be fully developed. Examining Fig. 3(b) and (d), one can see that the recovered image is smaller as if a scaling factor is applied. This is because the 1st order wavefront distortion cannot be corrected by our algorithm. The resolved image features shift inwards due to the different wavefront distortion in opposing directions. This 1st order distortion can be calculated using the reconstructed SOS distribution. Also, making use of global information can make the solution more stable. Further restrictions can be imposed on from a global perspective to avoid apparent inconsistency among different patches. To further improve the performance of the algorithm, the frequency domain technique using Eq. (9) can also be modified – blind deconvolution, and prior knowledge about such as nonnegativity [49] can be applied for better accuracy.

Although the focused ring array system was used throughout this study, the method can be extended to other array geometries with modified formulas. In the current system, the out-of-focus signals may contaminate the in-plane ones, making our model inaccurate and adversely affect the overall performance. So we expect even better performance when both PA signal detection and image reconstruction are extended to 3D.

Conclusion

Herein we present an image reconstruction algorithm in PACT with adaptive wavefront correction that bears a strong resemblance to AO in optical microscopy. Using the new method, PACT image aberrations induced by SOS heterogeneity can be remedied locally. The wavefront distortion is defined as a perturbation added on the ideal spherical wavefront. We derived that this distortion, which corresponds to a signal arrival time difference, will mathematically determine a distorted point spread function. In an isoplanatic patch where wavefront distortion is almost translational invariant, the distorted image can be regarded as the initial pressure convolved with a distorted PSF. To locally decode the wavefront information, a phase diversity technique is used: A series of local images within an isoplanatic patch is reconstructed, in which various delays are used in the conventional delay-and-sum image reconstruction process. These images originate from the same target but are blurred by different PSFs. From the Fourier transforms of these images, we can extract both the wavefront distortion and the aberration-free image. The patches are then stitched together to form the whole image. Then by solving a set of linear equations, we used the extracted wavefronts to calculate the SOS map. Because the PA image and SOS map are no longer simultaneously reconstructed, our method is more stable and is unlikely to fall into local extrema. We demonstrated the effectiveness of our method using in silico, phantom, and in vivo experiments. The method is developed in particular for the full ring array system. In the future, it can be further improved and extended to other PACT detection geometry. We also believe that the concept of APACT can benefit other imaging fields that involve wave propagation in heterogeneous media (e.g., optical microscopy and ultrasound imaging), for simultaneous image deblurring and mapping of the wave speed distribution.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.