figure table

# Generative Models for Reproducible Coronary Calcium Scoring

Sanne G.M. van Velzen  Department of Biomedical Engineering and Physics,
Amsterdam University Medical Centers - location AMC, the Netherlands  Image
Sciences Institute, Univerisity Medical Center Utrecht and Utrecht University,
Utrecht, the Netherlands  Bob D. de Vos  Department of Biomedical Engineering
and Physics, Amsterdam University Medical Centers - location AMC, the
Netherlands  Julia M.H. Noothout  Department of Biomedical Engineering and
Physics, Amsterdam University Medical Centers - location AMC, the Netherlands
Helena M. Verkooijen  Imaging Division, Univerisity Medical Center Utrecht,
Utrecht, the Netherlands  Max A. Viergever  Image Sciences Institute,
Univerisity Medical Center Utrecht and Utrecht University, Utrecht, the
Netherlands  Ivana Išgum  Department of Biomedical Engineering and Physics,
Amsterdam University Medical Centers - location AMC, the Netherlands
Department of Radiology and Nuclear Medicine, Amsterdam University Medical
Centers - location AMC, the Netherlands  Amsterdam Cardiovascular Sciences,
Amsterdam University Medical Centers, the Netherlands  Informatics Institite,
Faculty of Science, University of Amsterdam, the Netherlands

###### Abstract

Purpose: Coronary artery calcium (CAC) score, i.e. the amount of CAC
quantified in CT, is a strong and independent predictor of coronary heart
disease (CHD) events. However, CAC scoring suffers from limited interscan
reproducibility, which is mainly due to the clinical definition requiring
application of a fixed intensity level threshold for segmentation of
calcifications. This limitation is especially pronounced in non-ECG-
synchronized CT where lesions are more impacted by cardiac motion and partial
volume effects. Therefore, we propose a CAC quantification method that does
not require a threshold for segmentation of CAC. Approach: Our method utilizes
a generative adversarial network where a CT with CAC is decomposed into an
image without CAC and an image showing only CAC. The method, using a CycleGAN,
was trained using 626 low-dose chest CTs and 514 radiotherapy treatment
planning CTs. Interscan reproducibility was compared to clinical calcium
scoring in radiotherapy treatment planning CTs of 1,662 patients, each having
two scans. Results: A lower relative interscan difference in CAC mass was
achieved by the proposed method: 47% compared to 89% manual clinical calcium
scoring. The intraclass correlation coefficient of Agatston scores was 0.96
for the proposed method compared to 0.91 for automatic clinical calcium
scoring. Conclusions: The increased interscan reproducibility achieved by our
method may lead to increased reliability of CHD risk categorization and
improved accuracy of CHD event prediction.

###### keywords:

Calcium scoring, Reproducibility, CT, Generative models, CycleGAN

*Sanne G.M. van Velzen<EMAIL_ADDRESS>
##  1 Introduction

Coronary heart disease (CHD) is among the leading causes of death worldwide
[1]. A clear manifestation of CHD is calcification of the coronary arteries.
The amount of coronary artery calcifications (CAC), i.e. coronary calcium
score, is a strong and independent predictor of cardiac events, e.g.
myocardial infarction or heart failure [2, 3]. In CT, CAC is defined as a high
density area of $\geq$130 HU [4] in the artery. Using this definition CAC is
commonly quantified into CAC volume, mass or Agatston score[4]. The Agatston
score is used most because it directly indicates a patient’s risk of CHD
events.

In a clinical setting the amount of CAC is routinely quantified in non-
contrast-enhanced cardiac CT exams with ECG-synchronization. However, other CT
scans that show the heart, like chest CT, are increasingly used for
quantification of CAC and determination of CVD risk. Such scans are made about
ten times more often[5] than ECG-synchronized cardiac CT and moreover, the CAC
scores in these scans correlate well with the CAC scores in cardiac CT[6, 7].
Hence, given its clinical relevance, the guidelines of the Society of
Cardiovascular Computed Tomography and the Society of Thoracic Radiology [8]
recommend to quantify and report CAC on all non-contrast CTs showing the
heart, including non-ECG-synchronized scans. However, the interscan
reproducibility of calcium scores is limited. For example, the reported mean
interscan variability of Agatston scores in cardiac CT with ECG-
synchronization ranges from 15% to 41% [9, 10, 11, 12]. Hence, in non-ECG-
synchronized scans the limited interscan reproducibility is further emphasized
due to increased influence of partial volume effect and cardiac motion. For
instance, Jacobs et al.[13] reported an average interscan variability of 71%
for Agatston scores in non-triggered chest CT scans, that led to assignment of
different cardiovascular risk category in 24% of subjects. The compromised
insterscan reproducibility limits the applicability calcium quantification in
non-ECG-synchronized scans in clinic and for longitudinal studies, in which
reproducibility is of key importance to obtain reliable measurements.

Figure 1: An image containing coronary artery calcium (CAC) is decomposed into
an healthy tissue image without CAC and an image containing only CAC.

The clinical definition that requires the use of the aforementioned single
intensity level threshold at 130 HU has been identified as one of the major
causes of the limited reproducibility of CAC quantification[14], which may
lead to under- or overestimation of the amount of CAC. This is especially
pronounced in CT exams without ECG-synchronization where visualization of the
calcifications is often extremely affected by cardiac motion. As a
consequence, the calcifications are blurred or completely remain below the
threshold. Nevertheless, common to the vast majority of manual and automatic
calcium scoring methods is the use of the clinical definition that identifies
voxels in the coronary arteries above the 130 HU threshold as CAC in an
explicit manner (through segmentations) [15, 16] or implicit manner (through
regression of the CAC score)[17, 18]. To address this issue, automatic methods
have been proposed that either use automatic adaptive thresholds or omit
thresholding for segmentation of calcified lesions. Groen et al.[19] proposed
using an adaptive threshold dependent on the maximum intensity value of each
CAC lesion. Song et al. [20] proposed to adapt the threshold based on the
intensity of the background in the vicinity of lesions in combination with
deconvolution of the image, using a scanner-specific point spread function to
decrease partial volume effect. Saur et al. [21] omitted the intensity
threshold for segmentation by employing a mesh-based algorithm that segmented
CAC lesions by refining their boundary based on the intensity value profile. A
different method that circumvents thresholding for segmentation, was proposed
by Šprem et al. [22], who built on work by Dehmeshki et al. [23]. Both used an
expectation-maximization algorithm to determine the partial calcium content in
each voxel of a CAC lesion and its vicinity.

While the aforementioned methods all show improved interscan reproducibility
of CAC quantification, they use the clinically used CAC definition for initial
lesion detection, thereby missing lesions below the intensity level threshold.
This is particularly a problem with small lesions and scans with motion
artefacts, where complete lesions may remain under the detection threshold.
Although motion detection [24] in CAC lesions has been proposed, this does not
include correction that would allow quantification of such CAC lesions.

To address the limited interscan reproducibility of CAC detection and
quantification, we propose a method that does not depend on the current
clinically used definition of CAC that requires intensity-level thresholding.
Instead, we define CAC as the difference between an image visibly containing
CAC and a healthy tissue image. Our method takes an image with CAC and
decomposes it into an image without CAC and a CAC-map, i.e. an image showing
only CAC (Figure 1). The CAC-map comprises the location of lesions and the
grey values that indicate the attenuation in the CT scan due to CAC lesions.
Due to this new definition of CAC lesions, our method can detect lesions that
are completely - and partly - below the clinically used threshold.

Several earlier studies proposed to use an autoencoder for generating healthy
counterparts of images with pathology in an unsupervised manner[25, 26, 27].
Subsequently, the pixel-wise intensity difference between the original image
and the reconstruction was used to e.g. detect brain lesions in MRI scans
[25], brain CT [27], or detected anomalies in head CT [26]. Alternative
approaches used generative adversarial networks (GAN) for translating an image
containing pathologies to a healthy image. Baur et al. [28] and Seah et al.
[29] used conditional GANs to generate healthy images and used the difference
with the input image for multiple sclerosis lesion segmentation, and
visualizing features of congenital heart failure in chest X-rays,
respectively. Furthermore, Baumgartner et al. [30] proposed to use a
Wasserstein GAN to directly generate a difference map that indicates changes
in the brain related to Alzheimer’s disease. The map is added to the input
image to create a healthy brain image from an image containing Alzheimer’s
disease. Tang et al. [31] proposed to predict both an healthy image and a
difference map that indicates disease using a disentangled conditional GAN. To
tackle the instability and mode collapse problems that often arise during
training of GANs, Zhu et al. [32] proposed to use a cycle-consistency loss in
their CycleGAN for image-to-image translation using natural images. Sun et al.
[33] proposed to use a CycleGAN implementation for translating images
containing lesions to healthy images. Subsequently, the pixel-wise intensity
difference was used to detect brain lesions in MRI and liver tumors in CT
scans.

To decompose an image with CAC into an image without CAC and CAC-map, we
employ a CycleGAN to translate images between the containing CAC domain and
not containing CAC domain. To our knowledge, our method is the first to use a
generative model for accurate segmentation of CAC lesions. The method builds
on our preliminary work [34], in which we presented a semi-automatic method
for segmentation of CAC, that analyzed the vicinity of manually identified CAC
lesions. The method was trained using clinical calcium scores as reference and
hence, was not able to detect lesions below the threshold. In the current
work, we extend the method to perform fully automatic CAC detection and
quantification, including CAC lesions below the standard detection threshold.
The method is aimed at scans where the low interscan reproducibility is most
dire and therefore, experiments were performed in non-ECG-synchronized scans,
with and without breathing motion, that often contain excessive cardiac motion
artefacts and have a low spatial resolution. To simplify the challenging task
and allow the method to focus on the region of interest only, CAC
quantification was performed in the heart slices containing CAC. The method is
developed and evaluated using 4,038 CT scans of 2,276 breast cancer patients
and 626 CT scans of lung screening participants. Moreover, we show that the
here proposed method outperforms the gold standard for calcium scoring in the
clinic and a state-of-the-art automatic calcium scoring method in terms of
interscan reproducibility of CAC quantification and detection of visible
lesions.

##  2 Data

We confirm that all methods were carried out in accordance with relevant
guidelines and regulations. We included 4,038 radiotherapy treatment planning
(RTP) CTs of 2,276 breast cancer patients[35]. In this set 1,762 patients had
two scans made on the same day: one with and one without breath-hold (Figure
2, A and B). Therefore, this set allows evaluation of CAC scoring interscan
reproducibility. The remaining 514 patients had one CT scan made without
breath-hold. Scans were acquired in the University Medical Center Utrecht with
a Philips Brilliance Big Bore scanner. Intravenous (IV) contrast was not
induced and the acquisition was not ECG-synchronized. The CTs were acquired
with 120 kVp and reconstructed to 0.92-1.37 mm in-plane resolution and 3.0 mm
slice thickness and increment. This study was approved and the need for
informed consent was waived by the institutional review board of the
University Medical Center Utrecht, the Netherlands.

Additionally, we included 626 low-dose chest CT images of 626 subjects
included in the National Lung Screening Trail (NLST, Figure 2C). Scans were
acquired in 29 different hospitals, on scanners from all major vendors [15,
36]. Tube voltage was either set to 120 or 140 kVp, depending on the subjects
weight. IV contrast was not induced and acquisition was not ECG-synchronized.
Images were reconstructed to 0.49-0.98 mm in-plane resolution and 1.0-2.5 mm
slice thickness and 0.6-2.5 mm increment.

