# Ultrasound Speckle Suppression and Denoising using MRI-derived Normalizing
Flow Priors

Vincent van de Schaft and Ruud J.G. van Sloun  V. van de Schaft and R.J.G. van
Sloun are with the department of Electrical Engineering, Eindhoven University
of Technology, Eindhoven, The Netherlands. (e-mails:
<EMAIL_ADDRESS>r.j.g.v.sloun@tue.nl). R.J.G. van Sloun is also
with Philips Research, Eindhoven.

###### Abstract

Ultrasonography offers an inexpensive, widely-accessible and compact medical
imaging solution. However, compared to other imaging modalities such as CT and
MRI, ultrasound images notoriously suffer from strong speckle noise, which
originates from the random interference of sub-wavelength scattering. This
deteriorates ultrasound image quality and makes interpretation challenging. We
here propose a new unsupervised ultrasound speckle reduction and image
denoising method based on maximum-a-posteriori estimation with deep generative
priors that are learned from high-quality MRI images. To model the generative
tissue reflectivity prior, we exploit normalizing flows, which in recent years
have shown to be very powerful in modeling signal priors across a variety of
applications. To facilitate generaliation, we factorize the prior and train
our flow model on patches from the NYU fastMRI (fully-sampled) dataset. This
prior is then used for inference in an iterative denoising scheme. We first
validate the utility of our learned priors on noisy MRI data (no prior domain
shift), and then turn to evaluating performance on both simulated and in-vivo
ultrasound images from the Plane-wave Imaging Challenge in Medical UltraSound
(PICMUS) and Challenge on ultrasound beamforming with deep learning (CUBDL)
datasets. The results show that the method outperforms other (unsupervised)
ultrasound denoising methods (nonlocal means (NLM) and Optimized Bayesian
nonlocal means (OBNLM)) both quantitatively and qualitatively.

###### Index Terms:

ultrasound, deep generative models, normalizing flows, denoising, speckle, MRI

##  I Introduction

Ultrasound (US) imaging offers a low-cost and widely accessible diagnostic
imaging solution, with form factors becoming ever more compact and handheld.
This portability and cost-effectiveness enables point-of-care imaging at the
bedside, in emergency settings, rural clinics, and developing countries. US
sees increasing adoption across many medical specialties, ranging from
cardiology to obstetrics and oncology. Unfortunately, US image quality suffers
from many sources of image degradation, one of which being speckle noise. The
resulting pseudo-random fluctuations in intensity obscure important image
features and make US images hard to interpret.

Many previous works have found ways to combat speckle noise. The Lee filter
[1] uses local statistics to denoise individual pixels in the presence of
multiplicative and/or additive noise. Later methods include the Fost-filter,
the Kuan-filter, and speckle reducing anisotropic diffusion [2, 3, 4]. The
nonlocal-means algorithm uses statistics from patches that are similar to the
pixel to realize more smoothing while retaining sharp edges [5]. This approach
was later adapted specifically for combating ultrasound speckle with the
OBNLM-algorithm [6, 7].

More recently, inspired by the success of deep learning in computer vision,
deep neural network based approaches have been proposed for US denoising.
These data-driven methods act as regressors that learn a parameterized mapping
from images with speckle to corresponding clean images [8, 9], or from raw
channel data (i.e. before beamforming) to clean images [10]. Deep networks for
US image quality improvement are typically based on discriminative
convolutional architectures (such as the well-known U-Net), and often resort
to pairs of clean and artificially corrupted images obtained via simulation
for training [10, 11, 12]. Unfortunately, including real US data for improved
fitting of the measurement data distribution is challenging due to lack of
access to clean target US images. Some works therefore aim at acceleration of
existing yet computationally demanding methods that provide such targets [13,
14], or resort to domain adaptation approaches [15]. A second disadvantage of
discriminative methods (especially when trained on simulations) is their
sensitivity to out-of-distribution measurement statistics, such as changes in
the noise statistics or system point-spread-function characteristics. The
alternative, using deep generative models as image priors for US image quality
improvement has to the best of our knowledge thus far not been explored.

In this paper, we propose to use deep generative models based on Normalizing
Flows (NFs) in combination with Maximum-a-Posteriori (MAP) estimation for US
image denoising and speckle suppression. NFs [16, 17, 18, 19, 20] are
generative models which have been shown to be very effective in modeling
complex image priors for sample generation, solving inverse tasks such as
image inpainting and denoising [21, 22], as well as in density estimation
[23]. We here use its ability to explicitly model density functions to fit a
high-quality target image distribution, which in turn enables performing US
image despeckling/denoising via MAP optimization under a Gaussian likelihood
model. To train a high-image-quality anatomical prior, we fit our generative
model on patches of (by nature) speckle-free MRI data, sampled from the NYU
fastMRI dataset[24, 25]. Our method thus does not rely on pairs of noisy and
speckle-free US images to train a direct estimator, but instead performs
unsupervised reconstruction via an explicit generative prior trained a-priori.

