Benoit Dherin††thanks: equal contribution<EMAIL_ADDRESS>
Google, Dublin and Michael Munn11footnotemark: 1<EMAIL_ADDRESS>
Google, New York and David G.T. Barrett<EMAIL_ADDRESS>
DeepMind, London

# The Geometric Occam’s Razor Implicit in Deep Learning

###### Abstract

In over-parameterized deep neural networks there can be many possible
parameter configurations that fit the training data exactly. However, the
properties of these interpolating solutions are poorly understood. We argue
that over-parameterized neural networks trained with stochastic gradient
descent are subject to a Geometric Occam’s Razor; that is, these networks are
implicitly regularized by the geometric model complexity. For one-dimensional
regression, the geometric model complexity is simply given by the arc length
of the function. For higher-dimensional settings, the geometric model
complexity depends on the Dirichlet energy of the function. We explore the
relationship between this Geometric Occam’s Razor, the Dirichlet energy and
other known forms of implicit regularization. Finally, for ResNets trained on
CIFAR-10, we observe that Dirichlet energy measurements are consistent with
the action of this implicit Geometric Occam’s Razor.

##  1 Introduction

Naively, we might expect that over-parameterized models will overfit the
training data and that under-parameterized models will be better since they
have fewer degrees of freedom. However, it turns out that over-parameterized
models can find better solutions than the under-parameterized models - a
paradoxical phenomenon known as the double-descent curve [Belkin (2021),
Nakkiran et al. (2020)]. One possible explanation for this behaviour is that
over-parameterized models are subject to an Occam’s razor that filters out
unnecessarily complex solutions in favour of simpler solutions.

Typically we might expect an Occam’s razor to take the form of a complexity
measure on the number of model parameters or the size of the hypothesis space,
for instance [Zhang et al. (2017)]. However, for neural networks, the precise
form of this hypothesized Occam’s razor is not known, since it is not
explicitly enforced during training. There has been some progress recently to
identify sources of implicit regularization that may play a role here [Barrett
and Dherin (2021), Smith et al. (2021), Ma and Ying (2021), Blanc et al.
(2020)]. For instance, recent work has exposed a hidden from of regularization
in Stochastic Gradient Descent (SGD) called Implicit Gradient Regularization
(IGR) [Barrett and Dherin (2021), Smith et al. (2021)] which penalizes
learning trajectories that have large loss gradients.

For over-parameterized neural networks trained with SGD, we hypothesize that
the hidden Occam’s razor takes the form of a geometric complexity measure. Our
key contributions are as follows: (1) define this notion of geometric
complexity; (2) show that the Dirichlet energy can be used as a proxy for
geometric complexity; (3) show that the IGR mechanism from SGD puts a
regularization pressure on the geometric complexity; and (4) show that the
strength of this pressure increases with the size of the learning rate, which
we verify with numerical experiments.

##  2 The Geometric Occam’s razor in 1-dimensional regression

To build intuition, we begin with a simple 1-dimensional example. Consider a
ReLU neural network consisting of 3 layers with 300 units per layer, trained
using SGD, without any form of explicit regularization to perform
1-dimensional regression using only 10 data points. In this extreme setting,
we should expect the network to overfit the dataset, since the function space
described by that neural network is extremely large - consisting of piecewise
linear functions with thousands of linear pieces (Arora et al., 2018). Yet, if
we plot the learned function during training from the first step all the way
up to interpolation, as in Figure 1, we observe that the learned function is
the ‘simplest’ possible function, in some sense, among all functions with the
same training error.

But what do we mean by ‘simple’? Our key intuition in this example is that the
arc length of the learned function over the smallest interval containing the
data points provides our measure of model complexity. At the end of training,
we see that the arc length of the learned function is close to that of the
shortest possible path interpolating between the data points, suggesting that
this measure of geometric complexity is somehow optimised during training.

Figure 1: Training sequence for an over-parameterised neural network trained
using 10 data points plotted at iteration step 1, 20, 40, 50, 5000, 10000,
20000, and 30000 (from top left to bottom right, respectively). The blue
circles represent the data points while the red crosses represent the points
predicted by the neural network.

##  3 Dirichlet energy as a measure of function complexity

In the previous section, we used the arc length of the learnt 1-dimensional
function as a measure of its geometric complexity. What is the corresponding
notion for a function in a high-dimensional feature space $f:\rm
I\\!R^{d}\rightarrow\rm I\\!R$? In this case, we can define the geometric
complexity of a function $f$ as the volume $\Omega_{D}(f)$ of its graph

$\textrm{gr}(f_{X_{D}}):=\\{(x,f(x)):\>x\in X_{D}\\}\subset\rm
I\\!R^{d}\times\rm I\\!R$ (1)

restricted to the feature polytope $X_{D}$; that is, the polytope with
smallest volume containing all the feature points $x_{i}$ of the dataset
$D=\\{(x_{i},y_{i}):i=1,\dots,n\\}$. From differential geometry (do Carmo,
1976), for a smooth function $f:\rm I\\!R^{d}\rightarrow\rm I\\!R$, the graph
$\textrm{gr}(f_{X_{D}})$ is an $n$-dimensional smooth submanifold of $\rm
I\\!R^{n+1}$. Using the Riemannian metric on $\textrm{gr}(f_{X_{D}})$ induced
from the Euclidean metric on $\rm I\\!R^{n+1}$ and its corresponding
Riemannian volume form, the volume of the graph of $f$ can be expressed as

$\Omega_{D}(f)=\int_{X_{D}}\sqrt{1+\|\nabla_{x}f\|^{2}}~{}dx.$ (2)

This can in turn be approximated using a first-order Taylor series expansion
$\sqrt{1+z}\simeq 1+\frac{1}{2}z$ so that

$\Omega_{D}(f)\approx\int_{X_{D}}1+\frac{1}{2}\|\nabla_{x}f\|^{2}~{}dx=\textrm{Vol}(X_{D})+\textrm{DE}(f),$
(3)

where

$\textrm{DE}(f)=\frac{1}{2}\int_{X_{D}}\|\nabla_{x}f\|^{2}dx.$ (4)

is the Dirichlet energy  of the function $f$ over $X_{D}$. The computation
above suggests that both the function volume or its Dirichlet energy can be
used as a measure of a function’s geometric complexity.

One way to compute the Dirichlet energy numerically is to use a quadrature
formula summing up $\|\nabla_{x}f\|^{2}$ over a number of points in $X_{D}$
and multiplying the summands by the volume element of the point. So, if we use
the data points themselves for evaluation and $1/|D|$ as a proxy for the
volume element, we obtain a discrete version of the Dirichlet energy which we
call the discrete Dirichlet energy, denoted by $\widehat{\textrm{DE}}(f)$.
This provides an easily computable measure of a function’s geometric
complexity:

$\widehat{\textrm{DE}}(f)=\frac{1}{2|D|}\sum_{x\in D}\|\nabla_{x}f(x)\|^{2}.$
(5)

##  4 How neural networks tame model complexity

We now argue that the geometric complexity, as measured by the discrete
Dirichlet energy, is implicitly regularized during the training of neural nets
with vanilla SGD.

In recent work Smith et al. (2021), it was shown that the discrete steps of
SGD from epoch to epoch closely follow, on average, the gradient flow of a
modified loss of the form:

$\tilde{L}=L+\frac{h}{4|D|}\sum_{(x,y)\in
D}\|\nabla_{\theta}L(x,\theta)\|^{2},$

where $L$ is the original loss and $L(x,\theta):=E(f_{\theta}(x),y)$ is the
error $E$ between the prediction $f_{\theta}(x)$ and the true label $y$. This
means that during SGD the quantities $\|\nabla_{\theta}L(x,\theta)\|^{2}$ at
each data point $(x,y)\in D$, are implicitly regularized, with the learning
rate $h$ acting as an implicit regularization rate.

Now, for models whose losses come from the application of a maximum likelihood
estimation on a conditional probability distribution in the exponential family
such as the least-square loss or the cross-entropy loss, we obtain loss
gradients that have the following form:

$\nabla_{\theta}L(x,\theta)=\epsilon_{x}(\theta)\nabla_{\theta}f_{\theta}(x),$

where $\epsilon_{x}(\theta)=(f_{\theta}(x)-y)$ is the signed residual,
yielding

$\tilde{L}=L+\frac{h}{4|D|}\sum_{(x,y)\in
D}\epsilon_{x}(\theta)^{2}\|\nabla_{\theta}f_{\theta}(x)\|^{2}.$ (6)

From that last expression for $\tilde{L}$, we see that the terms
$\|\nabla_{\theta}f_{\theta}(x)\|^{2}$ are implicitly regularized at each data
point $(x,y)$ and even more so in the region where the residual errors are
large, such as the beginning of training.

We now argue that for neural networks, in particular, the regularization
pressure on the gradient of the network with respect to the parameters
$\|\nabla_{\theta}f_{\theta}(x)\|^{2}$ acts as a regularization pressure on
the gradient of the network with respect to the input
$\|\nabla_{x}f_{\theta}(x)\|^{2}$. Hence, this creates a pressure for the
Dirichlet energy to be implicitly regularized during training. In fact, this
follows from the fact that for neural networks their derivatives with respect
to the inputs and the parameters can be related as follows (proof in Appendix
A):

######  Theorem 4.1.

Consider a neural network with $l$ layers $f_{\theta}(x)=f_{1}\circ\dots\circ
f_{l}(x),$ where $f_{i}(z)=a_{i}(w_{i}z+b_{i})$ with
$\theta=(w_{1},b_{1},\dots,w_{l},b_{l})$ being the vector of layer weight
matrices $w_{i}$ and biases $b_{i}$ and the $a_{i}$’s are the layer activation
functions. Then we have that

$\|\nabla_{x}f_{\theta}(x)\|^{2}\left(\frac{1+\|h_{1}(x)\|^{2}}{\|w_{1}\|^{2}\|h_{1}^{\prime}(x)\|^{2}}+\cdots+\frac{1+\|h_{l}(x)\|^{2}}{\|w_{l}\|^{2}\|h_{l}^{\prime}(x)\|^{2}}\right)\leq\|\nabla_{\theta}f(x)\|^{2},$
(7)