Figure 2: Examples of two RTP scans of one breast cancer patient, one with
breath-hold (A) and one without breath-hold (B), and an NLST scan of one
participant (C).

Table 1 lists the used data sets, availability of reference annotations and
describes the purpose of each data set. For development and training of the
method, RTP CT scans of 514 patients with a single CT ($RTP_{train}$ and
$RTP_{trainHS}$) and all 626 NLST scans ($NLST_{train}$) were used. The latter
were used to enhance training with a diverse set of scans of heavy smokers
with high CAC burden. In method development, a validation set of 100 patients
with two available RTP scans (200 scans) was used ($RTP_{val}$). We divided
the test scans into two sets with patient scan pairs: $RTP_{test1}$, a set
consisting of scans of 119 patients totalling 238 RTP CTs with manually
segmented CAC allowing detailed evaluation of calcium scoring; and
$RTP_{test2}$, a set consisting of scans of 1,543 patients of in total 3,086
RTP CTs without reference CAC scores, used for assessment of interscan
reproducibility.

Figure 3: Schematic example of the weak training labels per axial slice
indicating whether CAC is present in the slice (A). Indicated region of the
lesion with assigned coronary artery label (B).  Table 1: A. Data sets,
consisting of low-dose CT scans from the National Lung Screening Trial (NLST)
or radiotherapy treatment planning (RTP) CT scans, used for development
($NLST_{train},RTP_{train},RTP_{trainHS},RTP_{val}$) and evaluation
($RTP_{testHS},RTP_{test1},RTP_{testlesions},RTP_{test2}$) of the presented
method. Scans are made with (BH) or without deep inspiration breath-hold
(nonBH). B. Availability of manual reference annotations are indicated per set
for heart segmentations (Heart Segm.), labels per axial slice indicating
presence of CAC as a binary decision assessmed visually (Figure 3A) and
availability of coarse segmentations of lesion regions (Lesion Regions, Figure
3B). C. The way the data is used in the experiments is listed as training
(Tr), validation (Val) and testing (Test) for the different networks (Heart
Segmentation, CAC Classification, Calcium CycleGAN).

[b]

A | Scan information
---|---
| Scan Type | Patients | BH | nonBH
Development |  |  |
| $NLST_{train}$ | 626 | 626 |
| $RTP_{train}$ | 514 |  | 514
| $RTP_{trainHS}$ | 120 |  | 120
| $RTP_{val}$ | 100 | 100 | 100
Evaluation |  |  |
| $RTP_{testHS}$ | 79 |  | 79
| $RTP_{test1}$ | 119 | 119 | 119
| $RTP_{testlesions}$ | 34 | 34 | 34
| $RTP_{test2}$ | 1,543 | 1,543 | 1,543
B | Reference Annotations
|  |  | Heart
---
Segmentation
| Per Slice
---
CAC
| Lesion
---
Regions
Development |  |  |
| $NLST_{train}$ |  | ✓ |
| $RTP_{train}$ |  | ✓ |
| $RTP_{trainHS}$ | ✓ |  |
| $RTP_{val}$ |  | ✓ |
Evaluation |  |  |
| $RTP_{testHS}$ | ✓ |  |
| $RTP_{test1}$ |  | ✓ |
| $RTP_{testlesions}$ |  | ✓ | ✓
| $RTP_{test2}$ |  |  |
C | Purpose
|  |  | Heart
---
Segmentation
| CAC
---
Classification
| Calcium
---
CycleGAN
|  | Tr | Val | Test | Tr | Val | Test | Tr | Val | Test
Development |  |  |  |  |  |  |  |  |
| $NLST_{train}$ |  |  |  | ✓ |  |  | ✓ |  |
| $RTP_{train}$ |  |  |  | ✓ |  |  | ✓ |  |
| $RTP_{trainHS}$ | ✓ | ✓ |  |  |  |  |  |  |
| $RTP_{val}$ |  |  |  |  | ✓ |  |  | ✓ |
Evaluation |  |  |  |  |  |  |  |  |
| $RTP_{testHS}$ |  |  | ✓ |  |  |  |  |  |
| $RTP_{test1}$ |  |  |  |  |  | ✓ |  |  | ✓
| $RTP_{testlesions}$ |  |  |  |  |  | ✓ |  |  | ✓
| $RTP_{test2}$ |  |  |  |  |  | ✓ |  |  | ✓

To allow training, all training CTs, validation CTs and CTs from $RTP_{test1}$
were manually labelled without using a threshold by an observer with $>$3
years experience in calcium scoring. Because voxelswise annotation without a
threshold is practically infeasible, as the exact outline of the CAC is not
known, weak labels were used to train the methods. Therefore, the observer
labelled axial slices of the scans as containing CAC or not containing CAC,
based on visual assessment (Figure 3A). Moreover, for a subset of 199 RTP
scans of the training set ($RTP_{trainHS}$ and $RTP_{testHS}$), manual
segmentations of the heart were available from clinical practice.

To define reference standard for CAC, lesions were annotated and assigned
coronary artery label in $RTP_{test1}$. Because accurate voxel-level manual
annotation without threshold is hardly feasible, regions containing lesions
were indicated in the scans, by rough manual voxel-painting without a
threshold, with a brush size of 1 cm (Figure 3B). This was only done in the 34
scan pairs in which both scans had CAC, because pairs without CAC and visibly
discordant pairs are not informative for measuring interscan reproducibility.
Visibly discordant pairs are pairs in which in one scan CAC is visible and in
the second scan it is not. Next, to allow comparison with calcium scoring used
in clinic, manual annotations following the clinical definition from a
previous study for $RTP_{test1}$ were available [16].

Figure 4: Schematic overview of the proposed CycleGAN that translates the
images from the domain of images with CAC (CAC domain) to the domain of images
without CAC (No CAC domain) and back. In the pathway that removes CAC, a CAC
map is predicted by the generator (GR) and subtracted from the image to obtain
a synthetic image without CAC. In the pathway that synthesizes CAC the CAC map
predicted by the generator (GS) is added to the image. Synthetic images are
compared to real examples by discriminators (D).

##  3 Method

We propose an automatic method that detects and segments CAC by decomposing a
CT slice into an image without CAC and an image showing a CAC map, i.e. only
CAC. By defining CAC as the difference between an image visibly containing CAC
and a healthy tissue image, the method allows segmentation of whole lesions,
including parts below the current clinical definition applying 130 HU
threshold, leading to increased interscan reproducibility of CAC
quantification.

The need for improved interscan reproducibility of CAC quantification is
highest in non-ECG-synchronized scans. However, non-ECG-synchronized scans are
typically made for other purposes than calcium scoring, thus, the field of
view in these scans is highly variable. Therefore, before generating the CAC-
map, we define the region of interest to simplify the CAC detection task by
initial heart segmentation, and subsequent identification of CT slices with
CAC.

Table 2: Architecture of the Heart Segmentation CNN, CAC Classification CNN
and the Generators in the proposed CycleGAN

[b]  Heart Segmentation CNN CAC Classification CNN CAC-map generators CycleGAN
General Dimensions 3D 2D 2D Input size 64x64x64 voxels 224x224 voxels 224x224
voxels Type of padding (layers) Zero (All) Reflection (All) Reflection (All)
Architecture (kernel, output channels, stride) 1st Layer Conv (7x7x7, 16) Conv
(7x7, 64) Conv (7x7, 64) 2nd Layer Conv (3x3x3, 32, 2) Conv (3x3, 128, 2) Conv
(3x3, 128, 2) 3rd Layer Conv (3x3x3, 64, 2) Conv (3x3, 256, 2) Conv (3x3, 256,
2) ResNet Blocks 9 Blocks (3x3x3, 64, 1) 6 Blocks (3x3, 256, 1) 6 Blocks (3x3,
256, 1) 3rd last Layer TConv (3x3x3, 32, 2) - TConv (3x3, 128, 2) 2nd last
Layer TConv (3x3x3, 16, 2) Global Average Pool TConv (3x3, 64, 2) Last Layer
Conv (7x7x7, 2) Dense Layer (256x2 channels) Conv (7x7x7, 1) Activation last
Layer Sigmoid Softmax Sigmoid

  * •

Conv = Convolution layer; TConv = Transposed convolution layer; All
convolution layers, except the last layer, are followed by a batch
normalization layer and rectifier linear unit activation. Resnet Blocks
consist of 2 convolution layer pairs with a residual connection.

Heart segmentation is performed with a 3D convolutional neural network (CNN),
which, for efficiency, analyses the full image in a patch based manner. The
network architecture is based on ResNet [37] and is detailed in Table 2, left.
The Dice loss is used during training and Adam optimization algorithm
($lr=0.001$) is used for optimization. During inference, overlapping patches
were analysed and the results were averaged for final prediction.

We classify each slice in the heart mask according to the presence of visible
CAC with a 2D CNN with ResNet [37] based architecture (Table 2, middle). The
CNN analyses axial image slices centered and cropped around the heart. For
cropping, the center of mass of the heart segmentation is used. The network is
trained with per slice labels indicating the presence of CAC, using the cross-
entropy loss and Adam optimization algorithm ($lr=0.001$).

For the identified slices, we generate a CAC-map. Manually identifying voxels
with visible CAC to directly predicts a CAC-map with a CNN, would be hardly
feasible. Hence, we propose to obtain a CAC-map using weakly supervised
labels, i.e. labels indicating whether CAC is present in the image. In our
approach a CAC-map is obtained by subtracting a generated healthy tissue image
from an image containing CAC. Because it is not possible to obtain paired CT
scans of the same patient, where one would contain CAC and the other not, a
paired approach for image generation is not feasible. Therefore, an unpaired
approach is used to generate synthetic $\hat{noCAC}$ images from images
containing CAC. For this, images are assigned to two domains based on weak
labels: one domain consisting of images containing visible CAC and the other
domain containing images without visible CAC. Hence, to generate a healthy
tissue image for subtraction, a 2D CycleGAN[32] for unpaired image-to-image
translation is used to translate images from one domain to another. The
CycleGAN is a type of GAN [38] that improves in training stability and
performance over the traditional and conditional GAN, by adding cycle-
consistency [32].

The CycleGAN consists of two generators $G$, and two discriminators $D$ that
compete with each other. The generators are CNNs that convert an input image
from either domain to a realistic image from the other domain. In our task the
images are divided into the CAC domain that contains images with CAC and the
noCAC domain that contains images without CAC. The discriminator CNNs compare
the synthesized images against real examples from the target domain. In
contrast to a conventional CycleGAN, that directly translates an input image
to an image of the target domain, we use an adjusted CycleGAN (Figure 4) that
adds an extra step to minimize data hallucination: in the removing pathway,
images with CAC are translated to the noCAC domain by predicting a CAC map and
subtracting it from the input image: $\hat{noCAC}=CAC-G_{R}(CAC)$. Hence, the
CAC map $G_{R}(CAC)$ comprises both the location of CAC lesions and grey
values that indicate the attenuation in the CT scan due to CAC lesions. In the
synthesizing pathway the CAC map $G_{S}(noCAC)$ is added to the input image to
generate a synthetic image containing CAC from an image without CAC:
$\hat{CAC}=noCAC+G_{S}(noCAC)$. In this pathway the generator learns to
synthesize realistic lesions based on probable lesion locations and severity
inferred from the training data. In this work the synthesizing pathway is
solely utilized for training stability. The generator in the pathway that
removes CAC, i.e. $G_{R}$, is used to obtain the CAC map and, subsequently,
the CAC map can be used to quantify the amount of CAC in the image.

