# Unsupervised Anomaly Detection using Aggregated Normative Diffusion

Alexander Frotscher Jaivardhan Kapoor∗ Thomas Wolfers Christian F. Baumgartner
University of Tübingen
Equal contribution.

###### Abstract

Early detection of anomalies in medical images such as brain MRI is highly
relevant for diagnosis and treatment of many conditions. Supervised machine
learning methods are limited to a small number of pathologies where there is
good availability of labeled data. In contrast, _unsupervised_ anomaly
detection (UAD) has the potential to identify a broader spectrum of anomalies
by spotting deviations from normal patterns. Our research demonstrates that
existing state-of-the-art UAD approaches do not generalise well to diverse
types of anomalies in realistic multi-modal MR data. To overcome this, we
introduce a new UAD method named Aggregated Normative Diffusion (ANDi)111Code
is available at https://github.com/alexanderfrotscher/ANDi.. ANDi operates by
aggregating differences between predicted denoising steps and ground truth
backwards transitions in Denoising Diffusion Probabilistic Models (DDPMs) that
have been trained on pyramidal Gaussian noise. We validate ANDi against three
recent UAD baselines, and across three diverse brain MRI datasets. We show
that ANDi, in some cases, substantially surpasses these baselines and shows
increased robustness to varying types of anomalies. Particularly in detecting
multiple sclerosis (MS) lesions, ANDi achieves improvements of up to 178% in
terms of AUPRC.

##  1 Introduction

Detection of anomalies in medical imaging data is highly important for early
diagnoses and timely treatment of diseases. In order to detect anomalies,
clinicians typically evaluate _multi-modal_ image data, consisting of a range
of different medical images of the same patient that provide different
information about potential pathologies. In neuroimaging, routinely, a range
of different magnetic resonance (MR) imaging sequences are acquired that
together offer a comprehensive view of the brain.

Figure 1: Aggregated Normative Diffusion (ANDi). First, we train a DDPM model
using our proposed Gaussian pyramidal noise on healthy brain slices to
approximate the normative distribution $q_{norm}(x_{0})$. In order to obtain
an anomaly map for a possibly anomalous image $x_{0}$, we first partially
noise it using the Gaussian forward process (indicated by the gray arrow). We
then calculate the pixel-wise Euclidean distance $d_{t}$ between the ground
truth backwards transition $\tilde{x}_{t}=\mu_{q}(x_{t},x_{0})$ and the
denoising step $\hat{x}_{t}=\mu_{\theta}(x_{t},t)$ for a partial range of $t$.
The denoising step towards $\hat{x}_{t}$ can be thought of as _normative
diffusion_ as it is taking one step towards the normative distribution.
Finally, using the geometric mean $A_{\text{gm}}$, we aggregate deviations
over $d_{t}$.

Due to the diversity of possible anomalies, their detection is very
challenging, time-consuming, and prone to errors. Therefore, there has been
substantial interest in automating anomaly detection. The majority of research
effort on automated anomaly detection has been focused on supervised learning
(_e.g_. [10, 12, 19, 17]). However, supervised algorithms must be tailored to
each type of anomaly and hinge on large-scale hand-annotated data. This limits
their clinical applicability and prevents their application to rare or
understudied diseases.

Recently, there has been an increased interest in unsupervised anomaly
detection (UAD) methods. UAD attempts to learn a concept of normality and
label anything that deviates from normality as an anomaly. Recent methods
often apply generative modeling to learn a _normative reference distribution_
from a large-scale population of healthy images. To this end, prior work has
employed techniques such generative adversarial networks (GANs) [28, 27],
autoencoders (AEs) [1, 26], variational autoencoders (VAEs) [4, 5, 8], and
most recently denoising diffusion probabilistic models (DDPMs) [31, 3, 23].
Most prior work then uses the trained networks to reconstruct a pseudo-healthy
image from an input image potentially containing an anomaly [5, 8, 14, 3, 23].
An approximate segmentation of the anomaly is then obtained via the
thresholded difference between the original and the pseudo-healthy image.
Recent DDPM-based approaches often perform the reconstruction in a multi-stage
fashion: First, a rough map anomaly map is computed, and then the highlighted
areas are in-painted with healthy-looking tissue [3, 23]. While these multi-
stage approaches may boost performance, the additional steps also add to the
complexity and execution time of the UAD pipeline.

Since anomaly maps are typically continuous, a key question in the UAD
literature is how to set the threshold to binarize a segmentation, and
consequently, whether a sample can be considered anomalous [2]. In practice,
many previous approaches directly tune the threshold for lesion detection on
labeled data from the target domain, which overestimates detection performance
and limits the approaches’ generalizability to other domains [23, 3, 14, 8, 2,
5]. Conceptually, such approaches cannot be considered truly unsupervised
since they require labeled data for determining the optimal operating point.

Knowledge about the target domain can also be introduced via the selection of
the modalities on which the algorithm is trained [31, 3, 8, 23]. For example,
tumors are easiest to detect on T2-weighted or FLAIR images, while chronic
post-stroke changes are best detected in T1-weighted images. Evaluating only
the modality on which the target anomaly is easiest to identify, again,
overestimates detection performance in a truly unsupervised scenario. Most
prior work is currently evaluated with a single image modality (_e.g_. [31, 3,
23, 8]), which makes it unclear how well these methods would generalize to the
clinically relevant multi-modal scenario.

Lastly, a number of techniques use implicit assumptions about the appearance
of the target anomalies. For instance, a recent approach based on denoising
autoencoders (DAEs) employs coarse noise during training that mimics the
structure of large tumors on which the approach is evaluated [14]. However, as
we show, the approach does not generalize to smaller multiple sclerosis (MS)
lesions.

In this paper, we propose Aggregated Normative Diffusion models ANDi, a novel
DDPM-based unsupervised anomaly detection algorithm that requires no
inpainting stage, is simple to use, and achieves state-of-the-art performance
on diverse multi-modal neuroimaging datasets in a realistic multi-modal MR
imaging setting. Anomaly maps are obtained by directly aggregating deviations
between the actual and the learned transitions in the denoising process (see
Fig. 1). Key to the performance of our proposed method is the use of multi-
resolution pyramidal (Gaussian) noise for training the DDPM. This
substantially boosts performance and allows the method to generalize to
different types of pathologies. We also introduce a novel mechanism for
aggregating deviations in the diffusion process based on the geometric mean
that leads to further improvements. Crucially, ANDi works with no a priori
modality selection and does not require threshold tuning on labeled data.

Our contributions are as follows:

  1. 1.

We demonstrate that current state-of-the-art UAD methods are not robust to
varying lesion types in the clinically relevant multi-modal MR scenario.

  2. 2.

We propose ANDi, a simple and fast DDPM-based UAD method based on aggregating
deviations between the actual and the learned denoising transition, a novel
pyramidal noise, and a novel aggregation mechanism.

  3. 3.

We introduce a novel domain-agnostic thresholding strategy based on Yen
thresholding [32] that can be applied to any method and allows evaluation of
segmentation accuracy without domain-specific threshold tuning.

  4. 4.

We show that our method consistently outperforms related work on three diverse
multi-modal brain anomaly datasets, both in a domain-agnostic evaluation and
in terms of the commonly used “ceil” Dice score.

##  2 Related Work

