Streak artifact suppression in photoacoustic computed tomography using adaptive back projection.

For photoacoustic computed tomography (PACT), an insufficient number of ultrasound detectors can cause serious streak-type artifacts. These artifacts get overlaid on top of image features, and thus locally jeopardize image quality and resolution. Here, a reconstruction algorithm, termed Contamination-Tracing Back-Projection (CTBP), is proposed for the mitigation of streak-type artifacts. During reconstruction, CTBP adaptively adjusts the back-projection weight, whose value is determined by the likelihood of contamination, to minimize the negative influences of strong absorbers. An iterative solution of the eikonal equation is implemented to accurately trace the time of flight of different pixels. Numerical, phantom and in vivo experiments demonstrate that CTBP can dramatically suppress streak artifacts in PACT and improve image quality.


Introduction
Photoacoustic tomography (PAT) can form images with contrast of optical absorption and spatial resolution of ultrasound imaging [1]. Due to multiple scattering of photons, traditional optical imaging has poor resolution in deep tissue. For PAT, the pressure waves, which have three orders of magnitude weaker scattering than light in biological tissues [1], are generated by absorption of short-pulse light and detected by ultrasound detectors. Therefore, PAT can break the optical diffusion limit [2].
In actual applications, the number of elements in the photoacoustic (PA) signal detection array is often insufficient, due to the highly complex PA waveform on one hand, and to system cost on the other. This results in 1) spatial sampling below the Nyquist criterion and 2) limited view (with detection solid angle of less than 4π), and can cause serious streak-type artifact, one of the main artifacts of PAT. These streak patterns are similar to the undersampling artifacts in X-ray CT [3] which are caused by spatial aliasing, and can be very strong and get overlaid on top of real image features, and thus locally jeopardize image quality and resolution [4].
Model based methods can introduce constraints to yield sharp image with weak streak artifact, at the cost of large computational workload and memory demand. For example, Deán-Ben et al. proposed a three-dimensional model-based method based on a discretization of the Poisson-type integral model, and reduced streak-type artifacts due to limited number of measuring points [4]. Sparsity constraints, including the L 1 norm of the image in sparse representation and total variation regularization, could be adopted to improve contrast and overcome the problem of data insufficiency [5][6][7]. Partially known support could further reduce the transducer element number [8]. Special sparse dictionary, which was a combination of different transforms, was used to capture features and suppress artifacts [9].
Recently, apart from the model-based methods, deep learning is introduced to eliminate undersampling artifacts. Antholzer et al. combined filtered back projection (BP) and U-net to reconstruct images from sparse data, with a quality comparable to state-of-the-art iterative approaches in numerical experiments [10]. This is a learning based post-processing method. Hauptmann et al. proposed a network to represent an iterative scheme and incorporated gradient information to suppress limited-view artifacts [11]. However, methods based on deep learning require large amount of data to train the network, and their robustness under in vivo conditions still needs to be verified.
Back projection (BP) based methods have the merits of low computational complexity and memory demand. Using multi-wavelength light and transducer displacement, some methods based on BP were developed for artifacts caused by acoustic reflection and out-of-plane signals [12,13]. As for artifact-suppression methods without a substantial increase of system complexity, Paltauf et al. designed special weights, but the algorithm only applies to the limited view problem where a 'detection region' is defined such that all boundaries of an object are visible by the detector array [14]. The delay multiply and sum method was suitable for point source imaging and was aimed at sidelobe suppression and resolution improvement, but was sensitive to high levels of noise [15,16]. Short-lag spatial coherence beamforming could enhance imaging contrast, at the expense of sacrificed lateral resolution [17,18].
In this paper, a reconstruction algorithm named contamination-tracing back-projection (CTBP) is proposed for the mitigation of streak-type artifacts. CTBP is a modified back projection method with high efficiency and good extensibility. In real practice, spoke-like streak artifacts always originate from bright features in the image. For instance, surface vessels or deep contrast agents with high absorptions and concentrations are candidates of such contamination sources. As a result of failed destructive interference due to a lack of information, such bright sources contaminate all points along a back-projection arc. In CTBP, special weight, which depends on the contamination degree of data, is added to minimize the negative influences of those contamination sources during back-projection reconstruction. Moreover, note that the accuracy of contamination tracing relies on the precision of the time of flight (TOF) estimation, and in our application acoustic speed heterogeneity correction based on iterative solution of the eikonal equation [19] is implemented for the purpose of accurate TOF calculation of contamination sources. Improvements of image quality can be observed in phantom and in vivo experiments, under the conditions of sparse sampling and limited view.