The generators used in the CycleGAN for generating the CAC map have a ResNet-
based architecture [37], with 6 ResNet blocks using ReLU activation and no
downsampling in between blocks (Table 2, right). To ensure the CAC maps only
contain positive values, the generators have a sigmoid activation in the last
layer. For the discriminator networks we use PatchGANs, which classify 70x70
overlapping image patches as real or synthetic, described in detail by [39].
The PatchGANs are able to capture high-frequency structure and have fewer
parameters than a full-image discriminator and are therefore, easier to train
[39]. The patch-level discriminators are applied to the image in a fully
convolutional manner, averaging all responses to provide the output of the
discriminator.

The system is trained with the loss terms on the input and synthetic images
that are traditionally used for training a CycleGAN[32]: Adversarial loss
(Equation 2), Cycle-consistency loss (Equation 3) and Identity loss (Equation
4). The Identity loss trains $G_{R}$ to predict an empty CAC-map when it is
presented with an image without CAC and vice versa for $G_{S}$. Additionally,
a Sparsity (L1) loss (Equation 5) on the CAC map is used to encourage the
generator to predict sparse CAC maps, giving the full objective:

$\textbf{L}=\textbf{L}_{Adversarial}+\lambda\textbf{L}_{Cycle}+\alpha\textbf{L}_{Identity}+\beta\textbf{L}_{Sparsity}$
(1)

Where

$\textbf{L}_{Adversarial}=\mathbb{E}_{CAC\sim
P_{data(CAC)}}[\log\textbf{D}_{CAC}(CAC)]+\mathbb{E}_{noCAC\sim
P_{data(noCAC)}}[\log(1-\textbf{D}_{CAC}(\hat{CAC}))]\\\
+\mathbb{E}_{noCAC\sim P_{data(noCAC)}}[\log\textbf{D}_{noCAC}(noCAC)]\\\
+\mathbb{E}_{CAC\sim P_{data(CAC)}}[\log(1-\textbf{D}_{noCAC}(\hat{noCAC}))]$
(2)

$\textbf{L}_{Cycle}=\mathbb{E}_{CAC\sim
P_{data(CAC)}}||\textbf{G\textsubscript{S}}(\textbf{G\textsubscript{R}}(CAC))-CAC||_{1}\\\
+\mathbb{E}_{noCAC\sim
P_{data(noCAC)}}||\textbf{G\textsubscript{R}}(\textbf{G\textsubscript{S}}(noCAC))-noCAC||_{1}$
(3)

$\textbf{L}_{Identity}=\mathbb{E}_{noCAC\sim
P_{data(noCAC)}}||\textbf{G\textsubscript{R}}(noCAC)||_{1}+\mathbb{E}_{CAC\sim
P_{data(CAC)}}||\textbf{G\textsubscript{S}}(CAC)||_{1}$ (4)

$\textbf{L}_{Sparsity}=\mathbb{E}_{CAC\sim
P_{data(CAC)}}||\textbf{G\textsubscript{R}}(CAC)||_{1}+\mathbb{E}_{noCAC\sim
P_{data(noCAC)}}||\textbf{G\textsubscript{S}}(noCAC)||_{1}$ (5)

To make the method more robust to noise, during training images are augmented
by reducing or amplifying noise. To reduce the noise, images are smoothed with
a Gaussian filter with $\sigma=0.5$ pixels. To amplify noise, a noise image is
obtained by subtracting the smoothed image from the original input image that
is, subsequently, added to the input image. Moreover, during training random
cropping and rotation are used for augmentation.

Despite the sparsity loss, the generated CAC-map may contain low levels of
noise, because of excessive image noise. Noise in the CAC map may lead to
unrealistically low values for soft tissue in the synthetic noCAC image.
Therefore, we mask voxels in the CAC-map that lead to a value of $<$ -10 HU in
the synthetic noCAC image, as those are likely noise voxels in the CAC-map.

Finally, CAC was quantified by CAC pseudo-mass and adjusted Agatston score.
The pseudo-mass was calculated by multiplying the sum of HU values of the
voxels in the lesions by the voxel volume. The adjusted Agatston score was
calculated by multiplying the area of a lesion with a density score (1: $<$199
HU, 2: 200–299 HU, 3: 300–399 HU, 4: $>$399 HU) determined by the maximum
density of the area, and summing over axial slices. Please note that the
adjusted Agatston score includes segmented (parts of) lesions below the
clinically used threshold. For comparison, the conventional Agatston score was
calculated based on segmentations obtained using the clinical definition,
which uses the following density scores: 0: $<$130 HU, 1: 130-199 HU, 2:
200–299 HU, 3: 300–399 HU, 4: $>$399 HU. Please note that the adjusted
Agatston score is equal to the conventional Agatston score when CAC lesions
are segmented using the 130 HU threshold according to the clinical definition
of CAC.

##  4 Evaluation

Detection of CAC lesions was evaluated using accuracy, sensitivity, rate of
false positive scans and F1-score, between the automatic and reference labels
that indicated visible CAC.

The reproducibility of the CAC score was evaluated in scan pairs, using
absolute relative differences of quantified CAC pseudo masses. To measure the
interscan agreement of adjusted Agatston scores, the two-way random intraclass
correlation coefficient (ICC) for absolute agreement was computed and Bland
Altman plots with 95% limits of agreement were examined[40]. Moreover,
patients were assigned one of four CHD risk categories (I: 0–10, II: 11–100,
III: 101–400, IV: $<$400) based on their adjusted Agatston score. Agreement
between risk categories of each patient determined from the patient’s two
scans was evaluated using Cohen’s linearly weighted kappa ($\kappa$).

##  5 Experiments and Results

We evaluate the detection performance of visible CAC and the interscan
reproducibility of the proposed method. Moreover, we compare the performance
of our method with two clinical calcium scoring methods that use the standard
clinically used CAC definition for segmentation: manual calcium scoring and
automatic calcium scoring with a validated, state-of-the-art deep learning
approach[16]. Finally, we perform an ablation study to evaluate the different
components of the method.

Prior to heart segmentation all scans were resampled to 1.5 mm in-plane
resolution for standardization across the data set and 3.0 mm slice spacing,
the slice spacing of manual reference annotations. The heart segmentation CNN
was trained using $RTP_{trainHS}$: a subset of 120 radiotherapy treatment
planning CT scans of $RTP_{train}$, for which reference heart segmentations
were available. Of this set, 100 scans were used for training and 20 scans for
validation. $RTP_{testHS}$ was used for evaluation (see Table 1). The CNN was
trained for 250,000 iterations with batches of 10 patches of 64x64x64 voxels.
Batches were balanced for presence of the heart during training. After
training, the heart segmentation CNN achieved a median Dice score of 0.95
(interquartile range 0.94-0.96), a median Haussdorf distance of 15 mm
(interquartile range 11 - 21 mm) and a median asymmetric surface distance of
1.6 mm (interquartile range 1.2 - 2.2 mm) over all test images in
$RTP_{testHS}$.

After the heart segmentation, images were resampled to 1.0 mm in-plane
resolution, which was the average resolution of the training set, and 1.5 mm
slice spacing, for comparison with clinical calcium scoring[16]. Moreover,
slices were clipped between -50 HU and 950 HU, for optimal contrast between
soft tissue and CAC, and scaled between 0 and 1. The training set of 626 low
dose chest CTs ($NLST_{train}$) and 315 radiotherapy treatment planning CTs
($RTP_{train}$) was used, with in total 65,675 2D slices containing the heart.
Given that the method is aimed for analysis of RTP CT, the validation set of
200 such scans of 100 breast cancer patients (100 pairs) was used
($RTP_{val}$, Table 1). During training manually defined reference labels
indicating the presence of CAC in axial CT slices were used. The slice
classification CNN was trained for 1,500,000 iterations using batches of 20
axial CT slices of 224x224 voxels. Batches were randomly sampled and balanced
for presence of visible CAC. After training, the slice classification CNN
achieved an per slice accuracy of 0.93, with a sensitivity of 0.78 for visible
CAC and a false positive rate of 0.06, on all 17,450 CT slices from
$RTP_{test1}$ containing the heart. On average 9.9% of the slices containing
the heart was analysed by the CycleGAN.

The adjusted CycleGAN for generating the CAC map was trained with the same
data set and preprocessing as the slice classification CNN. The adjusted
CycleGAN was trained for 375,000 iterations, after which the learning curves
converged, using batches of 4 axial CT slices of 224x224 voxels, which was the
maximum that fit in the available GPU memory. Batches were balanced for both
the presence of CAC and relative slice position in the heart, using the heart
segmentations. The training parameters were determined experimentally through
qualitative (visual) and quantitative assessment. Qualitative assessment
evaluated the appearance of the generated synthetic $\hat{noCAC}$ and
synthetic $\hat{CAC}$ images, and it led to setting the learning rate to
0.0001 and $\lambda$ to 10[32]. Thereafter, quantitative assessment was used
to define parameters $\alpha$ and $\beta$ through evaluation of the CAC
detection (F1-score) and CAC reproducibility (absolute relative difference) on
$RTP_{val}$. Figure 5 shows the evaluated parameters and the obtained results.
Based on these results, $\alpha=10$ and $\beta=0.001$ were used in further
experiments. Finally, the predicted CAC maps were resampled to the original
in-plane resolution prior to evaluation.

All experiments were performed using Python with Pytorch[41] on an Intel Xeon
E5-1620 3.60 GHz CPU with an NVIDIA Titan X GPU. The complete analysis takes
on average about 2 minutes per scan.

Figure 5: F1 score for per scan CAC detection and per scan absolute relative
difference for concordant positive scans on the validation set Figure 6: The
results of the ablation study for the method with and without heart
segmentation (HS -/+) and with and without slice classification (CL -/+). We
evaluated the detection performance (A) with sensitivity, false positive rate
(FPR) and F1-score. We evaluated the interscan reproducibility (B) with
absolute relative difference in CAC mass for all pairs ($\Delta_{R}$ All) and
concordant positive pairs ($\Delta_{R}$ Pos). For comparison, the performance
of clinical manual and clinical automatic calcium scoring are shown with
horizontal lines. Figure 7: Examples of lesions, visual assessment of the CAC
lesion area, calcium segmentation using the clinical definition, corresponding
synthetic images without CAC and segmentation of CAC with the proposed method
and following the clinical protocol. Lesions with excessive motion are
indicated with the blue box.

###  5.1 CAC detection per scan

In the clinic CAC is routinely detected on a scan-level and used for CVD risk
prediction. To provide insight in the detection performance of visible CAC on
a scan-level, Figure 7 illustrates the results and Figure 6A (proposed) shows
the sensitivity, rate of false positive scans (FPR) and F1-score in
$RTP_{test1}$ (Table 1). The proposed method achieved a sensitivity for
visible CAC of 0.83, with an FPR of 0.14, meaning 21 out of 150 scans without
visible CAC, and an F1-score of 0.80. In comparison, the sensitivity was 0.69
and 0.67 for manual and automatic clinical calcium scoring, respectively. The
F1-score of the clinical calcium scoring methods was comparable with our
method, due to a low FPR. Among our false positive cases, in 10 cases a lesion
in the aorta was erroneously detected (Figure 7, FP Aorta), among which 6 were
located in the ostia of the coronary arteries. In 9 scans noise was detected
as CAC by the proposed method. In 5 of these 9 scans, the false positive
lesion was located near the edge of the heart segmentation (Figure 7, FP
Noise).

