Department of Physics and Astronomy
Heidelberg University

Master thesis
in Physics
submitted by
Sarah Schuhegger
born in Traunstein
2021

Body Part Regression for CT Images

This Master thesis has been carried out by Sarah Schuhegger
at the
German Cancer Research Center
under the supervision of
Prof. Dr. Klaus H. Maier-Hein,
Dr. Fabian Isensee,
Lisa Kausch,
and
Prof. Dr. Carsten Rother
Visual Learning Lab

Body Part Regression for CT Images:

One of the greatest challenges in the medical imaging domain is to
successfully transfer deep learning models into clinical practice. Since
models are often trained on a specific body region, a robust transfer into the
clinic necessitates the selection of images with body regions that fit the
algorithm to avoid false-positive predictions in unknown regions. Due to the
insufficient and inaccurate nature of manually-defined imaging meta-data,
automated body part recognition is a key ingredient towards the broad and
reliable adoption of medical deep learning models. While some approaches to
this task have been presented in the past, building and evaluating robust
algorithms for fine-grained body part recognition remains challenging. So far,
no easy-to-use method exists to determine the scanned body range of medical
Computed Tomography (CT) volumes. In this thesis, a self-supervised body part
regression model for CT volumes is developed and trained on a heterogeneous
collection of CT studies. Furthermore, it is demonstrated how the algorithm
can contribute to the robust and reliable transfer of medical models into the
clinic. Finally, easy application of the developed method is ensured by
integrating it into the medical platform toolkit Kaapana and providing it as a
python package on GitHub.

Körperbereichs-Regression für CT Bilder:

Eine der größten Herausforderungen im Bereich der medizinischen Bildgebung
besteht darin, Deep-Learning-Modelle erfolgreich in die klinische Praxis zu
übertragen. Da Modelle häufig auf einer bestimmten Körperregion trainiert
werden, erfordert eine robuste Übertragung in die Klinik die Auswahl von
Bildern mit Körperregionen, die zum Algorithmus passen, um falsch-positive
Vorhersagen in unbekannten Regionen zu vermeiden. Aufgrund der typischerweise
unzureichenden Informationen des aufgenommenen Körperbereichs in den
Bildmetadaten, stellt die automatische Erkennung von Körperbereichen ein
Schlüsselelement für die breite und zuverlässige Einführung von medizinischen
Deep-Learning-Modellen dar. Obwohl in den letzten Jahren einige Ansätze für
diese Aufgabe vorgestellt wurden, bleibt das Entwickeln und die Evaluation
robuster Algorithmen für die genaue Körperbereichserkennung eine
Herausforderung. Bisher gibt es keine einfach zu handhabende Methode zur
Bestimmung des gescannten Körperbereichs von medizinischen Computertomografie
(CT)-Volumina. In dieser Arbeit wird ein selbstüberwachtes Körperbereichs-
Regressionsmodell für CT-Volumen entwickelt und auf einer heterogenen Sammlung
von CT-Studien trainiert. Weiterhin wird gezeigt, wie der Algorithmus zur
robusten und zuverlässigen Übertragung von medizinischen Modellen in die
klinische Praxis beitragen kann. Schließlich wird die Methode auf dem
medizinischen Plattform-Toolkit Kaapana und als einfach zu bedienendes Python-
Paket auf GitHub bereitgestellt.

###### Contents

  1. 1 Introduction
    1. 1.1 Motivation
    2. 1.2 Contributions
  2. 2 Background
    1. 2.1 Deep Learning
      1. 2.1.1 Learning Tasks
      2. 2.1.2 Learning Methodologies
      3. 2.1.3 Loss Functions
      4. 2.1.4 Regularization
    2. 2.2 Medical Background
      1. 2.2.1 Anatomical Foundations
      2. 2.2.2 Computed Tomography
    3. 2.3 Deployment
  3. 3 Related Work
    1. 3.1 Body Part Regression
      1. 3.1.1 Self-Supervised Body Part Regression
      2. 3.1.2 Blind Unsupervised-Supervision Network
    2. 3.2 Evaluation Metrics
    3. 3.3 Application and Deployment
  4. 4 Materials and Methods
    1. 4.1 Body Part Regression
      1. 4.1.1 Dataset Selection
      2. 4.1.2 Annotations
      3. 4.1.3 Learning Procedure
      4. 4.1.4 Order Loss
      5. 4.1.5 Data Augmentations
    2. 4.2 Evaluation Metrics
      1. 4.2.1 Landmark Mean Square Error
      2. 4.2.2 Accuracy
    3. 4.3 Application and Deployment
      1. 4.3.1 Estimate Examined Body Part
      2. 4.3.2 Known Region Cropping
      3. 4.3.3 Data Sanity Checks
      4. 4.3.4 Deployment
  5. 5 Experiments and Results
    1. 5.1 Body Part Regression
      1. 5.1.1 Loss Functions
      2. 5.1.2 Data Augmentations
      3. 5.1.3 Model Evaluation
      4. 5.1.4 State-of-the-Art Comparison
    2. 5.2 Application and Deployment
      1. 5.2.1 Estimate Examined Body Part
      2. 5.2.2 Known Region Cropping
      3. 5.2.3 Data Sanity Checks
      4. 5.2.4 Deployment
  6. 6 Discussion
    1. 6.1 Body Part Regression
    2. 6.2 Evaluation Metrics
    3. 6.3 Application and Deployment
  7. 7 Conclusion
    1. 7.1 Summary
    2. 7.2 Future Work
  8. 8 Appendix
  9. Acronyms
  10. Acknowledgement

##  Chapter 1 Introduction

###  1.1 Motivation

One of the greatest challenges in the medical imaging domain is to
successfully transfer deep learning models from homogeneous development
environments into clinical practice. Most algorithms are trained on a specific
body region, disease, or modality. If the models are then required to process
images that are out of their scope, one can no longer rely on their results to
be accurate. This can lead to potentially harmful and unintended consequences
for the clinical environment.

One important strategy to avoid these issues is to only present images to the
models close to their respective training distribution and reject images or
image regions that could be considered out-of-distribution. For example, in
the medical domain, the modality and scanned body part in the image meta-data
could be leveraged to decide whether an image fits a given algorithm.

The standard format for storing radiologic images and their meta-data is the
Digital Imaging and Communication in Medicine (DICOM) format. The meta-data is
stored in the DICOM header and includes, e.g., the patient’s name, examination
date, examined device settings, or modality. Furthermore, it includes fields
with the information of the scanned body part called Body Part Examined and
Anatomic Structure. While image meta-data can, to some extent, address the
above-outlined problem, it is inherently unreliable in practice
[gueld2002quality] and cannot solve cases where the desired anatomical region
is located somewhere within a larger image, such as a whole-body CT scan.

Automatic body part recognition can pave the way towards a broader application
of models in hospitals and can become an essential part of medical deep
learning software. Moreover, a better estimate of the examined body part than
provided by the meta-data can be obtained through automatic body part
recognition. This can lead to a significant time reduction in preprocessing,
sorting, and filtering large medical datasets.

So far, no easy-to-use method exists for automatic body part recognition in
Computed Tomography (CT) images. In addition, no easy applicable and fully
automated technique exists to crop CT images to the scope of medical deep
learning algorithms to provide robust end-to-end computer vision tools for the
clinic.

Previous publications delivered the proof of concept of self-supervised body
part regression for automatic body part recognition [yan2018deep,
yan2018unsupervised, tang2021body]. This thesis will investigate the
generalizability of the approach considering a heterogeneous collection of CT
studies and deploy the method into an easily applicable framework. Finally, it
will be shown how the method can be used as an essential pillar for robust
clinical end-to-end computer vision pipelines.

###  1.2 Contributions

Specifically, the work done in this thesis will focus on three main
contributions:

  1. C1:  Obtaining insights into the underlying optimization problem of robust body part regression to train a better-performing and generalizable body part regression model for CT images.

  2. C2:  Proposing a more fine-grained, thorough and robust evaluation strategy for cross-model performance comparison for body part regression models.

  3. C3:  Proposing three use cases of body part regression and simplify their application for the research community and the clinical environment. The presented use cases tackle the challenges of estimating the examined body part, cropping the known region from a CT volume for an algorithm, and carrying out basic data sanity checks.

##  Chapter 2 Background

In this thesis, a deep learning based body part regression model for CT
volumes is trained and deployed on GitHub and Kaapana. The foundations of deep
learning will be explained in section 2.1, and the medical background
regarding the human anatomy and Computed Tomography will be explained in
section 2.2. The software toolkit Kaapana which tries to build the bridge
between the deep learning research community and the clinical environment, is
explained in section 2.3.

###  2.1 Deep Learning

Deep learning is a subfield of machine learning that uses multi-layer neural
networks to solve optimization tasks numerically. Neural networks are
concatenations of linear and non-linear functions. They are a cascading of
many simple functions, which in the end learns complex representations of the
input [DeepVision, bishop2006pattern]. For example, in figure 2.1 a simplified
version of a standard fully connected feedforward neural network with one
input layer, three hidden layers, and one output layer can be seen. Each
hidden layer consists of several neurons, and each neuron consists of a series
connection of a linear function of the input values and a subsequent non-
linear function.

Figure 2.1: Simplified version of a fully connected feedforward neural
network. The Image was taken from [neuralnetwork].

Universal Function Approximation Theorem:  In general, what makes neural
networks so powerful, is that they are trainable parametric universal function
approximators. It means that they can learn any function with arbitrary
precision as long as the complexity of the neural network is high enough. This
holds because of the Universal Function Approximation Theorem, which states
that every feedforward neural network with at least one single hidden layer is
sufficient to represent any continuous function with arbitrary precision
[DeepVision]. Feedforward neural networks are networks where connections
between nodes do not form a cycle like, e.g., the network in figure 2.1. If
the network has just one hidden layer, it may be extremely big and fail to
learn. Therefore, from a practical point of view, it is better to use deeper
models. This leads to more complex representations of the input, which helps
the model to generalize better [krizhevsky2012imagenet].

Learning Framework:  From the Universal Function Approximation Theorem, we
know that a deep neural network $f_{\boldsymbol{\theta}}$ is, in theory, able
to represent any function with arbitrary precision, given a sufficient number
of linear and non-linear components. The main task that remains is to find
parameters $\boldsymbol{\theta}$ for the network to solve an optimization
problem reasonably good. The optimization problem is defined by an objective
function $\Phi$

$\Phi(\boldsymbol{\theta})=\frac{1}{N}\sum_{i=1}^{N}l(f_{\boldsymbol{\theta}}(\mathbf{X}_{i}),y_{i}),\qquad
f_{\boldsymbol{\theta}}(\mathbf{X}_{i})=\hat{y}_{i},$

where $N$ is the amount of available data pairs $(\mathbf{X}_{i},y_{i})$ and
$l$ the loss function per sample. For supervised tasks, the objective function
normally compares the predicted output from the neural network $\hat{y}_{i}$
and the ground truth label $y_{i}$. The loss is high if the prediction from
the neural network and the ground truth label differ and low if similar. The
objective function is numerically minimized by Gradient Descent and
Backpropagation. Gradient Descent is a numerical optimization method to find
extreme values of a function by following the negative gradient. The numerical
update step of classic gradient descent is given by

$\boldsymbol{\theta}=\boldsymbol{\theta}-\epsilon\nabla_{\boldsymbol{\theta}}\Phi(\boldsymbol{\theta}),$

where $\epsilon$ is the learning rate. The learning rate is a hyperparameter
that controls the step size of the parameter updates along the negative
gradient. If the learning rate is too low, the convergence may be slow. If the
learning rate is too high, the gradients might diverge. Backpropagation
describes how the gradient $\nabla_{\theta_{i}}L(\theta)$ for each parameter
$\theta_{i}$ can be found in a neural network through the chain rule
[DeepVision]. For training neural networks, a robust approximation of gradient
descents is used to save computation time. One common optimization strategy to
gain robust gradient estimates is called Adam [kingma2014adam]. The Adam
optimizer approximates the gradients by averaging over a mini-batch $B$ and
computing adaptive learning rates $\epsilon$ for stabilization
[kingma2014adam]. The adaptive learning rates are estimated from the first and
second moment of the gradients (for further details, refer to [DeepVision]). A
simplified version of the learning framework for a neural network
$f_{\boldsymbol{\theta}}(\mathbf{X}_{i})=\hat{y}_{i}$ with the mini-batch size
of $B=1$ can be seen in figure 2.2. For the optimization problem, it is
important to have many independent example data pairs $(\mathbf{X}_{i},y_{i})$
for training so that the neural network can find suitable parameters and
generalize on unseen new data. In general, before training neural networks the
total dataset of input-output pairs $(\mathbf{X}_{i},y_{i})$ is split into
three independent datasets: the training dataset, validation dataset and test
dataset. The training dataset is used to train the neural network. The
validation dataset is used to choose the best hyperparameters for the neural
network based on the network’s performance on the validation set. Finally, the
test set is the dataset with which the neural network’s generalization to new
data is estimated. The dataset is separated at the beginning to avoid
overestimating the generalizability of the model.

Figure 2.2: Neural network learning framework for mini-batch size of $B=1$ and
neural network $f_{\boldsymbol{\boldsymbol{\theta}}}$.

Convolutional Neural Networks: Different types of feedforward neural networks
$f_{\boldsymbol{\theta}}$ are distinguished based on their architecture and
input-output types. The state-of-the-art neural network architectures for
analyzing images are called Convolutional Neural Networks (CNNs). This thesis
concentrates on analyzing images. Therefore, we will focus on discussing the
architecture of CNNs in more detail. In the later part of the work, we will
come back to the VGG16 network, which is one of the mose common CNN
architectures introduced by Simonyan and Zisserman [simonyan2014very]. The
network architecture can be found in figure 2.3. Besides this, a standard
architecture of a CNN is shown in figure 2.4. In general, an image is feed
into multiple CNN building blocks, which are stacked together. One building
block consists of a convolutional layer, an activation layer, which is
normally the Rectified Linear Unit (ReLU) activation function, and a pooling
layer [DeepVision, bishop2006pattern]. The convolutional layer performs a
discrete convolution on the input image. The ReLU function is a simple non-
linear function given by

$\text{ReLU}(x)=\text{max(}0,x).$ (2.1)

Additionally, the pooling layer is a downsampling layer that reduces the size
of the input matrix. The output is flattened and used as input for a fully
connected neural network at the end of these building blocks. In summary, the
convolutional body transforms the input image into a dense feature map
encoding spatially separated semantic information. Finally, the fully
connected model head combines the semantic information from various regions to
produce a final prediction. For further details regarding the CNN architecture
and functionality, refer to [DeepVision].

Figure 2.3: Visualization of the VGG16 architecture. The image was taken from
[ferguson2017automatic]. Figure 2.4: Example of a CNN for predicting the
vehicle visible in an image. At the beginning, convolutional layers learn a
good feature representation for separating the classes. At the end fully
connected layers take over the classification similar to logistic regression.
The image was taken from [cnn].

####  2.1.1 Learning Tasks

We can divide deep learning tasks by the format of their target. This thesis
will concentrate on three main deep learning tasks: classification,
regression, and semantic segmentation.

Classification:  In classification tasks, categorical information is
predicted. For example, a neural network that predicts if a cat or a dog is
seen in an image would be a neural network for classification. Neural networks
for classification have a softmax function following the last layer, which
maps continuous values to values between zero and one. For a multi-class
classification problem, the softmax function $f$ is given by