Methods
Back projection [20] is a commonly used method for PA initial pressure reconstruction. For heterogeneous media, the calculation of TOF should take into account the speed of sound (SOS) distribution. The formula of BP can be modified as in which dΩ 0 is the solid angle for a detector at ì r 0 with respect to a reconstruction point at ì r. Under practical conditions, the detected signal can be directly adopted as the back-projection term b(ì r 0 , t) [21]. Eikonal equation is a non-linear partial differential equation that tracks advancing wavefront [19]. The eikonal equation is where t is TOF and v is SOS. Iterative solution of Eq. (2) [22] can provide a TOF map [19], which is used in Eq. (1). Because acoustic ray paths are reversible, during reconstruction each detector is regarded as a source and the TOFs of the image pixels are calculated. All the results in this paper has taken into account SOS correction with Eq. (2). To illustrate the principle and effectiveness of CTBP. Figure 1 shows the results of a representative phantom experiment. We used spatially sparse sampling to make the streak artifacts more prominent (details of the system and the phantom will be given in the Experiments and Results section). Figure 1(a) is the direct reconstruction result using Eq. (1). In this case, the 'T' shaped feature was the main contamination source, which caused serious streak artifacts in the entire object. Figure 1(b) is the result of CTBP. Figures 1(c) and (d) illustrate the procedure of CTBP. In Fig. 1(c), the white points are identified as the contamination sources, which can be obtained with simple thresholding and appropriate image dilation of Fig. 1(a). The black curves are partial iso-TOF curves of a specific detector, here "Iso-TOF" means on these curves the time-of-flight values are the same. Since the initial pressure values along the same iso-TOF curve were integrated into a single temporal value, conventional reconstruction (i.e., applying Eq. (1)) assigns that detected value uniformly to the iso-TOF curve. The number of contamination sources along an iso-TOF curve reflects the possibility of contamination, which can be taken into account to form the back-projection weight. Thus, the proposed CTBP is formulated as where w (CT) (ì r, ì r i ) is the suppression weight used for minimizing the influences of contamination sources, and ì r i is the location of the detector. A simple but effective design of w (CT) (ì r, ì r i ) is where n(ì r, ì r i ) is the number of contamination sources along the corresponding iso-TOF curve, a is a decay parameter, and w min is the minimum weight. It should be mentioned that Eq. (4) is only adopted for points not identified as contamination sources. w (CT) of any contamination source is always 1. With Eq. (4), if the number of contamination sources is bigger, the weight is smaller. If there is no contamination point on an iso-TOF curve, the weight is 1. An example of the suppression weights for a single detector element is given in Fig. 1(d). With CTBP (Eq. (3)), the streak-type artifacts are significantly reduced ( Fig. 1(b)), yielding a final result ( Fig. 1(b)) well resembling the real phantom ( Fig. 1(e)). A more quantitative comparison is given in Fig. 1 where the top-to-bottom profile fluctuation along a line in Fig. 1(a) is plotted along with the one in Fig. 1(b). The profile of Fig. 1(a) has larger fluctuation than that of Fig. 1 The CTBP algorithm can be summarized as follows: Algorithm Contamination-tracing back-projection Set parameters Set the decay parameter a, and the minimum weight w min .
3. Adaptive suppression weights (Eq. (4)) are generated based on the number of contamination sources along each iso-TOF curve. Perform contamination-tracing back-projection (Eq. (3)). The calculation of the suppression weight (Eq. (4)) nearly costs no time, so the CTBP computational time is only twice as much as the traditional BP (Eq. (1)). Even if SOS heterogeneity correction is applied without graphics processing unit (GPU) acceleration, the whole computational time for the phantom experiment with sparse sampling is only 178.4 s (128 channels, Matlab R2016b, Intel(R) Core(TM) i7-7700 CPU @ 3.6 GHz). The computational time can be greatly reduced with parallel computation.

