# Speckle2Speckle: Unsupervised Learning of Ultrasound Speckle Filtering
Without Clean Data

Rüdiger Göbl<EMAIL_ADDRESS>
Computer Aided Medical Procedures, Technical University Munich, Munich,
Germany
OneProjects Design Innovation GmbH, Munich, Germany Christoph Hennersperger
<EMAIL_ADDRESS>
Computer Aided Medical Procedures, Technical University Munich, Munich,
Germany
OneProjects Design Innovation GmbH, Munich, Germany Nassir Navab
<EMAIL_ADDRESS>
Computer Aided Medical Procedures, Technical University Munich, Munich,
Germany

###### Abstract

In ultrasound imaging the appearance of homogeneous regions of tissue is
subject to speckle, which for certain applications can make the detection of
tissue irregularities difficult. To cope with this, it is common practice to
apply speckle reduction filters to the images. Most conventional filtering
techniques are fairly hand-crafted and often need to be finely tuned to the
present hardware, imaging scheme and application. Learning based techniques on
the other hand suffer from the need for a target image for training (in case
of fully supervised techniques) or require narrow, complex physics-based
models of the speckle appearance that might not apply in all cases. With this
work we propose a deep-learning based method for speckle removal without these
limitations. To enable this, we make use of realistic ultrasound simulation
techniques that allow for instantiation of several independent speckle
realizations that represent the exact same tissue, thus allowing for the
application of image reconstruction techniques that work with pairs of
differently corrupted data. Compared to two other state-of-the-art approaches
(non-local means and the Optimized Bayesian non-local means filter) our method
performs favorably in qualitative comparisons and quantitative evaluation,
despite being trained on simulations alone, and is several orders of magnitude
faster. Our code is available at https://github.com/goeblr/Speckle2Speckle.

Keywords: CNN, Deep Learning, Speckle filtering, Ultrasound, Unsupervised

##  1 Introduction

Ultrasound imaging is a valuable tool to support non-invasive diagnosis and
for use during interventions, as it is portable, affordable and does not
subject the patient or staff to ionizing radiation. It is not, however,
without problems. One of the major differences between CT and MRI on the one
hand and ultrasound on the other hand is the presence of speckle, a grainy,
noiselike image degradation that even affects tissue regions that are
homogeneous in their microstructure such as liver or thyroid tissue (e.g. see
1(a)). This speckle can make interpretation of the images difficult—especially
for users with limited experience working with ultrasound. Consequently,
most—if not all—ultrasound systems apply speckle reduction techniques after
image formation, where the amount of ”graininess” that is retained usually can
be controlled by the user through a linear combination or blending of filtered
and unfiltered data, depending on the clinical task and user preference.
Figure 1 shows a thyroid image before and after filtering with our method
compared to a state-of-the-art method and outlines a way to preserve a
controlled amount of speckle.

(a) Unfiltered

(b) Filtered: OBNLM

(c) Filtered: Ours

(d) Blending

Figure 1: Comparison of image appearance of a thyroid. Unfiltered (a),
despeckling with a state-of-the-art method (OBNLM) (b), despeckling with our
method (c) and a linear combination of (a) and (c) in (d).

The speckle patterns that emerge in ultrasound images are directly influenced
by the location of sub-resolution scatterers in the medium. Being (spaced)
sub-resolution, the reflections caused by the individual scatterers cannot be
separated by the ultrasound imaging system. The result is interference between
the echoes of different amplitude and phase, leading to an apparent
echogenicity that changes depending on the spatial configuration of the
scatterers within each resolution cell and the orientation of transmitting and
receiving elements. In addition to speckle, there is also electronic noise
that degrades ultrasound images, which is not the focus of this work. Speckle,
while sometimes referred to as noise, is a deterministic phenomenon. In fact,
if it were possible to place an ultrasound transducer at the exact same
location when repeating a scan, the resulting speckle pattern would be exactly
the same.

##  2 Related Works and Contribution

Conventional despeckling approaches such as speckle-reducing anisotropic
diffusion (SRAD) (Yu and Acton, 2002) and the Optimized Bayesian non-local
means filter (OBNLM) (Coupé et al., 2009)—itself building upon non-local means
(NLM) (Buades et al., 2005)—have been studied extensively. We refer the
interested reader to the review from Mohd Sagheer and George (2020).
Especially the latter method OBNLM exhibits excellent despeckling performance,
but requires application-specific tuning and significant processing time.

Cammarasana et al. (2021) alleviate both these concerns by using a regression
CNN that is trained against images despeckled with a tuned, computationally
expensive conventional method, thereby allowing the CNN to replicate its
results while lowering the execution time.

In other areas different to ultrasound imaging, there exists a body of work on
the task of image denoising without requiring perfectly matched pairs of
corrupted and clean data. In Noise2Noise (Lehtinen et al., 2018), a CNN learns
to restore images from a number of different corruptions, by using a corrupted
image as input and a second image of the same ”scene” with a different
realization of the corruption as training target. The authors show that this
can reach the performance of methods training against clean images. Noise2Void
(Krull et al., 2019) takes this even further by requiring only unpaired
corrupted images, it requires however, strong assumptions on the corruption,
namely zero-mean and per-pixel independent noise.

Given the difficulties of acquiring clean images in many medical imaging
modalities, these general approaches have been further specialized. In optical
coherence imaging, a few well-chosen modifications to the acquisition setup
can allow for the acquisition of two independent speckle instances with
otherwise unchanged scanning setup (Yin et al., 2020), effectively allowing
for direct use of Noise2Noise. However, this method cannot be directly applied
to ultrasound imaging, as these modifications are not possible in ultrasound
scanning setups.

In this work, we build on the concept of exploiting two independent
_corrupted_ observations with the same scanning setup and geometry, and
transfer this to a generalizable method for ultrasound imaging. In this way,
we present a technique for learning a CNN-based ultrasound despeckle filter
only from images containing speckle. To do so, we apply the concept of
Noise2Noise (Lehtinen et al., 2018), combined with a scheme for ultrasound
simulation for data generation, and modifications to the training loss to
enhance the image appearance in a method we call Speckle2Speckle. By
simulating different observations of the same general anatomical geometry but
with different scatterer locations, we create images with uncorrelated speckle
patterns. In our qualitative comparison and quantitative evaluation with
state-of-the-art approaches, we show that the method can remove speckle from
real ultrasound data efficiently while not relying on manual data annotation
or acquisitions, while the application of the filter is several orders of
magnitude faster.

##  3 Method

In this section we first summarize the original formulation of Noise2Noise,
followed by a description how we create simulated ultrasound images to be used
with a denoising approach like this, and the modifications that were required
to account for the peculiarities of ultrasound. Figure 2 shows a comparison
between the concepts for conventional (clean-target) (Cammarasana et al.,
2021), Noise2Noise and Speckle2Speckle filter learning.

(a) Clean target

(b) Noise2Noise

(c) Ours

Figure 2: Overview of training a model with corrupted data input to clean data
output (a), Noise2Noise training with both corrupted data inputs and outputs
of different realizations (b), and Speckle2Speckle training with corrupted
realizations of ultrasound data simulated from the same tissue model (c).

###  3.1 Noise2Noise

Noise2Noise (Lehtinen et al., 2018) is a deep-learning based technique for
general image restoration that does not require clean, uncorrupted data for
training. In a conventional restoration approach, one would use corrupted
images as input and clean images of the same scene as the training target. In
their work however, Lehtinen et al. derive, from a statistical perspective,
that a network trained with infinite samples can learn to estimate the
expectation of the target samples. This allows the input images as well as the
target images to be corrupted with noise, like image noise observed with
digital image sensors. The authors show that the networks trained in this
manner perform on-par with networks trained against clean data.