We show that our MRI-derived image priors are very effective at denoising US
images across a series of experiments, ranging from simulations with
artificially-generated speckle noise to in-vivo US images, thereby
outperforming strong traditional US speckle reduction algorithms.

Our main contributions are as follows:

  * •

We formulate the US de-speckling problem as a MAP optimization problem under a
deep generative image prior.

  * •

We propose to learn this image prior from high-quality MRI data using NFs, and
factorize the distribution to model local tissue structure (patches) rather
than high-level anatomical semantics.

  * •

We propose a method for automatic sampling of structurally informative patches
from full MRI volumes to learn improved priors for downstream MAP estimation.

The remainder of this paper is organized as follows. Section II provides the
reader with some background on ultrasound speckle and normalizing flows for
deep generative modeling. In sections III and IV the methodology, algorithms
and the performed experiments are explained, respectively. The results are
presented in section V. Finally, in sections VI and VII, the results are
discussed and conclusions derived.

##  II Background

In this section we provide the reader with some background on the ultrasound
acquisition process, leading up to the origin of speckle noise and its
characteristics (II-A). We then turn to describing normalizing flow (NF) s and
their use for flexible generative density modelling (II-B).

###  II-A Ultrasound speckle

Brightness-mode (B-mode) US aims at imaging the spatial distribution of
backscatter intensities in tissue. This is achieved by first transmitting a
short ultrasonic pulse (typically in the range of 1-10MHz), probing the
tissue, and then beamsteering an array of reiver elements to each individual
pixel (pixel-based dynamic receive beamforming) to estimate the local tissue
reflectivity (or backscatter intensity). The spatial selectivity of the
beamsteering operation (i.e. spatial filtering), impacts the resolution (main-
lobe width and pulse length) and contrast (side- and grating-lobes) of the
system. The nature of this pulse-echo experiment causes so called speckle
noise in the resulting image: the result of constructive and destructive
interference between reflections originating from many sub-wavelength spaced
scatters at pseudo-random locations. This summation of phasors with random
phase and amplitude within one resolution cell yields a strong multiplicative
noise component in the magnitude domain, that is spatially correlated by the
system resolution. A simple convolutional model of speckle noise is thus [26]:

$\bm{y}(r_{x},r_{y})=\left\\{\bm{x}(r_{x},r_{y})\cdot\bm{n}(r_{x},r_{y})\right\\}*h(r_{x},r_{y}),$
(1)

where $\bm{y}$ is the speckle-corrupted image, $\bm{x}$ is the clean
backscatter magnitude image, $\bm{n}$ is an i.i.d. multiplicative (Rayleigh)
noise component, $*$ denotes a spatial convolution, and $h$ is the
radiofrequency(RF)-modulated point spread function. Speckle-corrupted RF US
images are then envelope-detected and log-compressed before display. Empirical
models for log-compressed, envelope-detected US data often assume speckle
noise models of the form
$u(r_{x},r_{y})=v(r_{x},r_{y})+\sqrt{v(r_{x},r_{y})}\eta(r_{x},r_{y})$, with
$u(r_{x},r_{y})$ the noise-corrupted image, $v(r_{x},r_{y})$ the clean image,
and $\eta(r_{x},r_{y})$ zero-mean Gaussian noise with some standard deviation
[7]. Importantly, speckle noise is thus signal dependent.

###  II-B Deep generative models and normalizing flows

Deep generative modeling aims to learn the structure of data by fitting
complex parametric distributions to training examples. Such distributions
should satisfy $\int_{\bm{x}}p_{X}(\bm{x})=1$. For arbitrary neural networks
this does not hold and enforcing this quickly becomes intractable for high
dimensional problems. This problem is addressed by either optimizing an
evidence lower bound of the likelihood under arbitrary neural architectures
(VAEs), implicitly optimizing likelihood through adversarial training (GANs),
or more recently through the use of architectures based on invertible
bijective mappings: normalizing flows. Normalizing flows are a class of models
that learn to transform a complex high-dimensional data distribution into a
simple base distribution through a series (flow) of such invertible
transformations. This base distribution is often chosen to be a Gaussian
distribution (hence the name normalizing flow), enabling straightforward
sampling. Deep generative models based on normalizing flows have the advantage
of enabling directly likelihood estimation of data samples,
To achieve this, normalizing flows leverage the change of variables formula.
Let $X\in\mathbb{R}^{N}$ be a random variable from which we wish to model the
probability density function $p_{X}(x)$. Let $f_{\theta}(x)$ be a
differentiable, invertible transformation parameterized by $\theta$.

$Z=f_{\theta}(X)$

where $X\sim p_{X}(\bm{x})$ and $Z\sim p_{Z}(\bm{z})$. The change of variables
formula can be used to determine the probability density function of random
variable $X$ as a function of the probability density function of $Z$ as

$p_{X}(\bm{x})=p_{Z}(f_{\theta}(\bm{x}))\cdot|\det
J_{\bm{x}}(f_{\theta}(\bm{x}))|$ (2)

