footnote

# Plane-Wave Ultrasound Beamforming Through
Independent Component Analysis

Sobhan Goudarzi, Amir Asif and Hassan Rivaz

###### Abstract

Beamforming in plane-wave imaging (PWI) is an essential step in creating
images with optimal quality. Adaptive methods estimate the apodization weights
from echo traces acquired by several transducer elements. Herein, we formulate
plane-wave beamforming as a blind source separation problem. The output of
each transducer element is considered as a non-independent observation of the
field. As such, beamforming can be formulated as the estimation of an
independent component out of the observations. We then adapt the independent
component analysis (ICA) algorithm to solve this problem and reconstruct the
final image. The proposed method is evaluated on a set of simulation, real
phantom, and in vivo data available from the PWI challenge in medical
ultrasound. The performance of the proposed beamforming approach is also
evaluated in different imaging settings. The proposed algorithm improves
lateral resolution by as much as $36.5\%$ and contrast by $10\%$ as compared
to the classical delay and sum. Moreover, results are compared with other
well-known adaptive methods. Finally, the robustness of the proposed method to
noise is investigated.

###### Index Terms:

Plane-wave imaging, adaptive beamforming, ICA, image quality.

##  I Introduction

Ultrasound imaging experienced a revolution with the introduction of Plane-
Wave Imaging (PWI) in which frame-rate can reach several thousands per second.
In contrast to other techniques, PWI fires all elements of the probe
simultaneously to form a flat wavefront and span the whole region of interest
in a single shot. This technique has been successfully applied to different
applications such as imaging of shear waves, contrast imaging, and Doppler
imaging of blood flow [1]. Having a non-focused transmitted beam, however,
leads to poor resolution and low contrast in PWI. This drawback was addressed
by coherent compounding of images obtained by several insonifications of
different angles [2]. Consequently, there is always a trade-off between image
quality and frame-rate. Hence, beamforming is witnessing a growing attention
in order to enhance the quality of images without sacrificing the frame-rate.

In PWI, beamforming mainly refers to the method of merging the outputs of
different crystal elements. More specifically, it is a weighting function
across the probe which is referred to as apodization. It can also be used
during transmission. Delay-and-sum (DAS) is a classical nonadaptive
beamforming method in which apodization weights for different pixels of the
image are assigned based on the F-number as well as a predefined window shape.
As known from spectral estimation, there is often a trade-off between the
width of main lobe and energy of side lobes of the apodization window. When
measured backscattered signals are directly used to optimize the apodization
weights, the beamforming method is considered adaptive.

Capon or minimum variance (MV) is a well-known adaptive method in which
apodization weights are estimated to minimize the variance of output while
preserving the unity gain in the steering direction [3]. Asl and Mahloojifar
[4] extended eigenspace-based MV (EMV) to better suppress off-axis signals.
The main issue with MV is the estimation of covariance matrix of data, which
is time consuming and limits its real-time applications [5, 6].

Matrone et al. [7] proposed Filtered-delay multiply and sum (F-DMAS). The
pipeline of algorithm contains the multiplication of delayed radio-frequency
(RF) signals followed by summation of signed square roots. Finally, the
beamformed signal is achieved by bandpass filtering to remove the DC
component.

A family of adaptive beamforming algorithms are based on phase coherence.
First, coherence factor (CF) is defined as the ratio between the coherent and
incoherent energy across the aperture [8] and then used as an adaptive weight
following the DAS beamformer to enhance the image quality [9]. CF was
generalized to be computed from Fourier spectra over the aperture of the
delayed channel data and in a range of low spatial frequency region [9].
Subsequently, Camacho et al. [10] used phase information of aperture data to
compute the adaptive correction weight and proposed phased CF (PCF). However,
the estimated correction weights of CF methods can be affected by speckle
noise.

More recently and more specifically for PWI, the MV approach was applied in
[11, 12]. Nguyen and Prager [13] proposed extensions to MV for coherent plane-
wave compounding (CPWC). Beamforming based on compressive sensing for PWI was
introduced in [14, 15, 16, 17]. Dei et al. [18, 19] investigated the
performance of their beamforming method entitled aperture domain model image
reconstruction (ADMIRE) on PWI. Beamforming in Fourier domain on PWI was
introduced in [20]. This was further developed to incorporate coherent
compounding and angular weighting in [21]. Beamforming as a regularized
inverse problem was introduced in [22] and applied at each depth separately.
This point of view was extended in [23] to solve inverse problem for all image
depths jointly. Recently, a statistical interpretation of beamforming entitle
iterative maximum-a-posteriori (iMAP) was introduced in [24].

Herein, we propose a new framework for adaptive plane-wave beamforming wherein
apodization weights are estimated through independent component analysis
(ICA). In the field of US imaging, ICA has been mainly used for clutter
filtering and noise suppression [25, 26, 27, 28, 29, 30]. Recently, ICA was
used as a dimensionality reduction technique to speed up ADMIRE beamforming
[31].