$\begin{split}&f:\quad\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}:\quad\mathbf{z}\mapsto\frac{\text{exp}(\mathbf{z})}{\sum_{j}\text{exp}(z_{j})}=\mathbf{\hat{p}},\\\
&\mathbf{z}=\mathbf{W}^{T}\mathbf{x}+\mathbf{b};\quad
z_{i}=(\mathbf{W}^{T}\mathbf{x}+\mathbf{b})_{i},\end{split}$

with

$\sum_{i=1}^{m}\hat{p}_{i}=1,$

where $m$ is the number of distinct classes and $\mathbf{W}$ the weights of
the last fully connected layer and $\mathbf{b}$ the bias. With the outputs
$\hat{p}_{i}\in(0,1)$ of the last fully connected layer, the loss is computed.
For classification tasks, the standard loss function is the Cross-Entropy loss
(see sec. 2.1.3).

Regression:  If the target variable is continuous, we have a regression
problem to solve. The network architecture of classification networks and
regression networks look quite similar. The regression network has a linear
layer at the end of the network instead of a softmax function. Moreover, the
standard loss function for regression models is the  Mean Square Error (MSE)
(see sec. 2.1.3).

Semantic Segmentation:  Classification on pixel level is called semantic
segmentation. Semantic segmentation is mostly related to organ and disease
(e.g., tumor) segmentation in the medical domain. Example images of organ
segmentation can be found in figure 2.5. The standard baseline model for
medical image segmentation is the  no new U-Net (nnU-Net). The nnU-Net was
developed by Isensee et al. [isensee2021nnu]. It is an out-of-the-box
segmentation framework for organ and tumor segmentation for radiological
images. High performing semantic segmentation algorithms can be easily trained
without manual interventions for a variety of different dataset settings
[isensee2021nnu].

(a) Liver segmentation

(b) Spleen segmentation

Figure 2.5: Example segmentations of liver and spleen on the SegTHOR dataset
[lambert2020segthor]. The visualization was done with MITK [wolf2005medical].

####  2.1.2 Learning Methodologies

There are many deep learning methodologies, such as supervised learning,
unsupervised learning, semi-supervised learning, reinforcement learning, and
self-supervised learning. In this thesis, the two most relevant learning
methodologies are supervised learning and self-supervised learning. Therefore,
these two methodologies will be explained in the following.

Supervised Learning:  In supervised learning, the neural network maps an input
to an output based on many annotated input-output data points. Examples of
supervised learning are:

  * •

Classifying color images with the ImageNet database. The ImageNet database has
many different targets as dog, sailboat, human, and jet [deng2009imagenet].

  * •

Organ segmentation, where the network is trained on lots of pairs of three-
dimensional volumes and manually annotated segmentation masks.

  * •

Predicting the brain age through a regression network and using a database
with pairs of MRI images and the brain age [jonsson2019brain].

Self-Supervised Learning:  In self-supervised learning, the basic idea is that
the neural network architecture generates itself a kind of supervisory signal
through a proxy objective function to solve a task. Self-supervised learning
is often used to find a good initialization for supervised learning tasks,
where only a small annotated dataset is available. Through self-supervised
learning, meaningful representations of the data should be learned without the
need for annotated data. The self-supervised pretraining task is referred to
as "pretext task", and the actual task to be solved is referred to as
"downstream task". Examples for self-supervised learning tasks are:

  * •

Transform gray scale images to color images [zhang2016colorful].

  * •

Ordering image patches correctly. Cut puzzle pieces out of the image, shuffle
them and predict the correct permutation [noroozi2016unsupervised].

  * •

Ordering movie frames correctly. Shuffle movie frames and predict the
chronological order [lee2017unsupervised].

####  2.1.3 Loss Functions

For training a neural network, the loss function is crucial. The main loss
functions for classification (Cross-Entropy loss) and for regression (Mean
Square Error) will be explained in this section.

Cross-Entropy Loss:  The Cross-Entropy loss is the main loss function for
classification tasks. For $M$ classes and data point $i$ it is given by

$l_{i}=\sum_{c=1}^{M}y_{ic}\cdot\text{ln}(\hat{p}_{ic}),$

where $y_{ic}\in[0,1]$ is the ground truth of data point $i$ and
$\hat{p}_{ic}$ is the predicted probability for $p(y_{ic}=1)$. For $N$ data
points, the corresponding Cross-Entropy loss function is respectively given by

$L=\sum_{i=1}^{N}l_{i}(\mathbf{y_{i}},\mathbf{\hat{p}_{i}}).$

To understand the idea and the derivation behind the Cross-Entropy loss
better, we will derive the Cross-Entropy loss for the example of a binary
image classification problem with the ground truth label $y_{i}\in[0,1]$ (e.g.
cat or dog). We assume that the class $y_{i}$ is Bernoulli distributed

$\begin{split}&y_{i}\sim\mathcal{B}(p_{i}),\\\ &\Rightarrow
P(y_{i}=1|\mathbf{X}_{i})=p_{i},\quad P(y_{i}=0|\mathbf{X}_{i})=1-p_{i},\\\
&\Rightarrow
P(y_{i}|\mathbf{X}_{i})=p_{i}^{y_{i}}\cdot(1-p_{i})^{1-y_{i}}.\end{split}$

The probability for the positive class $p$ will be estimated by a neural
network $f_{\boldsymbol{\theta}}$

$f_{\boldsymbol{\theta}}(\mathbf{X}_{i})=\hat{p}_{i},\quad\hat{p}_{i}\in(0,1),$

where $\mathbf{X}_{i}$ is the image with the ground truth label $y_{i}$. For
this problem the likelihood is given by

$\begin{split}\mathcal{L}(\\{\mathbf{X},y\\}_{i\in[0,...,N]}|\
\\{\hat{p}_{i}\\}_{i\in[0,...,N]})&=\prod_{i=1}^{N}P(y=y_{i}|\hat{p}_{i}),\\\
&=\prod_{i=1}^{N}\hat{p}_{i}^{y_{i}}\cdot(1-\hat{p}_{i})^{1-y_{i}}.\end{split}$

The corresponding log-likelihood is given by the natural logarithm of the
likelihood function

$\text{ln
}\mathcal{L}=\sum_{i=1}^{N}{y_{i}}\cdot\text{ln(}\hat{p}_{i})+(1-y_{i})\cdot\text{ln(}1-\hat{p}_{i}).$
(2.2)

Maximizing the log-likelihood of equation 2.2 is equivalent of minimizing the
Cross-Entropy given by

$L=-\sum_{i=1}^{N}y_{i}\cdot\text{ln(}\hat{p}_{i})-(1-y_{i})\cdot\text{ln(}1-\hat{p}_{i}),$
(2.3)

where $N$ is the number of observed data points $(\mathbf{X},\ y)_{i}$.

Mean Square Error:  The Mean Square Error (MSE) is the main loss function for
regression tasks. For $N$ samples and the ground truth $y_{i}\in\mathbb{R}$
the MSE is given by

$L=\frac{1}{N}\sum_{i=1}^{N}||y_{i}-f_{\boldsymbol{\theta}}(\mathbf{X}_{i})||_{2}^{2}.$
(2.4)

Minimizing the Mean Square Error is equivalent to maximizing the likelihood of
the regression problem, if a Gaussian error $\boldsymbol{\epsilon}$ is assumed
[DeepVision]:

$\mathbf{y}=f_{\boldsymbol{\theta}}(\mathbf{X})+\boldsymbol{\epsilon},\quad\boldsymbol{\epsilon}\sim
N(0,\ \sigma^{2}\mathds{1}).$

####  2.1.4 Regularization

One of the main pitfalls in training a neural network is to overfit or to
underfit the problem. A model that overfits the training data performs
excellent on the training data itself but fails if it sees new data. A model
which underfits on the training data is performing poorly on the training data
and unseen data. If a model underfits a problem, the model complexity is
chosen too low to capture the problem’s complexity (see fig. 2.6). A higher
model complexity will resolve the underfitting problem. If a model overfits on
the training data, the model complexity is too high (see fig. 2.6). If the
complexity is too high, the model learns the noise in the training data as
well. This leads to a poor generalization on new data. The key method to
combat overfitting is regularization. With regularization techniques, the
complexity of a model is reduced, and the model performance on unseen data can
be improved. Regularization sets additional constraints to the optimization
problem or adds additional noise to the training process to reduce the
complexity of the learned model. There are different ways of regularization as
for example, using data augmentations, a regularization loss, or dropout. For
further details regarding regularization techniques, refer to [DeepVision].

Besides preventing overfitting, regularization is also used for ill-posed
problems, which are problems with no unique solution. Through regularization,
the solution space shrinks, and the numerical optimization is more stable
[buhlmann2011statistics].

Figure 2.6: Visualization of the overfitting vs. underfitting problem on a
simple regression model. The problems are based on an inappropriate model
complexity. The image was taken from [badillo2020introduction].

###  2.2 Medical Background

For this thesis, CT scans are used as the data basis. Moreover, the annotation
of robust anatomical landmarks is an important pillar of this work. Therefore,
the anatomical foundations will be discussed in section 2.2.1 and CT will be
explained in section 2.2.2.

####  2.2.1 Anatomical Foundations

For a better orientation in the human body, anatomical planes are used. They
are hypothetical planes to better describe locations and areas in the human
body. A distinction is made between the following planes (see fig. 2.7):

  * •

The axial plane is parallel to the ground and separates the head from the
feet.

  * •

The sagittal plane is perpendicular to the ground and separates the left arm
from the right one.

  * •

The coronal plane is perpendicular to the ground and separates the chest from
the back.

In figure A.3 the appearance of different body parts in a CT scan is
visualized. For each anatomical region, three different axial CT slices from
different patients can be seen. This demonstration emphasizes the inter-
patient variability in each anatomical region. The inter-patient variability
of organs is much higher than the inter-patient variability of bones. One
reason for this is that organs can move, and the position can depend, e.g., on
the breathing cycle [schiebler2013anatomie]. Therefore, for robust inter-
patient anatomical landmarks, bone landmarks are beneficial compared to organ
landmarks. The central support system of the human body which is suitable to
define bone landmarks is the spine.

Figure 2.7: Visualization of axial (gray), sagittal (blue) and coronal (green)
plane.

Spine:  The spine consists of five regions: Cervical, Thoracic, Lumbar,
Sacral, and Coccygeal [schiebler2013anatomie]. The Cervical spine consists of
seven vertebrae, the Thoracic spine of twelve vertebrae, and the Lumbar spine
of five vertebrae (see fig. 2.8). Each Thoracic vertebra carries one rib. The
Lumbar vertebrae need to carry the most weight. Therefore, the volume of the
Lumbar vertebrae is bigger than the vertebrae of the Thoracic spines or the
Cervical spine [schiebler2013anatomie]. In return, the vertebrae of the
Cervical spine are more flexible [schiebler2013anatomie]. One of the primary
diseases of the spine is scoliosis, where the spine has a lateral curvature
[kuznia2020adolescent]. This can lead to a bad posture and can affect
breathing, and movement [weinstein2008adolescent, yang2016early]. It occurs in
about 3 % of humans [shakil2014scoliosis].

In figure 2.8 b) the projection of the organs to the trunk can be seen. The
Thoracic vertebrae protect most of the upper body and upper abdomen organs,
including the lungs, the spleen, and the liver. Therefore, the vertebrae
positions can be used as a reference to describe organ positions in the upper
body.

Figure 2.8: The anatomy of the spine is visible in figure a) and the surface
projection of the organs of the trunk is visible in figure b). The images were
taken from [spine, anatomy-spine-projection].

####  2.2.2 Computed Tomography

Computed Tomography (CT) is a quantitative medical imaging method, which
produces images in Hounsfield unit (HU). It determines the X-ray absorption
coefficient in a matrix of equal size volume elements named voxels. Figure 2.9
shows a clinical CT scanner. It consists of a patient table and a gantry,
including the detectors and at least one X-ray tube. The detectors and the
X-ray tubes rotate around the patient during measurement time, and the patient
table moves horizontally. The X-ray tubes emit X-rays, and the detectors
measure the attenuation of the X-rays after passing through the patient.

Figure 2.9: Clinical CT scanners consists usually of an adjustable patient
table and a gantry in which detectors and X-ray tube rotate around the
patient. The image was taken from [ctscanner].

The Beer-Lambert law can describe the attenuation of the radiation

$I=I_{0}\cdot e^{-\int_{0}^{d}\mu(\mathbf{r})dr},$

where $I_{0}$ is the original beam intensity, $I$ the final beam intensity,
$d$ the distance between sender and receiver and $\mu(\boldsymbol{r})$ the
absorption coefficient at the position $\mathbf{r}$. The CT scanner measures
the radiation attenuation of a variety of different directions. From the
ensemble of all data points, the distribution of $\mu(\mathbf{r})$ can be
reconstructed. Finally, the distribution of $\mu(\mathbf{r})$ is linear
transformed into Hounsfield units and visualized in gray values. A typical CT-
volume consists of multiple slices with a layer thickness of 0.5 - 10 mm
[alkadhi2011funktioniert]. Each slice typically has the size of 512 x 512
pixels, with a corresponding pixel spacing of $0.2\text{mm/pixel}$ to
$1\text{mm/pixel}$ [alkadhi2011funktioniert]. These settings lead to a
measured cross-section of 10 cm to 50 cm.

Hounsfield Unit:  From the X-ray absorption coefficients $\mu$, the Hounsfield
unit can be calculated by

$\text{HU-value}=1000\cdot\frac{\mu-\mu_{w}}{\mu_{w}},$

where $\mu_{w}$ is the X-ray absorption coefficient of water. This is given by
$\mu_{w}\approx\text{0.192/cm}$ [alkadhi2011funktioniert]. It follows that
water always has the CT-value of 0 HU. Air has a CT-value of about - 1000 HU,
because its absorption coefficient lies around zero. In table 2.1 typical HU
values for different tissue types and organs are shown. In a typical CT-
volume, the data is saved in an unsigned 12 bit format with an offset of 1024
HU [schlegel2018medizinische]. Therefore, the HU values range between - 1024
HU till 3071 HU. In the first approximation, the attenuation coefficients are
proportional to the density of the tissue. Therefore, the CT values can be
interpreted in good approximation as density values. The quantitative physical
interpretation of the intensity units of a CT image opens up the opportunity
of setting intensity ranges for visualization, called windowing. Windowing can
be beneficial to enhance the visual contrast for particular organs or
pathology and filtering useless information.

Table 2.1: Organs and tissue types and their typical HU-values
[schlegel2018medizinische].

.  Tissue HU-value Tissue HU-value Air \- 1000 HU Pancreas 20 HU to 50 HU Lung
-900 HU to -500 HU Kidney 20 HU to 40 HU Fat -100 HU to -70 HU Liver 40 HU to
70 HU Water 0 HU Bone (cancellous) 70 HU to 350 HU Blood 30 HU to 60 HU Bone
(cortical) 350 HU to 2000 HU