where $J_{x}(f_{\theta}(x))$ is the Jacobian of the transformation. The first
term, $p_{Z}(f_{\theta}(\bm{x}))$, is the probability density function of the
base distribution evaluated at $\bm{z}$. The second term, $|\det
J_{\bm{x}}(f_{\theta}(\bm{x}))|$, corrects for the change in volume caused by
the transformation. The corresponding negative-log-likelihood (NLL) is given
by

$-\log p_{X}(\bm{x})=-\log p_{Z}(\bm{z})-\log|\det
J_{\bm{x}}(f_{\theta}(\bm{x}))|,$ (3)

which can be directly optimized on a training dataset. Once the transformation
has been learned it provides a model of $p_{X}(\bm{x})$. To evaluate the
likelihood of a sample $p_{X}(\bm{x})$ one can compute $z=f_{\theta}(\bm{x})$
and evaluate the probability density function of the base distribution at
$\bm{z}$, combine it with the determinant of the Jacobian using (2). Sampling
from $p_{X}(\bm{x})$ is also straightforward. One can draw a sample $\bm{z}$
from the simple base distribution and transform it to the data distribution
using the inverse transformation

$\bm{x}=f_{\theta}^{-1}(\bm{z}).$ (4)

Figure 1: An overview of the proposed approach. We first (1) train a
normalizing flow prior on patches sampled from (by nature) speckle-free
anatomical MRI images. We then (2) use this prior for MAP inference on noisy
ultrasound measurements, by performing gradient ascent on the posterior in
latent $\bm{z}$ space, minimizing the $\ell_{2}$ distance between noisy
measurements and images lying on the learned manifold of clean anatomical
images.

##  III Methods

In this section the proposed method and experimental details are provided. The
denoising algorithm is explained in section III-A. Neural network architecture
of the prior model as well as the training data are discussed in section
III-C.

###  III-A MAP denoising using generative priors

We formulate the task of recovering the clean image $\bm{x}$ from a measured,
noisy, image $\bm{y}$ as a MAP optimization problem:

$\displaystyle\hat{\bm{x}}$
$\displaystyle=\arg\max_{\bm{x}}p_{X|Y}(\bm{x}\;|\;\bm{y})$
$\displaystyle=\arg\max_{\bm{x}}p_{Y|X}(\bm{y}\;|\;\bm{x})\cdot
p_{X}(\bm{x}).$

Under a given deep generative NF prior with normalized hidden space $\bm{z}$,
we then perform optimization in the log-$z$-domain, such that:

$\displaystyle\hat{\bm{z}}$ $\displaystyle=\arg\max_{\bm{z}}\log
p_{Z|Y}(\bm{z}\;|\;\bm{y})$ $\displaystyle=\arg\max_{\bm{z}}\left\\{\log
p_{Y|X}(\bm{y}\;|\;f_{\theta}^{-1}(\bm{z}))+\log p_{Z}(\bm{z})\right\\},$

where $f^{-1}_{\theta}(\bm{z})$ is the inverse transformation of
$f_{\theta}(\bm{z})$, i.e. the generative direction. We here approximate the
measurement process as a simple Gaussian likelihood model for the log-
compressed envelope detected US images, such that:

$\bm{y}=f^{-1}_{\theta}(\bm{z})+n,$ (5)

where $n$ follows an i.i.d Gaussian distribution. Since $p_{Z}(\bm{z})$ by
design also follows a normal distribution, we can rewrite the above
optimization problem as:

$\hat{\bm{z}}=\operatorname*{arg\,min\;}_{\bm{z}}||f^{-1}_{\theta}(\bm{z})-\bm{y}||_{2}^{2}+\lambda||\bm{z}||_{2}^{2},$
(6)

where $\lambda$ is a parameter that depends on the assumed noise variance of
$n$. Intuitively this means that the algorithm places more trust in the prior
when noise variance is high.

In this work we perform optimization using gradient descent with the Adam
optimizer, making use of the Pytorch Autograd functionality to automatically
compute gradients. The gradient step size is decayed each time the loss
reaches a plateau. We finally recover the clean image as
$\hat{\bm{x}}=f^{-1}_{\theta}(\hat{\bm{z}})$. Fig. 1 gives an overview of the
proposed approach.

###  III-B Prior factorization and MRI content-based patch sampling

The design of the prior distribution clearly $p_{X}(\bm{x})$ has a strong
impact on the effectiveness of our strategy, both in terms of noise/speckle
suppression as well as generalization to unseen anatomical instances. To
achieve the latter we factorize our prior distribution to only consider the
joint distribution of pixels within a local image patch of $14\times 14$ mm
($32\times 32$ pixels), retaining independence across larger spatial distances
and thereby per design not modeling high-level semantics that characterize
particular antomical (or even pathological) images. Note that for final image-
level denoising according to (6), we perform parallel inference (250 patches
in parallel) via gradient checkpointing [27] and then stitch the patches
together.