Experiments and results
A 256-element full-ring ultrasound detector array (Imasonic Inc.; 5.5 MHz central frequency; > 60% -6 dB bandwidth) was adopted in this study. Two data acquisition units (Analogic Corp., SonixDAQ; 40 MHz sampling rate; 12-bit ADC; 36-51 dB programmable gain) were employed for data acquisition. Details of the system can be found in [22]. Rotation of the array by 2π/512 doubles the sampling density. Partial or all samples of the 512-channel raw data were used for reconstruction. A 1,064 nm Nd:YAG laser (LOTIS, LS-2145-LT150, 10 Hz pulse repetition rate, 12 ns pulse duration) was employed for in-vivo mouse liver imaging with no other targets added. The maximum energy density on the sample (11 mJ/cm 2 )) was under the American National Standards Institute (ANSI) security limits (100 mJ/cm 2 at 1,064 nm, 10 Hz laser pulse repetition rate). The raw data of the normal liver experiments (Section 3.2) have been used in our previous paper [22]. An optical parametric oscillator (OPO) (SOLAR LP604, 680-1,064 nm tuning range, 10 Hz repetition, 10 ns temporal width) was employed for in-vivo mouse liver imaging with carbon rods (working at 1,050 nm). The maximum energy density on the sample (11mJ/cm 2 at 830nm) was below the ANSI limits. For all in vivo liver imaging experiments, light delivery was realized with a customized one-to-ten fiber bundle (Silica fiber, transmission band 500-2,400 nm).
For in-vivo brain imaging, the ten-branched fiber bundle was replaced by a single-ended one. The angle of the fiber head is continuously adjustable within a certain range in the vertical plane. The OPO was adopted and the light of 850 nm wavelength was used to balance optical penetration and energy density. The emergent light from the fiber head was delivered obliquely onto the mouse scalp. The laser fluence was approximately 50 mJ/cm 2 which was below the ANSI safety limits. For phantom experiments, the excitation wavelength was also 850 nm.

Phantom experiments
An agarose phantom was made ( Fig. 1(e)). The outer part of the phantom was composed of 4% w/w agar, 0.85% w/w diluted India ink (0.625% v/v), and 0.5% w/w intralipid. The dark features were made of 4% w/w agar, 66% w/w diluted ink and 0.5% w/w intralipid. The inner agarose rod (8.75 mm in diameter) was composed of 2.5% w/w agar and 0.5% w/w intralipid. Different SOS values were set for different parts of the phantom when reconstruction was carried out. Undersampling was performed by selecting 128 channels along the ring array uniformly, to generate the results shown in Fig. 1. More detailed discussion about the results of this phantom experiment was provided earlier in the Method section. Figure 2 shows the results corresponding to limited-view measurements. Figures 2(a)-(c) and Figs. 2(d)-(f) are the results due to 1/3 ring (120°) and full ring acquisitions, respectively. It should be emphasized that in our case, even the full ring measurement should be considered as imperfect: 1) the detector array had a limited view, because the detection solid angle was less than 4π; 2) the array itself was unideal, presenting uneven sensitivity among different elements, dead channels, and noise. The effect of imperfect measurement can be mitigated by moving average of the raw data, if the spatial sampling density is not too sparse. Figures 2(a)(d), (b)(e), and (c)(f) are the results of direct back-projection, back-projection after raw data moving average, and CTBP with moving average, respectively. After data filtering, streak artifacts can be suppressed slightly, while most of the obvious artifacts still remain. In contrast, most artifacts are effectively suppressed with CTBP (Figs. 2(c) and (f)). In the following section, moving average of raw data was carried out before CTBP reconstruction, if sampling density is not sparse, unless otherwise stated. The data filtering might make images a little blurred. There are shape distortions for the whole phantom in Figs. 2(a)-(c), due to the limited view problem. The edges, whose surface normal do not point to any detectors, tend to disappear in the reconstructed images. These kind of shape distortions are out of the scope of this paper and need to be studied in the future.