These relaxed requirements on the training targets allow for the application
of this technique to situations where the acquisition of clean targets is
either costly or impossible all-together. In their work, Noise2Noise using a
modified U-Net (following Ronneberger et al. (2015), see Figure 3) is applied
to a number of different image corruption scenarios ranging from photographs
subjected to Gaussian-, Poisson- and Bernoulli-noise or text-overlays, Monte-
Carlo rendered images, and even MRI acquisitions within the randomly sampled
k-space.

Figure 3: Network architecture of Noise2Noise and Speckle2Speckle. Numbers
left of each box show the sizes of the feature maps, exemplified with a 128 by
128 pixel input. The number of channels per feature map is given by the
numbers above the boxes. Different arrows represent different operations. The
rightmost convolution (marked with $*$) has linear activation.

###  3.2 Speckle2Speckle: Data

Ultrasound imaging is a perfect example for a situation where a noise-free
target image is difficult to obtain: since diffuse scattering is present for
any ultrasound acquisition, speckle patterns remain visible at any time for
acquired echo data. Thus, one can argue that in ultrasound imaging, there is
no image without speckle. The strong correlation between the location and
scattering behaviour of the tissue inhomogeneities causing the echos and the
locations of extrema in the observed reflections (i.e., the speckle pattern),
however, makes the acquisition of multiple images with independent corruption
(of the desired echogenicity measurement) difficult for real acquisitions.

However, we can instead turn to realistic ultrasound simulations, where we can
exploit the specific control over individual scatterers to generate an
arbitrary amount of corrupted datasets for a given setup. The training data
for Speckle2Speckle is generated with the simulation software Field II (Jensen
and Svendsen, 1992; Jensen, 1996). Images are generated in pairs, using
independent in-silico scatterer phantom datasets (collections of scatterers
with associated scattering strengths). These in-silico scatterer phantom
datasets are derived from the same geometric phantom, which is a collection of
3D shapes each of which with a number of parameters describing their
appearance in ultrasound imaging and scatterer parameters:

  * •

scatterer density

  * •

scatterer amplitude distribution

  * •

presence or absence of an interface

  * •

interface scatterer density

  * •

interface scatterer amplitude distribution

Based on the these properties, this approach allows for the generation of
scatterer phantom datasets with independent scatterer locations and thus
ideally suited for approaches such as Noise2Noise. Scatterer locations can be
generated randomly, but since the ”macro-scale” geometry (e.g. areas of higher
scattering) is the same, the ultrasound images created from the different
scatterer phantoms look structurally identical and even the speckle
distributions are identical. The main difference, however, is that the speckle
pattern ”instance” is different, for which an example is shown in 4(a) and
4(b).

(a) Instance 1

(b) Instance 2

(c) Average

(d) Interface map

Figure 4: Comparison of the different image types created through simulation:
Ultrasound images of two different scatterer instances from the same phantom
geometry (a) and (b), average of nine ultrasound images from the same geometry
(c), interface map derived from the phantom geometry (d).

While the training is performed with pairs of scatterer instances, we use a
different target for the testing and validation. In order to facilitate the
evaluation of speckle removal in the absence of noise-free images, we compare
the output of our method with an image that is composed by averaging multiple
ultrasound images of the same phantom geometry but different scatterer
instances. For the training and validation sets, we generated 10 instances per
phantom, one of which is used (randomly) as input, the other nine are
averaged. 4(c) shows the average of nine images corresponding to the same
geometry as in 4(a) and 4(b).

In addition to the ultrasound images, we also derive an interface map directly
from the geometric phantom. Being of the same size and spatial origin as the
simulated ultrasound images, it contains whether an interface of a shape was
present at each location. 4(d) shows the interface map associated with the
geometry used for the simulations of the ultrasound images in Figure 4. This
is later utilized to adapt the training loss to the challenges of ultrasound
imaging in order to achieve enhancements of prominent interfaces in the
output.

Following this general scheme, a dataset was generated with a training set
consisting of 1000 image pairs while the validation and test sets each contain
100 phantom geometries with 10 speckle instances per phantom geometry,
resulting in 4000 simulated images in total. Fundamentally, the size of the
dataset was limited by the simulation time, where each image required roughly
between one and two days on a single CPU core to simulate (depending on the
concrete number of scatterers and the CPU performance), although
parallelization between images was trivial.

In order to limit the simulation run-time, the scatterer densities values were
determined as high as required such that further increase did not cause a
change in the speckle patterns. The probabilities of geometries being anechoic
and the scatterer amplitude distributions were determined arbitrarily to
ensure visibility of the different regions. The specific parametrization is
listed in Table 1.

Region | Parameter | Value
---|---|---
Background | probability anechoic | 0.4
| scatterer density | $0.333\text{\times}{10}^{10}\text{\,}\mathrm{m}^{-3}$
| scatterer amplitude |  normal distributed $\mathcal{N}(0,\sigma_{\text{bg}}^{2})$; $\sigma_{\text{bg}}$ drawn per phantom from $\mathcal{N}(1,0.5^{2})$
Inclusion | random objects | spheroids & cuboids, uniformly distributed
| number per phantom | 100
| extents | drawn per axis uniformly from $[1,5]\text{\,}\mathrm{m}\mathrm{m}$
| probability anechoic | 0.4
| scatterer density | $0.333\text{\times}{10}^{10}\text{\,}\mathrm{m}^{-3}$
| scatterer amplitude |  normal distributed $\mathcal{N}(0,\sigma_{\text{inc}}^{2})$; $\sigma_{\text{inc}}$ drawn per inclusion from $\mathcal{N}(4,2^{2})$
Inclusion Interface | probability interface | 0.5
| scatterer density | $5\text{\times}{10}^{6}\text{\,}\mathrm{m}^{-2}$
| scatterer amplitude |  constant per inclusion, drawn per inclusion from $\mathcal{N}(14,2^{2})$
Table 1: Simulated phantom parameters and values

The simulated imaging protocol used focused imaging and it followed the same
restrictions as an actual focused scan using a CPLA12875 probe with a cQuest
Cicada (both Cephasonics, CA, USA) would be subject to. Table 2 contains the
detailed acquisition parameters.

Parameter | Value
---|---
transducer | CPLA12875
number elements | 128
element pitch | $0.30\text{\,}\mathrm{m}\mathrm{m}$
elevation focus | None
transmit frequency | $7\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$
transmit pulse | bipolar
imaging depth | $60\text{\,}\mathrm{m}\mathrm{m}$
focus depth | $30\text{\,}\mathrm{m}\mathrm{m}$
steering | $0\text{\,}\mathrm{\SIUnitSymbolDegree}$
number scanlines | 128
TX aperture size | 32
RX aperture size | 64
Table 2: Simulated acquisition parameters and values

The output of the simulations was ultrasound radio-frequency channel data,
which was subsequently processed to ultrasound images using SUPRA (Göbl et
al., 2018) with parameters as in Table 3. After beamforming, the resulting
images were converted to PNGs and interface indicator maps for inclusions with
an interface were created by examining the geometry at every pixel location
for the presence of an interface.

Parameter | Value
---|---
beamforming | DAS with hamming window
dynamic range | 70 dB
image resolution | 0.075 mm
image width |  $37.6\text{\,}\mathrm{m}\mathrm{m}$ $\widehat{=}$ 502 pixels
image height |  $60.0\text{\,}\mathrm{m}\mathrm{m}$ $\widehat{=}$ 801 pixels
Table 3: Image reconstruction parameters and values

###  3.3 Speckle2Speckle: Learning

The network architecture employed in Speckle2Speckle follows Noise2Noise (see
subsection 3.1).

