# An Empirical Evaluation of Time-Series Feature Sets

1st Trent Henderson  School of Physics
The University of Sydney
Sydney, Australia
<EMAIL_ADDRESS>2nd Ben D. Fulcher  School of Physics
The University of Sydney
Sydney, Australia
<EMAIL_ADDRESS>
###### Abstract

Solving time-series problems using informative features has been rising in
popularity due to the availability of numerous software packages for time-
series feature extraction. Feature-based time-series analysis can now be
performed using any one of a range of time-series feature sets, including
hctsa (7730 features: Matlab), feasts (42 features: R), tsfeatures (63
features: R), Kats (40 features: Python), tsfresh (up to 1558 features:
Python), TSFEL (390 features: Python), and the C-coded catch22 (22 features,
able to be run from Matlab, R, Python, and Julia). There is substantial
overlap in the types of time-series analysis methods included in these feature
sets (including properties of the autocorrelation function and Fourier power
spectrum, and distributional shape statistics), but they are yet to be
systematically compared. Here we compare these seven feature sets on their
computational speed, assess the redundancy of features contained in each set,
and evaluate the overlap and redundancy across different feature sets. We take
an empirical approach to measuring feature similarity, based on the similarity
of their outputs across a diverse set of real-world and model-simulated time
series. We find that feature sets vary across approximately three orders of
magnitude in their computation time per feature on a laptop for a 1000-sample
time series, from the fastest feature sets catch22 and TSFEL ($\sim 0.1$ ms
per feature) to tsfeatures ($\sim 3$ s per feature). Using PCA to evaluate
feature redundancy within each set, we find the highest within-set redundancy
for TSFEL and tsfresh. For example, in TSFEL, 90% of the variance across 390
features can be captured with just four principal components. Finally, we
introduce a metric for quantifying overlap between pairs of feature sets,
which indicates substantial overlap between the feature sets. We found that
the largest feature set, hctsa, is the most comprehensive, and that tsfresh is
the most distinctive, due to its incorporation of large numbers of Fourier
coefficients that are summarized at higher levels in the other sets. Our
results provide empirical understanding of the differences between existing
feature sets, information that can be used to better understand and tailor
feature sets to their applications.

###### Index Terms:

time-series analysis, time-series features.

##  I Introduction

Time series, or ordered measurements taken over time, are a ubiquitous form of
data in the sciences and industry. Applications of time-series analysis are
correspondingly diverse, from forecasting (the spread of disease through a
population or future rainfall) to classification (from inferring cord pH from
fetal heart-rate dynamics [1] to detecting abnormalities in silicon
semiconductors [2]). As a result, there is a wide range of quantitative tools
for analyzing time series.

In this work, we focus on methods for extracting informative summary
statistics, or features, from univariate time series (sampled uniformly in
time). Such methods are diverse, and capture properties of the time-series
distribution of values, autocorrelation structure, entropy, model-fit
statistics, nonlinear time-series analysis, stationarity, and many others [4,
3]. The first comprehensive effort to organize the interdisciplinary time-
series analysis literature, using an approach termed highly comparative time-
series analysis, encoded thousands of diverse analysis methods as features and
compared their behavior on a wide range of time-series data [3]. This
culminated in a software framework that includes implementations of over 7 000
time-series features, as the Matlab package hctsa [4, 5]. While hctsa is
comprehensive in its coverage of time-series analysis methods, it involves a
computational cost exceeding what is required of many real-world applications,
for which a smaller number of simpler features is often sufficient.
Additionally, hctsa requires access to the proprietary Matlab software, which
limits broader use. To address this limitation, catch22 is a reduced set of 22
canonical time-series features, that were derived from hctsa via a procedure
designed to both maximize classification performance across 93 tasks and
reduce redundancy among the top-performing features [6]. These 22 features
capture the main conceptual areas of time-series analysis covered in hctsa,
including properties of distributional shape, linear and nonlinear
autocorrelation, long-range scaling, symbolic transition rates, and others.
The catch22 feature set is coded in C, with wrappers for Python and Matlab,
and derived native packages in R (Rcatch22 [7]) and Julia (Catch22.jl [8]).

Other independently developed packages for computing time-series features also
exist in R, including feasts [9] and tsfeatures [10]. feasts is an updated and
refined iteration of tsfeatures, and both include sets of features tailored
largely to time series commonly encountered in econometrics, such as
seasonality and autoregressive processes. Examples of features included in
these sets include those associated with generalized autoregressive
conditional heteroscedasticity (GARCH) models, crossing points, seasonality,
and Seasonal and Trend decomposition using Loess (STL).

In Python, three popular libraries for time-series feature extraction are
tsfresh [11], TSFEL [12], and Kats [13]. tsfresh is a popular feature set that
contains implementations of up to 1 558 features, including autocorrelation
and Fourier decomposition methods, entropy, and distributional properties.
TSFEL includes features associated with autocorrelation properties and Fourier
transforms, spectral quantities and wavelet decompositions, and distributional
characteristics. The feature set included in Kats, developed by Facebook’s
Engineering team, was designed to be similar in composition to R’s existing
tsfeatures package and includes features associated with crossing points, STL
decomposition, sliding windows, autocorrelation properties, and Holt-Winters
algorithms.

There are major barriers to leveraging these various feature-extraction
packages, as they span programming languages and differ in required data
formats and syntax. For example, Kats requires a DateTime formatted time index
and a corresponding numeric vector of values, whereas other packages represent
uniformly sampled time-series data as an ordered vector of values. The R
package, theft: Tools for Handling Extraction of Features from Time series
[15], addresses these difficulties, providing a standardized computational
framework for time-series feature extraction, supporting the catch22, feasts,
tsfeatures, tsfresh, TSFEL, and Kats feature sets.

While the time-series analysis community now has ready access to multiple
time-series feature sets, how does one select an appropriate feature set to
use for a given application? Does the behavior of features in different
feature sets largely overlap, or are some distinctive? Different applications
may require different trade-offs between feature comprehensiveness and
computation speed—how do the feature-sets compare? While some sets contain
many features, are these features largely capturing independent properties of
the data, or are they highly correlated (and hence redundant), such that care
will need to be taken to statistically account for the non-independence (or
dimensionality reduction may be particularly useful). Addressing these
important questions requires systematic comparison of the behavior of the
existing tools for time-series feature extraction, which has not previously
been done. In this work we aim to address this shortcoming by comparing seven
time-series feature sets: hctsa (v1.06), catch22 (v0.1.12), feasts (v0.2.1),
tsfeatures (v1.0.2), tsfresh (v0.18.0), TSFEL (v0.1.4), and Kats (initial
release). We assess these feature sets in three different settings: (i) their
computation time and its scaling with time-series length; (ii) the redundancy
of features within a given feature set; and (iii) overlap and redundancy in
feature behavior between feature sets. We take a novel empirical approach of
assessing features using their behavior across a diverse set of real-world and
model-generated time series. Note that tsfeatures includes multiple feature
modules: we included all features here.

##  II Computation Time

In applications involving large datasets, where computational resources are
limited, smaller feature sets that are fast to compute are advantageous. But
for applications in which accuracy is prioritized, larger sets of complex
time-series features may be slower to compute but can yield higher
performance. For example, some feature sets, like hctsa, contain complex
nonlinear time-series analysis methods, and tsfeatures includes GARCH
modeling. These features are slower to compute than simple properties of the
distribution (e.g., mean and variance) or linear correlation statistics (e.g.,
Fourier coefficients or linear autocorrelations). It is thus important to
empirically benchmark the computation-time scaling of different time-series
feature sets with time-series length.

Here we estimated the computation-time scaling of each feature set across the
length range of 100–1000 samples (indicative of time-series lengths in many
applied problems). We computed the time taken to evaluate each of the seven
feature sets on a set of noise time series generated as independent samples
from a Gaussian distribution, across five lengths, $T=100,250,500,750$, and
$1000$. The median and interquartile range of computation time was returned
after repeating the process for ten independent noise time series of each
length. We confirmed that the results are not specific to our choice to use
white noise time series for benchmarking; we found quantitatively similar
results when repeating the analysis with noisy sinusoids (generated from:
$x_{t}=2\sin(2t)+\epsilon$, for i.i.d. $\epsilon\sim\mathcal{N}(0,1)$ and
$0\leq t\leq 3\pi$). Calculations were run on a 2019 MacBook Pro with an Intel
Core i7 2.6 GHz 6-Core CPU, and all feature sets were timed in their native
programming language (hctsa was run in parallel mode, allowing access to all
six CPU cores).

Feature sets varied considerably in the number of features they contain, from
catch22 (22 features), feasts (43 features), tsfeatures (62 features), Kats
(40 features), TSFEL (185 features for time series of length $T=100$, 260 for
length $T=250$, 285 for length $T=500$, and 390 for lengths $T>500$; feature
numbers scale approximately linearly with time-series length until $T>500$),
tsfresh (1 558 features), and hctsa (7 730 features). All features computed
successfully for all lengths for catch22, feasts, tsfeatures, and TSFEL. For
tsfresh, 392 features ($25.2$%) did not compute successfully for $T=100$, but
no issues were noted for any other time-series lengths. These features all
corresponded to fast Fourier transform (FFT) coefficients for real, imaginary,
absolute, and angle components of 51–99. For Kats, 3 features ($7.5$%) did not
compute successfully across all time-series lengths. These three features
related to level, trend, and seasonal parameters of a Holt–Winters model fit,
which understandably would not fit successfully on a white noise time series.
For hctsa, 706 features ($0.91$%) did not compute successfully for $T=100$,
469 features ($0.61$%) for $T=250$, 458 features ($0.59$%) for $T=500$, 470
features ($0.61$%) for $T=750$, and 443 features ($0.57$%) for $T=1000$. Some
features (such as stationarity on a certain fixed timescale) will not compute
on short time series, and some TISEAN-based features [14], requiring
compilation from C and Fortran code, were excluded from this analysis.

Median and interquartile ranges of computation times for each feature set are
plotted as a function of time-series length in Fig. 1A. We found that feature-
set computation times varied over approximately four orders of magnitude. For
example, for a 1000-sample time series, feature-extraction times varied from
$<10$ ms for the 22 C-coded features in catch22, to $16.5$ s for the $7300$
hctsa features. The other feature sets showed intermediate computation times,
from TSFEL ($0.03$ s) and Kats ($0.06$ s), to feasts (0.47 s), tsfresh (2.53
s), and tsfeatures (6.18 s). To more fairly compare feature sets of different
sizes (from the 22 features in catch22 to the $>7000$ features in hctsa), we
next investigated the computation time per successfully computed feature,
shown in Fig. 1B. We find that computation time per feature also varied over
several orders of magnitude. For a 1000-sample time series, feature
computation times were fastest for TSFEL ($<0.1$ ms per feature) and catch22
(0.3 ms per feature); intermediate for Kats (2 ms per feature), hctsa (2 ms
per feature), tsfresh (2 ms per feature), and feasts (11 ms per feature); and
slowest for tsfeatures (280 ms per feature).