where $h_{i}(x)$ is the sub-network from input $x$ to layer $i$ and
$\|w_{i}\|$ is the spectral norm of the weight matrix $w_{i}$.

From equation:gradient_relation, we see that the regularization pressure from
IGR translates into a regularization pressure on the discrete Dirichlet energy
when the positive quantities

$A_{i}(\theta,x):=\|w_{i}\|\|h_{i}^{\prime}(x)\|,\quad i=1,\dots,l$ (8)

remain small. Note that this is expected to happen at the beginning of
training when the spectral norms of the layers are close to zero, while they
tend to grow as the training progresses if no spectral regularization Miyato
et al. (2018) is applied. Furthermore, note here that preventing the $A_{i}$’s
from becoming too large during training may be an important consideration
which informs the choice of model architecture and layer regularization.

Experimental evidence: From Equation equation:resdidual_loss, since the
strength of IGR is a function of the learning rate, we should expect an
increased pressure on the Dirichlet energy as a result of Equation
equation:gradient_relation when training with higher learning rates. We verify
this prediction for a ResNet-18 trained to classify CIFAR-10 images. Measuring
the discrete Dirichlet energy at the time of maximal test accuracy for a range
of learning rates, we observe this predicted behaviour, consistent with our
theory; see Figure 2.

Figure 2: For a ResNet-18 trained on CIFAR-10, we observe that the Dirichlet
energy and the loss-surface slope (essentially the IGR term
$\|\nabla_{\theta}L\|^{2}$) decrease as the learning rate increases. Each dot
corresponds to a fully trained model at time of maximal test accuracy.

######  Remark 4.2.

Note also that for linear models (i.e., neural networks with a single linear
layer), the Dirichlet energy of the network coincides in this case with the
L2-norm of the parameters. Therefore, this results recovers the already known
fact that linear models trained with SGD have an inductive bias towards low
L2-norm solutions (see Zhang et al. (2017)). This also points toward the fact
that the Dirichlet energy may be the right generalization of the L2-norm for a
general network.

##  5 Related work, Future directions, and Discussion

Splines and connections to harmonic function theory: The Dirichlet energy
equation:DE is well-known in harmonic function theory (Axler et al., 2013)
where it can be shown using calculus of variations that harmonic functions
subject to a boundary condition minimize the Dirichlet energy over the space
of differentiable functions. This is known as Dirichlet’s principle. The
minimization of the Dirichlet energy itself is also related to the theory of
splines Jin and Montufar (2021). Our work seems to indicate that neural
networks are biased towards (a notion of) harmonic functions with the dataset
acting as the boundary condition. Complexity theory: The notion of geometric
complexity introduced has similarities to the Kolgomorov complexity
Schmidhuber (1997) as well as the minimum description length given in Hinton
and Zemel (1994). Smoothness regularization: The notion of geometric
complexity considered here is related to the notion of smoothness with respect
to the input as discussed in Rosca et al. (2020) as well as to the Sobolev
regularization effect of SGD discussed in Ma and Ying (2021), where
inequalities similar to equation:gradient_relation but involving only the
first layer are considered. In particular, various forms of gradient
penalties, reminiscent of the Dirichlet energy, have been devised to achieve
Lipschitz smoothness Elsayed et al. (2018); Gulrajani et al. (2017); Fedus et
al. (2018); Arbel et al. (2018); Kodali et al. (2018). It has been shown that
the discrete Dirichlet energy (evaluated at the data points) is a powerful
regularizer Hoffman et al. (2020); Varga et al. (2018) and in Rosca et al.
(2020) that it has advantages over other form of smoothness regularization
(such as spectral norm regularization Miyato et al. (2018); Lin et al.
(2021)). Our analysis shows that we can control this form of regularization
cheaply through the learning rate. In image processing, the Dirichlet energy
is also called the Rudin–Osher–Fatemi total variation and it as been
introduced as a powerful explicit regularizer for image denoising; see Rudin
et al. (1992) and Getreuer (2012). It may be possible that these various forms
of smoothness regularization are useful because they provide implicit control
over the model geometric complexity. Regularization through noise: The
discrete Dirichlet energy is reminiscent of the Tikhonov regularizer which is
implicitly regularized with added input noise Bishop (1995). The modified loss
in equation:resdidual_loss is also very reminiscent of the modified loss in
Blanc et al. (2020), which is argued to be implicitly minimized by SGD when a
random white noise is added to the labels. In Section 3 of Mescheder et al.
(2018), it is argued that explicit gradient regularization with respect to
input and noise instance produce similar types of regularization. Altogether,
this suggests that feature noise, label noise, and the optimization scheme all
conspire to implicitly tame the geometric complexity in the case of neural
networks trained with gradient-based optimization schemes. Regularization
through the number of layers: In Equation equation:gradient_relation, one sees
that each layer contributes an additional positive term, increasing the
pressure on the Dirichlet energy. This suggests that the pressure on the model
geometric complexity may increase with the neural network depth in a similar
spirit as Gao and Jojic (2016). Training of GANs: For GANs, explicit gradient
regularization both with respect to the input (Gulrajani et al., 2017; Arbel
et al., 2018; Fedus et al., 2018; Kodali et al., 2018; Miyato et al., 2018)
and the parameters Rosca et al. (2021); Qin et al. (2020); Balduzzi et al.
(2018); Mescheder et al. (2017); Nagarajan and Kolter (2017) has been proven
to be beneficial and related to smoothness. Our main theorem provides a way to
relate gradient penalties with respect to the input and with respect to the
parameters for neural networks, in a way where the spectral norm of the weight
matrices plays a key role. This points toward geometric complexity being a
useful notion to relate and understand these different forms of regularization
(including spectral normalization as in Miyato et al. (2018) and Lin et al.
(2021)).

##  6 Conclusion

In conclusion, we have found that neural networks trained with SGD are subject
to an implicit Geometric Occam’s razor, which selects parameter configurations
that have low geometric complexity ahead of configurations with high geometric
complexity. This geometric complexity is given by the arc length in
1-dimensinal regression; is linearly related to the Dirichlet energy in
higher-dimensional settings; and has many intriguing similarities to other
known quantities, including various forms of implicit and explicit
regularisation and model complexity. More generally, our work develops
promising new theoretical connections between optimization and the geometry of
over-parameterised neural networks.

### Acknowledgments

We would like to thank Mihaela Rosca, Maxim Neumann, Yan Wu, Samuel Smith, and
Soham De for helpful discussion and feedback. We would also like to thank
Patrick Cole and Shakir Mohamed for their support.