The training loss, however, was adjusted for the peculiarities of ultrasound
imaging. While the original Noise2Noise approach used different loss data
terms depending on the noise, we extended the data term $L_{2}$ with a term to
enhance the appearance of interfaces in the images, as interfaces retain key
information for interpreting ultrasound imaging data.

The interface term is based on the known locations of interfaces given by the
interface indicator map $I_{\text{i}}$. We chose
$I_{\text{i}G}=g_{\text{i}}*I_{\text{i}}$ as spatial weight for the interface
loss, where $g_{\text{i}}$ is a 2D Gaussian with standard deviation
$\sigma_{i}$. This serves to i) extend the range of the interface effects from
the potentially sharp peaks in the interface map, as well as ii) achieve a
smooth transition between the effects of the data term and the interface term.
$\tilde{I}_{\text{i}G}=1-I_{\text{i}G}$ is the inverse weight. With that, we
formulate the overall loss $L_{D}$ as

$\displaystyle
L_{D}(I_{\text{o}},I_{\text{t}},I_{\text{i}})=L_{2}(I_{\text{o}}\odot\tilde{I}_{\text{i}G},I_{\text{t}}\odot\tilde{I}_{\text{i}G})+\lambda
S_{D}(I_{\text{o}},I_{\text{t}},I_{\text{i}}),$ (1) $\displaystyle
S_{D}(I_{\text{o}},I_{\text{t}},I_{\text{i}})=\frac{1}{|\Omega|}\sum_{x\in\Omega}\left((I_{\text{o}}*g_{\text{psf}}-I_{\text{t}})\odot
I_{\text{i}G}\right)^{2},$ (2)

where $I_{\text{o}}$ is the network output, $I_{\text{t}}$ the target image,
$\lambda$ the weighting of the interface term, $\odot$ denotes the Hadamard
product (element-wise multiplication) and $g_{\text{psf}}$ a filter kernel of
the point spread function (PSF). The notion is: In the target image, the
interfaces are smoothed by the PSF by the process of imaging. Whereas in the
output image the goal is to reduce the blurring, i.e., increase the sharpness.
Consequently, the term $I_{\text{o}}*g_{\text{psf}}-I_{\text{t}}$ penalizes
deviations in $I_{\text{o}}$ from sharp peaks. The PSF in this case could be
the result of a per-image estimation, learned for one specific system,
manually determined for one system or treated as hyper-parameter. We assumed
the PSF to be an anisotropic Gaussian, with
$\sigma_{\text{psf}}^{x}=\text{FWHM}_{\text{lat}}/2$ and
$\sigma_{\text{psf}}^{y}=\text{FWHM}_{\text{ax}}/2$ determined experimentally
from a simulated point reflector, where we measured the full width half
maximum (FWHM) in lateral and axial directions. The hyper-parameters and their
values are listed in Table 4.

Parameter | Value
---|---
$\lambda$ | $500$
learning rate | $3\text{\times}{10}^{-5}$
epochs | $1000$
$\sigma_{\text{i}}$ (interface) |  $0.375\text{\,}\mathrm{m}\mathrm{m}$ $\widehat{=}$ 5 pixels
$\sigma_{\text{psf}}^{y}$ (axial) |  $0.075\text{\,}\mathrm{m}\mathrm{m}$ $\widehat{=}$ 1 pixel
$\sigma_{\text{psf}}^{x}$ (lateral) |  $0.525\text{\,}\mathrm{m}\mathrm{m}$ $\widehat{=}$ 7 pixels
Table 4: Training hyper-parameters and their values

##  4 Evaluation

In order to achieve a fair comparison we evaluate our approach to a group of
well-established and state-of-the-art speckle filters, and provide qualitative
comparisons and quantitative evaluations for the simulated, phantom, and in-
vivo datasets in the following. To this end, we compare the results of our
method against SRAD111implementation from Virginia Image and Video Analysis,
School of Engineering and Applied Science, University of Virginia (2006) (Yu
and Acton, 2002), median filter222implementation from MATLAB R2020b, bilateral
filter††footnotemark:  (Tomasi and Manduchi, 1998), non-local
means††footnotemark:  (NLM) (Buades et al., 2005), and optimized Bayesian non-
local means333implementation from Coupé (2011) (OBNLM) (Coupé et al., 2009).
Table 5 contains the parametrization used for the mentioned methods.

Method | Parameter | Value
---|---|---
SRAD | number of iterations | 200
| lambda | 0.1
Median | window size | 15
Bilateral | degree of smoothing | 0.05
| spatial sigma | 5
NLM | degree of smoothing | 0.075
| search windows size | 101
| comparison window size | 21
OBNLM | search area size | 101
| patch size | 45
| degree of smoothing | 1.05
Table 5: Comparison method parameters and their values after manual tuning on
the image in 7(a).

###  4.1 Qualitative comparison

(a) Original

(b) Average image

(c) SRAD

(d) Median

(e) Bilateral

(f) NLM

(g) OBNLM

(h) Ours

Figure 5:  Qualitative comparison on a simulated image. The original speckled
image (a), the average of nine different scatterer instantiations (b), output
of SRAD (c), median filter (d), bilateral filtering (e), NLM (f), OBNLM (g)
and ours (h).

_Simulated data._ In the simulated image from the validation set (see
subsection 3.2) shown in Figure 5, we can visually compare the result of
applying the different methods directly with the average image of nine speckle
instances. It is clear that SRAD (c), median filtering (d) and the bilateral
filter (e) mostly change the appearance of the speckle, but don not fully
remove it. NLM (f), OBNLM (g) and our method (h) all remove the speckle very
well. In the outlined area in the top of the image, the less dense speckle is
not fully removed by OBNLM. Another noteworthy region is the rectangular
object in the top left that has a bright interface in the input as well as the
average image. NLM and OBNLM do not retain this feature with high contrast,
whereas our method shows contrast close to the average image. The elliptic
object in the bottom right appears in our method not as clearly as it is
visible in the average image, we believe this is due to comparatively low
brightness of interface around it. The other methods, however, also struggle
to highlight the object.

(a) Original

(b) SRAD

(c) Median

(d) Bilateral

(e) NLM

(f) OBNLM

(g) Ours

Figure 6:  Qualitative comparison on an experimental phantom image. The
original speckled image (a), the output of SRAD (b), median filter (c),
bilateral filtering (d), NLM (e), OBNLM (f) and ours (g). The marked regions
in the inclusion and background are used for quantitative evaluation.

_Phantom data._ The experimental image in Figure 6 shows an ultrasound quality
assurance phantom (Model 040GSE, CIRS Inc., Norfolk, VA, USA). The image was
acquired with a CPLA12875 probe and a cQuest Cicada (both Cephasonics, CA,
USA), with equivalent imaging parameters as used for the simulations (see
Table 2), but with an imaging depth of $45\text{\,}\mathrm{m}\mathrm{m}$.
Being acquired from a physical phantom, there is no average image as
reference, but the distinct regions in the phantom are constructed as regions
of homogeneous scattering. The image shows a background region, several string
reflectors, as well as two different circular hyper-echoic inclusions. The
result of SRAD (b), median (c) and bilateral (d) show the same general
patterns as with the simulated data—not succeeding in completely removing the
speckle. NLM (e), OBNLM (f) and our method (g) are again close to that goal.
While NLM shows the smoothest image among the three methods and shows a sharp
edge of the bright inclusion, it suppresses the string responses more. The
result of OBNLM is almost as smooth as NLM overall, but it shows
inhomogeneities around the bright inclusion and its outline is not as well
defined as with the other methods. Our method—while smooth locally—shows a
certain inhomogeneity in the background, but does not suppress the point
reflectors as strongly and has a clear demarcation of both inclusions.