Figure 1:  Computation times of different feature sets vary across orders of
magnitude. A Median and interquartile range of computation time across ten
Gaussian noise time series is plotted as a function of time-series length on
logarithmic axes for lengths $T=100,250,500,750$, and $1000$. B Median
computation time per successfully computed feature is plotted for the same
data as in A. Results are shown for catch22 (dark green), feasts (orange),
hctsa (purple), Kats (pink), tsfeatures (light green), TSFEL (yellow), and
tsfresh (brown).

##  III Evaluation Dataset

Having evaluated the differences in computation time of different feature
sets, we next aimed to compare their behavior on real data, in order to assess
the similarity of the constituent features, both within and between feature
sets. To do this, we use the ‘Empirical 1000’ dataset (version 10) [16], which
consists of a diverse set of 1000 real-world and model-generated time series
[17]. Model-generated data spans a broad range of generative processes,
including chaotic and non-chaotic deterministic dynamical systems (both
discrete and continuous), and a range of stochastic processes. The real-world
time series span a broad range of domains, including hydrology, meteorology,
finance, astrophysics, medicine, and animal-sound audio. While the Empirical
1000 time series range in length from $T=208$ to $T=10\,000$ samples, we
analyzed the first 1000 samples for time series longer than 1000 points to
minimize the computational expense. To demonstrate the range of temporal
patterns contained in the Empirical 1000 dataset, six example time series
(after truncation) are plotted in Fig. 2. This diversity of dynamical patterns
reflects the corresponding diversity in the types of real-world and generative
models studied in science and industry [17].

Given the nature of some of the time series in the Empirical 1000, some
individual features failed to compute successfully. To minimize issues with
further analysis, we employed a filtering process to retain the largest subset
of features possible for each feature set while also removing consistently
problematic time series from the Empirical 1000. First, all features were
first normalized independently using a $z$-score transformation. We then
filtered each set to only the $z$-score transformed features that had $<10$%
missing values. This removed no features for catch22 or TSFEL, 2 features
($4.7$%) for feasts (41 remaining), 3 features ($7.5$%) for Kats (37
remaining), 34 features ($4.4$%) for tsfresh (745 remaining—a maximum of 779
features were computed on the Empirical 1000 compared to the 1 558 features
computed for the Gaussian noise computation time analysis), 3 features
($4.8$%) for tsfeatures (59 remaining), and 22 features ($0.3$%) for hctsa (7
708 remaining). After filtering for problematic features, we then filtered out
any time series in the Empirical 1000 with missing values across any feature
in any set. This removed 105 time series, leaving 895 time series remaining.
Of the 105 removed problematic time series, 99 ($94.3$%) were generated due to
features in TSFEL, ($3.8$%) were generated due to tsfresh, 1 ($0.6$%) was
generated due to feasts, and 1 ($0.6$%) was generated due to tsfeatures. We
use these 895 time series and the retained features for each set to
characterize feature behaviour throughout the remainder of this paper.

Figure 2:  We assess time-series features based on their behavior on a dataset
of 1000 diverse real-world and model-generated time series, a selection of
which are plotted here [16, 17]. Six example time series from the Empirical
1000 dataset [16, 17] are plotted (after truncation to the first 1000 samples,
as described in text): A audio of the Black-bellied Plover
(AS_s2.1_f4_b8_l8280_135971), B the Trans Polar Index normalized pressure
difference between Hobart and Stanley (CM_tpi_hob), C a simulation from the
dissipative forced Brusselator ODE [18] (FL_fbruss_L300_N10000_IC_0.3_2_y), D
a gait time series from the Parkinson’s disease gait dataset [19, 20]
(MD_gaitpdb_GaCo05_01l_SNIP_12562955), E river flow (ME_winter_rf), and F a
simulation from the dissipative predator-prey map [18]
(MP_predprey_L5000_IC_0.5_0.55_x). Time-series names are provided in
parentheses for reference to their corresponding data entries on www.comp-
engine.org [17].

##  IV Within Feature-Set Redundancy

In constructing a feature set, choices need to be made in selecting the range
of conceptually distinct types of time-series analysis algorithms to include,
their parametric variations, and how to represent the method’s outputs as one
or more real-valued features. For example, after computing the linear
autocorrelation function, one could simply extract the raw values of the
function at a chosen set of time-lags, $\tau$ (e.g., $\tau=1,\dots,40$), but
could also capture different properties of the shape of the function (e.g., an
exponential fit), or extract a timescale estimate, such as the function’s
first $1/e$ or zero crossing. Similar choices need to be made when extracting
features from the output of many other time-series analysis methods, including
a Fourier or wavelet decomposition, the Lyapunov exponent spectrum, and many
others. In this process of coding the outputs of time-series analysis methods
as real-valued features, some feature sets that take a comprehensive approach,
in which many features are included and thus use multiple features to capture
similar properties, are expected to have higher feature–feature redundancy.
Other feature sets may take a more restrained approach, aiming to include a
smaller set of features that are less redundant with each other. When doing
statistical analysis with a feature set, it is important to understand the
effective number of independent features it measures, especially when there is
a high degree of collinearity between the included features.

We therefore aimed to investigate feature-set redundancy by assessing the
range of distinct behaviors exhibited by features in each of the seven feature
sets. Within-set redundancy was measured by the ability of low-dimensional
feature components to capture variance in the full feature set when applied to
the Empirical 1000 dataset [16]. A small number of low-dimensional components
will explain a large amount of variance in a feature set that contains lots of
redundant features (that produce highly correlated outputs across a diverse
time-series dataset). But feature sets containing linearly independent
features that display distinctive behavior on time-series data will require
many more components to explain a high percentage of variance across the full
set of features. Here we used the linear dimensionality-reduction method,
principal components analysis (PCA), on the 895 time series $\times$ filtered
feature matrix obtained from performing feature extraction for a given feature
set on the 895 time series in the Empirical 1000 dataset. Extracted principal
components were ordered by their percentage of explained variance.

To put feature sets of different sizes on a comparable scale, we expressed the
number of principal components (PCs) as a proportion of the total number of
PCs computed on the dataset. This ranged from 22 for catch22, up to a maximum
of 895 (the number of time series in the dataset). We then summarized
redundancy in a feature set using a simple summary statistic: the proportion
of PCs required to explain 90% of the variance in the feature space. Results
for this with full cumulative variance curves are plotted in Fig. 3A and the
overall summary statistic is plotted in Fig. 3B. The results reveal a
substantial variation in feature redundancy across feature sets. We observed
the highest level of within-set redundancy for TSFEL and tsfresh. For TSFEL,
90% of the variance across all 390 features can be captured with just 4 PCs
(1.0%), and for tsfresh, 90% of the variance of all 779 features was captured
with 55 PCs (7.4%). To investigate this result, we inspected the features that
make up these feature sets. We found that a large proportion of features in
both feature sets were derived from fast Fourier transform (FFT) coefficients:
256/390 (or 66%) of TSFEL features, and 400/779 (or 51%) of tsfresh features.
This large proportion of features derived from the FFT, in different forms of
the raw coefficients, contributes to the high redundancy of features in TSFEL
and tsfresh.

The most diverse feature set for its size was catch22, for which 11 PCs (50%)
were required to capture 90% of the variance across its 22 features. This is
consistent with the construction of catch22 as a reduction from 1000s of
initial features that included a clustering procedure to reduce redundancy
[6]. An intermediate level of redundancy is seen in feasts, tsfeatures, and
Kats, which have overlaps between their constituent features (the Kats feature
set was based on tsfeatures, and feasts includes five of the same features as
tsfeatures). Despite its goal of being as extensive as possible in surveying
time-series analysis algorithms, hctsa exhibited an intermediate level of
redundancy, with 24.9% of the PCs required to capture 90.0% of the variance
across all features.

Figure 3:  Feature sets vary substantially in their within-set redundancy. We
assess within-set redundancy in terms of the ability of low-dimensional
principal components to capture large amounts of variance in the feature space
comprising a diverse set of 895 real-world and model-generated time series
[16]. A Cumulative variance explained is plotted as a function of the
percentage of the number of principal components included for the seven
labeled feature sets. The dashed horizontal line indicates the 90% cumulative
variance threshold used in B. B Percentage of the number of principal
components required to explain 90% of the variance in the feature space.
Results are plotted for the seven feature sets, ordered from least to most
redundancy.

##  V Between-Set Overlap and Redundancy

Feature sets vary both in their coverage of time-series analysis methods and
in their style of distilling these methods in the form of features. Do some
sets contain many unique features (that are not contained in other sets), or
do some feature sets reproduce the behavior of features contained in other
feature sets? In this section, we aimed to address these questions by
investigating the similarity between features contained in different feature
sets. We used the same set of well-behaved features and 895 time series that
were used above (see Sec. III): catch22 (22 features), TSFEL (390), feasts
(41), Kats (37), tsfresh (745), tsfeatures (59), and hctsa (7708).

We developed a method to quantify the redundancy of a feature set, given
another feature set. To achieve this for a given pair of feature sets, we
first labeled one as a benchmark, $B$. We then aimed to capture how similar
features in the second set, labeled test, $T$, were to the features in the
benchmark set. If all individual features in the test feature set, $T$, behave
similarly to a feature that is contained in the benchmark set, $B$, then we
conclude that $T$ is redundant given $B$. To quantify the similarity between
two features, we used the absolute value of the Spearman correlation
coefficient, $|\rho|$, between their outputs on 895 diverse time series from
the Empirical 1000 dataset [16]. In this way, features with highly correlated
outputs on these diverse data are judged as having similar behavior. We
defined an ordered similarity metric, $\mathcal{S}$, to judge the redundancy
of a feature set, $T$, given a benchmark set, $B$, based on the similarities
of their constituent features as

$\mathcal{S}(T|B)=\langle|\rho|^{\mathrm{max}}_{i}\rangle\,,$ (1)

where $|\rho|^{\mathrm{max}}_{i}=\max_{j}|\rho|_{ij}$ is the maximum Spearman
rank correlation between each feature, $i$, in $T$, and all features, $j$, in
$B$. The average, $\langle\cdot\rangle$, is then taken over all features, $i$,
in $T$, yielding $\mathcal{S}(T|B)$, the mean absolute maximum correlation
between features in $T$, given those that exist in the benchmark, $B$.

Importantly, $\mathcal{S}$ is asymmetric, allowing us to understand the
relative overlaps and redundancies between feature sets. For example, consider
a comprehensive feature set, $X$, containing many diverse and useful time-
series features, and another simple set, $Y$, that only contains
autocorrelation-related features. Then $\mathcal{S}(Y|X)$ will be high
(indicating that there is, on average, a feature in $X$ that is similar to
each feature in $Y$), but $\mathcal{S}(X|Y)$ will be low (indicating that the
simple features in $Y$ cannot capture the wider range of behavior of the
features in $X$). In the case that $X$ and $Y$ are sets containing very
similar features, then both $\mathcal{S}(Y|X)$ and $\mathcal{S}(X|Y)$ will be
high. We computed $\mathcal{S}(T|B)$ for all pairs of feature sets, treating
each feature set as both a benchmark and a test. The result was represented as
a pairwise matrix, with each benchmark feature set, $B$, shown as rows, and
each test feature set, $T$, shown as columns, as shown in Fig. 4.