As mentioned, we train our deep generative prior on patches sampled from high-
quality MRI images, obtained from the NYU fast-MRI dataset (fully-sampled
images of human knees) [24, 25]. During our initial experiments we however
noticed that merely sampling patches at random from the full MRI volumes did
not yield satisfactory results. The generative model was not sufficiently
incentivised to model samples with strong inter-patch diversity and important
features such as sharp tissue boundaries and fine structural details, since
large parts of the MRI volumes did not contain such features. To promote
structural fidelity and expressiveness of the model, we favour patches with
sharp features during sampling, measured using the Tenengrad focus measure
$\phi$ [28]:

$\phi=\sum_{(i,j)\in\Omega(x,y)}(G_{x}(i,j)^{2}+G_{y}(i,j)^{2}),$ (7)

where $G_{x}(i,j)$, and $G_{y}(i,j)$ are the horizontal and vertical gradients
computed by convolving the image with the horizontal and vertical Sobel
filters respectively [29]. Fig. 2 (left) shows a random batch of selected
patches. The final train dataset consisted of $10^{6}$ such image patches.
Note that for inference on simulations with many regions without texture (e.g.
an artificial cyst phantom) which we complemented this with $50.000$
synthesized uniform (no texture) patches with intensities drawn randomly
between $0$ and $0.5$. These additional patches are included to also enable
assigning high likelihood to image regions that contain no texture.

Figure 2: A set of patches sampled from the training data (left), and
generated using the adopted normalizing flow model (right)

Figure 3: PSNR and SSIM across various Gaussian noise levels for images from
both the knee and brain NYU fast-MRI datasets. Our method outperforms the
standard non-local means filter (NLM), a competitive baseline in ultrasound
image de-speckling that we optimized to maximize PSNR via a grid search.

.

Figure 4: An illustrative example of denoising across various noise levels,
comparing clean images, (optimized) NLM, and the proposed method. The
corresponding PSNR and SSIM with respect to the clean image are displayed
below the images.

.

###  III-C Normalizing flow architecture

Our parameterization for the prior distribution is based on the Glow
architecture [18]. Glow consists of a sequence of invertible steps, that each
transform and mix the input dimensions. Each step comprises three consecutive
layers: An actnorm layer, a 1x1-convolution layer, and an affine coupling
layer. The actnorm layer is performs an affine transformation that is
initialized to result in a zero mean, and unit variance output. After
initialization the parameters become free trainable parameters.

The 1x1-convolution layer can be seen as a generalisation of a permutation
operation of the input channels. It mixes together the channels of the input
by matrix-multiplying the vector of channel values of each pixel with a
learned weight matrix $W$. $W$ is initialized as an orthogonal matrix.

Finally the data is passed through an affine coupling layer. The affine
coupling layer splits up the input $\bm{a}$ into two equal parts,
$\bm{a}_{1}$, and $\bm{a}_{2}$, along the channel dimension. It then performs
an affine transformation on $\bm{a}_{1}$, where the applied scale and shift
are a function of $\bm{a}_{2}$. $\bm{a}_{2}$ passes through unchanged. The two
parts are then concatenated into output $\bm{b}$.

$\displaystyle\bm{a}_{1},\bm{a}_{2}=\text{split}(\bm{a})$ $\displaystyle(\log
s,t)=\text{NN}(\bm{a}_{2})$ $\displaystyle s=S(\log s+2)$
$\displaystyle\bm{b}_{1}=s\cdot\bm{a}_{1}+t$
$\displaystyle\bm{b}_{2}=\bm{a}_{2}$
$\displaystyle\bm{b}=\text{concat}(\bm{b}_{2},\bm{b}_{2}),$

where $S(\cdot)$ denotes the sigmoid function. The key advantage of this
scheme is that for the coupling layer to be invertible, the function of
$\bm{a}_{2}$, $NN(\bm{a}_{2})$, does not have to be invertible. It can thus be
an arbitrarily complex function or any neural network. Here we use the same
convolutional neural network with three convolutional layers that was used in
the Glow paper[18].

The model was implemented in Pytorch and was trained with a batch size of
$800$ for $10.000$ batches using the Adam optimizer. We use gradient
checkpointing to reduce its memory footprint during training and inference
[27]. During training we moreover clip the gradients and the norm of all
gradients to improve stability.

##  IV Experiments

We perform a series of experiments that incrementally close the gap between
the assumptions made in the MAP estimator on noise statistics and real data.
To that end we test performance on anatomical MRI samples with synthetic
additive Gaussian noise, MRI-based synthetic US speckle simulations, and
finally in-silico and in-vivo US data.

Figure 5: Speckle suppression on simulated ultrasound B-mode images from MRI-
derived tissue backscatter intensity maps, comparing NLM, OBNLM (both
qualitatively optimized through a grid search), and the proposed method.

.

###  IV-A Baselines

Popular techniques for edge-preserving ultrasound speckle-suppression include
anisotropic diffusion methods [26] and stochastic iterative techniques for
outlier removal [30]. More recently, non-local means methods have demonstrated
excellent performance in speckle reduction [31, 32]. We thus compare our
method with the standard Non-Local Means (NLM) [5] and Optimized Bayesian Non-
Local Means (OBNLM) [6] methods. Non-local means methods have been shown to
outperform other ultrasound de-speckling algorithms [7, 32], and are to date
competitive [33, 32]. For fair comparison we fine-tune the settings of these
baseline algorithms to the data in the particular experiment via an exhaustive
parameter grid search. We refer the reader to the original papers for more
details.