(a) Original

(b) SRAD

(c) Median

(d) Bilateral

(e) NLM

(f) OBNLM

(g) Ours

Figure 7:  Qualitative comparison on an in-vivo image. The original speckled
image (a), the output of SRAD (b), median filter (c), bilateral filtering (d),
NLM (e), OBNLM (f) and ours (g). The marked region in the thyroid is used for
quantitative evaluation.

_In-vivo data._ In the in-vivo image in Figure 7—a cross-sectional view of the
thyroid of a healthy volunteer acquired with the same configuration as the
phantom data with $35\text{\,}\mathrm{m}\mathrm{m}$ depth—we can see that SRAD
(b), median (c) and bilateral filtering (d) perform better than in the
simulated and experimental cases, although NLM (e), OBNLM (f) and our method
(g) exhibit a significantly more homogeneous appearance in the thyroid and
especially OBNLM and our method retain more sharpness in comparison. While
OBNLM better captures the muscle interfaces over the carotid artery on the
right of the image compared to our method, it is also less smooth in the
overall area of the carotid. Our method does not retain said interfaces
completely, but is not disturbed by the clutter within the artery.

###  4.2 Quantitative evaluation

_Simulated data._ For the quantitative evaluation in the simulated case we
utilize the average images as an approximation of the speckle-free image. We
perform this for all 100 simulated images in the validation set and the
corresponding average images. The input images are scaled to the range $[0,1]$
before application of any of the methods. Table 6 shows the mean and standard
deviation of the mean squared errors (MSE) and the mean absolute differences
(MAD) between the considered methods and the average images.

Method | MSE ($\times 10^{-3}$) | MAD ($\times 10^{-2}$)
---|---|---
Input | 5.87 |  $\pm$ 1.31 | 5.86 |  $\pm$ 0.87
SRAD | 2.33 |  $\pm$ 0.86 | 3.61 |  $\pm$ 0.76
Median | 1.82 |  $\pm$ 0.25 | 3.25 |  $\pm$ 0.33
Bilateral | 1.48 |  $\pm$ 0.20 | 2.94 |  $\pm$ 0.30
NLM | 1.27 |  $\pm$ 0.09 | 2.70 |  $\pm$ 0.13
OBNLM | 1.48 |  $\pm$ 0.15 | 2.87 |  $\pm$ 0.20
Ours | 1.20 |  $\pm$ 0.11 | 2.60 |  $\pm$ 0.20
Table 6: Quantitative comparison on 100 simulated images from the validation
set. It shows the mean and standard deviation of the MSE and the MAD between
images processed with each methods (and the unprocessed input) against the
average images of nine speckle instances. Lowest error mean and standard
deviation for both MSE and MAD are highlighted in bold without implying
significance.

All methods besides SRAD and median show a similar performance with respect to
MSE and MAD, where our method yields the best result. In combination with the
clearly varying image appearance of the different methods, this emphasizes the
trade-offs between the methods compared here.

(a) Background

(b) Hyperechoic Inclusion

Figure 8: Image intensity distribution within homogeneous background region
(a) and hyperechoic inclusion (b) of the experimental images. Figure 9: Image
intensity distribution within homogeneous thyroid region of the in-vivo
images.

_Phantom and in-vivo data._ In the experimental and in-vivo images it is not
possible to observe different speckle instantiations and thus no average image
can be computed. Consequently, we focus on homogeneous image regions, namely
the hyperechoic inclusion and the background for the phantom dataset (see
6(a)) as well as the thyroid in the in-vivo dataset (see 7(a)). Within those
regions of assumed homogeneous scatterer distribution the ideal speckle
removal filter would result in a flat response. Figure 8 shows violin plots of
the intensity distribution in the marked regions of the phantom images, while
Figure 9 shows those of the in-vivo image regions. The results of NLM, OBNLM
and our method on the experimental image show very similar distributions;
SRAD, median and bilateral-filter however are less tightly packed. This
generally is matched by the evaluation for the in-vivo image data, but the
distribution of our method is less pronounced than NLM and OBNLM.

In Table 7, we show the image intensity standard deviations within the
homogeneous image regions. As expected from the violin plots, SRAD, median,
and the bilateral-filter exhibit larger standard deviations in the three
regions than NLM, OBNLM and our method. While our method does not outperform
NLM and OBNLM—which share best performance within the three regions and were
tuned on the thyroid image—it shows how well our method translates to real
ultrasound data, despite not having been trained with real images.

Region | Experimental ($\times 10^{-2}$) | In-Vivo ($\times 10^{-2}$)
---|---|---
Method | Background | Inclusion | Thyroid
Input | 7.10 | 7.58 | 8.96
SRAD | 3.55 | 3.96 | 5.89
Median | 3.37 | 3.83 | 6.05
Bilateral | 2.39 | 2.69 | 5.42
NLM | 1.49 | 0.96 | 4.21
OBNLM | 1.32 | 1.60 | 4.20
Ours | 1.61 | 1.51 | 4.54
Table 7: Quantitative evaluation on the experimental and in-vivo images and
shown in 6(a) and 7(a) respectively. The table shows the image intensity
standard deviation in homogeneous regions in the background and foreground of
the experimental images, and within the thyroid in the in-vivo case. Lowest
standard deviation per region highlighted in bold without implying
significance.

##  5 Discussion

The qualitative and quantitative evaluations show that Speckle2Speckle works
well in the simulated images, performing on-par with or even outperforming NLM
and OBNLM, as can be seen in Figure 5 and Table 6. It also performs well in
real acquisitions, although no real images have been used for training. There
are shortcomings however, especially regarding the overall smoothness as in
the background area and the reproduction of the point scatterers in Figure 6.
Yet, it seems the highlighting of interfaces translates well to real images,
with the exception of the linear structures in the muscle layer above the
carotid artery in Figure 7.

The execution time of our method is significantly lower than for NLM and
OBNLM. While the used NLM and OBNLM implementations (CPU only) take in the
range of minutes to execute, our method can be performed within tens of
milliseconds on a recent GPU, making it suitable for real-time use.

##  6 Conclusion

We have shown with Speckle2Speckle a technique for training despeckling
filters purely on simulated data as well as the generation of the simulated
data. The concept is intriguing, since it allows for the generation of
multiple speckle realizations from the same underlying geometry. The learned
filter translates well to real acquisitions, performing similar to NLM and
OBNLM, but only requiring a fraction of the execution time.

Since no clean training data is required, the number of simulations is
drastically reduced compared to approaches that would use average images were
used as target images. It also enables creation of matching maps with high-
level info on the structures present in the image that can be used to guide
the outcome of the method during training—as we have done with interfaces in
this case.

Including single point scatterers and/or strings in the simulated images could
improve the reproduction of those structures. They could be included in the
interface map used for sharpening or placed in a separate map all-together,
only depending on the desired appearance. This aspect of using additional
geometric maps to tune the result for particular locations could possibly be
exploited beyond speckle removal, for example for clutter filtering in blood
vessels within the same filter pass.

Acknowledgments

We acknowledge the important conversations with our colleagues Maximilian
Baust, Maria Tirindelli, and Stefan Wörz.

Ethical Standards

The work follows appropriate ethical standards in conducting research and
writing the manuscript, following all applicable laws and regulations
regarding treatment of animals or human subjects.

Conflicts of Interest

We declare we do not have conflicts of interest.