The matrix reveals a range of interesting overlaps and redundancies between
feature sets. Overall, there are many high values, $\mathcal{S}$, indicating
substantial overlap between the types of features included in the feature
sets. Looking across a row, feature sets with many high values, $\mathcal{S}$,
are ‘comprehensive’ in the sense that they contain features that exhibit
similar behavior as those in the other feature sets. As expected, hctsa, the
package with the most numerous and comprehensive coverage of time-series
analysis algorithms, was the strongest benchmark set: on average, there exists
an feature in hctsa that is highly correlated to features in catch22
($\langle|\rho|^{\mathrm{max}}_{i}\rangle=0.96$), feasts
($\langle|\rho|^{\mathrm{max}}_{i}\rangle=0.84$), Kats
($\langle|\rho|^{\mathrm{max}}_{i}\rangle=0.84$), tsfeatures
($\langle|\rho|^{\mathrm{max}}_{i}\rangle=0.88$), and TSFEL
($\langle|\rho|^{\mathrm{max}}_{i}\rangle=0.88$). The exception is tsfresh,
for which $\langle|\rho|^{\mathrm{max}}_{i}\rangle=0.55$, indicating that, on
average, for a feature selected from hctsa, the most correlated tsfresh
feature will have $|\rho|=0.55$. Indeed, tsfresh is characterized by its
distinctiveness overall: when considered as a test set, $T$ (inspecting the
tsfresh column in Fig. 4), $\mathcal{S}(\texttt{tsfresh}|B)$ values are even
lower ($0.25-0.5$) for all other feature sets treated as benchmarks, $B$. This
indicates that tsfresh contains many unique time-series features that are not
present in the other feature sets.

To investigate this uniqueness, we examined the individual features in tsfresh
that were the least correlated to features in hctsa; in particular, those with
$|\rho|^{\mathrm{max}}<0.2$ (the maximum taken across all hctsa features). We
found 128 tsfresh features with $|\rho|^{\mathrm{max}}<0.2$ to features in
hctsa. Of these, 118 (or 92%) were Fourier coefficients (real and imaginary
components, and angles). While hctsa contains summary statistics derived from
the Fourier power spectrum, they do not output the set of all raw coefficients
from the fast Fourier transform (FFT). This appears to drive the unique
behavior of tsfresh, and also contributes to its relatively high within-set
redundancy, cf. Fig. 3. Although TSFEL computes the absolute value of a large
number of FFT coefficients, it does not also return the real and imaginary
components and angles as features, and thus has a low
$\mathcal{S}(\texttt{TSFEL}|\texttt{tsfresh})=0.5$.

As expected, given the history of their construction, tsfeatures and feasts
were strongly overlapping in their behavior (in both directions,
$\mathcal{S}\approx 0.82$), with high similarity scores also with the
similarly constructed Kats ($\mathcal{S}\approx 0.72$).

Figure 4:  Time-series feature sets show a wide variation of overlap and
redundancy. We assessed the redundancy of a given feature set by treating it
as a ‘test’ feature set, $T$, given another feature set, which is treated as a
‘benchmark’, $B$. The resulting directed similarity metric,
$\mathcal{S}(T|B)$, Eq. (1) (a mean of the maximum absolute Spearman
correlation coefficients between feature outputs across a diverse time-series
dataset), gives high values of $\mathcal{S}$ when features in $T$ reproduce
the behavior of features that exist in $B$, and low values when features in
$T$ are on average uncorrelated (distinctive) from the features in $B$. Here
we plot $\mathcal{S}(T|B)$ as a heatmap, where all feature sets have been
treated as both a benchmark feature set (rows) and a test feature set
(columns).

##  VI Discussion

The feature-based approach to time-series analysis, in which a range of
algorithms underpinned by theory can be used to distill a large time-series
dataset down to a set of interpretable summary statistics is a promising way
to address important time-series analysis problems in science and industry.
Here we characterized the behavior of, and overlaps between, the most
prominent existing software libraries for time-series feature extraction:
hctsa, catch22, tsfeatures, feasts, Kats, tsfresh, and TSFEL. We took an
empirical approach of characterizing features in terms of their pattern of
outputs across a set of over 800 diverse real-world and model-generated time
series, allowing us to distinguish ‘similar’ features as those with highly
correlated behavior on data. We found that TSFEL and catch22 were the fastest
to evaluate for time series in the range 100–1000 samples ($\sim 10^{-4}$ s
per feature). However, our low-dimensional analysis highlighted a large amount
of redundancy in particularly TSFEL and tsfresh, with a large proportion of
highly correlated features (including many raw Fourier coefficients) in both
feature sets. By contrast, catch22, a feature set designed to reduce within-
set redundancy, displayed the least within-set redundancy. We also analyzed
the overlap between feature sets by introducing a directed pairwise metric of
feature-set redundancy, that treats one set as a ‘benchmark’ and measures how
correlated features in a second ‘test’ feature set are to those in the
benchmark set. As expected, hctsa was distinguished as the most comprehensive
feature set, with features in catch22, feasts, Kats, tsfeatures, and TSFEL all
having a highly correlated match ($\langle\rho\rangle\geq 0.84$) to a feature
contained in hctsa. The exception was tsfresh, which was marked as the most
unique feature set, largely driven by a large number of features that capture
the raw outputs of a fast Fourier transform as real and imaginary components,
magnitudes, and phase angles (whereas many other packages extract summary
statistics from the Fourier power spectrum).

Our results also highlight different choices made in feature composition and
behavior across the seven time-series feature extraction libraries. It is
important to understand the characteristics of different feature sets, as
different applications of time-series feature extraction call for different
choices, such as prioritizing speed or accuracy. Different feature sets take
different approaches to feature-set construction, from the large hctsa feature
set, which was developed with the aim of representing a comprehensive survey
of the time-series analysis literature, to its distillation as the fast and
less correlated features in catch22. We also found a large amount of
redundancy in many of the feature sets, such as TSFEL and tsfresh, where the
variance across a large number of included features can be well captured using
a small proportion ($<10$%) of reduced components. This redundancy was driven
by large numbers of features quantifying similar types of time-series
structure (such as Fourier coefficients or quantiles of the distribution of
values). Feature comprehensiveness comes at the cost of redundancy, and could
be beneficial for some applications, e.g., involving large datasets and for
which raw Fourier coefficients (including phase information of specific
components) is relevant. When applying statistical models to data, it is also
important to be aware of and properly account for high collinearity of
features within a feature set, which is substantial in all feature sets
analyzed here.

While this work focused on the similarity of feature behavior on diverse time-
series data (from the Empirical 1000 dataset [16, 17]), it will be important
for future work to properly evaluate the relative performance of the different
feature sets on different types of problems involving time series. For
example, the most common application of feature-based time-series analysis has
been to time-series classification tasks—comparing the performance of
different feature sets for different types of classification (and other) tasks
will allow for clearer understanding of the relative strengths and weaknesses
of these feature sets. Other tasks which may yield informative differences
between feature sets include large-scale anomaly detection [21] and feature-
based forecasting [22]. Our results may also depend on our choice to quantify
feature similarity using a wide variety of data from the Empirical 1000
dataset [16]. While features that are highly correlated across such diverse
data can be considered to be measuring similar time-series properties in
similar ways, feature–feature relationships may differ on more specific
datasets, or on time series of lengths outside of the length range ($\leq
1000$) considered here. We note that the R package, theft [15], provides easy
access to the six Python and R feature sets analyzed here (excluding the
Matlab-based hctsa), allowing researchers to leverage combinations of features
from different feature sets. In the future, we hope that comparisons of time-
series analysis methods across disciplines and feature sets will highlight
fruitful directions for leveraging the combined strengths of diverse features,
towards optimising an analysis approach for a given task.

##  VII Available Code

Code to reproduce all analysis included in this paper (including downloading
and processing the Empirical 1000 dataset [16]) is available at
https://github.com/hendersontrent/feature-set-comp.

## References

  * [1] B. D. Fulcher, A. E. Georgieva, C. W. G. Redman, and N. S. Jones, “Highly comparative fetal heart rate analysis,” 34th Annual International Conference of the IEEE EMBS, 2012.
  * [2] B. D. Fulcher and N. S. Jones, “Highly comparative feature-Based time-series classification,” IEEE Transactions on Knowledge and Data Engineering, vol. 26, no. 12, pp. 3026–3037, 1 Dec. 2014.
  * [3] B. D. Fulcher, “Feature-based time-series analysis,” Feature engineering for machine learning and data analytics, pp. 87–-116, CRC Press, 2018.
  * [4] B. D. Fulcher, M. A. Little, and N. S. Jones, “Highly comparative time-series analysis: the empirical structure of time series and their methods,” J. Roy. Soc. Interface 10, 20130048, 2013.
  * [5] B. D. Fulcher and N. S. Jones. “hctsa: A computational framework for automated time-series phenotyping using massive feature extraction,” Cell Systems, 5, 527, 2017.
  * [6] C. H. Lubba, S. S. Sethi, P. Knaute, S. R. Schultz, B. D. Fulcher, and N. S. Jones, “catch22: CAnonical Time-series CHaracteristics.,” Data Min Knowl Disc, 33, pp. 1821-–1852, 2019.
  * [7] T. Henderson, “Rcatch22: Calculation of 22 CAnonical Time-Series CHaracteristics,” R package version 0.1.12, Zenodo, https://doi.org/10.5281/zenodo.4851882, 2021.
  * [8] B. Harris. “brendanjohnharris/Catch22.jl: v0.2.1,” Zenodo, https://doi.org/10.5281/zenodo.5030712, 2021.
  * [9] M. O’Hara-Wild, R. J. Hyndman, and E. Wang, “feasts: Feature Extraction and Statistics for Time Series,” R package version 0.2.1, 2021.
  * [10] R. J. Hyndman, Y. Kang, P. Montero-Manso, T. Talagala, E. Wang, Y. Yang, and M. O’Hara-Wild, “tsfeatures: Time series feature extraction,” R package version 1.0.2, 2020.
  * [11] M. Christ, N. Braun, J. Neuffer, and A. W. Kempa-Liehr, “Time Series FeatuRe Extraction on basis of Scalable Hypothesis tests (tsfresh – A Python package),” Neurocomputing, 307, pp.72-77, 2018.
  * [12] M. Barandas, D. Folgado, Fernandes, L, et al., “TSFEL: Time Series Feature Extraction Library,” SoftwareX 11, 2020.
  * [13] Facebook Infrastructure Data Science. “Kats,” https://facebookresearch.github.io/Kats/
  * [14] R. Hegger, H. Kantz, and T. Schreiber, “Practical implementation of nonlinear time series methods: The TISEAN package,” CHAOS, 9, 413, 1999.
  * [15] T. Henderson, “theft: Tools for Handling Extraction of Features from Time series,” R package version 0.1.32, Zenodo, https://doi.org/10.5281/zenodo.5220019, 2021.
  * [16] B. D. Fulcher, “1000 Empirical Time series,” figshare, Dataset, https://doi.org/10.6084/m9.figshare.5436136.v9, 2017.
  * [17] B. D. Fulcher, C. H. Lubba, S. S. Sethi, and N. S. Jones, “A self-organizing, living library of time-series data,” Scientific Data, 7, 213, 2020.
  * [18] J. C. Sprott, “Chaos and Time-Series Analysis,” Oxford University Press, Inc., USA, 2003.
  * [19] A. Goldberger, L. Amaral, L. Glass, J. Hausdorff, P. C. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C. K. Peng, and H. E. Stanley, “PhysioBank, PhysioToolkit, and PhysioNet: Components of a new research resource for complex physiologic signals,” Circulation [Online], 101 (23), pp. e215-–e220, 2000.
  * [20] J. Hausdorff, “Gait in Parkinson’s Disease v1.0.0,” https://doi.org/10.13026/C24H3N, 2008.
  * [21] R. J. Hyndman, E. Wang, and N. Laptev, “Large-Scale Unusual Time Series Detection,” 2015 IEEE International Conference on Data Mining Workshop (ICDMW), pp. 1616-1619, 2015.
  * [22] P. Montero-Manso, G. Athanasopolous, R. J. Hyndman, and T. S. Talagala, “FFORMA: Feature-based forecast model averaging,” International Journal of Forecasting, 36, pp. 86-92, 2020.