## References

  * Arbel et al. (2018) Michael Arbel, Danica J Sutherland, Mikolaj Binkowski, and Arthur Gretton.  On gradient regularizers for mmd gans.  In _Advances in Neural Information Processing Systems_ , volume 31, 2018.
  * Arora et al. (2018) Raman Arora, Amitabh Basu, Poorya Mianjy, and Anirbit Mukherjee.  Understanding deep neural networks with rectified linear units.  In _International Conference on Learning Representations_ , 2018.
  * Axler et al. (2013) Sheldon Axler, Paul Bourdon, and Ramey Wade.  _Harmonic Function Theory_ , volume 137.  Springer, 2013.
  * Balduzzi et al. (2018) David Balduzzi, Sebastien Racaniere, James Martens, Jakob Foerster, Karl Tuyls, and Thore Graepel.  The mechanics of n-player differentiable games.  In _International Conference on Machine Learning_ , 2018.
  * Barrett and Dherin (2021) David G.T. Barrett and Benoit Dherin.  Implicit gradient regularization.  In _International Conference on Learning Representations_ , 2021.
  * Belkin (2021) Mikhail Belkin.  Fit without fear: remarkable mathematical phenomena of deep learning through the prism of interpolation.  2021\.  URL https://arxiv.org/abs/2105.14368.
  * Bishop (1995) Chris M. Bishop.  Training with noise is equivalent to tikhonov regularization.  _Neural Computation_ , 7(1), 1995.
  * Blanc et al. (2020) Guy Blanc, Neha Gupta, Gregory Valiant, and Paul Valiant.  Implicit regularization for deep neural networks driven by an ornstein-uhlenbeck like process.  In _Conference on Learning Theory_ , 2020.
  * do Carmo (1976) Manfredo P. do Carmo.  _Differential geometry of curves and surfaces._ Prentice Hall, 1976.
  * Elsayed et al. (2018) Gamaleldin Elsayed, Dilip Krishnan, Hossein Mobahi, Kevin Regan, and Samy Bengio.  Large margin deep networks for classification.  In _Advances in Neural Information Processing Systems_ , volume 31, 2018.
  * Fedus et al. (2018) William Fedus, Mihaela Rosca, Balaji Lakshminarayanan, Andrew M. Dai, Shakir Mohamed, and Ian J. Goodfellow.  Many paths to equilibrium: Gans do not need to decrease a divergence at every step.  In _International Conference on Learning Representations_ , 2018.
  * Gao and Jojic (2016) Tianxiang Gao and Vladimir Jojic.  Degrees of freedom in deep neural networks.  In _Proceedings of the Thirty-Second Conference on Uncertainty in Artificial Intelligence_ , 2016.
  * Getreuer (2012) Pascal Getreuer.  Rudin-Osher-Fatemi Total Variation Denoising using Split Bregman.  _Image Processing On Line_ , 2012.
  * Gulrajani et al. (2017) Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron C Courville.  Improved training of wasserstein gans.  In _Advances in neural information processing systems_ , 2017.
  * Hinton and Zemel (1994) Geoffrey E Hinton and Richard Zemel.  Autoencoders, minimum description length and helmholtz free energy.  In J. Cowan, G. Tesauro, and J. Alspector, editors, _Advances in Neural Information Processing Systems_ , volume 6. Morgan-Kaufmann, 1994.
  * Hoffman et al. (2020) Judy Hoffman, Daniel A. Roberts, and Sho Yaida.  Robust learning with jacobian regularization.  2020\.  URL https://arxiv.org/abs/1908.02729.
  * Jin and Montufar (2021) Hui Jin and Guido Montufar.  Implicit bias of gradient descent for mean squared error regression with wide neural networks.  2021\.  URL https://arxiv.org/abs/2006.07356.
  * Kodali et al. (2018) Naveen Kodali, James Hays, Jacob Abernethy, and Zsolt Kira.  On convergence and stability of GANs.  In _Advances in Neural Information Processing Systems_ , 2018.
  * Lin et al. (2021) Zinan Lin, Vyas Sekar, and Giulia C. Fanti.  Why spectral normalization stabilizes gans: Analysis and improvements.  In _Advances in neural information processing systems_ , 2021.
  * Ma and Ying (2021) Chao Ma and Lexing Ying.  The sobolev regularization effect of stochastic gradient descent.  2021\.  URL https://arxiv.org/abs/2105.13462.
  * Mescheder et al. (2017) Lars Mescheder, Sebastian Nowozin, and Andreas Geiger.  The numerics of gans.  In _Advances in Neural Information Processing Systems_ , 2017.
  * Mescheder et al. (2018) Lars Mescheder, Andreas Geiger, and Sebastian Nowozin.  Which training methods for GANs do actually converge?  In _Proceedings of the 35th International Conference on Machine Learning_ , 2018.
  * Miyato et al. (2018) Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida.  Spectral normalization for generative adversarial networks.  In _International Conference on Learning Representations_ , 2018.
  * Nagarajan and Kolter (2017) Vaishnavh Nagarajan and J Zico Kolter.  Gradient descent gan optimization is locally stable.  In _Advances in neural information processing systems_ , 2017.
  * Nakkiran et al. (2020) Preetum Nakkiran, Gal Kaplun, Yamini Bansal, Tristan Yang, Boaz Barak, and Ilya Sutskever.  Deep double descent: Where bigger models and more data hurt.  In _International Conference on Learning Representations_ , 2020.
  * Qin et al. (2020) Chongli Qin, Yan Wu, Jost Tobias Springenberg, Andy Brock, Jeff Donahue, Timothy Lillicrap, and Pushmeet Kohli.  Training generative adversarial networks by solving ordinary differential equations.  In _Advances in Neural Information Processing Systems_ , volume 33, 2020.
  * Rosca et al. (2020) Mihaela Rosca, Theophane Weber, Arthur Gretton, and Shakir Mohamed.  A case for new neural network smoothness constraints.  In _Proceedings on ”I Can’t Believe It’s Not Better!” at NeurIPS Workshops_ , volume 137 of _Proceedings of Machine Learning Research_ , 2020\.
  * Rosca et al. (2021) Mihaela Rosca, Yan Wu, Benoit Dherin, and David G.T. Barrett.  Discretization drift in two-player games.  In _Proceedings of the 38th International Conference on Machine Learning_ , volume 139, 2021.
  * Rudin et al. (1992) Leonid Rudin, Stanley Osher, and Emad Fatemi.  Nonlinear total variation based noise removal algorithms.  _Physica D: Nonlinear Phenomena_ , 60(1), 1992.
  * Schmidhuber (1997) J. Schmidhuber.  Discovering neural nets with low kolmogorov complexity and high generalization capability.  _Neural networks : the official journal of the International Neural Network Society_ , 10 5, 1997.
  * Seong et al. (2018) Sihyeon Seong, Yegang Lee, Youngwook Kee, Dongyoon Han, and Junmo Kim.  Towards flatter loss surface via nonmonotonic learning rate scheduling.  In _Proceedings of the Thirty-Fourth Conference on Uncertainty in Artificial Intelligence_ , 2018.
  * Smith et al. (2021) Samuel L Smith, Benoit Dherin, David G.T. Barrett, and Soham De.  On the origin of implicit regularization in stochastic gradient descent.  In _International Conference on Learning Representations_ , 2021.
  * Varga et al. (2018) Dániel Varga, Adrián Csiszárik, and Zsolt Zombori.  Gradient regularization improves accuracy of discriminative models.  2018\.  URL https://arxiv.org/abs/1712.09936.
  * Zhang et al. (2017) Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.  Understanding deep learning requires rethinking generalization.  In _International Conference on Learning Representations_ , 2017.

##  Appendix A Proof of Theorem 1

Consider a neural network with $l$ layers $f_{\theta}(x)=f_{1}\circ\dots\circ
f_{l}(x),$ where $f_{i}(z)=a_{i}(w_{i}z+b_{i})$ with
$\theta=(w_{1},b_{1},\dots,w_{l},b_{l})$ being the vector of layer weight
matrices $w_{i}$ and biases $b_{i}$ and the $a_{i}$’s are the layer activation
functions. We will use the notation $f_{w_{i}}$ or $f_{b_{i}}$ instead of
$f_{\theta}$ when we want to consider $f_{\theta}$ as dependent on the
$w_{i}$’s or $b_{i}$’s only.

For this model structure and following Pythagoras, we have:

$\|\nabla_{\theta}f_{\theta}(x)\|^{2}=\|\nabla_{w_{1}}f_{\theta}(x)\|^{2}+\|\nabla_{b_{1}}f_{\theta}(x)\|^{2}+\cdots+\|\nabla_{w_{l}}f_{\theta}(x)\|^{2}+\|\nabla_{b_{l}}f_{\theta}(x)\|^{2}$
(9)

For each layer $i$, we can rewrite the network function as

$f_{w_{i}}(x)=g_{i}(w_{i}h_{i}(x)+b_{i}),$

where $g_{i}(z)$ consists of the deeper layers above $i$ and $h_{i}(x)$
consists of the shallower layers below $i$. The idea now, inspired from Seong
et al. (2018), is to show that a small perturbation $\delta x$ of the input is
equivalent to a small perturbation $u(\delta x)$ of the weights of layer $i$.
We will use this idea to prove the following two lemmas.

######  Lemma A.1.

In the notation above, for each layer $i$, we have:

$\|\nabla_{x}f_{\theta}(x)\|^{2}\left(\frac{\|h_{i}(x)\|}{\|w_{i}\|\|h_{i}^{\prime}(x)\|}\right)^{2}\leq\|\nabla_{w_{i}}f(x)\|^{2}.$
(10)

######  Proof A.2.

Consider a small perturbation $x+\delta x$ of the input $x$. We start by
showing that we can always find a corresponding perturbation $w_{i}+u(\delta
x)$ of the weight matrix in layer $i$ such that

$f_{w_{i}}(x+\delta x)=f_{w_{i}+u(\delta x)}(x).$ (11)

Namely, because $f_{w_{i}}(x)=g_{i}(w_{i}h_{i}(x)+b_{i})$, to show this, it is
enough to find $u(\delta x)$ such that

$w_{i}(h_{i}(x)+h_{i}^{\prime}(x)\delta x)+b_{i}=(w_{i}+u(\delta
x))h_{i}(x)+b_{i},$ (12)

where we identify $h_{i}(x+\delta x)$ with its linear approximation around $x$
for small $\delta x$. Then equation:condition is always satisfied if we set

$u(\delta x):=\frac{(w_{i}h_{i}^{\prime}(x)\delta
x)h_{i}^{T}(x)}{\|h_{i}(x)\|^{2}},$ (13)

since $h_{i}^{T}(x)h_{i}(x)=\|h_{i}(x)\|^{2}$. Now taking the derivative with
respect to $\delta x$ at $\delta x=0$ on both sides of Equation
equation:equivalence, and using the chain rule and that $u(\delta x)$ is
linear in $\delta x$, we obtain a relation between the network derivative with
respect to the weight matrices and w.r.t the input:

$\nabla_{x}f_{\theta}(x)=\nabla_{w_{i}}f(x)\frac{(w_{i}h_{i}^{\prime}(x))h^{T}(x)}{\|h_{i}(x)\|^{2}}.$
(14)

Taking the norm on both sides, squaring, and rearranging the terms yields
equation:weight_derivative.

Following the same strategy, we now prove a corresponding lemma for the biases
at each layer:

######  Lemma A.3.

In the notation above, for each layer $i$, we have:

$\|\nabla_{x}f_{\theta}(x)\|^{2}\left(\frac{1}{\|w_{i}\|\|h_{i}^{\prime}(x)\|}\right)^{2}\leq\|\nabla_{b_{i}}f(x)\|^{2}.$
(15)

######  Proof A.4.

Consider a small perturbation $x+\delta x$ of the input $x$. We start again by
showing that we can always find a corresponding perturbation $b_{i}+u(\delta
x)$ of the biases in layer $i$ such that

$f_{b_{i}}(x+\delta x)=f_{b_{i}+u(\delta x)}(x).$ (16)

Namely, because $f_{b_{i}}(x)=g_{i}(w_{i}h_{i}(x)+b_{i})$, to show this, it is
enough to find $u(\delta x)$ such that

$w_{i}(h_{i}(x)+h_{i}^{\prime}(x)\delta x)+b_{i}=w_{i}h_{i}(x)+b_{i}+u(\delta
x),$ (17)

where we again identify $h_{i}(x+\delta x)$ with its linear approximation
around $x$ for small $\delta x$. Then equation:condition2 is always satisfied
if we set this time

$u(\delta x):=w_{i}h_{i}^{\prime}(x)\delta x.$ (18)

Now taking the derivative with respect to $\delta x$ at $\delta x=0$ on both
sides of Equation equation:equivalence, and using the chain rule and that
$u(\delta x)$ is linear in $\delta x$, we obtain a relation between the
network derivative w.r.t. the biases and w.r.t the input:

$\nabla_{x}f_{\theta}(x)=\nabla_{w_{i}}f(x)w_{i}h_{i}^{\prime}(x).$ (19)

Taking the norm on both sides, squaring, and rearranging the terms yields
equation:bias_derivative.

Using this in quality for each layer in the decomposition given by Equation
equation:decomposition, we obtain that

$\|\nabla_{x}f_{\theta}(x)\|^{2}\left(\left(\frac{k_{1}(x)}{\|w_{1}\|}\right)^{2}+\cdots+\left(\frac{k_{l}(x)}{\|w_{l}\|}\right)^{2}\right)\leq\|\nabla_{\theta}f(x)\|^{2}$
(20)

Now if we use equation:weight_derivative and equation:bias_derivative in
equation:decomposition, we finally obtain that