As is shown in Figure 7, the proposed method also segments the parts of the
calcium lesions that remain below the threshold. Particularly in lesions
severely affected by motion, e.g. commonly lesions in the RCA, the difference
between the proposed method and manual clinical calcium scoring is especially
pronounced.

Table 3: Per artery detection performance (A) and interscan reproducibility
(B) of clinical manual (ClinicalM) and automatic (ClinicalA) calcium scoring,
and the proposed method in $RTP_{test1}$. CAC detection is shown for pairs and
for single scans. Reproducibility is evaluated with the absolute relative
difference in CAC pseudo-mass in scan-rescan artery pairs for concordant
positive artery pairs ($\Delta_{R}$Pos) as well as all artery pairs
($\Delta_{R}$All).

|  | (A) |  | (B)
---|---|---|---|---
|  | Detection |  | Interscan Reproducibility
|  | Single (%) | Pairs (%) |  |  $\Delta_{R}$ All |  $\Delta_{R}$ Pos$\dagger$
Proposed |  |  |  |
| CAC | 81 (79) | 34 (67) |  | 0.95 | 0.55
| LAD | 57 (92) | 26 (84) |  | 0.84 | 0.61
| RCA | 16 (57) | 5 (36) |  | 1.25 | 0.35
| LCX | 8 (67) | 3 (50) |  | 0.99 | 0.33
ClinicalM |  |  |  |
| CAC | 70 (69) | 30 (59) |  | 1.05 | 0.77
| LAD | 53 (85) | 23 (74) |  | 1.00 | 0.75
| RCA | 10 (36) | 4 (29) |  | 1.42 | 1.14
| LCX | 7 (58) | 3 (50) |  | 0.84 | 0.46
ClinicalA |  |  |  |
| CAC | 64 (63) | 26 (51) |  | 1.15 | 0.76
| LAD | 48 (77) | 21 (68) |  | 1.07 | 0.81
| RCA | 9 (32) | 2 (14) |  | 1.63 | 0.74
| LCX | 7 (58) | 3 (50) |  | 0.86 | 0.48

  * $\dagger$

Pairs in which CAC was detected in the artery in both scans by the respective
method.

###  5.2 CAC detection per artery

Because the impact of cardiac motion is different per location, we
additionally evaluate the detection performance per artery. A subset of 34
concordant scan pairs from $RTP_{test1}$ that contain CAC with per artery CAC
annotations ($RTP_{testlesions}$) were used for evaluation. Overall, the
proposed method detected CAC in 81 out of 102 arteries with lesions (Table 3A,
left). Out of 51 lesion pairs, the method detected both lesions in 34 pairs
(Table 3A, right). In comparison, manual clinical calcium scoring detected 70
out of 102 single lesions, and automatic clinical calcium scoring detected 64
out of 102. Manual clinical calcium scoring detected CAC in both lesions in 30
out of 51 lesion pairs and automatic clinical calcium scoring detected 26
pairs out of 51. The achieved improvement of our method over the clinical
methods was largest in the RCA, in which 64% of the visible lesions remained
below the 130 HU threshold due to blurring caused by cardiac motion. This is
also illustrated in Figure 7, where all lesions are affected by partial volume
effect due to low resolution and especially lesions in the RCA are subject to
excessive motion artefacts.

Table 4: Reproducibility of CAC pseudo mass and Agatston scores for automatic
clinical calcium scoring and the proposed method in $RTP_{test2}$

.  All pairs Concordant positive pairs CAC mass AG score CAC mass AG score
1543 pairs $\Delta_{R}$ ICC ICC n (%) $\Delta_{R}$ ICC ICC Proposed 0.88 0.97
(0.96-0.97) 0.96 (0.96-0.96) 515 (34%) 0.44 0.96 (0.96-0.97) 0.95(0.95-0.96)
Clinical 0.92 0.93 (0.92-0.94) 0.91 (0.91-0.92) 375 (24%) 0.58 0.92
(0.90-0.93) 0.90 (0.88-0.92)

###  5.3 Interscan reproducibility of CAC quantification

To make CAC quantification in non-ECG-synchronized CT scans useful for
clinical application the interscan reproducibility ideally should be high.
Therefore, the interscan reproducibility of CAC pseudo mass was evaluated in
$RTP_{test1}$. First, we evaluated the absolute relative difference in CAC
pseudo mass in 51 artery pairs using per artery annotations of lesions
described in Section 2. The interscan difference was 95% on average in all
artery pairs (Table 3B, left). In concordant positive artery pairs, i.e. pairs
in which the method detected CAC in both arteries, the interscan difference
was 55% using the proposed method (Table 3B, right). Second, we evaluated the
reproducibility of total CAC pseudo mass in scan pairs. The interscan
difference was 94% on average in all scan pairs (Figure 6B, left). For
concordant positive scan pairs the absolute relative difference was 47%
(Figure 6B, right). The higher absolute relative difference in CAC pseudo mass
(i.e. lower reproducibility) in all scan pairs compared to positive pairs, is
partly explained by the presence of 21 discordant pairs. In these pairs the
method shows a high absolute relative difference, since a CAC lesion is only
visible in one of the two scans. Additionally, we compared our results to
clinical calcium scoring. Both clinical methods show a lower interscan
reproducibility than our proposed method in artery pairs as well as in scan
pairs. For manual clinical calcium scoring the interscan difference was 77% in
concordant positive artery pairs and 89% in concordant positive scan pairs.
For automatic clinical calcium scoring the interscan difference was 76% and
114% in artery pairs and scan pairs, respectively. For similar reasons as for
our method, the relative interscan difference is high when all scan pairs are
considered.

Figure 8: Bland-Altman plots for adjusted Agatston scores of the proposed and
automatic clinical method. For comparison purposes the adjusted Agatston
scores of the proposed method are scaled to the range of scores of the
clinical method. Please note that, since the same linear scaling is used for
all scores, this does not influence the agreement. 95% limits of agreement are
indicated by dashed lines. Because the errors tend to increase with increasing
CAC, regression for nonuniform differences was used to model the variation of
the absolute differences in scan pairs. To calculate the 95% limits of
agreement, the predicted absolute differences were multiplied by
$1.96\times(\pi/2)^{1/2}$, because the absolute differences have a half-normal
distribution [40].

Moreover, the reproducibility of the proposed method is evaluated in
$RTP_{test2}$, in which manual annotations are not available. Hence, the
results of our method are compared with automatic clinical calcium scoring.
The results are listed in Table 4. The relative difference in CAC mass for
positive pairs with detected CAC in both scans was 44%. For all scan pairs
this was 88%. Moreover, the interscan agreement of CAC mass was high with an
ICC of 0.96 (95% CI 0.96-0.97) in positive pairs and 0.97 (95% CI 0.96-0.97)
in all pairs. In comparison, the relative difference was 58% in concordant
positive pairs and 92% in all pairs when using the automatic clinical method.
Moreover, the ICC for the quantified CAC mass was significantly lower for the
clinical automatic method than for the proposed method (Table 4).

Table 5: Interscan agreement of risk categorization based on (adjusted)
Agatston score in $RTP_{test1}$ and $RTP_{test2}$. The risk categories are I:
0-10, II: 11-100, III: 101-400, IV: $\geq$400.

[b]  Proposed$\dagger$ ClinicalM ClinicalA $RTP_{test1}$ $\kappa$ 0.82
(0.71-0.94) 0.63 (0.42-0.84) 0.64 (0.41-0.87) Accuracy 0.92 0.92 0.95
$RTP_{test2}$ All pairs $\kappa$ 0.84 (0.82-0.87) - 0.81 (0.78-0.84) Accuracy
0.90 - 0.92 Concordant positive pairs$\ddagger$ $\kappa$ 0.77 (0.72-0.82) -
0.67 (0.61-0.72) Accuracy 0.79 - 0.68

  * $\dagger$

For comparison purposes, the adjusted Agatston scores obtained with the
proposed method were linearly scaled to the range of the scores obtained with
automatic clinical (ClinicalA) calcium scoring, using the mean score of both
methods.

  * $\ddagger$

Pairs in which the ClinicalA method found CAC in both scans.

The proposed method achieved an interscan agreement for adjusted Agatston
scores with an ICC of 0.96 (95% CI 0.96-0.96). This is significantly better
than for the automatic clinical method that achieved an ICC of 0.91 (95% CI
0.91-0.92). This improvement is clearly visible in the Bland-Altman plots
(Figure 8), where the 95% limits of agreement are narrower for the proposed
method than for automatic clinical calcium scoring.

###  5.4 Ablation Study

To investigate the benefits of prior heart segmentation and slice
classification, we evaluated performance using the following settings: 1)
using the settings as proposed, 2) the proposed method without slice
classification, where all CT slices of the heart were analysed, 3) the
proposed method without heart segmentation, 4) the proposed method without
heart segmentation and without slice classification. When heart segmentation
was not used, we standardized the field of view over the data sets to include
the heart and chest wall. CAC detection, FPR and absolute relative interscan
difference of the quantified CAC mass was evaluated on $RTP_{test1}$. Figure 6
shows the results. The proposed method achieved the best performance, with an
F1-score of 0.80 and an absolute relative difference in CAC mass of 47%. While
both approaches without slice classification show a high sensitivity, the
false positive rate is high: with settings 2 the method found false positive
lesions in all negative scans and with settings 4 in 57% of the negative
scans. The effect of defining the region of interest by heart segmentation is
not only reflected in a better detection, but also in a higher interscan
reproducibility of CAC mass, i.e. lower absolute relative difference.

###  5.5 Risk category assignment

Because CVD risk stratification based on the Agatston score is clinically
relevant, we evaluate the interscan reproducibility of risk categorization.
The interscan agreement of CVD risk categories is shown in Table 5. The
agreement in $RTP_{test1}$ was higher using scores derived with the proposed
($\kappa$ = 0.82) method than with manual of automatic clinical calcium
scoring, with $\kappa$ of 0.63 and 0.64 respectively. In concordant positive
pairs of $RTP_{test2}$ the $\kappa$ was 0.77 for the proposed method. In
contrast, for automatic clinical calcium scoring this was 0.67.

##  6 Discussion

We presented a method for calcium scoring that does not depend on the standard
clinically used definition of CAC that applies an intensity value threshold of
130 HU for segmentation of CAC lesions. Instead, our method defines CAC as the
difference between an image visibly containing CAC and a healthy tissue image.
Hence, the method can identify CAC - complete or parts of lesions - that may
remain undetected by the clinical threshold due to e.g. partial volume effect
in small lesions or blurring due to cardiac motion. We achieve this by
separating an image containing CAC into an image without CAC and an image only
containing CAC. For this we exploit a CycleGAN that translates images between
the CAC and noCAC domain. We have shown that segmentation of lesions using our
method enables increased interscan reproducibility of CAC quantification,
compared to clinically used manual and automatic CAC scoring. Increased
interscan reproducibility may lead to more reliable risk estimation and enable
longitudinal studies in non-ECG-synchronized CT.