# DDPET-3D: Dose-aware Diffusion Model for 3D Ultra Low-dose PET Imaging

Huidong Xie∗
Yale University
<EMAIL_ADDRESS>Weijie Gan∗
Washington University in St. Louis
<EMAIL_ADDRESS>Bo Zhou
Yale University
<EMAIL_ADDRESS>Xiongchao Chen
Yale University
<EMAIL_ADDRESS>Qiong Liu
Yale University
<EMAIL_ADDRESS>Xueqi Guo
Yale University
<EMAIL_ADDRESS>Liang Guo
Yale University
<EMAIL_ADDRESS>Hongyu An
Washington University in St. Louis
<EMAIL_ADDRESS>Ulugbek S. Kamilov
Washington University in St. Louis
<EMAIL_ADDRESS>Ge Wang
Rensselaer Polytechnic Institute
<EMAIL_ADDRESS>Chi Liu
Yale University
<EMAIL_ADDRESS>
###### Abstract

As PET imaging is accompanied by substantial radiation exposure and cancer
risk, reducing radiation dose in PET scans is an important topic. Recently,
diffusion models have emerged as the new state-of-the-art generative model to
generate high-quality samples and have demonstrated strong potential for
various tasks in medical imaging. However, it is difficult to extend diffusion
models for 3D image reconstructions due to the memory burden. Directly
stacking 2D slices together to create 3D image volumes would results in severe
inconsistencies between slices. Previous works tried to either apply a penalty
term along the z-axis to remove inconsistencies or reconstruct the 3D image
volumes with 2 pre-trained perpendicular 2D diffusion models. Nonetheless,
these previous methods failed to produce satisfactory results in challenging
cases for PET image denoising. In addition to administered dose, the noise
levels in PET images are affected by several other factors in clinical
settings, e.g. scan time, medical history, patient size, and weight, etc.
Therefore, a method to simultaneously denoise PET images with different noise-
levels is needed. Here, we proposed a Dose-aware Diffusion model for 3D low-
dose PET imaging (DDPET-3D) to address these challenges. We extensively
evaluated DDPET-3D on 100 patients with 6 different low-dose levels (a total
of 600 testing studies), and demonstrated superior performance over previous
diffusion models for 3D imaging problems as well as previous noise-aware
medical image denoising models. The code is available at:
https://github.com/xxx/xxx.

††footnotetext: ∗ Equal Contributions

##  1 Introduction

Positron Emission Tomography (PET) is a functional imaging modality widely
used in oncology, cardiology, and neurology studies [19, 23, 4]. Given the
growing concern about radiation exposure and cancer risks accompanied with PET
scans, reducing the PET injection dose is desirable [18]. However, PET image
quality is negatively affected by the reduced injection dose and may affect
diagnostic performance such as the identification of low-contrast lesions
[22]. Therefore, reconstructing high-quality images from noisy input is an
important topic.

Iterative methods such as Maximum Likelihood Expectation Maximization (MLEM)
[24] and Ordered Subset Expectation Maximization (OSEM) [11] are commonly used
for PET reconstructions. However, they are vulnerable to noise in low-dose PET
data. Iterative methods incorporating a Total Variation (TV) regularization
could be used for noise reduction [27]. However, they are time-consuming and
often fail to produce high-quality reconstructions.

Deep learning has emerged as a new reconstruction algorithm for medical
imaging tasks [28], and many different deep-learning methods were proposed for
low-dose PET image reconstructions [30, 35, 17, 31, 33, 32, 34, 15]. Recently,
diffusion models have become the new state-of-the-art generative models [5,
13]. They are capable of generating high-quality samples from Gaussian noise
input, and have demonstrated strong potential for low-dose PET imaging. For
example, using the DDPM framework [10], Gong et al. proposed to perform PET
image denoising with MRI as prior information for improved image quality [8].
However, PET-MR systems are not widely available and it is difficult to obtain
paired PET-MR images in reality for network training and image
reconstructions. Jiang et al.[12] adopted the latent diffusion model [20] for
unsupervised PET denoising. Moreover, diffusion models have also been proposed
for other imaging modalities, such as CT [7] and MRI [9, 2].

However, these previous diffusion papers focus on 2D and do not address the 3D
imaging problem. This is particularly important for PET imaging as PET is
intrinsically a 3D imaging modality. A 3D diffusion model is desirable,
however, due to the hardware memory limit, directly extending the diffusion
model to 3D would be difficult. There are a few previous works that aim to
address the 3D imaging problems of diffusion models. For example, Chung et
al.[3] proposed to apply a TV penalty term along the z-axis to remove
inconsistencies within each reverse sampling step in the diffusion model. On
the other hand, Lee et al.[14] utilizes 2 pre-trained perpendicular 2D
diffusion models to remove inconsistencies between slices. Nonetheless, these
previous methods failed to produce satisfactory results in challenging cases
for 3D PET image denoising.

Another challenge with PET image denoising is the high variation of image
noise. The final noise level in PET images can be affected by different
factors: _(1)_ Variations of acquisition start time. _(2)_ Occasional tracer
injection infiltration in some patients, causing the tracer stuck in the arm,
resulting in high image noise in the body. _(3)_ The variation of injection
dose based on the patient’s medical history, weight, and size. _(4)_ Different
hospitals have different standard scan times for patients, resulting in
variability of noise levels in the reconstructed images. Because of these
reasons, a method to simultaneously denoise images with varying noise levels
is desirable. However, previously proposed methods mentioned above have
limited generalizability to different noise levels [29]. A network trained on
one noise level fails to produce high-quality reconstructions on other noise
levels. To address this problem, Xie et al.  proposed to combine multiple
U-net-based [21] sub-networks with varying denoising power to generate optimal
results for any input noise levels [29]. However, training multiple sub-
networks requires tedious data pre-processing and long training time. The
testing time also linearly increases with the number of sub-networks. Also,
paired training data with different low-count levels may not be readily
available.

Moreover, different from other imaging modalities, PET was developed as a
quantitative tool for clinical diagnosis. PET quantitative characteristics are
increasingly being recognized as providing an objective, and more accurate
measure for prognosis and response monitoring purposes than visual inspection
alone [1]. However, our experimental results showed that, although standard
diffusion models (DDPM [10], DDIM[25]) produce visually appealing
reconstructions, they typically failed to maintain accurate quantification.
For example, the total activities change after diffusion models. This is
particularly true for whole-body PET scans since the tracer uptake could vary
significantly in different organs. This problem was not able to address by
simple normalization of the images.

In this work, we developed a dose-aware diffusion model for 3D PET Imaging
(DDPET-3D) to address these limitations. The main contributions of the
proposed DDPET-3D framework are: _(1)_ We proposed a dose-embedding strategy
that allows noise-aware denoising. DDPET-3D can simultaneously denoise 3D PET
images with varying low-dose/count levels. _(2)_ We proposed a 2.5D diffusion
strategy with multiple fixed noise variables to address the 3D inconsistency
issue between slices. DDPET-3D maintains a similar memory burden to 2D
diffusion models while achieving high-quality reconstructions. _(3)_ We
proposed to use a denoised prior in DDPET-3D, allowing it to converge within
25 sampling steps. DDPET-3D can denoise 3D PET images within a reasonable time
constraint ($15$ mins on a single GPU). Previous diffusion methods using DDPM
sampling would require approximately $6$ hrs. Evaluated on real-world ultra-
low-dose 3D PET data, DDPET-3D demonstrates superior quantitative and
qualitative results as compared with previous baseline methods. External
validation of data from another country also shows our trained model can be
reasonably generalized to other sites.

##  2 Methods

###  2.1 Diffusion Models

The general idea of diffusion models is to learn the target data distribution
$q({\bm{x}}_{0})$ (_i.e.,_ full-dose PET images in our case) using neural
network. Once the distribution is learned, we can synthesize a new sample from
it. Diffusion models consist of two Markov chains: the forward diffusion
process and the learned reverse diffusion process. The forward diffusion
process $q$ gradually adds small amount of Gaussian noise to ${\bm{x}}_{0}\sim
q({\bm{x}}_{0})$ in each step, until the original image signal is completely
destroyed. As defined in the DDPM paper [10]

$\begin{split}q({\bm{x}}_{1:T}|{\bm{x}}_{0})&\coloneqq\prod_{t=1}^{T}q({\bm{x}}_{t}|{\bm{x}}_{t-1})\
,\\\
\text{where}\;q({\bm{x}}_{t}|{\bm{x}}_{t-1})&\coloneqq\mathcal{N}({\bm{x}}_{t};\sqrt{1-\beta_{t}}{\bm{x}}_{t-1},\beta_{t}\mathbf{I})\
.\end{split}$ (1)

One property of the diffusion process is that, one can sample ${\bm{x}}_{t}$
for any arbitrary time-step $t$ without gradually adding noise to
${\bm{x}}_{0}$. By denoting $\alpha_{t}\coloneqq 1-\beta_{t}$ and
$\bar{\alpha}_{t}\coloneqq\prod_{s=1}^{t}\alpha_{s}$, we have

$\begin{split}q({\bm{x}}_{t}|{\bm{x}}_{0})&=\mathcal{N}({\bm{x}}_{t};\sqrt{(\bar{\alpha}_{t}}){\bm{x}}_{0},(1-\bar{\alpha}_{t})\mathbf{I})\
,\\\ \text{and}\
{\bm{x}}_{t}&=\sqrt{\bar{\alpha}_{t}}{\bm{x}}_{0}+\sqrt{1-\bar{\alpha}_{t}}{\bm{\epsilon}}\
,\end{split}$ (2)

where ${\bm{\epsilon}}\sim\mathcal{N}(0,\mathbf{I})$. The latent
${\bm{x}}_{T}$ is nearly an isotropic Gaussian distribution for a properly
designed $\beta_{t}$ schedule. Therefore, one can easily generate a new
${\bm{x}}_{T}$ and then synthesize a ${\bm{x}}_{0}$ by progressively sampling
from the reverse posterior $q({\bm{x}}_{t-1}|{\bm{x}}_{t})$. However, this
reverse posterior is tractable only if ${\bm{x}}_{0}$ is known

