Image improvement in linear-array photoacoustic imaging using high resolution coherence factor weighting technique

Abstract

Background

In Photoacoustic imaging (PAI), the most prevalent beamforming algorithm is delay-and-sum (DAS) due to its simple implementation. However, it results in a low quality image affected by the high level of sidelobes. Coherence factor (CF) can be used to address the sidelobes in the reconstructed images by DAS, but the resolution improvement is not good enough, compared to the high resolution beamformers such as minimum variance (MV). In this paper, it is proposed to use high-resolution-CF (HRCF) weighting technique in which MV is used instead of the existing DAS in the formula of the conventional CF.

Results

The higher performance of HRCF is proved numerically and experimentally. The quantitative results obtained with the simulations show that at the depth of 40 mm, in comparison with DAS+CF and MV+CF, HRCF improves the full-width-half-maximum of about 91% and 15% and the signal-to-noise ratio about 40% and 14%, respectively.

Conclusion

Proposed method provides a high resolution along with a low level of sidelobes for PAI.

Background

In photoacoustic imaging (PAI), a short electromagnetic pulse, i.e. laser or radio frequency (RF), illuminates the target of imaging, and Ultrasound (US) waves are generated based on the thermoelastic effect [1, 2]. In comparison with other imaging modalities, PAI has multiple advantages leading to many investigations [3, 4]. The main incentive in PAI is having the merits of the US imaging spatial resolution and the optical imaging contrast in one imaging modality [5]. PAI can be used in different fields of study such as tumor detection [6, 7], ocular imaging [8] and functional imaging [9, 10]. Moreover, contrast agents and nanoparticles play a significant role in PAI [11, 12]. PAI can be separated into two fields: photoacoustic tomography (PAT) and photoacoustic microscopy (PAM) [13, 14]. PAT, for the first time, was successfully used as in vivo functional and structural brain imaging modality in small animals [15]. In PAT, an array of elements may be formed in linear, arc or circular shape, and mathematical reconstruction algorithms are used to obtain the optical absorption distribution map of the tissue [16–18]. Most of the used reconstruction algorithms for image formation in PAI are based on some assumptions leading to artifacts and disturbing effects on the formed photoacoustic (PA) images. One of the challenges in PA image formation is related to reduction of these effects for different number of transducers and properties of imaging media [19–21].

Some modifications should be considered if an algorithm in US imaging is going to be used in PAI. These modifications have led using different hardware to implement an integrated US-PA imaging device [22, 23]. DAS is the most commonly used beamforming algorithm in PAI. However, it leads to a low quality image, having a wide mainlobe and high level of sidelobes [24]. Adaptive beamformers, commonly employed in radar, have the ability of weighting the aperture based on the characteristics of detected signals, providing a high quality image with a wide range of off-axis signals rejection. MV can be treated as one of the commonly used adaptive methods in medical US imaging [25, 26]. Vast variety of modifications have been investigated on MV such as complexity reduction [27, 28], shadowing suppression [29], using eigenstructure to enhance MV performance [30, 31], and combination of MV and multi-line transmission (MLT) technique [32]. Matrone et al.proposed a new algorithm namely delay-multiply-and-sum (DMAS), as a beamforming technique for medical US imaging [33]. Double stage DMAS (DS-DMAS), outperforming DMAS in the terms of contrast and sidelobes, was introduced for the linear-array US and PAI [34–36]. Minimum variance-based DMAS has been proposed for resolution improvement in DMAS while the level of sidelobes would be retained [37, 38]. Coherence factor (CF) can be mentioned as one of the prevalent weighting methods in beamforming field [39]. The performance of CF has been investigated for US imaging and PAI in [40] and [41], respectively. Short-lag spatial coherence beamforming was also used to enhance the visualization of prostate brachytherapy seeds [42, 43]. Moreover, a high resolution CF (HRCF) has been investigated for high-frame rate US imaging [44]. Recently, a modified version of the CF has been reported by the authors where the aim was to achieve a higher contrast compared to the conventional CF [45].

