## 1.

## Introduction

Photoacoustic imaging (PAI) is a promising medical imaging modality that uses a short electromagnetic pulse to generate ultrasound (US) waves based on the thermoelastic effect.^{1} Having the merits of the US imaging spatial resolution and the optical imaging contrast in one imaging modality is the main motivation of using PAI.^{2} Unlike the x-ray that uses an ionizing radiation, PAI uses nonionizing waves, i.e., short laser or radio frequency pulses. In comparison with other imaging modalities, PAI has multiple advantages leading to many investigations.^{3}^{,}^{4} PAI is a multiscale imaging modality that has been used in different cases of study such as tumor detection,^{5}^{,}^{6} cancer detection and staging,^{7} ocular imaging,^{8} monitoring oxygenation in blood vessels,^{9} and functional imaging.^{10}^{,}^{11} There are two techniques of PAI: photoacoustic tomography (PAT) and photoacoustic microscopy.^{12}^{,}^{13} In 2002, for the first time, PAT was successfully used as *in vivo* functional and structural brain imaging modality in small animals.^{14} 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 a tissue.^{15}^{–}^{17} Most of the reconstruction algorithms are defined under an ideal imaging condition and full-view array of elements. Also, the noise of the measurement system is not considered as a parameter in the reconstruction procedure. Thus, photoacoustic (PA) reconstructed images contain inherent artifacts caused by imperfect reconstruction algorithms. Reducing these artifacts has become a crucial challenge in PA image reconstruction for different number of transducers and different properties of imaging media.^{18}^{,}^{19}

Since there is a high similarity between US and PA detected signals, many of beamforming algorithms used in US imaging can be used in PAI. Moreover, integrating these two imaging modalities has been a challenge.^{20}^{,}^{21} Common US beamforming algorithms such as delay-and-sum (DAS) and minimum variance (MV) can be used in PA beamforming with modifications.^{22} These modifications in algorithms have led to use different hardware to implement an integrated US-PA imaging device. There are many studies focused on developing one beamforming technique for US and PA image formation in order to reduce the cost of imaging system.^{23}^{,}^{24} Although DAS is the most common beamforming method in linear-array imaging, it is a blind beamformer. Consequently, DAS causes a wide mainlobe and high level of sidelobes.^{25} Adaptive beamformers are commonly employed in radar and have the ability of weighting the aperture based on the characteristics of detected signals. Apart from that, these beamformers form a high quality image with a wide range of off-axis signals rejection. MV can be considered as one of the commonly used adaptive methods in medical imaging.^{26}^{–}^{28} Over time, a vast variety of modifications have been investigated on MV such as complexity reduction,^{29}^{,}^{30} shadowing suppression,^{31} using of the eigenstructure,^{32}^{,}^{33} and combination of MV and multiline transmission technique.^{34} Matrone et al. proposed, in Ref. 35, a new beamforming algorithm namely delay-multiply-and-sum (DMAS) as a beamforming technique, used in medical US imaging. This algorithm, introduced by Lim et al., was initially used in confocal microwave imaging for breast cancer detection.^{36} In addition, DMAS was used in synthetic aperture imaging.^{37} Double stage DMAS (DS-DMAS), in which two stages of DMAS is used in order to achieve a higher contrast and resolution compared to DMAS, was proposed for linear-array US and PAI.^{38}^{–}^{40} In addition, a modified version of coherence factor and a high-resolution CF were used to provide a higher contrast and resolution in linear-array PAI, respectively, compared to conventional CF.^{41}^{,}^{42}

In this paper, a beamforming algorithm, namely minimum variance-based DMAS (MVB-DMAS), is introduced. The expansion of the DMAS algorithm is used, and it is shown that in each term of the expansion, there is a DAS algebra. Since the DAS algorithm is a nonadaptive beamformer and leads to low-resolution images, we proposed to use MV instead of the existing DAS in DMAS algebra expansion. It is shown that using MVB-DMAS results in resolution improvement and sidelobe levels reduction at the expense of higher computational burden. A preliminary version of this work and its eigenspace version have been reported before.^{43}^{–}^{45} However, in this paper, we are going to present a highly more complete description of this approach and evaluate, numerically and experimentally, its performance and the effects of its parameters.

The rest of the paper is organized as follows. In Sec. 2, the DMAS and MV beamforming algorithms are presented. In Sec. 3, the proposed method and the necessary modifications are explained. The numerical and experimental results are presented in Secs. 4 and 5, respectively. The advantages and disadvantages of proposed method are discussed in Sec. 6, and the conclusion is presented in Sec. 7.

## 2.

## Background

## 2.1.

### Beamforming

When PA signals are detected by a linear array of US transducer, beamforming algorithms such as DAS can be used to reconstruct the image using the