Previous research shows that interscan reproducibility in non-ECG synchronized
CT scans is lower than in dedicated cardiac CT [9, 10, 11, 12]. Because of
limitations in the acquisition of radiotherapy treatment planning CTs, our
method does not exceed the in literature reported reproducibility in dedicated
cardiac CT. However, CAC scoring is increasingly performed in non-dedicated CT
scans, hence, the need for more reproducible CAC quantification is growing,
especially for quantification in challenging scans like the ones used in this
work. Our proposed method shows potential for increasing the interscan
reproducibility in non-ECG synchronized scans to the level of dedicated
cardiac CT.

In contrast to clinical methods, the proposed method is able to detect visible
CAC lesions that remain below the threshold. Owing to a more accurate
segmentation that includes pairs of lesions below the threshold, our method
outperforms clinical CAC scoring methods in interscan reproducibility of CAC
quantification. The results of per artery evaluation show that this
improvement was especially pronounced for CAC in the RCA. This is probably due
to the ability of the proposed method to detect parts of lesions that are
heavily affected by motion artefacts, which occurs more for the lesions in the
RCA than in the LAD and LCX.

Errors were made in a number of scans. False positive detections mostly
consisted of noise in the proximity of an artery or other types of
calcifications than CAC, such as calcifications in the aorta or cardiac
valves. Despite noise augmentation during training and post processing, the
CAC map contained noise in a few cases. Such errors were often located close
to the apex and near the RCA, where the noise level was most severe. Noise
reduction strategies described in e.g. Wolterink et al. [42] could offer a
solution. False positives in the aorta or valves were often close to the ostia
of the coronaries or occurred in the mitral valve. Since in scans without ECG
triggering it is often difficult to distinguish CAC from calcifications in the
thoracic aorta or mitral valve calcifications, these type of errors are also
not uncommon for human observers [43] and are also present in other automatic
methods [16, 15, 17]. False negative lesions were typically small with very
low voxel intensities, making them difficult to distinguish from noise and
soft tissue, also for experts.

Similarly to our preliminary work [34], the method without heart segmentation
and slice classification suffered from a large amount of false positives in
the CAC map, which made fully automatic quantification of CAC infeasible. In
our ablation experiment, we showed that adding slice classification with
respect to CAC presence and heart segmentation solves this issue. Using
classification to identify the slices containing CAC, prevents false positive
errors in slices that do not contain CAC. By adding heart segmentation, we
prevent anatomical structures outside the heart in the data from influencing
the CAC segmentation, resulting in more reproducible quantification. Likewise,
in previous methods for CAC quantification in scans with high noise levels, a
false positive reduction step was used to decrease the number of false
positive findings[16, 15]. Possibly, if a very large amount of training data
would be available, the method would be able to learn to avoid these errors.
However, in the medical field this is often challenging to obtain. In this
work, heart segmentation and slice classification are performed independently.
Although they are relatively straightforward, future work will investigate
merging these steps into a single network.

The heart segmentation network showed good performance in terms of Dice and
asymmetric surface distance. However, a relatively high Hausdorff distance was
found. Visual evaluation revealed that errors causing high Hausdorff distance
were mostly located in the basal slices of the heart. This is probably due to
high interobserver variability in reference segmentations regarding the
difficult to outline superior boundaries of the heart. Given that the
segmented superior boundaries are typically well above the ostia of the
coronaries, we expect the influence of variability in the superior boundaries
of the heart on CAC classification and quantification to be minimal.

During experiments the impact of the different loss terms for training the
Calcium CycleGAN on the CAC detection and reproducibility was evaluated using
relatively large differences in the values of the tested parameters. Possibly
a performance gain can be achieved by further optimizing $\alpha$ and $\beta$.
Nevertheless, given the obtained results (Figure 5), we expect the benefit
would be limited.

Other approaches for increasing the reproducibility of CAC quantification have
previously been proposed[19, 23, 22, 20, 21]. Since the true amount of CAC is
not available for patients, often a phantom is used for development and
evaluation of the method. This makes direct comparison of our method
infeasible. Šprem et al.[22] additionally evaluated their approach in a set of
293 subject CT scans and showed an improved reproducibility of CAC volume
between different CT images. However, different reconstructions of the same
acquisition were used, excluding interscan differences, in contrast to two
separate acquisitions used in our study. Moreover, in contrast to our approach
a major drawback of previously proposed methods is the use of the standard 130
HU intensity threshold for detection of CAC prior to partial volume
correction, inevitably missing lesions below the threshold.

Risk categorization of subjects is typically done using total CAC Agatston
scores. However, the proposed method also includes parts of lesions below the
clinical threshold, and thus, the Agatston score tends to be (substantially)
higher than for the clinical methods. Therefore, assigning risk categories
based on the Agatston score derived with the proposed method would necessarily
lead to different definitions of CVD risk categories. In order to estimate the
performance of risk categorization of the proposed method, the obtained
Agatston scores were scaled to the range of clinical Agatston scores. The
results showed an improved agreement for the proposed method, compared to
automatic clinical calcium scoring. However, to establish the clinical
relevance, further research relating the adjusted Agatston scores with CHD
events is warranted.

Future research on evaluation of the proposed method in non-ECG-synchronized
against calcium scoring in dedicated calcium scoring CT is needed. However,
clinical calcium scoring in dedicated cardiac CT also suffers from low
reproducibility[9, 10, 11, 12]. A non-invasive imaging-based reference that
indicates the true amount of CAC does not exist. Hence, to define the clinical
relevance, future research investigating whether more reproducible CAC
quantification with the proposed method leads to better prediction of CVD
events is warranted. However, since CVD related outcomes were not available,
investigating the relation to CVD events is out of scope of this work.

##  7 Conclusion

We have proposed a method for detection and quantification of CAC in non ECG-
synchronized chest CT scans that that does not require the application of the
clinically used intensity level threshold. Hence, the method allows detection
of visible CAC that may not exceed the standardly used 130 HU threshold due to
e.g. partial volume effect or blurring due to cardiac motion. This leads to
more reproducible quantification of CAC and increased interscan
reproducibility of risk categorization compared to the current clinical
method.

## Disclosures

The authors declare no competing interests.

## Acknowledgment

The authors gratefully acknowledge the Dutch Cancer Society for the financial
support (NCT03206333). The authors thank the National Cancer Institute for
access to NCI’s data collected by the National Lung Screening Trial. The
statements contained herein are solely those of the authors and do not
represent or imply concurrence or endorsement by NCI.