In this paper, the performance of HRCF is investigated for linear-array PAI. The concept of this technique indicates that a high resolution image, obtained with an algorithm such as MV, can be used to weight the calculated samples instead of the formed image by DAS. It is shown that the proposed weighting algorithm (used with DAS) outperforms the DAS and MV (with/without CF) in the terms of resolution, sidelobes and contrast.

Numerical results and performance assessment

In this section, numerical results are presented to illustrate the performance of the proposed technique for PA image formation in comparison with DAS, DAS+CF, MV and MV+CF.

Simulated point targets

Simulation setup

The K-wave Matlab toolbox was used to simulate the numerical study [46]. Imaging region was 20 mm in the lateral axis and 80 mm in the vertical axis. A linear-array having M=128 elements operating at 7 MHz central frequency and 77% fractional bandwidth was used to detect the PA signals generated from the defined initial pressures. The schematic of the designed simulation is shown in Fig. 1. The speed of sound was assumed to be 1540 m/s during the simulations. The sampling frequency was 50 MHz, subarray length L=M/2, K=3 and Δ=1/100L for all the simulations.

Fig. 1

Schematic of the simulation study

Qualitative and quantitative evaluation

The reconstructed images are shown in Fig. 2, along with a zoomed version at the depth of 40 mm (shown in Fig. 3) for a better evaluation. As can be seen, the reconstructed image using DAS have a low quality along with high sidelobes. MV improves the resolution significantly, but the sidelobes still affect the image. Using CF combined with DAS or MV results in sidelobes reduction and image quality enhancement. Even though the image reconstructed by MV+CF, shown in Fig. 2d, has a high resolution, but the negative effects of the sidelobes still degrade the image quality. In Fig. 2e, it can be seen that the sidelobes are reduced compared to Fig. 2d while the resolution is retained. To assess in more details, the lateral variations of the reconstructed images shown in Fig. 2 are shown at four imaging depths in Fig. 4. As can be seen, DAS+HRCF results in lower level of sidelobes and narrower width of mainlobe compared to other beamformers. Moreover, the lateral valleys between the targets have the lowest levels using the proposed method. Consider, for instance, the depth of 50 mm where the level of sidelobes are of about -36 dB, -69 dB, -45 dB, -79 dB and -88 dB for DAS, DAS+CF, MV, MV+CF and DAS+HRCF, respectively. Thus, the proposed method leads to lowest sidelobes in comparison with other beamformers. Moreover, the levels of the lateral valleys for DAS, DAS+CF, MV, MV+CF and DAS+HRCF are about -29 dB, -61 dB, -37 dB, -70 dB and -80 dB, respectively. It indicates the higher separability of the proposed method. To evaluate the proposed method quantitatively, the full-width-half-maximum (FWHM) in -6 dB and signal-to-noise ratio (SNR) metrics are calculated and presented in Table 1 and Table 2, respectively. SNRs are calculated using the formula explained in [35]. For the axial FWHM, at the depth of 25 mm, DAS, DAS+CF, MV, MV+CF and DAS+HRCF leads to 442.7 μm, 441.0 μm, 433.2 μm, 431.7 μm and 425.1 μm, respectively.

Fig. 2

Reconstructed images using the simulated data. a DAS, b DAS+CF, c MV, d MV+CF and e DAS+HRCF. A linear-array and point phantom were used for the numerical design. All images are shown with a dynamic range of 60 dB. Noise was added to the detected signals having a SNR of 40 dB

Fig. 3

A close view of the reconstructed images shown in Fig. 1

Fig. 4

Lateral variations of the reconstructed images shown in Fig. 2 at the depths of a 20 mm, b 40 mm, c 55 mm and d 70 mm

Table 1 FWHM(μm) in -6 dB values at the different depths

Table 2 SNR (dB) values at the different depths