In vivo mouse liver imaging
The cross-section of the liver region of a normal nude mouse (8-week-old Crl: NU-Foxn 1 nu mouse) was imaged using the photoacoustic computed tomography (PACT) system. All animal experiments were conducted in conformity with the regulations of the Laboratory Animal Research Center at Tsinghua University, Beijing, China. Because the liver nearly occupied the entire cross-section, uniform SOS inside the mouse body was assumed. Figures 3(a)

Mouse liver imaging with exogenous absorbers
To further demonstrate the robustness of CTBP, carbon rods were added to the imaged scene. The design of these experiments mimicked the situations where contrast agents with strong absorption were used, or highly-absorptive melanomas existed. Two experiments were performed. In the first experiment, a normal nude mouse (6-week-old Crl: NU-Foxn 1 nu mouse) was scanned across its liver with a carbon rod (10 cm in length, 0.7 mm in diameter) placed nearby. Figure 4 shows the results. Figure 4(a) is the schematic diagram of the experimental setup. The carbon rod has strong absorption and it is closer to the illumination source, thus its initial pressure is obviously higher than that of the mouse. When conventional back-projection was carried out, streak artifacts deteriorated the image (Fig. 4(b)). After the rod was identified as the contamination source, CTBP successfully suppressed the streak patterns inside the imaged object (Fig. 4(c)). The corresponding image accorded with the image in the control experiment where no rods were added (Fig. 4(d)).
In the second experiment, a normal nude mouse (6-week-old Crl: NU-Foxn 1 nu mouse) was euthanized. A carbon rod (10 cm in length, 0.7 mm in diameter) was inserted into the mouse along the esophagus to the stomach. Then the mouse was scanned at its liver (near its lung). It is worth mentioning that in this case only 256 channel full-ring data were used to make the radial artifacts more serious. Figure 5(a) is a photograph of the mouse's cross-section corresponding to the imaged region. The artifacts caused by the rod in Fig. 5(b) can be effectively cancelled out with CTBP, as shown in Fig. 5(c). The result agreed well with the rod-free image in Fig. 5(d). Figures 4 and 5 were obtained at different axial positions along the mouse trunk, which explains why Fig. 4 shows more vascular features.

In vivo mouse brain imaging
A normal mouse (CD-1;10-week-old,30g body weight) was adopted for in-vivo brain imaging. Hair on the mouse scalp was first removed using depilatory cream avoiding acoustic scattering. 1% isoflurane anesthesia was used while the mouse head was mounted vertically in the center of transducer during the experiment. 512-channel data were acquired and the acoustic sectioning was along the coronal plane.
The results of brain imaging are shown in Fig. 6. As shown in Fig. 6(a), the brain was illuminated on one side. The blue arrow in Fig. 6(b) denotes the blood-rich part with high signal level. The inner region of the brain has weak signals, due to optical fluence attenuation, and the corresponding features have low contrast. Acoustic impedance mismatch at the skull further attenuated the intracranial signals. Under these conditions, the radial artifacts due to the strong absorber (blue arrow in Fig. 6(b)) can adversely affect any precise analysis of the deep brain. Dramatic suppression of the artifacts is demonstrated using CTBP, as shown in Fig. 6(c). This particular example demonstrates the capability of CTBP in PACT based brain imaging applications. CTBP combined with multispectral imaging can be used for functional and molecular imaging of the brain. Potential applications in whole-brain calcium or voltage imaging in rodents make PACT a powerful tool for neural science, however, the weak functional signals are highly susceptible to noises and artifacts. The apparent reduction of the artifacts shown in Fig. 6(c) shows great potential to advance PACT toward practical neural imaging applications.
In the above demonstrations, to implement CTBP, the parameters "a" and "w min " in Eq. (4) are usually set within 0.2-0.4 and 0-0.3, respectively. If streak artifacts are serious, "a" with high value and "w min " with low value should be adopted. For example, in the mouse experiment with carbon rod near the fur, "a" and "w min " are 0.4 and 0, respectively.