What motivates us to make use of ICA is a principal physical limitation which
is brought about by wave propagation. More specifically, the backscattered
waves from sources at the same time, leading to a single sensory data at the
resulting RF signal. Fortunately, this distance is not the same for other
elements and, therefore, a group of distinct sources are indistinguishable
from the output of each piezoelectric element. Our approach considers the
output of all piezoelectric elements as a set of non-independent observations
of the field and then uses beamforming to extract independent components. As
such, we reconstruct the whole image simultaneously. The performance of the
proposed adaptive plane-wave beamforming is evaluated and compared to other
methods on a set of simulation, phantom, and in vivo data provided by PWI
challenge in medical ultrasound (PICMUS) 2016 [32]. Moreover, we run the
algorithm on another set of publicly available phantom data recorded with an
Alpinion ultrasound machine (Bothell, WA) to test its performance on images
recorded with different settings. Finally, the robustness of proposed approach
is compared to other methods in the presence of noisy observation.

##  II METHODS

Figure 1: Geometrical illustration of the PWI. Notation $\alpha$ is the
steering angle of the incidence wave.

Consider a linear array of $n$ elements, symmetrically distributed on the
$x$-axis, transmitting along the positive $z$-axis (Fig. 1). Let us assume a
plane-wave with angle $\alpha$ spans the domain with a sound speed of $c$. The
backscattered signals received by element $i$ located at $x_{i}$ is denoted by
$\mathbf{h}_{i}(t)$. Without any loss of generality,
$zcos(\alpha)+xsin(\alpha)$ is the transmission distance $d_{t}$ from the
origin of the transmitted plane-wave to an arbitrary point $(x,z)$ in the
region-of-interest (ROI) and $\sqrt{(x-x_{i})^{2}+z^{2}}$ is the receiving
distance $d_{r}$ from $(x,z)$ to the location of element $i$. The RF data
corresponding to $(x,z)$ in $\mathbf{h}_{i}(t)$ can be found by applying the
associated propagation delay as follows (hereafter, capital and bold font
variables represent matrices and vectors, respectively):

$\tau=\frac{d_{t}+d_{r}}{c}\Longrightarrow
R_{i}(x,z)=\mathbf{h}_{i}(t)\mid_{t=\tau},$ (1)

where $R$ is the matrix containing RF data received by the entire transducer
array. Each point $(x,z)$ of the final image can be obtained through the
following weighted summation:

$R(x,z)=\sum_{i=1}^{n}\mathbf{w}(i)R_{i}(x,z),$ (2)

where $\mathbf{w}$ is the apodization window of length $n$. In practice,
however, we utilize dynamic beamforming where the F-number (F) is fixed for
the entire image. Therefore, the number of elements considered for the
reconstruction of each depth of the image $l$ is calculated as follows:

$F=z/l.$ (3)

Our goal is to estimate the apodization window using ICA.

###  II-A Independent Component Analysis

ICA is a framework used to separate signal components mixed in observations
recorded at different transducer elements [33]. Assuming an $n$-dimensional
signal space, i.e., an $n$-dimensional observed data $\mathbf{x}$,
$n$-dimensional independent sources $\mathbf{s}$, and a square transformation
matrix $W$ of size $n\times n$, the mixing model can be written as follows
[33]:

$\mathbf{s}=W\mathbf{x}.$ (4)

With the assumption of having independent and non-Gaussian sources (at the
most one independent Gaussian source is allowed), both of $W$ and $\mathbf{s}$
can be estimated using the ICA algorithm. In practice, the objective function
for ICA estimation can be formulated using different measures of non-
Gaussianity such as kurtosis, negentropy, and mutual information. Moreover, it
is very useful to center and whiten the observations before applying ICA. One
of the most famous algorithms of ICA implementation is FastICA, where a unit
vector $w$ is computed such that the dot product $\mathbf{w}^{T}\mathbf{x}$
maximizes negentropy. FastICA algorithm can be summerized in four steps as
follows [33]:

  1. 1.

Random initialization of vector $\mathbf{w}$.

  2. 2.

$\mathbf{w}_{new}=E\\{\mathbf{x}g(\mathbf{w}^{T}\mathbf{x})\\}-E\\{g^{\prime}(\mathbf{w}^{T}\mathbf{x})\\}\mathbf{w}$

  3. 3.

$\mathbf{w}=\mathbf{w}_{new}/\|\mathbf{w}_{new}\|$

  4. 4.

Return to step $2$ until the direction of $\mathbf{w}$ does not change.

where $E$ denotes the expectation operation. $g$ and $g^{\prime}$ are first
and second derivatives of a nonquadratic nonlinear function $f$, respectively.
It was shown that either of the two functions $f$ is robust for negentropy
estimation [34]:
$f(u)=\frac{1}{a_{1}}\log\cosh a_{1}u$ or $f(u)=-\exp(-u^{2}/2)$
where $1\leq a_{1}\leq 2$.

As discussed in [33], ICA can be considered as a variant of the projection
pursuit algorithm [35], which enables one-by-one estimation of the independent
components. This is an important feature since the computational load of the
method substantially decreases in our case wherein only one of the independent
components is required. The aforementioned algorithm can be easily extended to
estimate several independent components. More details regarding the FastICA
algorithm can be found in [33].