As demonstrated in Table 1, the proposed method for PA image reconstruction results in a narrower width of mainlobe in -6 dB compared to other beamformers in the all depth of imaging. Of note, there is no significant improvement compared to MV and MV+CF. Consider, in particular, the depth of 45 mm where DAS, DAS+CF, MV, MV+CF and DAS+HRCF leads to a FWHM of 2284 μm, 1388 μm, 131 μm, 131 μm and 103 μm, respectively. In comparison with a high resolution method such as MV, the proposed method leads to 28 μm FWHM improvement. As shown in Table 2, the proposed method results in a higher SNR in comparison with other reconstruction methods at the all depths of imaging. Consider, for instance, the depth of 55 mm where DAS, DAS+CF, MV, MV+CF and DAS+HRCF results in a SNR of 37.8 dB, 55.4 dB, 47.8 dB, 68.3 dB and 77.7 dB, respectively. In other words, DAS+HRCF improves the SNR for about 9 dB and 22 dB compared to MV+CF and DAS+CF, respectively, proving its superiority for linear-array PAI.

The proposed method is evaluated at the presence of high level of imaging noise. Eleven 0.1 mm radius spherical absorbers as initial pressure were positioned along the vertical axis every 5 mm beginning 25 mm from transducer surface. Noise was added to the detected signals having a SNR of 0 dB. The reconstructed images are shown in Fig. 5. As can be seen, the presence of the noise in the reconstructed images using DAS and MV degrade the images. CF results in the higher noise suppression and higher image quality, as shown in Fig. 5b and Fig. 5d. As shown in Fig. 5e, the sidelobes are better reduced using DAS+HRCF. The lateral variations for images shown in Fig. 5, at two depths of imaging, are shown in Fig. 6. As can be seen, the proposed method results in lower level of sidelobes and narrower width of mainlobe.

Fig. 5

Reconstructed images using the simulated data. a DAS, b DAS+CF, c MV, d MV+CF and e DAS+HRCF. A linear-array and point phantom were used for the numerical design. All images are shown with a dynamic range of 60 dB. Noise was added to the detected signals having a SNR of 0 dB

Fig. 6

Lateral variations of the reconstructed images shown in Fig. 5 at the depths of (a) 45 mm and b 55 mm

Experimental results

Experimental setup

To further evaluate the proposed weighting method and its effects on enhancing the PA images, phantom experiments were performed in which a phantom consists of 2 light absorbing wires with diameter of 150 μm were placed 1 mm apart from each other, as seen in Fig. 7. In this experiment, we utilized a Nd:YAG pulsed laser (Phocus core system, OPOTEK Inc, Carlsbad, CA, USA), with the pulse repletion rate of 30 Hz at wavelengths of 532 nm. A programmable digital ultrasound scanner (System Vantage 128,Verasonics, Inc., Kirkland, WA, USA), equipped with a linear array transducer (L11-4v, Verasonics, Inc., Kirkland, WA, USA) operating at frequency range between 4 to 9 MHz was utilized to acquire the PA RF data. A high speed FPGA was used to synchronize the light excitation and PA signal acquisition.

Fig. 7

The schematic of the setup used for the experimental PAI

Qualitative and quantitative evaluation

The reconstructed images are shown in Fig. 8. As can be seen, the artifact and noise affect the reconstructed image by DAS while the CF improves the image quality by suppressing them. As shown in Fig. 8c, MV results in an image having a high resolution, but the presence of the noise highly affects the image. As can be seen in Fig. 8e, HRCF results in a high resolution while the sidelobes are degraded, and the presence of the noise is clearly lower than other methods, comparing the background of the Fig. 8e with other images shown in Fig. 8. To assess the images in details, the lateral variations of the two wire targets are shown in Fig. 9. As can be seen, the HRCF outperforms the conventional CF combined with DAS and MV and results in a narrower width of mainlobe and lower level of sidelobes. Consider, for instance, the depth of 24 mm where DAS+CF, MV+CF and DAS+HRCF result in -36 dB, -47 dB and -60 dB sidelobes, respectively. In other words, DAS+HRCF improves the sidelobes for about 24 dB and 13 dB compared to DAS+CF and MV+CF, respectively. To compare the experimental images quantitatively, SNRs for all the methods are calculated and presented in Table 3 where the proposed weighting method leads to a higher SNR, for both imaging targets, compared to other methods, indicating the superiority of HRCF weighting method.