## References

  * [1] Wang, H., Naghavi, M., Allen, C., Barber, R. M., Bhutta, Z. A., Carter, A., Casey, D. C., Charlson, F. J., Chen, A. Z., Coates, M. M., et al., “Global, regional, and national life expectancy, all-cause mortality, and cause-specific mortality for 249 causes of death, 1980–2015: a systematic analysis for the global burden of disease study 2015,” The Lancet 388(10053), 1459–1544 (2016).
  * [2] Hecht, H. S., “Coronary artery calcium scanning: past, present, and future,” JACC: Cardiovascular Imaging 8(5), 579–596 (2015).
  * [3] Budoff, M. J., Nasir, K., McClelland, R. L., Detrano, R., Wong, N., Blumenthal, R. S., Kondos, G., and Kronmal, R. A., “Coronary calcium predicts events better with absolute calcium scores than age-sex-race/ethnicity percentiles: MESA (Multi-Ethnic Study of Atherosclerosis),” Journal of the American College of Cardiology 53(4), 345–352 (2009).
  * [4] Agatston, A. S., Janowitz, W. R., Hildner, F. J., Zusmer, N. R., Viamonte, M., and Detrano, R., “Quantification of coronary artery calcium using ultrafast computed tomography,” Journal of the American College of Cardiology 15(4), 827–832 (1990).
  * [5] De González, A. B., Mahesh, M., Kim, K.-P., Bhargavan, M., Lewis, R., Mettler, F., and Land, C., “Projected cancer risks from computed tomographic scans performed in the united states in 2007,” Archives of internal medicine 169(22), 2071–2077 (2009).
  * [6] Wu, M.-T., Yang, P., Huang, Y.-L., Chen, J.-S., Chuo, C.-C., Yeh, C., and Chang, R.-S., “Coronary arterial calcification on low-dose ungated mdct for lung cancer screening: concordance study with dedicated cardiac ct,” American Journal of Roentgenology 190(4), 923–928 (2008).
  * [7] Budoff, M. J., Nasir, K., Kinney, G. L., Hokanson, J. E., Barr, R. G., Steiner, R., Nath, H., Lopez-Garcia, C., Black-Shinn, J., and Casaburi, R., “Coronary artery and thoracic calcium on noncontrast thoracic ct scans: comparison of ungated and gated examinations in patients from the copd gene cohort,” Journal of cardiovascular computed tomography 5(2), 113–118 (2011).
  * [8] Hecht, H. S., Cronin, P., Blaha, M. J., Budoff, M. J., Kazerooni, E. A., Narula, J., Yankelevitz, D., and Abbara, S., “2016 SCCT/STR guidelines for coronary artery calcium scoring of noncontrast noncardiac chest CT scans: A report of the society of cardiovascular computed tomography and society of thoracic radiology,” Journal of Cardiovascular Computed Tomography 11(1), 74–84 (2017).
  * [9] Detrano, R. C., Anderson, M., Nelson, J., Wong, N. D., Carr, J. J., McNitt-Gray, M., and Bild, D. E., “Coronary calcium measurements: effect of CT scanner type and calcium measure on rescan reproducibility: MESA study,” Radiology 236(2), 477–484 (2005).
  * [10] Mao, S., Bakhsheshi, H., Lu, B., Liu, S. C., Oudiz, R. J., and Budoff, M. J., “Effect of electrocardiogram triggering on reproducibility of coronary artery calcium scoring,” Radiology 220(3), 707–711 (2001).
  * [11] Hoffmann, U., Siebert, U., Bull-Stewart, A., Achenbach, S., Ferencik, M., Moselewski, F., Brady, T. J., Massaro, J. M., and O’Donnell, C. J., “Evidence for lower variability of coronary artery calcium mineral mass measurements by multi-detector computed tomography in a community-based cohort—consequences for progression studies,” European Journal of Radiology 57(3), 396–402 (2006).
  * [12] Van Hoe, L. R., De Meerleer, K. G., Leyman, P. P., and Vanhoenacker, P. K., “Coronary artery calcium scoring using ecg-gated multidetector CT: effect of individually optimized image-reconstruction windows on image quality and measurement reproducibility,” American Journal of Roentgenology 181(4), 1093–1100 (2003).
  * [13] Jacobs, P. C., Isgum, I., Gondrie, M. J., Mali, W. P. T. M., van Ginneken, B., Prokop, M., and van der Graaf, Y., “Coronary artery calcification scoring in low-dose ungated CT screening for lung cancer: interscan agreement,” American Journal of Roentgenology 194(5), 1244–1249 (2010).
  * [14] Hong, C., Bae, K. T., and Pilgram, T. K., “Coronary artery calcium: accuracy and reproducibility of measurements with multi–detector row CT—assessment of effects of different thresholds and quantification methods,” Radiology 227(3), 795–801 (2003).
  * [15] Lessmann, N., van Ginneken, B., Zreik, M., de Jong, P. A., de Vos, B. D., Viergever, M. A., and Išgum, I., “Automatic calcium scoring in low-dose chest CT using deep neural networks with dilated convolutions,” IEEE Transactions on Medical Imaging 37(2), 615–625 (2017).
  * [16] van Velzen, S. G., Lessmann, N., Velthuis, B. K., Bank, I. E., van den Bongard, D. H., Leiner, T., de Jong, P. A., Veldhuis, W. B., Correa, A., Terry, J. G., et al., “Deep learning for automatic calcium scoring in CT: validation using multiple cardiac CT and chest CT protocols,” Radiology 295(1), 66–79 (2020).
  * [17] de Vos, B. D., Wolterink, J. M., Leiner, T., de Jong, P. A., Lessmann, N., and Išgum, I., “Direct automatic coronary calcium scoring in cardiac and chest CT,” IEEE Transactions on Medical Imaging (2019).
  * [18] Cano-Espinosa, C., González, G., Washko, G. R., Cazorla, M., and Estépar, R. S. J., “Automated agatston score computation in non-ecg gated CT scans using deep learning,” in [Medical Imaging 2018: Image Processing ], 10574, 105742K, International Society for Optics and Photonics (2018).
  * [19] Groen, J., Dijkstra, H., Greuter, M., and Oudkerk, M., “Threshold adjusted calcium scoring using CT is less susceptible to cardiac motion and more accurate,” Medical Physics 36(2), 438–446 (2009).
  * [20] Song, Y., Eck, B. L., Levi, J., and Wilson, D. L., “Improved reproducibility of calcium mass score using deconvolution and partial volume correction,” in [Medical Imaging 2019: Biomedical Applications in Molecular, Structural, and Functional Imaging ], 10953, 109531O, International Society for Optics and Photonics (2019).
  * [21] Saur, S. C., Alkadhi, H., Desbiolles, L., Székely, G., and Cattin, P. C., “ACCURATUM: improved calcium volume scoring using a mesh-based algorithm—a phantom study,” European Radiology 19(3), 591–598 (2009).
  * [22] Šprem, J., De Vos, B. D., Lessmann, N., Van Hamersvelt, R. W., Greuter, M. J., De Jong, P. A., Leiner, T., Viergever, M. A., and Išgum, I., “Coronary calcium scoring with partial volume correction in anthropomorphic thorax phantom and screening chest CT images,” PloS One 13(12), e0209318 (2018).
  * [23] Dehmeshki, J., Ye, X., Amin, H., Abaei, M., Lin, X., and Qanadli, S. D., “Volumetric quantification of atherosclerotic plaque in CT considering partial volume effect,” IEEE Transactions on Medical Imaging 26(3), 273–282 (2007).
  * [24] Šprem, J., de Vos, B. D., Lessmann, N., de Jong, P. A., Viergever, M. A., and Išgum, I., “Impact of automatically detected motion artifacts on coronary calcium scoring in chest computed tomography,” Journal of Medical Imaging 5(4), 044007 (2018).
  * [25] Chen, X. and Konukoglu, E., “Unsupervised detection of lesions in brain mri using constrained adversarial auto-encoders,” in [International Conference on Medical Imaging with Deep Learning ], (2018).
  * [26] Sato, D., Hanaoka, S., Nomura, Y., Takenaga, T., Miki, S., Yoshikawa, T., Hayashi, N., and Abe, O., “A primitive study on unsupervised anomaly detection with an autoencoder in emergency head CT volumes,” in [Medical Imaging 2018: Computer-Aided Diagnosis ], 10575, 105751P, International Society for Optics and Photonics (2018).
  * [27] Pawlowski, N., Lee, M. C., Rajchl, M., McDonagh, S., Ferrante, E., Kamnitsas, K., Cooke, S., Stevenson, S., Khetani, A., Newman, T., et al., “Unsupervised lesion detection in brain CT using bayesian convolutional autoencoders,” in [International Conference on Medical Imaging with Deep Learning ], (2018).
  * [28] Baur, C., Wiestler, B., Albarqouni, S., and Navab, N., “Deep autoencoding models for unsupervised anomaly segmentation in brain MR images,” in [International MICCAI Brainlesion Workshop ], 161–169, Springer (2018).
  * [29] Seah, J. C., Tang, J. S., Kitchen, A., Gaillard, F., and Dixon, A. F., “Chest radiographs in congestive heart failure: visualizing neural network learning,” Radiology 290(2), 514–522 (2019).
  * [30] Baumgartner, C. F., Koch, L. M., Can Tezcan, K., Xi Ang, J., and Konukoglu, E., “Visual feature attribution using wasserstein gans,” in [Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition ], 8309–8319 (2018).
  * [31] Tang, Y., Tang, Y., Zhu, Y., Xiao, J., and Summers, R. M., “A disentangled generative model for disease decomposition in chest x-rays via normal image synthesis,” Medical Image Analysis 67, 101839 (2021).
  * [32] Zhu, J.-Y., Park, T., Isola, P., and Efros, A. A., “Unpaired image-to-image translation using cycle-consistent adversarial networks,” in [Proceedings of the IEEE International Conference on Computer Vision ], 2223–2232 (2017).
  * [33] Sun, L., Wang, J., Huang, Y., Ding, X., Greenspan, H., and Paisley, J., “An adversarial learning approach to medical image synthesis for lesion detection,” IEEE Journal of Biomedical and Health Informatics (2020).
  * [34] van Velzen, S. G., de Vos, B. D., Verkooijen, H. M., Leiner, T., Viergever, M. A., and Išgum, I., “Coronary artery calcium scoring: Can we do better?,” in [Medical Imaging 2020: Image Processing ], 11313, 113130G, International Society for Optics and Photonics (2020).
  * [35] Emaus, M. J., Išgum, I., van Velzen, S. G., van den Bongard, H. D., Gernaat, S. A., Lessmann, N., Sattler, M. G., Teske, A. J., Penninkhof, J., Meijer, H., et al., “Bragatston study protocol: a multicentre cohort study on automated quantification of cardiovascular calcifications on radiotherapy planning CT scans for cardiovascular risk prediction in patients with breast cancer,” BMJ Open 9(7), e028752 (2019).
  * [36] National Lung Screening Trial Research Team, “Reduced lung-cancer mortality with low-dose computed tomographic screening,” New England Journal of Medicine 365(5), 395–409 (2011).
  * [37] He, K., Zhang, X., Ren, S., and Sun, J., “Deep residual learning for image recognition,” in [Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition ], 770–778 (2016).
  * [38] Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y., “Generative adversarial nets,” in [Advances in Neural Information Processing Systems ], 2672–2680 (2014).
  * [39] Isola, P., Zhu, J.-Y., Zhou, T., and Efros, A. A., “Image-to-image translation with conditional adversarial networks,” in [Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition ], 1125–1134 (2017).
  * [40] Sevrukov, A. B., Bland, J. M., and Kondos, G. T., “Serial electron beam CT measurements of coronary artery calcium: has your patient’s calcium score actually changed?,” American Journal of Roentgenology 185(6), 1546–1553 (2005).
  * [41] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S., “Pytorch: An imperative style, high-performance deep learning library,” in [Advances in Neural Information Processing Systems 32 ], Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., eds., 8024–8035, Curran Associates, Inc. (2019).
  * [42] Wolterink, J. M., Leiner, T., Viergever, M. A., and Išgum, I., “Generative adversarial networks for noise reduction in low-dose CT,” IEEE Transactions on Medical Imaging 36(12), 2536–2545 (2017).
  * [43] van Hamersvelt, R. W., Willemink, M. J., Takx, R. A., Eikendal, A. L., Budde, R. P., Leiner, T., Mol, C. P., Isgum, I., and de Jong, P. A., “Cardiac valve calcifications on low-dose unenhanced ungated chest computed tomography: inter-observer and inter-examination reliability, agreement and variability,” European Radiology 24(7), 1557–1564 (2014).


# Using In-Service Train Vibration for Detecting Railway Maintenance Needs

Irene Alisjahbana
Stanford University
<EMAIL_ADDRESS>
###### Abstract

The need for the maintenance of railway track systems have been increasing.
Traditional methods that are currently being used are either inaccurate, labor
and time intensive, or does not enable continuous monitoring of the system. As
a result, in-service train vibrations have been shown to be a cheaper
alternative for monitoring of railway track systems. In this paper, a method
is proposed to detect different maintenance needs of railway track systems
using a single pass of train direction. The DR-Train dataset that is publicly
available was used. Results show that by using a simple classifier such as the
k-nearest neighbor (k-NN) algorithm, the signal energy features of the
acceleration data can achieve 76% accuracy on two types of maintenance needs,
tamping and surfacing. The results show that the transverse direction is able
to more accurately detect maintenance needs, and triaxial accelerometer can
give further information on the maintenance needs. Furthermore, this paper
demonstrates the use of multi-label classification to detect multiple types of
maintenance needs simultaneously. The results show multi-label classification
performs only slightly worse than the simple binary classification (72%
accuracy) and that this can be a simple method that can easily be deployed in
areas that have a history of many maintenance issues.

##  1 Introduction

As most of the railway infrastructure in the US is aging and nearing the end
of its useful life, the need to monitor the health of the railway system is
increasing. In fact, the spending on operation and maintenance have recently
surpassed the spending on the infrastructure itself [1]. Today, most railway
infrastructure monitoring system still rely on two methods: manual visual
assessment or using a special track geometry car. However, visual inspection
is known to be inaccurate in addition to being labor and time intensive. Using
the track geometry car on the other hand, though very accurate, usually
requires that parts of the railway be closed for maintenance.

As a result, newer and more cost efficient methods to monitor railway
infrastructure health have been developed. Some examples include sensor-based
monitoring systems such as optical fiber [10] and strain gauges and image-
based systems such as using SAR images [2] or employing Unmanned Aerial
Vehicles (UAV) [3]. Nevertheless, there are still many limitations to the
aforementioned methods. For example, extremely large amount of optical fibers
must be employed on the railway track to accurately monitor the whole railway
infrastructure system. Image-based systems, though provides information on an
interval basis, cannot provide continuous information on the condition of the
railway infrastructure.

Researchers have found that the health of the railway infrastructure can be
indirectly monitored through the dynamic response of a train that is in
service itself. In addition to being cheaper and more reliable than visual
inspection, the use of sensors that are placed on operational trains enable
the continuous monitoring of the railway health. Furthermore, the availability
of these sensors on the train have now provided researchers with a lot of
data, thus enabling more data-driven approaches. For example, the recently
released DR-Train dataset [7] provides the dynamic response of operational
trains through vibration data and includes other corresponding data including
GPS position and environmental factors. Nonetheless, past work using train
vibration for railway track monitoring have mostly focused on detecting
changes of the railway track and not the underlying reasons causing them.

To that end, this paper aims to study how the different maintenance needs of
railway tracks can be detected using a single pass of train vibration. In
particular, this paper addresses the following questions: (1) which features
and acceleration direction that is best at capturing the characteristics of
the different maintenance needs and (2) how to detect multiple maintenance
needs simultaneously.

##  2 Related Work

In 2017, Lederman et al. [4] used uniaxial acceleration data obtained from
sensors on board operational trains to detect 2 different changes on the
railway track, particularly from track replacement and tamping, two common
maintenance procedures of railway track systems. In their study, they found
that the two main challenges of using train acceleration data are that the
speed of the train varies from pass to pass, and that the position of the
train is not known precisely as a result of the uncertainties of the GPS
system. As a result, Lederman et al. aligned the raw acceleration spatially to
reduce the uncertainty caused by the variability in train speed. Using this
spatially aligned acceleration data, they observed that the signal-energy best
represents track change which is visible after multiple passes.