where ${y}_{\mathrm{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 ${\mathrm{\Delta}}_{i}$ are detected signals and corresponding time delay for detector $i$, respectively. DAS is a simple algorithm and can be used for real-time PA and US imaging. However, in the linear-array transducer only a few numbers of detection angles are available. In other words, a low-quality image is formed due to the limited available angles in linear-array transducers. DMAS was introduced in Ref. 35 to improve the image quality. DMAS calculates corresponding sample for each element of the array, the same as DAS, but before summation, samples are combinatorially coupled and multiplied. The DMAS formula is given as## Eq. (2)

$${y}_{\mathrm{DMAS}}(k)=\sum _{i=1}^{M-1}\sum _{j=i+1}^{M}{x}_{i}(k-{\mathrm{\Delta}}_{i}){x}_{j}(k-{\mathrm{\Delta}}_{j}).$$^{35}

## Eq. (3)

$${\widehat{x}}_{ij}(k)=\mathrm{sign}[{x}_{i}(k-{\mathrm{\Delta}}_{i}){x}_{j}(k-{\mathrm{\Delta}}_{j})]\sqrt{|{x}_{i}(k-{\mathrm{\Delta}}_{i}){x}_{j}(k-{\mathrm{\Delta}}_{j})|}.$$## 2.2.

### Minimum Variance

The output of MV adaptive beamformer is given as

## Eq. (5)

$$y(k)={\mathbf{W}}^{H}(k){\mathbf{X}}_{d}(k)=\sum _{i=1}^{M}{w}_{i}(k){x}_{i}(k-{\mathrm{\Delta}}_{i}),$$## Eq. (6)

$$\mathbf{X}(k)=\mathbf{s}(k)+\mathbf{i}(k)+\mathbf{n}(k)=\mathbf{s}(k)\mathbf{a}+\mathbf{i}(k)+\mathbf{n}(k),$$^{28}To acquire the optimal weights, signal-to-interference-plus-noise ratio (SINR) needs to be maximized

^{46}

## Eq. (7)

$$\mathrm{SINR}=\frac{{\sigma}_{s}^{2}{|{\mathbf{W}}^{H}\mathbf{a}|}^{2}}{{\mathbf{W}}^{H}{R}_{i+n}\mathbf{W}},$$## Eq. (8)

$$\underset{\mathbf{W}}{\mathrm{min}}\text{\hspace{0.17em}\hspace{0.17em}}{\mathbf{W}}^{H}{\mathit{R}}_{i+n}\mathbf{W},\phantom{\rule[-0.0ex]{1em}{0.0ex}}s.t.\text{\hspace{0.17em}\hspace{0.17em}}{v}^{H}\mathbf{a}=1.$$^{47}

## Eq. (9)

$${\mathit{W}}_{\mathrm{opt}}=\frac{{R}_{i+n}^{-1}\mathbf{a}}{{\mathbf{a}}^{H}{R}_{i+n}^{-1}\mathbf{a}}.$$In practical application, the interference-plus-noise covariance matrix is unavailable. Consequently, the sample covariance matrix is used instead of unavailable covariance matrix using $N$ recently received samples and is given as

Using MV in medical US imaging encounters some problems that are addressed in Ref. 27, and we briefly review it here. It should be noticed that by applying delays on each element of the array, the steering vector $\mathbf{a}$ for each signal waveform becomes a vector of ones.^{26}^{–}^{28} The subarray-averaging or spatial-smoothing method can be used to achieve a good estimation of covariance matrix using decorrelation of the coherent signals received by array elements. The covariance matrix estimation using spatial-smoothing can be written as

## Eq. (11)

$$\widehat{R}(k)=\frac{1}{M-L+1}\sum _{l=1}^{M-L+1}{\mathbf{X}}_{d}^{l}(k){\mathbf{X}}_{d}^{l}{(k)}^{H},$$## Eq. (13)

$$\widehat{R}(k)=\frac{1}{(2K+1)(M-L+1)}\times \sum _{n=-K}^{K}\sum _{l=1}^{M-L+1}{\mathbf{X}}_{d}^{l}(k+n){\mathbf{X}}_{d}^{l}{(k+n)}^{H},$$## Eq. (14)

$$\widehat{y}(k)=\frac{1}{M-L+1}\sum _{l=1}^{M-L+1}{\mathbf{W}}_{*}^{H}(k){\mathbf{X}}_{d}^{l}(k),$$## 3.

## Proposed Method

In this paper, it is proposed to use the MV adaptive beamformer instead of the existing DAS algebra inside DMAS mathematical expansion. To illustrate this, consider the expansion of the DMAS algorithm, which can be written as

## Eq. (15)

$${y}_{\mathrm{DMAS}}(k)=\sum _{i=1}^{M-1}\sum _{j=i+1}^{M}{x}_{id}(k){x}_{jd}(k)=[{x}_{1d}(k){x}_{2d}(k)+{x}_{1d}(k){x}_{3d}(k)+\dots +{x}_{1d}(k){x}_{Md}(k)]\phantom{\rule{0ex}{0ex}}+[{x}_{2d}(k){x}_{3d}(k)+{x}_{2d}(k){x}_{4d}(k)+\dots +{x}_{2d}(k){x}_{Md}(k)]+\dots +[{x}_{(M-2)d}(k){x}_{(M-1)d}(k)+{x}_{(M-2)d}(k){x}_{Md}(k)]+[{x}_{(M-1)d}(k){x}_{Md}(k)],$$## Eq. (16)

$${y}_{\mathrm{DMAS}}(k)=\sum _{i=1}^{M-1}\sum _{j=i+1}^{M}{x}_{id}(k){x}_{jd}(k)={x}_{1d}(k)\genfrac{}{}{0ex}{}{\underset{\u23df}{[{x}_{2d}(k)+{x}_{3d}(k)+{x}_{4d}(k)+\dots +{x}_{Md}(k)]}}{\text{first term}}+{x}_{2d}(k)\genfrac{}{}{0ex}{}{\underset{\u23df}{[{x}_{3d}(k)+{x}_{4d}(k)+\dots +{x}_{Md}(k)]}}{\text{second term}}\phantom{\rule{0ex}{0ex}}+\dots +{x}_{(M-2)d}(k)\genfrac{}{}{0ex}{}{\underset{\u23df}{[{x}_{(M-1)d}(k)+{x}_{Md}(k)]}}{(M-2)\text{'}\mathrm{th}\text{\hspace{0.17em}}\text{term}}+\genfrac{}{}{0ex}{}{\underset{\u23df}{[{x}_{(M-1)d}(k).{x}_{Md}(k)]}}{(M-1)\text{'}\mathrm{th}\text{\hspace{0.17em}}\text{term}}.$$## 3.1.

### Modified DMAS

It should be noticed that the quality of covariance matrix estimation in MV highly depends on the selected length of subarray. The upper boundary is limited to $M/2$ and the lower boundary to 1. Choosing $L=M/2$ leads to resolution enhancement at the cost of robustness, and $L=1$ leads to resolution reduction and robustness increment. In Eq. (16), each term can be considered as a DAS algorithm with different number of elements of array. In other words, the number of samples of elements contributing in the existing DAS is different in each term, which results from the nature of the DMAS algorithm. To illustrate this, consider the length of array $M$ and $L=M/2$. There will be $M-1$ terms in DMAS expansion, while first term contains $M-1$ entries, second term contains $M-2$ entries, and finally the last term contains only one entry. Limited number of entries in each term causes problem for MV algorithm due to the limited length of the subarray. This problem can be addressed by adding the unavailable elements in each term in order to acquire large enough number of available elements and consequently high quality covariance matrix estimation. The extra terms, needed to address the problem, are given as

## Eq. (17)

$${y}_{\text{extra}}(k)=\sum _{i=M-2}^{2}\sum _{j=i-1}^{1}{x}_{id}(k){x}_{jd}(k)+{y}_{{\text{extra}}^{*}}\phantom{\rule{0ex}{0ex}}={x}_{(M-2)d}(k)[{x}_{(M-3)d}(k)+{x}_{(M-4)d}(k)+\dots +{x}_{2d}(k)+{x}_{1d}(k)]+{x}_{(M-3)d}(k)[{x}_{(M-4)d}(k)+{x}_{(M-5)d}(k)+\dots +{x}_{2d}(k)+{x}_{1d}(k)]+\dots +{x}_{3d}(k)\phantom{\rule{0ex}{0ex}}[{x}_{2d}(k)+{x}_{1d}(k)]+{x}_{2d}(k){x}_{1d}(k)+{y}_{{\text{extra}}^{*}}(k),$$## Eq. (18)

$${y}_{{\text{extra}}^{*}}(k)={x}_{Md}(k)[{x}_{(M-1)d}(k)+{x}_{(M-2)d}(k)+\dots +{x}_{3d}(k)+{x}_{2d}(k)+{x}_{1d}(k)].$$## Eq. (19)

$${y}_{\mathrm{MDMAS}}(k)={y}_{\mathrm{DMAS}}(k)+{y}_{\text{extra}}(k)\phantom{\rule{0ex}{0ex}}=\sum _{i=1}^{M}\sum _{j=1,j\ne i}^{M}{x}_{id}(k){x}_{jd}(k)=\phantom{\rule{0ex}{0ex}}={x}_{1d}(k)\genfrac{}{}{0ex}{}{\underset{\u23df}{[{x}_{2d}(k)+{x}_{3d}(k)+\dots +{x}_{(M-1)d}(k)+{x}_{Md}(k)]}}{\text{first term}}\phantom{\rule{0ex}{0ex}}+{x}_{2d}(k)\genfrac{}{}{0ex}{}{\underset{\u23df}{[{x}_{1d}(k)+{x}_{3d}(k)+\dots +{x}_{(M-1)d}(k)+{x}_{Md}(k)]}}{\text{second term}}+\dots +{x}_{(M-1)d}(k)\genfrac{}{}{0ex}{}{\underset{\u23df}{[{x}_{1d}(k)+{x}_{2d}(k)+\dots +{x}_{(M-2)d}(k)+{x}_{Md}(k)]}}{(M-1)\text{'}\mathrm{th}}\phantom{\rule{0ex}{0ex}}+{x}_{Md}(k)\genfrac{}{}{0ex}{}{\underset{\u23df}{[{x}_{1d}(k)+{x}_{2d}(k)+\dots +{x}_{(M-2)d}(k)+{x}_{(M-1)d}(k)]}}{M\text{'}\mathrm{th}\text{\hspace{0.17em}}\text{term}}.$$## Eq. (20)

$${y}_{\mathrm{MV}\text{-}\mathrm{DMAS}}(k)=\sum _{i=1}^{M}{x}_{id}(k)[{\mathbf{W}}_{i,M-1}^{H}(k){\mathbf{X}}_{id,M-1}(k)]=\sum _{i=1}^{M}{x}_{id}(k)[\sum _{j=1,j\ne i}^{M}{w}_{j}(k){x}_{jd}(k)]=\sum _{i=1}^{M}{x}_{id}(k)[{\mathbf{W}}^{H}(k){\mathbf{X}}_{d}(k)-{w}_{i}(k){x}_{id}(k)]=\sum _{i=1}^{M}{x}_{id}(k)[\sum _{j=1}^{M}{w}_{j}(k){x}_{jd}(k)-{w}_{i}(k){x}_{id}(k)]=\sum _{i=1}^{M}{x}_{id}(k)\genfrac{}{}{0ex}{}{\underset{\u23df}{[\sum _{j=1}^{M}{w}_{j}(k){x}_{jd}(k)]}}{\mathrm{MV}}-\sum _{i=1}^{M}{x}_{id}(k)[{w}_{i}(k){x}_{id}(k)],$$## Eq. (21)

$${y}_{\mathrm{MV}\text{-}\mathrm{DMAS}}(k)=\sum _{i=1}^{M}{x}_{id}(k)\genfrac{}{}{0ex}{}{\underset{\u23df}{[\sum _{j=1}^{M}{w}_{j}(k){x}_{jd}(k)]}}{MV}-\sum _{i=1}^{M}{x}_{id}(k)[{w}_{i}(k){x}_{id}(k)]=\sum _{i=1}^{M}\genfrac{}{}{0ex}{}{\underset{\u23df}{{x}_{id}(k)\genfrac{}{}{0ex}{}{\underset{\u23df}{[\sum _{j=1}^{M}{w}_{j}(k){x}_{jd}(k)]}}{MV}-{w}_{i}(k){x}_{id}^{2}(k)}}{i\text{'}\mathrm{th}\text{\hspace{0.17em}}\text{term}}.$$## Eq. (22)

$${y}_{\mathrm{MVB}\text{-}\mathrm{DMAS}}(k)=\sum _{i=1}^{M}{w}_{i,\text{new}}\genfrac{}{}{0ex}{}{\underset{\u23df}{\{{x}_{id}(k)[\sum _{j=1}^{M}{w}_{j}(k){x}_{jd}(k)]-{w}_{i}(k){x}_{id}^{2}(k)\}}}{i\text{'}\mathrm{th}\text{\hspace{0.17em}}\text{term}},$$## 4.

## Numerical Results and Performance Assessment

In this section, numerical results are presented to illustrate the performance of the proposed algorithm in comparison with DAS, DMAS, and MV.

## 4.1.

### Simulated Point Target

## 4.1.1.

#### Simulation setup

The $K$-wave MATLAB^{®} toolbox was used to simulate the numerical study.^{48} Eleven 0.1-mm radius spherical absorbers as initial pressure were positioned along the vertical axis every 5 mm beginning 25 mm from the transducer surface. The imaging region was 20 mm in lateral axis and 80 mm in vertical axis. A linear-array having $M=128$ elements operating at 5 MHz central frequency and 77% fractional bandwidth was used to detect the PA signals generated from the defined initial pressures. The speed of sound was assumed to be $1540\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}/\mathrm{s}$ during simulations. The sampling frequency was 50 MHz, subarray length $L=M/2$, $K=5$, and $\mathrm{\Delta}=1/100L$ for all the simulations. Also, a band-pass filter was applied by a Tukey window ($\alpha =0.5$) to the beamformed signal spectra, covering 6 to 16 MHz, to pass the necessary information.

## 4.1.2.

#### Qualitative evaluation

Figures 1(a)–1(d) show the output of DAS, MV, DMAS, and MVB-DMAS beamformers, respectively. It is clear that DAS and DMAS result in low-resolution images, and at the high depths of imaging, both algorithms lead to a wide mainlobe. However, DMAS leads to lower level of sidelobes and a higher resolution. In Fig. 1(b), it can be seen that MV results in high resolution, but the high level of sidelobes affect the reconstructed image. Formed image using MVB-DMAS is shown in Fig. 1(d) where the resolution of MV beamformer is maintained and the level of a sidelobe is highly degraded compared to MV.

To assess the different beamforming algorithms in detail, lateral variations of the formed images are shown in Fig. 2. Lateral variations at the depth of 50 mm are shown in Fig. 2(c) where DAS, MV, DMAS, and MVB-DMAS result in about $-40$, $-50$, $-55$, and $-65\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$ sidelobes, respectively. On the other hand, the width of a mainlobe can be regarded as a parameter, indicating the resolution metric. It can be seen that MV and MVB-DMAS result in significant higher resolution in comparison with DAS and DMAS.

## 4.1.3.

#### Quantitative evaluation

To quantitatively compare the performance of the beamformers, the full-width-half-maximum (FWHM) in $-6\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$ and signal-to-noise ratio (SNR) are calculated in all imaging depths using point targets in the reconstructed images. The results for FWHM and SNR are shown in Tables 1 and 2, respectively. As can be seen in Table 1, MVB-DMAS results in the narrowest $-6\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$ width of a mainlobe in all imaging depths compared to other beamformers. In particular, consider depth of 50 mm where FWHM for DAS, DMAS, MV, and MVB-DMAS is about 3565, 2355, 172, and $95\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, respectively. More importantly, the FWHM differentiation of the first and last imaging depth indicates that MVB-DMAS and MV techniques variate 158 and $378\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, respectively, while DAS and DMAS variate 5425 and $3394\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, respectively. As a result, FWHM is more stabilized using MVB-DMAS and MV in comparison with DAS and DMAS. The represented SNRs in Table 2 are calculated using

## Eq. (23)

$$\mathrm{SNR}=20\text{\hspace{0.17em}}{\mathrm{log}}_{10}{P}_{\text{signal}}/{P}_{\text{noise}},$$^{39}As can be seen in Table 2, MVB-DMAS outperforms other beamformers in SNR. Consider, in particular, the depth of 50 mm where SNR for DAS, DMAS, MV, and MVB-DMAS is 21.7, 32.0, 33.2 and 45.0 dB, respectively.

## Table 1

FWHM (μm) values (in −6 dB) at the different depths.

Depth (mm) | Beamformer | |||
---|---|---|---|---|

DAS | DMAS | MV | MVB-DMAS | |

25 | 1200 | 850 | 93 | 54 |

30 | 1476 | 1059 | 99 | 54 |

35 | 1842 | 1286 | 115 | 61 |

40 | 2277 | 1584 | 133 | 72 |

45 | 2710 | 1862 | 143 | 82 |

50 | 3565 | 2355 | 172 | 95 |

55 | 3800 | 2535 | 187 | 102 |

60 | 4400 | 2937 | 226 | 113 |

65 | 4967 | 3273 | 288 | 123 |

70 | 5512 | 3639 | 305 | 146 |

75 | 6625 | 4244 | 471 | 212 |

## Table 2

SNR (dB) values at the different depths.

Depth (mm) | Beamformer | |||
---|---|---|---|---|

DAS | DMAS | MV | MVB-DMAS | |

25 | 30.7 | 51.1 | 42.7 | 57.7 |

30 | 28.6 | 46.9 | 39.3 | 53.2 |

35 | 26.8 | 43.6 | 37.7 | 50.8 |

40 | 25.2 | 40.2 | 35.3 | 48.2 |

45 | 23.5 | 36.2 | 34.4 | 46.5 |

50 | 21.7 | 32.0 | 33.2 | 45.0 |

55 | 20.1 | 28.2 | 32.4 | 50.8 |

60 | 18.5 | 25.0 | 31.3 | 48.2 |

65 | 17.2 | 22.6 | 30.4 | 46.5 |

70 | 16.3 | 20.7 | 30.1 | 45.0 |

75 | 15.5 | 19.1 | 29.0 | 43.3 |

## 4.2.

### Sensitivity to Sound Velocity Inhomogeneities

In this section, the proposed method is evaluated in the term of robustness against the sound velocity errors resulting from medium inhomogeneities, which are inevitable in practical imaging. The simulation design for Fig. 1 is used in order to investigate the robustness, except that the sound velocity is overestimated by 5%, which covers and may be more than the typical estimation error.^{26}^{,}^{27} It should be noticed that in the simulation, we have intentionally overestimated the sound velocity by 5% to evaluate the proposed method. This is important since we usually face this phenomenon in the practical situations where the medium is inhomogeneous, but the images are reconstructed assuming a sound velocity of $1540\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}/\mathrm{s}$. As can be seen in Fig. 3(b), MV leads to higher resolution compared to DAS, but the high levels of sidelobe and negative effects of overestimated sound velocity still affect the reconstructed image. It should be noticed that the appeared noise and artifacts in Fig. 3(b) are due to the overestimated sound velocity and the temporal averaging with $K=5$. DMAS, in Fig. 3(c), reduces these negative effects, but the resolution is not well enough. As can be seen in Fig. 3(d), MVB-DMAS results in the high negative effects reduction of DMAS and the high resolution of MV. However, the reconstructed image using MVB-DMAS contains more artifacts compared to DMAS, which is mainly as a result of the lower SNR of MVB-DMAS compared to DMAS. Figure 4 shows the lateral variation of the reconstructed images in Fig. 3. As can be seen, MVB-DMAS detects the peak amplitude of point target as well as DAS. The resolution of the formed image using MVB-DMAS is improved in comparison with DAS and DMAS. Moreover, the levels of sidelobe using MVB-DMAS is reduced in comparison with other mentioned beamformers.

## 4.3.

### Effects of Varying L

To evaluate the effects of varying $L$, the proposed method has been implemented using $L=64$, $L=45$, $L=32$, and $L=16$. The lateral variations of the formed images at the depth of 45 mm are presented in Fig. 5. Clearly, increasing $L$ results in a higher resolution and lower level of sidelobes. Moreover, the SNR, in two imaging depths, is presented in Table 3. It is shown that SNR does not significantly vary for different $L$. However, $L=45$ results in higher SNR. In addition, Table 4 shows the calculated FWHM for different amounts of $L$, and it proves that FWHM is reduced with higher $L$.

## Table 3

SNR (dB) values of MVB-DMAS for the different amounts of L.

Depth (mm) | Number of L | |||
---|---|---|---|---|

16 | 32 | 45 | 64 | |

45 | 67.5 | 67.2 | 67.3 | 66.0 |

65 | 61.8 | 61.6 | 62.4 | 61.0 |

## Table 4

FWHM (μm) values of MVB-DMAS (in −6 dB) for the different amounts of L.

Depth (mm) | Number of L | |||
---|---|---|---|---|

16 | 32 | 45 | 64 | |

45 | 480 | 199 | 95 | 59 |

65 | 735 | 259 | 161 | 97 |

## 4.4.

### Effects of Coherence Weighting

The proposed algorithm is also evaluated when it has been combined with CF. To have a fair comparison, all the beamformers are combined with CF weighting. The reconstructed images along with the corresponding lateral variations are shown in Figs. 6 and 7, respectively. As can be seen in Fig. 6, MVB-DMAS + CF results in higher resolution and lower sidelobes compared to other beamformers. The higher resolution of MVB-DMAS + CF is visible compared to DMAS + CF, especially at high depths of imaging. In addition, it is clear that MVB-DMAS + CF reduces the sidelobes compared to MV + CF and improves the target detectability. To have a better comparison, consider Fig. 7 where MVB-DMAS + CF outperforms other beamformers in the terms of the width of a mainlobe in $-6\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$ and sidelobes. It is shown that sidelobes for DAS, DMAS, MV, and MVB-DMAS, when all of them are combined with CF, are about $-120$, $-120$, $-111$, and $-160\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$, respectively, showing the superiority of MVB-DMAS+CF.

## 5.

## Experimental Results

To evaluate the MVB-DMAS algorithm, in this section results of designed experiments are presented.

## 5.1.

### Experimental Setup

A linear-array of PAI system was used to detect the PA waves and the major components of system include an US data acquisition system, Vantage 128 Verasonics (Verasonics, Inc., Redmond, Washington), a Q-switched Nd:YAG laser (EverGreen Laser, double-pulse Nd: YAG system) with a pulse repetition rate of 25 Hz, wavelength 532 nm, and a pulse width of 10 ns. A transducer array (L7-4, Philips Healthcare) with 128 elements and 5.2 MHz central frequency was used as a receiver. A function generator is used to synchronize all operations (i.e., laser firings and PA signal recording). The data sampling rate was 20.8320 MHz. The schematic of the designed system is presented in Fig. 8, and a gelatin-based phantom used as imaging target is shown in Fig. 9, including two blood inclusions to provide optoacoustic properties. The experimental setup for PA linear-array imaging is shown in Fig. 10 where two parallels wire are used as phantom for another experiment. It should be noticed that in all the experiments, surface of the transducer is perpendicular to the imaging targets. Thus, it is expected to see a cross section of the targets. A band-pass filter was applied by a Tukey window ($\alpha =0.5$) to the beamformed signal spectra, covering 6 to 13 MHz, to pass the necessary information.

## 5.2.

### Qualitative Evaluation

The reconstructed images using the phantom shown in Fig. 9 are presented in Fig. 11. Clearly, there are three structures seen in the reconstructed images, Fig. 11, which two of them are blood inclusions, and the first one is because of the small fracture on the upper part of the phantom shown in Fig. 9. As is demonstrated, DAS leads to a low-resolution image having a high level of sidelobe, especially the target at the depth of 35 mm. MV leads to a higher resolution in comparison with DAS, but negative effects of the high level of sidelobes are obvious Fig. 11(b), and the background of the reconstructed image are affected by noise. DMAS enhances the image in the terms of sidelobes and artifacts but still provides a low-resolution image. MVB-DMAS leads to a higher resolution image having lower sidelobes compared to DAS, DMAS, and MV. It is clear that MVB-DMAS provides the high resolution of MV and low sidelobes of DMAS. The line above the targets is due to the PA signal generation at the surface of the phantom (due to the top illumination) and is supposed to look like what is seen in Fig. 11(d), considering the area of illumination, its location, and laser beam profile. The phantom we used was a bit old, and its surface was slightly dried. Hence, when we added the US gel on the top of the sample, there was a rather large impedance mismatch created at the interface between the surface of the sample and the US gel. In our other similar tests (Fig. 12), we observed the similar artifact but much weaker. In Figs. 11(a) and 11(c), the extended version of the line is seen, which is considered as an artifact. DAS and DMAS stretch imaging targets [see the two imaging targets in Figs. 11(a) and 11(c)]. On the other hand, the stretch in MVB-DMAS is much less due to the correlation process and two stages of MV. The reconstructed images for the designed experiment shown in Fig. 10 are shown in Fig. 12. Since the surface of the transducer is perpendicular to the wires, it is expected to see the targets such as points. As is demonstrated in Fig. 12(a), DAS results in low-resolution points, along with high levels of artifacts, especially at the depth of about 30 mm. In Fig. 12(b), MV leads to resolution improvement while the image is still suffers from high level of sidelobes. The reconstructed image using DMAS, shown in Fig. 12(c), contains low level of sidelobes, but the resolution is low. Finally, MVB-DMAS provides an image with characteristics of DMAS and MV, which are reduced sidelobes and high resolution, respectively. Figure 13 demonstrates the lateral variations of the beamformers at the depth of 33 mm of Fig. 12. As shown in the green circle, MVB-DMAS results in a narrower width of a mainlobe and lower sidelobes (see the arrows).

## 5.3.

### Quantitative Evaluation

To compare the experimental results quantitatively, SNR and contrast ratio (CR) metrics are used. Tables 5 and 6 show the calculated SNR and CR for the two targets in the Fig. 11. CR formula is explained in Ref. 35. As can be seen, the calculated metrics show that MVB-DMAS outperforms other beamformers. In other words, it leads to higher SNR and CR.

## Table 5

SNR (dB) values at the different depths using the targets in Fig. 11.

Depth (mm) | Beamformer | |||
---|---|---|---|---|

DAS | DMAS | MV | MVB-DMAS | |

35 | 28.6 | 46.96 | 39.32 | 53.21 |

55 | 24.20 | 45.87 | 37.91 | 50.11 |

## Table 6

CR (dB) values at the different depths using the targets in Fig. 11.

Depth (mm) | Beamformer | |||
---|---|---|---|---|

DAS | DMAS | MV | MVB-DMAS | |

35 | 30.4 | 33.8 | 27 | 36.3 |

55 | 26.3 | 30.1 | 26.2 | 32.5 |

## 5.4.

### In Vivo Imaging

We imaged median antebrachial vein of a 30-year-old Middle Eastern male *in vivo* (see Fig. 14). The illumination was done from side using a single large diameter poly(methyl methacrylate) (PMMA) fiber (10 mm). $\mathrm{L}22\text{-}14v$ transducer was used to collect the PA signals while it was positioned perpendicular to the vein. The institutional review board at Wayne State University (Independent Investigational Review Board, Detroit, Michigan) approved the study protocol, and informed consent was obtained from the individual before enrolment in the study. In the reconstructed images, as shown in Fig. 15, the top and bottom of the antebrachial vein showed up. It can be seen that the reconstructed image using MVB-DMAS has lower sidelobes, noise, and artifacts compared to other methods, and the cross sections (top and bottom) of the vein are more detectable.

## 6.

## Discussion

The main improvement gained by the introduced method is that the high resolution of the MV beamforming algorithm is retained while the level of sidelobes is reduced. PA images reconstructed by DAS beamformer have a low quality, along with high effects of off-axis signals and high sidelobes. This is mainly due to the blindness of DAS. In fact, the DAS algorithm is a procedure in which all contributing samples are treated identically. On the other hand, DMAS beamformer is a nonlinear algorithm and leads to a high level of off-axis signals rejection due to its correlation process. In DMAS beamformer, all the calculated samples are weighted using a linear combination of the received signals. This procedure makes DMAS a nonblind beamforming algorithm, which results in lower effects of off-axis signals and higher contrast reconstructed images compared to DAS. However, the resolution improvement by DMAS is not good enough in comparison with the MV algorithm. In MV beamformer, samples are weighted adaptively resulting significant resolution improvement. However, it leads to a high level of sidelobes. Therefore, we face two types of beamformers, in which one of them (DMAS) results in sidelobes improvement, and the other one (MV) leads to significant resolution enhancement.

The expansion of DMAS algebra shows there are multiple terms and each of them can be interpreted as a DAS with different lengths of array. This could be the source of the low resolution of the DMAS algorithm, and using MV instead of these terms can be an appropriate choice to improve the resolution. However, as shown in Eq. (16), the number of contributing samples in each term of the expansion is different. The length of the subarray in the spatial smoothing highly affects the performance of MV algorithm, and in Eq. (16) there are some terms representing a low length of array and subarray. To address this problem, necessary terms are added to each term, and then MV algorithm is applied on it. The superiority of MV has been proved compared to DAS, and it is expected to have resolution improvement using MV instead of the existing DAS inside the expansion. This method has been used in the introduced algorithm twice to suppress the artifacts and sidelobes of MV. In other words, there are two MV algorithms inside the proposed method, one on the delayed signals and one on the $i$’th term of Eq. (21). The MV implemented on the delayed signals improves the resolution, but since there is another summation procedure interpreting as DAS, shown in Eq. (16), the level of sidelobes and artifacts reduces the image quality. Second, MV is implemented on the $i$’th term of Eq. (21) to use the properties of MV algorithm in order to improve the image quality. It should be noticed that since the expansion of DMAS is used to integrate the MV algorithm for resolution improvement, there are multiplication operations in the introduced algorithm. The same as DMAS, a band-pass filter is needed to only pass the necessary information.^{35} The proposed algorithm adaptively calculates the weights for each samples, which improves the resolution. Since the correlation procedure of DMAS contributes in the proposed method, the sidelobe level of MV is reduced while the resolution is retained due to the existence of MV in the proposed method. MVB-DMAS has been evaluated numerically and experimentally. It should be noticed that the processing time of the proposed method is higher than other mentioned beamformers. Table 7 shows the order of beamformers computations and corresponding processing time. The correlation process of DMAS needs more time compared to DAS, and MV needs time to adaptive calculation of the weights. MVB-DMAS uses two stages of MV algorithm and a correlation procedure, so it is expected to result in higher processing time compared to MV and DMAS. The computational complexity for calculating the weighting coefficients in MVB-DMAS is in the order of $O({L}^{3})$. Considering the fact that $L$ supposed to be a fraction of $M$, the computational complexity is a function of ${M}^{3}$. Given the weighting coefficient, the computational complexity of the reconstruction procedure is a function of $M$, so the bottle neck of the computational burden is ${M}^{3}$, which is the same as a regular MV algorithm. Note that, the complexity of DMAS and DAS are $O({M}^{2})$ and $O(M)$, respectively. Since the MV algorithm is used in the proposed method, twice, the effects of length of $L$ have been investigated, and the results showed that it affects MVB-DMAS the same as it affects MV. The proposed algorithm significantly outperforms DMAS and MV in the terms of resolution and level of sidelobes, respectively, mainly due to having the specifications of DMAS and MV at the same time. In fact, MVB-DMAS uses the correlation process of DMAS to suppress the artifacts and noise, and adaptive weighting of MV to improve the resolution.

## Table 7

Computational operation and processing time (s).

Metric | Beamformer | |||
---|---|---|---|---|

DAS | DMAS | MV | MVB-DMAS | |

Order | $M$ | ${M}^{2}$ | ${M}^{3}$ | ${M}^{3}$ |

Processing time | 1.1 | 10.8 | 90.9 | 187.1 |

## 7.

## Conclusion

In PAI, DAS beamformer is a common beamforming algorithm, capable of real-time imaging due to its simple implementation. However, it suffers from poor resolution and high level of sidelobes. To overcome these limitations, a DMAS algorithm was used. Expanding DMAS formula leads to multiple terms of DAS. In this paper, we introduced a beamforming algorithm based on the combination of MV and DMAS algorithms, called MVB-DMAS. This algorithm was established based on the existing DAS in the expansion of DMAS algebra, and it was proposed to use MV beamforming instead of the existing DAS. Introduced algorithm was evaluated numerically and experimentally. It was shown that MVB-DMAS beamformer reduces the level of sidelobes and improves the resolution in comparison with DAS, DMAS, and MV, at the expense of higher computational burden. Qualitative results showed that MVB-DMAS has the capabilities of DMAS and MV concurrently. Quantitative comparisons of the experimental results demonstrated that the MVB-DMAS algorithm improves CR about 20%, 9%, and 33%, and enhances SNR 89%, 15%, and 35%, with respect to DAS, DMAS, and MV.

## Disclosures

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors. The authors have no potential conflicts of interest to disclose.

## References

## Biography

**Moein Mozaffarzadeh** received his BSc degree in electronics engineering from Babol Noshirvani University of Technology, Mazandaran, Iran, in 2015, and his MSc degree in biomedical engineering from Tarbiat Modares University, Tehran, Iran, in 2017. He is currently a research assistant at the Research Center for Biomedical Technologies and Robotics, Institute for Advanced Medical Technologies, Tehran, Iran. His current research interests include photoacoustic image reconstruction, ultrasound beamforming, and biomedical imaging.

**Ali Mahloojifar** received his BSc degree in electronic engineering from Tehran University, Tehran, Iran, in 1988, and he received his MSc degree in digital electronics from Sharif University of Technology, Tehran, Iran, in 1991. He obtained his PhD in biomedical instrumentation from the University of Manchester, Manchester, United Kingdom, in 1995. He joined the Biomedical Engineering Group at Tarbiat Modares University, Tehran, Iran, in 1996. His research interests include biomedical imaging and instrumentation.

**Mahdi Orooji** received his BSc degree in electrical engineering from the University of Tehran, Iran, and his MSc and PhD degrees in communication systems from the Louisiana State University, Louisiana, USA. He was a postdoctoral fellow in the Case Western Reserve University, Cleveland, Ohio, USA. Currently, he is an assistant professor of bioelectrical engineering at the Tarbiat Modares University, Tehran, Iran. His research interests include the computer aided diagnostics System and machine vision in medical images.

**Karl Kratkiewicz** graduated summa cum laude with his BS degree in biomedical physics from Wayne State University’s Honors College, Detroit, Michigan, USA, in 2017. He is now a biomedical engineering PhD student at Wayne State University. His area of research includes multimodal (ultrasound/photoacoustic/elastography) imaging and photoacoustic neonatal brain imaging.

**Saba Adabi** currently is a research fellow in medical imaging in Mayo Clinic. She received her PhD in Applied Electronics Department in biomedical engineering from Roma Tre University in Rome, Italy, in 2017. She joined Wayne State University as a scholar researcher in 2016 to 2017. Her main research interests include optical coherence tomography, photoacoustic/ultrasound imaging. She is a member of OSA and SPIE.

**Mohammadreza Nasiriavanaki** received his PhD with outstanding achievement in medical optical imaging and computing from the University of Kent in the United Kingdom. His bachelor's and master's degrees with honors are in electronics engineering. He is currently an assistant professor in the Biomedical Engineering, Dermatology and Neurology Departments of Wayne State University and scientific member of Karmanos Cancer Institute. His area of expertise is designing medical devices for skin and brain diseases diagnosis using photoacoustic imaging and optical coherence tomography technologies.