The first methods tackling anomaly detection in medical images originated in
the field of neuroimaging. The earliest approaches rely on building a
population-based brain atlas reflecting an “average” normal subject and
registering potentially anomalous images to that atlas. Unmatchable areas are
then used to identify potential pathologies [6, 29, 24]. A limitation of these
approaches is that the registration step may fail in the presence of
anomalies. While this issue can be partially mitigated with registration
frameworks that are more robust to anomalies (e.g.[18]), handling anomalies
that deviate significantly from the atlas remains a problem. Moreover, the
creation of an atlas can be subjective and may not represent all patient
populations equally.

More recently, Meissen _et al_. [20] proposed a simple semi-supervised
baseline approach for tumor and white matter lesion segmentation in which the
pre-processed FLAIR modality is simply thresholded in image space. The optimal
threshold is determined with a labeled validation set. This naive approach
works surprisingly well for tumors and white matter lesion detection, as they
appear bright compared to the surrounding tissue in FLAIR images. However,
this approach is data-specific and will not generalize to more complex types
of anomalies. In this paper, we use this approach as a semi-supervised
baseline to our proposed fully unsupervised method.

###  2.1 Deep Learning-based UAD in Medical Imaging

The emergence of deep learning-powered generative models such as VAEs [16], or
GANs [11] has spurred a new generation of UAD methods. In one of the earliest
approaches, Sato _et al_. [26] use an AE framework to detect anomalies in
computed tomography images of the head. This work was followed by a series of
contributions based on AEs [1], VAEs [7, 8, 4, 5], and GANs [28, 27]. For the
most part, these techniques rely on generating a pseudo-healthy image from a
potentially anomalous image.

Kascenas _et al_. [14] proposed to train a denoising autoencoder (DAE) on
multi-modal brain MR images with a custom “coarse” noise. While simple, this
approach reached state-of-the-art performance and, to date, was not
conclusively outperformed by newer approaches. Indeed, in more recent work,
Kascenas _et al_. [15] showed that the technique performs competitively when
compared to recent DDPM-based approaches. Moreover, this technique is one of
the few to consider the realistic multi-modal scenario where the network is
trained and evaluated with four diagnostically important MR sequences per
patient.

Recently, DDPMs [13] have been shown to reach unprecedented performance for
modeling high-dimensional probability distributions, and a number of recent
works build on this approach for UAD [31, 3, 23].

Similar to our proposed work, Pinaya _et al_. [23] aggregate the Euclidian
distances between the ground truth backwards transition and the denoising step
of the diffusion model. In contrast to our work, these aggregated distances
are used to obtain a preliminary anomaly mask, which is then inpainted via
another pass through the DDPM model. Interestingly, Kascenas _et al_. [15]
showed that omitting the inpainting step leads to superior performance,
however, the method still does not outperform the simpler DAE [14].

In closely related work, Wyatt _et al_. [31] introduce AnoDDPM, a pseudo-
healthy generation method based on DDPMs. Similar to our work, the authors
employed a multi-scale noise in the training of the DDPM in order to capture
anomalies at different scales. Specifically, the authors use multi-scale
simplex noise [22]. However, generating simplex noise incurs a comparatively
large computational cost. In contrast, our proposed pyramidal Gaussian noise
is much faster and easier to generate, and empirically leads to a better
performance within our framework.

Very recently, Bercea _et al_. [3] proposed a DDPM-based UAD framework based
also on pseudo-healthy image generation. The pipeline uses the DDPM to produce
an initial anomaly likelihood map, which is then in-painted in an iterative
process. The method was shown to outperform the earlier AnoDDPM [31] on stroke
lesions. However, as we will show in our experiments, the method is currently
not outperforming the simple DAE approach [14].

##  3 Denoising Diffusion Probabilistic Models

Algorithm 1 ANDi for unsupervised anomaly detection.

Input: Initial volume $v\in\mathbb{R}^{H\times W\times D\times C}$, $T_{l}$,
$T_{u}$, noise distribution $N_{\epsilon}$, thresholding $h:a\rightarrow s$.

1:$a_{m}\leftarrow\\{\\}$ $\triangleright$ empty multi-modal anomaly map

2:for each $x_{0}\in\mathbb{R}^{H\times W\times C}$ from $v$ do

3: for each $t\in\\{T_{l},\ldots,T_{u}\\}$ do

4: $\epsilon_{t}\sim N_{\epsilon}$ $\triangleright$ Sample noise

5: $x_{t}\leftarrow\sqrt{\bar{\alpha}_{t}}\cdot
x_{0}+\sqrt{1-\bar{\alpha}_{t}}\cdot\epsilon_{t}$$\triangleright$ Noise the
image

6: $\tilde{x}_{t-1}\leftarrow\mu_{q}(x_{t},x_{0})$ $\triangleright$ Expected
transition

7: $\hat{x}_{t-1}\leftarrow\mu_{\theta}(x_{t},t)$ $\triangleright$ Predicted
transition

8: $d_{t}\leftarrow\left(\tilde{x}_{t-1}-\hat{x}_{t-1}\right)^{2}$
$\triangleright$ Get deviation (Eq. 6)

9: end for

10: $i_{m}=A_{\text{gm}}(\\{d_{t}\\})$$\triangleright$ Aggregate per modality
(Eq. 8)

11: $a_{m}\oplus i_{m}$ $\triangleright$ Concatenate volume slice

12:end for

13:$a=\max_{m}a_{m}$ $\triangleright$ Max over modalities (Sec. 4.5)

14:$s\leftarrow h(a)$ $\triangleright$ Thresholding (Sec. 4.5)

Output: Binarized anomaly map $s\in\\{0,1\\}^{H\times W\times D}$.

DDPMs [13] are a class of generative models that aim to approximate a data
distribution $q(x_{0})$ given some training data from that distribution. A
DDPM consists of a fixed forward diffusion (noising) process and a learned
reverse denoising process. The forward process gradually adds Gaussian noise
$\epsilon$ to an input image $x_{0}$ in a fixed forward process
($x_{t-1}\xrightarrow{\text{noise}}x_{t}$), transforming it into Gaussian
noise over $T$ steps. A neural network then learns to reverse this process
iteratively ($x_{t}\xrightarrow{\text{denoise}}x_{t-1}$) to be able to
generate samples from the original data distribution.

The forward noise step from the original image $x_{0}$ to the image $x_{t}$ at
timestep $t$ is given by the Gaussian transition

$q(x_{t}|x_{0})=\mathcal{N}\left(\sqrt{\bar{\alpha}_{t}}x_{0},(1-\bar{\alpha}_{t})I\right),$
(1)

where $\alpha_{1}\ldots\alpha_{T}$ are scaling factors that follow a fixed
schedule, and $\bar{\alpha}_{t}=\prod_{k=1}^{t}\alpha_{k}$. To generate an
image from noise, we require sampling from the reverse transition
$q(x_{t-1}|x_{t})$, which is intractable. However, as described in [13], we
can apply Bayes rule to compute the closed-form expression for the reverse
transitions conditioned on $x_{0}$ as
$q\left(x_{t-1}|x_{t},x_{0}\right)=\mathcal{N}\left(\mu_{q}(x_{t},x_{0}),\Sigma_{q}(t)\right)$,
where