## References

  * Buades et al. (2005) Antoni Buades, Bartomeu Coll, and J-M Morel.  A non-local algorithm for image denoising.  In _2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05)_ , volume 2, pages 60–65. IEEE, 2005.
  * Cammarasana et al. (2021) Simone Cammarasana, Paolo Nicolardi, and Giuseppe Patanè.  A universal deep learning framework for real-time denoising of ultrasound images.  _arXiv preprint arXiv:2101.09122_ , 2021.
  * Coupé (2011) Pierrick Coupé.  Optimized bayesian non-local means filter, matlab package.  https://web.archive.org/web/20201015011836/https://sites.google.com/site/pierrickcoupe/softwares/denoising-for-medical-imaging/speckle-reduction/obnlm-package, 2011\.  [Online; accessed 7-December-2021].
  * Coupé et al. (2009) Pierrick Coupé, Pierre Hellier, Charles Kervrann, and Christian Barillot.  Nonlocal means-based speckle filtering for ultrasound images.  _IEEE Transactions on Image Processing_ , 18(10):2221–2229, 2009.
  * Göbl et al. (2018) Rüdiger Göbl, Nassir Navab, and Christoph Hennersperger.  Supra: open-source software-defined ultrasound processing for real-time applications.  _International Journal of Computer Assisted Radiology and Surgery_ , Mar 2018.  ISSN 1861-6429.  doi: 10.1007/s11548-018-1750-6.
  * Jensen (1996) Jørgen Arendt Jensen.  Field: A program for simulating ultrasound systems.  In _Proc. of the 6th Nordic Conf. on Human-Computer Interaction_ , Medical & Biological Engineering & Computing, pages 351–353. Springer, 1996.
  * Jensen and Svendsen (1992) Jørgen Arendt Jensen and Niels Bruun Svendsen.  Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers.  _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , 39(2):262–267, 1992.
  * Krull et al. (2019) Alexander Krull, Tim-Oliver Buchholz, and Florian Jug.  Noise2void-learning denoising from single noisy images.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 2129–2137, 2019.
  * Lehtinen et al. (2018) Jaakko Lehtinen, Jacob Munkberg, Jon Hasselgren, Samuli Laine, Tero Karras, Miika Aittala, and Timo Aila.  Noise2Noise: Learning image restoration without clean data.  In Jennifer Dy and Andreas Krause, editors, _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pages 2965–2974, Stockholmsmässan, Stockholm Sweden, 10–15 Jul 2018. PMLR.
  * Mohd Sagheer and George (2020) Sameera V. Mohd Sagheer and Sudhish N. George.  A review on medical image denoising algorithms.  _Biomedical Signal Processing and Control_ , 61:102036, 2020\.  ISSN 1746-8094.  doi: https://doi.org/10.1016/j.bspc.2020.102036.
  * Ronneberger et al. (2015) Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI’15)_ , pages 234–241. Springer, 2015.
  * Tomasi and Manduchi (1998) Carlo Tomasi and Roberto Manduchi.  Bilateral filtering for gray and color images.  In _Sixth International Conference on Computer Vision (IEEE Cat. No. 98CH36271)_ , pages 839–846. IEEE, 1998.
  * Virginia Image and Video Analysis, School of Engineering and Applied Science, University of Virginia (2006) Virginia Image and Video Analysis, School of Engineering and Applied Science, University of Virginia.  MATLAB and C implementation of speckle reducing anisotropic diffusion (SRAD).  https://web.archive.org/web/20171213062541/http://viva-lab.ece.virginia.edu/downloads.html, 2006\.  [Online; accessed 7-December-2021].
  * Yin et al. (2020) Da Yin, Zhongzheng Gu, Yanran Zhang, Fengyan Gu, Shouping Nie, Shaotong Feng, Jun Ma, and Caojin Yuan.  Speckle noise reduction in coherent imaging based on deep learning without clean data.  _Optics and Lasers in Engineering_ , 133:106151, 2020.
  * Yu and Acton (2002) Yongjian Yu and Scott T Acton.  Speckle reducing anisotropic diffusion.  _IEEE Transactions on image processing_ , 11(11):1260–1270, 2002.


Directly applying it to sfl, i.e., each satellite contacts the ps directly, leads to very slow convergence speed, especially in scenarios with terrestrial orchestration <cit.>. The question is whether the usage of isl resolves the connectivity bottleneck sufficiently to make synchronous orchestration feasible.

To this end, we measure the test accuracy with respect to the wall clock time for synchronous terrestrial orchestration in the and constellations, with and without isl. We complement this by an experiment with a ps in leo. Results for non-iid MNIST and CIFAR-10 are displayed in <ref>.
All scenarios show the distinct step function behavior first observed in <cit.> for the non-isl scenarios. This is caused by the ps being forced to wait for all results before updating the global model. As expected, the usage of isl, as proposed in <ref>, shortens the convergence time significantly in all scenarios.

With MNIST training usually showing very fast convergence, it can be easily observed that isl reduce the time to convergence in all scenarios by 4 hours. This leads to almost instantaneous convergence in the leo scenario, at least in relation to the training duration in conventional sfl, and reduces the convergence time by up to 50 for terrestrial orchestration.

[MNIST dataset]

xlabel=Time [$\si{\hour}$],
ylabel=Top-1 Accuracy,
ylabel near ticks,
minor x tick num = 4,
minor y tick num = 1,

xmin = 0,
ymin = 0,
ymax = .93,
xmax = 15,
xtick = 0,5,...,15,
legend pos=south east,
legend cell align=left,
cycle list name=default,
legend image post style=very thick,
no markers,

legend entries = isl, PS: LEO, , non-isl, PS: LEO, , isl, PS: GS, , non-isl, PS: GS, , isl, PS: GS, ,   non-isl, PS: GS,  ,
legend to name = fig:sync_mnist_cifar10_data:legend,
legend columns = 2,
transpose legend,
legend style=/tikz/every even column/.append style=column sep = 2mm,

[color=blue,line width=1.1pt] table [x=Time_FedISL_Sync_MEO, y=Acc_FedISL_Sync_MEO, col sep=comma] images/FedISL_FedNonISL_MNIST_NonIID.csv;

[color=blue,line width=1.3pt, dashed] table [x=Time_FedNonISL_Sync_MEO, y=Acc_FedNonISL_Sync_MEO, col sep=comma] images/FedISL_FedNonISL_MNIST_NonIID.csv;

[color=black!20!red,line width=1.3pt] table [x=Time_FedISL_Sync_Bremen_walker_star, y=Acc_FedISL_Sync_Bremen_walker_star, col sep=comma] images/FedISL_FedNonISL_MNIST_NonIID_walker_star.csv;

[color=black!20!red,line width=1.3pt, dashed] table [x=Time_FedNonISL_Sync_Bremen_walker_star, y=Acc_FedNonISL_Sync_Bremen_walker_star, col sep=comma] images/FedISL_FedNonISL_MNIST_NonIID_walker_star.csv;

[color=black!40!green,line width=1.1pt] table [x=Time_FedISL_Sync_Bremen, y=Acc_FedISL_Sync_Bremen, col sep=comma] images/FedISL_FedNonISL_MNIST_NonIID.csv;

[color=black!40!green,line width=1.3pt, dashed] table [x=Time_FedNonISL_Sync_Bremen, y=Acc_FedNonISL_Sync_Bremen, col sep=comma] images/FedISL_FedNonISL_MNIST_NonIID.csv;

[CIFAR-10 dataset]

xlabel=Time [$\si{\hour}$],
ylabel=Top-1 Accuracy,
ylabel near ticks,
minor x tick num = 4,
minor y tick num = 1,

xmin = 0,
ymin = 0,
ymax = .8,
xmax = 15,
xtick = 0,5,...,15,
legend pos=south east,
legend cell align=left,
cycle list name=default,
legend image post style=very thick,
no markers,

