Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Streak artifact suppression in photoacoustic computed tomography using adaptive back projection

Open Access Open Access

Abstract

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.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. 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 L1 norm of the image in sparse representation and total variation regularization, could be adopted to improve contrast and overcome the problem of data insufficiency [57]. 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.

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

$$p_0^{(b)}(\vec{r}) = \int_{{\Omega _0}} {b({{\vec{r}}_0},t = TOF(\vec{r},{{\vec{r}}_0}))d{\Omega _0}/{\Omega _0}} ,$$
in which $d{{\Omega }_0}$ is the solid angle for a detector at ${\vec{r}_0}$ with respect to a reconstruction point at $\vec{r}$. Under practical conditions, the detected signal can be directly adopted as the back-projection term $b({{{\vec{r}}_0},\textrm{t}} )$ [21]. Eikonal equation is a non-linear partial differential equation that tracks advancing wavefront [19]. The eikonal equation is
$$|{\nabla t(\vec{r})} |= \frac{1}{{v(\vec{r})}},$$
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

$$p_0^{(CT)}(\vec{r}) = \frac{1}{{{\Omega _0}}}\sum\limits_i {{w^{(CT)}}(\vec{r},{{\vec{r}}_i})b({{\vec{r}}_i},t = TOF(\vec{r},{{\vec{r}}_i}))\Delta {\Omega _0}(\vec{r},{{\vec{r}}_i})} ,$$
where ${w^{({CT} )}}({\vec{r},{{\vec{r}}_i}} )$ is the suppression weight used for minimizing the influences of contamination sources, and ${\vec{r}_i}$ is the location of the detector. A simple but effective design of ${w^{({CT} )}}({\vec{r},{{\vec{r}}_i}} )$ is
$${w^{(CT)}}(\vec{r},{\vec{r}_i}) = (1 - {w_{\min }}){e^{ - a \cdot n(\vec{r},{{\vec{r}}_i})}} + {w_{\min }},$$
where $n({\vec{r},{{\vec{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(f), 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(b).

 figure: Fig. 1.

Fig. 1. Phantom experiment with sparse sampling. (a) The direct reconstruction result and (b) the result of CTBP, under the condition of sparse sampling density. (c) The TOF calculation. The black curves denote the iso-TOF curves of a specific detector. White points are the contamination sources, which are generated with (a). (d) The weight map used in CTBP. (e) The photograph of phantom. (f) The reconstructed profiles correspond to the dotted lines in (a) and (b) (from the top to bottom). Scale bar: 5 mm.

Download Full Size | PDF

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}}$.

Run CTBP

  • 1. Rough reconstruction of initial pressure $p_0^{(b )}$ with Eq. (1).
  • 2. Extraction of pixels with high initial pressure, i.e., contamination sources.
  • 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 Core i7-7700 CPU @ 3.6 GHz). The computational time can be greatly reduced with parallel computation.

3. 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/cm2)) was under the American National Standards Institute (ANSI) security limits (100 mJ/cm2 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 ($11\textrm{mJ}/\textrm{c}{\textrm{m}^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 \,{\textrm{mJ}}/\textrm{c}{\textrm{m}^2}$ which was below the ANSI safety limits. For phantom experiments, the excitation wavelength was also 850 nm.

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

 figure: Fig. 2.

Fig. 2. Phantom experiment with limited view. (a)–(c) and (d)–(f) are the results of 1/3 ring and full ring detection, respectively. (a)(d), (b)(e), and (c)(f) are the results of direct back-projection, back-projection after raw data moving average, and CTBP with data moving average, respectively. Scale bar: 5 mm.

Download Full Size | PDF

3.2 In vivo mouse liver imaging

The cross-section of the liver region of a normal nude mouse (8-week-old Crl: NU-Foxn 1nu 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)(b) and Figs. 3(c)(d) are the results of 1/3 ring and full ring acquisitions, respectively. Due to light fluence attenuation along the depth direction, superficial features usually have stronger initial pressure. Therefore, in the results of direct back-projection (Figs. 3(a) and (c)), the vertical surface vessels (i.e., bright dots on the periphery of the body) have high contrast. If their signals could not be averaged out elsewhere along the back-projection paths, obvious streak artifacts appear near them. Most of those artifacts are suppressed with CTBP (Figs. 3(b) and (d)).

 figure: Fig. 3.

Fig. 3. In vivo imaging of the mouse liver cross-section. (a)(b) and (c)(d) are the results of 1/3 ring and full ring acquisitions, respectively. (a)(c) and (b)(d) are the results of direct back-projection and CTBP, respectively. Scale bar: 5 mm.

Download Full Size | PDF

3.3 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 1nu 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)).

 figure: Fig. 4.

Fig. 4. Imaging of the mouse liver region with external carbon rod. (a) The schematic diagram of experimental setup. OF: optical fiber; UA: ultrasound array; WT: water; CR: carbon rod. (b) The result of direct back-projection. The blue arrow denotes the rod. (c) The result of CTBP. (d) Back-projection result of the corresponding slice without the rod, as gold standard. Scale bar: 5 mm.

Download Full Size | PDF

In the second experiment, a normal nude mouse (6-week-old Crl: NU-Foxn 1nu 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.

 figure: Fig. 5.

Fig. 5. Imaging of the mouse liver region (near the lungs) with inserted carbon rod. (a) Photo of the frozen section corresponding to the region imaged by PA. IVC: inferior vena cava; LV: liver; CR: carbon rod; AA: abdominal aorta; L: lung; SC: spinal cord. (b) The result of direct back-projection. The blue arrow denotes the rod. (c) The result of CTBP. (d) Back-projection result of the corresponding slice without the rod, as gold standard. Scale bar: 5 mm.

Download Full Size | PDF

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

 figure: Fig. 6.

Fig. 6. In vivo mouse brain imaging along the coronal plane. (a) Experimental setup. OF: optical fiber; AF: acoustic focus; UA: ultrasound array; WT: water. (b) The result of direct back-projection. (c) The result of CTBP. The blue arrow denotes the part with strong absorption. Scale bar: 5 mm.

Download Full Size | PDF

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.

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

In Fig. 7, the reconstruction results of direct BP (Figs. 7(a)–(d)) and CTBP (Figs. 7(e)–(h)) are given. Figures 7(a), (b), (e), and (f) are the results of sparse sampling, where detecting elements are evenly distributed along the ring transducer array. Figures 7(a) and (e) correspond to 128 elements, while Figs. 7(b) and (f) correspond to 64. Sparse sampling results in streak artifacts with uniform distribution in this numerical phantom. Artifacts appear throughout the whole object in Figs. 7(a) and (b). With CTBP, the artifacts can be suppressed (Figs. 7(e) and (f)). If the number of elements is too small, weak background fluctuations still exist but at a dramatically reduced amplitude (Fig. 7(f)). Figures 7(c), (d), (g) and (h) show the reconstruction results where partial arrays were employed—1/2 ring (180° angular spread, 256 consecutive elements) for Figs. 7(c) and (g) and 1/3 ring (120° angular spread, 170 consecutive elements) for Figs. 7(d) and (h). The artifacts in these scenarios appear drastically different from those in sparse sampling, and are mainly distributed close to the contamination sources (Figs. 7(c) and (d)). They can also be effectively suppressed with CTBP (Figs. 7(g) and (h)).

 figure: Fig. 7.

Fig. 7. Results of numerical experiments. (a)–(d) and (e)–(h) are the results of direct back-projection and CTBP, respectively. (a)(b)(e)(f) are the results of sparse sampling, where elements are evenly distributed on a full ring array. The total number of elements for (a)(e) and (b)(f) are 128 and 64, respectively. (c)(d)(g)(h) are the results of limited view, where (c)(g) and (d)(h) are the results of 1/2 ring and 1/3 ring, respectively. Red box denotes the part used to perform quantitative evaluation. Scale bar: 5 mm.

Download Full Size | PDF

An indicator ${M_G}$ based on image gradient is defined to quantitatively evaluate the severity of artifacts,

$${M_G} = \sum\limits_i {||\nabla {I_i}|{|_2}} ,$$
in which $\nabla {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. For both sparse sampling and limited view, the fewer elements are used, the higher the image gradient magnitude (${M_G}$) can be obtained. Compared to those obtained by standard BP, the values of MG are nearly reduced by half if CTBP is used. CTBP performs especially well when sampling is made very sparse. When the number of elements is 32, the ${M_G}$ value of CTBP is only 40.7% of that of BP.

Tables Icon

Table 1. Quantitative evaluation of image gradient (MG)

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

References

1. L. V. Wang, “Multiscale photoacoustic microscopy and computed tomography,” Nat. Photonics 3(9), 503–509 (2009). [CrossRef]  

2. L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science 335(6075), 1458–1462 (2012). [CrossRef]  

3. J. F. Barrett and N. Keat, “Artifacts in CT: recognition and avoidance,” RadioGraphics 24(6), 1679–1691 (2004). [CrossRef]  

4. X. L. Deán-Ben, A. Buehler, V. Ntziachristos, and D. Razansky, “Accurate model-based reconstruction algorithm for three-dimensional optoacoustic tomography,” IEEE Trans. Med. Imag. 31(10), 1922–1928 (2012). [CrossRef]  

5. Y. Han, S. Tzoumas, A. Nunes, and V. Ntziachristos, “Sparsity-based acoustic inversion in cross-sectional multiscale optoacoustic imaging,” Med. Phys. 42(9), 5444–5452 (2015). [CrossRef]  

6. Y. Han, L. Ding, X. L. Deán-Ben, D. Razansky, J. Prakash, and V. Ntziachristos, “Three-dimensional optoacoustic reconstruction using fast sparse representation,” Opt. Lett. 42(5), 979–982 (2017). [CrossRef]  

7. Y. Zhang, Y. Wang, and C. Zhang, “Total variation based gradient descent algorithm for sparse-view photoacoustic image reconstruction,” Ultrasonics 52(8), 1046–1055 (2012). [CrossRef]  

8. J. Meng, L. V. Wang, L. Ying, D. Liang, and L. Song, “Compressed-sensing photoacoustic computed tomography in vivo with partially known support,” Opt. Express 20(15), 16510–16523 (2012). [CrossRef]  

9. P. Omidi, M. Zafar, M. Mozaffarzadeh, A. Hariri, X. Haung, M. Orooji, and M. Nasiriavanaki, “A novel dictionary-based image reconstruction for photoacoustic computed tomography,” Appl. Sci. 8(9), 1570 (2018). [CrossRef]  

10. S. Antholzer, M. Haltmeier, and J. Schwab, “Deep learning for photoacoustic tomography from sparse data,” Inverse Probl. Sci. Eng. 27(7), 987–1005 (2019). [CrossRef]  

11. F. Hauptmann, M. Lucka, N. Betcke, J. Huynh, B. Adler, P. Cox, S. Beard, S. Ourselin, and Arridge, “Model-based learning for accelerated, limited-view 3-d photoacoustic tomography,” IEEE Trans. Med. Imag. 37(6), 1382–1393 (2018). [CrossRef]  

12. H. N. Y. Nguyen, A. Hussain, and W. Steenbergen, “Reflection artifact identification in photoacoustic imaging using multi-wavelength excitation,” Biomed. Opt. Express 9(10), 4613–4630 (2018). [CrossRef]  

13. H. N. Y. Nguyen and W. Steenbergen, “Reducing artifacts in photoacoustic imaging by using multi-wavelength excitation and transducer displacement,” Biomed. Opt. Express 10(7), 3124–3138 (2019). [CrossRef]  

14. G. Paltauf, R. Nuster, and P. Burgholzer, “Weight factors for limited angle photoacoustic tomography,” Phys. Med. Biol. 54(11), 3303–3314 (2009). [CrossRef]  

15. G. Matrone, A. S. Savoia, G. Caliano, and G. Magenes, “The delay multiply and sum beamforming algorithm in ultrasound B-mode medical imaging,” IEEE Trans. Med. Imag. 34(4), 940–949 (2015). [CrossRef]  

16. M. Mozaffarzadeh, A. Mahloojifar, M. Orooji, S. Adabi, and M. Nasiriavanaki, “Double-stage delay multiply and sum beamforming algorithm: Application to linear-array photoacoustic imaging,” IEEE Trans. Biomed. Eng. 65(1), 31–42 (2018). [CrossRef]  

17. M. A. L. Bell, N. Kuo, D. Y. Song, and E. M. Boctor, “Short-lag spatial coherence beamforming of photoacoustic images for enhanced visualization of prostate brachytherapy seeds,” Biomed. Opt. Express 4(10), 1964–1977 (2013). [CrossRef]  

18. M. A. Lediju, G. E. Trahey, B. C. Byram, and J. J. Dahl, “Short-lag spatial coherence of backscattered echoes: Imaging characteristics,” IEEE Trans. Sonics Ultrason. 58(7), 1377–1388 (2011). [CrossRef]  

19. M. S. Hassouna and A. A. Farag, “Multi-stencils fast marching methods: A highly accurate solution to the eikonal equation on cartesian domains,” IEEE Trans. Pattern Anal. Mach. Intell. 29(9), 1563–1574 (2007). [CrossRef]  

20. M. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E 71(1), 016706 (2005). [CrossRef]  

21. M. Pramanik, “Improving tangential resolution with a modified delay-and-sum reconstruction algorithm in photoacoustic and thermoacoustic tomography,” J. Opt. Soc. Am. A 31(3), 621–627 (2014). [CrossRef]  

22. C. Cai, X. Wang, K. Si, J. Qian, J. Luo, and C. Ma, “Feature coupling photoacoustic computed tomography for joint reconstruction of initial pressure and sound speed in vivo,” Biomed. Opt. Express 10(7), 3447–3462 (2019). [CrossRef]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (7)

Fig. 1.
Fig. 1. Phantom experiment with sparse sampling. (a) The direct reconstruction result and (b) the result of CTBP, under the condition of sparse sampling density. (c) The TOF calculation. The black curves denote the iso-TOF curves of a specific detector. White points are the contamination sources, which are generated with (a). (d) The weight map used in CTBP. (e) The photograph of phantom. (f) The reconstructed profiles correspond to the dotted lines in (a) and (b) (from the top to bottom). Scale bar: 5 mm.
Fig. 2.
Fig. 2. Phantom experiment with limited view. (a)–(c) and (d)–(f) are the results of 1/3 ring and full ring detection, respectively. (a)(d), (b)(e), and (c)(f) are the results of direct back-projection, back-projection after raw data moving average, and CTBP with data moving average, respectively. Scale bar: 5 mm.
Fig. 3.
Fig. 3. In vivo imaging of the mouse liver cross-section. (a)(b) and (c)(d) are the results of 1/3 ring and full ring acquisitions, respectively. (a)(c) and (b)(d) are the results of direct back-projection and CTBP, respectively. Scale bar: 5 mm.
Fig. 4.
Fig. 4. Imaging of the mouse liver region with external carbon rod. (a) The schematic diagram of experimental setup. OF: optical fiber; UA: ultrasound array; WT: water; CR: carbon rod. (b) The result of direct back-projection. The blue arrow denotes the rod. (c) The result of CTBP. (d) Back-projection result of the corresponding slice without the rod, as gold standard. Scale bar: 5 mm.
Fig. 5.
Fig. 5. Imaging of the mouse liver region (near the lungs) with inserted carbon rod. (a) Photo of the frozen section corresponding to the region imaged by PA. IVC: inferior vena cava; LV: liver; CR: carbon rod; AA: abdominal aorta; L: lung; SC: spinal cord. (b) The result of direct back-projection. The blue arrow denotes the rod. (c) The result of CTBP. (d) Back-projection result of the corresponding slice without the rod, as gold standard. Scale bar: 5 mm.
Fig. 6.
Fig. 6. In vivo mouse brain imaging along the coronal plane. (a) Experimental setup. OF: optical fiber; AF: acoustic focus; UA: ultrasound array; WT: water. (b) The result of direct back-projection. (c) The result of CTBP. The blue arrow denotes the part with strong absorption. Scale bar: 5 mm.
Fig. 7.
Fig. 7. Results of numerical experiments. (a)–(d) and (e)–(h) are the results of direct back-projection and CTBP, respectively. (a)(b)(e)(f) are the results of sparse sampling, where elements are evenly distributed on a full ring array. The total number of elements for (a)(e) and (b)(f) are 128 and 64, respectively. (c)(d)(g)(h) are the results of limited view, where (c)(g) and (d)(h) are the results of 1/2 ring and 1/3 ring, respectively. Red box denotes the part used to perform quantitative evaluation. Scale bar: 5 mm.

Tables (1)

Tables Icon

Table 1. Quantitative evaluation of image gradient (MG)

Equations (5)

Equations on this page are rendered with MathJax. Learn more.

p 0 ( b ) ( r ) = Ω 0 b ( r 0 , t = T O F ( r , r 0 ) ) d Ω 0 / Ω 0 ,
| t ( r ) | = 1 v ( r ) ,
p 0 ( C T ) ( r ) = 1 Ω 0 i w ( C T ) ( r , r i ) b ( r i , t = T O F ( r , r i ) ) Δ Ω 0 ( r , r i ) ,
w ( C T ) ( r , r i ) = ( 1 w min ) e a n ( r , r i ) + w min ,
M G = i | | I i | | 2 ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.