###  II-B Beamforming Using ICA

Figure 2: Steps 1 and 2 of the proposed ICA beamforming. (a) From all crystal
elements, RF data of each pixel ($R_{i}(x,z)$) is extracted by applying the
propagation delays to received backscattered signals ($\mathbf{h}_{i}(t)$).
(b) Central pixels of the image around which the crystals are symmetric are
cropped from each $R_{i}(x,z)$. (c) The matrix of observation $X$ is
constructed by stacking $\bar{\mathbf{r}_{i}}$.

In general, our goal is to reconstruct a high-quality spatial map of target
echogenicity. More specifically, what we do is discretizing the map of
scatterers that leads to pixels. Each pixel corresponds to an averaged tissue
reflectivity function over the extent of the pixel. According to the central
limit theorem, such a pixel has a Gaussian distribution if the number of
diffusion scatters in the medium is significantly high [36]. Therefore, ICA
may appear not suited because it finds independent components that maximize
non-Gaussianity. However, it is shown in [33] that the ICA can still be used
even if only one of the independent components is Gaussian. Therefore, we look
for the discretized tissue reflectivity function (i.e. the map of scatterers)
as the only independent Gaussian source.

When the RF data corresponding to each pixel of the final image is extracted
from the output of each element, there is also the possibility of source
ambiguity. More specifically, the backscattered waves of at least two
different scatterers arrive simultaneously. Without loss of generality, when
$\alpha=0$, the backscattered waves of two distinct scatterers (with indices 1
and 2) arrive at the same time in element $i$ if and only if they have the
same propagation delay $\tau$. Form Eq. (1) and if the first scatterer is in
the lateral position $x_{i}$, it can be written that:

$\tau_{1}=\tau_{2}\Longrightarrow
2z_{1}=z_{2}+\sqrt{(x_{2}-x_{i})^{2}+z_{2}^{2}}.$ (5)

Simplifying (5) gives:

$z_{2}=z_{1}-\frac{(x_{2}-x_{i})^{2}}{4z_{1}}.$ (6)

So, for $z_{2}<z_{1}$, there are a bunch of scatterers located on an ellipse,
whose reflections arrive at the same time with scatterer $1$. In the
continuous case, this problem is fully addressable. In the discrete case,
however, there is the error of quantization as well. Although this problem was
shown for the specific case of $\alpha=0$, the concept can be extended for
different angles.

As can be seen from Eq. (6), the group of scatterers from whom their
reflections arrive simultaneously are not the same for each crystal element.
This point provides the opportunity of source separation. Herein, the task is
to extract an independent RF signal out of a set of non-independent
observations. In the ideal case, the desired independent RF signal contains
the trace of only one distinct scatterer in each sample. It has to be
mentioned that in practice, the axial and lateral resolution are based on the
sampling frequency of the system, the center frequency of the transmitted
wave, and the transducer design. So, when we refer to one scatterer, it means
one pixel of the final image.

In ultrasound beamforming, an issue is that the apodization window is not
fixed throughout the image. More specifically, ICA works with a fixed
transformation matrix $W$ in Eq. (4). In ultrasound images, however, the
apodization weight is not spatially invariant, rendering a different set of
values of $W$ for different pixels. Two points make the apodization weight
spatially variant. First, for pixels lying at the two lateral ends of the
image, there are crystal elements predominately lying along one side. Second,
as explained in section II, pixels of different depths of the image are
reconstructed using a different number of elements to hold F-number fixed in
the entire image. Hence, if we do not consider these points, ICA fails to
estimate the source and apodization windows, leading to images that are even
lower in quality than DAS.

To solve the aforementioned problem, we consider only the central pixels of
the image around which the crystals are symmetric as the input to the ICA
algorithm. In this way, the cropped portion of $R_{i}$ (Eq. 1) is used to
construct the observation matrix $X$. In our ICA formulation, therefore, the
observations are RF data corresponding to central pixels of final image that
are recorded by all crystal elements of the probe.

Our proposed beamforming algorithm for PWI using ICA includes the following
steps:

  1. 1.

The propagation delay is applied to all channels of data in order to generate
$R_{i},i=0,...,n-1$.

  2. 2.

Each $R_{i}$ is considered as an observation of the field (Fig. 2,). First, it
is cropped (denoted by $\bar{R}_{i}(\bar{x},z)$) and then vectorized in a row
vector ($\bar{\mathbf{r}_{i}}$). Consequently, the observation matrix $X$ is
constructed by stacking the row vectors ($\bar{\mathbf{r}_{i}}$). Finally,
matrix $X$ is centralized and whitened before running ICA.

  3. 3.

The independent source and corresponding mixing vector are estimated using the
FastICA algorithm [37] by maximizing Negentropy as the measure of non-
Gaussianity.

  4. 4.

The apodization window (the estimated transformation matrix in the last step)
is applied throughout the image based on a predefined F-number.