legend entries = ,
yticklabel pos = right,

[color=blue,line width=1.1pt] table [x=Time_FedISL_Sync_MEO, y=Acc_FedISL_Sync_MEO, col sep=comma] images/FedISL_FedNonISL_CIFAR.csv;

[color=blue,line width=1.3pt, dashed] table [x=Time_FedNonISL_Sync_MEO, y=Acc_FedNonISL_Sync_MEO, col sep=comma]images/FedISL_FedNonISL_CIFAR.csv;

[color=black!40!green,line width=1.1pt] table [x=Time_FedISL_Sync_Bremen, y=Acc_FedISL_Sync_Bremen, col sep=comma]images/FedISL_FedNonISL_CIFAR.csv;

[color=black!40!green,line width=1.3pt, dashed] table [x=Time_FedNonISL_Sync_Bremen, y=Acc_FedNonISL_Sync_Bremen, col sep=comma]images/FedISL_FedNonISL_CIFAR.csv;

[color=black!20!red,line width=1.3pt] table [x=Time_FedISL_Sync_Bremen_star, y=Acc_FedISL_Sync_Bremen_star, col sep=comma] images/FedISL_FedNonISL_CIFAR_star.csv;

[color=black!20!red,line width=1.3pt, dashed] table [x=Time_FedNonISL_Sync_Bremen_star, y=Acc_FedNonISL_Sync_Bremen_star, col sep=comma] images/FedISL_FedNonISL_CIFAR_star.csv;

Test accuracy with respect to wall-clock time. Synchronous orchestration for the MNIST and CIFAR-10 datasets with non-iid distributions, considering both terrestrial and non-terrestrial PS. i.e., the GS in Bremen and an LEO satellite. Note that isl and non-isl stand for the FedAvg with and without isl algorithms respectively.

The CIFAR-10 experiment was chosen to evaluate the convergence behavior of a more involved ml model. As such, convergence requires a large number of global iterations. From <ref>, we observe a more nuanced convergence behavior as compared to the MNIST experiment. In particular, leo-based orchestration benefits most from isl, while profits only in the initial learning phase. Surprisingly, it basically shows the same convergence speed as without isl, simply shifted in time by a few hours. The reason behind this is best understood from <ref>, which displays exactly the connectivity patterns of these two scenarios. Connectivity in the leo ps scenario changes from a sporadic towards a near-persistent connectivity pattern. Instead, for terrestrial orchestration, the connectivity windows grow significantly in duration, but the overall connectivity pattern is still sporadic. That is, the offline periods of individual orbital groups still lead to blocking in synchronous aggregation algorithms. However, this behavior is not universal to terrestrial orchestration. It depends on the constellation design and gs location, as can be seen by the improved convergence speed in the scenario. Moreover, combining synchronous orchestration and isl with the scheduling approach from <cit.> is expected to alleviate this problem in training scenarios with limited computational complexity.

Following the discussion in <ref> and <cit.>, a viable alternative to synchronous orchestration in sporadic connectivity scenarios is asynchronous ps operation. This is evaluated in <ref> for with terrestrial orchestration, with maximum time between updates $T_u$ set to 147 minutes. Convergence with asynchronous aggregation is improved over synchronous operation, but still involves at least a ten-fold increase in convergence time. Especially initial convergence speed, i.e., within the first few hours of training, is greatly increased. This leads to the conclusion that asynchronous aggregation is a potential solution if the system design suffers from sporadic connectivity. However, if changing the ps location is an option, it might be preferable to asynchronous orchestration.

[iid distribution]

yminorgrids = true,
legend entries = isl Async., isl Sync., non-isl Sync.,
xlabel=Time [h],
ylabel=Test Accuracy,
minor x tick num = 4,
minor y tick num = 1,
extra y ticks = .772,
xmin = 0,
xmax = 110,
ymin = 0,
ymax = 0.79,
grid = major,

legend cell align=left,

legend pos=south east,

legend style=font=

[color=blue, line width = 1.1pt] table [x=Time_FedISL_Async_Bremen_IID, y=Acc_FedISL_Async_Bremen_IID, col sep=comma]images/FedISL_Sync_Async_CIFAR10.csv;

[color=black!40!green,line width=1.1pt] table [x=Time_FedISL_Sync_Bremen_IID, y=Acc_FedISL_Sync_Bremen_IID, col sep=comma]images/FedISL_Sync_Async_CIFAR10.csv;

[color=black!20!red,line width=1.1pt] table [x=Time_FedNonISL_Sync_Bremen_IID, y=Acc_FedNonISL_Sync_Bremen_IID, col sep=comma]images/FedISL_Sync_Async_CIFAR10.csv;

[non-iid distribution]

yminorgrids = true,
legend entries = isl Async., isl Sync., non-isl Sync.,
xlabel=Time [h],
ylabel=Test Accuracy,
minor x tick num = 4,
minor y tick num = 1,
extra y ticks = .772,
xmin = 0,
xmax = 110,
ymin = 0,
ymax = 0.79,
grid = major,

legend cell align=left,

legend pos=south east,

legend style=font=,
yticklabel pos = right,

[color=blue, line width = 1.1pt] table [x=Time_FedISL_Async_Bremen_NIID, y=Acc_FedISL_Async_Bremen_NIID, col sep=comma]images/FedISL_Sync_Async_CIFAR10.csv;

[color=black!40!green,line width=1.1pt] table [x=Time_FedISL_Sync_Bremen_NIID, y=Acc_FedISL_Sync_Bremen_NIID, col sep=comma]images/FedISL_Sync_Async_CIFAR10.csv;

[color=black!20!red, line width = 1.1pt] table [x=Time_FedNonISL_Sync_Bremen_NIID, y=Acc_FedNonISL_Sync_Bremen_NIID, col sep=comma]images/FedISL_Sync_Async_CIFAR10.csv;

Test accuracy with respect to wall-clock time for Constellation. Synchronous and asynchronous orchestrations for the CIFAR-10 dataset with iid and non-iid distributions are included, where the PS is located in Bremen.

§.§ Gradient Sparsification and Transmission Load

<ref> shows the transmission load required for collecting the gradients in one orbit and sending them toward the PS, within a single global iteration. The results are based on the MNIST dataset with non-iid distribution. The figure depicts the total transmitted data as a function of the number of satellites in an orbit.
The total transmitted data in FedAvg with isl algorithm are compared for two different cases, with and without incremental aggregation (IA), respectively. When IA is used, each satellite aggregates its gradient updates with the received ones, according to the method described in <ref>. Without IA, each satellite only forwards its gradients and the received ones to the other satellite or the ps, without performing any aggregation. We have evaluated the total transmitted data for different sparsification ratios $q=0.01$, $q=0.1$, as well as without sparsification, which corresponds to $q=1$.

<ref> shows the number of the transmitted elements for different sparsification ratios, normalized with the length of the sparsified gradient vector $n_a$. Clearly, this comparison is only meaningful if the IA does not affect the model accuracy and we have verified that this is the case since when IA is applied, satellites only perform the aggregation operations, which is equivalent to the operations performed by the ps without IA. The final parameter values obtained, with or without IA, remain unchanged. By comparing the cases with and without IA, we can observe a significant reduction in the transmission load when the number of satellites increases, $>10$ times with $40$ satellites and no sparsification. This reduction is particularly valuable as the number of global iterations increases since it consistently reduces the transmitted data in each iteration, thereby addressing the problem of overall communication cost in fl. In addition to the simulation curves, <ref> also depicts and verifies the theoretical upper bound for the transmission load, derived in <ref>.