Quantitative evaluation of the proposed method
Numerical experiments were carried out to quantitatively evaluate the performance of CTBP. In the simulations, a 512-element ring array with a radius of 5 cm was used, while some of the elements were selectively disabled in various testing conditions. The SOS values of the target (2.85 cm in diameter) and water were 1,560 m/s and 1,480 m/s, respectively. Inside the imaged object, there were eight strong absorbers and two ellipses. Sparse sampling and partial array (limited view) were studied quantitatively. To ensure fairness, data moving average was not used, which means that the improvement of image quality was only due to CTBP. An indicator M G based on image gradient is defined to quantitatively evaluate the severity of artifacts, in which ∇I i is the gradient of the ith image pixel's value. We chose a fixed feature-free region in the reconstructed image (highlighted by the red box in Fig. 7(a)) to calculate M G . The results are given in Table 1

Discussion and conclusion
Insufficient detectors often lead to serious streak artifacts in PAT, especially when strong absorbers exist. Compared to deep learning or model based reconstruction methods, the proposed CTBP algorithm is a back-projection-based method, which is free of heavy data training, excessive computational burden and memory demand. Moreover, SOS heterogeneity correction is applied to guarantee good image quality and accurate tracing of the iso-TOF curves. Various phantom and mice experiments, under the conditions of sparse spatial sampling or limited detection view, were carried out. It is demonstrated that CTBP can effectively suppress the streak artifacts and recover the real details. If there are strong exogenous absorbers (such as those shown in Fig. 4), streak artifacts can deteriorate image quality significantly. The problem can be effectively mitigated by CTBP. In brain imaging (Fig. 6), deep region usually has weak signal, due to light fluence attenuation and acoustic impedance mismatch. Under these conditions, bright superficial features (such as the superior longitudinal sinuses) can contaminate the deep regions, resulting in failure of quantitative analysis. Therefore, CTBP is envisioned to help photoacoustic brain imaging (as demonstrated in Fig. 6), especially when delicate functional and molecular studies are concerned. The determination of the contamination sources, i.e., the strong absorbers, is based on an image reconstructed using BP. If the amplitudes of the sources are not significantly higher than the rest of the regions, simple thresholding is not effective to identify them. If the threshold is set too high, some contamination sources will be missed; however, if the threshold is set too low, false positives start to appear. So in this case, manual selection has to be carried out to remove them. In the future, automatic selection method for contamination sources will be studied. Besides, the limit view problem can cause edge disappearance (Fig. 2), which is out of the scope of this paper but needs to be studied in the future. Some parts of the edge of the inner rod in the phantom experiments were also suppressed together with the artifacts, posing another problem to be addressed in the future. Moreover, the strength of the contamination sources, together with the number of them, can be taken into account in the future.
In conclusion, CTBP is an effective method for suppression of streak artifacts in PACT. The proposed algorithm can also be used to other measurement geometries, such as linear arrays, planar arrays, and so on. As such, we envision that PAT image qualities will benefit from this simple, effective and generic method.

Funding
National Natural Science Foundation of China (61735016, 61871251); Tsinghua National Laboratory for Information Science and Technology (Youth Innovation Fund).

Disclosures
The authors declare that there are no conflicts of interest related to this article.