Figure 3: Beamforming results on the single $0^{\circ}$ plane wave. Columns
indicate different image data sets while rows correspond to beamforming
methods.

##  III EXPERIMENTS

###  III-A Dataset

Herein, we use a publicly available PWI dataset, entitled PICMUS, which was
provided by the IEEE International Ultrasonics Symposiun (IUS 2016) in order
to benchmark novel beamforming methods [32]. The PICMUS data utilized in this
work include:

  1. 1.

Simulation resolution (SR): a simulated ultrasound image containing point
targets distributed vertically and horizontally over an anechoic background
designed to assess the performance of beamforming methods in terms of spatial
resolution.

  2. 2.

Simulation contrast (SC): a simulated ultrasound image containing anechoic
cysts distributed vertically and horizontally over fully developed speckle
designed to assess the performance of beamforming methods in terms of
contrast.

  3. 3.

Experimental Resolution (ER): an experimental ultrasound image was recorded on
a CIRS Multi-Purpose Ultrasound Phantom (Model 040GSE) in the regions
containing several wires against speckle background to assess the performance
of beamforming methods in terms of spatial resolution.

  4. 4.

Experimental contrast (EC): an experimental ultrasound image was recorded on
the same phantom as ER but in the regions containing two anechoic cysts
against speckle background to assess the performance of beamforming methods in
terms of contrast.

In addition, PICMUS dataset also contains two in vivo ultrasound images,
showing cross-sectional and longitudinal views, recorded on the carotid artery
of a volunteer. All of the phantom and in vivo data were collected using a
Verasonics Vantage 256 research scanner and a L11 probe (Verasonics Inc.,
Redmond, WA). The simulation settings were selected to be as similar as
possible to the experimental setup.

For each of mentioned groups, a collection of received prebeamformed data
corresponding to 75 steered Plane-Waves covering the angle span from
$-16^{\circ}$ to $16^{\circ}$ was provided. Both RF and IQ (phase quadrature)
formats of data were provided. The proposed algorithm works on the RF version
of data. More details regarding PICMUS dataset can be found in [32].

Another publicly available dataset on PWI is used to investigate the
performance of the proposed method on data collected with other imaging
settings. In this database, two datasets of CPWC recorded using an Alpinion
scanner with a L3-8 probe from a CIRS phantom are used. These data sets are
available through the ultrasound toolbox [38], a MATLAB toolbox for processing
ultrasonic signals. The first dataset was recorded from regions containing
hyperechoic cyst and points scatterers, and the second one recorded from
regions containing hypoechoic cyst. Imaging setting such as probe geometry,
transmit frequency, sampling frequency, and etc. are completely different from
PICMUS. More details regarding this data can be found in https://www.ustb.no/.

###  III-B Evaluation Metrics

Contrast and resolution are considered for the sake of evaluation. More
specifically, resolution is estimated as the full width at half maximum (FWHM)
both in axial and lateral directions. The average value of FWHM among all
scatterers in the image is reported. As for contrast, the contrast to noise
ratio (CNR) is calculated as follows [32]:

$\text{CNR}=20\log_{10}(\frac{\mid\mu_{in}-\mu_{out}\mid}{\sqrt{(\sigma_{in}^{2}+\sigma_{out}^{2})/2}}),$
(7)

where $\mu$ and $\sigma$ are the mean gray level and the gray level standard
deviation. Subscribes $._{in}$ and $._{out}$ refer to inside and outside of
the anechoic cystic region, respectively. Indexes are calculated on B-Mode
images. In order to unify the comparison, we use the codes provided by PICMUS
to compute the indexes.

###  III-C Implementation Details

Figure 4: ICA beamforming using 1, 11, and 75 plane waves. Columns indicate
different image data sets and rows correspond to the number of transmitted
plane waves.

As explained in Section LABEL:sec:sec23, the FastICA algorithm is used to
estimate the apodization window. The maximum number of iterations is set to
$100$ and the stopping criterion is set to be $\epsilon=10^{-6}$. The weights
are initialized with random numbers extracted from standard distribution. The
reduction of dimension through PCA is not used and the best results which are
most reproducible are attained by considering all eigenvalues in the
estimation procedure. We use the Matlab implementation of the Fixed point ICA,
the main algorithm of FastICA, which is publicly available online
http://research.ics.aalto.fi/ica/fastica/code/dlcode.shtml.

Throughout the results section, we consider $F=1.75$ and use Tukey (tapered
cosine) window with constant parameter set to $0.25$ for DAS and other
adaptive methods on top of DAS.

It is not possible to theoretically prove the convergence of FastICA algorithm
with the mentioned parameters. In practice, however, we set the maximum number
of iterations equal to $100$ and observe that for all of data sets, the
algorithm converges in a lower number of iterations.