Fig. 8

Reconstructed images using the experimental detected data. a DAS, b DAS+CF, c MV, d MV+CF and e DAS+HRCF. A linear-array and wire target phantom were used for the experimental design. All images are shown with a dynamic range of 60 dB

Fig. 9

Lateral variations of the reconstructed images obtained with DAS+CF, MV+CF and DAS+HRCF, shown in Fig 8, at the depths of a 22 mm and b 24 mm

Table 3 SNR (dB) values at the different depths for images shown in Fig. 7

Discussion

The main improvement gained by HRCF is having a high resolution and low sidelobes at the same time. DAS is the most commonly used beamformer in PA and US imaging which is mainly as a result of its simple implementation. Moreover, it provides a real-time imaging. However, it results in a low quality image having a low resolution and high sidelobes due to its blindness and non-adaptiveness. To put it more simply, DAS considers all the calculated samples the same as each other, and there is just a summation process. On the other hand, adaptive beamformers, such as MV, provides a higher image quality compared to DAS, especially in the term of resolution. However, in MV, sidelobes affect the reconstructed image and degrade the image quality. CF is a weighting method that can be used with beamformers, such as DAS or MV, for sidelobes reduction. However, conventional CF weighting does not improve the resolution and the width of mainlobe significantly, compared to beamformers such as MV. It can be seen that in (2), the numerator of the formula of CF is the output of DAS. While CF reduces the sidelobes, the performance of CF is not high in the term of resolution, which is mainly due to the existence of DAS on the numerator of the formula of CF. Using MV instead of the exiting DAS in the (2) can improve the resolution gained by the conventional CF (15). The proposed method, HRCF, is a weighting method which can be applied on any beamforming algorithm (DAS was used in this paper). The reconstructed images (Fig. 2 and Fig. 3) show that the HRCF outperforms CF combined with DAS and MV. As shown in Fig. 2 and Fig. 3, the point targets are better distinguished and detectable using HRCF weighting procedure, and the sidelobes are better reduced. The proposed method was evaluated in the term of the presence of high level of noise, and the reconstructed images were shown in Fig. 5. As can be seen, the HRCF reduces the negative effects of the added noise, and it provides a higher robustness compared to other methods. The images have been evaluated using the lateral variations shown in Fig. 4 and Fig. 6, and all the results indicate the superiority of HRCF in the terms of sidelobes, lateral valley and the width of mainlobe. Tables 1, 2 and 3 show the quantitative evaluation of the proposed weighting method. They indicate that the HRCF reduces the presence of the noise and results in the narrower width of mainlobe. Despite all the evaluation with the simulations, the algorithm should be evaluated using experimental data. The generated experimental images are shown in Fig. 7, and the superiority of HRCF can be clearly seen. The lateral variations of the experimental images are shown in Fig. 9, proving the higher performance of HRCF. Compared to similar method proposed in [45], modified CF (MCF), it should be noticed that the proposed method in this paper focuses on a high resolution (better than MV), providing a contrast slightly better than CF. However, MCF provides almost 50 dB lower sidelobes compared to CF, with a slightly improved resolution. It should be mentioned that the higher performance of the HRCF is obtained at the expense of the higher computational burden where replacing DAS by MV on the numerator of the formula of CF would result the order of complexity change from O(M) to O(L³). The computational time to generate images shown in Fig. 7a-e was 0.96 s, 1.09 s, 144.16 s, 145.05 s and 145.11 s, respectively. All the results indicate that the HRCF can be an effective weighting method for image formation in linear-array PAI, and it provides a higher contrast and resolution compared to DAS and MV combined with conventional CF.

Conclusions

In this paper, the HRCF was proposed as a weighting method in linear-array PAI. It was shown that there is a DAS on the numerator of the formula of CF, and it can be replaced with MV beamformer. The proposed method (HRCF) was evaluated numerically and experimentally, and it was shown that it leads to a higher image quality compared to MV and DAS (with/without CF). The quantitative results show that at the depth of 55 mm, compared to DAS+CF and MV+CF, HRCF improves the SNR of about 9 dB and 22 dB, respectively, and reduces the FWHM of about 1752 μm and 44 μm, respectively.