Although recent speckle suppression based on end-to-end supervised training of
deep neural networks has been proposed [33], we here restrict ourselves to the
unsupervised setting, i.e. the setting where no pairs of speckled and speckle-
free images are available.

Figure 6: Comparison of the proposed method (B,E) against an optimized NLM
baseline (C,F) on an artificial cyst simulation (top) and a human in-vivo
carotid artery (bottom). Note that we here mix back a small amount of speckle
(20%) for improved perceptual quality (as is common practice in clinical
systems). Qualitatively, NLM produces patchy images in vivo, giving a
“posterization” effect, while the proposed method qualitatively preserves a
more natural appearance an better reveals details. This is also apparent from
the zoomed inset.

###  IV-B Denoising simulations

We first test feasibility on a set of $16$ artificially noise-corrupted MRI
images of $(64\times 64)$ pixels, extracted from a hold-out test subset of the
NYU fastMRI database, containing knees and brains. Note that our NF prior is
only trained on patches from images of knees, so the latter is useful for
assessing performance on out-of-distribution anatomy. We apply additive
Gaussian white noise having standard deviations
$\sigma_{n}\in\\{0.1,0.2,0.3,0.4\\}$. For denoising we use a
$\lambda=0.75\sigma^{2}$. Testing the algorithm for a range of noise
intensities is important since speckle noise intensity in the image domain is
known to be signal dependent (see Sec. II). We compare ourselves to the NLM-
algorithm by evaluating the PSNR and SSIM statistics [34]. We do not compare
to OBNLM here since it is not designed for additive Gaussian noise.

###  IV-C Speckle suppression simulations

To synthesize images with more realistic speckle noise from clean (MRI-based
tissue contrast) images, we follow the method by Yu et al. in [26]. Given an
oversampled and normalized MRI tissue contrast image $\bm{x}$, we generate a
synthetic image with speckle $\bm{y}$ through:

$\bm{y}(r_{x},r_{y})=\left\\{\bm{x}(r_{x},r_{y})\cdot\bm{n}(r_{x},r_{y})\right\\}*h(r_{x},r_{y}),$
(8)

with $\bm{n}(r_{x},r_{y})$ being Rayleigh-distributed noise and where

$h(r_{x},r_{y})=\sin\left[\frac{2\pi
f_{0}r_{x}}{c}\right]\exp\left[-\frac{r_{x}^{2}}{2\sigma_{r_{x}}^{2}}\right]\exp\left[-\frac{r_{y}^{2}}{2\sigma_{r_{y}}^{2}}\right],$
(9)

with $f_{0}$ being the RF pulse center frequency, $c$ the speed of sound in
the medium, and $\sigma_{r_{x}}$ and $\sigma_{r_{y}}$ determining the axial
and lateral resolution, respectively, of the virtual system. In our
experiments we use $f_{0}=7$ MHz, $c=1540$ m/s, $\sigma_{r_{x}}=0.25$ mm and
$\sigma_{r_{y}}=0.3$ mm. After this we perform envelope detection via the
Hilbert transform, log compress the images, and subsequently downscale to the
original image pixel dimensions. We again compare with the NLM algorithm, and
additionaly with the OBNLM algorithm based on parts of knee MRI images that
were withheld from the training set.

###  IV-D PICMUS and CUBDL datasets

We then evaluate performance on various US images from the PICMUS [35] and
CUBDL datasets [36, 37, 38]. From the former we include images of a cyst
phantom and a carotid artery; from the latter images of the heart. Before
processing, the beamformed images were re-sampled and interpolated to obtain
images having equal axial and lateral resolution and a pixel size that is
about same as that of the MRI patches the generative model was trained on. For
the cardiac images we performed denoising in the polar domain, i.e. before
scan conversion. Finally the images were converted to decibels and normalized
between zero and one. For denoising we used a fixed $\lambda=10$.

The NLM algorithm was applied with $h=0.8\sigma$, a patch size of $10$, and a
patch distance of $10$. The OBNLM algorithm was applied with $h=\sigma$, a
patch size of $3$, and a similarity window of size $10$.

##  V Results

###  V-A Denoising simulations

We compare performance of the proposed method against our baseline at various
noise levels in Fig. 3. The proposed method outperforms the NLM algorithm
across all evaluated noise levels. Noteably, it also generalizes well to
images from the brain dataset of the NYU fast-MRI challenge. This supports our
hypothesis that by choosing a prior factorization that only models joint pixel
dependencies across a relatively small local neighborhood of $1.4\times 1.4$
cm, the model generalizes across anatomies.

Fig. 4 qualitatively shows how performance is impacted by the overall noise
level for the two methods in a typical example. Across the noise levels, the
proposed method better captures the details and structure of the original
image. Although we notice that our method produces slightly more blurred
features compared to the clean input, it visually outperforms the NLM
algorithm.

###  V-B Speckle suppression simulations