$\|\nabla_{x}f_{\theta}(x)\|^{2}\left(\frac{1+\|h_{1}(x)\|^{2}}{\|w_{1}\|^{2}\|h_{1}^{\prime}(x)\|^{2}}+\cdots+\frac{1+\|h_{l}(x)\|^{2}}{\|w_{l}\|^{2}\|h_{l}^{\prime}(x)\|^{2}}\right)\leq\|\nabla_{\theta}f(x)\|^{2}.$
(21)


# Comments on the NEMA NU 4-2008 Standard on Performance Measurement of Small
Animal Positron Emission Tomographs

Patrick Hallen1, David Schug1,2, Volkmar Schulz1,2,3,4 1Department of Physics
of Molecular Imaging Systems, Institute for Experimental Molecular Imaging,
RWTH Aachen University, Aachen, Germany  2Hyperion Hybrid Imaging Systems
GmbH, Aachen, Germany  3III. Physikalisches Institut B, RWTH Aachen University
4Fraunhofer Institute for Digital Medicine MEVIS, Aachen, Germany
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The National Electrical Manufacturers Association’s (NEMA) NU 4-2008 standard
specifies methodology for evaluating the performance of small-animal PET
scanners. The standard’s goal is to enable comparison of different PET
scanners over a wide range of technologies and geometries used. In this work,
we discuss if the NEMA standard meets these goals and we point out potential
flaws and improvements to the standard.

For the evaluation of spatial resolution, the NEMA standard mandates the use
of filtered backprojection reconstruction. This reconstruction method can
introduce star-like artifacts for detectors with an anisotropic spatial
resolution, usually caused by parallax error. These artifacts can then cause a
strong dependence of the resulting spatial resolution on the size of the
projection window in image space, whose size is not fully specified in the
NEMA standard. If the PET ring has detectors which are perpendicular to a
Cartesian axis, then the resolution along this axis will typically improve
with larger projection windows.

We show that the standard’s equations for the estimation of the random rate
for PET systems with intrinsic radioactivity are circular and not satisfiable.
However, a modified version can still be used to determine an approximation of
the random rates under the assumption of negligible random rates for small
activities and a constant scatter fraction. We compare the resulting estimated
random rates to random rates obtained using a delayed coincidence window and
two methods based on the singles rates. While these methods give similar
esimates, the estimation method based on the NEMA equations overestimates the
random rates.

In the NEMA standard’s protocol for the evaluation of the sensitivity, the
standard specifies to axially step a point source through the scanner and to
take a different scan for each source position. Later, in the data analysis
section, the standard does not specify clearly how he different scans have to
be incorporated into the analysis, which can lead to unclear interpretations
of publicizeded results.

The standard’s definition of the recovery coefficients in the image quality
phantom includes the maximum activity in a region of interest, which causes a
positive correlation of noise and recovery coefficients. This leads to an
unintended trade-off between desired uniformity, which is negatively
correlated with variance (i.e. noise), and recovery.

With this work, we want to start a discussion on possible improvements in a
next version of the NEMA NU-4 standard.

##  1 Introduction

The National Electrical Manufacturers Association’s (NEMA) NU 4-2008 standard
on “Performance Measurements of Small Animal Positron Emission Tomography”
specifies “standardized methodology for evaluating the performance of positron
emission tomographs (PET) designed for animal imaging” [1]. The standard’s
goal is to enable comparison of the performance of different PET systems over
a wide range of technologies and geometries used. Thus, the methods specified
in the standard should not artificially favor or disfavor certain choices in
scanner geometry and technology and the performance results should indicate
the expected performance in real-world applications as closely as possible.
Virtually all commercial small-animal PET systems and most research prototype
PET systems have published performance evaluations based on the NEMA standard
and Goertzen et al. [2] have published a review comparing small-animal PET
systems based on the respective NEMA performance publications. These
publications are an essential benchmark in the development of new PET systems
and an important tool for the purchase decisions of potential buyers.

The NEMA standard specifies 5 measurements with respective analysis:
evaluation of spatial resolution; evaluation of total, true, scattered, random
and noise-equivalent count rates; evaluation of system sensitivity; and
quantitativ evaluation of image quality in a standardized imaging situation
using a hot-rod phantom.

The standard was devised over 10 years ago, so it does not incorporate newer
technological developments and paradigm shifts. For instance, the use of data
acquisition into sinograms and filtered backprojection reconstruction mandated
in the standard was more widespread than it is today. Nowadays, these methods
are often only implemented to evaluate the PET performance based on NEMA but
never actually used for real-world applications

In this work, we examine if the NEMA standard meets its goals to enable a fair
comparison of PET systems and we point out potential flaws and improvements in
the standard. In our opinion, the standard is underspecified in several parts,
limiting the comparability of different systems, since the investigators
performing the performance evaluations are still free to choose parameters
which significantly influence the results. The methods specified for
evaluation of the spatial resolution disadvantages certain system geometries,
where those geometries don’t exhibit the same reduction in spatial resolution
in real-world applications. The definition of random rates is circular and
allows the use of very different other methods generating different results.
The chapter on sensitivity is ambiguous, leading to publications using
different or even unclear methods for the measurement of sensitivity, creating
ambiguity in the interpretion of sensitivity of different PET systems.

If applicable, we demonstrate the claimed issues with simple simulation
studies. All discussions in this work should be universally applicable to any
PET system. However, it is still helpful and instructive to support the claims
in this work with real-world data. This is done using data obtained with the
Hyperion II${}^{\text{D}}$ PET/MRI scanner, which was developed by our group
[3]. Using the same data, we already have published a performance evaluation
based on the NEMA standard [4].

The goal of this work is to start a discussion on a revised version of the
NEMA standard and to provide input for this discussion.

##  2 Spatial Resolution

To evaluate the spatial resolution, the NEMA standard mandates the use of
point source scans which are reconstructed using filtered backprojection.
However, basically all modern PET scanners instead use an iterative maximum
likelihood expectation maximization (MLEM) algorithm for reconstruction [5, 6,
7, 8, 9, 10, 4, 11, 12, 13, 14, 15], so a scanner’s spatial resolution using
filtered backprojection is not necessarily indicative of its spatial
resolution for applications. While the mandated filtered backprojection is
intended to benchmark the detector performance alone, we will demonstrate in
the following that it disadvantages certain scanner geometries. Furthermore,
the NEMA standard specifies that the spatial resolution must be determined
using the projections of the reconstructed point sources inside a window in
image space, without strictly specifying the size of this projection window.
We will demonstrate that this can lead to an ambiguous spatial resolution
which depends on the size of the projection window and allows for artificially
enhancing the spatial resolution by choosing a particularly large projection
window for certain scanner geometries.

Figure 1: Ring geometry that was used for the simulations and the measurement.
The blue bands show the parallax error of LORs, which increases approximately
proportional to the angle $\varphi$ to the normal of the block detector.

The main disadvantage of filtered backprojection is that it doesn’t include
any model of the detector and assumes an ideal, ring-like PET scanner, while
the detectors in real-world PET scanners are usually in a block geometry with
anisotropic spatial resolutions. Line of responses (LORs) perpendicular to the
detectors front face are detected with the highest resolution, while tilted
LORs have a parallax error in the detected position, which increases with more
tilt of the LORs relative to the detector’s front faces as illustrated in
Figure 1. In principle, this effect can be reduced by detectors which are able
to determine the depth of interaction (DOI) of the gamma interaction, but in
practice most PET system don’t employ detectors with DOI determination [5, 6,
4, 11, 12, 14, 16]. Additionally, PET rings have gaps between the detector,
where no LORs are detected at all.

These issues with filtered backprojection will lead to artifacts in the
reconstructed activity. For instance, each angle where the PET ring has an
enhanced spatial resolution creates an excess in reconstructed activity along
the line connecting this position with the point source and each angle with
degraded spatial resolution creates a reduction in reconstructed activity
along the respective line. Similarly, gaps between the detector create a lack
of reconstructed activity along these lines.

To understand and demonstrate this behaviour, it is instructive to look at
these effects in sinogram space. In sinogram space, the enhanced spatial
resolution of perpendicular LORs manifests as hot-spots or rather peaks in the
center of each detector modules as Figure LABEL:fig:measurement_sinogram
shows. With increasing distance from the center of the detector module the
spatial resolution degrades, blurring the line of the point source in sinogram
space. This is equivalent to the convolution of the sinogram of a Gaussian
point source and the parallax error of the detector. The parallax error of the
detector stack can be modeled as the shape of two triangles, connected at
their tips as shown in Figure LABEL:fig:psf_sinogram. The parallelax error is
proportional to $\sin\varphi$, where $\varphi$ is the angle to the normal of
the block detector as defined in Figure 1. The parallax error shown in Figure
LABEL:fig:psf_sinogram is a small-angle approximation of this.

In addition to the inherent problems of mandating the use of filtered-
backprojection in the NEMA standard, the standard additionally mandates
projecting the reconstructed three-dimensional activity onto different one-
dimensional axes using a projection window. However, the size of the
projection window is not fully specified: ”The response function is formed by
summing all one-dimensional profiles that are parallel to the direction of
measurement and within at least two times the FWHM of the orthogonal
direction” [1, p. 7]. The first issue is that this definition is circular,
since the minimal size of the projection window to determine the FWHM is
defined using the FWHM itself. One can easily fix this problem, either using a
sufficiently large projection window in the first place, or by reducing the
size of the projection window iteratively in dependence of the determined FWHM
in the previous iteration. However, the much bigger problem is that the size
of the projection window can strongly influence the resulting spatial
resolution. The cause of this is the integration of the star-like artifacts
created by the anisotropic spatial resolution, as we demonstrate with the
following simulation, shown in Figure LABEL:fig:spatial_resolution.

We created the activity distribution of an ideally reconstructed point source
by assuming a rotationally symmetric two-dimensional normal distribution,
shown in Figure LABEL:fig:gaussian_point_source. The position of the point
source is off-center at a radial offset of 10 mm. To investigate the essence
of the effects, we don’t include noise in our simulation. From this ideally
reconstructed point source we create a sinogram by forward projection (i.e. by
applying a Radon transformation). The resulting sinogram is shown in Figure
LABEL:fig:ideal_sinogram.