Variability:  Although the intensity values of a CT scan have a physical
interpretation, variability across different CT scans from the same subject
exists. One reason for the variability in CT scans are artifacts. Artifacts
can be distinguished into physics-based, hardware-based, patient-based and
setup-related artifacts [boas2012ct, barrett2004artifacts]. They can cause
rings, noise, stripes and blurring in the CT scan [boas2012ct,
barrett2004artifacts]. Variability in contrast and intensity exists as well
across different scans of the same scanner based on different acquisition
parameters [berenguer2018radiomics]. Moreover, there is variability across
different CT scanners. Inter-scanner variability can be traced back to
different image reconstruction techniques, different scanner designs, and
initial scanner settings, but are not fully understood yet
[mackin2015measuring, berenguer2018radiomics, choe2019deep].

DICOM:  The standard format for storing and exchanging CT scans and other
radiological images is the Digital Imaging and Communication in Medicine
(DICOM) format. DICOM files contain additional meta information, such as
information about the scanner, the patient, and the imaged body part. In
hospitals, the DICOM files are saved in the Picture Archiving And
Communication Systems (PACS), the hospital’s internal medical database for
patient-related data.

###  2.3 Deployment

One main challenge for the medical computer vision research community is to
provide the developed algorithms in the clinical environment. The software
toolkit Kaapana aims to build the bridge between the deep learning research
community and the clinical environment. Kaapana is an open-source toolkit for
platform provisioning in the field of radiological data analysis [kaapana]. It
is not a ready-to-use software but rather a toolkit that can be used to build
customized medical platforms for data analysis and AI applications. Because
Kaapana builds upon Docker containers, users are enabled to add and remove
tools from the platform. Docker containers are isolated from one another,
including their own libraries, configuration files, and software. For further
details to Docker container refer to [docker]. Processing pipelines are build
in Airflow, where a workflow is called  Directed Acyclic Graph (DAG)
[airflow]. The bricks of the pipeline are called operators. Each operator
triggers a particular Docker container that runs a specific task. Kubernetes
is used to manage the Docker containers and Kibana to trigger the DAG
workflows (for further details, refer to [kubernetes, kibana]). To implement
new image analysis pipelines in Kaapana, a DAG needs to be built.

Figure 2.10 depicts the infrastructure of the Joint Imaging Platform (JIP)
developed by the German Cancer Research Center. The JIP is a medical platform
for data analysis that was built on Kaapana and is used in several German
university hospitals [scherer2020joint, JIP].

Figure 2.10: Infrastructure architecture of the Joint Imaging Platform (JIP),
which is build on Kaapana [scherer2020joint]. The image was taken from [JIP].

##  Chapter 3 Related Work

This work will refer to previously trained body part recognition models
explained in section 3.1 and their evaluation methods explained in section
3.2. The previously stated use cases for body part regression are presented in
section 3.3.

###  3.1 Body Part Regression

For automatic body part recognition, several approaches exist. For example,
through slice-wise body part classification, each slice in a body part gets
predicted to a particular region. The classification approach was studied
already in several studies [roth2015anatomy, yan2016multi, de20162d,
de2017convnet]. For example, Roth et al. [roth2015anatomy] implemented a CNN-
based body part classifier, which takes two-dimensional axial CT slices as
input and distinguishes between 5 classes (legs, pelvis, liver, lungs & neck).
They have shown that data augmentation is essential for training a robust body
part classifier and achieved an accuracy of 94.1 %. Most misclassifications
appeared between two neighboring classes because of ambiguous class labels.
Yan et al. [yan2016multi] extended the work of Roth et al. [roth2015anatomy]
and trained a CNN-based classifier with more classes to gain a fine-grained
prediction of the imaged body parts. In total, they defined 12 different
classes (nose, teeth, neck, shoulder, lung apex, sternal, aorta arch, cardiac,
liver upper, liver middle, abdomen/kidney, and ilium). A multi-stage deep
learning framework was proposed and applied for body part recognition. In the
first stage, the network discovers local discriminative regions, and in the
second stage, it learns a slice-level classification on the local regions.
With this approach, an $F_{1}-$score of about 92.23 % was reached. For the
slice-wise body part classification approach, an extensive dataset with
annotated images is needed. Furthermore, ambiguities will exist for slices
that lie in between two anatomic regions. This deliberation leads to the
assumption that anatomy-specific classification may not be the best choice to
convert human anatomy into a machine-readable form. The human body is a whole
continuous object instead of multiple disjoint classes. This consideration
strengthens the approach of a regression model to represent the natural order
in the human body.

There exist multiple approaches to create a body part regression model with
deep learning methods [yan2018unsupervised, yan2018deep, zhang2017self,
criminisi2010regression]. Most of these approaches use the intrinsic order
information in three-dimensional radiological volumes to train the models in a
self-supervised manner. Self-supervised methods have a crucial advantage in
their independence to manually annotated data. As the acquisition of annotated
data is usually time-consuming and expensive, self-supervised approaches shine
when applied to large datasets for which no manual annotations are available.
Therefore, in the following sections, the Self-Supervised Body Part Regression
(SSBR) model proposed by Yan et al. [yan2018unsupervised, yan2018deep] and the
Blind Unsupervised-Supervision Network (BUSN) model proposed by Tang et al.
[tang2021body] which are the most state-of-the-art self-supervised body part
regression models are explained. This thesis builds upon the SSBR model.

####  3.1.1 Self-Supervised Body Part Regression

The SSBR model was proposed by Yan et al.[yan2018unsupervised, yan2018deep].
The main goal of the SSBR model is to map every axial two-dimensional slice
$\mathbf{X}$ of a CT-image to a continuous score $s$

$f:\mathbb{R}^{n\times n}\rightarrow\mathbb{R},\quad\mathbf{X}\mapsto s,$
(3.1)

where $n$ denotes the width of the axial slices. The slice score $s_{i}$
should increase monotonously with slice index. Moreover, similar anatomical
regions should map to similar scores. The slice index itself can not be used
to gain body part related information, because CT images often have different
scan ranges and therefore different anatomical start and endpoints or inter-
slice intervals. Additionally, the size of people varies naturally.

Linear Relation: By the intrinsic information of the two-dimensional slice
ordering in a three-dimensional volume, the network should learn the mapping
from equation 3.1 by itself without any guidance through manual labels. Every
predicted slice score difference $\Delta s$ should be approximately
proportional to the spatial distance between the slices [yan2018unsupervised].
Formally, the desired target function corresponds to a linear mapping $g$
between slice index $j$ and slice score $s_{ij}$

$g_{i}:\mathbb{R}\rightarrow\mathbb{R},\quad j\mapsto s_{ij}\qquad\Rightarrow
g_{i}(j)=s_{ij},$

with:

$g_{i}(j)\approx m_{i}\cdot j+s_{i0},$ (3.2)

where the $s_{i0}$ indicates the start region of the CT-scan and the slope
$m_{i}$ depends on the height of the patient and the inter-slice distance of
the volume $\mathbf{X}_{i}$.

Loss Functions: One data item for training the SSBR model consists of $m$
equidistant two-dimensional axial CT slices, as it can be seen in figure 3.1.

Figure 3.1: One data item for training the SSBR model with the number of
slices $m=3$, the start index $j$ and the slice index difference $k$.

Random cropping is applied to every slice. After that, the score $s$ is
calculated by the network. With the help of a defined loss function, the
network obtains feedback and can learn through backpropagation. The loss
function consists of an order loss and a distance loss

$L=L_{\text{dist}}+L_{\text{order}},$

where the order loss $L_{\text{order}}$ is defined as

$L_{\text{order}}=-\sum_{i=1}^{B}\sum_{j=1}^{m-1}\text{ln}\ \sigma(\Delta
s_{ij}),\quad\Delta s_{ij}=s_{i,j+1}-s_{ij},$ (3.3)

where $\sigma$ is the sigmoid activation function, $B$ the batch size, $m$ the
sampled slices per volume and $s_{ij}$ dedicates the predicted slice score of
volume $i$ and slice index $j$ of the sampled data item. The sigmoid function
is given by

$\sigma(x)=\frac{1}{1+e^{-x}}.$ (3.4)

The order loss is big if the slice scores do not monotonously increase with
slice index. The distance loss is defined as

$L_{\text{dist}}=\sum_{i=1}^{B}\sum_{j=1}^{m-2}l_{1}(\Delta s_{i,j+1}-\Delta
s_{ij}),$ (3.5)

where $l_{1}$ is the smooth L1-loss [ren2016faster]. The smooth L1-loss is
given by

$l_{1}(\Delta s_{i,j+1}-\Delta s_{ij})=\begin{cases}0.5(\Delta
s_{i,j+1}-\Delta s_{ij})^{2}&|\Delta s_{i,j+1}-\Delta s_{ij}|<1\\\ |\Delta
s_{i,j+1}-\Delta s_{ij}|-0.5&\text{otherwise}.\end{cases}$ (3.6)

The distance loss is big if consecutive slice differences $\Delta s_{ij}$ are
particular different. The distance loss enforces the network to maintain a
linear relation between slice index and slice score, as seen in equation 3.2.

Figure 3.2: Proposed framework of the SSBR model with random two-dimensional
cropping as data augmentation step. The image was taken from [yan2018deep].

Model Architecture:  Figure 3.2 shows the structure of the SSBR model. For the
first convolutional layers, the Imagenet pre-trained VGG16 [simonyan2014very]
model was adopted. Then, an additional convolutional layer with 512 x 1 x 1
filters and a stride of one was added. After that, a ReLU layer, a global
average pooling layer, and a fully connected layer follow to predict the slice
score. The SSBR model was trained on 800 CT volumes for body ranges between
pelvis to chest [yan2018deep]. The model was evaluated on a 3-class
classification problem with the classes: pelvis, abdomen, and chest. The SSBR
model was able to achieve an accuracy of $95.99\%$ on the test set.

Research Gap:  The theory behind the loss functions of the SSBR model is not
fully understood until now. Moreover, finding good hyperparameters and
comparing SSBR models with different hyperparameter settings is challenging.
The accuracy metric is not sensitive to minor deviations in prediction and
makes it difficult to decide which model performs best. The scope of the SSBR
model is the pelvis to the chest body region. One of the remaining tasks is to
expand the algorithm’s scope to further body regions and to use more CT
studies for training to improve generalization towards different institutions
and scanner protocols.

####  3.1.2 Blind Unsupervised-Supervision Network

Recently, Tang et al. [tang2021body] proposed an additional modification to
the SSBR model. It is based on a two-stage approach, where in the first stage,
the SSBR model is trained. The predicted slice scores are corrected by a
robust linear regression algorithm (RANSAC). In the second stage, a 2D U-Net
model uses these corrected scores to train in a supervised fashion. An
additional Neighbor Message Passing (NMP) method was introduced to improve the
scores of a single slice by taking the information of the scores from the
neighbor slices into account [tang2021body]. The proposed modified model is
referred to as BUSN model. In figure 3.3 an overview of the BUSN model can be
found. The BUSN model was trained with 1030 CT volumes and was evaluated on
about 100 CT scans of the BTCV dataset [dataSynapse]. As evaluation metric,
the $R^{2}$-score was used (see sec. 3.2).

Research Gap: Due to the novelty of the BUSN model, it was not yet been used
and compared in other papers. The authors compared the proposed BUSN model
with the original SSBR model and reported relatively poor performance of the
baseline approach. From the paper [tang2021body], it is not clear where the
deviation in performance compared to the original SSBR paper
[yan2018unsupervised] comes from.

Figure 3.3: (A) The SSBR model proposed by Yan et al. [yan2018deep]. (B) Deep
supervised network using the refined prediction scores from the robust
regression. During test-time, only the second model part is required to
perform body part regression. The image was taken from [tang2021body].

###  3.2 Evaluation Metrics

Evaluating a body part regression model is non-trivial because no ground truth
labels are available. In the literature, different approaches for evaluating a
body part regression model exist.

Accuracy:  In the paper from Yan et al. [yan2018unsupervised] the primary
evaluation measure was the accuracy. Three classes were defined: pelvis,
abdomen, and chest. For testing and validation, each slice was labeled
manually. The abdomen class starts from the upper border of the ilium bone and
ends at the upper border of the liver. With the help of two extra validation
CT volumes, two thresholds to separate the three-body zones were identified.
Based on the predicted slice score, each slice is assigned to a class
$\hat{c}$. In the end, the accuracy $\psi$ compares the predicted classes
$\hat{c}$ with the annotated ground truth classes $c$ by calculating

$\begin{split}\psi=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{M_{i}}\sum_{j=1}^{M_{i}}\chi(\hat{c}_{ij},c_{ij}),\\\
\end{split}$ (3.7)

with

$\chi(\hat{c}_{ij},c_{ij})=\begin{cases}1&\text{if }\hat{c}_{ij}=c_{ij},\\\
0&\text{otherwise},\end{cases}$

where $N$ is the number of CT volumes in the validation set and $M_{i}$ the
number of slices in volume $i$. The continuous output slice score is converted
to discrete classes, implying that deviations of the scores within one class
are neglected. Thus, one limitation of the accuracy metric is that it is not
particularly sensitive.

$R^{2}$-Metric:  The $R^{2}$-metric is used by Tang et al. [tang2021body] to
validate their two-stage body part regression approach. For each slice score
curve, a linear curve $\hat{s}$ was fitted to the slice scores $s$. For a
volume $X$ with $M$ slices, the $R^{2}$-metric is given by

$R^{2}_{X}=\frac{\sum_{i=1}^{M}(s_{i}-\bar{s})^{2}}{\sum_{i=1}^{M}(\hat{s}_{i}-\bar{s})^{2}}=1-\frac{\sum_{i=1}^{M}(\hat{s}_{i}-s_{i})^{2}}{\sum_{i=1}^{M}(\hat{s}_{i}-\bar{s})^{2}},$
(3.8)

where $s_{i}$ is the predicted slice score on index $i$, $\hat{s}_{i}$ is the
linear fit at index $i$ and $\bar{s}$ the mean slice score from the linear
fit. On the one hand, an advantage of the $R^{2}$-metric is that no manually
annotated data is needed. On the other hand, the assumption of linear
correlation between body height and slice score does not consider inter-
patient anatomical variation. Therefore, one limitation of the $R^{2}$-metric
is that it does not measure if the same anatomical regions are mapped to the
identical scores across patients.

Research Gap:  Up to now, no fine-grained metric was proposed for regression
models that precisely measure if the same anatomical regions are predicted to
similar scores. Moreover, there exists no evaluation metric for body part
regression models on which the research community agreed on. Instead, each
research group uses a different evaluation metric, and no comparison of the
proposed evaluation metrics exists to agree on a universal metric for cross-
model comparison.

###  3.3 Application and Deployment

Due to the fine-grained recognition of body regions, body part regression
models have a broad range of possible use cases for medical image analysis
tasks. Some of them were already investigated in previous literature and will
be introduced in this section.

Annoation:  One evident use case is to obtain body part region annotations
through a body part regression model. For example, Yan et al. annotated the
position of lesions in the human body with a body part regression model to
structure a large lesion dataset [yan2018deep]. Further, diverse papers showed
that a body part regression model can predict rough estimates of organ
positions [tang2021body, proskurov2021fast], and that it is capable of
cropping certain body regions as the abdomen [tang2021high, tang2020learning,
xu2020outlier] or the pelvis [ivanovska2021efficient] from a CT scan.

Anomaly Detection:  A further use case is anomaly detection in radiologic
volumes by visualizing the predicted slice scores against the slice index
[yan2018unsupervised]. Yan et al. showed that for some volumes, the slice
score curves exhibit discontinuities where the appearance of the CT slices is
abnormal due to scanning artifacts or large lesions [yan2018unsupervised].