TABLE I: Quantitative results in terms of CNR and FWHM indexes for simulation and real phantom experiments. dataset |  | SR | ER | SC | EC
---|---|---|---|---|---
index |  |  FWHMA FWHML |  FWHMA FWHML | CNR | CNR
1 PW |  DAS ICA |  0.4 0.82 0.39 0.52 |  0.57 0.88 0.57 0.81 |  9.95 10.67 |  8.15 8.9
11 PW |  DAS ICA |  0.4 0.54 0.4 0.41 |  0.56 0.54 0.56 0.51 |  12.48 12.6 |  11.25 11.4
75 PW |  DAS ICA |  0.4 0.56 0.4 0.42 |  0.56 0.56 0.56 0.53 |  15.55 15.96 |  12 12.1

##  IV RESULTS

###  IV-A Simulated and Experimental Data

The results of DAS beamforming versus our proposed method based on ICA on a
single $0^{\circ}$ plane wave of simulated and experimental data are
illustrated in Fig. 3. As seen from this figure, the proposed beamforming
method outperforms DAS and improves the resolution as well as contrast on both
simulated and experimental phantom data. In order to better investigate the
amount of improvement, the quantitative indices are reported in Table I. What
causes the improvement is the window used for apodization. So, as Table I
confirms, improvement in resolution can only be acquired in the lateral
direction. The highest improvement in lateral FWHM is $36.5\%$ on simulated
plane-wave data of only one single transmission. In terms of CNR,
approximately $10\%$ of improvement is achieved on the simulated cyst data of
a single transmission while boarders of the cyst are also sharper.

Figure 5: Comparison of beamforming results using different number of plane
waves. Rows indicate beamforming methods while columns correspond to the
number of transmitted plane waves.

In order to investigate the effect of CPWC, the results of the proposed method
on higher number of plane waves are illustrated in Fig. 4. The indexes of
Table I as well as Fig. 4 confirm that CPWC improves the image quality. As for
CPWC, we do not repeat beamforming for each angle and use the apodization
weights of $0^{\circ}$ plane wave for the remaining angles. Moreover, we do
not apply any angular apodization weights to limit the sources of achieved
improvement. In fact, our main focus is on beamforming of the received
signals.

To better understand the effect of proposed method, Fig. 6 shows a comparison
between Tukey25 window used in DAS and the apodization weights estimated by
ICA on ER dataset. The estimated window has a lower leakage factor as well as
a relative side lobe attenuation while its main lobe is wider. The estimated
window is of a different shape which can not be found among predefined common
windows. So, this point confirms the necessity of estimating the apodization
window from the received data.

Fig. 5 demonstrates qualitative improvements with ICA using fewer angles than
needed with DAS. As can be seen in Fig. 5, the proposed approach achieves
similar image quality with $11$ plane waves compared with DAS using $75$
angles. Therefore, it is possible to reduce the number of plane wave transmits
needed to achieve image quality similar to a fully sampled transmit.

Figure 6: Comparison of apodization window estimated using ICA (second row)
and Tukey25 (first row) used in DAS. Windows are shown in both space and
frequency domains.

###  IV-B In vivo Data

In real ultrasound tissues, there are more sources of degradation in image
quality. In order to make sure that the proposed method also works on in vivo
data, the results of beamforming on real carotid images of PICMUS dataset are
provided in Fig. 7. Visual comparison of beamformed images with different
number of angles reveals that ICA outperforms classical DAS in both cross as
well as longitudinal sections. Furthermore, ICA results in sharper images with
a better contrast.

Figure 7: Beamforming results on in vivo data using 1, 11, and 75 plane waves.
Two columns in left indicate cross-section images while left two columns
correspond to longitudinal-section. Rows denote different number of
transmitted plane waves used in beamforming.

Figure 8: Beamforming results on the single $0^{\circ}$ plane wave data
collected with Alpinion Scanner. Columns indicate different image data sets
while rows correspond to beamforming methods.

Figure 9: Results of other adaptive beamforming methods on the single
$0^{\circ}$ plane wave. Columns indicate different image data sets while rows
correspond to different adaptive beamforming methods.

###  IV-C Results on Alpinion Scanner

We investigated the performance of our method on data collected with Alpinion
Scanner thanks to ultrasound toolbox [38]. Fig. 8 illustrates the result of
proposed method on images of cyst regions and points scatterers recorded on an
Alpinion scanner with a L3-8 probe from a CIRS phantom. The ICA results in
Fig. 8 obviously have a better quality as compared to DAS results in terms of
resolution and contrast. As illustrated in Fig. 8, point targets in the image
beamformed with ICA are finer in ER dataset, and contrast of cyst regions in
EC dataset is noticeably improved. As we do not have the exact location of
cysts and point targets in this dataset, quantitative comparison is not
possible.

###  IV-D Comparison with Other Adaptive Methods