Methods

In this section, the concept of image reconstruction in linear-array PAI, along with the concerned algorithms in this paper, are discussed.

Beamforming In linear-array PAI, a laser illuminates the target of imaging. Then, PA signals are recorded using an US transducer. The detected signals can be used for the image formation using a beamforming algorithm. The most common beamforming algorithm in linear-array PAI is DAS. Its formula is as follows:

$$ y_{DAS}\left(k\right)=\sum_{i=1}^{M}x_{i}\left(k-\Delta_{i}\right), $$

(1)

where y_DAS(k) is the output of the beamformer, k is the time index, M is the number of elements of array, and x_i(k) and Δ_i are the detected signals and the corresponding time delay for the detector i, respectively. To have a more efficient beamformer and improve the reconstructed image, CF can be combined with DAS [40]. The combination of DAS and CF results in sidelobes reduction and contrast enhancement. CF, as an effective weighting process, is given by:

$$ CF(k)=\frac{\left|\sum\limits_{i=1}^{M}x_{i}(k-\Delta_{i})\right|^{2}}{M\sum\limits_{i=1}^{M}\left|x_{i}(k-\Delta_{i})\right|^{2}}. $$

(2)

As can be seen in (2), the argument inside the squared absolute value is the output of DAS algorithm. (1) can be simply implemented and provides a real-time PAI. However, due to the low range of the off-axis signals rejection, it leads to low quality images. The combination of DAS and CF can be written as follows:

$$ y_{DAS+CF}(k)=CF(k)*y_{DAS}(k). $$

(3)

MV can be chosen as an algorithm which provides a high resolution in PAI [47]. However, sidelobes caused by MV highly affect the image quality and degrade the contrast of the reconstructed image. The output of MV adaptive beamformer is given by:

$$ y(k)=\boldsymbol W^{H}(k)\boldsymbol X_{d}(k)=\sum_{i=1}^{M}w_{i}(k)x_{i}(k-\Delta_{i}), $$

(4)

where X_d(k) is the time-delayed array detected signals X_d(k)=[x₁(k), x₂(k),...,x_M(k)]^T, W(k)=[w₁(k),w₂(k),...,w_M(k)]^T is the beamformer weights, and (.)^T and (.)^H represent the transpose and the conjugate transpose, respectively. The detected array signals can be written as follows:

$$ \boldsymbol X(k)=\boldsymbol s(k)+\boldsymbol i(k)+\boldsymbol n(k)= s(k)\boldsymbol a+\boldsymbol i(k)+\boldsymbol n(k), $$

(5)

where s(k),i(k) and n(k) are the desired signal, interference and noise components received by the transducer, respectively. Parameters s(k) and a are the signal waveform and the related steering vector, respectively. MV bemaformer can be used to adaptively weight the calculated samples. Its goal is to achieve the optimal weights in order to estimate the desired signal as accurately as possible. The superiority of MV algorithm has been evaluated in comparison with static windows, such as Hamming window [26]. To acquire the optimal weights, signal-to-interference-plus-noise ratio (SINR) needs to be maximized:

$$ SINR=\frac{\sigma_{s}^{2}|\boldsymbol W^{H}\boldsymbol a|^{2}}{\boldsymbol W^{H}\boldsymbol R_{i+n}\boldsymbol W}, $$

(6)

where R_i+n and $\sigma _{s}^{2}$ are the M×M interference-plus-noise covariance matrix and the signal power, respectively. The maximization of SINR can be gained by minimizing the output interference-plus-noise power while maintaining a distortionless response to the desired signal using following equation:

$$ \min_{\boldsymbol W}\boldsymbol W^{H}\boldsymbol R_{i+n}\boldsymbol W, \quad s.t. \quad \boldsymbol W^{H}\boldsymbol a=1. $$

(7)