Preprocessing for Registration:  The body part regression model could also be
used as pre-alignment in z-direction for registration purposes. Heine and
Hering pointed out that preliminary results confirm that a fast and robust
pre-alignment of CT volumes through body part regression is possible
[meine2019efficient].

Preprocessing for Segmentation:  Tang et al. pointed out that it is possible
to enhance the quality of nnU-Net organ segmentation models in terms of dice
score by using body part regression [tang2021body]. Through introducing an
additional preprocessing step, the volumes from the training, validation, and
test set were spatially normalized with the body part regression model by
automatically cropping the abdomen region [tang2021body]. This preprocessing
step led to a significant increase in segmentation performances
[tang2021body].

Research Gap:  Although several use cases were already proposed, a broad field
of use cases remains not yet analyzed. For example, up to now, no one has
analyzed how the body part regression model can be used for false-positive
detection in medical algorithms. Medical algorithms are typically trained on
images from a specific body part and lead to false-positive predictions in
unknown regions. The body part regression model opens up new opportunities to
detect false-positive predictions in unknown regions, which can lead to a
significant increase in model performance of medical algorithms, such as,
e.g., tumor segmentation.

Yan et al. pointed out that finding outliers in CT volumes with the body part
regression model through visual inspection is possible [yan2018unsupervised].
Nevertheless, up to our knowledge, no automatic approach of finding outliers
within CT volumes or CT datasets due to, e.g., improper data loading or
incorrectly assigned modality exists.

Moreover, up to now, the possibility of gaining a better estimate for the
imaged body region than the BodyPartExamined tag in the DICOM meta-data
through the body part regression model was not analyzed in the previous
literature. Gueld et al. emphasized that the DICOM body part meta-data usually
is not sufficient to describe the examined body part correctly
[gueld2002quality]. Therefore, the body part regression model can bring
significant benefits in predicting robust and precise body part examined tags
for radiological images.

Regarding the deployment of a body part regression model, no user-friendly
model exists with which the body part regression model can be easily used in
the research community or the clinical environment.

##  Chapter 4 Materials and Methods

In this chapter, the theory behind the trained body part regression models is
explained (sec. 4.1) and the evaluation methods for cross-model comparison are
introduced (sec. 4.2). Building on this, in section 4.3 several use cases for
the model, are presented. Further, the model deployment is explained to make
the model publicly available for the clinical and research communities.

###  4.1 Body Part Regression

A self-supervised body part regression model is presented that automatically
maps similar anatomical regions to the same score across different patients.
For training a model with good generalization capabilities, a diverse dataset
was selected (sec. 4.1.1). For evaluating the model performance, robust
anatomical landmarks were defined and annotated in selected CT volumes (sec.
4.1.2). Additionally, the network architecture of the body part regression
model is described (sec. 4.1.3), different order loss terms for the model are
introduced (sec. 4.1.4) and the used data augmentation techniques are
explained (sec. 4.1.5).

####  4.1.1 Dataset Selection

In order to train a body part regression model with good generalization
properties, it is essential to use a heterogeneous dataset with CT volumes
from different hospitals, CT devices, and various diseases. The dataset used
for this thesis was compiled from studies of The Cancer Imaging Archive (TCIA)
[clark2013cancer] and different computer vision challenges (see tab. 4.1 and
tab. 4.2). Table 4.1 summarizes the 12 studies used for training and
validation, resulting in a dataset size of 2192 CT scans for training and 100
CT scans for validation. In table 4.2 an overview of the compiled test dataset
is shown, which consists of 10 different studies and in total 100 CT scans.
The pixel spacing in the z-direction varies between 0.45 mm and 6 mm. Only
images with a $z$-range greater than 12 cm and 42 or more slices were used for
training. The original pixel spacing in the axial plane in x- and y-direction
vary between 0.3 mm and 7.5 mm.

To avoid overfitting on studies and certain body regions, CT volumes were
limited to 150 CT images for the studies Task007 and LIDC-IDRI. For the
datasets, HNSCC, QIN-HEADNECK, ACRIN-NSCLC-FDG-PET, the limit was increased to
400 CT scans since those studies cover almost the whole body. By selecting the
400 volumes from the greater population, care was taken to favor volumes
showing rare body regions like the head or the thigh. Nearly all studies are
publicly available. Only the "whole body CT" study is private. This study
includes the knees and the feet as well. In nearly no other study, these body
areas are visible. Adding the "whole body CT" study to the dataset has two
advantages: on the one hand, the model sees more volumes in rare regions, and
on the other hand, it could be interesting to investigate the model behavior
on whole body CT volumes during validation. Nevertheless, the valid scope of
the trained body part regression model will be the lower pelvis until the end
of the head because there are too few volumes from the legs available to
extend the valid scope of the model. Moreover, due to a lack of data, body
regions of different modalities, children, or pregnant women are not in the
scope of the trained model.

Table 4.1: Overview of studies used for training and validation. For a subset of the training data and the whole validation set landmarks were annotated (see sec. 4.1.2). Some studies are quite heterogeneous in terms of the scanned body areas. For these datasets, the minimum and maximum body part of the scanned volumes is documented. Study name | Body target |  # Train. data |  # Annotated train. data |  # Val. data | Source
---|---|---|---|---|---
Task003 | liver | 177 | 4 | 9 | [simpson2019large, bilic2019liver]
Task007 | pancreas | 150 | 4 | 9 | [simpson2019large, attiyeh2018survival, attiyeh2019preoperative, chakraborty2018ct]
Task051 | upper body | 42 | 4 | 8 | [dataGrandChallange]
Task062 | pancreas | 74 | 4 | 8 | [clark2013cancer, dataPancreas, roth2015deeporgan]
whole body CT | foot - head | 45 | 4 | 6 | private
CT COLONOGRAPHY | colon | 74 | 4 | 8 | [clark2013cancer, johnson2008accuracy, dataCOLONOGRAPHY]
COVID-19-AR | lung | 113 | 5 | 9 | [clark2013cancer, jenjaroenpun2021two, desai2020chest, dataCOVID]
CT LYMPH NODES | pelvis - neck | 167 | 4 | 9 |  [clark2013cancer, dataLYMPH, roth2014new] [seff20142d, seff2015leveraging]
LIDC-IDRI | lung | 150 | 4 | 9 | [clark2013cancer, dataLIDC, armato2011lung]
HNSCC | pelvis - head | 400 | 5 | 9 | [clark2013cancer, dataHNSCC, grossberg2018imaging, elhalawani2017matched]
QIN-HEADNECK | pelvis - head | 400 | 4 | 8 | [clark2013cancer, dataQIN, fedorov2016dicom]
ACRIN-NSCLC- FDG-PET | pelvis - head | 400 | 4 | 8 | [clark2013cancer, dataACRIN, machtay2013prediction]
Sum |  | 2192 | 50 | 100 |
Table 4.2: Overview of test dataset. For each study, 10 CT volumes were randomly selected and the positions of the defined landmarks were annotated (see sec. 4.1.2). Study name | Body target | # Annotated data | Source
---|---|---|---
Task006 | lung | 10 | [simpson2019large, bakr2018radiogenomic, napel2014nsclc, gevaert2012non]
Task008 | abdomen | 10 | [simpson2019large, pak2018quantitative, simpson2015texture, simpson2017computed, zheng2017preoperative]
Task009 | spleen | 10 | [simpson2019large, simpson2015chemotherapy]
Task010 | colon | 10 | [simpson2019large]
Task017 | abdomen | 10 | [dataSynapse]
Task018 | pelvis | 10 | [dataSynapse]
Task046 | abdomen | 10 | [gibson2018multi]
Task049 | head-neck | 10 | [dataGrandChallange]
Task055 | thorax | 10 | [lambert2020segthor]
Task064 | kidneys | 10 | [heller2019kits19]
Sum |  | 100 |

Data Preprocessing:  The intensity values in Hounsfield unit were clipped
between -1000 HU and 1500 HU and rescaled to -1 to 1. Axial slices were
resampled to a pixel spacing of 3.5 mm. All images were resized to 128 px
$\times$ 128 px by zero-padding or center cropping. The pixel spacing was
chosen to assure that the resulting field-of-view is 44.8 cm $\times$ 44.8 cm
and covers the entire cross-section of an average person.

Before down-sampling, a Gaussian filter was used to avoid aliasing artifacts
[rao2018digital, engelberg2008digital]. For a downsampling factor of $a=0.25$,
the smoothing standard deviation $\sigma=0.8$ for the Gaussian kernel was
experimentally determined by visual inspection of resized images. This
Gaussian filter configuration reduces aliasing effects with small edge
information loss. Based on the original pixel spacing of the axial slices, we
gain different downsampling factors $a^{\prime}$. The most suitable smoothing
standard deviation $\sigma^{\prime}$ for the Gaussian kernel is dependent on
the downsampling factor. For bigger downsampling factors, we would assume
stronger smoothing and vise versa. Therefore we used the following relation to
gain the custom smoothing standard deviation $\sigma^{\prime}$ for a
downsampling factor of $a^{\prime}$:

$\sigma^{\prime}=\frac{a}{a^{\prime}}\sigma.$

####  4.1.2 Annotations

For evaluating the model performance, anatomical landmarks were defined. Table
4.3 summarizes all defined 35 landmarks and also contains a short description
of each landmark. All landmarks except the vertebrae landmarks were annotated
manually. For annotating the centroid of the vertebrae, a nnU-Net
[isensee2021nnu] vertebrae segmentation algorithm was used, which was kindly
provided by Fabian Isensee. The segmentation algorithm was trained on the data
from the Verse2020 segmentation challenge data [sekuboyina2020verse,
loffler2020vertebral]. With this segmentation algorithm, all vertebrae of the
Cervical spine, Thoracic spine, and Lumbar spine are segmented. In addition,
automatic vertebra landmark annotations were derived by computing the centroid
of the predicted vertebra masks. Segmenting the vertebrae of one volume took
about 20 minutes on an Nvidia V100 SXM2 GPU with 32 GB GPU memory.

Figure A.4 shows an example slice for each landmark. The positions of the
landmarks for 100 volumes of the validation set, 50 volumes of the training
set, and 100 volumes of the test set were annotated. It was taken care to
annotate the volumes uniformly across different studies. Therefore, in the
test dataset for every study 10 volumes were annotated, in the validation set
6 to 9 volumes and in the training set 4 to 5 volumes per study (see tab. 4.1
and tab. 4.2). For validation, only six volumes were annotated from the "whole
body CT" study to use more whole body CT volumes for training. A few misplaced
landmarks were identified in a second pass through the annotated slices. The
corresponding annotations were removed from the database. A summary of the
number of annotated landmarks in the training, validation, and test dataset
can be found in table A.1.

The anatomical variation of the different landmarks is highly diverse. In
general, bone landmarks are anatomically more robust than organ landmarks. For
the later evaluation of the models, only a subset of the defined bone
landmarks was used. The evaluation landmarks are pelvis-start, femur-end, L5,
L3, L1, Th11, Th8, Th5, Th2, C6, C1, and eyes-end. The evaluation landmarks
will be referred to as $l_{k}$ with $k\in\\{1,...,12\\}$. It was taken care of
the distance between evaluation landmarks to ensure that they have
approximately the same distance to each other. A detailed definition for each
landmark can be found in table 4.3. Figure 4.1 visualizes in two example
images where the evaluation landmarks are located in the human body.

Figure 4.1: Two example coronal CT slices from the validation set with highlighted evaluation landmark positions.  Table 4.3: Used anatomical landmarks and defined bone landmarks for evaluation. Landmark name | Description |  Eval. landmark | Notation
---|---|---|---
pelvis-start | Slice before ischium bone starts | ✓ | $l_{1}$
femur-end | Slice after the femur bone ends | ✓ | $l_{2}$
pelvis-end | Slice after the end of the illium bone |  |
kidneys | First slice where both kidneys can be seen well |  |
lung-start | Slice before lung starts |  |
liver-end | Slice after end of liver |  |
lung-end | Slice after end of lung |  |
teeth | Middle slice of lower jaw |  |
nose | Slice where tip of nose can be seen the first time |  |
eyes-end | Slice after the end of the eye socket | ✓ | $l_{3}$
head-end | Last slice where bone and brain can be seen |  |
L5 | Centroid of 5. Lumbar vertebrae | ✓ | $l_{4}$
L4 | Centroid of 4. Lumbar vertebrae |  |
L3 | Centroid of 3. Lumbar vertebrae | ✓ | $l_{5}$
L2 | Centroid of 2. Lumbar vertebrae |  |
L1 | Centroid of 1. Lumbar vertebrae | ✓ | $l_{6}$
Th12 | Centroid of 12. Thoracic vertebrae |  |
Th11 | Centroid of 11. Thoracic vertebrae | ✓ | $l_{7}$
Th10 | Centroid of 10. Thoracic vertebrae |  |
Th9 | Centroid of 9. Thoracic vertebrae |  |
Th8 | Centroid of 8. Thoracic vertebrae | ✓ | $l_{8}$
Th7 | Centroid of 7. Thoracic vertebrae |  |
Th6 | Centroid of 6. Thoracic vertebrae |  |
Th5 | Centroid of 5. Thoracic vertebrae | ✓ | $l_{9}$
Th4 | Centroid of 4. Thoracic vertebrae |  |
Th3 | Centroid of 3. Thoracic vertebrae |  |
Th2 | Centroid of 2. Thoracic vertebrae | ✓ | $l_{10}$
Th1 | Centroid of 1. Thoracic vertebrae |  |
C7 | Centroid of 7. Cervical vertebrae |  |
C6 | Centroid of 6. Cervical vertebrae | ✓ | $l_{11}$
C5 | Centroid of 5. Cervical vertebrae |  |
C4 | Centroid of 4. Cervical vertebrae |  |
C3 | Centroid of 3. Cervical vertebrae |  |
C2 | Centroid of Axis vertebrae |  |
C1 | Centroid of Atlas vertebrae | ✓ | $l_{12}$

####  4.1.3 Learning Procedure

In this thesis, a body part regression model is trained. The learning
procedure is based on the SSBR model proposed by Yan et al.
[yan2018unsupervised, yan2018deep]. At the end of this section, a detailed
overview of the differences to the original method will be given.

The final body part regression model $f_{\boldsymbol{\theta}}$ with parameters
$\boldsymbol{\theta}$ should map CT slices $\mathbf{X}_{ij}$ of volume $i$ and
slice index $j$ to anatomical slice scores $s_{ij}$:

$f_{\boldsymbol{\theta}}:\mathbb{R}^{n\times n}\rightarrow\mathbb{R}:\quad
f_{\boldsymbol{\theta}}(\mathbf{X}_{ij})=s_{ij}.$

The same anatomical region should map to the same slice score, independently
of patient and study. This means that our function $f_{\boldsymbol{\theta}}$
should be invariant under certain transformations, which do not change the
visible anatomy but may change study- or patient-related features. Moreover,
each body part should be assigned to a unique slice score. Mathematically
speaking, we want to learn an injective function which maps anatomical regions
to slice scores. Therefore, we need a function that produces slice scores that
increase strictly monotonic with slice index. For example, the patient’s
pelvis should map to a lower slice score than the head of a patient. Finally,
we have two main conditions for our mapping function
$f_{\boldsymbol{\theta}}$:

  1. 1\. Independence condition: The same anatomical region should map to the same score, independently of patient and study.

  2. 2\. Monotony condition: The anatomical slice scores should increase monotonically with height.