The normalized transmitted data increases as the sparsification ratio $q$ increases. To explain this, recall that, when IA is applied under the assumption of sparsification, the satellites transmit the non-zero elements along with their respective indices.
When the indices of non-zero elements differ among satellites, each satellite must transmit its own elements and the elements from the other satellites, along with their indices. This leads to an increase in the number of elements transmitted. However, as the sparsification ratio $q$ increases, it becomes more likely that the same indices will have the non-zero entries at different satellites. Each index that is present at multiple satellites needs to be transmitted only once, along with the aggregated value from all the satellites, which reduces the transmitted data.

<ref> shows the absolute value, rather than the normalized, of the total transmitted data. The storage size of each element $\omega$ is set to 32 bits. For the case with IA, with sparsification (q>1), a total number of 45 bits is required to transmit each element, which is the sum of the storage space for each element and its index. In the considered model, the index storage space is 13 bits. The difference between the schemes with and without IA becomes substantial as $q$ increases, owing to the increased correlation among the indices of non-zero elements across the satellites. Additionally, when the number of satellites increases, the reduction in data size achieved with IA is significant to the extent that the case without sparsification $q=1$ outperforms the case without IA and a sparsification ratio of $q=0.1$.

font = ,
ylabel style = align=center,
xlabel=Satellites (orbital plane),

ylabel=Total transmitted data [normalized],

minor x tick num = 1,
minor y tick num = 1,
no markers,
legend columns = 3,
legend cell align=left,
legend pos = north west,
legend style = font=,
legend image post style = very thick,

[col sep=comma]images/IA_vs_q.dat

(b001) at (axis cs:35,325.5);
(s001) at (axis cs:35,304.4);
(b01) at (axis cs:35,200.8);
(s01) at (axis cs:35,166.5);
(b1) at (axis cs:35,35);

[line width = .4pt, densely dotted] let 1 = (b001), 2 = (s001), 1 = veclen((1-2, 1-2)) in ($(b001)!.5!(s001)$) ellipse [y radius = 1, x radius = 3pt];
[line width = .4pt, densely dotted] let 1 = (b01), 2 = (s01), 2 = veclen((1-2, 1-2)), 1 = 4pt in ($(b01)!.5!(s01)$) ellipse [y radius = 2, x radius = 1];
[line width = .4pt, densely dotted](b1) ellipse [y radius = 6pt, x radius = 3pt];

legend image with text=simulated

legend image with text=bound

empty legend

[] table [x=sats, y=noIA] ;

empty legend

empty legend
no IA

[black!20!red] table [x=sats, y=simulated_0.01] ;

[densely dashed, black!20!red] table [x=sats, y=bound_0.01] ;

empty legend
$q = 0.01$

[black!40!green] table [x=sats, y=simulated_0.1] ;

[densely dashed, black!40!green] table [x=sats, y=bound_0.1] ;

empty legend
$q = 0.1$

[blue] table [x=sats, y=IA] ;

empty legend

empty legend
no spars.

[line cap = round] let 1 = (b001), 2 = (s001), 2 = veclen((1-2, 1-2)), 1 = 3pt in ($(b001)!.5!(s001)$) + (76:1 and 2) arc (76:348:1 and 2) node[anchor=north west,font=, inner sep = 1pt, fill = white] $q = 0.01$;
[line cap = round] let 1 = (b01), 2 = (s01), 2 = veclen((1-2, 1-2)), 1 = 4pt in ($(b01)!.5!(s01)$) + (50:1 and 2) arc (50:330:1 and 2) node[anchor=north west,font=, inner sep = 1pt, fill = white] $q = 0.1$;
[line cap = round] (b1) + (25:3pt and 6pt) arc (25:340:3pt and 6pt) node[anchor=north west,font=, inner sep = 1pt, fill = white] $q = 1$;

font = ,
xlabel=Satellites (orbital plane),
ylabel=Total transmitted data [/global iteration],
ylabel near ticks,
minor x tick num = 1,
minor y tick num = 1,
scaled y ticks = base 10:-6,
ytick scale label code/.code=,
no markers,
legend columns = 3,
legend cell align=left,
legend pos = north west,
legend style = font=,
legend image post style = very thick,
yticklabel pos = right,

[col sep=comma]images/totaldata_IA.dat

legend image with text=IA

legend image with text=no IA

empty legend

[blue] table [x=sats, y=simulated_1] ;

[densely dashed, blue] table [x=sats, y=no_IA_1] ;

empty legend
no spars.

[black!40!green] table [x=sats, y=simulated_0.1] ;

[densely dashed, black!40!green] table [x=sats, y=no_IA_0.1] ;

empty legend
$q = 0.1$

[black!20!red] table [x=sats, y=simulated_0.01] ;

[densely dashed, black!20!red] table [x=sats, y=no_IA_0.01] ;

empty legend
$q = 0.01$

Total transmission data in one global iteration for three different sparsification ratios; a) $q=0.01$ and b) $q=0.1$, and c) without sparsification (q=1). IA and no IA stand for the proposed algorithm with and without incremental aggregation respectively.

The proposed FedAvg with isl effectively offloads communication from gsls to isls. This is especially relevant due to the fact that the isls are increasingly implemented as  high-throughput optical links <cit.>.
When FedAvg with isl is used, the number of required gsls is independent of the number of satellites in each orbit and remains constant at 2 for each global iteration, one for uplink and one for downlink, respectively. However, when FedAvg does not use isls, the number of gsls increases linearly with the number of satellites, as each satellite has to transmit and receive the parameters on its own.

§ DISCUSSION AND CONCLUSIONS

We have designed a clustered fl system tailored to distributed ml in modern satellite megaconstellations. It efficiently uses isl to avoid connectivity bottlenecks that impair convergence speed. Owing to a fl-specific in-network aggregation strategy, this is achieved without increasing the total transmitted amount of data. Indeed, it can be expected that this method actually decreases the communication cost due to relying on cheaper-to-operate isl instead of long distance gsl. Moreover, the communication load is evenly spread among the network instead of focusing all communications on direct worker-to-ps links.
This strategy requires a
careful predictive route planning
in order to operate successfully in a dynamic network topology imposed by orbital mechanics. We have developed the necessary routing algorithms for this in conjunction with a rigorous distributed system design that inflicts a low overhead for synchronization and consensus finding.
The proposed system is compatible with a wide variety of fl algorithms, supports gradient sparsification, and includes facilities for asynchronous clustered aggregation. While not explicitly mentioned, the system's extension to incorporate worker scheduling is rather straightforward.

We have evaluated the performance of the proposed system for a few carefully chosen examples. The results highlight the major increase in convergence speed to to isl and the bandwidth effectiveness of our in-network aggregation approach.

Regarding future work, it is interesting to revise the assumption adopted in this paper that
all participants in the fl process are trustworthy, under complete control, and the software operates nominally. For instance, if the satellites belong to different operators, the use of ISLs may incur a certain cost and/or a group for satellites may decide not to participate in the FL process, which would require an extension of the proposed algorithms.

§ STUFF THAT NEEDS TO GO SOMEWHERE ELSE

§.§ From the communications model

A free-space environment is considered for the ISL connections which is
concerned with the free-space path loss and the thermal noise power.
The maximum achievable data rate for the transmission between satellites $k$ and $i$ with the Gaussian channels is

\begin{equation} \label{eq: rate}
r(k,i) = B \log_2 \left(1+\mathrm{SNR}(k,i)\right),
\end{equation}

where $B$ is the channel bandwidth, and $\rm{SNR}(k,i)$ is the signal to noise ratio (SNR)  which can be expressed as <cit.>

\begin{equation} \label{eq: SNR}
\mathrm {SNR}(k,i) = \frac{P_t G_k(i) G_i(k) }{N_0  L(k,i)}
\end{equation}