As mentioned before, our focus in current study is on beamforming of the
received signals. So, comparison with other adaptive approaches is of crucial
importance. In this way, we present the results of five well-known approaches,
namely MV [3], EMV [4], CF [9], generalized CF (GCF) [9], and PCF [10]. The
comparison with these methods was not possible without using codes provided by
Rindal et al. [39] in ultrasound toolbox repository (
http://www.ustb.no/publications/dynamic_range/). As for comparison with F-DMAS
method, the sampling frequency of the PICMUS data was not large enough to
apply the method with filtering around 2 times the center frequency. Fig. 9
shows the result of different adaptive beamforming algorithms on a single
$0^{\circ}$ plane wave of simulated and experimental data. The quantitative
comparison is provided in Table II. The EMV method outperforms all other
methods, even our proposed method, in terms of indices. However, methods based
on the MV are very time consuming and are not practical for online
applications. More specifically, the run time of our proposed method is in
order of millisecond while EMV takes a few minutes. Furthermore, as it can be
seen in Fig. 9, the EMV method destroys the image texture and its results are
not visually appealing. Other approaches based on CF outperform the proposed
approach in terms of FWHM index while are worse in terms of contrast.

###  IV-E Robustness to Noise

In order to further demonstrate the superiority of our proposed method, the
quality of beamforming approaches is investigated when some of the receiving
channels are noisy. More specifically, we add independent Gaussian noise with
different levels of signal to noise ratio (SNR) to different number of
channels of ER dataset to see how noise affects the performance of different
methods.

The first three rows of Fig. 10 illustrate the results for additive Gaussian
noise with $-10$, $-20$, and $-40dB$ SNR levels with noise added to channel
$64$ (the total number of channels is $128$). As seen in Fig. 10, other
adaptive methods and DAS fail to reconstruct the affected region with noise
while the proposed beamforming algorithm is completely robust. The fourth and
fifth rows of Fig. 10 correspondingly show the results when $-40dB$ additive
noise is added to 3 and 5 of crystal elements. As the number of noisy channels
with a bad SNR goes up, the quality of ICA starts to decrease. However, it
still outperforms other approaches.

Figure 10: Results of beamforming methods on the single $0^{\circ}$ plane wave of noisy ER dataset. Columns indicate different beamforming methods while three first rows correspond to different levels of noise added to channel $64$. Rows 4 and 5 correspondingly depict the results when we have 3 and 5 noisy channels with $-40dB$ SNR. TABLE II: Quantitative results of other adaptive beamforming methods in terms of CNR and FWHM indexes for simulation and real phantom experiments. dataset |  | SR | ER | SC | EC
---|---|---|---|---|---
index |  |  FWHMA FWHML |  FWHMA FWHML | CNR | CNR
1 PW |  MV EMV CF GCF PCF |  0.41 0.1 0.39 0.09 0.32 0.44 0.32 0.43 0.29 0.37 |  0.59 0.43 0.58 0.33 0.48 0.47 0.48 0.47 0.46 0.41 |  11.1 14 8.2 8.1 6.9 |  7.95 10.5 6.3 6.3 5.2
11 PW |  MV EMV CF GCF PCF |  0.43 0.1 0.4 0.09 0.37 0.37 0.38 0.36 0.37 0.3 |  0.59 0.29 0.56 0.28 0.55 0.37 0.55 0.37 0.55 0.31 |  11.4 15.2 11.9 11.8 11 |  9.8 11.5 10.2 10.2 9.05
75 PW |  MV EMV CF GCF PCF |  0.43 0.1 0.4 0.09 0.4 0.38 0.4 0.38 0.39 0.29 |  0.58 0.31 0.56 0.29 0.56 0.38 0.56 0.38 0.56 0.32 |  14.7 17 14.05 13.9 14.13 |  11 10.4 10 10 10.3

##  V DISCUSSION

Using a part of samples of each channel which only correspond to the middle
part of the final image is important from two aspects. First, it removes the
effect of incomplete data of borders on the ICA performance. Second, the
FastICA algorithm converges faster as it works with a lower amount of data.
Note that the estimated apodization weights and its specifications such as
width of main lobe or the amount of side lobe attenuation in each dataset are
different. So, it can be concluded that there is not a unique solution that
works for all data.

The algorithm can be applied for each angle separately. However, the estimated
weights for different insonification angles are not totally different and the
improvement is negligible while the processing time is multiplied
corresponding to the number of angles. The angular apodization can also be
performed using ICA for CPWC. However, the main focus of this study was
apodization of the received signals. The angular weights are not used to limit
the sources of improvement, which make the comparison with other approaches
possible.

As each adaptive beamforming method alters the dynamic range of the final
image differently, we normalized the range of final outputs to plot the
results in a uniform manner and make the visual comparison possible.

Robustness of the proposed method to noise can be explained as the ability of
ICA in extracting independent components and discarding noise. This aspect of
the proposed beamforming algorithm is very important in practice when we have
a missing or noisy channel. Although a malfunctioning element affects both
transmit and receive, the available data does not give us the chance of a
perfect simulation of this case. In future, we plan to test the performance of
the proposed method for the case of having broken channels in both transmit
and receive.

##  VI CONCLUSIONS

We have proposed a new beamforming approach for ultrasound plane-wave imaging
based on ICA. Beamforming has been formalized as the estimation of one
independent image out of several non-independent observation and the
apodization weights have been estimated based on collected data. The images of
one single plane-wave transmission as well as multiangle plane-wave
acquisitions have been successfully reconstructed. Moreover, the performance
of the algorithm has been demonstrated on different imaging settings. Results
show that the proposed method simultaneously improves the resolution and
contrast while the resulting image is also visually appealing. Furthermore,
the proposed algorithm is strongly robust when some of channels are noisy.

## Acknowledgement

This project was funded by NSERC Discovery grants RGPIN-04136 and
RGPIN2017-06629. The authors would like to thank the organizers of the PICMUS
challenge as well as the ultrasound toolbox for providing publicly available
codes and data.

## References

  * [1] M. Tanter and M. Fink,  “Ultrafast imaging in biomedical ultrasound,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 61, no. 1, pp. 102–119, January 2014.
  * [2] G. Montaldo, M. Tanter, J. Bercoff, N. Benech, and M. Fink,  “Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 56, no. 3, pp. 489–506, March 2009.
  * [3] J. Synnevag, A. Austeng, and S. Holm,  “Benefits of minimum-variance beamforming in medical ultrasound imaging,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 56, no. 9, pp. 1868–1879, Sep. 2009.
  * [4] B. M. Asl and A. Mahloojifar,  “Eigenspace-based minimum variance beamforming applied to medical ultrasound imaging,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 57, no. 11, pp. 2381–2390, November 2010.
  * [5] J. F. Synnevag, A. Austeng, and S. Holm,  “Adaptive beamforming applied to medical ultrasound imaging,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 54, no. 8, pp. 1606–1613, August 2007.
  * [6] A. C. Jensen and A. Austeng,  “An approach to multibeam covariance matrices for adaptive beamforming in ultrasonography,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 59, no. 6, pp. 1139–1148, June 2012.
  * [7] G. Matrone, A. S. Savoia, G. Caliano, and G. Magenes,  “The delay multiply and sum beamforming algorithm in ultrasound b-mode medical imaging,”  IEEE Transactions on Medical Imaging, vol. 34, no. 4, pp. 940–949, April 2015.
  * [8] Raoul Mallart and Mathias Fink,  “Adaptive focusing in scattering media through sound‐speed inhomogeneities: The van cittert zernike approach and focusing criterion,”  The Journal of the Acoustical Society of America, vol. 96, no. 6, pp. 3721–3732, 1994.
  * [9] Pai-Chi Li and Meng-Lin Li,  “Adaptive imaging using the generalized coherence factor,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 50, no. 2, pp. 128–141, Feb 2003.
  * [10] J. Camacho, M. Parrilla, and C. Fritsch,  “Phase coherence imaging,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 56, no. 5, pp. 958–974, May 2009.
  * [11] Nghia Q. Nguyen and Richard W. Prager,  “Minimum variance beamformers for coherent plane-wave compounding,”  in Medical Imaging 2017: Ultrasonic Imaging and Tomography, Neb Duric and Brecht Heyde, Eds. International Society for Optics and Photonics, 2017, vol. 10139, pp. 254 – 263, SPIE.
  * [12] O. M. H. Rindal and A. Austeng,  “Double adaptive plane-wave imaging,”  in 2016 IEEE International Ultrasonics Symposium (IUS), Sep. 2016, pp. 1–4.
  * [13] N. Q. Nguyen and R. W. Prager,  “A spatial coherence approach to minimum variance beamforming for plane-wave compounding,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 65, no. 4, pp. 522–534, April 2018.
  * [14] Congzhi Wang, Xi Peng, Dong Liang, and Hairong Zheng,  “Plane-wave ultrasound imaging based on compressive sensing with low memory occupation,”  in 2015 IEEE International Ultrasonics Symposium (IUS), Oct 2015, pp. 1–4.
  * [15] Guillaume David, Jean-luc Robert, Bo Zhang, and Andrew F Laine,  “Time domain compressive beam forming of ultrasound signals,”  The Journal of the Acoustical Society of America, vol. 137, no. 5, pp. 2773–2784, 2015.
  * [16] A. Besson, M. Zhang, F. Varray, H. Liebgott, D. Friboulet, Y. Wiaux, J. Thiran, R. E. Carrillo, and O. Bernard,  “A sparse reconstruction framework for fourier-based plane-wave imaging,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 63, no. 12, pp. 2092–2106, 2016.
  * [17] A. Besson, D. Perdios, F. Martinez, Z. Chen, R. E. Carrillo, M. Arditi, Y. Wiaux, and J. Thiran,  “Ultrafast ultrasound imaging as an inverse problem: Matrix-free sparse image reconstruction,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 65, no. 3, pp. 339–355, 2018.
  * [18] Kazuyuki Dei, Jaime Tierney, and Brett Byram,  “Aperture domain model image reconstruction (ADMIRE) with plane wave synthesis,”  in Medical Imaging 2017: Ultrasonic Imaging and Tomography, Neb Duric and Brecht Heyde, Eds. International Society for Optics and Photonics, 2017, vol. 10139, pp. 244 – 253, SPIE.
  * [19] Kazuyuki Dei, Jaime E. Tierney, and Brett C. Byram,  “Model-based beamforming with plane wave synthesis in medical ultrasound,”  Journal of Medical Imaging, vol. 5, no. 2, pp. 1 – 12, 2018.
  * [20] M. Albulayli and D. Rakhmatov,  “Fourier domain depth migration for plane-wave ultrasound imaging,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 65, no. 8, pp. 1321–1333, Aug 2018.
  * [21] C. Chen, G. A. G. M. Hendriks, R. J. G. van Sloun, H. H. G. Hansen, and C. L. de Korte,  “Improved plane-wave ultrasound beamforming by incorporating angular weighting and coherent compounding in fourier domain,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 65, no. 5, pp. 749–765, May 2018.
  * [22] T. Szasz, A. Basarab, and D. Kouamé,  “Beamforming through regularized inverse problems in ultrasound medical imaging,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 63, no. 12, pp. 2031–2044, Dec 2016.
  * [23] E. Ozkan, V. Vishnevsky, and O. Goksel,  “Inverse problem of ultrasound beamforming with sparsity constraints and regularization,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 65, no. 3, pp. 356–365, March 2018.
  * [24] T. Chernyakova, D. Cohen, M. Shoham, and Y. C. Eldar,  “imap beamforming for high-quality high frame rate imaging,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 66, no. 12, pp. 1830–1844, 2019.
  * [25] R. R. Wildeboer, F. Sammali, R. J. G. Van Sloun, Y. Huang, P. Chen, M. Bruce, C. Rabotti, S. Shulepov, G. Salomon, B. C. Schoot, H. Wijkstra, and M. Mischi,  “Blind source separation for clutter and noise suppression in ultrasound imaging: Review for different applications,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, pp. 1–1, 2020.
  * [26] A. C. H. Yu and L. Lovstakken,  “Eigen-based clutter filter design for ultrasound color flow imaging: a review,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 57, no. 5, pp. 1096–1111, 2010.
  * [27] Jaime E Tierney, Don M Wilkes, and Brett C Byram,  “Independent component analysis-based tissue clutter filtering for plane wave perfusion ultrasound imaging,”  in Medical Imaging 2019: Ultrasonic Imaging and Tomography. International Society for Optics and Photonics, 2019, vol. 10955, p. 1095503.
  * [28] Caterina M Gallippi, Kathryn R Nightingale, and Gregg E Trahey,  “Bss-based filtering of physiological and arfi-induced tissue and blood motion,”  Ultrasound in medicine & biology, vol. 29, no. 11, pp. 1583–1592, 2003.
  * [29] P. Li, X. Yang, D. Zhang, and Z. Bian,  “Adaptive clutter filtering based on sparse component analysis in ultrasound color flow imaging,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 55, no. 7, pp. 1582–1596, 2008.
  * [30] Caterina M Gallippi and Gregg E Trahey,  “Adaptive clutter filtering via blind source separation for two-dimensional ultrasonic blood velocity measurement,”  Ultrasonic imaging, vol. 24, no. 4, pp. 193–214, 2002.
  * [31] K. Dei, S. Schlunk, and B. Byram,  “Computationally efficient implementation of aperture domain model image reconstruction,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 66, no. 10, pp. 1546–1559, 2019.
  * [32] H. Liebgott, A. Rodriguez-Molares, F. Cervenansky, J. A. Jensen, and O. Bernard,  “Plane-wave imaging challenge in medical ultrasound,”  in 2016 IEEE International Ultrasonics Symposium (IUS), Sep. 2016, pp. 1–4.
  * [33] A. Hyvärinen and E. Oja,  “Independent component analysis: algorithms and applications,”  Neural Networks, vol. 13, no. 4, pp. 411 – 430, 2000.
  * [34] Aapo Hyvärinen,  “New approximations of differential entropy for independent component analysis and projection pursuit,”  in Advances in neural information processing systems, 1998, pp. 273–279.
  * [35] J. H. Friedman and J. W. Tukey,  “A projection pursuit algorithm for exploratory data analysis,”  IEEE Transactions on Computers, vol. C-23, no. 9, pp. 881–890, Sep. 1974.
  * [36] Robert F Wagner,  “Statistics of speckle in ultrasound b-scans,”  IEEE Trans. Sonics & Ultrason., vol. 30, no. 3, pp. 156–163, 1983\.
  * [37] A. Hyvarinen,  “Fast and robust fixed-point algorithms for independent component analysis,”  IEEE Transactions on Neural Networks, vol. 10, no. 3, pp. 626–634, May 1999.
  * [38] A. Rodriguez-Molares, O. M. H. Rindal, O. Bernard, A. Nair, M. A. L. Bell, H. Liebgott, and A. Austeng,  “The ultrasound toolbox,”  in 2017 IEEE International Ultrasonics Symposium (IUS), Sep. 2017, pp. 1–4.
  * [39] O. M. H. Rindal, A. Austeng, A. Fatemi, and A. Rodriguez-Molares,  “The effect of dynamic range alterations in the estimation of contrast,”  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 66, no. 7, pp. 1198–1208, July 2019.