In a subsequent study, Lederman et al. [6] proposed a novel technique to
reduce the errors caused by the uncertainty of the GPS position, a major
challenge that was identified in their earlier study. Here, Lederman et al.
used the roughness and bumps that are found on the track profile to align the
measurements. The inverse problem was solved by modelling the train as a
simple damped oscillator and constraining the problem to have a limited number
of bumps.

Lederman et al. [5] also presented a technique to combine features obtained
through multiple sensors in order to achieve better detection accuracy. They
used a two-step approach in which the data is first aligned based on the GPS
position before using a novel adaptive Kalman filter that weighs the data
based on its reliability to fuse the data together. By using this method,
their results indicated that fusing data from multiple trains enable the
change to be detected earlier.

Liu et al. [8] proposed an alternative method to reduce the errors caused by
the uncertainty of the GPS position. Given the ground truth track coordinates
obtained from the authorities, the raw data were aligned through the iterative
closest point (ICP) algorithm to reduce the error of of the GPS positions.
Using a variational autoencoder (VAE), they used the aligned data to detect
anomalies in the signal energy of slopes of the track geometry.

Even so, the studies done by Lederman et al. [4, 6, 5] and Liu et al. [8] both
focused on detecting the change itself and not the underlying reasons causing
it, particularly whether the changes in the signal correspond to a type of
damage that occurred or some maintenance need that needs to occur. In
addition, only a single type of maintenance work or anomaly can be detected at
a single time.

##  3 Dataset

###  3.1 Dataset Description

The DR-Train dataset is a public dataset first released by Liu et al. [7]. The
dataset contains measurements obtained from 2 Light Rail Vehicles in
Pittsburgh, Pennsylvania. In addition to vibration data collected by
accelerometers, the dataset also includes the GPS positions, corresponding
environmental conditions (temperature and weather) as well as ground truth
maintenance records provided by the Port Authority of Allegheny County (PAAC).

As explained in [7], the complete data management system consists of the
sensing module, data-acquisition module and data-storage module and data-
processing module.

The sensing module that was used to obtain the measurements are as follows:

  * •

LRV 4306: two uni-axial accelerometers inside the cabin of the train
(VibraMetrics 5102) and a tri-axial accelerometer (PCB 354C03) on the central
wheel truck

  * •

LRV 4313: two uni-axial accelerometers (VibraMetrics 5102) and two tri-axial
accelerometers (PCB 354C03)

  * •

Both LRV’s were instrumented with a BU 353 GPS antenna

Detailed placement of the sensors on the two LRV’s can be seen in [7].

For the data-acquisition module, the National Instruments LabView was used.
The system samples acceleration signals at 1.6kHz whereas the GPS position was
sampled at 1 Hz. The data was downloaded manually to an external hard drive
every two weeks. A sample of the raw acceleration data obtained from sensor 1
of LRV 4313 can be seen in Figure 1.

Figure 1: Sample raw acceleration data from sensor 1 of LRV 4313

###  3.2 Location of Interest

One major issue encountered in this study is related to the ground truth
maintenance logs. As explained in Section 3.1, the raw acceleration data also
contains the corresponding GPS coordinates. The ground truth data however,
only contains the descriptive locations and not the exact GPS coordinates (for
example, Panhandle Bridge 1A). Furthermore, the location of interest had to be
subjected to multiple types of maintenance work. However, the maintenance
record shows that only a few locations have multiple type of maintenance work
performed.

As a result, the locations provided in the maintenance record had to be
manually matched with Google Maps. Taking into account the number of
maintenance types performed, the location of interest selected for this
project is a location with the description: Denise crossover and Glenbury
crossover. This location has a recorded maintenance work of tamping and
surfacing. The location of interest with the corresponding Google Maps for
reference can be seen in Figure 2. A 250m track portion near that area was
used to extract the sample data.

Figure 2: Portion of track chosen for samples as compared to Google Maps
location

###  3.3 Maintenance Record Ground Truth

Figure 3 shows the dates of the specific maintenance work that was performed
on both the outbound and inbound track as specified by the maintenance
records.

Surfacing refers to the maintenance of the tracks in order to have the correct
surface elevation, such as raising or lowering the track structure. On the
other hand, tamping refers to the work that is done to pack the track ballast
underneath the track itself.

Train passes that occurred within a two-month range before the maintenance
work occurred were chosen for the data samples. For each sensor, a total of 17
passes were obtained for tamping, 14 passes were obtained for surfacing, 29
passes for both tamping and surfacing, and 27 passes for no maintenance work.

Figure 3: Selection of passes with respect to ground truth maintenance work

##  4 Methods

###  4.1 Pre-Processing

For the data pre-processing, the passes are first aligned with the ground
truth tack coordinates using the ICP algorithm according to the method
proposed by Liu et al. [8]. Passes that have incorrect GPS positions are
removed. The results of the GPS alignment can be seen in Figure 4.

The acceleration data are then aligned spatially to reduce the uncertainties
caused by the variability in train speed, similar to what was done in Lederman
et al. [4].

For all experiments, the dataset is split into training and validation sets
with 80% to 20% ratio.

Figure 4: The passes before alignment and after alignment using the iterative
closest point (ICP) algorithm. Passes that have bad GPS data are removed.

###  4.2 Feature Extraction

For the individual sensor directions, the signal-energy that was first
proposed in [4] is used as the features. It has been shown that this feature
is more robust to position uncertainty for different types of maintenance
works. The feature can be written as follows:

$f_{n}=\ddot{u}_{n}^{2}|_{x_{\epsilon}}$

where $\ddot{u}_{n}^{2}$ is the raw acceleration data and ${x_{\epsilon}}$
indicates the acceleration of the data given the position of the train
(spatially interpolated as opposed to temporally interpolated).

Following the methods done by Lederman et al, the signal is first averaged
using a 25m window of track. Then, only the 50 most discriminative indices are
selected. This was done by taking the indices with the greatest difference
with the mean signal.

Given the limited number of data samples in the training and test sets,
principal component analysis (PCA) is then performed in order to reduce the
dimensionality of the features. The explained variance for the number of
components included can be seen in Figure 5. It can be seen that using only 2
components, over 95% of the data sample variance can be explained, and
therefore only the first 2 components are used as features.

For the tri-axial accelerometer data, two types of features are defined:
concatenation and major voting. In concatenation, the 50 most discriminative
indices of each direction (x-, y-, and z-axis) are first concatenated in order
to obtain samples with 150 features each. PCA is then performed on the
concatenated features to obtain the first 2 components. In major voting, PCA
is performed for each individual sensor direction (x-, y-, and z-axis). The
final classification is obtained through the major voting result of each
individual sensor.

Figure 5: PCA results of the explained variation vs number of components for
the signal energy feature

###  4.3 Machine Learning Algorithms

####  4.3.1 Binary Classification

To address the first objective of the paper in which features and acceleration
direction are studied to detect maintenance needs, binary classification (no
maintenance vs maintenance) is performed. Classification models such as
k-Nearest Neighbours, Logistic Regression and Support Vector Machines (SVM)
are implemented.

k-NN: k-NN is a simple classification model in which the data are group
according to the notion of distance. A data point is classified into a
specific class by taking the ”k” nearest neighbours and taking a majority vote
based on the class of its neighboring data points. The k-NN model works well
for smaller data sets and can produce decision boundaries that are non-linear
and complex. It is also a non-parametric model and does not make any
assumptions of the underlying sample distribution. Given the limited number of
data samples, k-NN is an ideal starting point for the classification of
maintenance needs. k-NN does not perform as well in higher dimensions, but
given that the features of the acceleration data has been reduced to just two
components, k-NN is still applicable.

Logistic Regression: Logistic regression uses a logistic function to model a
binary dependent variable. Unlike k-NN, the logistic regression is only able
to produce linear boundary decisions. In this study, the L2 norm is used as
the penalization. Furthermore, the value C, or the regularization strength, is
used as the hyperparameter.

SVM: The SVM is another supervised classification algorithm that uses a subset
of training points in the decision function. The advantage of SVM is that it
is effective even in high dimensional cases, and more robust to outliers.
Similar to logistic regression, the regularization strength C is a
hyperparameter. In this study, the radial basis function (RBF) kernel is used
to produce more non-linear boundaries. The kernel coefficient $\gamma$ is
another hyperparameter that needs to be tuned through cross-validation.

####  4.3.2 Multi-Label Classification

In practice, railway tracks can be subjected to multiple maintenance needs
simultaneously. As a result, multi-label classification can be implemented to
capture this scenario. Unlike regular classification in which each data point
has one corresponding outcome label, in multi-label classification, a single
data point can have multiple outcome labels (in this case tamping and
surfacing). As a result, combinations of the maintenance need can be used to
detect no maintenance, tamping, surfacing, or both (Figure 6).

Figure 6: Data samples and corresponding outcome labels for multi-label
classification

Two different methods were used to implement multi-label classification:
problem transformation and adaptive algorithm. As its name suggests, in
problem transformation, the classification problem is transformed from a
multi-label classification into regular classification with one outcome value.
This can be done through several methods, which are binary relevance,
classifier chains and label powerset. In binary relevance, each label is
treated as a seperate single class classification problem, and multiple models
are trained for each label. In chain classifiers, the next classifier is
trained on the input space and all the previous classifiers in the chain.
Unlike binary relevance, chain classifier can take into account correlations
between the label. Label powerset on the other hand, treats the outcomes as
unique label combinations. The k-NN model is used as the basis model for all
the methods.

Adaptive algorithm on the other hand does not transfrom the problem, but
instead derives regular classification methods for multi-label classification.
Using the base model of k-NN, the ML-kNN [11], first proposed in 2007, is a
classifier that is used to predict multi-label data. The ML-kNN model was
implemented using the scikit-multilearn package in Python [9].

##  5 Experiments

###  5.1 Evaluation Metrics

The evaluation metrics used to reflect the overall performance of the binary
classification model is accuracy, which is the ratio of number of correct
predictions to the total number of input samples. It can also be written in
the following equation:

$Accuracy=\frac{TP+TN}{TP+TN+FP+FN}$ (1)

where TP = True Positives, TN = True Negatives, FP = False Positives, and FN =
False Negatives.

On the other hand, the error metric for the multi-label classification is
Exact Matio Ratio (subset accuracy) and Hamming loss. The exact match ratio is
the most strict metric and calculates the percentage of samples that have all
their labels classified correctly. It ignores partially correct matches. The
exact match ratio can be written as follows:

$ExactMatchRatio=\frac{1}{n}\sum_{i=1}^{n}I(Y_{i}=Z_{i})$ (2)

where $Y_{i}$ is the target label and $Z_{i}$ is the predicted label.

The Hamming loss measures the fractions of labels that are incorrectly
predicted and can be written as follows:

$HammingLoss=\frac{1}{|N|.|L|}\sum_{i=1}^{|N|}\sum_{j=1}^{|L|}xor(y_{i,j},z_{i,j})$
(3)

where $y_{i,j}$ is the target and $z_{i,j}$ is the prediction.

The metrics are calculated for both the training set and testing set.

###  5.2 Binary Classification