if $k$ and $i$ are visible and zero otherwise. The average antenna gain of satellite $j$ towards satellite $l$ is
$G_j(l)$, $N_0 = k_B T B$ is the total noise power with $k_B = 1.380649 \times
10^{-23}\,\si{\joule/\kelvin}$ being the Boltzmann constant, $T$ is the receiver noise temperature. The free space path loss $L(k,i)$ is

\begin{equation} \label{eq: LOSS}
L(k,i) = \left(\frac{4\pi f_c d(k,i)}{c}\right)^2,
\end{equation}

where $f_c$ and $c$ are the carrier frequency and the speed of light, respectively. The rate between satellites and the GS is derived based on the longest distance possible for the connections in each orbit.
Time to transmit the effective gradient parameters from satellite $k$ to satellite $i$ is the sum of propagation and transmission times as

\begin{equation} \label{eq: Time for transmitting}
t_c(k, i)= \underbrace{\frac{S(\vec{g}_k(\vec{w}_k^{n,I}))}{r(k,i)}}_{\mathclap{\text{Transmission}}} + \underbrace{\frac{d(k,i)}{\mathstrut c}}_{\mathclap{\text{Propagation}}},
\end{equation}

where $S(\vec{g}_k(\vec{w}_k^{n,I}))$ is the size of $\vec{g}_k(\vec{w}_k^{n,I})$ in bits. For transmitting the model parameters, $S(\vec{g}_k(\vec{w}_k^{n,I}))$ is substituted with $S(\vec{w}_k^{n,I})$ and we assume $S(\vec{g}_k(\vec{w}_k^{n,I})) = S(\vec{w}_k^{n,I})$. Moreover, note that the size of $\vec{g}_k(\vec{w}_k^{n,I})$ is model-dependent and remains fixed during the whole training process. Likewise, this model can be used for ground-to-satellite links if necessary.

§.§ visiting pattern

somewhere we should have a visiting pattern

Visiting patterns between 12 satellites and a PS located in Bremen over 12 hours, where $k_{p,i}$ denotes the $i^{th}$ satellite in orbit $p$.TODO (Nasrin): Please use exactly the same style as in the magazine paper. This way, we kind of create a “visual identity” that makes our work easy to recognize across papers.

§.§ Async update rules

* Satellite $k$ gets weight at $n_0$

* Satellite $k$ computes updates $\vec w_k^{n_0,0}, \vec w_k^{n_0,1}, \vec w_k^{n_0,2}, \dots, \vec w_k^{n_0,i}$

* Satellite connects at epoch $n$ and delivers update $\vec w_k^{n_0,i}$

from the letter:

\begin{equation}
\vec w^{n+1} = \vec w^n - \frac{D_k}{D} \left( \vec w_k^{n_0,0} - \vec w_k^{n_0,i} \right)
\end{equation}

\begin{equation}
\vec w^{n+1} = \vec w^n - \frac{D_k}{D} \left( \vec w_k^{n_k-1,I} - \vec w_k^{n_k,I} \right)
\end{equation}

transform this:

\begin{align}
\vec w^{n+1} &= \vec w^n - \frac{D_k}{D} \left( \vec w_k^{n_0,0} - \vec w_k^{n_0,i} \right) \\
&= \vec w^n + \frac{D_k}{D} \left( \vec w_k^{n_0,i} - \vec w_k^{n_0,0} \right) \\
&= \vec w^n + \frac{1}{D} D_k \vec g_k(\vec w_k^{n_0})
\end{align}

\begin{align}
\vec w^{n+1} &= \vec w^n - \frac{D_k}{D} \left( \vec w_k^{n_k-1,I} - \vec w_k^{n_k,I} \right)  \\
&= \vec w^n - \frac{D_k}{D} \left( \vec w_k^{n_k-1,I} - \vec w_k^{n_k-1,0} + \vec w_k^{n_k-1,0} - \vec w_k^{n_k,I} + \vec w_k^{n_k,0} - \vec w_k^{n_k,0} \right)  \\
&=  \vec w^n - \frac{D_k}{D} \left( \vec g_k(\vec w_k^{n_k-1,I}) + \vec w_k^{n_k-1,0} -\vec g_k( \vec w_k^{n_k,I}) - \vec w_k^{n_k,0} \right) \\ & = \vec w^n - \frac{D_k}{D}  \vec w_k^{n_k-1,0} + \frac{D_k}{D} \vec w_k^{n_k,0} -\frac{D_k}{D} \left( \vec g_k(\vec w_k^{n_k-1,I})  -\vec g_k( \vec w_k^{n_k,I}) \right)
\end{align}

\begin{align}
\vec w^{n+1} &= \vec w^n - \frac{D_k}{D} \left( \vec w_k^{n_k,I} - \vec w_k^{n_k+1,I} \right)  \\
&= \vec w^n - \frac{D_k}{D} \left( \vec w_k^{n_k,I} - \vec w_k^{n_k,0} + \vec w_k^{n_k,0} - \vec w_k^{n_k+1,I} + \vec w_k^{n_k+1,0} - \vec w_k^{n_k+1,0} \right)  \\
&=  \vec w^n - \frac{D_k}{D} \left( \vec g_k(\vec w_k^{n_k,I}) + \vec w_k^{n_k,0} -\vec g_k( \vec w_k^{n_k+1,I}) - \vec w_k^{n_k+1,0} \right) \\ & = \vec w^n - \frac{D_k}{D}  \vec w_k^{n_k,0} + \frac{D_k}{D} \vec w_k^{n_k+1,0} -\frac{D_k}{D} \left( \vec g_k(\vec w_k^{n_k,I})  -\vec g_k( \vec w_k^{n_k+1,I}) \right)
\end{align}

Note that this is the same update rule as in the synchronous case, but without the sum. This is in stark contrast to the rule you state in equation <ref>. Please check where I'm wrong!
Maybe, we can hold the same formulation with $w$ for the asynchronous update due to the extra terms that is extracted and make it a bit complex .
Can we define the $w$ part in equation (18) as $w^{'}$ and then say that satellite just transmit the $g_k(w_k^{n_k-1})$ and the PS makes the new version.

For the global model parameter in an asynchronous version, we can write:

\begin{align}
\vec w^{n}=\sum_{k=1}^{K} \frac{D_k}{D} \vec g_k(w_k^{n_k,I})+\sum_{k=1}^{K} \frac{D_k}{D}\vec w_k^{n_k,0}
\end{align}

Then, for update with a new received gradients from one satellite:

\begin{align}
\vec w^{n+1}= \vec w^{n} + \frac{D_k}{D}((\vec g_k(\vec w_k^{n_k+1,I}) - \vec g_k(w_k^{n_k,I}))+(\vec w_k^{n_k+1,0} - \vec w_k^{n_k,0}))
\end{align}

For the satellite side, it only needs to transmit $\vec g_k(\vec w_k^{n_k+1,I})$. For the PS side, it needs to put some memory for $g_k(\vec w_k^{n_k+1,I})$, $w_k^{n_k+1,I}$, and $w_k^{n_k,0}$.

And if we want to have a simple notation for the equation, we can write:

\begin{align}
\vec w^{n+1}= \vec w^{n} + \frac{D_k}{D}(\vec \Delta_ {\vec g_k}+ \vec \Delta_{\vec w_k})
\end{align}

where $\vec \Delta_ {\vec g_k} = \vec g_k(\vec w_k^{n_k+1,I}) - \vec g_k(w_k^{n_k,I})$ and $\vec \Delta_{\vec w_k}=\vec w_k^{n_k+1,0} - \vec w_k^{n_k,0}$.