$q({\bm{x}}_{t-1}|{\bm{x}}_{t},{\bm{x}}_{0})={\mathcal{N}}\Big{(}{\bm{x}}_{t-1};\mu_{q}({\bm{x}}_{t},{\bm{x}}_{0}),\frac{\beta_{t}(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_{t}}\textbf{I}\Big{)}\
,$ (3)

where

$\mu_{q}({\bm{x}}_{t},{\bm{x}}_{0})=\frac{\sqrt{\alpha_{t}}(1-\bar{\alpha}_{t-1}){\bm{x}}_{t}+\sqrt{\bar{\alpha}_{t-1}}(1-\alpha_{t}){\bm{x}}_{0}}{1-\bar{\alpha}_{t}}\
.$ (4)

Note that
$q({\bm{x}}_{t-1}|{\bm{x}}_{t}):=q({\bm{x}}_{t-1}|{\bm{x}}_{t},{\bm{x}}_{0})$,
where the extra conditioning term ${\bm{x}}_{0}$ is superfluous due to the
Markov property. DDPM thus proposes to learn a parameterized Gaussian
transitions $p_{\bm{\theta}}({\bm{x}}_{t-1}|{\bm{x}}_{t})$ to approximate the
reverse diffusion posterior (3)

$p_{\bm{\theta}}({\bm{x}}_{t-1}|{\bm{x}}_{t})={\mathcal{N}}\Big{(}{\bm{x}}_{t-1};\mu_{\bm{\theta}}({\bm{x}}_{t},t),\sigma^{2}_{t}\textbf{I}\Big{)}\
,$ (5)

where

$\mu_{\bm{\theta}}({\bm{x}}_{t},t)=\frac{1}{\sqrt{\alpha_{t}}}\Big{(}{\bm{x}}_{t}-\frac{1-\alpha_{t}}{\sqrt{1-\bar{\alpha}_{t}}}\epsilon_{\bm{\theta}}({\bm{x}}_{t},t)\Big{)}\
.$ (6)

Here, $\epsilon_{\bm{\theta}}$ denotes a neural network. Through some
derivations detailed in the DDPM paper [10], the training objective of
$\epsilon_{\bm{\theta}}({\bm{x}}_{t},t)$ can be formulated as follow

$\mathbb{E}_{{\bm{x}},{\bm{\epsilon}},t\sim[1,T]}\big{[}\left\lVert{\bm{\epsilon}}-{\epsilon}_{\bm{\theta}}({\bm{x}}_{t},t)\right\rVert^{2}\big{]}\
.$ (7)

It is worth noting that the original DDPM set $\sigma_{t}$ to a fixed constant
value based on the $\beta_{t}$ schedule. Recent studies [6, 16] have shown the
improved performance by using the learned variance
$\sigma_{t}^{2}\coloneqq\sigma_{\bm{\theta}}^{2}({\bm{x}}_{t},t)$. We also
adopted this approach. To be specific, we have
$\sigma_{\bm{\theta}}({\bm{x}}_{t},t)\coloneqq\exp(v\log\beta_{t}+(1-v)\log\tilde{\beta}_{t})$,
where $\tilde{\beta}_{t}$ refers to the lower bounds for the reverse diffusion
posterior variances [10], and $v$ denotes the network output. We used a single
neural network with two seperate output channels to estimate the mean and the
variance of (6) jointly. Based on the learned reverse posterior
$p_{\bm{\theta}}({\bm{x}}_{t-1}|{\bm{x}}_{t})$, the iteration of obtaining a
${\bm{x}}_{0}$ from a ${\bm{x}}_{T}$ can be formulated as follow

${\bm{x}}_{t-1}=\mu_{\bm{\theta}}({\bm{x}}_{t},t)+\sigma_{t}{\bm{z}},\text{
where }{\bm{z}}\sim{\mathcal{N}}(0,\textbf{I})\ .$ (8)

Figure 1: The training (_top_) and the sampling (_bottom_) pipeline of the
proposed DDPET-3D method. The DDPET-3D uses multiple neighboring slices as
additional inputs to predict the the central slice, allowing network to
observe 3D information during training and testing. The DDPET-3D is also
conditioned on the injected dose in order to accommodate inputs with varying
noise levels. In sampling, we propose to fix the Gaussian latent for each
slice in the 3D volume to address the inconsistency along z-axis. DDPET-3D
also use a pre-trained denoiser prior to ensure accurate quantification in
sampling.

###  2.2 Conditional PET Image Denoising

The framework described above only allows unconditional sampling. For the
purpose of low-dose PET denoising, instead of generating new samples, the
network needs to denoise the images based on input noisy counterparts. This
can be achieved by adding additional condition to the neural network.
$\epsilon_{\bm{\theta}}{({\bm{x}}_{t},t)}$ becomes
$\epsilon_{\bm{\theta}}{({\bm{x}}_{t},t,{\bm{x}}_{\mathrm{noisy}})}$, where
${\bm{x}}_{\mathrm{noisy}}$ denotes the noisy input PET images. Specifically,
the input to the neural network becomes a 2-channel input, one is
${\bm{x}}_{t}$, the other one is ${\bm{x}}_{\mathrm{noisy}}$. However, we
noticed that, this technique results in severe inconsistencies in the
reconstructed 3D image volumes. One example is presented in Fig. 2
(DDIM50+2.5D results).

###  2.3 Proposed DDPET-3D

The proposed framework is depicted in Fig. 1. The proposed network can observe
3D information from adjacent slices during the training process. This is
achieved by using the $n$ neighboring slices as input to predict the central
slice. Ablation studies showed that the performance converges at $n=31$, which
was used in this work. One may use 3D convolutional layers to replace the 2D
convolutional layers in the diffusion model for an enlarged receptive field
and to allow the network to observe 3D structural information. However, we
found that it significantly increases memory burden and makes the network
difficult to optimize. To alleviate the memory burden and allow faster
convergence, we embed neighboring slices in the channel dimension.
Specifically, when using only one 2D slice as conditional information ($n=1$),
the input dimension is $N_{b}\times W\times W\times 2$, where the last
dimension is the channel dimension, $N_{b}$ is input batch size, and $W$ is
the width of the images. One channel is conditional 2D slice, and the other
one is ${\bm{x}}_{t-1}$. When $n=31$, the input dimension becomes $N_{b}\times
W\times W\times 32$.

Such technique allows the network to observe neighboring slices for 3D image
reconstruction with only incremental increase in memory burden. Also, ablation
studies showed that such technique also noticeably improves inconsistencies
problems for 3D imaging. This is consistent with our expectation. With $n=1$,
when the network tries to predict the next slice, we observe inconsistent
reconstructions because both the conditional information and the starting
Gaussian noise $p({\bm{x}}_{T})$ change. When $n=31$, the starting Gaussian
noise $p({\bm{x}}_{T})$ still changes, but all the 31-channel conditional
information only shifts a little. The network should produce more consistent
output with only subtle changes in the input.

Another reason for the inconsistency issue is because the starting Gaussian
noise of the reverse process is different for different slices. The diffusion
sampling strategy starts from a random location in the high-dimensional space
(random Gaussian noise), and approximates the data distribution
$q({\bm{x}}_{0})$ based on the trained neural network, or the score function
as described in [26], through many iterations. The denoising problem is ill-
posed, and we could generate an infinite number of different denoised images
given the same low-count input with different starting Gaussian noise. This is
beneficial for generative models to produce a wide variety of different
images. However, the stochastic nature of the diffusion model would be
problematic for medical image reconstruction problems because we expect
neighboring slices to be consistent with each other in the 3D image volume. To
address this problem, we proposed to fix the starting Gaussian noise when
reconstructing all the slices in the entire 3D volume. Specifically, the
starting Gaussian noise at the last time step $T$, ${\bm{x}}_{T}$ is fixed for
all the slices during the sampling process, only the conditional low-count PET
images change for different slices.

This approach produced more consistent reconstructions along the z-axis.
However, we noticed that by fixing $p({\bm{x}}_{T})$, noise-dependent
artifacts would propagate to all the slices along the z-axis. One example is
presented in the first sub-figure in Fig. 6 (single $\epsilon$ image). To
address this issue, as presented in Fig. 1, we proposed to initialize 2
different noise variables ${\bm{\epsilon}}_{0}^{a}$ and
${\bm{\epsilon}}_{0}^{b}$, and we have 2 different ${\bm{x}}_{t-1}$ in the
reverse process. Fortunately, having 2 different noise variables does not
significantly increase sampling time since we only need to average them at the
first reverse step. This technique effectively addresses this issue.

As mentioned in the introduction section, despite the visually appealing
results, the diffusion model usually produced images with inaccurate
quantification. We noticed that, even though U-net-based methods failed to
produce satisfactory results, especially for extremely low-dose settings, they
typically maintained overall image quantification much better than diffusion
models. To take advantage of that, instead of starting from random Gaussian
noise, we first generate a denoised prior using a pre-trained U-net-based
network, and then add Gaussian noise to it based on Equation 2. This added-
noise denoised prior were used as the starting point of the reverse process.
We adapted the Unified Noise-aware Network (UNN) proposed by Xie et al. [29]
to generate the denoised priors.

For injected dose embedding, the values of the administered dose were first
converted to Becquerel (Bq). Both $\sin$ and $\cos$ functions were then used
to encode the injected dose and add with the diffusion time-steps (i.e.,
$t+\sin(\mathrm{dose})+\cos(\mathrm{dose})$, where $t$ is the diffusion time-
step). Time-step $t$ was also encoded using $\sin$ and $\cos$ functions. The
encoded values were then fed into 2 linear layers to generate the embedding
with Sigmoid linear unit (SiLU) in between.

For diffusion sampling, DDIM sampling enables faster convergence. But we
noticed that using DDIM alone would tend to produce over-smoothed images. To
alleviate this issue, adapted from the DDIM paper [25], we proposed to have an
interpolated sampling between DDIM and DDPM. We inserted DDPM samplings every
5 steps in between DDIM samplings to prevent the images from becoming over-
smoothed. Algorithm 1 displays the complete training procedure of the proposed
method. Algorithm 2 displays the complete sampling procedure of the proposed
method.

Algorithm 1 Training

Repeat

${\bm{x}}_{0}\sim q({\bm{x}}_{0})$ $\triangleright$ Sample single-slice full-
dose data

${\bm{x}}_{\mathrm{noisy}}\sim q({\bm{x}}_{\mathrm{noisy}})$ $\triangleright$
Sample multi-slice low-dose data

$\mathrm{dose}_{{\bm{x}}_{\mathrm{noisy}}}$ $\triangleright$ Injected dose of
${\bm{x}}_{\mathrm{noisy}}$

$t\sim\mathrm{Uniform({1,...,T})}$ $\triangleright$ Sample diffusion time-step

${\bm{\epsilon}}\sim\mathcal{N}(0,\mathbf{I})$

$\nabla_{\theta}||{\bm{\epsilon}}-\epsilon_{\bm{\theta}}(\sqrt{\bar{\alpha}_{t}}{\bm{x}}_{0}+\sqrt{1-\bar{\alpha}_{t}}{\bm{\epsilon}},t,{\bm{x}}_{\mathrm{noisy}},\mathrm{dose}_{{\bm{x}}_{\mathrm{noise}}})$

$\triangleright$ Take gradient descent step

Until convergence

Algorithm 2 Testing

${\bm{x}}_{\mathrm{noisy}}\sim q({\bm{x}}_{\mathrm{noisy}})$ $\triangleright$
Get entire low-dose image volume

$\mathrm{dose}_{{\bm{x}}_{\mathrm{noisy}}}$ $\triangleright$ Injected dose of
${\bm{x}}_{\mathrm{noisy}}$

${\bm{x}}_{\mathrm{prior}}=\mathrm{UNN}({\bm{x}}_{\mathrm{noisy}})$
$\triangleright$ Generate denoised prior using

a pre-trained network (UNN) [29]

${\bm{x}}_{\mathrm{recon}}\leftarrow\mathrm{zeros}({\bm{x}}_{\mathrm{noisy}})$
$\triangleright$ Initialize ${\bm{x}}_{\mathrm{recon}}$

${\bm{\epsilon}}_{0}^{a}\sim\mathcal{N}(0,\mathbf{I})$,
${\bm{\epsilon}}_{0}^{b}\sim\mathcal{N}(0,\mathbf{I})$ $\triangleright$ Obtain
2 noise variables

while $s=1,...,S$ do $\triangleright$ Total number of slices

${\bm{x}}_{\mathrm{noisy}}[s]$ $\triangleright$ Conditioned multi-slice low-
dose input

${\bm{x}}_{\mathrm{prior}}[s]$ $\triangleright$ Conditioned single-slice
denoised prior

${\bm{x}}_{T^{\prime}}^{a}=\sqrt{\bar{\alpha}_{T^{\prime}}}{\bm{x}}_{\mathrm{prior}}[s]+\sqrt{1-\bar{\alpha}_{T^{\prime}}}{\bm{\epsilon}}_{0}^{a}$

${\bm{x}}_{T^{\prime}}^{b}=\sqrt{\bar{\alpha}_{T^{\prime}}}{\bm{x}}_{\mathrm{prior}}[s]+\sqrt{1-\bar{\alpha}_{T^{\prime}}}{\bm{\epsilon}}_{0}^{b}$

while $t=T^{\prime},...,1$ do

$z\sim\mathcal{N}(0,\mathbf{I})$

${\bm{x}}_{t-1}^{a}=\mathrm{Sampler}({\bm{x}}_{t}^{a},{\bm{x}}_{\mathrm{noisy}}[s],\mathrm{dose}_{{\bm{x}}_{\mathrm{noisy}}})$

if $t==T^{\prime}$ then

${\bm{x}}_{t-1}^{b}=\mathrm{Sampler}({\bm{x}}_{t}^{b},{\bm{x}}_{\mathrm{noisy}}[s],\mathrm{dose}_{{\bm{x}}_{\mathrm{noisy}}})$

${\bm{x}}_{t-1}^{a}=({\bm{x}}_{t-1}^{a}+{\bm{x}}_{t-1}^{b})/2$

end if

end while

${\bm{x}}_{\mathrm{recon}}[s]={\bm{x}}_{1}^{a}$

end while

Return ${\bm{x}}_{\mathrm{recon}}$

###  2.4 Low-count PET Data

The PET dataset used in this study was collected at XX. 320 subjects with
18F-FDG tracer were acquired using a United Imaging uExplorer PET/CT system.
The reconstruction matrix size is $673\times 360\times 360$ with a $2.89\times
1.67\times 1.67$\mathrm{m}\mathrm{m}$^{3}$ voxel size. We randomly selected
210 subjects for training, 10 subjects for validation, and 100 subjects for
testing. We downsampled the PET list-mode data to 1%, 2%, 5%, 10%, 25%, and
50% to simulate low-dose settings. All the images were reconstructed using
vendors’ software from United Imaging to simulate the clinical reality. In the
rest of this paper, instead of “low-dose”, the term “low-count” will be used
to be technically more accurate in a PET setting.

##  3 Results

Table 1: Quantitative assessment for different methods. The measurements were
obtained by averaging the values on the testing human studies. The proposed
method consistently produced promising denoised results regardless of input
count levels. $p<0.01$ was observed in all groups when comparing the proposed
method with all other methods (except for 50% SSIM values). The best results
among different low-count levels are marked in red. Black voxels were removed
for calculations.

United Imaging uExplorer Scanner (20 patients $\times$ 6 = 120 studies)
---
PSNR$\uparrow$/NRMSE$\downarrow$/SSIM$\uparrow$ | 1% Count Input | 2% Count Input | 5% Count Input | 10% Count Input | 25% Count Input | 50% Count Input
Input | 44.667 / 0.682 / 0.788 | 49.260 / 0.404 / 0.895 | 53.575 / 0.249 / 0.957 | 56.085 / 0.189 / 0.977 | 59.346 / 0.131 / 0.991 | 62.253 / 0.093 / 0.996
UNN | 52.034 / 0.293 / 0.953 | 54.044 / 0.235 0.973 | 55.838 / 0.194 / 0.982 | 56.963 / 0.173 / 0.987 | 58.570 / 0.146 / 0.991 | 59.885 / 0.126 / 0.994
DDIM50 | 42.550 / 0.852 / 0.904 | 42.718 / 0.836 / 0.926 | 42.846 / 0.824 / 0.938 | 42.942 / 0.815 / 0.944 | 43.065 / 0.804 / 0.952 | 43.145 / 0.797 / 0.957
DDIM50+2.5D | 43.636 / 0.754 0.910 | 43.851 / 0.736 / 0.924 | 43.856 / 0.735 / 0.930 | 43.829 / 0.738 / 0.932 | 43.805 / 0.740 / 0.936 | 43.820 / 0.738 / 0.938
DiffusionMBIR | 42.590 / 0.848 / 0.913 | 42.747 / 0.833 / 0.934 | 42.868 / 0.822 / 0.944 | 42.960 / 0.814 / 0.950 | 43.079 / 0.803 / 0.957 | 43.156 / 0.796 / 0.961
DiffusionMBIR+2.5D | 42.615 / 0.846 / 0.913 | 42.834 / 0.825 / 0.935 | 42.892 / 0.819 / 0.945 | 42.911 / 0.818 / 0.948 | 42.931 / 0.816 / 0.952 | 42.950 / 0.814 / 0.955
TPDM | 42.691 / 0.839 / 0.906 | 42.810 / 0.827 / 0.929 | 42.939 / 0.815 / 0.940 | 43.051 / 0.805 / 0.946 | 43.202 / 0.791 / 0.954 | 43.291 / 0.783 / 0.959
TPDM+2.5D | 42.334 / 0.873 / 0.887 | 42.458 / 0.861 / 0.912 | 42.483 / 0.859 / 0.923 | 42.499 / 0.857 / 0.927 | 42.528 / 0.854 / 0.932 | 42.536 / 0.853 / 0.935
DDPET-3D (proposed) | 52.899 / 0.267 / 0.965 | 54.937 / 0.215 / 0.977 | 57.119 / 0.171 / 0.985 | 58.551 / 0.148 / 0.989 | 60.916 / 0.117 / 0.993 | 63.804 / 0.088 / 0.996

Table 2: Quantitative assessment for different methods on all the testing
patients from United Imaging scanner. The measurements were obtained by
averaging the values on the testing human studies. The proposed method
consistently produced promising denoised results regardless of input count
levels. $p<0.01$ was observed in all groups when comparing the proposed method
with all other methods (except for 50% SSIM values). The best results among
different low-count levels are marked in red. Black voxels were removed for
calculations.

United Imaging uExplorer Scanner (100 patients $\times$ 6 = 600 studies)
---
PSNR$\uparrow$/NRMSE$\downarrow$/SSIM$\uparrow$ | 1% Count Input | 2% Count Input | 5% Count Input | 10% Count Input | 25% Count Input | 50% Count Input
Input | 43.141 / 0.706 / 0.780 | 47.727 / 0.418 / 0.889 | 52.052 / 0.254 / 0.954 | 54.535 / 0.191 / 0.975 | 57.481 / 0.135 / 0.990 | 59.805 / 0.103 / 0.995
UNN | 50.369 / 0.304 / 0.950 | 52.441 / 0.239 / 0.970 | 54.013 / 0.201 / 0.979 | 55.063 / 0.179 / 0.984 | 56.298 / 0.155 / 0.989 | 57.327 / 0.138 / 0.992
DDPET-3D (proposed) | 51.629 / 0.263 / 0.962 | 53.776 / 0.207 / 0.976 | 55.911 / 0.163 / 0.984 | 57.370 / 0.140 / 0.988 | 59.690 / 0.108 / 0.993 | 62.516 / 0.080 / 0.997

Figure 2: Low-count denoised images reconstructed using different methods.
Note that all the images were normalized by the total injected activities.
“2.5D” indicates multi-slice condition. Blue arrows point to a possible lesion
that is better recovered by DDPET-3D. Compared to other baseline methods,
DDPET-3D produced images with better details and maintain consistency along
the slice dimension. Dashed blue line indicate the location of the transverse
slice.

Figure 3: Low-count denoised images reconstructed using the proposed DDPET-3D
with different low-count levels. 6 different patients at 6 different low-count
levels are shown. DDPET-3D can produce consistent denoised results for a wide
range of different noise-levels.

Figure 4: Images reconstructed using DDPET-3D with different numbers of
conditional slices. Note that all the images were normalized by the total
injected activities. Blue arrows point to a possible lesion that is better
recovered by DDPET-3D with more conditional slices.

Figure 5: Visual comparison of images reconstructed using the proposed
DDPET-3D methods with and without denoised priors. Note that all the images
were normalized by the total injected activities. Note the tracer activity
differences in different organs. For the brain transverse slice, due to
quantification bias, image intensity in the ”DDPET-3D (no prior)” method was
manually tuned for visual comparison. Many image details were lost with no
prior. Dashed blue line indicates the location of the transverse slice.

Figure 6: Images reconstructed with different settings of noise variables.
Note that all the images were normalized by the total injected activities.
Note the undesired artifacts along the z-direction by using only one noise
variable $\epsilon$. Also note the inconsistencies of difference slices
without fixing $\epsilon$. Despite similar quantitative values, proposed
network settings produced images with more consistent 3D reconstructions and
less artifacts.

###  3.1 Baseline Comparison

The proposed method was compared with DiffusionMBIR [3], and TPDM [14]. Both
of them were proposed to address the diffusion 3D inconsistencies problem in
medical imaging. DiffusionMBIR applies a TV penalty term along the z-axis, and
TPDM samples the images using 2 pre-trained perpendicular diffusion models. We
also compared with the UNN method [29], which was proposed for noise-aware PET
denoising. Lastly, we compared with the standard DDIM sampling.

We extended the standard DDIM sampling method using conditional information
from neighboring slices (denoted as DDIM50+2.5D). For sampling strategy, we
used DDIM sampling for all diffusion models. Since our proposed method used 2
noise variables, for a fair comparison, all the comparison diffusion model
used 50 sampling steps, while our proposed method only used 25 sampling steps
(note that the proposed method runs 26 steps with 2 noise variables, since it
only averages the 2 outputs from the first reverse step).

As presented in Fig. 2, even though we used neighboring slices as conditional
information for DDIM sampling, both DiffusionMBIR and TPDM still produced more
consistent reconstructions along the z-axis. However, the inconsistent issue
still exist in both images. Proposed DDPET-3D method produced noticeable more
consistent reconstructions by simply fixing the noise variables in the reverse
process. Moreover, in the transverse plane, despite the promising visual
results, all other comparison diffusion models produced images with distorted
organs (e.g., myocardium in the transverse slice), which is unacceptable in
clinical settings. The U-net-based method (UNN), on the other hand, still
maintained the overall structure of the heart.

Another major issue of diffusion model is inaccurate image quantification.
Note that all the presented images were already normalized by the total
injected activities. As shown in Fig. 2, the tracer activities in certain
organs are completely wrong in other comparison diffusion models. For example,
the activities in the brain are noticeably lower in DDIM, DiffusionMBIR, and
TPDM results. Such difference may affect certain diagnostic tasks such as
lesion detection [22]. Even though UNN produced over-smoothed reconstructions,
it does not alter the overall tracer activities in different organs. Using UNN
output as denoised prior, the proposed DDPET-3D maintained overall image
quantification. Also, using UNN as a denoised prior allows the proposed method
to recover some subtle features that are almost invisible in the low-count
input (blue arrows in Fig. 2).

Images were quantitatively evaluated using SSIM (Structural Similarity Index),
PSNR (Peak signal-to-noise ratio), and RMSE (Root-mean-square error). To
facilitate the testing process, we tested all the comparison methods using 20
patients from the entire testing dataset. There are 6 different low-count
levels of all the patients, resulting in a total of 120 testing studies.
Quantitative results are presented in Table 1. We also extended both
DiffusionMBIR and TPDM, and re-trained them using neighboring slices as
conditional information (denoted as DiffusionMBIR+2.5D and TPDM+2.5D in Table
1).

As presented in Table 1, the proposed method outperformed other methods in all
the 6 low-count levels. When the input count-level increases, the performance
of all the methods gradually improves. Both the proposed method and UNN can
achieve noise-aware denoising. However, at higher count-levels (25% and 50%),
UNN results were even worse than the input. In contrast, the proposed method
produced optimal results at all count-levels. Also, due to inaccurate image
quantification, other diffusion models produced images with even worse
quantitative results compared with low-count inputs at higher count-levels.

It is worth noting that, simply adding neighboring slices as conditional
information does not necessarily lead to better performance for DiffusionMBIR
and TPDM. Visual comparisons for DiffusionMBIR+2.5D, TPDM+2.5D, and other
methods are presented in Supplemental Fig. 7.

As presented in Table 1, UNN is the second-best method. To comprehensively
evaluate the proposed method, We compared the proposed method with UNN for all
the 100 testing patients. As presented in Table 2, DDPET-3D consistently
outperformed UNN on different low-count levels.

To demonstrate that DDPET-3D achieves dose-aware denoising, denoised results
using inputs with different low-count levels are presented in Fig. 3. DDPET-3D
consistently produced promising denoised results regardless of input count-
levels. The reconstruction results gradually improve when the input count-
level increases. Some subtle structures are better visualized in denoised
results with higher count-levels input images.

###  3.2 Ablation Studies

We performed several ablated experiments to demonstrate the effectiveness of
different proposed components in DDPET-3D.

Impact of Number of Conditional Slices: We evaluated the DDPET-3D with
different numbers of neighboring conditional slices. Three variants of
DDPET-3D networks were trained using 9, 21, and 41 neighboring slices as
conditional information. These networks are denoted as ”DDPET-3D ($n$ slice)”,
where $n$ is the number of conditioned neighboring slices. Results showed that
$n=31$ performed the best in most quantitative metrics. Therefore, $n=31$ was
used in this paper. Quantitative results are presented in Supplemental Table
3. As presented in Fig. 4, using more conditioned slices helped the network to
recover some subtle details in the images (blue arrows in Fig. 4).

Impact of Denoised Prior: We tested the DDPET-3D without the denoised prior
during the sampling process to demonstrate the improvement in image
quantification using this denoised prior. This network is denoted as ”DDPET-3D
(no prior)” in Supplemental Table 3. As presented in Fig. 5, we can see that
diffusion models produced images with inaccurate tracer activities in
different organs without the proposed denoised prior (especially in the brain
and liver). Note that the images were already normalized by the total injected
activities of the entire 3D volume. Also, without the denoised prior, many
details in the images were not able to recovered, as presented in the brain
slice in Fig. 5.

Impact of Fixing Noise Variables: We analyzed the DDPET-3D with and without
fixing the 2 noise variables to demonstrate its effectiveness in producing
consistent 3D reconstructions. Specifically, $\epsilon_{0}^{a}$ and
$\epsilon_{0}^{b}$ are sampled from ${\mathcal{N}}(0,\textbf{I})$ for every
slice in the 3D image volume, instead of fixing them for all slices. This
network is denoted as ”DDPET-3D (no fix $\epsilon$)” in Supplemental Table 3.
Although the differences in quantitative evaluations with and without fixing
noise variables are small, and not fixing noise even led to better
quantitative measurements in certain cases, we noticed significant
improvements in visual quality, as presented in Fig. 6. Not fixing noise
variables produced images with inconsistent slices and unclear organ
boundaries, which are unfavorable in clinical settings.

Impact of Using Multiple Noise Variables: We tested the DDPET-3D method with
only one noise variable. Specifically, we only initialized one noise variable
instead of 2 ($\epsilon_{0}^{a}$ and $\epsilon_{0}^{b}$). This method is
denoted as ”DDPET-3D (single $\epsilon$)” in Supplemental Table 3. Using 2
noise variables produced images with better quantitative assessments. Also,
Fig. 6 shows that using only one noise variable would result in undesired
artifacts along the z-direction.

Impact of Dose Embedding: We tested the DDPET-3D without dose embedding, which
is denoted as ”DDPET-3D (no dose)” in Supplemental Table 3. Experimental
results showed that DDPET-3D with dose embedding produced images with better
quantitative results.

Generalizability Test: To evaluate DDPET-3D’s generalizability, we applied the
trained model to patient data acquired using a different scanner from a
different hospital in another country. 20 patient studies were acquired using
the Siemens Vision Quadra scanner at XX in this experiment. Results are shown
in Supplemental Fig. 11 and Supplemental Table 4.

##  4 Discussion and Conclusion

We proposed DDPET-3D, a dose-aware diffusion model for 3D ultra-low-dose PET
imaging. Compared with previous methods, DDPET-3D has the following
contributions: _(1)_ 3D reconstructions with fine details and address the 3D
inconsistency issue in previous diffusion models. This was achieved by using
multiple neighboring slices as conditional information and enforcing the same
Gaussian latent for all slices in sampling. _(2)_ DDPET-3D maintains accurate
image quantification by using a pre-trained denoised prior in sampling. _(3)_
DDPET-3D achieves dose-aware denoising. It can be generalized to different
low-count levels, ranging from 1% to 50%. _(4)_ With all the proposed
strategies, DDPET-3D converges within 25 diffusion steps, allowing fast
reconstructions. It takes roughly 15 minutes to reconstruct the entire 3D
volume, while DDPM sampling takes about 6 hours. The experimental results
showed that DDPET-3D achieved SOTA performance compared with previous methods.
In addition, for the first time, we showed that diffusion model can gain
promising results on 1% ultra-low-dose PET denoising problem, both
quantitatively and qualitatively. While this study focuses on PET image
denoising, we believe that the proposed method could be easily extended to
other 3D reconstruction tasks for different imaging modalities.

## References

  * [1] Ronald Boellaard.  Standards for PET image acquisition and quantitative data analysis.  _Journal of Nuclear Medicine_ , 50:11S–20S.
  * [2] Hyungjin Chung and Jong Chul Ye.  Score-based diffusion models for accelerated MRI.  _Medical Image Analysis_ , 80:102479.
  * Chung et al. [2023] Hyungjin Chung, Dohoon Ryu, Michael T. McCann, Marc L. Klasky, and Jong Chul Ye.  Solving 3d inverse problems using pre-trained 2d diffusion models.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 22542–22551, 2023.
  * Clark et al. [2012] Christopher M Clark, Michael J Pontecorvo, Thomas G Beach, Barry J Bedell, R Edward Coleman, P Murali Doraiswamy, Adam S Fleisher, Eric M Reiman, Marwan N Sabbagh, Carl H Sadowsky, et al.  Cerebral pet with florbetapir compared with neuropathology at autopsy for detection of neuritic amyloid-$\beta$ plaques: a prospective cohort study.  _The Lancet Neurology_ , 11(8):669–678, 2012.
  * Croitoru et al. [2023] Florinel-Alin Croitoru, Vlad Hondru, Radu Tudor Ionescu, and Mubarak Shah.  Diffusion models in vision: A survey.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2023.
  * Dhariwal and Nichol [2021] Prafulla Dhariwal and Alexander Nichol.  Diffusion models beat gans on image synthesis.  _Advances in neural information processing systems_ , 34:8780–8794, 2021.
  * [7] Qi Gao, Zilong Li, Junping Zhang, Yi Zhang, and Hongming Shan.  CoreDiff: Contextual error-modulated generalized diffusion model for low-dose CT denoising and generalization.  _IEEE Transactions on Medical Imaging_ , pages 1–1.
  * [8] Kuang Gong, Keith Johnson, Georges El Fakhri, Quanzheng Li, and Tinsu Pan.  PET image denoising based on denoising diffusion probabilistic model.  _European Journal of Nuclear Medicine and Molecular Imaging_.
  * [9] Alper Güngör, Salman UH Dar, Şaban Öztürk, Yilmaz Korkmaz, Hasan A. Bedel, Gokberk Elmas, Muzaffer Ozbey, and Tolga Çukur.  Adaptive diffusion priors for accelerated MRI reconstruction.  _Medical Image Analysis_ , 88:102872.
  * Ho et al. [2020] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  _arXiv preprint arxiv:2006.11239_ , 2020.
  * Hudson and Larkin [1994] H.M. Hudson and R.S. Larkin.  Accelerated image reconstruction using ordered subsets of projection data.  _IEEE Transactions on Medical Imaging_ , 13(4):601–609, 1994.
  * [12] Caiwen Jiang, Yongsheng Pan, Mianxin Liu, Lei Ma, Xiao Zhang, Jiameng Liu, Xiaosong Xiong, and Dinggang Shen.  PET-diffusion: Unsupervised PET enhancement based on the latent diffusion model.  In _Medical Image Computing and Computer Assisted Intervention – MICCAI 2023_ , pages 3–12. Springer Nature Switzerland.
  * Kazerouni et al. [2022] Amirhossein Kazerouni, Ehsan Khodapanah Aghdam, Moein Heidari, Reza Azad, Mohsen Fayyaz, Ilker Hacihaliloglu, and Dorit Merhof.  Diffusion models for medical image analysis: A comprehensive survey.  _arXiv preprint arXiv:2211.07804_ , 2022.
  * Lee et al. [2023] Suhyeon Lee, Hyungjin Chung, Minyoung Park, Jonghyuk Park, Wi-Sun Ryu, and Jong Chul Ye.  Improving 3d imaging with pre-trained perpendicular 2d diffusion models, 2023.
  * Liu et al. [2022] Qiong Liu, Hui Liu, Niloufar Mirian, Sijin Ren, Varsha Viswanath, Joel Karp, Suleman Surti, and Chi Liu.  A personalized deep learning denoising strategy for low-count pet images.  _Physics in Medicine & Biology_, 67(14):145014, 2022.
  * Nichol and Dhariwal [2021] Alexander Quinn Nichol and Prafulla Dhariwal.  Improved denoising diffusion probabilistic models.  In _International Conference on Machine Learning_ , pages 8162–8171. PMLR, 2021.
  * Ouyang et al. [2019] Jiahong Ouyang, Kevin T. Chen, Enhao Gong, John Pauly, and Greg Zaharchuk.  Ultra-low-dose PET reconstruction using generative adversarial network with feature matching and task-specific perceptual loss.  _Medical Physics_ , 46(8):3555–3564, 2019.
  * Robbins [2008] Elizabeth Robbins.  Radiation risks from imaging studies in children with cancer.  _Pediatric blood & cancer_, 51(4):453–457, 2008.
  * Rohren et al. [2004] Eric M Rohren, Timothy G Turkington, and R Edward Coleman.  Clinical applications of pet in oncology.  _Radiology_ , 231(2):305–332, 2004.
  * Rombach et al. [2022] Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer.  High-resolution image synthesis with latent diffusion models.  _arXiv preprint arXiv:2112.10752_ , 2022.
  * Ronneberger et al. [2015] O Ronneberger, P Fischer, and T Brox.  U-Net: Convolutional Networks for Biomedical Image Segmentation.  In _Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015_ , pages 234–241. Springer International Publishing, 2015.
  * Schaefferkoetter et al. [2015] Joshua D Schaefferkoetter, Jianhua Yan, David W Townsend, and Maurizio Conti.  Initial assessment of image quality for low-dose pet: evaluation of lesion detectability.  _Physics in Medicine & Biology_, 60(14):5543, 2015.
  * Schwaiger et al. [2005] Markus Schwaiger, Sibylle Ziegler, and Stephan G Nekolla.  Pet/ct: challenge for nuclear cardiology.  _Journal of Nuclear Medicine_ , 46(10):1664–1678, 2005.
  * Shepp and Vardi [1982] LA Shepp and Y Vardi.  Maximum Likelihood Reconstruction for Emission Tomography.  _IEEE Transactions on Medical Imaging_ , 1(2):113–122, 1982.
  * Song et al. [2022] Jiaming Song, Chenlin Meng, and Stefano Ermon.  Denoising diffusion implicit models.  _arXiv preprint arXiv:2010.02502_ , 2022.
  * Song et al. [2021] Yang Song, Jascha Sohl-Dickstein, Diederik P. Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  _arXiv preprint arXiv:2011.13456_ , 2021.
  * Wang et al. [2014] Chenye Wang, Zhenghui Hu, Pengcheng Shi, and Huafeng Liu.  Low dose PET reconstruction with total variation regularization.  In _2014 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society_ , pages 1917–1920, 2014.  ISSN: 1558-4615.
  * Wang [2016] G. Wang.  A Perspective on Deep Imaging.  _IEEE Access_ , 4:8914–8924, 2016.
  * Xie et al. [2023] Huidong Xie, Qiong Liu, Bo Zhou, Xiongchao Chen, Xueqi Guo, and Chi Liu.  Unified noise-aware network for low-count PET denoising.  _arXiv preprint arXiv:2304.14900_ , 2023.
  * Xu et al. [2017] Junshen Xu, Enhao Gong, John Pauly, and Greg Zaharchuk.  200x low-dose PET reconstruction using deep learning.  _arXiv preprint arXiv:1712.04119_ , 2017.
  * Zhou et al. [2021] Bo Zhou, Yu-Jung Tsai, Xiongchao Chen, James S Duncan, and Chi Liu.  Mdpet: A unified motion correction and denoising adversarial network for low-dose gated pet.  _IEEE Transactions on Medical Imaging_ , 2021.
  * Zhou et al. [2022] Bo Zhou, Tianshun Miao, Niloufar Mirian, Xiongchao Chen, Huidong Xie, Zhicheng Feng, Xueqi Guo, Xiaoxiao Li, S Kevin Zhou, James S Duncan, et al.  Federated transfer learning for low-dose pet denoising: A pilot study with simulated heterogeneous data.  _IEEE Transactions on Radiation and Plasma Medical Sciences_ , 2022.
  * Zhou et al. [2023a] Bo Zhou, Yu-Jung Tsai, Jiazhen Zhang, Xueqi Guo, Huidong Xie, Xiongchao Chen, Tianshun Miao, Yihuan Lu, James S Duncan, and Chi Liu.  Fast-mc-pet: A novel deep learning-aided motion correction and reconstruction framework for accelerated pet.  _arXiv preprint arXiv:2302.07135_ , 2023a.
  * Zhou et al. [2023b] Bo Zhou, Huidong Xie, Qiong Liu, Xiongchao Chen, Xueqi Guo, Zhicheng Feng, Jun Hou, S Kevin Zhou, Biao Li, Axel Rominger, et al.  Fedftn: Personalized federated learning with deep feature transformation network for multi-institutional low-count pet denoising.  _Medical image analysis_ , 90:102993, 2023b.
  * Zhou et al. [2020] Long Zhou, Joshua D Schaefferkoetter, Ivan WK Tham, Gang Huang, and Jianhua Yan.  Supervised learning with cyclegan for low-dose fdg pet image denoising.  _Medical image analysis_ , 65:101770, 2020.

Supplementary Material

Fig. 7: we extended DiffusionMBIR and TPDM methods using neighboring slices as
conditional information (denoted as DiffusionMBIR+2.5D and TPDM+2.5D).

Fig. 9: we added coronal view in addition to Fig. 2 in the main text.

Fig. 8: Images reconstructed using DDPET-3D with inputs with different low-
count levels.

Fig. 10: Same patients presented in Fig. 3 with larger field-of-view.

Fig. 11: to demonstrate that the proposed DDPET-3D method has superior
generalizability, we directly apply the trained model on patient data acquired
on a different scanner from another hospital. 20 patient studies were used for
this experiment.

Table 3: quantitative assessments of different ablation studies described in
the main text.

Table 4: quantitative assessments of the generalizability test described in
the main text.

Figure 7: Visual comparison of low-count denoised images reconstructed using
different methods. Note that all the images were normalized by the total
injected activities. ”2.5D” indicates multi-slice condition.

Figure 8: Visual comparison of low-count denoised images reconstructed using
the proposed DDPET-3D method with different low-count levels. Blue arrows
point to subtle details that are gradually better recovered with higher-count
inputs. DDPET-3D can produce consistent denoised results for a wide range of
different noise-levels.

Figure 9: Visual comparison of low-count denoised images reconstructed using
different methods. Note that all the images were normalized by the total
injected activities. “2.5D” indicates multi-slice condition. Blue arrows point
to subtle details that are better recovered by the proposed method. It is
worth highlighting that, compared to other baseline methods, DDPET can
reconstruct images with fine-details and maintain consistency along the slice
dimension.

Figure 10: Low-count denoised images reconstructed using the proposed DDPET-3D
method with different low-count levels. DDPET-3D can produce consistent
denoised results for a wide range of different noise-levels.

Figure 11: Patient data acquired on a Siemens Vision Quadra scanner at Country
YY reconstructed using the proposed DDPET-3D. Specifically, the DDPET-3D was
trained using data acquired on the United Imaging uExplorer scanner at Country
XX. The presented results demonstrate the strong generalizability of the
proposed method. As pointed out by the blue arrow, the aorta can be recovered
from the 1% low-count input. Table 3: Quantitative assessment for different
ablation studies. The best results among different low-count levels are marked
in red. The second best results are marked in blue. Black voxels were removed
for calculations.

United Imaging uExplorer Scanner (20 patients $\times$ 6 = 120 studies)
---
PSNR$\uparrow$/NRMSE$\downarrow$/SSIM$\uparrow$ | 1% Count Input | 2% Count Input | 5% Count Input | 10% Count Input | 25% Count Input | 50% Count Input
DDPET-3D (9 slices) | 52.563 / 0.277 / 0.964 |  54.666 / 0.215 / 0.977 |  56.194 / 0.162 / 0.984 |  57.688 / 0.137 / 0.988 |  60.011 / 0.106 / 0.993 |  62.267 / 0.082 / 0.995
DDPET-3D (21 slices) | 52.502 / 0.278 / 0.960 | 54.509 / 0.224 / 0.974 | 56.574 / 0.181 / 0.981 | 57.949 / 0.157 / 0.985 | 60.137 / 0.125 / 0.990 | 62.338 / 0.100 / 0.994
DDPET-3D (41 slices) |  51.866 / 0.262 / 0.964 |  53.926 / 0.208 / 0.977 |  56.147 / 0.162 / 0.985 |  57.618 / 0.138 / 0.989 |  59.948 / 0.107 / 0.993 |  62.438 / 0.081 / 0.995
DDPET-3D (no prior) | 41.497 / 0.961 / 0.868 | 41.638 / 0.945 / 0.898 | 41.704 / 0.938 / 0.913 | 41.715 / 0.937 / 0.920 | 41.747 / 0.934 / 0.927 | 41.767 / 0.932 / 0.932
DDPET-3D (no fix $\epsilon$) |  52.991 / 0.264 / 0.961 |  55.049 / 0.212 / 0.973 |  57.225 / 0.170 / 0.980 |  58.695 / 0.146 / 0.983 |  61.026 / 0.116 / 0.988 |  63.916 / 0.087 / 0.991
DDPET-3D (no dose) | 52.817 / 0.269 / 0.960 | 54.822 / 0.217 / 0.972 | 56.903 / 0.175 / 0.979 | 58.259 / 0.152 / 0.983 | 60.360 / 0.123 / 0.987 | 62.871 / 0.096 / 0.990
DDPET-3D (single $\epsilon$) | 52.289 / 0.285 / 0.959 | 54.261 / 0.231 / 0.974 | 56.324 / 0.185 / 0.983 | 57.701 / 0.161 / 0.987 | 59.917 / 0.128 / 0.992 | 62.548 / 0.098 / 0.995
DDPET-3D (proposed) |  52.899 / 0.267 / 0.965 |  54.937 / 0.215 / 0.977 |  57.119 / 0.171 / 0.985 |  58.551 / 0.148 / 0.989 |  60.916 / 0.117 / 0.993 |  63.804 / 0.088 / 0.996

Table 4: Quantitative assessment for the generalizability test. The model was
trained using the United Imaging uExplorere scanner at Country XX and directly
applied to 20 patient studies acquired using a Siemens Vision Quadra scanner
at Country YY. The best results among different low-count levels are marked in
red. 50% count images were not available for these patient studies. Black
voxels were removed for calculations.

Siemens Vision Quadra Scanner (20 patients $\times$ 5 = 100 studies)
---
PSNR$\uparrow$/NRMSE$\downarrow$/SSIM$\uparrow$ | 1% Count Input | 2% Count Input | 5% Count Input | 10% Count Input | 25% Count Input
Input | 50.964 / 0.385 / 0.904 | 53.737 / 0.275 / 0.946 | 56.993 / 0.188 / 0.976 | 59.273 / 0.145 / 0.987 | 62.708 / 0.098 / 0.995
DDPET-3D (proposed) | 56.818 / 0.192 / 0.981 | 58.130 / 0.165 / 0.986 | 59.604 / 0.139 / 0.991 | 60.692 / 0.123 / 0.993 | 62.850 / 0.096 / 0.996