k-NN classifies a data point based on its closest neighbors, therefore, the
parameter ”k” will highly determine the performance of the model. Thus, the
number of neighbour was tuned through using cross-validation. This is done by
dividing the training set into k-folds, training it on (k-1) folds and holding
1 fold for the validation set. The accuracy is calculated for all k folds, and
the mean is taken to obtain the misclassification error for the chosen
hyperparameter. A sample result for the cross validation using 3-folds can be
seen in Figure 7. The number of folds for the cross validation was chosen by
taking into account the total number of samples in the training set.

Figure 7: Cross validation results to determine number of k

The results for the different accelerometer directions can be seen in Table 1.
It can be seen that the transverse direction, the direction pointing
perpendicular to the railway actually performs slightly better than the
vertical direction - the direction that has always been used in past studies.
This however, makes sense as tamping and surfacing are maintenance works that
are performed to make a smoother ride. The more unnecessary movements
occurring in any other direction indicates that maintenance is needed.

Sensors | Optimal k | Train Acc. | Test Acc.
---|---|---|---
z-axis (vertical) | 8 | 0.76 | 0.73
x-axis (longitudinal) | 7 | 0.84 | 0.70
y-axis (transverse) | 7 | 0.85 | 0.74
Table 1: Accelerometer direction binary classification results

Table 2 shows the classification results for the different accelerometer
types, particularly uniaxial and triaxial accelerometers. As described in
section 4.2, two methods were implemented for the triaxial accelerometer:
concatenation and major voting. It can be seen that either method performs
better than the uniaxial accelerometer (only vertical direction), which
further proves that the directions other than the vertical can provide more
information to the maintenance needs of the railway track.

Sensors | Optimal k | Train Acc. | Test Acc.
---|---|---|---
Uniaxial | 3 | 0.84 | 0.71
Triaxial (Concat) | 5 | 0.80 | 0.76
Triaxial (Voting) | varies | 0.82 | 0.76
Table 2: Accelerometer type binary classification results Figure 8: Training
and test set with the decision boundaries for the triaxial accelerometer using
k-NN model Figure 9: Confusion matrix for triaxial accelerometer using k-NN
model

Figure 8 shows the decision boundary of the k-NN model. It can be seen that
the decision boundary is quite linear, with the ”no maintenance” sample
located on the left side of the plots. As a result, other classification
models such as logistic regression and SVM with a radial basis function (RBF)
kernel were implemented to improve the accuracy of the models. Similarly to
the k-NN model, cross validation with 5 folds was performed to obtain the
optimal parameters.

For logistic regression model, the L2 penalty was used and the hyperparameter
that had to be tuned is C, or the regularization strength. Cross validation
results show that the optimal C value is 0.62. On the other hand, for the SVM
model, the hyperparameters to be tuned is the C and the $\gamma$. Cross
validation results show that the optimal C value is 10, with a $\gamma$ value
of 0.01.

Comparison results for the three methods can be seen in Table 3. All the
method performs just as well in the test set, though the k-NN classifier had a
higher train accuracy. This can be seen through the decision boundaries for
the logistic regression (Figure 10) and SVM (Figure 11) that are similar to
that of the k-NN model (Figure 8. The higher training accuracy of the k-NN
classifier can just be caused by some overfitting to the training set.

Model | Train Acc. | Test Acc.
---|---|---
k-NN | 0.80 | 0.76
Logistic Regression | 0.72 | 0.76
SVM RBF Kernel | 0.74 | 0.76
Table 3: Binary classification results using different classification models
Figure 10: Training and test set with the decision boundaries for the triaxial
accelerometer using Logistic Regression Figure 11: Training and test set with
the decision boundaries for the triaxial accelerometer using SVM with RBF
kernel

Given the uncertainty of the location of the maintenance work itself, a
sensitivity analysis on the window size or length of track portion to take
into account was also performed. From Table 4 it can be seen that by
increasing the length of track, the accuracy is reduced, and is most likely
due to more noise being included. When thinking about how the method would be
implemented in real life, the length of track to take into account will
greatly determine the detection accuracy.

Train Pass Length | Train Acc. | Test Acc.
---|---|---
250m | 0.80 | 0.76
500m | 0.71 | 0.72
1000m | 0.72 | 0.72
1500m | 0.76 | 0.50
Table 4: Binary classification results using varying window sizes (length of
track portion)

.

###  5.3 Multi-Label Classification

Results of the multi-label classification with the Exact Match Ratio error
metric can be seen in Table 5, whereas results using Hamming loss can be seen
in Table 6. Though the Exact Matio Ratio is a very strict error metric that
does not take into account partial correctness, it can be seen that the best
performing model was able to achieve a 72% accuracy, only 4% lower than
ordinary binary classification. This shows that the signal-energy feature of
the train vibration is able to be used for multiple maintenance needs
detection.

Chain classifier performs slightly better overall, as it took into
consideration label correlation. This would mean that a higher probability
that the signal required tamping, means it also needed surfacing and vice
versa. Label Powerset and ML-kNN did not perform as well probably due to the
limited number of samples.

Method | Train Acc. | Test Acc.
---|---|---
Binary Relevance | 0.67 | 0.72
Chain Classifier | 0.70 | 0.72
Label Powerset | 0.70 | 0.56
ML-kNN | 0.62 | 0.62
Table 5: Multi-label classification results using triaxial accelerometer and Exact Match Ratio error metric Method | Train Loss | Test Loss
---|---|---
Binary Relevance | 0.22 | 0.22
Chain Classifier | 0.21 | 0.22
Label Powerset | 0.22 | 0.31
ML-kNN | 0.28 | 0.31
Table 6: Multi-label classification results using triaxial accelerometer and
Hamming Loss error metric

##  6 Limitations and Future Improvements

In reality, the methods that have been described in this paper will be easily
applicable and trainable to new areas. However, the methods will not be able
to perform real-time classification on the track signals. This is partly
because of the underlying machine learning method used (kNN), which is a lazy
classifier, but also because the features that was used in the model. Though
the features can be easily extracted from the raw acceleration and GPS
signals, the features must first be aligned based on the GPS coordinates as
well as the speed of the train. This prevents the ability of real-time
monitoring.

In line with real-time monitoring, the methods described would be unable to
accurately classify the signals during the transition periods. In other words,
there is a grey time period between when the railway track does not need
maintenance and when the track does need maintenance. This is because a
specific time range (around 2 months) prior to the maintenance work was
labeled as ”need maintenance” and anything before that was labeled as ”no
maintenance”. However, in reality, the threshold would be more gradual, and no
clear distinction would be made. A deeper study on the effects of the time
window on the accuracy of the models can be performed. Perhaps one would use
passes from multiple days instead of a single day to see whether the model
would consistently predict ”need maintenance”. Thus, the confidence levels of
the classification would increase if the track portion is consistently
classified as ”need maintenance”.

Furthermore, only one location, the Denise Crossover - Glenbury Crossover, was
used as a case study in this study. Therefore, the ability for the model to
work well on other locations is still unknown. The model will most likely have
to be retrained to be implemented in other locations. Finally, only two types
of maintenance work, tamping and surfacing was taken into account in this
study. In reality, there are many other types of maintenance work that a
railway track can be subjected to. Including more locations and types of
maintenance needs can be an extension to this study.

As described in Section 2, one of the main challenges of this study is the
unavailability of GPS coordinates corresponding to the maintenance work.
Obtaining the corresponding GPS coordinates associated with each location of
the maintenance work would be extremely beneficial to generalizing the models.
This might mean some additional agreement with the PAAC authorities as exact
track damage locations might not be available to the public. However, if this
information were available, more location and maintenance need types can be
included into the training and testing dataset.

##  7 Conclusion

In this study, we demonstrated how detecting maintenance needs, particularly
tamping and surfacing, of railway tracks can be performed using a single pass
of in-service train vibration. The results show that the transverse direction
was able to provide more information on the maintenance needs of railway
tracks as compared to the vertical direction, the direction that has mostly
been used in past studies. This was further supported by the highest accuracy
obtained by the triaxial accelerometer with an accuracy of 76% using a simple
k-NN classifier. The results also show that increasing the length of track
taken into consideration will drastically reduce the classification accuracy.

Finally, this study showed that multiple maintenance needs can be detected
simultaneously through multi-label classification while only reducing the
overall accuracy by 4%.

Though this study only demonstrated the method on a single location and only
two maintenance needs, it is a step towards an early warning detection system
for detecting maintenance needs of railway track systems. The methods
presented in this study is easily trainable and applicable to any location
that have a history of many maintenance issues.

##  8 Acknowledgments

I would like to acknowledge the Stanford University teaching team of CEE 286
(Structural Health Monitoring) particularly Professor Haeyoung Noh and Mostafa
Mirshekari for the constant feedback on the project.

## References

  * [1] The state of u.s. infrastructure.
  * [2] L. Chang, R. P. Dollevoet, and R. F. Hanssen.  Nationwide Railway Monitoring Using Satellite SAR Interferometry.  IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 10(2):596–604, 2017.
  * [3] F. Flammini, C. Pragliola, and G. Smarra.  Railway infrastructure monitoring by drones.  2016 International Conference on Electrical Systems for Aircraft, Railway, Ship Propulsion and Road Vehicles and International Transportation Electrification Conference, ESARS-ITEC 2016, pages 1–6, 2017.
  * [4] G. Lederman, S. Chen, J. Garrett, J. Kovačević, H. Y. Noh, and J. Bielak.  Track-monitoring from the dynamic response of an operational train.  Mechanical Systems and Signal Processing, 87:1–16, 2017.
  * [5] G. Lederman, S. Chen, J. H. Garrett, J. Kovačević, H. Y. Noh, and J. Bielak.  A data fusion approach for track monitoring from multiple in-service trains.  Mechanical Systems and Signal Processing, 95:363–379, 2017.
  * [6] G. Lederman, S. Chen, J. H. Garrett, J. Kovačević, H. Y. Noh, and J. Bielak.  Track monitoring from the dynamic response of a passing train: a sparse approach.  Mechanical Systems and Signal Processing, 90:141–153, 2017.
  * [7] J. Liu, S. Chen, G. Lederman, D. B. Kramer, H. Y. Noh, J. Bielak, J. H. Garrett, J. Kovačević, and M. Bergés.  Dynamic responses, gps positions and environmental conditions of two light rail vehicles in pittsburgh.  Scientific data, 6(1):1–11, 2019.
  * [8] J. Liu, Y. Wei, M. Bergés, J. Bielak, J. H. Garrett Jr, and H. Noh.  Detecting anomalies in longitudinal elevation of track geometry using train dynamic responses via a variational autoencoder.  In Sensors and Smart Structures Technologies for Civil, Mechanical, and Aerospace Systems 2019, volume 10970, page 109701B. International Society for Optics and Photonics, 2019.
  * [9] P. Szymański and T. Kajdanowicz.  A scikit-based Python environment for performing multi-label classification.  ArXiv e-prints, Feb. 2017.
  * [10] X. W. Ye, Y. H. Su, and J. P. Han.  Structural health monitoring of civil infrastructure using optical fiber sensing technology: A comprehensive review.  Scientific World Journal, 2014, 2014.
  * [11] M.-L. Zhang and Z.-H. Zhou.  Ml-knn: A lazy learning approach to multi-label learning.  Pattern recognition, 40(7):2038–2048, 2007.