Fig. 5 shows the results on two speckle simulations generated using intensity
maps derived from MRI images of (segments of) human knees. Comparing the
resulting denoised images by NLM, OBNLM and the proposed methods with the
intensity maps (and its features), we observe that both NLM and the proposed
method provide strong denoising. OBNLM preserves more details than NLM, but
also retains more noise. The proposed method preserves details well, without
sacrificing denoising capabilities.

###  V-C PICMUS and CUBDL datasets

The results from experiments on the PICMUS dataset are shown in Fig. 6. The
first observation is that the proposed method translates well to the real in-
vivo ultrasound data, showing the feasibility of using MRI-based image priors
for ultrasound denoising and speckle suppression. While the MRI database did
not include perfect circle-shaped low-intensity regions, the proposed method
does not distort the shapes and remains true to the geometry, showing the
strong generalization of our prior. Moreover, comparing to the NLM method, we
observe that it produces more natural and less patchy denoising, as is common
with NLM methods. Although this effect is not very noticeable on the cyst
simulations since these artificial constellations indeed only consist of two
levels of scattering intensity, it is clear on the more intricate in-vivo data
of a human carotid artery. From the zoomed areas we can also appreciate better
preservation of detail with the proposed approach.

Fig. 7 offers a side by side comparison including OBNLM. OBNLM has been
designed to cope with the intensity-dependent noise statistics, and produces
less patchy results. Yet, OBNLM qualitatively also remains much more noisy.

We show some additional denoising results on cardiac images from the CUBDL
database in Fig. 8. This indicates how the method generalizes across domains,
and effectively leverages the image priors derived from MRI for various
applications.

Figure 7: Comparison among NLM (B), OBNLM (C), and the proposed method (D) in-
vivo. Qualitatively, OBNLM is less patchy than NLM, but also more noisy. The
proposed method retains a more natural appearance, yet with strong noise
suppression.

##  VI Discussion

In this paper, we demonstrated that patch-based iterative denoising with a
deep generative normalizing-flow image prior can produce sharp and detailed
images from ultrasound acquisitions distorted by severe noise and speckle. Our
experiments show that the proposed method numerically outperforms a fine-tuned
NLM algorithm in terms of both SSIM, and PSNR across a range of noise powers.
The method qualitatively outperforms both NLM and OBNLM on ultrasound speckle
simulations and real in-vivo data from the PICMUS challenge. Notably, it
produces less patchy and perceptually more natural outputs compared to the
other algorithms, which yield an almost segmented output, reminiscent of
“posterization”.

Although the results are promising, the required processing time used for
iterative MAP inference is currently too high for real-time application.
Although not the focus of this work, one avenue to explore is the use of deep
unfolding [39], in which the iterations of the algorithm are unfolded as a
feed-forward deep network [22]. In other applications this has yielded major
speedups (as high as a factor 100) [40, 41], potentially enabling real-time
inference. Another means of acceleration would be to exploit the temporal
structure/persistence of ultrasound sequences, initializing the MAP solver for
the next frame with the with the solution for the previous frame. This way the
algorithm “tracks” the statistics. We leave this for future work.

To further boost performance, one could explore improvements in the adopted
data likelihood model. While we here assume white Gaussian error statistics,
this is in practice naive, since speckle noise is (1) spatially correlated,
and (2) signal-dependent. Formulating a more accurate model could lead in the
future lead to additional performance gains. Note that such models can also be
formulated to jointly perform denoising/despeckling and deconvolution, by
including a (convolutional) point-spread function model in Eqn. 6.

More broadly, we foresee application of the proposed MRI-derived image priors
to ultrasound compressed sensing settings such as sub-Nyquist acquisitions
[42] and sparse arrays [43], e.g. overcoming the challenges of analytically
formulating an adequate sparsifying basis with learned priors.

Finally, a more in-depth comparison with a large set of real ultrasound images
is required, including clinical assessment, which we leave to future work.

##  VII Conclusion

We present a new method for ultrasound image denoising and speckle suppression
that uses MRI-based structural priors learned with deep generative models to
perform maximum-a-posteriori estimation of clean images. To learn intricate
image priors from MRI data that generalize to unseen anatomies, we devise an
invertible normalizing flow architecture that models a factorized density
function based on local image patches. Notably, the proposed deep learning
method thus requires no supervision for its inference (no pairs of clean and
noisy data), and outperforms strong baselines based on non-local-means
quantitatively and qualitatively.

Figure 8: Additional denoising results on cardiac images from the CUBDL
dataset [38], demonstrating how the method generalizes across various
anatomies.