$\displaystyle\mu_{q}(x_{t},x_{0})=\frac{\sqrt{\alpha_{t}}(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_{t}}x_{t}+\frac{\sqrt{\bar{\alpha}_{t-1}}(1-\alpha_{t})}{1-\bar{\alpha}_{t}}x_{0}.$
(2)

We train the model to approximate the true denoising distribution
$q(x_{t-1}|x_{t})$ by the variational posterior

$p_{\theta}(x_{t-1}|x_{t})=\mathcal{N}\left(\mu_{\theta}(x_{t},t),\Sigma_{\theta}(x_{t},t)\right),$
(3)

for each step $t$, where a neural network is used to estimate the mean
$\mu_{\theta}(x_{t},t)$ while the variance $\Sigma_{\theta}(x_{t},t)$ is set
to
$\Sigma_{q}(t)=\frac{(1-\alpha_{t})(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_{t}}I$.

The model parameters $\theta$ are learned by maximizing the Evidence Lower
Bound on the log-likelihood of the joint distribution $q(x_{0}\ldots x_{T})$,
which can be decomposed into a sum of KL divergence terms of the form
$D_{KL}\left[q(x_{t-1}|x_{t},x_{0})\|p_{\theta}(x_{t-1}|x_{t})\right]$ for
each denoising step $t$. The training objective can thus be written as

$\displaystyle\underset{\theta}{\arg\min}\leavevmode\nobreak\
\mathbb{E}_{x_{0\ldots
T}}\sum_{t=1}^{T}\Sigma^{-1}_{q}(t)\|\mu_{\theta}(x_{t},t)-\mu_{q}(x_{t},x_{0})\|^{2}.$
(4)

In practice, a simpler training objective is used, where the scaling terms
$\Sigma^{-1}_{q}(t)$ are omitted and the predicted and target means are
reparameterized in terms of the noise $\epsilon$, i.e.

$\displaystyle\underset{\theta}{\arg\min}\leavevmode\nobreak\
\mathbb{E}_{x_{0},\epsilon_{0}}\sum_{t=1}^{T}\|\epsilon_{\theta}(x_{t},t)-\epsilon_{0}\|^{2},\text{with}$
(5)

$\epsilon_{0}$ being the noise added to $x_{0}$ to get $x_{t}$. Regardless of
reparameterization, this objective essentially consists of a sum of the
discrepancies between the predicted and target means. In the following
section, we demonstrate how the form of the optimization objective can be used
to create anomaly maps by aggregating these discrepancies.

##  4 Aggregated Normative Diffusion ANDi

Given a potentially anomalous multi-modal volume $v\in\mathbb{R}^{H\times
W\times D\times C}$, the aim of UAD is to obtain a pixel-wise anomaly map $a$
indicating the likelihood of an anomaly being present, where $C$ is the number
of image modalities, and $H,W,D$ are the image dimensions. An anomaly is
anything that deviates from the anatomy of healthy subjects. In the multi-
modal setting, we initially obtain an anomaly map $a_{m}$ for each channel,
which then needs to be reduced to a pixel-wise single-channel anomaly map $a$.
From the anomaly map $a$, a binary anomaly segmentation mask
$s\in\\{0,1\\}^{H\times W\times D}$ can then be generated using a thresholding
mechanism. In practice, $s$ could, for example, be the segmentation of a tumor
or an MS lesion.

In the following, we first describe how we train a DDPM to model the normative
reference distribution (Sec. 4.1), and how that model can be used to obtain
anomaly maps via aggregation of differences between the ground truth backwards
transitions and the denoising step (Sec. 4.2). In Sec. 4.3, we describe a
novel approach for aggregating those differences, which leads to more robust
anomaly detection, and in Sec. 4.4 we describe a pyramidal noise type that can
be used in DDPM training and further improves performance. Lastly, in Sec. 4.5
we describe how multi-channel anomaly maps $a_{m}$ are reduced to single-
channel anomaly maps $a$, and how we obtain final anomaly segmentations from
them. An overview of our final anomaly detection algorithm can be found in
Algorithm 1.

###  4.1 Learning the Normative Distribution using DDPMs

Similar to recent work [31, 23, 3], we use a DDPM [13] to learn the normative
distribution $q_{norm}(x)$ of normal image _slices_ $x\in\mathbb{R}^{H\times
W\times C}$ given a training set of healthy image slices
$\\{x^{(0)},\ldots,x^{(M)}\\}$.

We empirically observed that training with standard Gaussian noise transitions
given in Eq. 3 do not allow our final model to capture anomalies at different
scales. Instead, we train our DDPM with pyramidal Gaussian noise, which will
be introduced in Sec. 4.4.

###  4.2 Generating Anomaly Maps via Aggregation of Denoising Errors

After training, we use the learned denoising step along with the ground truth
backwards process to obtain an anomaly map. The denoising step can be thought
of as _normative diffusion_ as it is taking one step towards the normative
distribution, while the ground truth backwards transition is going back
towards the anomalous image. The discrepancy between the two is therefore
indicative of potential anomalies.

Specifically, for a potentially anomalous multi-modal input image $x_{0}$, we
first obtain the noised images $x_{t}$ for all $t\in[T_{l},T_{u}]$ using the
forward process (see Eq. 1). We use $T_{l}=75$, $T_{u}=200$ for all
experiments. Next, for each time step in that range, we obtain the KL-
divergence term (Eq. 4) which amounts to the Euclidean distance $d_{t}$
between the ground truth backwards transition
$\tilde{x}_{t}=\mu_{q}(x_{t},x_{0})$ and the learned denoising step
$\hat{x}_{t}=\mu_{\theta}(x_{t},t)$,

$d_{t}=\left(\mu_{q}(x_{t},x_{0})-\mu_{\theta}(x_{t},t)\right)^{2}.$ (6)

Importantly, each $d_{t}$ term can be computed independently. Thus, the entire
computation can be parallelized across time steps, which has the potential for
significant speedups compared to DDPM approaches performing pseudo-healthy
image generation typically requiring step-by-step denoising (_e.g_. [31, 3]).

In order, to obtain the initial multi-channel anomaly maps $a_{m}$, we
aggregate the individual $\\{d_{t}\\}$ with $t\in[T_{l},T_{u}]$ using some
aggregation function $A$. In this work, we explore a simple aggregation via
the arithmetic mean over time steps,

$A_{\text{am}}(\\{d_{t}\\})=\frac{1}{T_{u}-T_{l}}\sum_{t={T_{l}}}^{T_{u}}d_{t},$
(7)

as well as a novel aggregation mechanism based on the geometric mean which we
will describe next.

###  4.3 Geometric Mean for Error Aggregation

The arithmetic mean in Eq. 7 can be sensitive to outliers when aggregating
across timesteps. To address this, we propose to use the more robust
_geometric mean_ for aggregation:

$A_{\text{gm}}(\\{d_{t}\\})=\exp\left(\frac{1}{T_{u}-T_{l}}\sum_{t={T_{l}}}^{T_{u}}\log
d_{t}\right).$ (8)

This approach dampens the effect of extreme values across timesteps. We
empirically show that the geometric mean leads to better results for anomaly
detection compared to the arithmetic mean used in [23].

Figure 2: An illustration of pyramidal Gaussian noise.

###  4.4 Robust Anomaly Detection via Pyramidal Noise

In order to make the denoising transitions more sensitive to anomalies at
different scales (e.g. different lesion sizes), we propose to train our model
using a multi-resolution pyramidal Gaussian noise in which noise is sampled at
increasing image resolutions, bilinearly upsampled to the original image
resolution, and then added. This type of noise was previously explored for
improved image generation in the following blog post [30]. An illustration of
pyramidal noise is shown in Fig. 2.

Mathematically, we define pyramidal noise as

$\epsilon=\sum\limits_{i=1}^{10}c^{i}\cdot U(\epsilon^{(i)};H,W),$ (9)

where $c$ is a constant scaling factor ($c=0.8$ for all experiments),
$\epsilon^{(i)}\sim\mathcal{N}(0,I_{h_{i}\times w_{i}\times C})$ is
independent standard Gaussian noise of size $h_{i}\times w_{i}\times C$, and
$U$ is a bilinear upsampling operation that upscales the images to size
$H\times W\times C$. The Gaussian noise height and width $h_{i},w_{i}$ are
calculated by successively dividing the base width and height in half for each
level. Following [30], a slight jitter is added to these values at each level,
resulting in the following equations

$\displaystyle h_{i}=$ $\displaystyle\lceil H/(r_{i})^{i-1}\rceil$ (10)
$\displaystyle w_{i}=$ $\displaystyle\lceil W/(r_{i})^{i-1}\rceil$
$\displaystyle\text{with}\ \ r_{i}=2+\text{Unif}(0,2).$

Finally, we normalize the pyramidal noise to unit variance.

Similar to the multi-scale simplex noise previously proposed in AnoDDPM [31],
the pyramid noise allows the learned backwards transitions to denoise at
multiple spatial scales and thus increases robustness across different anomaly
sizes. However, in contrast to multi-scale simplex noise, pyramidal noise is
very simple to generate and incurs much smaller computational costs. In our
experiments, we empirically show that we obtain superior performance compared
to the AnoDDPM which uses simplex noise.

The pyramidal noise may be used for both the training (Sec. 4.1) and anomaly
map generation (Sec. 4.2). However, we found that training with pyramidal
noise, but generating anomaly maps with normal Gaussian noise consistently
improved performance (see also ablation study in Tab. 2). This can be
explained by the following intuition: During training, the denoising function
$\mu_{\theta}(x_{t},t)$ in Eq. 1 learns to remove pyramidal noise. During
evaluation, we only noise the images with Gaussian noise. This means that the
denoising function tries to remove the simple Gaussian noise plus the
additional lower-frequency components. Since anomalies may resemble lower-
frequency noise, the backwards transition is more likely to remove them. This
is similar to the intuition of the DAE approach in [14], where the denoising
autoencoder is trained with coarse noise and evaluated with no noise.

###  4.5 From Anomaly Scores to Segmentation Masks

Given a potentially anomalous multi-channel input volume $v$, we first slice
it into axial multi-channel slices. Then, for each slice, we obtain the multi-
channel anomaly maps, which are then stacked back into an anomaly volume
$a_{m}$. In order to reduce the channel-wise anomaly map to a pixel-wise
anomaly map $a$, prior work has used the arithmetic mean across the channel
dimension [14]. In this work, we additionally investigated a $max$ operation
across channels. We found that this led to slightly better performance for all
examined methods, including the DAE approach in [14], and thus used the $max$
operation for all experiments.

Following related work (_e.g_. [2]), we then apply a 3-dimensional median
filtering to the anomaly maps. We explore kernel sizes of $3\times 3\times 3$,
and $5\times 5\times 5$. We also apply median filtering to all baselines and
find that it improves all examined techniques.

In order to obtain a binarized segmentation mask or a decision about the
abnormality of a sample, thresholding must be performed. In reality, it is
non-trivial to find a good threshold without making use of any labeled data.
Many prior UAD works indeed do use labels to obtain the binarized segmentation
masks [23, 3, 14, 8, 2, 5]. Additionally, the ceil Dice is a widely used
metric in UAD. This metric is derived by creating a segmentation mask $s$
using an optimal threshold that is determined based on the test data, and then
comparing this mask to the ground truth segmentation. While this approach
provides an upper limit of performance, it does not reflect the true
unsupervised performance in scenarios where labeled data for tuning is
unavailable.

As an additional domain-agnostic thresholding strategy, we propose to
automatically threshold the anomaly maps using the Yen technique [32]. This
approach does not require access to the test dataset and allows to obtain a
unique threshold for every subject. In that sense, it allows evaluation in a
truly unsupervised setting with no knowledge of the target application. We
observed that Yen thresholding leads to consistent underestimation of the
anomaly size. Therefore, we additionally post-process the resulting
segmentation with a 3D dilation step. We additionally report the Dice score
obtained in this manner ($\text{Dice}_{\text{Yen}}$) for all methods examined
in this paper.

Figure 3: Comparative analysis of inferred anomalies and binarized
segmentations using Yen thresholding[32] across baseline methods and ANDi.
These anomaly maps, enhanced through modality max-pooling and median filtering
(kernel size 5), highlight the limitations of calibration-dependent FLAIR
thresholding[20] and DAE’s[14] inefficiency with small lesions. Among methods
based on diffusion models, both AnoDDPM[31] and ANDi demonstrate proficiency
in isolating anomalies of varying sizes, with ANDi showing superior
performance in unsupervised thresholding. The anomaly map brightness is
volume-scaled, potentially appearing darker in individual slices.

##  5 Experiments & Results

###  5.1 Data

We evaluated ANDi on data from three multi-modal brain MR datasets. Namely, we
used images from the brain tumor segmentation (BraTS’21) challenge [10] as
well as the MSSEG [9] and Ljubljana MS lesion datasets, which are both part of
the 2021 Shifts challenge [19]. Each subject in these datasets comprises the
following four MR scans: native T1-weighted, post-contrast T1 gadolinium
weighted (T1Gd), T2-weighted (T2), and Fluid Attenuated Inversion Recovery
(FLAIR).

###  5.2 Preprocessing

All datasets have been skull-stripped and resampled to a unified resolution.
We removed one subject from the Ljubljana dataset where the skull-stripping
has failed. The volumes and segmentation masks from the Shifts data were
coregistered to the SRI24 atlas [25]. After registration, the segmentation
masks were binarized using a threshold of 0.5. Afterward, we applied pairwise
histogram matching [21] between all subjects from the Shifts datasets to a
randomly chosen subject of the BraTS’21 dataset in order to reduce domain
shifts resulting from different acquisition protocols.

All volumes were individually divided by the 99th percentile foreground voxel
intensity for each subject and modality. The axial slices were downsampled to
a size of 128x128 using bilinear downsampling with anti-aliasing. The ground
truth segmentation masks were downsampled to the same size using nearest
neighbor interpolation.

###  5.3 Experiments

Similar to [14], we used the healthy slices in the BraTS’21 dataset to train
all our models. The BraTS’21 dataset is comparatively large consisting of MR
images of 1251 subjects, and contains considerable diversity in contrasts.
Using the data splits previously proposed by [14], we divided the dataset into
938 subjects for the training, 62 subjects for the validation, and 251
subjects for the testing. In the training split, we then identified all slices
that did not contain any tumors using the ground truth segmentations. The
training slices that do contain tumor tissue are discarded. The healthy
BraTS’21 training slices were then used to train our proposed as well as all
baseline methods.

We evaluated the different UAD algorithms examined in this study on all of the
test subjects of the BraTS’21 dataset, the 52 subjects of the MSSEG dataset,
and the 24 subjects of the Ljubljana dataset.

Table 1: Anomaly detection performance measured in AUPRC ($\uparrow$), $\lceil\text{Dice}\rceil$ ($\uparrow$) and $\text{Dice}_{\text{Yen}}$ ($\uparrow$) on the test dataset of BraTS’21 (251 subjects) and Shifts Challenge (MSSEG – 52 subjects, Ljubljana – 24 subjects). We highlight the first, second, and third best-performing unsupervised methods for each metric. We exclude the semi-supervised thresholding technique from the ranking as it requires labeled data for segmentation threshold tuning, and relies on an a priori modality selection. |  | BraTS21 | Shifts – MSSEG | Shifts – Ljubljana
---|---|---|---|---
Method | AUPRC | $\lceil\text{Dice}\rceil$ | $\text{Dice}_{\text{Yen}}$ | AUPRC | $\lceil\text{Dice}\rceil$ | $\text{Dice}_{\text{Yen}}$ | AUPRC | $\lceil\text{Dice}\rceil$ | $\text{Dice}_{\text{Yen}}$
Semi-sup. Thresholding (FLAIR) [20] |  | 0.741 | 0.674 | 0.106 | 0.430 | 0.331 | 0.018 | 0.425 | 0.345 | 0.027
w/ MF size 3 |  | 0.795 | 0.712 | 0.106 | 0.465 | 0.303 | 0.018 | 0.403 | 0.302 | 0.027
w/ MF size 5 |  | 0.818 | 0.726 | 0.106 | 0.386 | 0.231 | 0.018 | 0.304 | 0.226 | 0.023
Unsupervised DAE [14] |  | 0.810 | 0.724 | 0.700 | 0.022 | 0.117 | 0.037 | 0.138 | 0.140 | 0.077
\+ Modality max pooling |  | 0.817 | 0.733 | 0.699 | 0.049 | 0.126 | 0.043 | 0.059 | 0.100 | 0.041
w/ MF size 3 |  | 0.835 | 0.752 | 0.740 | 0.048 | 0.126 | 0.054 | 0.066 | 0.099 | 0.069
w/ MF size 5 |  | 0.845 | 0.763 | 0.759 | 0.041 | 0.120 | 0.058 | 0.070 | 0.096 | 0.070
AnoDDPM[31] |  | 0.600 | 0.532 | 0.488 | 0.104 | 0.140 | 0.150 | 0.070 | 0.105 | 0.062
w/ MF size 3 |  | 0.731 | 0.650 | 0.609 | 0.187 | 0.194 | 0.214 | 0.137 | 0.139 | 0.128
w/ MF size 5 |  | 0.769 | 0.687 | 0.657 | 0.156 | 0.169 | 0.167 | 0.152 | 0.131 | 0.129
AutoDDPM [3] |  | 0.543 | 0.511 | 0.437 | 0.084 | 0.144 | 0.152 | 0.103 | 0.192 | 0.133
w/ MF size 3 |  | 0.686 | 0.630 | 0.558 | 0.134 | 0.178 | 0.202 | 0.128 | 0.210 | 0.213
w/ MF size 5 |  | 0.731 | 0.674 | 0.606 | 0.108 | 0.127 | 0.150 | 0.076 | 0.155 | 0.155
ANDi (ours) |  | 0.654 | 0.587 | 0.388 | 0.219 | 0.274 | 0.290 | 0.325 | 0.331 | 0.228
w/ MF size 3 |  | 0.829 | 0.747 | 0.701 | 0.321 | 0.311 | 0.380 | 0.422 | 0.347 | 0.356
w/ MF size 5 |  | 0.858 | 0.781 | 0.766 | 0.174 | 0.201 | 0.240 | 0.293 | 0.242 | 0.253

###  5.4 Metrics

For all methods, we report the following three metrics: Area under the
precision-recall curve (AUPRC), ceil Dice ($\lceil\text{Dice}\rceil$), and the
fully unsupervised $\text{Dice}_{\text{Yen}}$ values, where no prior domain
knowledge is used to tune the segmentation threshold (see Sec. 4.5). To
calculate the AUPRC, the anomaly scores are normalized to be in the range of
$[0,1]$. AUPRC and $\text{Dice}_{\text{Yen}}$ evaluate the performance of the
method without assumptions about the target domain, while
$\lceil\text{Dice}\rceil$ uses the test dataset to determine an optimal
threshold and gives an upper-bound on the achievable performance.

###  5.5 Baselines

In addition to our proposed ANDi approach, we also evaluate 3 state-of-the-art
UAD baselines: DAEs [14], AnoDDPM [31], and AutoDDPM [3]. As an additional
baseline, we consider the semi-supervised thresholding technique by Meissen
_et al_. [20] where we used the FLAIR modality.

Out of the UAD baselines, only the DAE [14] was specifically designed for
multi-modal data. The original method uses the arithmetic mean to combine the
anomaly maps for the different input modalities. We additionally report
results with our proposed max modality pooling (see Sec. 4.5), showing that it
improves results in 2 out of 3 test datasets. Since the other methods were all
proposed for a single modality scenario, we by default extend them by the max
modality pooling operation to make them multi-modality compatible. AutoDDPM
[3] requires a masking threshold for their first mask calculation that was
obtained using the 95th percentile of the residuals. In our experiments we
found that the 99th percentile worked better on the validation dataset of
BraTS’21, and we used this threshold for all our analyses.

Lastly, we report all results with median filtering with a $3\times 3\times 3$
(MF size 3) and a $5\times 5\times 5$ (MF size 5) kernel, as we found the
results to be particularly sensitive to this hyperparameter choice.

###  5.6 Implementation Details

In line with prior work, we used a U-Net to implement our DDPM model [31, 3,
23, 13] as well as all DDPM-based baselines in this paper. We trained the
models with the AdamW optimizer with standard hyperparameters. The learning
rate is changed per step and follows a linear warmup with a cosine decay. A
description of all hyperparamters can be found in Appendix A.

Unlike usual DDPM implementations, we did not clip the predictions at each
timestep, as clipping may remove important outlier values from the anomaly
maps. We have found minor improvements when neglecting the clipping for a
Gaussian and pyramidal noise schedule.

##  6 Results

Our main results for all baselines on the three datasets are reported in Tab.
1 with qualitative results shown in Fig. 3. We find that our proposed ANDi
consistently and, on some datasets, substantially outperforms all unsupervised
baselines. Importantly, ANDi is the first UAD algorithm to surpass the strong
baseline of DAEs [14], which had never been outperformed on the BraTS’21
dataset. On the more challenging MS lesion datasets (MSSEG & Ljubliana), all
methods perform significantly worse. However, ANDi clearly outperforms the
baselines with and without median filtering. On the AUPRC metric, the best-
performing configuration of our method outperforms the next best method by 72%
on the MSSEG dataset and by 178% on the Ljubliana dataset. While the strongest
DAE [14] configuration performs very closely to our strongest configuration on
BraTS’21, it performs comparatively poorly on the Shifts datasets. This
suggests that the DAE approach may be specifically tailored to brain tumor
detection due to its use of coarse noise during training, which resembles
large anomalies.

For small lesions, we find that MF size 5 degrades the performance of all
unsupervised methods since the smoothing operation leads to small lesions
being discarded from the anomaly map. We also observe that the domain-agnostic
Yen thresholding achieves Dice scores similar to $\lceil\text{Dice}\rceil$
when median filtering is applied. This implies that the Yen thresholding
yields results close to the theoretical Dice upper bound when using the test
set to tune the threshold.

We observe that for BraTS’21 our method and the DAE approach [14] outperform
the semi-supervised thresholding baseline [20]. We emphasize that thresholding
is not directly comparable to the unsupervised methods because it can, by
definition, only detect hyperintense lesions, and requires labeled data to
determine a threshold for binary segmentation. The fact that it would not
perform well in a realistic scenario is reflected by the low
$\text{Dice}_{\text{Yen}}$ values.

In order to better understand the contribution of each of our proposed
novelties, we also performed an ablation study on the BraTS’21 validation
shown in Tab. 2. We started with the ANDi${}_{\text{base}}$ algorithm, which
is a DDPM trained with standard Gaussian noise, and anomaly map generation via
aggregation of KL-terms using the arithmetic mean as described in Sec. 4.2. We
then sequentially added the modality max pooling (Sec. 4.5), the geometric
mean aggregation (Sec. 4.3), training and testing with pyramidal noise (Sec.
4.4), median filtering with kernel size 5 (Sec. 4.5), and finally training
with pyramidal noise but evaluating with standard Gaussian noise (Sec. 4.4).
We find that each of these contributions is crucial for achieving the final
performance.

Table 2: Ablation study of our individual contributions on the validation dataset of BraTS’21 (62 subjects) measured in ($\uparrow$), $\lceil\text{Dice}\rceil$ ($\uparrow$) and $\text{Dice}_{\text{Yen}}$ ($\uparrow$). Each row includes the cumulative additions from the previous rows. |  BraTS21${}_{\text{val}}$
---|---
Method | AUPRC | $\lceil\text{Dice}\rceil$ | $\text{Dice}_{\text{Yen}}$
ANDi${}_{\text{base}}$ | 0.221 | 0.270 | 0.114
\+ Max pooling | 0.230 | 0.283 | 0.114
\+ Geometric mean | 0.268 | 0.319 | 0.123
\+ Pyramidal training | 0.674 | 0.609 | 0.385
\+ MF size 5 | 0.853 | 0.784 | 0.760
\+ Eval. with Gauss | 0.887 | 0.807 | 0.784

##  7 Discussion & Conclusion

In this work, we proposed a novel DDPM-based algorithm for unsupervised
anomaly detection in medical imaging data. Our proposed approach functions by
directly aggregating deviations of the learned denoising step from the ground
truth backwards transition. In contrast to recent related work, our method
does not require an additional in-painting or pseudo-healthy image generation
step, making it simple and fast to use. In order to achieve state-of-the-art
performance we introduced crucial methodological contributions. Most notably,
we use pyramidal noise during training of the DDPM, as well as geometric mean
for aggregation of the above deviations.

We further proposed an evaluation strategy based on Yen thresholding of
individual anomaly maps ($\text{Dice}_{\text{Yen}}$), which does not require
any knowledge of, or data from the target domain. While the established AUPRC
evaluation also does not require setting a specific threshold, the advantage
of the Yen thresholding is that it can be used to obtain a binary segmentation
without knowledge of the target domain.

When evaluating our approach in a clinically relevant multi-modal brain MR
image scenario on three challenging and diverse datasets (BraTS’21, MSSEG,
Ljubljana), we clearly outperform two recent diffusion model-based approaches
(AnoDDPM & AutoDDPM), as well as the strong DAE baseline proposed by Kascenas
_et al_. [14]. In particular, the detection of MS lesions detection is
substantially improved, exceeding the next best method by up to 178% in terms
of AUPRC.

Combined, the developments presented in this paper are a step in the direction
of achieving more clinical utility for diffusion model-based UAD. However, it
is also clear that while the work is promising, UAD methods have not yet
achieved the performance required for clinical utility. In future work, we aim
to investigate the interplay of DDPMs with different noise types, which
appears to be of great importance for improving the performance of our anomaly
detection method. While our method is currently trained on 2D slices, in
subsequent work, we aim to better incorporate the inherently three-dimensional
nature of brain MRI scans.

##  8 Acknowledgements

We thank Jakob H. Macke for his helpful comments. This work was funded by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy – EXC 2064/1, Project number 390727645, with
additional support from the DFG Emmy Noether 513851350 (TW) and the German
Federal Ministry of Education and Research (BMBF) via the Tübingen AI Center,
FKZ: 01IS18039A. We also utilized the Tübingen Machine Learning Cloud,
supported by DFG FKZ INST 37/1057-1 FUGG. The authors thank the International
Max Planck Research School for Intelligent Systems (IMPRS-IS) for their
support of Jaivardhan Kapoor, Thomas Wolfers, and Christian F. Baumgartner.

## References

  * Baur et al. [2020] Christoph Baur, Benedikt Wiestler, Shadi Albarqouni, and Nassir Navab.  Bayesian skip-autoencoders for unsupervised hyperintense anomaly detection in high resolution brain mri.  pages 1905–1909, 2020.
  * Baur et al. [2021] Christoph Baur, Stefan Denner, Benedikt Wiestler, Nassir Navab, and Shadi Albarqouni.  Autoencoders for unsupervised anomaly segmentation in brain mr images: a comparative study.  _Medical Image Analysis_ , 69:101952, 2021.
  * Bercea et al. [2023a] Cosmin I. Bercea, Michael Neumayr, Daniel Rueckert, and Julia A Schnabel.  Mask, stitch, and re-sample: Enhancing robustness and generalizability in anomaly detection through automatic diffusion models.  In _ICML 3rd Workshop on Interpretable Machine Learning in Healthcare (IMLH)_ , 2023a.
  * Bercea et al. [2023b] Cosmin I Bercea, Benedikt Wiestler, Daniel Rueckert, and Julia A Schnabel.  Generalizing unsupervised anomaly detection: Towards unbiased pathology screening.  In _Medical Imaging with Deep Learning_ , 2023b.
  * Bercea et al. [2023c] Cosmin I Bercea, Benedikt Wiestler, Daniel Rueckert, and Julia A Schnabel.  Reversing the abnormal: Pseudo-healthy generative networks for anomaly detection.  _arXiv preprint arXiv:2303.08452_ , 2023c.
  * Cabezas et al. [2011] Mariano Cabezas, Arnau Oliver, Xavier Lladó, Jordi Freixenet, and Meritxell Bach Cuadra.  A review of atlas-based segmentation for magnetic resonance brain images.  _Computer methods and programs in biomedicine_ , 104(3):e158–e177, 2011.
  * Chen and Konukoglu [2018] Xiaoran Chen and Ender Konukoglu.  Unsupervised detection of lesions in brain mri using constrained adversarial auto-encoders.  _arXiv preprint arXiv:1806.04972_ , 2018.
  * Chen et al. [2021] Xiaoran Chen, Nick Pawlowski, Ben Glocker, and Ender Konukoglu.  Normative ascent with local gaussians for unsupervised lesion detection.  _Medical Image Analysis_ , 74:102208, 2021.
  * Commowick et al. [2018] Olivier Commowick, Audrey Istace, Michael Kain, Baptiste Laurent, Florent Leray, Mathieu Simon, Sorina Camarasu Pop, Pascal Girard, Roxana Ameli, Jean-Christophe Ferré, et al.  Objective evaluation of multiple sclerosis lesion segmentation using a data management and processing infrastructure.  _Scientific reports_ , 8(1):13650, 2018.
  * et. al. [2021] Ujjwal Baid et. al.  The rsna-asnr-miccai brats 2021 benchmark on brain tumor segmentation and radiogenomic classification.  _arXiv preprint arXiv: 2107.02314_ , 2021.
  * Goodfellow et al. [2014] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  _Advances in neural information processing systems_ , 27, 2014.
  * Hernandez Petzsche et al. [2022] Moritz R Hernandez Petzsche, Ezequiel de la Rosa, Uta Hanning, Roland Wiest, Waldo Valenzuela, Mauricio Reyes, Maria Meyer, Sook-Lei Liew, Florian Kofler, Ivan Ezhov, et al.  Isles 2022: A multi-center magnetic resonance imaging stroke lesion segmentation dataset.  _Scientific data_ , 9(1):762, 2022.
  * Ho et al. [2020] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  In _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020.
  * Kascenas et al. [2022] Antanas Kascenas, Nicolas Pugeault, and Alison Q O’Neil.  Denoising autoencoders for unsupervised anomaly detection in brain mri.  In _International Conference on Medical Imaging with Deep Learning_ , pages 653–664. PMLR, 2022.
  * Kascenas et al. [2023] Antanas Kascenas, Pedro Sanchez, Patrick Schrempf, Chaoyang Wang, William Clackett, Shadia S. Mikhael, Jeremy P. Voisey, Keith Goatman, Alexander Weir, Nicolas Pugeault, Sotirios A. Tsaftaris, and Alison Q. O’Neil.  The role of noise in denoising models for anomaly detection in medical images.  _Medical Image Analysis_ , 90:102963, 2023.
  * Kingma and Welling [2013] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Liew et al. [2022] Sook-Lei Liew, Bethany P Lo, Miranda R Donnelly, Artemis Zavaliangos-Petropulu, Jessica N Jeong, Giuseppe Barisano, Alexandre Hutton, Julia P Simon, Julia M Juliano, Anisha Suri, et al.  A large, curated, open-source stroke neuroimaging dataset to improve lesion segmentation algorithms.  _Scientific data_ , 9(1):320, 2022.
  * Liu et al. [2014] Xiaoxiao Liu, Marc Niethammer, Roland Kwitt, Matthew McCormick, and Stephen Aylward.  Low-rank to the rescue–atlas-based analyses in the presence of pathologies.  In _Medical Image Computing and Computer-Assisted Intervention–MICCAI 2014: 17th International Conference, Boston, MA, USA, September 14-18, 2014, Proceedings, Part III 17_ , pages 97–104. Springer, 2014.
  * Malinin et al. [2021] Andrey Malinin, Neil Band, Alexander Ganshin, German Chesnokov, Yarin Gal, Mark J. F. Gales, Alexey Noskov, Andrey Ploskonosov, Liudmila Prokhorenkova, Ivan Provilkov, Vatsal Raina, Vyas Raina, Denis Roginskiy, Mariya Shmatova, Panos Tigar, and Boris Yangel.  Shifts: A dataset of real distributional shift across multiple large-scale tasks.  _arXiv preprint arXiv:2107.07455_ , 2021.
  * Meissen et al. [2021] Felix Meissen, Georgios Kaissis, and Daniel Rueckert.  Challenging current semi-supervised anomaly segmentation methods for brain mri.  In _International MICCAI brainlesion workshop_ , pages 63–74. Springer, 2021.
  * Nyúl et al. [2000] László G Nyúl, Jayaram K Udupa, and Xuan Zhang.  New variants of a method of mri scale standardization.  _IEEE transactions on medical imaging_ , 19(2):143–150, 2000.
  * Perlin [2002] Ken Perlin.  Improving noise.  In _Proceedings of the 29th annual conference on Computer graphics and interactive techniques_ , pages 681–682, 2002.
  * Pinaya et al. [2022] W. H. Pinaya, M. Graham, R. Gray, P. F. D. Costa, Petru-Daniel Tudosiu, P. Wright, Y. Mah, A. Mackinnon, J. Teo, H. Jäger, D. Werring, Geraint Rees, P. Nachev, S. Ourselin, and M. Cardoso.  Fast unsupervised brain anomaly detection and segmentation with diffusion models.  _MICCAI_ , 2022.
  * Prastawa et al. [2004] Marcel Prastawa, Elizabeth Bullitt, Sean Ho, and Guido Gerig.  A brain tumor segmentation framework based on outlier detection.  _Medical image analysis_ , 8(3):275–283, 2004.
  * Rohlfing et al. [2010] Torsten Rohlfing, Natalie M Zahr, Edith V Sullivan, and Adolf Pfefferbaum.  The sri24 multichannel atlas of normal adult human brain structure.  _Human brain mapping_ , 31(5):798–819, 2010.
  * Sato et al. [2018] Daisuke Sato, Shouhei Hanaoka, Yukihiro Nomura, Tomomi Takenaga, Soichiro Miki, Takeharu Yoshikawa, Naoto Hayashi, and Osamu Abe.  A primitive study on unsupervised anomaly detection with an autoencoder in emergency head ct volumes.  In _Medical Imaging 2018: Computer-Aided Diagnosis_ , pages 388–393. SPIE, 2018.
  * Schlegl et al. [2017] Thomas Schlegl, Philipp Seeböck, Sebastian M Waldstein, Ursula Schmidt-Erfurth, and Georg Langs.  Unsupervised anomaly detection with generative adversarial networks to guide marker discovery.  In _International conference on information processing in medical imaging_ , pages 146–157. Springer, 2017.
  * Schlegl et al. [2019] Thomas Schlegl, Philipp Seeböck, Sebastian M Waldstein, Georg Langs, and Ursula Schmidt-Erfurth.  f-anogan: Fast unsupervised anomaly detection with generative adversarial networks.  _Medical image analysis_ , 54:30–44, 2019.
  * Van Leemput et al. [2001] Koen Van Leemput, Frederik Maes, Dirk Vandermeulen, Alan Colchester, and Paul Suetens.  Automated segmentation of multiple sclerosis lesions by model outlier detection.  _IEEE transactions on medical imaging_ , 20(8):677–688, 2001.
  * Whitaker [2023] Jonathan Whitaker.  Multi-Resolution-Noise-for-Diffusion-Model-Training.  https://wandb.ai/johnowhitaker/multires_noise/reports/Multi-Resolution-Noise-for-Diffusion-Model-Training--VmlldzozNjYyOTU2, 2023.  Accessed: 2023-4-21.
  * Wyatt et al. [2022] Julian Wyatt, Adam Leach, Sebastian M Schmon, and Chris G Willcocks.  Anoddpm: Anomaly detection with denoising diffusion probabilistic models using simplex noise.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 650–656, 2022.
  * Yen et al. [1995] Jui-Cheng Yen, Fu-Juay Chang, and Shyang Chang.  A new criterion for automatic multilevel thresholding.  _IEEE Transactions on Image Processing_ , 4(3):370–378, 1995.

Supplementary Material

Figure 4: Wallclock time for processing 1 volume ($128\times 128\times 155$)
for the unsupervised baselines and ANDi. We report the mean $\pm$ standard
deviation over 10 runs. ANDi is substantially faster than the other DDPM-based
methods, and is only outperformed by the simpler DAE.

##  Appendix A Training and Evaluation Hyperparameters

All DDPM models were parameterized by the same U-Net architecture. ANDi was
trained with pyramidal noise, AutoDDPM with Gaussian noise, and AnoDDPM with
simplex noise. All hyperparameters for the simplex noise are equal to the
original publication with a starting frequency of $\nu=2^{-6}$ , octave of
$N=6$ and a decay of $\gamma=0.8$. When generating the simplex noise, the seed
was changed before every calculation, and a slice $s$ from the 3-dimensional
noise function was extracted. We trained ANDi with a total number of diffusion
timesteps $T=1000$, a linear noise schedule between
$\beta_{1}=10^{-4},\beta_{T}=0.02$ and a batch size of $128$ over $232$
epochs. Starting with a learning rate of $2\cdot 10^{-5}$, we linearly
increased it to $1\cdot 10^{-4}$ over 5% of the total training steps and then
decayed it back to $2\cdot 10^{-5}$ with a cosine schedule.

All diffusion model-based methods require setting an upper ($T_{u}$) and a
lower bound ($T_{l}$) for the denoising operations during evaluation. For ANDi
we empirically determined $T_{u}=200$ to $T_{l}=75$ to work well, though we
show that ANDi is not sensitive to the exact choice of time range in Appendix
C. For AutoDDPM and AnoDDPM we followed the original publications when setting
the time range. AutoDDPM and AnoDDPM require $T_{l}=0$ as they need denoising
back to the image space fully. The initial mask generation step of AutoDDPM
uses $T_{u}=200$ and AnoDDPM uses $T_{u}=250$. AutoDDPM requires an additional
thresholding hyperparameter during evaluation to obtain the binary mask for
in-painting. In contrast to the original publication, this threshold was set
to the 99th percentile of the residuals obtained on the validation dataset of
BraTS’21. This choice improved the results for AutoDDPM. The subsequent
inpainting process was run for the last $50$ denoising steps with $5$
resamplings per step.

Figure 5: Influence of the chosen time step interval ($T_{u}-T_{l}$) on the
performance measured in AUPRC for ANDi and the Diffusion Model baselines. ANDi
is the most robust method regarding this hyperparameter choice and reaches
optimal performance with the smallest amount of denoising steps.

##  Appendix B Wallclock Analysis of Methods

In order to better understand the computational time required for anomaly
detection using the different methods, we performed a wallclock analysis for
one volume.

Fig. 4 shows the wallclock time for processing the anomaly maps for each of
the unsupervised baselines, as well as ANDi. This analysis includes
postprocessing steps over the anomaly maps, such as maxpool over modalities,
median filtering, and thresholding. These steps are common across all methods
and, therefore, contribute a constant offset $1.08\pm 0.02s$ to the wallclock
time. We used an NVIDIA-V100 GPU on a 16-core machine for all evaluations.

ANDi is the fastest among the DDPM-based methods and can process a volume
almost twice as fast as the next best method in that category, AnoDDPM. ANDi
is only outperformed by the DAE baseline, which predicts the anomaly map in a
single pass of the UNet.

The speedups of ANDi with respect to the other DDPM-based methods can largely
be explained due to the number of denoising steps considered in the
evaluation. ANDi requires fewer timesteps than AnoDDPM and AutoDDPM to achieve
maximum performance (see also Fig. 5). Further, the gains with respect to
AutoDDPM can be explained by the fact that AutoDDPM requires two stages for
anomaly detection, an initial anomaly mask creation and a subsequent
inpainting stage, while ANDi requires only a single stage.

In our evaluations of AnoDDPM, in order to reduce the computational
complexity, we reused the same simplex noise sample for each slice of a test
volume, but augmented it with random rotations and flipping. Computing
independent simplex noise for each slice of a volume would take significantly
longer due to the need to sample 3D simplex noise to reduce corner artifacts.
The noise samples are still independent for each step of the diffusion
process. If implemented with slice-wise independent simplex noise, AnoDDPM
would be substantially slower, and would likely fall behind AutoDDPM in terms
of execution time.

We note that ANDi has the potential to be parallelized over denoising
timesteps, given adequate GPU memory and FLOPs. While we didn’t explore this
in the present paper, parallelization over denoising steps would lead to
further substantial speedups. In contrast, AnoDDPM and AutoDDPM are sequential
by nature and hence cannot be parallelized in this way.

##  Appendix C Sensitivity Analysis on Number of Denoising Timesteps

As mentioned earlier, all DDPM-based algorithms required setting an upper
($T_{u}$) and a lower bound ($T_{l}$) for the denoising operations during
evaluation. The resulting time range $T_{u}-T_{l}$ is an important
hyperparameter for all methods.

AnoDDPM and AutoDDPM require setting $T_{l}=0$ because they need to denoise
back to the input space. ANDi, on the other hand, is not required to reach the
end of the diffusion chain because no pseudo-healthy image is not calculated.
Therefore, we can choose any $T_{l}$ as an endpoint for the denoising steps.

In order to better understand the DDPM-based algorithms’ sensitivity to the
time range $T_{u}-T_{l}$, we performed a parameter exploration on the
validation dataset of BraTS’21. The results are shown in Fig. 5. We observe
that ANDi is significantly more robust to the choice of $T_{u}$ and $T_{l}$.
Moreover, we note that all ANDi configurations within our exploration range
outperform the best AnoDDPM and AutoDDPM configurations. It can further be
observed that ANDi obtains optimal anomaly detection with a smaller number of
total denoising steps compared to the other two methods. This partially
explains the speedups we measured in Appendix B.

##  Appendix D Additional Randomly Sampled Qualitative Results

In the following we present additional qualitative anomaly detection results
for all methods on the BraTS’21 (Fig. 6), MSSEG (Fig. 7), and Ljubljana
datasets (Fig. 8). As in the main article, we show the anomaly maps and the
binarized segmentations obtained using Yen thresholding. The shown images were
randomly sampled from the test datasets.

The results shown use median filtering with a filter size of $5\times 5\times
5$ for the BraTS21 dataset and $3\times 3\times 3$ for the two Shifts
datasets. As mentioned in the main text, a larger kernel size hurts
performance for small anomalies, such as MS lesions in the Shifts dataset.

For MS lesions, the DAE and the Thresholding approaches often detect the
complete slice as anomalous when using Yen thresholding. This happens when the
volume is segmented into foreground and background due to the anomaly maps not
being informative enough to find an appropriate threshold.

Figure 6: Randomly selected qualitative results for anomaly maps and Yen
thresholding on the BraTS’21 test dataset. The anomaly maps are post-processed
using median filtering (kernel size 5). Figure 7: Randomly selected
qualitative results for anomaly maps and Yen thresholding on the Shifts-MSSEG
test dataset. The anomaly maps are post-processed using median filtering
(kernel size 3). Figure 8: Randomly selected qualitative results for anomaly
maps and Yen thresholding on the Shifts-Ljubljana test dataset. The anomaly
maps are post-processed using median filtering (kernel size 3).