Open Access
Add to BibSonomy
Abstract
Delay-and-sum (DAS) is the most common algorithm used in photoacoustic (PA) image formation. However, this algorithm results in a reconstructed image with a wide mainlobe and high level of sidelobes. Minimum variance (MV), as an adaptive beamformer, overcomes these limitations and improves the image resolution and contrast. In this paper, a novel algorithm, named Modified-Sparse-MV (MS-MV), is proposed in which a ℓ1-norm constraint is added to the MV minimization problem after some modifications, in order to suppress the sidelobes more efficiently, compared to MV. The added constraint can be interpreted as the sparsity of the output of the MV beamformed signals. Since the final minimization problem is convex, it can be solved efficiently using a simple iterative algorithm. The numerical results show that the proposed method, MS-MV beamformer, improves the signal-to-noise (SNR) about 19.48 dB, in average, compared to MV. Also, the experimental results, using a wire-target phantom, show that MS-MV leads to SNR improvement of about 2.64 dB in comparison with the MV.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Photoacoustic imaging (PAI) is a non-invasive hybrid medical imaging modality that combines the physics of ultrasound (US) and optical imaging modalities [1,2]. It uses short pulses to illuminate the imaging medium and consequently, generate the acoustic waves. The major benefit of PAI, in comparison with optical imaging, is its higher penetration depth while a high contrast is achieved. It has several applications such as functional and structural imaging [3], ocular imaging [4], tumor detection [5] and imaging the whole body of small animal [6]. PAI is divided into two categories: Photoacoustic microscopy (PAM) and Photoacoustic tomography (PAT) which is suitable for in vivo applications [7]. In PAT, an array sensor which is formed in linear, circular or arc shape, is used to detect the propagated photoacoustic (PA) waves from the surface of the medium [8]. Therefore, a reconstruction algorithm is necessary to construct the absorption distribution of the medium [9]. Due to the high similarity between the US and PAI, beamforming algorithms used in US can be applied to the PAI with some modifications [10]. Delay-and-Sum (DAS) is the most common algorithm used to reconstruct the images in US and PAI due to its low computational complexity. However, the reconstructed images obtained from this algorithm suffers from wide mainlobe and high level of sidelobes. Delay-Multiply-and-Sum (DMAS) is proposed to enhance the image quality compared to DAS [11]. In order to improve the reconstructed image in the term of contrast, in comparison with DMAS, some other algorithms have been proposed [12,13]. Adaptive beamformers, which have a wide range of applications in Radar and Sonar, weight the received PA waves proportional to the characteristics of the signals. Therefore, adaptive beamformers would overcome the limitations exist in non-adaptive beamformers and improve the image quality. One of the adaptive beamformers that has been improved in medical image reconstruction, is minimum variance (MV) algorithm [14], in which the weights are calculated in a way that the interference-plus noise power of the output is minimized while the unit gain at the focal point is retained [15]. The reconstructed images obtained from this algorithm is significantly improved in terms of mainlobe width and sidelobe levels compared to non-adaptive beamformers. Some modifications have been applied to the MV beamformer to improve the image quality. Coherence factor (CF) has been applied to MV to improve the resolution and contrast [16,17]. A modification has been applied to CF in order to suppress the sidelobes and improve the resolution more efficiently compared to the conventional CF [18]. A novel beamforming algorithm, called MV-Based DMAS (MVB-DMAS) is introduced in which the MV beamformer is applied inside the expansion of the non-adaptive DMAS beamformer [19,20]. Temporal averaging has been employed to improve the contrast while the resolution is retained [21]. Also, Eigenspace-Based Minimum Variance (EIBMV) and its combination with DMAS have been employed in order to improve the image quality [22–24].
Recently, sparse-based algorithms have been used for image formation, especially in tomographic image reconstruction where the difference between the measured signals and the predicted sparse beamformed data is considered to be minimized through a least square model, referred as the forward problem. In addition, some regularization terms are combined with the mentioned minimization problem in order to suppress the artifacts caused by the limited-view and improve the image quality [25–28]. The effect of the regularization term is investigated in different domains, such as frequency domain and wavelet domain, using the sparsifying basis [29, 30]. However, the reconstructed images obtained by these algorithms still suffer from the wide mainlobe compared to MV-based algorithms. It should be noted that the sparse-based algorithms are expensive in terms of memory requirement.
Dispersion Distortionless Response (MDDR) is another beamforming method in which the dispersion term is interpreted as a generalized form of the variance. This algorithm is suitable for non-Gaussian signals and noise, where the necessary information is obtained from higher-order and lower-order statistics of the output. In MDDR minimization problem, the weight is estimated in a way that the expectation of the ℓp-norm (p ≥ 1) of the output is minimized, while the unit gain at the focal point is retained unity [31]. It is obvious that MDDR would be reduced to MV beamformer for p = 2, where the second-order statistics of the output is exploited. Higher order statistics can be exploited for p ≥ 2 and lower order statistics can be exploited for p ≤ 2. Generally, the value of p is chosen depending on the characteristics of the noise added to the signals [32].
In order to suppress the sidelobes more efficiently, several beamforming methods have been developed in Radar and Sonar applications by adding a new constraint to the existing optimization problem; In [33], a sparse beamforming method, named ℓ1-regularized minimum absolute distortionless response (ℓ1-MADR), is introduced in which ℓ1-regularization term is added to MDDR beamformer with the assumption of p = 1. It has been shown that this robust algorithm suppresses the sidelobe level more efficiently compared to MV. In [34], a new constraint is added to the MV criterion in order to maximize the Mainlobe-to-Sidelobe Power Ratio (MSPR). Sparse capon (SC) beamforming is another method in which a ℓp-regularization term (p ≤ 1) is added to the MV beamforming constraint in order to force the sparsity to the beampattern and improve the reconstructed images in terms of sidelobes compared to the pure MV beamformer [35,36].
In this paper, it is proposed to add a new ℓ1-norm sparse constraint to the existing MV beamformer after some modifications in order to enhance the image quality and further suppress the sidelobes. The new added constraint can be interpreted as the sparsity of the output which is forced to the entire beampattern. The optimum weight would be obtained from a simple iterative algorithm used in [33]. The results show that the proposed algorithm, named modified-sparse-MV (MS-MV), leads to a significant image improvement in the term of sidelobe levels compared to MV.
The rest of the paper is organized as follows. In section 2, MV beamforming method is described. The proposed method is presented in section 3. The results obtained from the numerical data are presented in section 4. Also, the experimental results are evaluated in section 5. Finally, the discussion and the conclusion are reported in section 6 and 7, respectively.
2. Background
2.1. Minimum variance beamformer
In MV beamforming method, the output of the beamformed signal, y(k), is defined as below:
where M is the number of elements in the linear array sensor, xm(k) is the generated PA signal received by mth element, k is the time index, Δm is the time delay which is calculated proportional to the distance between the point target and mth element and applied to the received signals, and wm(k) is the calculated weight corresponding to the delayed received signals. Equation (1) can be written as below: where W(k) = [w1(k), w2(k), · · · , wm(k)]T is the array of weights and X(k) = [x1 (k − Δ1(k)), x2 (k − Δ2(k)), · · · , xM (k − ΔM(k))]T is the array of delayed received signals. The weight is calculated in a way that the interference-plus-noise power of the output is minimized, while the unit gain is retained at the focal point [37]: where R(k) is the spatial covariance matrix of the interference-plus-noise, a is the steering vector and (.)H denotes the conjugate transpose operator. It should be noted that a is considered as a vector of ones since the received signals are delayed. The optimum weight would be obtained using Lagrange multiplier method: The optimum weight is obtained in a way that the large sidelobes are allowed in directions where the received energy is low. Therefore, the off-axis signals would be suppressed [15]. Note that due to the non-stationary property of the signals in medical US and PA, obtaining the exact spatial covariance matrix is not possible. Therefore, the spatial covariance matrix should be estimated in the way that the array sensor is divided into M − L + 1 overlapping subarrays with length L. Then, R(k) is estimated by averaging the calculated covariance matrixes of each subarray, referred to as spatial smoothing technique. In addition to the spatial smoothing, temporal averaging over 2K + 1 samples also can be used in spatial covariance matrix estimation to achieve a better resolution, while the contrast is retained [21]. The estimated covariance matrix is written as below: where X̄l(k) = [xl(k), xl+1(k), · · · , xl+L−1(k)]T is the delayed received signals of lth subarray. Note that the smaller L increases the number of subarrays, and consequently a more robust estimation would be achieved. However, the resolution would be decreased due to the reduction of employed signals in each subarray. Therefore, it can be concluded that there is a trade off between the subarray length and the resolution. Also, the diagonal loading (DL) can be applied to the spatial covariance matrix to obtain a more robust estimation. To apply DL, Δ.trace{R} should be added to the diagonal of the estimated covariance matrix before the weight is calculated. Δ is a constant parameter that corresponds to the subarray length and is much smaller than [38]. The output beamformed signal would be obtained after the weight is estimated using R̂(k) instead of (k) in (4) from the following equation:3. Proposed method
3.1. Sparse MV beamforming
In SC beamforming method, which has a wide range of applications in the field of Radar and Sonar [35], a sparse constraint is applied to the MV minimization problem in order to force the sparsity to the hole beampattern and suppress the sidelobes more efficiently compared to the MV beamformer. The optimization problem is defined as follows:
where α is the regularization parameter, C is a matrix which covers K steering vectors sampled from −90° to 90°. Note that in PAI, the steering vector is considered as a vector of ones in all directions, as mentioned before; therefore, the C would be only a matrix of ones with a dimension of L × K. ‖.‖p represents the ℓp-norm of a vector. Considering p = 1, this minimization problem is convex, and therefore, it can be solved efficiently using CVX MATLAB toolbox [39]. In this paper, it is suggested to add a new ℓ1-norm sparse constraint to this minimization problem and generate sparse-MV (S-MV) beamformer in order to enhance the image quality and suppress the sidelobes more efficiently. The new optimization problem is defined as follows: where XT = [X̄1, X̄2, · · · , X̄M+L−1]. The parameter β is considered as a weighting factor which determines the balance between the MV constraint and the new added sparse constraint. In the other words, the sparsity of the output of the MV beamformed signal is forced to the beampattern by adding the new ℓ1-norm constraint. This new added constraint leads to significantly image improvement in terms of sidelobe suppression. Combining the two ℓ1-norm constraints, (8) can be expressed as follows: where B = [βXT, αC]H. In this work, the optimum weight is obtained using the iterative algorithm used in [33]. In the following, the problem solving procedure is explained.3.1.1. Problem solving
Equation (9) can be rewritten as below:
where Φ1 is a diagonal weighting matrix: The motivation behind this scheme, is to convert the ℓp-norm form into the ℓ2-norm form, and therefore make the optimization problem solvable [31]. From this new form of the minimization problem, the optimum weight would be obtained using the Lagrange multiplier method; First, the objective function f(WS, γ) is created as bellow: where By taking the derivation of the objective function with respect to the WS and γ, separately, we have: Multiplying two sides of (14) to (aH),we have: and finally, the optimum weight would be obtained using (16) and (14): where (.)k indicates kth(k = 0, 1, · · ·) step of the iteration procedure. Since this iterative algorithm is not sensitive to the initial weight, the MV beamformer or the Bartlett beamformer can be used to obtain the initial weight. The iteration procedure continues until the calculated weights are converged, according to the following criteria: The results show that the optimum weight obtained from (17) leads to image enhancement in terms of sidelobe levels and noise suppression compared to MV.3.2. MS-MV beamforming method
The proposed S-MV minimization problem, consists of two sparse constraints; the first sparse constraint which comes from the SC minimization problem, and the new added constraint which forces the sparsity of the output to the entire beampattern. Analysing the first constraint, ‖CHWS‖1, we have concluded that this constraint does not have any effect on the sparsity of the beampattern; the steering vector, which is considered as a vector of ones in PA and US, forces the added constraint to be a constant parameter. In other words, a constant parameter is added to the spatial covariance matrix for all time indexes using this constraint, which can not make any changes to the calculated weight compared to MV (see Appendix). In fact, the improvement of the reconstructed images obtained from S-MV, shown in (8), is mainly due to the second new added constraint. As a result, the first constraint could be omitted while the quality of the reconstructed image remains unchanged. The proposed minimization problem, called modified-sparse-MV (MS-MV), finally would be achieved as follows:
The optimum weight (WMS) would be achieved similar to the previous minimization problem and the final proposed iterative algorithm would be expressed as below: where , and The results in the next section show that the proposed MS-MV beamforming method would suppress the sidelobes more efficiently compared to MV beamformer.4. Numerical result
4.1. Imaging setup
The simulation is performed on a PC with 3.3 GHz core(TM) i7-5820k CPU and 64 GB memory. The k-wave MATLAB toolbox is used to simulate the array sensor and the absorbers [40]. An imaging region is designed with the vertical and lateral axis of 70 mm and 20 mm, respectively. Ten 0.1 mm radius spherical absorbers, as initial pressures, are centered on the lateral axis. The absorbers are located along the vertical axis. 5 mm distance is considered between each two absorbers, starting from 20 mm of the array sensor. A linear array sensor including 128 elements is used to detect the propagated PA signals with the central frequency of 5 MHZ and 77% bandwidth. The speed of sound is assumed to be 1540 m/s. In all simulation cases, the length of the subarrays is considered L = M/2. DL is applied to the spatial covariance matrix with the assumption of δ = 1/100L. Also, temporal averaging over 5 samples (K = 2) is applied. Gaussian noise is added to the received signals. Finally, Hilbert transform, normalization and log-compression procedure are performed after applying the reconstruction algorithm.
4.2. Effect of varying α and β
Before the qualitative evaluation of the reconstructed images using different algorithms, the effect of varying α and β needs to be discussed in order to fix them in the procedure of the simulations. Note that α is a parameter that determines the balance between the first added constraint ‖CHW‖1 and MV constraint, as mentioned before. It can be concluded that different values of α would not make any changes to the reconstructed images due to the effectiveness of this added constraint to S-MV minimization problem. Therefore, only the different values of β needs to be investigated. Consider the reconstructed images of a single point target located at the depth of 250 mm using MS-MV for three different values of β (≥ 0). From the lateral variations shown in Fig. 1, it can be seen that the lowest level of sidelobes would be obtained using β = 1. Note that the smaller value of β causes the algorithm to behave like MV. In this work, the simulation is performed with the consideration of β = 1. Due to the effectiveness of the first added constraint in S-MV, and therefore, the equality of S-MV and MS-MV, the reconstruction procedure is performed using MS-MV in the next part of this section.
4.3. Qualitative evaluation
Figure 2 shows the reconstructed images of the simulated point targets. The reconstructed images using DAS, MV and DS-MV are depicted in Fig. 2(a), Fig. 2(b) and Fig. 2(c), respectively. Considering the images, it is obvious that the non-adaptive DAS beamformer results in the reconstructed image with the lowest quality. MV enhances the image resolution and contrast compared to DAS. Comparing Fig. 2(b) and Fig. 2(c), it can be seen that the proposed DS-MV algorithm improves the image contrast compared to MV, as expected. To have a better evaluation, the lateral variations of the reconstructed images are presented in Fig. 3, at four different depths of imaging. It is demonstrated that the proposed DS-MV leads to image improvement in terms of sidelobe suppression more effectively compared to MV and the non-adaptive DAS beamformers, at all the depths of imaging.
Fig. 2 The numerical reconstructed PA images of 10 point targets using (a) DAS, (b) MV, (c) DS-MV. Noise is added to the received signals having a SNR of 50 dB.
Fig. 3 The lateral variations of the images shown in Fig. 2 at the depth of (a) 30 mm, (b) 40 mm, (c) 55 mm and (d) 65 mm. Noise is added to the received signals having a SNR of 50 dB.
4.4. Quantitative evaluation
One of the most common metrics to evaluate the reconstructed images quantitatively, is signal-to-noise ratio (SNR), which is calculated as below:
where Psignal is the difference between the maximum and the minimum intensity of the reconstructed image, and Pnoise is the standard deviation of the reconstructed image [12, 18]. The SNRs of the reconstructed images shown in Fig. 2 are calculated and presented in Table 1. As shown, the DAS leads to the lowest value of SNR due to the existence of high level of sidelobes, as expected. Comparing the calculated SNRs corresponding to MV and MS-MV, it can be concluded that the MS-MV beamformer significantly improves the SNR due to the added sparse constraint to the MV criterion. Full-width-half-maximum (FWHM) is another metrics to evaluate the performance of the beamformers. This metric is used in order to estimate the spatial resolution of the reconstructed images. The calculated FWHM is presented in Table 2 in all depths of positioning point targets. As can be seen, the highest value of calculated FWHM belongs to DAS. Therefore, it can be concluded that DAS leads to a reconstructed image with the lowest resolution, compared to other beamformers, as expected. Also, It can be concluded that no significant resolution improvement occurred using the proposed MS-MV in comparison with MV since the calculated FWHMs are almost similar. It worth to mention that the calculated FWHMs are almost constant in all depths of positioning points using MV and MS-MV, as shown in Table 2; the difference between the calculated FWHMs at the first and the last depth of positioning points is 2.8 mm, 0.04 mm and 0.01 mm using DAS, MV and MS-MV, respectively. As a result, it can be concluded that the resolution remains almost unchanged using MV and MS-MV in all depths of imaging compared to the non-adaptive DAS beamformer.
Table 1. The calculated SNR (dB) in different depths.
Table 2. The calculated FWHM (mm) in different depths.
4.5. Imaging at the presence of high level of noise
To evaluate the performance of the beamformers at the presence of high level of noise, a Gaussian noise is added to the received signals, resulting the SNR of 10 dB. The reconstructed images are shown in Fig. 4. The reconstructed images using DAS, MV and DS-MV are depicted in Fig. 4(a), Fig. 4(b) and Fig. 4(c), respectively. Also, the lateral variations in two different depths of positioning points are shown in Fig. 5. The results show that the reconstructed image using DAS is more affected by the presence of noise. Also, it can be concluded that the reconstructed image using DS-MV is significantly improved in terms of noise reduction compared to other beamformers.
Fig. 4 The reconstructed PA images of 10 point targets using (a) DAS, (b) MV, (c) DS-MV. Noise is added to the received signals having a SNR of 10 dB.
Fig. 5 The lateral variations of the images shown in Fig. 4 at the depth of (a) 40 mm and (b) 50 mm . Noise is added to the received signals having a SNR of 10 dB.
5. Experimental results
5.1. Wire-target phantom
To further evaluate the introduced algorithm and its effect on PA image improvement, phantom experiments were performed in which a phantom consists of 1 and 2 light absorbing wires with diameter of 150 μm were placed 1 mm apart from each other in a water tank. The schematic of the experimental setup is shown in Fig. 6. In the wire-target experiments, we utilized a Nd:YAG pulsed laser, with the pulse repetition rate of 30 Hz at wavelengths of 532 nm. A programmable digital ultrasound scanner (Verasonics Vantage 128), equipped with a linear array transducer (L11-4v) 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.
5.1.1. Qualitative and quantitative evaluation
First, a single wire target was used as the phantom of imaging. The reconstructed images using the concerned beamformers are shown in Fig. 7 where the higher resolution of the MV-based algorithms are clear, compared to the DAS. As demonstrated in Fig. 7(c), a reduced background noise is obtained, in comparison with the image shown in Fig. 7(b), indicating the higher SNR of the MS-MV. To quantitatively evaluate the experimental results, SNR was used where DAS, MV and MS-MV leads to 39.2 dB, 40.5 dB, 43 dB, respectively, which shows the superiority of the proposed method. A more complex target (using two wires as the target) was utilized, and the results are presented in Fig. 8. As demonstrated, DAS results in a low resolution image along with high sidelobes. MV improves the resolution, but the sidelobes still affect the reconstructed image. On the other hand, MS-MV leads to a high resolution while the sidelobes are degraded, and the presence of noise is suppressed in the image, as can be seen in Fig. 8(c). To evaluate the images in detail, lateral variations are presented in Fig. 9 where MS-MV leads to a better resolution compared to DAS while the sidelobes are degraded compared to MV (see the arrows and circles). It can be seen that at the depth of 22 mm, using MS-MV, the sidelobes are degraded about 8 dB, compared to MV. CF, as a weighting factor, can be added to a beamforming algorithm to eliminate the limited-view artifacts and improve the image quality, as mentioned before. It is worth to apply CF to MV (MV+CF) and the proposed MS-MV (MS-MV+CF) algorithms in order to evaluate and compare its performance on these two beamformers. The beamforming results are depicted in Fig. 10. The lateral variations are presented in Fig. 11. Considering Figs. 10 and 11, it can be seen that MS-MV+CF leads to improvement in image quality compared to MV+CF. The SNR is improved by almost 2.85 dB and 6.01 dB at the depth of 22 mm and 24 mm, respectively. The FWHM metric to compare various beamforming algorithms are presented in TABLE 3. It can be concluded that the largest value of FWHM belongs to DAS, and therefore, the resolution obtained by DAS is not well enough, as expected. MS-MV leads to an enhanced image in term of resolution, compared to MV. Comparing MV+CF and MS-MV+CF, it can be seen that MS-MV+CF improves the FWHM value in both depths of imaging wire targets.
Fig. 7 Reconstructed experimental PA images using (a) DAS, (b) MV and (c) MS-MV (β = 1). All the images are shown with a dynamic range of 50 dB. A wire target was used as the imaging target.
Fig. 8 Reconstructed experimental PA images using (a) DAS, (b) MV and (c) MS-MV (β = 1). All the images are shown with a dynamic range of 50 dB. Two wires were used as the imaging target.
Fig. 9 The lateral variations for the reconstructed experimental PA images shown in Fig. 8. Arrows and circles demonstrate the improvement caused by MS-MV algorithm.
Fig. 10 Reconstructed experimental PA images using (a) MV, (b) MV+CF, (c) MS-MV and (d) MS-MV+CF.Two wires were used as the imaging target.
Table 3. The calculated FWHM (mm) for wire phantom shown in Fig. 10.
5.2. Ex vivo imaging
In this study, an ex vivo experimental tissue study have been designed to evaluate the performance of the proposed algorithm. A piece of a chicken breast tissue (about 4 cm × 4 cm × 3 cm) is extracted. Two pencil leads with a diameter of 0.5 mm are embedded inside the breast tissue, having an axial distance of about 5 mm. Fig 12 shows the photographs of the imaged tissue. The PA signals are collected with a combined linear US/PA imaging probe. The reconstructed images are presented in Fig. 13. As demonstrated, DAS leads to a low resolution while MV provides a higher resolution. On the other hand, the background of the ex vivo image obtained by MS-MV is darker which indicates a reduced background noise, in comparison with MV, which indicates the higher SNR of the proposed method. It should be noticed that the sidelobes are lower using MS-MV while the resolution of MV is retained. The lateral variations shown in Fig. 14 clearly proves the lower sidelobes of MS-MV (15 dB), in comparison with MV. Moreover, the calculated FWHM of MV and the proposed MS-MV method at this depth is 0.69 mm and 0.34 mm, respectively. This indicates the improvement of image quality in term of resolution using MS-MV algorithm.
Fig. 13 Reconstructed ex vivo images using (a) DAS, (b) MV and (c) MS-MV. A linear-array and the phantom shown in Fig. 12 were used for the experimental design. All the images are shown with a dynamic range of 50 dB.
6. Discussion
The major benefits of the proposed DS-MV algorithm, is the improvement of the reconstructed images in terms of sidelobe levels and noise reduction compared to DAS and MV beamformers. DAS is a non-adaptive algorithm in which the weights are predefined, and therefore, it can not reject the of-axis signals contribution as well as adaptive beamformers. As a results, the reconstructed images using DAS leads to a low quality reconstructed image in terms of resolution and contrast. In the adaptive MV beamformer, the calculated weights are changed proportional to the characteristics of the received signals. Therefore, the reconstructed images would be improved in terms of mainlobe width and sidelobe levels. By adding the new ℓ1-norm constraint to the MV criterion, the sidelobes are suppressed more efficiently compared to MV. Consider Fig. 2 and Fig. 3, for instance, where the reconstructed images of the wire phantom and the corresponding lateral variations are shown, respectively. It can be seen that the sparse added constraint leads to a better image quality and more noise reduction compared to MV. Furthermore, the proposed algorithm results in the reconstructed image with a better quality at the presence of high level of noise, as shown in Fig. 4, where 10 dB noise is added to the received signals. In fact, the proposed constraint added to the minimization problem in DS-MV, ‖XHW‖1, forces the sparsity of the output to the whole beampattern; it takes the smaller values to the signals with large amplitudes. It should be noted that most of the large amplitudes exist in mainlobe, but the mainlobe is not affected by this added sparse constraint since the sparse constraint is combined with the other constraint, s.t. aHW = 1. It means that the sparse constraint is performed on the beampattern, while the focal point remains unity. It is worth to mention that an accurate metric of the sparsity is ℓ0-norm, where the non-zero elements of the vector is counted. However, ℓ0-norm constraint leads to a non-convex optimization problem. Therefore, ℓ1-norm is used to make the minimization problem convex and solvable. Note that in US and PA, the steering vector is considered as a vector of ones in all direction, as mentioned before. That’s why the constraint ‖CHW‖1 does not make any changes to the reconstructed images compared to MV (see APPENDIX for more information). In contrast to US and PA, in Radar and Sonar applications the steering vector varies in each direction, and therefore, this will cause the mentioned constraint to be considered as an efficient constraint which results in sidelobe suppression more efficiently compared to MV. In the other words, in US and PA, there is not any difference between the proposed S-MV and MS-MV due to the effectiveness of the first added constraint, ‖CHW‖1, exists in S-MV beamformer. However, these two algorithms would behave differently in Radar and Sonar applications. Finally, it should be noted that the computational complexity of DAS, MV and the proposed DS-MV is O(M) and O(M3) and O(Niter M3), respectively, where Niter is the number of iteration which needs to satisfy the convergence criterion expressed in (18). Obviously, DAS algorithm can reconstruct the images faster than other beamformers. Comparing MV and DS-MV, it can be concluded that the computational complexity of the proposed DS-MV is comparable with MV.
7. Conclusion
One of the most common algorithms used in PAI, is the non-adaptive DAS beamformer due to its simple implementation. However, the quality of the reconstructed images obtained from this algorithm suffers from the wide mainlobe and high sidelobes. Adaptive beamformers, such as MV, overcomes these limitation by calculating the weights based on the characteristics of the signals and results in an improved reconstructed images in terms of mainlobe width and sidelobe levels. In this paper, it was proposed to add a new constraint to the MV minimization criterion in order to further suppress the sidelobes. The new ℓ1-norm added constraint can be interpreted as the sparsity of the output that is forced to the beampattern. The numerical and experimental results showed that the proposed DS-MV algorithm outperforms MV in terms of sidelobe suppression and noise reduction. The quantitative results obtained from DS-MV showed that the calculated SNR was improved about 19.48 dB and 2.64 dB for the simulated point targets and the designed wire phantom, respectively, compared to MV.
Appendix
7.1. Analysis of the sparse constraint (‖CHW‖1)
In US and PAI applications, the first sparse constraint in S-MV minimization problem, ‖CHW‖1, does not have any effect on the reconstructed images compared to MV beamformer; the steering vector is considered as a vector of ones in all directions, as mentioned before. This property makes the added constraint to be a constant parameter which is added to the covariance matrix. It can be seen that adding a constant parameter would not make any changes to the estimated weight compared to the estimated weight obtained from MV. For better understanding, consider SC minimization problem in (7), in which the desired constraint is added to MV minimization problem. The optimum weight is obtained from the same method mentioned before:
where D(W) = diag {|CHW(1)|−1, · · · , |CHW(N)|−1}. D(W) can be written as below: It should be noted that each elements in diagonal of the matrix in (24), equals to the second constraint combined with the first constraint in SC beamformer, WHa = 1. Therefore, D(W) would be a diagonal matrix of ones with dimensions K × K. By expanding the term added to the covariance matrix shown in (23), we have: It is obvious from (25) that the term added to the covariance matrix would be a N × N matrix where each elements of this matrix is a constant parameter, αN. As a result, (23) would be rewritten as below: in which the constant parameter, αN, is added to the covariance matrix. In order to show the equality of the calculated weights obtained from (4) and (26), consider the covariance matrix with dimension 2 × 2 (N = 2) for simplicity, as below: The inverse of (27) would be as follows: Referring to (4), the estimated weight obtained from MV is as follows: Adding the constant term to the covariance matrix, we have: The inverse of (30) would be as below: Referring to (26), the optimum weight obtained from SC beamformer is estimated as follows: Comparing (29) and (31), it can be seen that the weigh obtained from MV equals to the weight obtained from SC beamformer in which a sparse constraint is added to the existing minimization problem (wMV = wSC). In the other words, it can be concluded that the mentioned sparse constraint does not make any changes to the calculated weight compared to MV in US and PA image reconstruction.Acknowledgments
This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.
Disclosures
The authors declare that there are no conflicts of interest related to this article.
References and links
1. B. Wang, L. Xiang, M. S. Jiang, J. Yang, Q. Zhang, P. R. Carney, and H. Jiang, “Photoacoustic tomography system for noninvasive real-time three-dimensional imaging of epilepsy,” Biomed. Opt. Express 3, 1427–1432 (2012). [CrossRef] [PubMed]
2. M. Mehrmohammadi, S. Joon Yoon, D. Yeager, and S. Y Emelianov, “Photoacoustic imaging for cancer detection and staging,” Curr. Mol. Imaging 2, 89–105 (2013). [CrossRef] [PubMed]
3. J. Yao and L. V. Wang, “Sensitivity of photoacoustic microscopy,” Photoacoustics 2, 87–101 (2014). [CrossRef] [PubMed]
4. S. Y. Nam and S. Y. Emelianov, “Array-based real-time ultrasound and photoacoustic ocular imaging,” J. Opt. Soc. Korea 18, 151–155 (2014). [CrossRef]
5. M. Heijblom, W. Steenbergen, and S. Manohar, “Clinical photoacoustic breast imaging: the twente experience,” IEEE Pulse 6, 42–46 (2015). [CrossRef] [PubMed]
6. J. Xia and L. V. Wang, “Small-animal whole-body photoacoustic tomography: a review,” IEEE Trans. Biomed. Eng. 61, 1380–1389 (2014). [CrossRef]
7. L. V. Wang, “Multiscale photoacoustic microscopy and computed tomography,” Nat. Photon. 3, 503–509 (2009). [CrossRef]
8. M. Nasiriavanaki, J. Xia, H. Wan, A. Q. Bauer, J. P. Culver, and L. V. Wang, “High-resolution photoacoustic tomography of resting-state functional connectivity in the mouse brain,” Proc. Natl. Acad. Sci. 111, 21–26 (2014). [CrossRef]
9. M. Xu, Y. Xu, and L. V. Wang, “Time-domain reconstruction algorithms and numerical simulations for thermoacoustic tomography in various geometries,” IEEE Trans. Biomed. Eng. 50, 1086–1099 (2003). [CrossRef] [PubMed]
10. C.-W. Wei, T.-M. Nguyen, J. Xia, B. Arnal, E. Y. Wong, I. M. Pelivanov, and M. O’Donnell, “Real-time integrated photoacoustic and ultrasound (paus) imaging system to guide interventional procedures: ex vivo study,” IEEE Trans. Ultrason., Ferroelectr., Freq. Control. 62, 319–328 (2015). [CrossRef]
11. 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, 940–949 (2015). [CrossRef]
12. 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, 31–42 (2018). [CrossRef]
13. M. Mozaffarzadeh, A. Mahloojifar, and M. Orooji, “Image enhancement and noise reduction using modified delay-multiply-and-sum beamformer: Application to medical photoacoustic imaging,” in “Iranian Conference on Electrical Engineering (ICEE), 2017,” (IEEE, 2017), pp. 65–69.
14. J. Capon, “High-resolution frequency-wavenumber spectrum analysis,” Proc. IEEE 57, 1408–1418 (1969). [CrossRef]
15. J.-F. Synnevag, A. Austeng, and S. Holm, “Benefits of minimum-variance beamforming in medical ultrasound imaging,” IEEE Trans. Ultrason., Ferroelectr., Freq. Control. 56, 1263 (2009). [CrossRef]
16. S.-L. Wang and P.-C. Li, “MVDR-based coherence weighting for high-frame-rate adaptive imaging,” IEEE Trans. Ultrason., Ferroelectr., Freq. Control. 56, 10943182 (2009).
17. S. Park, A. B. Karpiouk, S. R. Aglyamov, and S. Y. Emelianov, “Adaptive beamforming for photoacoustic imaging using linear array transducer,” in “Ultrasonics Symposium, 2008. IUS 2008. IEEE,” (IEEE, 2008), pp. 1088–1091.
18. M. Mozaffarzadeh, Y. Yan, M. Mehrmohammadi, and B. Makkiabadi, “Enhanced linear-array photoacoustic beamforming using modified coherence factor,” J. Biomed. Opt. 23, 026005 (2018). [CrossRef]
19. M. Mozaffarzadeh, A. Mahloojifar, and M. Orooji, “Medical photoacoustic beamforming using minimum variance-based delay multiply and sum,” in “Digital Optical Technologies 2017,”, vol. 10335 (International Society for Optics and Photonics, 2017), vol. 10335, p. 1033522.
20. M. Mozaffarzadeh, A. Mahloojifar, M. Orooji, K. Kratkiewicz, S. Adabi, and M. Nasiriavanaki, “Linear-array photoacoustic imaging using minimum variance-based delay multiply and sum adaptive beamforming algorithm,” J. Biomed. Opt. 23, 026002 (2018). [CrossRef]
21. J.-F. Synnevåg, C.-I. Nilsen, and S. Holm, “P2b-13 speckle statistics in adaptive beamforming,” in “Ultrasonics Symposium, 2007. IEEE,” (IEEE, 2007), pp. 1545–1548.
22. B. M. Asl and A. Mahloojifar, “Contrast enhancement of adaptive ultrasound imaging using eigenspace-based minimum variance beamfoming,” in “Ultrasonics Symposium (IUS), 2009 IEEE International,” (IEEE, 2009), pp. 349–352.
23. S. Mehdizadeh, A. Austeng, T. F. Johansen, and S. Holm, “Eigenspace based minimum variance beamforming applied to ultrasound imaging of acoustically hard tissues,” IEEE Trans. Med. Imag. 31, 1912–1921 (2012). [CrossRef]
24. M. Mozaffarzadeh, A. Mahloojifar, M. Nasiriavanaki, and M. Orooji, “Eigenspace-based minimum variance adaptive beamformer combined with delay multiply and sum: experimental study,” in “Photonics in Dermatology and Plastic Surgery 2018,”, vol. 10467 (International Society for Optics and Photonics, 2018), vol. 10467, p. 1046717.
25. R. Jirik, I. Peterlik, N. Ruiter, J. Fousek, R. Dapp, M. Zapf, and J. Jan, “Sound-speed image reconstruction in sparse-aperture 3-d ultrasound transmission tomography,” IEEE Trans. Ultrason., Ferroelectr., Freq. Control. 59, 2185 (2012). [CrossRef]
26. Y. Zhang, Y. Wang, and C. Zhang, “Total variation based gradient descent algorithm for sparse-view photoacoustic image reconstruction,” Ultrasonics. 52, 1046–1055 (2012). [CrossRef] [PubMed]
27. R. Lavarello, F. Kamalabadi, and W. D. O’Brien, “A regularized inverse approach to ultrasonic pulse-echo imaging,” IEEE Trans. Med. Imag. 25, 712–722 (2006). [CrossRef]
28. E. Ozkan, V. Vishnevsky, and O. Goksel, “Inverse problem of ultrasound beamforming with sparsity constraints and regularization,” IEEE Trans. Ultrason., Ferroelectr., Freq. Control. 65, 356–365 (2018). [CrossRef]
29. Y. Zhang, Y. Wang, and C. Zhang, “Efficient discrete cosine transform model–based algorithm for photoacoustic image reconstruction,” J. Biomed. Opt. 18, 066008 (2013). [CrossRef]
30. A. Rosenthal, T. Jetzfellner, D. Razansky, and V. Ntziachristos, “Efficient framework for model-based tomographic image reconstruction using wavelet packets,” IEEE Trans. Med. Imag. 31, 1346–1357 (2012). [CrossRef]
31. X. Jiang, W.-J. Zeng, A. Yasotharan, H.-C. So, and T. Kirubarajan, “Minimum dispersion beamforming for non-gaussian signals,” IEEE Trans. Signal Process. 62, 1879–1893 (2014). [CrossRef]
32. P. Tsakalides and C. L. Nikias, “Robust adaptive beamforming in alpha-stable noise environments,” 1996 IEEE International Conference on Acoustics, Speech, and Signal Processing, vol. 5 (IEEE, 1996), vol. 5, pp. 2884–2887.
33. X. Jiang, A. Yasotharan, and T. Kirubarajan, “Robust beamforming with sidelobe suppression for impulsive signals,” IEEE Signal Process. Lett. 22, 346–350 (2015). [CrossRef]
34. Y. Liu and Q. Wan, “Sidelobe suppression for robust capon beamforming with mainlobe-to-sidelobe power ratio maximization,” IEEE Antennas Wirel. Propag. Lett. 11, 1218–1221 (2012). [CrossRef]
35. Y. Liu, “Robust capon beamforming via shaping beam pattern,” arXiv preprint arXiv:1302.6173 (2013).
36. Y. Zhang, B. Ng, and Q. Wan, “Sidelobe suppression for adaptive beamforming with sparse constraint on beam pattern,” Electron. Lett. 44, 615–616 (2008). [CrossRef]
37. R. A. Monzingo and T. W. Miller, Introduction to adaptive arrays (Scitech publishing, 1980).
38. J. F. Synnevag, A. Austeng, and S. Holm, “Adaptive beamforming applied to medical ultrasound imaging,” IEEE Trans. Ultrason., Ferroelectr., Freq. Control. 54, 431 (2007). [CrossRef]
39. J. F. Sturm, “Using sedumi 1.02, a matlab toolbox for optimization over symmetric cones,” Optim. Methods Softw. 11, 625–653 (1999). [CrossRef]
40. B. E. Treeby and B. T. Cox, “k-wave: Matlab toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. 15, 021314 (2010). [CrossRef] [PubMed]