## References

  * [1] J.-S. Lee, “Digital image enhancement and noise filtering by use of local statistics,” 1977.
  * [2] V. S. Frost, J. A. Stiles, K. S. Shanmugan, and J. C. Holtzman, “A model for radar images and its application to adaptive digital filtering of multiplicative noise,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. PAMI-4, pp. 157–166, 1982.
  * [3] D. T. Kuan, A. A. Sawchuk, T. C. Strand, and P. Chavel, “Adaptive noise smoothing filter for images with signal-dependent noise.,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. PAMI-7, pp. 165–177, 1985.
  * [4] Y. Yu and S. T. Acton, “Speckle reducing anisotropic diffusion,” IEEE Transactions on Image Processing, vol. 11, pp. 1260–1270, 11 2002.
  * [5] A. Buades, B. Coll, and J. M. Morel, “A non-local algorithm for image denoising,” vol. II, pp. 60–65, IEEE Computer Society, 2005.
  * [6] Y. Zhan, M. Ding, L. Wu, and X. Zhang, “Nonlocal means method using weight refining for despeckling of ultrasound images,” Signal Processing, vol. 103, pp. 201–213, 2014.
  * [7] P. Coupé, P. Hellier, C. Kervrann, and C. Barillot, “Bayesian non local means-based speckle filtering,” in 2008 5th IEEE International Symposium on Biomedical Imaging: From Nano to Macro, pp. 1291–1294, IEEE, 2008\.
  * [8] D. Li, W. Yu, K. Wang, D. Jiang, and Q. Jin, “Speckle noise removal based on structural convolutional neural networks with feature fusion for medical image,” Signal Processing: Image Communication, vol. 99, p. 116500, 11 2021\.
  * [9] O. Karaoğlu, H. Şakir Bilge, and İhsan Uluer, “Removal of speckle noises from ultrasound images using five different deep learning networks,” Engineering Science and Technology, an International Journal, 2021.
  * [10] D. Hyun, L. L. Brickson, K. T. Looby, and J. J. Dahl, “Beamforming and speckle reduction using neural networks,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 66, pp. 898–910, 5 2019.
  * [11] S. Vedula, O. Senouf, A. M. Bronstein, O. V. Michailovich, and M. Zibulevsky, “Towards ct-quality ultrasound imaging using deep learning,” arXiv preprint arXiv:1710.06304, 2017.
  * [12] X. Feng, Q. Huang, and X. Li, “Ultrasound image de-speckling by a hybrid deep network with transferred filtering and structural prior,” Neurocomputing, vol. 414, pp. 346–355, 2020.
  * [13] B. Luijten, R. Cohen, F. J. de Bruijn, H. A. Schmeitz, M. Mischi, Y. C. Eldar, and R. J. van Sloun, “Adaptive ultrasound beamforming using deep learning,” IEEE Transactions on Medical Imaging, vol. 39, no. 12, pp. 3967–3978, 2020\.
  * [14] F. Dietrichson, E. Smistad, A. Ostvik, and L. Lovstakken, “Ultrasound speckle reduction using generative adversial networks,” in 2018 IEEE International Ultrasonics Symposium (IUS), pp. 1–4, IEEE, 2018.
  * [15] J. Tierney, A. Luchies, C. Khan, B. Byram, and M. Berger, “Domain adaptation for ultrasound beamforming,” in International Conference on Medical Image Computing and Computer-Assisted Intervention, pp. 410–420, Springer, 2020\.
  * [16] L. Dinh, J. Sohl-Dickstein, and S. Bengio, “Density estimation using real nvp,” 5 2016.
  * [17] L. Dinh, D. Krueger, and Y. Bengio, “Nice: Non-linear independent components estimation,” 10 2014.
  * [18] D. P. Kingma and P. Dhariwal, “Glow: Generative flow with invertible 1x1 convolutions,” 7 2018.
  * [19] J. Ho, X. Chen, A. Srinivas, Y. Duan, and P. Abbeel, “Flow++: Improving flow-based generative models with variational dequantization and architecture design,” 2 2019.
  * [20] I. Kobyzev, S. J. D. Prince, and M. A. Brubaker, “Normalizing flows: An introduction and review of current methods,” 8 2019.
  * [21] M. Asim, M. Daniels, O. Leong, A. Ahmed, and P. Hand, “Invertible generative models for inverse problems: mitigating representation error and dataset bias,” in International Conference on Machine Learning, pp. 399–409, PMLR, 2020.
  * [22] X. Wei, H. van Gorp, L. G. Carabarin, D. Freedman, Y. Eldar, and R. van Sloun, “Deep unfolding with normalizing flow priors for inverse problems,” arXiv preprint arXiv:2107.02848, 2021.
  * [23] G. Papamakarios, T. Pavlakou, and I. Murray, “Masked autoregressive flow for density estimation,” 5 2017.
  * [24] J. Zbontar, F. Knoll, A. Sriram, T. Murrell, Z. Huang, M. J. Muckley, A. Defazio, R. Stern, P. Johnson, M. Bruno, M. Parente, K. J. Geras, J. Katsnelson, H. Chandarana, Z. Zhang, M. Drozdzal, A. Romero, M. Rabbat, P. Vincent, N. Yakubova, J. Pinkerton, D. Wang, E. Owens, C. L. Zitnick, M. P. Recht, D. K. Sodickson, and Y. W. Lui, “fastmri: An open dataset and benchmarks for accelerated mri,” 11 2018.
  * [25] F. Knoll, J. Zbontar, A. Sriram, M. J. Muckley, M. Bruno, A. Defazio, M. Parente, K. J. Geras, J. Katsnelson, H. Chandarana, Z. Zhang, M. Drozdzalv, A. Romero, M. Rabbat, P. Vincent, J. Pinkerton, D. Wang, N. Yakubova, E. Owens, C. L. Zitnick, M. P. Recht, D. K. Sodickson, and Y. W. Lui, “fastmri: A publicly available raw k-space and dicom dataset of knee images for accelerated mr image reconstruction using machine learning,” Radiology: Artificial Intelligence, vol. 2, p. e190007, 1 2020.
  * [26] Y. Yu and S. T. Acton, “Speckle reducing anisotropic diffusion,” IEEE Transactions on image processing, vol. 11, no. 11, pp. 1260–1270, 2002.
  * [27] T. Chen, B. Xu, C. Zhang, and C. Guestrin, “Training deep nets with sublinear memory cost,” 4 2016.
  * [28] S. Pertuz, D. Puig, and M. A. Garcia, “Analysis of focus measure operators for shape-from-focus,” Pattern Recognition, vol. 46, pp. 1415–1432, 5 2013\.
  * [29] I. Sobel, “An isotropic 3x3 image gradient operator,” Presentation at Stanford A.I. Project 1968, 11 2014.
  * [30] P. C. Tay, C. D. Garson, S. T. Acton, and J. A. Hossack, “Ultrasound despeckling for contrast enhancement,” IEEE Transactions on Image Processing, vol. 19, no. 7, pp. 1847–1860, 2010.
  * [31] L. H. Breivik, S. R. Snare, E. N. Steen, and A. H. S. Solberg, “Real-time nonlocal means-based despeckling,” IEEE transactions on ultrasonics, ferroelectrics, and frequency control, vol. 64, no. 6, pp. 959–977, 2017.
  * [32] M. Ambrosanio, B. Kanoun, and F. Baselice, “wksr-nlm: An ultrasound despeckling filter based on patch ratio and statistical similarity,” IEEE Access, vol. 8, pp. 150773–150783, 2020.
  * [33] D. Hyun, L. L. Brickson, K. T. Looby, and J. J. Dahl, “Beamforming and speckle reduction using neural networks,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 66, pp. 898–910, 5 2019.
  * [34] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” IEEE Transactions on Image Processing, vol. 13, pp. 600–612, 4 2004.
  * [35] H. Liebgott, A. Rodriguez-Molares, F. Cervenansky, J. A. Jensen, and O. Bernard, “Plane-wave imaging challenge in medical ultrasound,” vol. 2016-November, IEEE Computer Society, 11 2016.
  * [36] D. Hyun, A. Wiacek, S. Goudarzi, S. Rothlubbers, A. Asif, K. Eickel, Y. C. Eldar, J. Huang, M. Mischi, H. Rivaz, D. Sinden, R. J. V. Sloun, H. Strohm, and M. A. L. Bell, “Deep learning for ultrasound image formation: Cubdl evaluation framework and open datasets,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 12 2021.
  * [37] O. M. H. Rindal, S. Aakhus, S. Holm, and A. Austeng, “Hypothesis of improved visualization of microstructures in the interventricular septum with ultrasound and adaptive beamforming,” Ultrasound in Medicine & Biology, vol. 43, no. 10, pp. 2494–2499, 2017.
  * [38] M. A. L. Bell, J. Huang, A. Wiacek, P. Gong, S. Chen, A. Ramalli, P. Tortoli, B. Luijten, M. Mischi, O. M. H. Rindal, V. Perrot, H. Liebgott, X. Zhang, J. Luo, E. Oluyemi, and E. Ambinder, “Challenge on Ultrasound Beamforming with Deep Learning (CUBDL) Datasets.”
  * [39] V. Monga, Y. Li, and Y. C. Eldar, “Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing,” IEEE Signal Processing Magazine, vol. 38, no. 2, pp. 18–44, 2021.
  * [40] O. Solomon, R. Cohen, Y. Zhang, Y. Yang, Q. He, J. Luo, R. J. van Sloun, and Y. C. Eldar, “Deep unfolded robust pca with application to clutter suppression in ultrasound,” IEEE transactions on medical imaging, vol. 39, no. 4, pp. 1051–1063, 2019.
  * [41] R. J. Van Sloun, R. Cohen, and Y. C. Eldar, “Deep learning in ultrasound imaging,” Proceedings of the IEEE, vol. 108, no. 1, pp. 11–29, 2019.
  * [42] T. Chernyakova and Y. C. Eldar, “Fourier-domain beamforming: the path to compressed ultrasound imaging,” IEEE transactions on ultrasonics, ferroelectrics, and frequency control, vol. 61, no. 8, pp. 1252–1267, 2014.
  * [43] I. A. Huijben, B. S. Veeling, K. Janse, M. Mischi, and R. J. van Sloun, “Learning sub-sampling and signal recovery with applications in ultrasound imaging,” IEEE Transactions on Medical Imaging, vol. 39, no. 12, pp. 3955–3966, 2020.