With these two conditions, we can now construct an objective function that
should be minimized during training. For training a body part regression
model, a self-supervised approach is used. Therefore, we do not need ground
truth labels, and we can use a large and diverse dataset. In the self-
supervised approach, we gain a learning incentive from the intrinsic order of
consecutive axial slices in a three-dimensional CT volume. Through the
relative relation to other slices in the human body, the model can learn the
relation between anatomy and slice score on its own without labeled data.

In the following, we will try to find a general approach for a possible
objective function, which considers the monotony condition and the
independence condition. As a first approach for the objective function, the
following equation can be stated

$\begin{split}&\Phi=\sum_{i=1}^{N}\sum_{j=1}^{M_{i}-1}L_{\text{order}}(\Delta
s_{ij}),\quad\Delta s_{ij}=s_{ij+1}-s_{ij}\end{split},$

where $N$ is the number of available CT-volumes and $M_{i}$ the number of
slices of volume $i$. The order loss $L_{\text{order}}$ should be designed
such that it is high if $\Delta s_{ij}<0$ and low if $\Delta s_{ij}>0$. It
makes sure that the slice scores increase monotonically with slice height. In
section 4.1.4 two possible order loss functions will be derived. This
objective function $\Phi$ can be seen as a potential, which we want to
minimize to obtain good parameters $\boldsymbol{\theta}$ for our model
$f_{\boldsymbol{\theta}}$, which will be a neural network in this work.
Stochastic gradient descent with a mini-batch of size $B$ can be used to
minimize the objective function. Assuming we are choosing $m$ ($m\geq 2$)
equidistant random slices per volume and $B$ random CT-volumes per update
step, we need to calculate the following loss function, which is an
expectation of our objective function, to calculate the gradients:

$L=\frac{1}{B}\sum_{i=1}^{B}\frac{1}{m-1}\sum_{j=1}^{m-1}L_{\text{order}}(\Delta
s_{ij}).$

The $L_{\text{order}}$ objective function is motivated by the monotony
condition. If the chosen neural network model has enough parameters
$\boldsymbol{\theta}$, an easy way for the neural network to reduce the loss
is to simply learn the slice score for each volume separately. Then the
independence condition would be violated and anatomical similar regions from
different patients would not necessarily be mapped to the same slice score. We
need additional constraints to incorporate the second condition into the model
as well. Loosely speaking, we want that for any transformation $g$ that does
not change the assigned anatomical region of an image, the slice score mapping
function should produce the same output. Let $G$ denote the family of
transformations $g$ for which this assumption is true. We want that our slice
to slice score mapping function $f$ is invariant under transformations $g\in
G$:

$f_{\theta}(g(\mathbf{X}_{ij}))=f_{\theta}(\mathbf{X}_{ij})=s_{ij}.$

In the machine learning setup, we can try to approximately implement this
constraint by using transformations $g\in G$ as data augmentation. These data
augmentations have to be applied independently to each individual sampled
slice $\mathbf{X}_{ij}$. With slice-wise data augmentation, we try to prevent
the network from overfitting.

An example of the used network structure can be seen in figure 4.2. The
training procedure is as follows: first, $m$ equidistant slices are sampled
from a random three-dimensional volume. The distance between the slices
$\Delta h$ is sampled randomly from a $z$-range, e.g., from 5 mm to 100 mm.
After sampling $\Delta h$, the distance in mm is converted to a difference in
slice indices $k$ with the help of the z-spacing $z$ by the equation

$k=\text{round}\left(\frac{\Delta h}{z}\right).$ (4.1)

For each of the $m$ two-dimensional slices, data augmentation techniques are
applied independently. After that, the slice scores are computed for each
slice by a neural network. Then, based on the stack of predicted slice scores,
the loss function is computed. Finally, the total loss is backpropagated to
the neural network so that the model can learn from its mistakes, and the
procedure starts all over again.

Figure 4.2: Overview of the learning framework for the self-supervised body
part regression model. In this simplified illustration, the mini-batch size is
equal to one. The abbreviation DA stands for diverse data augmentation
techniques. The variable $k$ is the randomly sampled slice distance between
the slices. The network architecture was adapted from the SSBR model
[yan2018deep, yan2018unsupervised].

The network architecture was adapted from the SSBR model [yan2018deep,
yan2018unsupervised] with the following modifications:

  1. 1.

In the SSBR model, the slice-index difference $k$ is sampled directly. We will
sample the distance $\Delta h$ in mm from the physical space between $[\Delta
h_{\text{min}},\Delta h_{\text{max}}]$. From this sampled distance $\Delta h$,
the slice-index difference $k$ is computed. Because the distribution of the
z-spacing in the training data is heterogeneous, the two approaches lead to
different samplings strategies. The sampling method was changed to obtain a
more natural and realistic slice sampling and to be less dependent on the
distribution of the z-spacing in the dataset.

  2. 2.

In the SSBR model, random cropping was used as the only data augmentation
method [yan2018deep]. This work will apply far more data augmentation
techniques to prevent the model from overfitting on studies or patients. The
used data augmentation techniques can be find in section 4.1.5.

  3. 3.

The original SSBR model was only trained on body regions from the lower pelvis
to the upper chest [yan2018unsupervised]. In this work, the investigated body
region was extended to the end of the head by using publicly available data.
This makes the method applicable to more CT images.

  4. 4.

In the SSBR model, the classification order loss was used (see sec. 4.1.4). We
will use the heuristic order loss. The derivations of the loss functions and
the use of the heuristic order loss makes the used hyperparameters better
interpretable (see sec. 4.1.4).

  5. 5.

The original SSBR model is validated by transferring the regression problem to
a 3-class classification problem and measuring the accuracy. In this work, a
metric is proposed that allows a more accurate evaluation of different models
without simplifying the problem to a classification problem (see sec. 4.2).

####  4.1.4 Order Loss

For self-supervised approaches, the loss function is crucial. The loss
function of the body part regression model consists of an order loss and
optionally an additional distance loss. The order loss ensures the monotony
condition, which states that slice scores should increase with height. In this
section, two versions of the order loss are introduced. The classification
order loss used by Yan et al. [yan2018unsupervised] is mathematically derived,
and a new variant, called heuristic order loss, is proposed. The
classification order loss is based on transferring the order problem to a
classification problem. For the heuristic order loss, we additionally
incorporate the physical distance between sampled slices into the loss
function.

Figure 4.3: Sampled CT slices $\mathbf{X}_{v},\mathbf{X}_{w}$, corresponding
slice scores $s_{v},s_{w}$ and physical distance $\Delta h$ of CT volume with
z-spacing $z$.

For the following calculations, we assume that we sample two two-dimensional
slices $\mathbf{X}_{v}$ and $\mathbf{X}_{w}$ from a three-dimensional human
CT-volume with the corresponding physical distance $\Delta h$, as it can be
seen in figure 4.3.

Classification Order Loss:  Because no mathematical derivation for the loss
function was given in the literature, the loss function will be derived in
this subsection. The key idea behind this order loss function is to transfer
the order problem to a classification problem. Therefore, the derivation will
be quite similar to the derivation of the Cross-Entropy loss (see sec. 2.1.3).
The ordering problem is converted to a classification problem by defining two
classes $c\in\\{0,1\\}$:

$c_{vw}=\begin{cases}1&\text{if }v>w,\\\ 0&\text{if }v\leq w.\\\ \end{cases}$

If slice $\mathbf{X}_{v}$ lies above slice $\mathbf{X}_{w}$, we assign the
class 1. Otherwise, the class 0 is assigned. The slices $\mathbf{X}_{v}$ and
$\mathbf{X}_{w}$ are randomly sampled from one patient. We assume that the
class labels $c$ are Bernoulli distributed

$\begin{split}&c_{vw}\sim\mathcal{B}(\pi_{vw}),\\\
&P(c=1|\mathbf{X}_{v},\mathbf{X}_{w})=\pi_{vw},\\\
&P(c=0|\mathbf{X}_{v},\mathbf{X}_{w})=1-\pi_{vw},\\\ &\Rightarrow
P(c|\mathbf{X}_{v},\mathbf{X}_{w})=\pi_{vw}^{c_{vw}}\cdot(1-\pi_{vw})^{1-c_{vw}}.\end{split}$

Provided that we have a perfect body part regression model, which fulfills the
monotony condition (see sec. 4.1.3), we can model the probability $\pi_{vw}$
based on the predicted scores $s_{v}$ and $s_{w}$. We assume that slice
$\mathbf{X}_{v}$ lies above slice $\mathbf{X}_{w}$, if the score difference
$\Delta s_{vw}=s_{v}-s_{w}$ is greater than zero. If the score difference
$\Delta s_{vw}$ is less than zero, we assume, that slice $\mathbf{X}_{v}$ lies
under slice $\mathbf{X}_{w}$. For $\Delta s_{vw}=0$, we expect a chance of 50
% that $\mathbf{X}_{v}$ lies above slice $\mathbf{X}_{w}$. The probability
$\pi_{vw}$ can be represented with the help of the sigmoid function $\sigma$:

$\pi_{vw}=\sigma(\Delta s_{vw}).$ (4.2)

The sigmoid function $\sigma$ is given by equation 3.4. The greater $s_{v}$
compared to $s_{w}$, the higher is the probability $P(c_{vw}=1)$. To train a
model, which fulfills the monotony condition, we can use the maximum
likelihood approach. With this approach, we find the model parameters
$\boldsymbol{\theta}$ under which the observed slice orderings are most
likely.

We can influence the observed slice orderings by the sampling strategy of the
input volumes. By sampling a random start index $j$ and a random index
difference $k$, we can define

$v=j+k,\quad w=j.$

Therefore, we know that slice $\mathbf{X}_{v}$ always lies above slice
$\mathbf{X}_{w}$ and that the observed class label is always one. If we sample
from $N$ independent volumes, the observed data is given by
$(\textbf{X}_{i,j},\textbf{X}_{i,j+k},c_{i}=1)_{i=1...N}$. Moreover, the
likelihood function is equal to

$\mathcal{L}({\boldsymbol{\theta}})=p(\\{c_{ij}=1\\}_{i=1...N}|\boldsymbol{\theta})=\prod_{i=1}^{N}\sigma(f_{\boldsymbol{\theta}}(\mathbf{X}_{i,j+k})-f_{\boldsymbol{\theta}}(\mathbf{X}_{ij}))=\prod_{i=1}^{N}\sigma(\Delta
s_{ij}),$ (4.3)

where the slice score difference $\Delta s_{ij}$ is defined by

$\Delta s_{ij}=s_{i,j+k}-s_{ij},\quad
f_{\boldsymbol{\theta}}(\mathbf{X}_{ij})=s_{ij},\quad
f_{\boldsymbol{\theta}}(\mathbf{X}_{i,j+k})=s_{i,j+k}.$

By taking the logarithm of equation 4.3 we obtain the log-likelihood

$\text{ln }\mathcal{L}({\theta})=\sum_{i=1}^{N}\text{ln}\left(\sigma(\Delta
s_{ij})\right).$

The order loss is given by an estimate of the normalized negative log-
likelihood d

$L_{\text{order}}=-\frac{1}{B}\sum_{i=1}^{B}\text{ln}\left(\sigma(\Delta
s_{ij})\right),$ (4.4)

where $B$ is the batch size. This is equivalent to the classification order
loss from the SSBR model [yan2018unsupervised] for $m=2$. For $m>2$, we sample
$m$ slices from a volume by choosing a random slice difference $k$ and a
random start index. We obtain the following normalized order loss function of
the SSBR model

$L_{\text{order}}^{c}=-\frac{1}{B}\frac{1}{m-1}\sum_{i=1}^{B}\sum_{j=1}^{m-1}\text{ln}\
\sigma(\Delta s_{ij}),\qquad\Delta s_{ij}=s_{i,j+1}-s_{ij}.$ (4.5)

Figure 4.4 shows that the loss and the gradients are small if $\Delta s$ is
greater than zero, and it is high if $\Delta s$ is smaller than zero. One
weakness of this loss is, that for slice pairs with small physical distance
$\Delta h<$ 5 mm we would assume that the slice scores are nearly the same,
because the images are looking quite similar. Although, we expect
$s_{ij}\simeq s_{i,j+1}$, the gradient of the loss at $\Delta s=0$ is still
high. This weakness can be overcome by introducing a minimum physical slice
distance $\Delta h_{\text{min}}$ for sampling $\mathbf{X}_{i,j}$ and
$\mathbf{X}_{i,j+1}$.

Figure 4.4: Trend of the classification order loss $L_{\text{order}}(\Delta
s)$ from equation 4.4 for $B=1$ and $m=2$. For $\Delta s>>0$, the loss
converges to zero. For $\Delta s<0$ the loss and the gradient is high, the
monotony condition is not met.

Heuristic Order Loss:  In the previous paragraph, we derived the
classification order loss used by Yan et al. to train the SSBR model
[yan2018unsupervised]. We will now investigate a completely different and new
approach for handling the order problem. The main reason behind the
introduction of the heuristic order loss is to incorporate the physical
distance $\Delta h$ between slices in a meaningful way. In the classification
order loss, the physical distance has no impact on the loss function. If two
slices lie physically close to each other, they would be treated in the same
way as two slices that lie far apart from each other. In the heuristic order
loss, we concentrate on building a heuristic for the following probability:

$p(c_{vw}=1|\Delta s_{vw})=\pi_{vw},\quad\Delta s_{vw}=s_{v}-s_{w}.$

The probability $\pi_{vw}$ is the probability that slice $\mathbf{X}_{v}$ lies
above slice $\mathbf{X}_{w}$, based on the predicted slice scores $s_{v}$ and
$s_{w}$. We can model the predicted probability $\hat{\pi}_{vw}$ through the
network with the sigmoid function as already derived in equation 4.2:

$\sigma(\Delta s_{vw})=\hat{\pi}_{vw},$

In general, the probability $\hat{\pi}_{vw}$ is high if the difference between
the predicted slice scores is as well high. Next, we want to construct a
heuristic for the target probability $\pi_{vw}$. To build this pseudo-label,
we assume that it should be easier to decide whether slice $\mathbf{X}_{v}$
lies above slice $\mathbf{X}_{w}$ if the physical slice distance $\Delta h$
between both slices is high. Let us consider two scenarios in which the goal
is to order two slices from a human CT image:

  1. 1.

One slice lies in the lung and the other one in the abdomen region. The
physical slice difference between both slices is equal to $\Delta h=$ 200 mm.

  2. 2.

The two slices are consecutive slices from the abdomen region, with a physical
slice difference of $\Delta h=$1 mm.

A human can order the slices from scenario 1 with nearly a hundred percent
certainty, whereas in scenario 2, a human would be less sure about the order
of the slices. Based on these considerations, we propose as a proxy for
$\pi_{vw}$

$\sigma(\beta\Delta h_{vw})=\pi_{vw}^{\star},\quad\Delta h_{vw}=h_{v}-h_{w},$

where $\beta$ is a fixed hyperparameter. We assume high certainty for large
$\Delta h$ (around 100 %) and low certainty for small $\Delta h$ (around 50
%). Moreover, we assume that the scores are continuous. For training the
model, we can try to minimize the distance between $\pi_{vw}^{\star}$ and
$\hat{\pi}_{vw}$ through the least square error method
[durstewitz2017advanced]. If we chose for every volume $i$ in a batch a random
start index $j$ and a random slice index difference $k$ with the corresponding
physical slice distance $\Delta h_{i}$, we gain the following least square
error