We include the gaps between the detector stacks in our simulation by creating
a sensitivity sinogram, where all bins corresponding to gaps are 0 and bins
corresponding to sensitive detector area are 1 shown in Figure
LABEL:fig:gaps_sinogram. The simulated geometry is depicted in Figure 1 and
follows the geometry of the Hyperion II${}^{\text{D}}$ scanner to allow a
comparison between simulation and measurement. When we include this model of
gaps in our simulation by multiplying the sensitivity sinogram with our point-
source sinogram (Figure LABEL:fig:simulation_gaps) and then performing a
filtered backprojection (i.e. an inverse Radon transformation), we get a
reconstructed point source with slight artifacts, shown in Figure
LABEL:fig:simulation_gaps_FBP. As stated above, the artifacts are a lack of
reconstructed activity along the lines connecting the gaps and the point
source. When analyzing the spatial resolution of the filtered backprojection
with gaps we observe little influence of the gaps compared to the filtered
backprojection of an ideal sinogram without gaps. More importantly, the
resulting spatial resolution of 1.2 mm FWHM is stable to changes in the size
of the projection window, as shown in Figure
LABEL:fig:simulation_gaps_projection. Thus, gaps between the detectors are not
the cause of severe artifacts and only have a very minor influence on the
resulting spatial resolution with the usually small gaps of PET scanners.

When we additionally include the effect of the anisotropic detector resolution
due to parallax errors by convolving the point-source sinogram and the point
spread function in Figure LABEL:fig:psf_sinogram, the resulting filtered
backprojection in Figure LABEL:fig:simulation_gaps_doi_FBP exhibits a star-
like artifact, i.e. the lines connecting the center of each detector stack and
the point source exhibit a visible excess in activity.

Figure 2: Line profile of the reconstructed point source projected onto the
y-axis. The star-like artifact which is aligned with the x-axis creates an
excess in activity at the peak of the profile boosting the spatial resolution.

If one of these excesses aligns with one of the Cartesion projection axis, and
with the simulated geometry they do so for the x-axis, the projection onto the
axis perpendicular to this axis will result in a peaked excess at the maximum
of the line profile, as shown in Figure 2. A scanner’s spatial resolution is
defined by the FWHM and FWTM of this profile, which depends strongly on the
height of the maximum. Therefore, a peaked excess of the maximum will
significantly enhance the resulting spatial resolution. For our geometry, this
enhancement is only observed for the y-axis, because only the x-axis has an
excess in activity aligned with it, as there aren’t any detector stack which
are perpendicular to the y-axis. This difference between the resolution in x
and y is essentially an artifact and basically non-existant in real-world
applications using an iterative maximum-likelihood expectation maximization
(MLEM) reconstruction. Even worse, the extent of this effect depends strongly
on the size of the projection window as demonstrated in Figure
LABEL:fig:simulation_gaps_projection. Increasing the size of the projection
window enhances the resulting spatial resolution in y (i.e. decreases FWHM and
FWTM) while degrading the spatial resolution in x. This makes comparison of
the spatial resolution of different PET system difficult and maybe even
impossible, as the NEMA standard does neither specify a clear projection
window size nor does it mandata that the used window size should be reported.
Thus, most publications do not state the used projection window. [5, 16, 7,
14]. Other geometries may not exhibit this behavior at all, favoring or
disfavoring systems which have detectors perpendicular to a Cartesian axis.
One cannot even say that systems exhibiting this behavior have strictly worse
results, as such results can be sold in abstracts and conclusions (and
marketing brochures) as ”spatial resolution up to”, cherry picking the
artificially inflated spatial resolution along one of the axes.

The measurement and filtered backprojection reconstruction of point sources
with the Hyperion II${}^{\text{D}}$ scanner shown in Figure
LABEL:fig:measurement_sinogram and LABEL:fig:measurement_FBP look very similar
to the simulation which includes parallax error and gaps: The sinogram has the
same hot spots at the angles where the line of responses are perpendicular to
the detector surface and the reconstruction exhibits the same star-like
artifact. The analysis of the reconstruction yields the same observed
difference in spatial resolution between the x- and y-axis. Additionally, we
observe the same strong dependence on the size of the projection window, shown
in Figure LABEL:fig:measurement_projection.

An extreme example of a scanner geometry affected by this issue would be a box
geometry instead of the conventional ring geometry, i.e. a PET scanner with 4
large perpendicular detector modules without DOI capabilities. With such a
geometry, the filtered backprojection artifact would have the shape of a
cross, with both lines of excessive activity aligned with the x- and y-axis.
Thus, the artifact would enhance the resolution along both x- and y-axis by
boosting the maximum of both projections. This scenario is not solely
hypothetical, as small-animal PET scanners with the described box-like
geometry exist such as PETbox 4 [17]. In PETbox’s NEMA NU-4 performance
evaluation they state that using FBP was not possible ”since a FBP algorithm
specific for the PETbox4 system with the unconventional geometry has not been
developed” [17, p. 3797].

Other examples of published performance evaluation which have omitted the
filtered backprojection altogether when evaluating the spatial resolution are
[8, 18]. This is an indication that these groups don’t find the results based
on filtered backprojection not indicative for the performance of their system.

Fixing the issues of this method and proposing a better method to evaluate the
spatial resolution is challenging. The NEMA standards committee surely knew
many of these issues and we believe most of the PET community will be aware of
issues with filtered backprojection, as well. However, so far, none of the
performance publications based on NEMA discussed the issues presented here, so
we believe it is worthwhile to state them to start a discussion.

One obvious solution would be to simply not use filtered backprojection and to
perform the reconstruction with a modern iterative reconstruction instead, as
all real-world measurement would be performed with this iterative product
reconstruction anyway. However, modern iterative reconstructions include
features like resolution recovery, which, in theory, are able to reconstruct a
point source as a perfect point source in the limit of infinite statistics.
Thus, the reconstruction of a point source would mostly be a benchmark of the
reconstruction and not of the underlying detector performance. We suspect that
these arguments were the main reason why the NEMA standards committee chose
filtered backprojection instead.

One alternative could be the evaluation of spatial resolution using a Derenzo
hot-rod phantom. The standard could specify the geometry of such a phantom,
specify the activity and scan time, allow the use of the reconstruction method
supplied by the manufacturer and then define a quantitative analysis method.
The Derenzo phantom is already well-established in the community as a
benchmark to evaluate the spatial resolution. For instance, several NEMA
performance publications already include such a measurement as a benchmark of
spatial resolution [5, 7, 12, 15]. However, these results are not easily
comparable, as there currently isn’t a standardized quantitative analysis
method to determine the spatial resolution from a measurement of a Derenzo
phantom. Usually, the spatial resolution is estimated by making a qualitative
judgement at which distance the hot rods are still discernible. In principle,
such a definition of spatial resolution based on the ability to resolve to
close points is very reasonable and commonly used as a definition of spatial
or angular resolution for telescopes and microscopes [19, 20]. However, for a
quantitative definition of spatial resolution there must be a standardized
limit of the peak-to-valley ratio between two resolvable point sources, i.e.
how much the intensity between the two peaks must dip to make them just
resolvable. In optics, there are two commonly used criteria: The Rayleigh
criterion with an intensity dip of 26.5% and the Dawes criterion with an
intensity dip of 5% [21, p. 409]. In a standardized definition of PET spatial
resolution, the PET community could follow a similar criterion with either the
same intensity dip of 26.5% for consistency, or standardize a new arbitrary
limit.

(a) Reconstruction of Derenzo phantom scan. The red lines show an example of
profile lines which would be used for determination of valley-to-peak ratios
to evaluate the spatial resolution. (b) Distribution of valley-to-peak ratios
for the region with a rod distance of 0.9 mm. All ratios are below 0.735,
which is marked with a red vertical line.

Figure 3: Evaluation of spatial resolution using a Derenzo phantom

For the scan of a Derenzo phantom, such a resolvability criterion would
require to determine the valley-to-peak ratios of the profile lines over the
different regions of the phantom. To include anisotropies in the spatial
resolution, the profile lines should be defined over multiple angles as
demonstrated in Figure 3(a). Figure 3(b) shows the resulting distribution of
valley-to-peak ratios for the phantom’s 0.9-mm region. We would recommend that
the spatial resolution is defined as the hot-rod distance in the region where
at least 90% of the peak-to-valley ratios are below 0.735, i.e the valley dips
are above 26.5% for consistency with the Rayleigh criterion. Alternatively,
one could define a limit based on the average peak-to-valley ratio of a region
or using a different percentile than the suggested 90%. As shown in Figure
3(b), the region with distances of 0.9 mm has 100% of the valley-to-peak
ratios below 0.735. For the 0.8 mm-region, over half of the valley-to-peak
ratios would be above 0.735 in our measurement. Thus, the resulting spatial
resolution would be 0.9 mm.

To prevent arbitrary selection of peaks and valleys in a noisy reconstruction,
the standard could specify a limit for the allowed deviation from the physical
hot-rod distances when selecting the position of peak and valleys in the
profiles of the Derenzo region.

To evaluate the influence of radial and axial offsets on the spatial
resolution, the standard could specify different radial distances at which the
Derenzo phantom should be placed. Similarly, the standard could also specify
additional measurements of the rotated phantom to investigate the isotropy of
the spatial resolution.

In our opinion, such a method would depend much less on the system’s geometry
and technology and would provide a much more realistic benchmark, closely
mirroring real-world use of the system. As one of the disadvantages, the
precision of this method would be limited by the differences in hot-rod
distances between the phantom’s region. However, with commonly used Derenzo
phantoms one would achieve a precision in the determination of the spatial
resolution of 0.1 mm, which is more than adequate to assess the scanner’s
viability for intended applications.

As another alternative, Lodge et al. [22] have recently proposed a novel
method for the measurement of clinical PET spatial resolution using a
homogeneous cylinder phantom at an oblique angel. Such methods should also be
taken into consideration for an updated version of the standard.

##  3 Scatter Fraction, Count Losses, and Random Coincidence Measurements