The solution of (7) is given by [48]:

$$ \boldsymbol W_{opt}=\frac{\boldsymbol R_{i+n}^{-1}\boldsymbol a}{\boldsymbol a^{H}\boldsymbol R_{i+n}^{-1}\boldsymbol a}. $$

(8)

In practice, the interference-plus-noise covariance matrix is unavailable. Consequently, the sample covariance matrix is used instead of the unavailable covariance matrix using N recently received samples and is given by:

$$ \hat{\boldsymbol R}=\frac{1}{N}\sum_{n=1}^{N}\boldsymbol X_{d}(n)\boldsymbol X_{d}(n)^{H}. $$

(9)

Using MV in medical US imaging encounters some problems which are addressed and discussed in reference [40]. The subarray-averaging or the spatial-smoothing method can be used to achieve a better estimation of the covariance matrix using decorrelation of the coherent signals received by the array. The covariance matrix estimation using the spatial-smoothing can be written as:

$$ \hat{\boldsymbol R}(k)=\frac{1}{M-L+1}\sum_{l=1}^{M-L+1}\boldsymbol X_{d}^{l}(k)\boldsymbol X_{d}^{l}(k)^{H}, $$

(10)

where L is the subarray length, and $\boldsymbol X_{d}^{l}(k)~=~\left [x_{d}^{l}(k),x_{d}^{l+1} (k),...,x_{d}^{l+L-1}(k)\right ]$ is the delayed input signal for the l_th subarray. Due to the limited statistical information, only a few temporal samples are used to estimate the covariance matrix. Therefore, to obtain a stable covariance matrix, the diagonal loading (DL) technique is used. This method leads to replacing $\hat {\boldsymbol R}$ by the loaded sample covariance matrix, $\hat {\boldsymbol R_{l}}=\hat {\boldsymbol R}+\gamma \boldsymbol I$, where γ is the loading factor:

$$ \gamma= \Delta. trace\left\{\hat{\boldsymbol R}(k)\right\}, $$

(11)

where Δ is a constant related to the subarray length. Also, the temporal averaging method can be applied along with the spatial averaging to gain resolution enhancement while the contrast is retained. The estimation of the covariance matrix using both temporal averaging and spatial smoothing in given by:

$$ \begin{aligned} &\hat{\boldsymbol R}(k)=\frac{1}{(2K+1)(M-L+1)}*\\ &\sum_{n=-K}^{K}\sum_{l=1}^{M-L+1}\boldsymbol X_{d}^{l}(k+n)\boldsymbol X_{d}^{l}(k+n)^{H}, \end{aligned}\ $$

(12)

where the temporal averaging is performed over (2K+1) samples. After estimation of the covariance matrix, the optimal weights are calculated by (8) and (12). Finally, the output of MV beamformer is given by:

$$ \hat{y}_{MV}(k)=\frac{1}{M-L+1}\sum_{l=1}^{M-L+1}\boldsymbol W_{*}^{H}(k)\boldsymbol X_{d}^{l}(k). $$

(13)

where W_∗(k)=[w₁(k),w₂(k),...,w_L(k)]^T. Considering (2), it can be seen that the numerator of the fraction is the output of DAS beamformer, and this is why the output of the combination of DAS and CF does not have a high resolution. To put it more simply, the combination of DAS and CF does not provide a high resolution because DAS is weighted using a procedure in which DAS plays a significant role. On the other hand, using MV combined with CF weighting is a good alternative. However, as will be shown in the next section, the output of the combination of DAS and MV can be further improved using HRCF weighting procedure combined with DAS. Its formula is as follows [44, 49]:

$$ y_{DAS+HRCF}(k)=HRCF(k)*y_{DAS}(k), $$

(14)

where

$$ HRCF(k)=\frac{M\left| \hat{y}_{MV}(k)\right|^{2}}{\sum\limits_{i=1}^{M}\left|x_{i}(k-\Delta_{i})\right|^{2}}. $$

(15)

In the next section, the results of the proposed method (the combination of DAS and HRCF) for PA image reconstruction is evaluated.