$L_{\text{order}}^{h}=\frac{1}{B}\sum_{i=1}^{B}||\pi_{ij}^{\star}-\hat{\pi}_{ij}||^{2}_{2}=\frac{1}{B}\sum_{i=1}^{B}||\sigma(\beta\Delta
h_{i})-\sigma(\Delta s_{ij})||^{2}_{2},$ (4.6)

where $B$ is the mini batch size.

Figure 4.5: Order loss from equation 4.6 for different $\Delta h$, $B=1,$ and
$\beta=\text{0.01/mm}$.

For $m>2$, we sample $m$ slices from a volume by choosing a random start index
and slice index difference $k$. We obtain the following loss function:

$L^{h}_{order}=\frac{1}{B(m-1)}\sum_{i=1}^{B}\sum_{j=1}^{m-1}||\sigma(\beta\Delta
h_{i})-\sigma(\Delta s_{ij})||^{2}_{2},\quad\Delta s_{ij}=s_{i,j+1}-s_{ij}.$
(4.7)

The main advantage of this loss function is a meaningful integration of
$\Delta h$. Figure 4.5 shows that the network can now learn something if
$\Delta s$ is too high for small $\Delta h$: For small $\Delta h$, we want
$\Delta s$ to be small as well because the visible slices look similar. If
$\Delta s$ is too big, the network obtains the feedback that the predicted
slice scores are too far apart. While we assume a linear relation between
$\Delta s$ and $\Delta h$ in the range of small $\beta\Delta h$, the relation
disappears for large $\beta\Delta h$, and the loss becomes more similar to the
classification order loss. The order loss vanishes at the minimum for $\Delta
s_{ij}$ equal to

$L_{\text{order}}=0\quad\Rightarrow\quad\sigma(\beta\Delta
h_{i})=\sigma(\Delta s_{ij})\quad\Rightarrow\quad\Delta s_{ij}=\beta\Delta
h_{i}.$ (4.8)

To understand the influence of parameter $\beta$ better, it is helpful to look
at the curvature of the minimum. At the minimum $\Delta s_{ij}=\beta\Delta
h_{i}$ for $B=1$ and $m=2$, the curvature is given by

$L_{\text{order}}^{{}^{\prime\prime}}|_{\Delta s_{ij}=\beta\Delta
h_{i}}=\frac{2e^{-2\beta\Delta h_{i}}}{(1+e^{-\beta\Delta h_{i}})^{4}}.$ (4.9)

The curvature is the second derivative of the loss function with respect to
the score difference $\Delta s_{ij}$. If the curvature is high, the change in
the gradients around the minimum is big, and the linear constraint is enforced
more strongly by gradient descent. For $\beta\Delta h$ around zero, the
curvature lies near $\frac{1}{8}$ and for $\beta\Delta h\rightarrow\infty$ the
curvature goes to zero. This illustrates that the linear relation between
$\Delta h$ and $\Delta s$ is only enforced by gradient descent for small
$\beta\Delta h$ and disappears for large $\beta\Delta h$. For large
$\beta\Delta h$, the curvature around the minimum is approximately zero, and
we are in a flat region. If $\Delta s$ is large enough, the gradients are
nearly zero, and the neural network has only a particular weak incentive to
find the analytical minimum.

Figure 4.6 b) visualizes the curvature at the minimum for different $\beta$.
For $\Delta h=0\ \text{mm}$, the curvature is equal for all $\beta$. For
$\Delta h>0\ \text{mm}$, the $\beta$ parameter influences how quickly the
curvature decreases. In figure 4.6 a), we observe as well that for $\Delta
h=0\ \text{mm}$ the different $\beta$ have no impact. For $\Delta h=2\
\text{cm}$, the curvature is smaller for bigger $\beta$. Therefore, we can say
that the $\beta$ parameter influences how long the linearity condition $\Delta
s\propto\Delta h$ persists.

Figure 4.6: Analysis of parameter $\beta$ in heuristic order loss. In figure
a), the heuristic order loss for different $\beta$ is visible. The graph was
shifted to the left by the minimum $\Delta s=\beta\Delta h$ so that the
minimum lies at $\Delta s=0$ for all curves. In figure b), the curvature of
the heuristic order loss for different $\beta$ at the minimum $\Delta
s=\beta\Delta h$ is visible.

####  4.1.5 Data Augmentations

Data augmentations, applied to each slice independently, play a key role in
the learning procedure. The data augmentation techniques take care that the
independence condition is met (see sec. 4.1.3). Table 4.4 summarizes the used
data augmentations and parameters. In the following, the data augmentation
methods applied throughout the experiments are explained. The hyperparameters
for the data augmentations were set manually to gain more variability in the
image quality, appearance, and intensity values in the input image without
receiving unrealistic CT slices.

Gaussian Noise:  Gaussian noise was randomly applied to the images. The
maximum possible standard deviation of the Gaussian noise $\sigma_{max}=0.04$
corresponds to $50\text{HU}$, because the input volumes are rescaled from
(-1000 HU, 1500 HU) to (-1,1):

$0.04\cdot 2500\text{ HU}/2=50\text{ HU}.$ (4.10)

This data augmentation was applied randomly with a chance of $p=50\%$ to the
images in the dataset.

Gaussian Blur:  Gaussian blur was randomly applied with a chance of $p=50\%$
to the images. The maximum kernel size for blurring the input image is
randomly chosen from three to seven. The maximum standard deviation for the
blurring filter is set to 0.5.

Flip, Scale and Rotate:  In 50 % of the cases, the image is randomly flipped
horizontally or vertically. The images are transposed with a probability of 50
%. By these two operations, mirroring and rotating the image by 90 degrees are
covered. Furthermore, the images are randomly scaled by a scale factor between
80 % and 120 % and randomly rotated between zero and ten degrees.

Contrast and Brightness:  The contrast is changed by linearly scaling the
pixel intensities by the scaling factor, which is randomly chosen between 80 %
and 120%. If the image extends the minimum or maximum possible value of -1 and
1, the extreme values are clipped. The brightness of the image is changed by
randomly adding a constant between zero and 0.08. The maximum shift limit is
equal to $100\text{ HU}=0.08\cdot 2500\text{ HU}/2$ (see eq. 4.10). Contrast
and brightness adaption are applied randomly with a chance of $p=50\%$ to the
images in the dataset.

Figure 4.7: Example of Add Frame data augmentation. In the left figure, the
original slice is visible. In the figure on the right, the Add Frame data
augmentation technique was applied by adding a circular mask to the image. As
a result, the CT table and the arms are not anymore visible.

Add Frame:  In the Add Frame data augmentation technique, a circular mask is
added to the input image. The transformation is applied randomly with a
probability of $p=25\%$ to the images in the dataset. The circle has a radius
of about 48 pixels. Figure 4.7 visualizes an example of this data augmentation
method based on a chest CT slice. Moreover, this figure visualizes that the
data augmentation method helps ignore arms and CT tables in CT slices. Some
patients have their arms positioned next to their thighs and others over their
heads. In addition, the CT table appearance differs from device to device.
With this data augmentation, the arms and the CT table are sometimes cut off,
which helps the network to ignore these features for predicting the slice
scores.

Table 4.4: Documentation of used data augmentation techniques and corresponding parameters for training a body part regression model. Several augmentation techniques were adopted from the python package Albumentation [albumentation].  Data augmentation | Parameters | Implementation
---|---|---
GaussianNoise | std_min=0, std_max=0.04, p=0.5 | custom
GaussianBlur | blur_limit=(3,7), sigma_limit=0.5, p=0.5 | Albumentation [albumentation]
ShiftScaleRotate |  shift_limit=0, scale_limit=0.2, rotate_limit=10, border_mode=reflect | Albumentation [albumentation]
Flip | p=0.5 | Albumentation [albumentation]
Transpose | p=0.5 | Albumentation [albumentation]
Contrast | scale_delta=0.2, p=0.5 | custom
Brightness | shift_limit=0.08, p=0.5 | custom
Add Frame | diameter=0.75, p=0.25 | custom

###  4.2 Evaluation Metrics

Evaluating a body part regression model is non-trivial. Because of a lack of
ground truth labels, it is not easy to compare model performances. The mapping
between anatomy and score is learned implicitly in an unsupervised manner,
resulting in different slice score scales for each model. For example, one
model may map the beginning of the pelvis to -5 and the end of the head to 6,
and another model may map the beginning of the pelvis to -100 and the end of
the head to 90. Even for a human, it is hard to define a slice-wise anatomy-
to-score mapping that fulfills the monotony and independence condition from
sec. 4.1.3. One reason for this is the inter-patient anatomical variability.
For example, the liver can start at the end of the kidneys or the end of the
pelvis bones. The lungs can be bigger or smaller based on the breathing cycle
[vauthey2002body, abdalla2004total]. A good evaluation metric is necessary to
measure how well a body part regression model fulfills the independence
condition while allowing for cross-model comparison. In the following, two
different metrics are described, which are used to evaluate the model
performance of the body part regression model. The proposed Landmark Mean
Square Error (LMSE) is introduced in section 4.2.1. Moreover, the setup for
the accuracy metric, which was already used for evaluation by Yan et al.
[yan2018unsupervised], is explained in section 4.2.2.

Before calculating the defined metrics, we have to compute the model specific
slice score reference table, which establishes a correspondence between each
evaluation landmark defined in table 4.3 and its slice score pseudo-label
$\bar{s}$. The pseudo-labels $\bar{s}$ are calculated for each of the twelve
defined evaluation landmarks by averaging across all annotated landmark
positions of the training set

$\bar{s}_{l_{k}}=\frac{1}{N_{l_{k}}}\sum_{i=1}^{N_{l_{k}}}s_{il_{k}},\qquad\forall
l_{k}\in\\{l_{1},...,l_{12}\\},$ (4.11)

where $N_{l_{k}}$ indicates the number of CT volumes where the landmark
$l_{k}$ was annotated. The pseudo-labels are calculated with an annotated
subset of the training dataset to be independent of the validation and test
dataset.

####  4.2.1 Landmark Mean Square Error

The Landmark Mean Square Error (LMSE) $\bar{\phi}$ is used as the primary
metric to compare different models. It is a normalized version of a Mean
Square Error (MSE) with the difference that we do not have ground truth labels
and therefore create pseudo-labels based on the annotated landmarks. As
pseudo-labels $\bar{s}_{l_{k}}$, we are using the values of the slice score
reference table. Further, an additional normalization is applied to enable
comparing different models with different slice score scales. The landmark
scores $s$ are transformed into $s^{\star}$ such that the pseudo-label
$\bar{s}_{l_{1}}$ of landmark $l_{1}$, which refers to the start of the
pelvis, maps to 0 and the pseudo-label $\bar{s}_{l_{12}}$ of landmark
$l_{12}$, which refers to the end of the eye socket, maps to 100. The
transformation is:

$s\mapsto
s^{\star}=\frac{100\cdot(s-\bar{s}_{l_{1}})}{\bar{s}_{l_{12}}-\bar{s}_{l_{1}}}=\frac{(s-\bar{s}_{l_{1}})}{d},\qquad
d=\frac{1}{100}(\bar{s}_{l_{12}}-\bar{s}_{l_{1}}).$ (4.12)

In the following, the constant $d$, will be referred to as the normalization
constant. The Landmark Mean Square Error per volume is defined as

$\begin{split}\phi_{i}&=\frac{1}{|L_{i}|}\sum_{l_{k}\in
L_{i}}(\bar{s}^{\star}_{l_{k}}-s^{\star}_{il_{k}})^{2}\\\
&=\frac{1}{|L_{i}|}\sum_{l_{k}\in
L_{i}}\left(\frac{\bar{s}_{l_{k}}-\bar{s}_{l_{1}}-(s_{il_{k}}-\bar{s}_{l_{1}})}{d}\right)^{2}\\\
&=\frac{1}{|L_{i}|}\sum_{l_{k}\in
L_{i}}\left(\frac{\bar{s}_{l_{k}}-s_{il_{k}}}{d}\right)^{2},\end{split}$
(4.13)

where $L_{i}$ is the set of annotated landmark positions for volume $i$. The
Landmark Mean Square Error for the whole dataset $\bar{\phi}$ is the mean of
the LMSEs per volume $\phi_{i}$

$\begin{split}\bar{\phi}=\frac{1}{N}\sum_{i=1}^{N}\phi_{i},\end{split}$ (4.14)

where $N$ is the number of volumes in the annotated dataset. The average over
the LMSEs per volume is taken, because we are interested in how well the model
generalizes to new volumes [holland2020drawing]. The standard error of the
LMSE $\sigma_{\bar{\phi}}$ can be estimated by

$\sigma_{\bar{\phi}}=\frac{\sigma_{\phi}}{\sqrt{N}}=\frac{1}{\sqrt{N}}\cdot\sqrt{\frac{1}{N-1}\sum_{i=1}^{N}(\phi_{i}-\bar{\phi})^{2}}.$
(4.15)

The LMSE is specialized in evaluating the independence condition, but it is
insensitive to the monotony condition. For example, the LMSE can be zero if
all slices are mapped to the same slice score. Therefore, we need to look at
the predicted scores of a model to be sure that the model also fulfills the
monotony condition.

Comparing Landmark Mean Square Errors: In the following, the choice of
significance test to compare different model performances based on the LMSE
will be explained. If the number of annotated volumes $N$ is high enough, a
normal distributed LMSE value can be assumed
$\mathcal{N}(\bar{\phi},\sigma_{\bar{\phi}})$ (based on the central limit
theorem [DeepVision]). Because of that, a standard two-sided z-test to compare
model $x$ with model $y$ can be performed [durstewitz2017advanced]. To test if
model $x$ performs significantly different to model $y$ we can state the
following null hypothesis $H_{0}$ and alternative hypothesis $H_{1}$:

  1. $H_{0}$:

$\mu_{\phi_{x}}=\mu_{\phi_{y}}$ ,

  2. $H_{1}$:

$\mu_{\phi_{x}}\neq\mu_{\phi_{y}}$.

In this thesis, the significance level will be set for all z-tests to 5%. The
test statistic $t$ is computed by calculating the $\sigma$-deviation to zero
of the difference between the performances $\bar{\phi}_{x}$ and
$\bar{\phi}_{y}$ [durstewitz2017advanced]:

$t=\frac{\bar{\phi}_{x}-\bar{\phi}_{y}}{\sqrt{\sigma_{\bar{\phi}_{x}}^{2}+\sigma_{\bar{\phi}_{y}}^{2}}}\sim\mathcal{N}(0,1).$
(4.16)

The test statistic is standard normal distributed if the null hypothesis
$H_{0}$ is true. If the absolute value of the test statistic is greater than
1.96, the null hypothesis $H_{0}$ can be rejected based on the significance
level, and we assume that model $x$ and model $y$ perform significantly
different.

####  4.2.2 Accuracy

Because of the easy interpretability of the accuracy and its usage in previous
work, we will also use the accuracy $\psi$ as an evaluation metric. For
calculating the accuracy, five classes were defined: pelvis, abdomen, chest,
neck, and head. Table 4.5 lists the class boundaries based on the annotated
landmarks. For evaluation, the body range between the lowest landmark pelvis-
start and the highest landmark eyes-end is considered. The accuracy $\psi$
compares the predicted classes $\hat{c}$ with the ground truth classes $c$ and
is defined in equation 3.7. The ground truth class $c_{ij}$ for volume $i$ and
slice index $j$ is determined by finding the start and end landmarks from
table 4.5 that are surrounding slice $j$. Only if both boundary slice indices
of a class are defined within a CT volume, the ground truth label was
assigned. Otherwise, the class assignment would not be unique.