The definitions of the random rate, scatter rate and scatter fraction are
circular and not satisfiable, and thus ill-defined for systems employing
detector material containing intrinsic radioactivity, such as LYSO or LSO
scintillators, as most modern PET systems do.

To explain this issue, we give a brief summary of the NEMA standard for the
measurement of the scatter fraction, count losses and random coincidence rate
in the following. The measurement is specified as a scan of an FDG-filled line
source inside a scatter phantom consisting of polyethylene. The rows of the
measured sinogram are centered at their maxima and the sum of all rows is
calculated. In the resulting radial profile of the phantom scan, the NEMA
standard specifies a signal window of 7 mm around the maximum. All event
counts outside this signal window are regarded as either scatter or randoms.
It is assumed that the sum of scatter and random event counts is at the same
level inside the signal window as on the edges of the signal window. The sum
of random and scatter event counts is denoted as $C_{r+s}$, and the sum of all
event counts are denoted as the total event count $C_{\text{TOT}}$.

For systems without intrinsic radioactivity, the scatter fraction is supposed
to be determined by assuming that the contribution of the random rate to the
combined scatter and random counts $C_{r+s}$ is negligible for measurements at
a low activity. Then, the random rate is determined from the total event rate
$R_{\text{TOT}}$ and true event rate $R_{t}$.

For systems with intrinsic radioactivity the sum of random and scatter event
counts also includes the random event counts produced by the intrinsic
radioactivity and this contribution of the intrinsic random rate cannot be
neglected at low measured activities. The NEMA standard therefore states: ”For
systems employing detector material containing intrinsic radioactivity, the
scatter fraction shall be evaluated by first evaluating the scattered event
counting rate (see section 4.4.5 below).” [1, p. 13] However, section 4.4.5
gives the following formula for the scattered event counting rate $R_{s}$,
which already includes the random rate $R_{r}$ [1, p. 14]

$R_{s}=R_{\text{TOT}}-R_{t}-R_{r}-R_{\text{int}}$ (1)

The formula for the random rate is given above, in section 4.4.4 and it
includes the scatter fraction $SF$

$R^{\text{NEMA}}_{r}=R_{\text{TOT}}-\left(\frac{R_{t}}{1-SF}\right)$ (2)

The scatter fraction $SF$, which is defined in the mentioned section 4.4.5, in
turn includes the scattered count rate

$SF=\frac{R_{s}}{R_{t}+R_{s}}$ (3)

Therefore, the definition of scattered rate $R_{s}$, scatter fraction $SF$ and
random rate $R_{r}$ are circular when reading the NEMA standard verbatim.

Next, we show that the definition is not only circular but also not
satisfiable. We insert the definition of $SF$ (i.e equation 3) into the
definition of $R_{r}$ (i.e. equation 2:

$\displaystyle R^{\text{NEMA}}_{r}$
$\displaystyle=R_{\text{TOT}}-\left(\frac{R_{t}}{1-\frac{R_{s}}{R_{t}+R_{s}}}\right)$
$\displaystyle=R_{\text{TOT}}-\left(\frac{R_{t}}{\frac{R_{t}+R_{s}-R_{s}}{R_{t}+R_{s}}}\right)$
$\displaystyle=R_{\text{TOT}}-R_{t}-R_{s}$

This is inserted into the definition of $R_{s}$ (i.e. equation 1):

$\displaystyle R_{s}$
$\displaystyle=R_{\text{TOT}}-R_{t}-\left(R_{\text{TOT}}-R_{t}-R_{s}\right)-R_{\text{int}}$
$\displaystyle=R_{s}-R_{\text{int}}\qquad\mbox{\lightning\;for}\;R_{\text{int}}\neq
0$

This is a contradiction, because by definition it is true that
$R_{\text{int}}\neq 0$, since the standard specifies these definitions of
$R_{r}$ and $R_{s}$ for scanners with intrinsic radioactivity.

We can speculate on the intended meaning of the NEMA standard’s definitions.
One sensible and probably originally intended modification to the definitions
would be to neglect the influence of the random rate $R_{r}$ (i.e. assume
$R_{r}=0$) in equation 1 for measurements at low activities to determine
$R_{s}$ and $SF$. We can then assume that $SF$ is approximately constant with
increasing activity and use $SF$ determined at a low activity to calculate the
random rates $R_{r}$ and scatter rates $R_{s}$ at higher activities.

This definition looks as if one could re-evaluate equations 1, 3 and 1
iteratively, starting with the assumption that $R_{r}=0$ at low activities, to
iteratively obtain a more accurate estimates of $R_{r}$, $R_{s}$ and $SF$.
However, this iteration diverges which is another indicator that the NEMA
standard’s definitions are not satisfiable. Nevertheless, this approach allows
to determine estimates of $R_{r}$, $R_{s}$ and $SF$ if we stop after one
iteration, i.e. we use the $SF$ determined at a low activity to calculate
$R_{r}$ and $R_{r}$ and then re-evaluate $SF$ at higher activities.

The NEMA standard specifies the following lower activity threshold: ”For
scanner employing, radioactive scintillator material, measurements shall be
performed until the single event rate is equal to twice intrinsic single event
rate” [1, p. 11]. Our scanner has an intrinsic single event rate of 80 kcps
and we reach a single event rate of 160 kcps at 430 kBq. Thus, we use this
activity to determine the scatter rate $R_{s}$ using equation 1 while
neglecting the random rate. This scatter rate is then used with equation 3 to
determine the scatter fraction $SF$. This scatter fraction is assumed to be
constant with varying activity and we use this with equation 2 to determine
the random rates $R_{r}$ at different activities. With these random rates we
can evaluate equation 1 and 3 again to determine the scatter rates and
fractions at higher activities without neglecting the random rates.

Alternatively, the NEMA standard allows the usage of a randoms estimate
supplied directly by the scanner. Such estimates usually use one of two
techniques: one using a delayed coincidence window (DCW) and one based on the
singles rates [23]. The singles rate (SR) method infers the randoms rate
$R_{ij}$ between to detector element $ij$ from the single rates $S_{i}$ and
$S_{j}$ using the formula

$R^{\text{SR}}_{ij}=2\tau S_{i}S_{j}$ (4)

with the time coincidence window $\tau$. However, this method systematically
overestimates the random randoms rate [24, 25]. Oliver et al. [26] proposed an
improved method ”Singles Prompt” (SP) which includes corrections based on the
coincidence rate (or prompt rate) $P_{i}$ to account for the contribution of
true coincidences and pile-up events:

$R^{\text{SP}}_{ij}=\frac{2\tau
e^{-(\lambda+S)\tau}}{(1-2\lambda\tau)^{2}}(S_{i}-e^{(\lambda+S)\tau}P_{i})(S_{j}-e^{(\lambda+S)\tau}P_{j}),$
(5)

where $\lambda$ is the solution of the equation

$2\tau\lambda^{2}-\lambda+S-P\,e^{(\lambda+S)\tau}=0.$ (6)

We have implemented these methods with the Hyperion II${}^{\text{D}}$ scanner
and can compare them empirically with the modified method the NEMA standard
suggests. The NEMA standard specifies a cylindrical signal window of 8 mm
around the phantom (i.e. a total diameter of 41 mm) in sinogram space. We
applied an equivalent cylindrical signal window, i.e. we only determined the
random rate for the pairs of detector elements whose line of responses
intersect with the cylindrical signal window.

Figure 4: Comparison of different methods for the determination of random
event rates: NEMA means a method based on the NEMA standard using equation 1,
DCW uses a delayed coincidence window, SR is based on the singles rate using
equation 4, and SP incorporates additional corrections using equation 5.

Figure 4 shows the total randoms rates as a function of activity inside the
scanner for the four different methods: NEMA, DCW, SR, and SP.

As expected, the SR estimate is larger than the SP estimate:
$R^{\text{SR}}\geq R^{\text{SP}}$. The DCW estimate is similar to the SP
estimate, and the adapted NEMA estimate is similar to the SR estimate, which
which is known to be the less precise than DCW and SP [26].

Oliver et al. [26] showed that random estimates $R^{\text{DW}}$ using a
delayed coincidence window (DW) are larger or equal to the SP estimates:
$R^{\text{SR}}\geq R^{\text{DW}}\geq R^{\text{SP}}$. There are many
publications investigating the correctness of these methods, providing
evidence from theory, simulations and measurements. For the NEMA method, on
the other hand, there doesn’t exist a single publication investigating the
correctness of the method, to our knowledge. Additionally, if the NEMA
standard is taken verbatim, the described random estimation method is not
implementable and the described estimation method has to be adapted as
described here, instead. Alternatively, the standard also allows the use of
any randoms estimation method instead. Altogether, the current standard does
not provide objective comparability of measured random rates of different
scanners.

All of these points apply also to the scatter rate $R_{S}$ defined in equation
1 and the noise-equivalent count rate, as the definitions of these observables
depend on the randoms rate. Therefore, objective comparison of these three
performance observables between different scanners based on the NEMA standard
is problematic.

##  4 Sensitivity

We think the NEMA standard’s protocol for the evaluation of the sensitivity is
unclear. Section 5.3 of the NEMA standard specifies to axially step a point
source through the scanner. Further, section 5.3.4 implies that a different
scan for each source position should be acquired. In section 5.4, all of the
data analysis is specified for single sinogram slices $i$. For instance, the
sensitivity is defined as

$S_{i}=\frac{R_{i}-R_{B,i}}{A_{\text{cal}}}$ (7)

with the counting rate $R_{i}$ and the background rate $R_{B,i}$ of sinogram
slice $i$. However, the NEMA standard only ever references sinogram slices and
never different measurements. We have one measurement per source position and
each of these measurements has many sinogram slices. In other words, there are
many measurements for each axial sinogram slice. Whenever the NEMA standard
refers to sinogram slices, it remains unclear which measurement to consider.
One possible intention could be to calculate the sum of all measurements;
however, this is never explicitly stated. This would effectively create a
sensitivity measurement with a virtual line source of activity $n\cdot A$,
where $n$ is the number of measurements. Such a line source would be similar
to the source distribution specified in the sensitivity protocol in the
clinical NEMA NU 2-2012 standard. However, the sensitivity $S_{i}$ is defined
by the activity $A_{\text{cal}}$ in equation 7, not a virtual activity $n\,A$
of the combined measurements. Unfortunately, the NEMA standard does not define
$A_{\text{cal}}$ in this equation, the only definition of $A_{\text{cal}}$ is
in section 1.2 as ”activity at time $T_{\text{cal}}$”. In conclusion, if this
interpretation were the intention of the NEMA standard, multiple required
instructions would be missing.

Another possible interpretation could be to take the slice $i$ of the
measurement where the point source is located at the center position of the
slice. However, this interpretation is not consistent with the formulas given
for the total system sensitivity

$S_{\text{tot}}=\sum_{\text{all}\,i}S_{i},$ (8)

which lack a normalization for the total number of slices. With a
normalization with the total number of slices, this would effectively be an
additional axial signal window around the point source, However, the size of
this axial signal window would depend on the scanner’s slice thickness, giving
an unfair disadvantage to high-resolution scanners. For instance, with a slice
thickness of 1 mm, this axial signal window would cut into the point source.
Additionally, this interpretation would not be realistic in the context of
real-world applications, where the sensitivity is supposed to be an indicator
of how many true coincidences one can expect for a given activity inside the
scanner’s FoV.

In summary, the NEMA standard does not include any instructions on how to
analyze the data of the multiple measurements it instructs to take. It only
defines the sensitivity of sinogram slices without specifying the relationship
of the sinogram slices and measurements with different source positions.

One consistent alternative definition of sensitivity could simply sum all
sinogram slices and then divide the total coincidence counts by the
acquisition time and activity for each measurement (i.e. source position). The
sensitivity profile would consequently be defined as this total sensitivity as
a function of the source position. To calculate the mouse- and rat-equivalent
sensitivities, one would average this sensitivity profile inside the central 7
cm or 15 cm. Because the NEMA standard specifies a transversal signal window
with a width of 20 mm in sinogram space, it would be consistent to apply the
same signal window around the point source in axial direction.

The ambiguity of the NEMA standard leads to unclear and incomparable results
in performance publication based on NEMA. Most publications seem to more or
less ignore the NEMA standard and simply measure the sensitivity for different
source positions as the total sum of all counts in this measurement. However,
the exact methods and definitions used stay usually unclear, impeding an
objective comparison of different sensitivity

For instance, Prasad et al. [13] seem to follow the formulas given by NEMA
quite closely, without clearly specifying how the data of the different
measurements at different source measurement is used in the data analysis. The
reported sensitivity profile has data points above 1 cps/Bq, i.e. an
impossible sensitivity larger than 100% for the central slices. They claim a
total absolute sensitivity of 12.74%, which is implausibly large compared to
the expected geometric sensitivity of 12.9%, We calculated this ideal
geometric sensitivity using their scanner’s diameter, axial length and crystal
thicknesses with the simple geometric model explained in [4]. The usual ratio
between measured peak sensitivity and geometric sensitivity is between 0.3 and
0.5 [4].

Other publications seem to more or less ignore the NEMA standard and simply
measure the sensitivity for different source positions as the total sum of all
counts in this measurement. However, basically all publications leave their
exact method unclear, preventing an objective comparison of different
sensitivity results.

##  5 Image Quality, Accuracy of Attenuation, and Scatter Corrections

The NEMA standard defines several observables for quantitative analysis of the
image quality phantom. The uniformity is defined as the relative standard
deviation of all voxels in a large cylindrical volume of interest over the
uniform region in the image quality phantom. For determination of the recovery
coefficients, the image slices along the central 10 mm of the hot rods are
averaged. Then, the recovery coefficients are defined as the maximum values in
a circular region of interest around the hot rods with different diameters,
divided by the mean activity in the volume of interest over the uniform
region. The issue with this definition is that the recovery coefficients are
correlated with the uniformity: The maximum value of a randomly distributed
sample increases with variance, even if the mean value of the distribution is
constant. Thus, this definition of the recovery coefficients does not measure
the mean recovery in the hot rods, but measures a combination of recovery and
variance. With a high variance and a good recovery the recovery coefficients
can even reach values larger than 1.

Figure 5: Simulated recovery coefficient of the 5-mm rod as a function of the
uniformity. The ground truth for the recovery was 0.95. The simulated recovery
coefficients are always larger than the ground truth and increase with
increasing variance (i.e. larger uniformity values).

We can demonstrate this behavior in a simple Monte Carlo simulation, where we
assume that the reconstructed activity in a voxel follows a normal
distribution with the standard deviation given by the uniformity. The
simulated geometry is the NEMA image quality phantom. Figure 5 shows the
simulated recovery coefficients of the 5-mm rod as a function of the
uniformity. The ground truth for the recovery coefficient for the activity in
the rod was 0.95. The data analysis follows the NEMA standard, i.e. the
recovery coefficient is defined by the maximum activity in the region of
interest. The drawn errors are calculated from the errors on the mean of the
averaged pixels in the region of interest. The simulation demonstrates that
the recovery coefficient is always overestimated compared to the ground truth
and increases with increasing variance (i.e. larger uniformity values).

Thus, the NEMA standard’s definition of the recovery coefficients hampers an
easy comparison of different scanner’s recovery performance, because the
recovery and uniformity must be compared at the same time. In other words, the
same scanner can achieve different recovery performance at different
uniformity points. The user can influence the uniformity with parameters such
as the amount of filtering during reconstruction. Figure 6 shows how the
measured uniformity and recovery coefficients changes with different widths of
a Gaussian kernel used during reconstruction for a scan of the image quality
phantom. We used the maximum likelihood expectation maximization
reconstruction described in [27]. As predicted by the Monte Carlo simulation,
the recovery coefficients are correlated with the relative standard deviation
in the uniformity region: Both values decrease with large filter width, i.e.
reduced variance in the image. Of course, it is not unexpected that the
recovery decreases with stronger filtering during reconstruction. However, the
observed effect is on top of the expected decrease in recovery due to
filtering. Using the NEMA standard’s observables, improving the uniformity
performance will always lead to a loss in observed recovery, regardless
whether the actual true recovery degraded or not. When conducting a NEMA
performance evaluation one has to chose an arbitrary point on the uniformity
and recovery curve resulting in one of many possible results, which are
difficult to compare with the results of other scanners.

Figure 6: Measured uniformity and recovery coefficients as a function of
filter width used during reconstruction. A larger filter reduces variance and
therefore increases uniformity (i.e. decreasing relative standard deviation).
The recovery coefficients decrease at the same time, so overall image quality
performance is a trade-off between uniformity and recovery.

As another minor issue, the NEMA standard derives the standard deviation of
the recovery coefficients from the standard deviations of the line profiles
along axial directions and the standard deviations of the uniform regions
using Gaussian error propagation. This is not the correct standard deviation
of the recovery coefficient, because the standard deviation of the maximum
value of a randomly distributed value is not the standard deviation of the
underlying distribution.

Fixing the definition of the recovery coefficients is not trivial. The NEMA
standard probably uses the maximum due to the small diameters of the hot rods.
For the very small rods, very few, if any, voxels lie clearly in the center of
the rods. Alternative definitions using the mean in a volume of interest will
therefore be biased by the smaller reconstructed activity in the border
regions of the rods. However, with today’s high-resolution PET scanners, we
believe it would be possible for most scanners to define volume of interest
(VoI) inside the hot rods and then define the recovery coefficients using the
mean reconstructed activity inside the VoI. Even if these VoIs would partially
include the border regions of the rods, it would still at least be a
comparable measure of recovery for every scanner. For the larger rods it
should not be any problem to define VoI which are well inside the hot rods
with a sufficient number of voxels. It is these larger rods where the current
definition of recovery coefficients leads to basically a recovery of 1 or
larger for all current scanners, hindering a differentiation of subtle
differences in recovery between the scanners.

##  6 General points

The NEMA standard does not explicitly mandate the use of the same settings for
each measurement. Most scanners offer a multitude of settings for measurements
and data processing, such as trigger settings, coincidence and energy window
sizes and quality filters for gamma interactions. The choice of setting
parameters requires often a trade-off for different performance parameters.
For example, the sensitivity benefits from wide energy and coincidence windows
and no quality filters, while the image quality and spatial resolution
benefits from narrow windows and strict quality filters. One could report very
misleading performance results by optimizing the settings for each performance
measurement separately, thus achieving performance results which are
unattainable at the same time in real-world applications.

While following the standard, many performance publications based on NEMA do
neither state if they used the same settings for every measurement explicitly
nor report all used settings for each measurement. For example, Nagy et al.
[5] use wide energy windows for the sensitivity and count rates measurements
and a narrow energy window for the measurement of spatial resolution. They do
not report any settings for the image quality measurement.

Another issue is the mandated use of sinograms. The data analysis for every
measurement except the image quality measurement are described on sinograms.
However, most modern scanners store their data in listmode format and often
only implement sinogram support to conduct the NEMA measurements. To our
knowledge, all NEMA NU-8 measurements published in in the last 5 years had to
convert their native listmode files to sinograms after the measurements [5, 6,
7, 8, 9, 10, 4, 11, 12, 13, 14, 15]. Spinks et al. [8] even mentions that the
calculation of scatter fractions were omitted due to missing sinogram support,
so this performance evaluation did apparently only use listmode data for the
data analysis. The number of scintillator crystals or equivalent positioning
bins in monolithic scintillators is usually above 30 000 in modern small-
animal PET systems, so that full sinograms have a file size of multiple
gigabytes even for very very short measurements. Listmode files on the other
hands are usually much smaller, making sinograms much more unwieldy.

The data analyses in the NEMA standard could be specified without the use of
sinograms, since most of the specified cuts in the sinograms could be
specified as cylindrical cuts in the scanner’s field of view. The standard
could still allow the use of sinograms as one possibility to implement the
specified geometric cuts for backwards compatibility.

##  7 Conclusion

Eleven years after the publication of the NEMA NU-4 standard, we believe it’s
time for a revision of the standard. In this work, we have pointed out several
flaws in the standard which should be addressed in the next version.
Additionally, the new technological developments in the last decade would
warrant discussing an updated version in itself. With this publication, we
would like to open this discution.

## Bibliography

## References

  * [1] National Electrical Manufacturers Association (NEMA): NEMA NU4-2008: Performance Measurements of Small Animal Positron Emission Tomographs (2008)
  * [2] Goertzen, A.L., Bao, Q., Bergeron, M., Blankemeyer, E., Blinder, S., Cañadas, M., Chatziioannou, A.F., Dinelle, K., Elhami, E., Jans, H.-S., Lage, E., Lecomte, R., Sossi, V., Surti, S., Tai, Y.-C., Vaquero, J.J., Vicente, E., Williams, D.A., Laforest, R.: NEMA NU 4-2008 Comparison of Preclinical PET Imaging Systems. Journal of Nuclear Medicine 53(8), 1300–1309 (2012). doi:10.2967/jnumed.111.099382. Accessed 2015-07-29
  * [3] Weissler, B., Gebhardt, P., Dueppenbecker, P.M., Wehner, J., Schug, D., Lerche, C.W., Goldschmidt, B., Salomon, A., Verel, I., Heijman, E., Perkuhn, M., Heberling, D., Botnar, R.M., Kiessling, F., Schulz, V.: A Digital Preclinical PET/MRI Insert and Initial Results. IEEE transactions on medical imaging 34(11), 2258–2270 (2015). doi:10.1109/TMI.2015.2427993
  * [4] Hallen, P., Schug, D., Weissler, B., Gebhardt, P., Salomon, A., Kiessling, F., Schulz, V.: PET performance evaluation of the small-animal Hyperion IID PET/MRI insert based on the NEMA NU-4 standard. Biomedical Physics & Engineering Express 4(6), 065027 (2018). doi:10.1088/2057-1976/aae6c2. Accessed 2019-01-23
  * [5] Nagy, K., Tóth, M., Major, P., Patay, G., Egri, G., Häggkvist, J., Varrone, A., Farde, L., Halldin, C., Gulyás, B.: Performance Evaluation of the Small-Animal nanoScan PET/MRI System. Journal of Nuclear Medicine 54(10), 1825–1832 (2013). doi:10.2967/jnumed.112.119065. Accessed 2016-03-02
  * [6] Ko, G.B., Kim, K.Y., Yoon, H.S., Lee, M.S., Son, J.-W., Im, H.-J., Lee, J.S.: Evaluation of a silicon photomultiplier PET insert for simultaneous PET and MR imaging. Medical Physics 43(1), 72–83 (2015). doi:10.1118/1.4937784. Accessed 2018-06-07
  * [7] Omidvari, N., Cabello, J., Topping, G., Schneider, F.R., Paul, S., Markus Schwaiger, Ziegler, S.I.: PET performance evaluation of MADPET4: a small animal PET insert for a 7 T MRI scanner. Physics in Medicine & Biology 62(22), 8671 (2017). doi:10.1088/1361-6560/aa910d. Accessed 2018-06-07
  * [8] Spinks, T.J., Karia, D., Leach, M.O., Flux, G.: Quantitative PET and SPECT performance characteristics of the Albira Trimodal pre-clinical tomograph. Physics in Medicine and Biology 59(3), 715 (2014). doi:10.1088/0031-9155/59/3/715. Accessed 2017-02-01
  * [9] Miyake, K.K., Matsumoto, K., Inoue, M., Nakamoto, Y., Kanao, S., Oishi, T., Kawase, S., Kitamura, K., Yamakawa, Y., Akazawa, A., Kobayashi, T., Ohi, J., Togashi, K.: Performance Evaluation of a New Dedicated Breast PET Scanner Using NEMA NU4-2008 Standards. Journal of Nuclear Medicine 55(7), 1198–1203 (2014). doi:10.2967/jnumed.113.131565. Accessed 2017-02-01
  * [10] Sato, K., Shidahara, M., Watabe, H., Watanuki, S., Ishikawa, Y., Arakawa, Y., Nai, Y.H., Furumoto, S., Tashiro, M., Shoji, T., Yanai, K., Gonda, K.: Performance evaluation of the small-animal PET scanner ClairvivoPET using NEMA NU 4-2008 Standards. Physics in Medicine and Biology 61(2), 696–711 (2016). doi:10.1088/0031-9155/61/2/696
  * [11] Mackewn, J.E., Lerche, C.W., Weissler, B., Sunassee, K., de Rosales, R.T.M., Phinikaridou, A., Salomon, A., Ayres, R., Tsoumpas, C., Soultanidis, G.M., Gebhardt, P., Schaeffter, T., Marsden, P.K., Schulz, V.: PET Performance Evaluation of a Pre-Clinical SiPM-Based MR-Compatible PET Scanner. IEEE Transactions on Nuclear Science 62(3), 784–790 (2015). doi:10.1109/TNS.2015.2392560
  * [12] Wong, W.-H., Li, H., Baghaei, H., Zhang, Y., Ramirez, R.A., Liu, S., Wang, C., An, S.: Engineering and Performance (NEMA and Animal) of a Lower-Cost Higher-Resolution Animal PET/CT Scanner Using Photomultiplier-Quadrant-Sharing Detectors. Journal of Nuclear Medicine 53(11), 1786–1793 (2012). doi:10.2967/jnumed.112.103507. Accessed 2015-05-04
  * [13] Prasad, R., Ratib, O., Zaidi, H.: NEMA NU-04-based performance characteristics of the LabPET-8™ small animal PET scanner. Physics in Medicine and Biology 56(20), 6649 (2011). doi:10.1088/0031-9155/56/20/009. Accessed 2015-01-21
  * [14] Szanda, I., Mackewn, J., Patay, G., Major, P., Sunassee, K., Mullen, G.E., Nemeth, G., Haemisch, Y., Blower, P.J., Marsden, P.K.: National Electrical Manufacturers Association NU-4 Performance Evaluation of the PET Component of the NanoPET/CT Preclinical PET/CT Scanner. Journal of Nuclear Medicine 52(11), 1741–1747 (2011). doi:10.2967/jnumed.111.088260. Accessed 2015-05-05
  * [15] Krishnamoorthy, S., Blankemeyer, E., Mollet, P., Surti, S., Holen, R.V., Karp, J.S.: Performance evaluation of the MOLECUBES BETA-CUBE—a high spatial resolution and high sensitivity small animal PET scanner utilizing monolithic LYSO scintillation detectors. Physics in Medicine & Biology 63(15), 155013 (2018). doi:10.1088/1361-6560/aacec3. Accessed 2018-08-13
  * [16] Bao, Q., Newport, D., Chen, M., Stout, D.B., Chatziioannou, A.F.: Performance Evaluation of the Inveon Dedicated PET Preclinical Tomograph Based on the NEMA NU-4 Standards. Journal of Nuclear Medicine 50(3), 401–408 (2009). doi:10.2967/jnumed.108.056374. Accessed 2015-04-22
  * [17] Gu, Z., Taschereau, R., Vu, N.T., Wang, H., Prout, D.L., Silverman, R.W., Bai, B., Stout, D.B., Phelps, M.E., Chatziioannou, A.F.: NEMA NU-4 performance evaluation of PETbox4, a high sensitivity dedicated PET preclinical tomograph. Physics in Medicine and Biology 58(11), 3791 (2013). doi:10.1088/0031-9155/58/11/3791. Accessed 2017-01-27
  * [18] España, S., Marcinkowski, R., Keereman, V., Vandenberghe, S., Holen, R.V.: DigiPET: sub-millimeter spatial resolution small-animal PET imaging using thin monolithic scintillators. Physics in Medicine and Biology 59(13), 3405 (2014). doi:10.1088/0031-9155/59/13/3405. Accessed 2014-08-22
  * [19] Lord Rayleigh, F.R.S.: XXXI. Investigations in optics, with special reference to the spectroscope. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 8(49), 261–274 (1879). doi:10.1080/14786447908639684. Accessed 2018-08-16
  * [20] Born, M., Wolf, E.: Principles of Optics, 7th edn. Cambridge University Press (1999)
  * [21] McKechnie, T.S.: Telescope Resolution and Optical Tolerance Specifications. In: General Theory of Light Propagation and Imaging Through The Atmosphere. Springer Series in Optical Sciences, pp. 405–464. Springer, ??? (2016). 10.1007/978-3-319-18209-4_14
  * [22] Lodge, M.A., Leal, J.P., Rahmim, A., Sunderland, J.J., Frey, E.C.: Measuring PET Spatial Resolution Using a Cylinder Phantom Positioned at an Oblique Angle. Journal of Nuclear Medicine 59(11), 1768–1775 (2018). doi:10.2967/jnumed.118.209593. Accessed 2019-02-01
  * [23] Hoffman, E.J., Huang, S.C., Phelps, M.E., Kuhl, D.E.: Quantitation in positron emission computed tomography: 4. Effect of accidental coincidences. Journal of Computer Assisted Tomography 5(3), 391–400 (1981)
  * [24] Rafecas, M., Torres, I., Spanoudaki, V., McElroy, D.P., Ziegler, S.I.: Estimating accidental coincidences for pixelated PET detectors and singles list-mode acquisition. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 571(1), 285–288 (2007). doi:10.1016/j.nima.2006.10.084. Accessed 2018-03-08
  * [25] Oliver, J.F., Rafecas, M.: Improving the singles rate method for modeling accidental coincidences in high-resolution PET. Physics in Medicine & Biology 55(22), 6951 (2010). doi:10.1088/0031-9155/55/22/022. Accessed 2018-03-08
  * [26] Oliver, J.F., Rafecas, M.: Modelling Random Coincidences in Positron Emission Tomography by Using Singles and Prompts: A Comparison Study. PLoS One 11(9), 0162096 (2016). doi:10.1371/journal.pone.0162096. Accessed 2016-09-14
  * [27] Salomon, A., Goldschmidt, B., Botnar, R., Kiessling, F., Schulz, V.: A Self-Normalization Reconstruction Technique for PET Scans Using the Positron Emission Data. IEEE Transactions on Medical Imaging 31(12), 2234–2240 (2012). doi:10.1109/TMI.2012.2213827