The predicted class label $\hat{c}_{ij}$ for volume $i$ and slice index $j$ is
determined by using the predicted slice score $s_{ij}$ of the slice $j$. With
the help of the slice score reference table, the surrounding landmark slice
score boundaries from table 4.5 are found and the corresponding body part
class $\hat{c}_{ij}$ is assigned.

Table 4.5: Definition of five body parts to determine the accuracy. Body part | Start landmark | End landmark
---|---|---
pelvis | pelvis-start | L5
abdomen | L5 | Th11
chest | Th11 | Th2
neck | Th2 | C1
head | C1 | eyes-end

###  4.3 Application and Deployment

In this thesis, three use cases (UC) for the body part regression model will
be discussed:

  1. UC1

Estimate Examined Body Part: the body part regression model can be used for a
better estimate of the BodyPartExamined DICOM tag [fedorov2016dicom]. This
information is beneficial for sorting and filtering medical datasets
[yan2018deep, gueld2002quality]. This use case will be discussed in section
4.3.1.

  2. UC2

Known Region Cropping: Deep learning methods in the medical imaging domain are
typically designed for one specific body region. For a robust transfer of deep
learning models into the clinic, it is important to only select images which
fit the scope of the algorithms. Selecting and cropping images to the scope of
an algorithm can be done by the body part regression model. Moreover, the body
part regression model can be used as a post-processing step to detect false-
positive predictions of a medical algorithm, which lie out of scope. This
application use case is demonstrated by using medical segmentation algorithms
as a representative for medical computer vision algorithms. The use case is
explained in section 4.3.2.

  3. UC3

Data Sanity Checks:  The body part regression model can be used for basic data
sanity checks, such as validating the z-spacing and the ordering of the
z-axis. This use case will be discussed in section 4.3.3.

In this thesis, the trained body part regression model is specialized on CT
images. Therefore, the use cases will be evaluated on CT images. To make the
use cases easier applicable, the methods were deployed through a python
package and integrated into Kaapana (see sec. 4.3.4).

Predicted slice scores are post-processed by filtering empty slices, linearly
transforming, and smoothing scores, and removing outliers. The different steps
are summarized in figure 4.8. To obtain the cleaned slice scores from figure
4.8, the following steps need to be taken:

  1. 1.

Filter $s_{0}$: For the body part regression model, the slice score prediction
$s_{0}$ of an empty slice needs to be determined in advance. In empty slices,
the pixels have all the HU value, which is equal to $-1024$ HU. No human and
no CT table is visible. Score predictions of empty slices $s=s_{0}$ are
filtered out.

  2. 2.

Transform:  As described in section 4.2 the slice score reference table is
employed to transform predicted slice scores to a common scale. The
calculation of the table is done as in section 4.2 with the difference that
the data foundation are the annotated volumes from the training and the
validation set to obtain a better estimate of the expected slice scores for
each landmark. With the help of the computed slice score reference table, the
slice scores are linearly transformed so that the slice of the pelvis-start
landmark maps to 0 and the landmark from eyes-end maps to 100 (see eq. 4.12).
This makes working with the slice scores more standardized and manageable.

  3. 3.

Smooth:  The slice scores are smoothed to remove local discontinuities. The
smoothing is done by applying a discrete convolution with a Gaussian filter
kernel

$\begin{split}&s(h)\mapsto(f\ast g)(h),\\\
&g(h)=\frac{1}{\sqrt{2\pi}\cdot\sigma}\cdot
e^{\frac{-h^{2}}{2\sigma^{2}}},\end{split}$

where $f(h)$ is our neural network which maps the slice at height $h$ of a
patient to a slice score and where $g(h)$ is the Gaussian filter kernel. The
standard deviation $\sigma$ of the Gaussian filter is set manually to
$\sigma=\text{10 mm = 1 cm}$. This means that approximately a range of
$\pm\text{1 cm}$ around a score has an influence on the smoothed score and
should be in a comparable slice score range. Through this step, we take
advantage of the assumption that the function should be continuous.

  4. 4.

Possibly remove tails:  Remove tails of the slice score curve before 0 or
after 100 if tangential slopes $m_{t}$ are atypical. The tangential slope at
index $i$ is calculated by

$m_{ti}=\frac{s_{i+1}-s_{i}}{h_{i+1}-h_{i}}.$ (4.17)

For the body part regression model, a minimum valid tangential slope
$m_{t,min}$ and a maximum valid tangential slope $m_{t,max}$ are set. If a
tangential slope $m_{t}$ before 0 or after 100 does not lie in the valid
range, the tails until this score are removed from the slice score curve. To
obtain $m_{t,min}$ and $m_{t,max}$ for a given model the distribution of the
tangential slopes $m_{t}$ for the training dataset was calculated. For this,
the score curve for each volume from the training dataset is transformed and
smoothed as described in step 1 till step 3, step 4 is ignored. The 0.5 %
quantile and the 99.5 % quantile of the distribution was set as $m_{t,min}$
and $m_{t,max}$.

For the data sanity checks, the mean slope $\bar{m}_{s}$ of the slice score
curves needs to be found. The difference between the previously defined
tangential slope $m_{t}$ and the slice score curve slope $m_{s}$ is that
$m_{s}$ represents the overall slope of a slice score curve for one volume
whereas the $m_{t}$ is the slope on one point. To obtain the mean slice score
curve slope $\bar{m}_{s}$, we are iterating through all volumes from the
training data. For each volume, the cleaned slice scores were predicted
through steps 1 till 4. To each slice score curves, a linear line is fitted

$s_{i}=m_{s,i}\cdot z\cdot i+b,$

where $z$ is the z-spacing of the volume. The slope $m_{i,s}$ is saved. In the
end, the mean $\bar{m}_{s}$ of the slopes $m_{i,s}$ across all training
volumes $i$ is used for the data sanity checks.

Figure 4.8: Slice score post-processing steps before using the slice scores
for application purposes.

####  4.3.1 Estimate Examined Body Part

To estimate the examined body part of a CT volume $X$ with $N$ slices, the
cleaned slice scores and the slice score reference table are used. The
estimated examined body part can be specified as a dictionary, called body
part examined dictionary $d$, or it can be specified as a single tag as in the
meta-data entry BodyPartExamined in a DICOM file. This tag will be referred to
as body part examined tag $t$ in the following section.

Body Part Examined Dictionary:  Six different core body regions were defined
to categorize the slices within a CT volume: legs, pelvis, abdomen, chest,
neck-and-shoulder, and head. In table 4.6, the defined body parts with the
corresponding start and end landmarks can be seen. The start and end landmarks
were manually set through visual inspection of the results. Some body regions
are overlapping, as for example, abdomen and chest. In the abdomen body part,
the abdominal organs should be visible as the liver or the kidneys. The chest
body part was defined from the lowest Thoracic spine to the highest Thoracic
spine. The body region leg has no start landmark and the body region head has
no end landmark. If head slices are in an image, the slices until the end of
the volume are assigned as head body region. If leg slices lie inside a CT
volume, all slices from index zero to the score $\bar{s}_{pelvis-start}$ are
assigned as leg body region. The body part examined dictionary $d$ summarizes
the slice indices for each body region:

$\begin{split}d=\\{&\text{legs: }\\{i|\infty<s_{i}<\bar{s}_{\text{pelvis-
start}}\wedge i\in[0,..,N]\\},\\\ &\text{pelvis: }\\{i|\bar{s}_{\text{pelvis-
start}}<s_{i}<\bar{s}_{\text{pelvis-end}}\wedge i\in[0,..,N]\\},\\\
&\text{abdomen: }\\{i|\bar{s}_{\text{L5}}<s_{i}<\bar{s}_{\text{Th8}}\wedge
i\in[0,..,N]\\},\\\ &\text{chest:
}\\{i|\bar{s}_{\text{Th12}}<s_{i}<\bar{s}_{\text{Th1}}\wedge
i\in[0,..,N]\\},\\\ &\text{neck-and-shoulder:
}\\{i|\bar{s}_{\text{Th3}}<s_{i}<\bar{s}_{\text{C2}}\wedge i\in[0,..,N]\\},\\\
&\text{head: }\\{i|\bar{s}_{\text{C5}}<s_{i}<\infty\wedge
i\in[0,..,N]\\}\\}.\end{split}$ (4.18)

With the body part examined dictionary, a detailed overview of the imaged body
parts is easily accessible.

Table 4.6: Defined start and end landmarks for each body region. The expected start and end slice scores $\bar{s}$ are used from the slice score reference table. Body part |  Start landmark |  End landmark |  Start slice score |  End slice score
---|---|---|---|---
legs | - | pelvis-start | $-\infty$ | $\bar{s}_{pelvis-start}$
pelvis | pelvis-start | pelvis-end | $\bar{s}_{pelvis-start}$ | $\bar{s}_{pelvis-end}$
abdomen | L5 | Th8 | $\bar{s}_{L5}$ | $\bar{s}_{Th8}$
chest | Th12 | Th1 | $\bar{s}_{Th12}$ | $\bar{s}_{Th1}$
neck-and-shoulder | Th3 | C2 | $\bar{s}_{Th3}$ | $\bar{s}_{C2}$
head | C5 | - | $\bar{s}_{C5}$ | $+\infty$

Body Part Examined Tag: For the body part examined tag $t$, one tag per CT
volume is assigned. The assigned body regions are PELVIS, ABDOMEN, CHEST,
NECK, and HEAD. If multiple body regions are visible, the tags are
concatenated by a hyphen. The legs are excluded from the tag because they are
out of the scope of the algorithm. If the data sanity check of section 4.3.3
outputs that the CT volume seems to be invalid based on the slope of the slice
score curve, the resulting tag is NONE. For all other volumes, we
differentiate between volumes with a small z-range of less than 10 cm and
volumes with a bigger z-range of more than 10 cm. An overview of the decision
tree can be found in figure 4.9.

For volumes with small z-ranges (z-range < 10 cm), no body region will likely
be completely visible. Therefore, for every slice, a body region will be
assigned based on the predicted cleaned slice score. After that, the most
frequent body tag will be assigned to the whole CT volume. The slice score
reference table is used to obtain the score boundaries for each class. The
defined landmark boundaries for each body part class can be found in figure
4.9.

For volumes with bigger z-ranges, it is more likely that multiple body ranges
were scanned. For each body part, a list of landmarks was defined through
visual inspection. If a volume includes at least a manually defined minimum
number of landmarks, the body part is declared as visible. A landmark $l_{k}$
is defined as visible if we find values in the cleaned slice scores which are
bigger and smaller than the expected landmark slice score $\bar{s}_{l_{k}}$.
For each body part, the related landmarks and the minimum number of required
landmarks can be found in figure 4.9. In the end, the visible body parts are
concatenated in descending order to one tag $t$ (see fig. 4.9).

Figure 4.9: Decision tree to obtain the body part examined tag $t$.

####  4.3.2 Known Region Cropping

Because of a limitation in labeled data, medical image algorithms are often
trained on a small dataset. These datasets normally only show the relevant CT
body part for the prediction task. For example, a CT kidney segmentation
network is only trained on abdomen CT images and has never seen head CT scans
during training [hellerkits19]. This can lead to miss-classifications in head
regions during test time. Therefore, clinical deep learning software can
benefit from a basic understanding of the visible anatomy in a CT volume to
filter false-positive predictions in irrelevant body areas.

The body part regression model can be used as a preprocessing step of clinical
end-to-end deep learning software to avoid false-positive predictions and to
reduce the run-time. Moreover, known region cropping can be used as a post-
processing step of a deep learning algorithm to catch false-positive
predictions. The known region cropping process is independent of the used
clinical algorithm. The different applications of known region cropping as
pre- or post-processing step are visualized in figure 4.10. In the figure, a
lung segmentation algorithm is used as an exemplary clinical computer vision
task.

Before known region cropping can be applied, we need to identify the known
region of a medical algorithm. The known region is defined by the minimum
known score $s_{\text{min}}$ and the maximum known score $s_{\text{max}}$.
There are two main approaches to identify the score boundaries:

  1. 1.

Through the data-driven approach, the known range of an algorithm is
determined with the used training data of the algorithm. In this approach, for
each volume of the training dataset, the cleaned slice score curve is
predicted (see fig. 4.8), and the minimum and maximum slice scores are
collected in two distinct lists. In the end, the $q_{\text{min}}$-quantile of
the list with the minimum slice-scores is used as the minimum known score
$s_{\text{min}}$. The $q_{\text{max}}$-quantile of the list with the maximum
slice scores is used as the maximum known score $s_{\text{max}}$. The minimum
quantile $q_{\text{min}}$ and the maximum quantile $q_{\text{max}}$-quantile
are hyperparameters, which are manually set.

  2. 2.

Through the manual approach, the known range of an algorithm is determined
manually based on the medical prediction task of the algorithm. It can be used
if the training data of a deep learning algorithm is unknown. The minimum
known score and maximum known score can be chosen based on the slice score
reference table of the model and the prediction task. For example, for a
kidney segmentation algorithm, the minimum known score can be set to the
expected score at the landmark femur-end, and the maximum known score can be
set to the expected score of the landmark Th12.

If we chose the difference $\Delta s=s_{\text{max}}-s_{\text{min}}$ too high,
it would lead to fewer detected false-positives. If we chose it too small, it
would lead to the truncation of correct predictions. We need to find a trade
of between the behaviors. In conclusion, for preparing known region cropping
for a model, the hyperparameters $q_{\text{min}}$ and $q_{\text{max}}$ needs
to be set if the data-driven approach is used. Additionally, we need to find
$s_{\text{min}}$ and $s_{\text{max}}$ for the model with the manual approach
or the data-driven approach. To apply the known region cropping method to a CT
volume $X$, we need to follow the steps:

  1. 1.

Predict slice scores of CT volume $X$.

  2. 2.

Calculate cleaned slice scores (see fig. 4.8).

  3. 3.

Crop prediction or volume $X$, where the score is outside the known region
defined by $s_{\text{min}}$ and $s_{\text{max}}$ (see fig. 4.10).

Figure 4.10: Visualization of known region cropping based on a lung
segmentation algorithm. False-positive predictions can arise if the algorithm
needs to process images that are out of the algorithm’s scope.

####  4.3.3 Data Sanity Checks

The height of women and men are approximately normal distributed with
$h_{\text{women}}=165\pm 7$ cm and $h_{\text{men}}=184\pm 8$ cm [heights].
Therefore, the mean slopes $\bar{m}_{s}$ of the slice score curves should be
as well approximately normal distributed. As discussed in section 4.3, the
mean slope $\bar{m}_{s}$ is estimated from the training data. With the help of
the mean slope, $\bar{m}_{s}$ basic data sanity checks can be applied. On the
one hand, it can be checked if the z-spacing is meaningful. On the other hand,
it can be checked if the z-axis order is reversed. For the validation of the
z-spacing a new variable, the slope ratio, $r_{m}$ and the expected relative
error $\hat{\gamma}$ is defined. The slope ratio $r_{m}$ is the fraction
between the observed slope $m_{X}$ for volume $X$ and the expected slope
$\bar{m}_{s}$:

$r_{m}=\frac{m_{X}}{\bar{m}_{s}}.$ (4.19)

The expected relative error is given by

$\hat{\gamma}=|1-|r_{m}||.$ (4.20)

The expected relative error $\hat{\gamma}$ is, therefore, the absolute value
of the relative error between the observed slope $m_{X}$ to the expected slope
$\bar{m}_{s}$. To decide if a slope ratio is valid, a threshold $\theta$ needs
to be defined. The threshold will be found data-driven based on the validation
set. The largest threshold $\theta$, for which for all volumes of the
validation set $\hat{\gamma}<\theta$ is chosen as threshold.

The expected z-spacing $\hat{z}$ can be calculated based on the slope of the
slice score curve and the expected slope $\bar{m}_{s}$. For a volume $X$, we
distinguish between two slopes. The slope $m^{\prime}_{X}$ depending on the
slice index and the slope $m_{X}$ depending on the height $h$ in mm. The
slopes are defined by fitting a linear line to the cleaned slice score curve.
Both slopes can be transformed into each other through

$m_{X}=\frac{1}{z_{X}}\cdot m^{\prime}_{X}.$

The expected slope $\bar{m}_{s}$ is declared in dependence of the slice height
$h$. Therefore, the expected z-spacing $\hat{z}_{X}$ for volume $X$ can be
calculated via the equation

$\hat{z}_{X}=\frac{m_{X}^{\prime}}{\bar{m}_{s}}.$ (4.21)

Figure 4.11 shows for an exemplary volume $X$ slice score curves of different
z-spacings. Moreover, the linear line for the expected z-spacing is visible.
It can be observed that the slope of the expected z-spacing and the observed
z-spacing are quite similar. The slopes of the slice score curves with twice
or half of the z-spacing can be visually detected as outliers. The slope of
the slice score curve is an indicator of the meaningfulness of the z-spacing.

Figure 4.11: Slice score slopes for an exemplary volume and different
z-spacings. The relation between z-spacing and the slope is visible. The
dotted line has the expected slope $\bar{m}_{s}$. Relative errors in the
z-spacing lead to a deviation in slope.

To apply both data sanity checks for an arbitrary CT-volumes $X$ during test-
time, we need to follow the steps:

  1. 1.

Predict slice scores for CT-volume $X$ with body part regression model

  2. 2.

Obtain cleaned slice scores (see fig. 4.8).

  3. 3.

Fit linear line to cleaned slice score curve.

  4. 4.

Check if slope is negative

$\begin{split}&\text{if: }m_{X}>0\Rightarrow\text{Natural z-ordering. Indices
increase with patient height.}\\\ &\text{if: }m_{X}<0\Rightarrow\text{Reverse
z-ordering. Indices decrease with patient height.}\\\ \end{split}$

  5. 5.

Check if the expected relative error $\hat{\gamma}_{X}$ is valid

$\begin{split}&\text{if: }\hat{\gamma}_{X}>\theta\Rightarrow\text{z-spacing is
invalid}\\\ &\text{if: }\hat{\gamma}_{X}<\theta\Rightarrow\text{z-spacing is
valid}\end{split}$

####  4.3.4 Deployment

To make the use cases UC1, UC2, and UC3 better applicable for the clinical
community and the research community, a public python package was deployed,
and the model was integrated into the medical platform toolkit Kaapana (see
sec. 2.3).

Main Functionality:  For the deployment, a main function for inference was
implemented. The function creates and saves a body part meta-data file for CT
images in the Neuroimaging Informatics Technology Initiative (NIfTI) data
format. The meta-data files are saved in the JavaScript Object Notation (JSON)
data format. Table 4.7 summarizes the information inside the JSON-file. The
threshold $\theta$ for the data sanity check can be customized as well as the
body part boundaries for calculating the body part examined dictionary and the
body part examined tag.

Table 4.7: Explanation of the tags in the JSON body part meta-data file. Tag | Format | Description
---|---|---
cleaned slice scores | list | cleaned slice scores from figure 4.8
unprocessed slice scores | list | transformed unprocessed slice scores
z | list | height in mm
body part examined | dictionary | body part examined dict from section 4.3.1
body part examined tag | string | body part examined tag from section 4.3.1
look-up table | dictionary | slice score reference table from section 4.3
reverse z-ordering | bool | data sanity check from section 4.3.3
valid z-spacing | bool | data sanity check from section 4.3.3
expected slope | float | expected slope $\bar{m}_{s}$ from section 4.3.3
observed slope | float | observed slope $m_{X}$ from section 4.3.3
slope ratio | float |  $r_{m}$ from equation 4.19
expected z-spacing | float |  $\hat{z}_{X}$ from equation 4.21
z-spacing | float | z-spacing of volume
settings | dictionary | used settings to generate the JSON-file

Public Model: To publish the body part regression model the model was trained
on a subset of the original training data from section 4.1.1. The data from
the whole body CT study and study Task051 was excluded based on the user
restriction for the data. The validation data and the test data stayed the
same. In total, 2091 volumes were used for training, 100 volumes for
validation and 100 volumes for testing. A subset of 42 volumes from the
training data was annotated.

Python Package: For the python package a README was created to describe how to
install and use the package. Moreover, jupyter notebook tutorials were
provided to explain how to train a body part regression model with the help of
the python package and how to use the package for inference of CT volumes in
different data formats.

Kaapana:  For Kaapana, the public model was used. The main function of
creating a JSON meta-data file for a NIfTI-volume was built into a Docker
container with the help of a Dockerfile. To integrate the Docker container
into the Kaapana toolkit, an operator was build upon the Docker container.
Additionally, a DAG was defined, which concatenates the following operators
(see sec. 2.3):

  1. 1.

get-input-data

  2. 2.

dcm-converter

  3. 3.

bodpartregression

  4. 4.

workflow-cleaner

All operators, apart from the bodypartregression operator, already existed in
Kaapana and are used to load the DICOM files, convert the DICOM files to
NIfTI-files, and to clean the temporary data of the workflow. The DAG can be
triggered through the Kaapana dashboard. The workflow uses the body parts
predicted from the DICOM volumes to extend the meta-data database of the
platform. This way, even more fine-grained exploration and filtering of data
can be achieved.

##  Chapter 5 Experiments and Results

In section 5.1, we investigate the effect of different loss functions and data
augmentation techniques on the model performance. Moreover, the best-
performing body part regression model is presented. In section 5.2, the
usefulness of the trained body part regression model is demonstrated through
several use cases. Additionally, the deployment of the body part regression
model is explained.

###  5.1 Body Part Regression

Multiple experiments were carried out to choose the best loss function (see
sec. 5.1.1) and to understand the model behavior better. All experiments were
conducted on the same dataset. The dataset consists of 2192 CT volumes for
training, 100 CT volumes for validation, and further 100 CT volumes for
testing (see sec. 4.1.1). Hyperparameters were optimized based on the model
performance on the validation set. Final statements about the generalization
to unseen data were made based on the model performance on the test set. After
running all experiments, the model with the best model performance on the
validation set was chosen and carefully evaluated.

During the experiments, the data augmentations and some parameters were fixed
to compare the experiments better and keep the search space small. The
selected fixed hyperparameter for the data augmentations can be found in table
4.4. For optimization, the Adam optimizer was used, which runs quite robustly
[kingma2014adam]. Additionally, a learning rate of $10^{-4}$ was used.
Moreover, the sampled slices per batch and the total sampled slices per CT
volume were fixed so that every model sees the same amount of two-dimensional
slices per batch and during the whole training procedure. Each mini-batch
contains about 256 axial slices. The sampled slices per batch are given by

$\text{sampled slices per batch}=\text{batch-size}\cdot m,$

where $m$ is the number of sampled slices per volume. From each CT volume
about 1920 slices are randomly sampled over the whole training period. The
total sampled slices per volume are given by

$\text{total sampled slices per volume}=m\cdot\text{epochs}.$

The random seed was fixed to zero. Further, the range from which $\Delta h$
was sampled, was fixed to [5 mm, 100 mm]. Table 5.1 summarizes the selected
hyperparameter settings. Training a body part regression model took between
1.5 and 2.5 hours on a Nvidia RTX 2080 ti GPU with 11 GB memory.

Table 5.1: Fixed parameters for the body part regression experiments.  Fixed parameter | Value
---|---
learning rate | $10^{-4}$
total sampled slices per volume | 1920
sampled slices per batch | 256
$\Delta h$-range | [5mm, 100mm]
random seed | 0

####  5.1.1 Loss Functions

In the following, the classification order loss (see sec. 4.1.4) and the
heuristic order loss (see sec. 4.1.4) will be compared. Further, we will
investigate the distance loss. To be able to compare both order losses, first
a good parameter $\beta$ for the newly proposed heuristic order loss needs to
be found. The optimal $\beta$ was searched on a logarithmic grid. Table 5.2
contains the search range, along with the other hyperparameters for this
experiment.

Table 5.2: Experiment setup for finding a good parameter $\beta$ for the heuristic order loss. Parameter | Description | Value[s]
---|---|---
m | sampled slices per volume | 4
batch-size | data-items in one mini-batch | 64
epochs | training epochs | 480
$\beta$ | parameter for heuristic order loss | 0.0001, 0.001, 0.01, 0.1, 1

As discussed in section 4.1.4, the parameter $\beta$ is an additional
constraint parameter, which influences how long the linearity condition
between $\Delta h$ and $\Delta s$ holds. Thus, a greater $\beta$ leads to a
weaker constraint. Figure 5.1 a) shows the Landmark Mean Square Error (LMSE)
curves during training for the different values of $\beta$. For
$\beta=1/\text{mm}$, the LMSE deteriorates, the constraint is too low.
Additionally, the absolute slice score predictions are drifting apart during
training. For $\beta=0.0001/\text{mm}$ the linear constraint is too high,
therefore the LMSE converges rather slowly. This would result in longer
training times. Moreover, the parameter $\beta$ is the proportionality factor
between $\Delta s$ and $\Delta h$ (see eq. 4.8). Therefore, we can observe in
figure 5.1 b) that the normalization constant $d$ increases with increasing
$\beta$.

Figure 5.1: Results for the heuristic order loss experiment. In figure a) the
influence of $\beta$ on the training convergence can be seen.The training
convergence is measured by the LMSE on the validation set. In figure b) the
procedure of the normalization constant
$d=\frac{1}{100}(\bar{s}_{l_{12}}-\bar{s}_{l_{1}})$ during training is
visible.

Table 5.3 shows a summary of the evaluation metrics for all trained models. It
can be seen that $\beta=0.01/\text{mm}$ performs best for the analyzed
metrics. This parameter is used for further experiments.

Table 5.3:  Influence of the hyperparameter $\beta$ on the model performance evaluated with the LMSE metric and the accuracy defined in section 4.2. $\beta$ | $\bar{\phi}_{\text{val}}\pm\sigma_{\bar{\phi}}$ | Accuracy in %
---|---|---
0.0001 | 6.3 $\pm$ 0.7 | 89.4
0.001 | 3.01 $\pm$ 0.27 | 93.9
0.01 | 2.39 $\pm$ 0.23 | 94.0
0.1 | 6.6 $\pm$ 1.3 | 91.7
1 | – | –

Next, we compare the heuristic order loss with the classification order loss.
The heuristic order loss and classification order loss are defined in equation
4.7 and equation 4.5.

Instability of the Classification Order Loss: During initial experiments, an
instability of the classification order loss was observed. The gradients and
the loss sometimes became "NaN". The reason for the instability is that the
classification order loss is not defined if $\sigma(\Delta s)=0$.
Mathematically, this can not happen, because the sigmoid function can not map
to zero. But numerically $\sigma(\Delta s)$ can be rounded to zero, if $\Delta
s\ll 0$. To stabilize the training, $\sigma(\Delta s)=0$ was ignored and the
updated order loss function is given by

$\begin{split}L_{\text{order}}^{c}&=-\frac{1}{N^{\star}}\sum_{i=1}^{B}\sum_{j=1,\text{if
}\sigma(\Delta s_{ij})\neq 0}^{m-1}\text{ln}\ \sigma(\Delta s_{ij}),\\\
N^{\star}&=\sum_{i=1}^{B}\sum_{j=1}^{m-1}[\sigma(\Delta s_{ij})\neq
0].\end{split}$ (5.1)

For comparing both loss functions, the following loss construction was used

$L=L_{\text{order}}+\alpha\cdot L_{\text{dist}},$ (5.2)

with the classification order loss (eq. 5.1) or the heuristic order loss (eq.
4.7) as order loss function and the distance loss function (eq. 3.5) as
additional constraint.

Experimental Setup:  To ensure a fair comparison of the methods, we will first
tune the hyperparameters. Table 5.4 describes the experimental setup. Based on
the fixed total sampled slices, slices per batch, and chosen $m$, the
corresponding batch size and epoch were calculated (see tab. 5.1). For the
classification order loss, Yan et al. used $\alpha=1$ for their analysis
[yan2018deep], therefore apart from the experiment with $\alpha=0$, the search
space was centered around one. For the sampled slices per data instance $m$,
Yan et al. chose $m=8$. Therefore, the search space for $m$ was centered
around 8. An additional experiment was carried out without any order loss and
only with the distance loss to understand the impact of the distance loss
better.

Table 5.4: Experimental setup for comparing the performances between the classification order loss and the heuristic order loss. The parameters $\alpha$, which controls the distance loss and the sampled slices per updated step $m$ were varied.  Order loss | Parameters | Values
---|---|---
heuristic order loss | (m, batch-size, epoch) | (4, 64, 480), (8, 32, 240), (12, 21, 160)
| $\alpha$ | 0, 0.01, 0.1, 1
| $\beta$ | 0.01
classification order loss | (m, batch-size, epoch) | (4, 64, 480), (8, 32, 240), (12, 21, 160)
| $\alpha$ | 0, 0.8, 1.0, 1.2
no order loss | (m, batch-size, epoch) | (8, 32, 240)
| $\alpha$ | 1

Experimental Results: In table 5.7 the LMSE and the accuracy for each trained
model is shown. The model which only uses the $L_{\text{dist}}$ as loss
function is by far the worst. Figure 5.2 shows the learning procedure of this
model. Based on figure 5.2 b), we can observe that the model predicts for all
landmarks nearly the same value. Moreover, in figure 5.2 c) we observe that
the LMSE has a high variance and does not decrease in the course of training,
although the validation loss $L_{\text{dist}}$ decreases. We can conclude that
the model performance is so bad because an easy minimum for $L_{\text{dist}}$
is to predict all slices to the same value. Without an additional order loss,
nothing prevents the model from converging to this state. The distance loss
$L_{\text{dist}}$ needs to be seen as an additional constraint, which
encourages a linear behavior between slice index and slice scores. However, as
an independent loss, it is not suitable for body part regression. It is
interesting to notice that the best model for the heuristic order loss
includes no additional distance loss $\alpha=0$. This supports the hypothesis
that the heuristic order loss has the linear constraint already incorporated,
and therefore the additional constraint through the distance loss is useless.

Figure 5.2: Training behavior of model with distance loss as only loss
function and with the hyperparameters $\alpha=1$ and $m=8$. The normalization
constant $d$ is given by $d=\frac{1}{100}(\bar{s}_{l_{12}}-\bar{s}_{l_{1}})$.

Based on the complete hyperparameter tuning results in table 5.7, the
following research questions are addressed: