# CNN Based Segmentation of Infarcted Regions in Acute Cerebral Stroke
Patients From Computed Tomography Perfusion Imaging

Luca Tomasetti<EMAIL_ADDRESS>,  Kjersti Engan<EMAIL_ADDRESS>,
Mahdieh Khanmohammadi<EMAIL_ADDRESS>University of Stavanger,
Department of Electrical Engineering and Computer Science, BMDLab
StavangerNorway4013 and  Kathinka Dæhli Kurz  Stavanger University Hospital,
Stavanger medical imaging laboratory (SMIL) University of
StavangerStavangerNorway4013<EMAIL_ADDRESS>
(2020)

###### Abstract.

More than 13 million people suffer from ischemic cerebral stroke worldwide
each year. Thrombolytic treatment can reduce brain damage but has a narrow
treatment window. Computed Tomography Perfusion imaging is a commonly used
primary assessment tool for stroke patients, and typically the radiologists
will evaluate resulting parametric maps to estimate the affected areas, dead
tissue (core), and the surrounding tissue at risk (penumbra), to decide
further treatments. Different work has been reported, suggesting thresholds,
and semi-automated methods, and in later years deep neural networks, for
segmenting infarction areas based on the parametric maps. However, there is no
consensus in terms of which thresholds to use, or how to combine the
information from the parametric maps, and the presented methods all have
limitations in terms of both accuracy and reproducibility.

We propose a fully automated convolutional neural network based segmentation
method that uses the full four-dimensional computed tomography perfusion
dataset as input, rather than the pre-filtered parametric maps. The suggested
network is tested on an available dataset as a proof-of-concept, with very
encouraging results. Cross-validated results show averaged Dice score of 0.78
and 0.53, and an area under the receiver operating characteristic curve of
0.97 and 0.94 for penumbra and core respectively.

††journalyear: 2020††copyright: acmcopyright††conference: Proceedings of the
11th ACM International Conference on Bioinformatics, Computational Biology and
Health Informatics; September 21–24, 2020; Virtual Event, USA††booktitle:
Proceedings of the 11th ACM International Conference on Bioinformatics,
Computational Biology and Health Informatics (BCB ’20), September 21–24, 2020,
Virtual Event, USA††price: 15.00††doi: 10.1145/3388440.3412470††isbn:
978-1-4503-7964-9/20/09††copyright: none Figure 1. A general overview of the
various steps involved in the proposed approach (green panels) that bypasses
the current methods used by the majority of the radiologists and state-of-the-
art approaches based on thresholding and semi-automatic approaches (red
panel).

##  1\. Introduction

A cerebral stroke is the second most common cause of death among adults
worldwide (Wang et al., 2016). An ischemic stroke can occur if the flow of
oxygen-rich blood in an artery to a portion of the brain is occluded. Ischemic
areas are usually heterogeneous with areas that are irreversibly damaged
(infarct core) and areas where the tissue is still vital, but critically
hypoperfused. This area is called the penumbra, and is the target for therapy,
as restoring blood flow for the penumbra will preserve neurological function
for the patient. This tissue will turn into infarct core if blood flow is not
restored timely, and therefore, there is a “time is brain” concept (Kurz et
al., 2016). A patient can lose up to 1.9 million neurons, 14 billion synapses,
and 12 km nerve fibers every minute from the time the stroke happened (Saver,
2006); the best treatment window for thrombolysis is estimated to be from 3 up
to 4.5 hours from symptom onset (Hacke et al., 2008).

Computed Tomography Perfusion (CTP) is a very fast and valuable tool to
accurately predict the prognosis of patients in the early stage of the
treatment (Campbell et al., 2013). Thus, Computed Tomography (CT) is the most
commonly used initial imaging method for acute stroke patients.

Figure 2. The 4D CTP volume is seen as 3D slices through the brain over
different time points. To the right, an example of a parametric map (TTP) is
seen, calculated from the entire time-span.

The four-dimensional (4D) examination of CTP is formed by acquiring three-
dimensional (3D) CT scans of the brain at many time points during the passage
of a contrast agent from the arteries through the capillaries to the veins and
then into the venous sinuses (Kurz et al., 2016; Campbell et al., 2013).
Parametric color-coded maps are calculated and all the different time points
are included in the analyses, describing the changes of blood perfusion and
other important factors over time.

Various parametric color-coded maps exist describing precise phenomena in the
brain: time to peak (TTP), cerebral blood volume (CBV), relative CBV (rCBV),
cerebral blood flow (CBF), relative CBF (rCBF), mean transit time (MTT) and
time to maximum $(T_{\text{Max}})$ (Kurz et al., 2016). A meticulous
assessment of these parametric maps, performed by a medical specialist, is
needed to identify the ischemic regions of a stroke. These findings will guide
the decision on who needs immediate thrombolytic treatment and/or
interventional thrombectomy to prevent larger cerebral stroke. Generation of
these parametric maps and their inspections by radiologists to make a final
decision about possible patient treatment is time-consuming.

Extensive research has proposed different methods to extract the core and the
penumbra based on distinct parametric maps (Kasasbeh et al., 2019; Lucas et
al., 2018; Campbell et al., 2012; Cereda et al., 2016; Ma et al., 2019; Lin et
al., 2014; Wintermark et al., 2006). Campbell et al. (Campbell et al., 2012),
Cereda et al. (Cereda et al., 2016), Ma et al. (Ma et al., 2019), Lin et al.
(Lin et al., 2014), and Wintermark et al.(Wintermark et al., 2006) are all
examples of researches that use extensive and different thresholding
approaches to define the infarct regions; hence, a global consensus in the CTP
parameters does not exist to identify the infarcted core and penumbra.

In all the mentioned studies, the gold-standard was found from Magnetic
Resonance Imaging (MRI) using the Diffusion-weighted Imaging (DWI) sequences
generated some hours after the acquisition of CTP images. DWI sequences
highlight the irreversibly ischemic areas; however, the delay in the
acquisition of these images may cause a mismatch between what is the infarcted
region of the gold-standard, and what was the actual hypoperfused region at
the time of the CTP acquisition, which may create an important obstacle in the
diagnostic work-up of these patients (Cereda et al., 2016). An improved
evaluation of these regions may lead to a better selection of patients for
treatment in the acute phase.

Image segmentation approaches using convolutional neural networks (CNN)
methods, as the backbone, have become a leading research field in the past
years. They have shown to provide promising results in many biomedical imaging
applications, such as the famous U-net model (Ronneberger et al., 2015) for
two-dimensional (2D) segmentation of neuronal structures in electron
microscopic stacks, the V-Net structure (Milletari et al., 2016) introduced to
segment prostate MRI volumes, the 3D U-Net network for volumetric segmentation
by Çiçek et al. (Çiçek et al., 2016), and the 3D CNN architecture for brain
lesion segmentation suggested by Kamnitsas et al. (Kamnitsas et al., 2017).
Specifically, in the context of stroke patients, Vargas et al. (Vargas et al.,
2019) proposed the use of CNN to predict the presence of perfusion deficit;
Kasasbeh et al. (Kasasbeh et al., 2019) used a semi-automatic approach based
on a CNN with the parametric maps as input. Lucas et al. (Lucas et al., 2018)
introduced a CNN based on the 2D U-Net with multi-scale information to segment
the infarct core inside the brain, also using a set of parametric maps as
input and other information. Also, CNN has been explored on the use of
Computed Tomography Angiography (CTA) images, enhancing occlusions in the
blood vessels (Öman et al., 2019; Barman et al., 2019; Sheth et al., 2019).

In this work, we propose a fully automatic CNN-based method to segment both
the infarct core and penumbra regions, using the complete 4D set of CTP scans
over the injection period. Our method contributes to the following:

  1. (1)

It uses CTP slices over time as input since this is normally the first
investigation a patient suspected to stroke goes through (Kurz et al., 2016;
Committee et al., 2008).

  2. (2)

The 4D set of CTP scans contains original information compared to the
calculated parametric maps.

  3. (3)

Being independent from parametric maps saves essential time during the
investigation of a patient’s condition.

  4. (4)

Our approach is based on a 3D CNN in a bottleneck model inspired by the U-Net
(Ronneberger et al., 2015) using as gold-standard manually annotated images
generated by an experienced radiologist after studying the parametric maps.

To the best of our knowledge, this is the first attempt that uses directly CTP
brain slices, and segments automatically the infarct core and penumbra regions
together, both essential to guide the decision of the patient’s treatment.

##  2\. Data Material

The data material consists of 4230 CTP slices, stored as DICOM files, from 10
pseudonymized patients with large vessel occlusions collected at Stavanger
University Hospital between 2014 and 2015. The mean age was 70.4 years (SD,
8.9), 60% were female and 50% of the patients had an acute stroke in the left
hemisphere. The median baseline National Institutes of Health Stroke Scale
(NIHSS) score was 14.5.

Fig. 2 presents an overview of a 4D CTP brain; slices of size $(512\times
512\times S)$, where $S$ is the number of slices $(s_{1},s_{2},\ldots,s_{S})$,
during the acquisition time $t_{1},t_{2},\ldots,t_{T}$ after injecting
contrast agent in a cubital vein. Using the information on a sample section
over time, from $t_{1}$ to $t_{T}$, the parametric maps are generated.

###  2.1. Imaging protocol

CTP images were acquired during injection of 40 ml iodine containing contrast
agent (Omnipaque 350 mg/ml) and 40 ml isotonic saline in a cubital vein with a
flow rate of 6 ml/s. The delay for the first slice acquisition was 4 seconds.
A number between 13 and 22 slices $(S)$, covering all brain, is acquired from
each patient with a slice thickness of 5 mm. Each slice was captured 30 times
($t_{1},\ldots,t_{T}$, where $T=30$) over the injection period. The total
number of images per patient is between 390 and 660, depending on the number
of slices acquired. Eight patients were subjected to 416 mA for the x-ray tube
current, while two patients were subjected to 350 mA; all of them were exposed
to 80 kVp.

The ground truth for each brain slice $s_{i}$ is an image with manually
annotated infarcted areas realized by an experienced radiologist from the
Stavanger University Hospital after a detailed study of the corresponding
parametric maps: TTP, CBV, CBF, and $T_{\text{Max}}$. The annotated areas in a
ground truth image (rightmost column of Fig. 3) are the brain without stroke
symptoms (black), the penumbra (dark gray), the core (light gray) and the
background (white), with a target of 0, 76, 150, 255 respectively in grayscale
pixel value.

Figure 3. A brain slice $s_{i}$ and the related ground truth image
$Y_{GT_{i}}(m,n)$ are processed through a sliding window technique to create a
series of 3D matrices $\tilde{V}_{s_{i}}(t,m_{j},n_{j})$ for each tile $j$ and
the relative output $Y_{GT_{i}}(m_{j},n_{j})$.

##  3\. Proposed Method

We propose a method for fully automatic segmentation of the core area of the
stroke as well as the penumbra area from the 4D CTP dataset with a CNN founded
method.

Let a CTP data be defined as, $V\in Z^{4}$, 4D signal of dimension $(T\times
M\times N\times S)$. $T$ is the number of time points $(T=30)$, $M\times N$ is
the dimension of each 2D image slice and $S$ the number of slices in the
brain. The image slices are of size $512\times 512$ with a resolution of
0.4258 mm/pixel. The slice thickness is 5 mm, and as such, the volume has a
much lower resolution in the $z$ dimension. At this stage we suggest doing a
slice by slice segmentation, but using the full-time series of images of each
slice.

Let $V_{s_{i}}(t,m,n)$ denote the 3D data for slice $s_{i}$ where $t$ denotes
the time index, and $(m,n)$ the spatial image coordinates. After pre-
processing steps (see Sec. 3.1) it is called $\tilde{V}_{s_{i}}(t,m,n)$. We
define $A_{j}$ as a set of image coordinate pairs, $(m,n)$, corresponding to a
tile number $j$ of size $16\times 16$ pixel, which extracted from a slice
$s_{i}$. Thus $\tilde{V}_{s_{i}}(t,m,n)_{(m,n)\in A_{j}}$ gives a volume of
size $(30\times 16\times 16)$ from slice $s_{i}$ to be used as the input of
the neural network. For notation simplicitiy we write:
${\tilde{V}_{s_{i}}(t,m,n)}_{(m,n)\in A_{j}}=\tilde{V}_{s_{i}}(t,m_{j},n_{j})$
in the reminder of the paper. For the 3D output and for the correlated ground
truth, the notation becomes $Y_{s_{i}}(m,n)$ and $Y_{GT_{i}}(m,n)$ for the
entire slice $s_{i}$; while for tile $j$ the notation becomes
$Y_{s_{i}}(m_{j},n_{j})$ and $Y_{GT_{i}}(m_{j},n_{j})$ respectively. Note that
the time dimension is no longer there, as the data is transformed from 4D to
3D through the analysis.

Due to the limited size of the dataset, a sliding window technique is applied
to each slice $s_{i}$ over the spatial coordinates during training to increase
the sample size and create a complete set of 3D matrices composed of $16\times
16$ tiles. During prediction using the learning model, tiles can be
overlapping, have different strides, or be non-overlapping, this can be taken
into account when composing the final output.

An example of one input and the relative ground truth image of the proposed
network is given in Fig. 3: a brain slice $s_{i}$ is processed through the
time series $(t_{1},t_{2}\ldots,t_{T})$ to generate a 3D matrix
$\tilde{V}_{s_{i}}(t,m_{j},n_{j})$, after pre-processing steps. The ground
truth, described as $Y_{GT_{i}}(m_{j},n_{j})$, represents tile $j$ in the
manually annotated brain slice (rightmost column of Fig. 3).

###  3.1. Pre-processing steps

A series of pre-processing steps are done before the CNN module and before
applying the overlapping sliding window technique to the dataset. Firstly
$V_{s_{i}}(t,m,n)$ is pairwise registered over the time sequence
$(t_{1},\ldots,t_{T})$ by a robust similarity-based method (Goshtasby, 2005).
Subsequently, the skull is removed using a watershed segmentation algorithm
(Vincent and Soille, 1991). Let $\hat{V}_{s_{i}}(t,m,n)$ denote the 3D
(2D+time) CTP data for slice $s_{i}$, after all these steps.

At last, $\hat{V}_{s_{i}}(t,m,n)$ is given as input to a contrast enhancement
algorithm, which involves three steps:

  1. (1)

image intensity adjustment to increase the contrast by saturating 1% of the
image to low and high intensities;

  2. (2)

histogram equalization (Kim, 1997) to distribute the intensities;

  3. (3)

image normalization to harmonize the contrast among all the images.

The input after the pre-processing steps is defined as
$\tilde{V}_{s_{i}}(t,m,n)$.

###  3.2. Architectures

Due to a limited amount of ground truth labeled data to train on, we pre-test,
as the first setup, three different CNN-based architectures followed by fully-
connected layers for classifier network, where an image block is classified as
belonging to the class “core”, “penumbra”, “brain”, or “background”. The
compressor part of the best architecture for classification is used as the
compressor part in a segmentation network, where a decompressor part is added,
but producing a 3D output from the 4D input. All the architectures are
implemented and trained with the same preprocessed input; they also produce
outputs interpolated to the same size.

Figure 4. Pipeline overview of the pre-test of different architectures
together with the details of their layers. (1) shows the input structure, (2)
defines the layers of the three architectures, (3) is the probability vector
output, (4) displays the tile expansion for the class with the highest
probability value and (5) represents the final image $Y_{s_{i}}(m,n)$.

A description of the various CNN based classifier networks is given in Fig. 4,
including the pipeline involved in the pre-test process. Arch_1, 2 and 3 are
using the same input described in Fig. 3 (marked with (1) in Fig. 4). The
neural networks generate a probability vector of four elements
$[y_{1},y_{2},y_{3},y_{4}]$ for each input tile $j$ $(T\times 16\times 16)$
(3). The highest prediction value is selected and expanded in the $16\times
16$ tile $j$ with the corresponding labeled grayscale value
$Y_{s_{i}}(m_{j},n_{j})$ (4). The final outcome is a composed pixelated image
$Y_{s_{i}}(m,n)$ of all the generated $16\times 16$ tiles (5).

Table 1. Comparison of the three classifiers based on the average accuracy,
precision, recall, and computational time for each tested patient.

Architecture | Acc.(%) | Prec. (%) | Recall (%) | F1 Score (%) | Time (s)
---|---|---|---|---|---
Arch_1 | 92.34 $\pm$ 4.35 | 83.6 | 75.1 | 79.12 |  $\sim$ 600
Arch_2 | 91.29 $\pm$ 4.71 | 82.4 | 73.8 | 77.86 |  $\sim$ 600
Arch_3 | 89.67 $\pm$ 4.44 | 77.1 | 70.4 | 73.60 | $\sim 700$

Table 2. Comparison of the three pre-tested architectures and the implemented
network (mJ-Net): conv3D represents a 3D convolution layer, max_pool and
avg_pool are a max and an average pooling layer respectively, fully_conn is
the fully connected layer and concat & transpose is the upsamples layer. $T$
is the number of time points in the 3D matrix (30).

| Arch_1 | Arch_2 | Arch_3 | mJ-Net
---|---|---|---|---
Input | 3D matrices | 3D matrices | 3D matrices | 3D matrices
Output |  | Probability
---
vector 4$\times$1
| Probability
---
vector 4$\times$1
| Probability
---
vector 4$\times$1
| $16\times 16$
---
labelled tile
| # classes
---
per tile
1 | 1 | 1 |  $\geq$ 1
| Problem
---
Approach
Classification | Classification | Classification | Segmentation
# Layers | 9 | 10 | 8 | 19
# Parameters | 203,320 | 773,384 | 63,312 | 981,553
Layers |  | conv3D +
---
avg_pool +
(2 $\times$ conv3D) +
max_pool +
conv3D +
max_pool +
(2 $\times$ fully_conn)
| 2 $\times$ [ (2 $\times$ conv3D) +
---
max_pool ] +
conv3D +
max_pool +
(2 $\times$ fully_conn)
|  $3\times$ (conv3D +
---
max_pool) +
(2 $\times$ fully_conn)
| conv3D +
---
avg_pool +
2 $\times$ (2 $\times$ conv3D +
max_pool) +
2 $\times$ conv3D +
2 $\times$ [ (concat & transpose) +
(2 $\times$ conv3D) +
max_pool ] +
conv3D
Kernel Size |  | (T,3,3) first layer
---
(1,3,3) other layers
(3,3,3) | (T,3,3) |  | (T,3,3) first layer
---
(1,3,3) other layers
(1,2,2) concat & transpose layers
| Pooling
---
Window
| (T,2,2) second layer
---
(1,2,2) other layers
(2,2,2) | (2,2,2) |  | (T,1,1) avg_pool
---
(1,2,2) first two max_pool
(2,1,1) last two max_pool
| # of filters
---
per conv3D
| 16, 32, 32, 64
---
| 16, 32, 32, 32, 64
---
| 16, 32, 64
---
| 16, 32, 64, 64, 128,
---
128, 256, 128, 64, 32, 32

The various architectures are illustrated in Fig. 4. Arch_1 reduces the third
dimension of its input into a 2D matrix after the avg_pool layer to average
the entire time dimension of the 3D tile into a 2D matrix. The implementation
of Arch_1 was mainly done to decrease the learning parameters in each layer
for reducing computational time while at the same time have useful information
for training purposes. On the contrary, Arch_2 and 3 work always with a 3D
matrix in their layers. Arch_2 uses a kernel size in each convolution layer of
$(3,3,3)$ while Arch_3 adopts a kernel size of $(T,3,3)$; thus, in each
convolutional layer, Arch_3 convolves a higher number of time points $t$ in
the $\tilde{V}_{s_{i}}(t,m_{j},n_{j})$ matrix than Arch_2.

Our final objective is to find a fast and automatic segmentation approach.
Following this motivation, the classifier model that yields the best results
(Arch_1, as displayed in Table 1) was further expanded to a segmentation
model, inspired both by the U-Net (Ronneberger et al., 2015) and the V-Net
(Milletari et al., 2016). We call it “mirror J-Net” (mJ-Net) because of its
particular shape. Table 2 presents a comparison of the three classifier
structures and the derived mJ-Net architecture, showing in detail the design
of all the layers involved in each model.

A detailed view of the mJ-Net model is given in Fig. 5. The output from mJ-Net
is a matrix of dimension $16\times 16$. Each value in the matrix corresponds
to a probability number which is multiplied by $255$ to represent a grayscale
pixel. A post-processing step is performed to combine all outputs and generate
the final segmented image, which is represented as a grayscale labeled image
with the same dimension as the initial CTP images ($512\times 512$ pixels) and
the labels “core”, “penumbra”, “brain”, “background” (Fig. 7(b)).

Figure 5. Detailed overview of the mirror J-Net (mJ-Net).

(a) Manual annotated image

(b) Automatically predicted images

(c) Histogram of one predicted image.

Figure 7. Example of three brain slices $s_{i}$ comparison between (a) their
ground truth, and (b) the deep neural network output using the best mJ-Net
model; (c) is the histogram used to divide the various classes.

The mJ-Net network is formed by a convolutional and a de-convolutional part.
The convolutional part (left side on the Fig. 5) presents a very similar
structure as the Arch_1: all the layers are the same with the only exception
of the last fully-connected layers, which are substituted by a new
convolutional layer useful for the second part of the network; moreover, the
average pooling operation in the second layer flats the time dimension of the
input into a 2D matrix in the same way as the Arch_1. The de-convolutional
part (right side in Fig. 5) is composed of concatenations of a transpose layer
and the output of one of the previous convolution layers with the same
dimension; convolutional layers and max-pooling operations are performed after
the creation of these layers to increase the resolution of the output
$Y_{s_{i}}(m_{j},n_{j})$. The structure of the mJ-Net is inspired by the U-Net
and the V-Net architectures, and has some similarities with these structures
and their layers. Nevertheless, while the U-Net and the V-Net preserve the
dimension of the input and output (respectively 2D for the U-Net and 3D for
the V-Net), the mJ-Net takes a 3D volume as input and produces a 2D image as
output.

##  4\. Experiments and Results

Due to the limited amount of training data, a leave-one-patient-out cross-
validation is performed. Also, the images are tiled into $(T\times 16\times
16)$ regions with an overlapping sliding window technique. The “core” class is
the most important infarcted region to detect, and the number of tiles in this
region is a very small portion of the entire dataset. Therefore a series of
data augmentation techniques were implemented to extend the number of core
class tiles and improve the balance between the classes during training. These
techniques consist of rotation or mirroring of the time-series tiles only if
the corresponding ground truth section is classified as a “core” area. The
total number of volumes $\tilde{V}_{s_{i}}(t,m_{j},n_{j})$ generated after
augmentation is 1,086,030, divided between the following classes: 28.3%
background, 54.9% brain, 11.7% penumbra and 5.1% core.

During the learning of the networks, two different optimizer functions were
tested: the Adaptive moment estimation (Adam) (Kingma and Ba, 2014), and the
stochastic gradient descent (SGD) optimizer function. The best results were
achieved with the SGD optimizer using a learning rate sets to 0.01 and the
Nesterov momentum sets to 0.9 (Nesterov, 2013). The statistical results are
shown in Table 3.

All experiments are performed with a GPU Tesla V100-PCIE (32GB) and 100 epochs
per training with early stopping if the loss does not decrease for a fixed
number of epochs (10). Keras 2.3.1 was used for the implementation detail of
the networks, with Tensorflow 1.14 as backend. For all the models, leave-one-
patient-out cross-validation was performed, using 9 patients (including
$\approx$4000 images) for training, and testing on the last patient (about
360-660 images) in each fold. The training of the mJ-Net, for each fold in the
leave-one-patient-out cross-validation, needs approximately 22 hours on
average to complete. Less than 20 seconds were necessary, on average, to
generate a brain slice image for a testing patient.

We choose a modified version of the Dice coefficient metric for the loss
function, called the soft Dice coefficient. This formulation was proven to
generate better outcomes for a dataset with unbalanced classes without the
need to assign weights to samples to balance the dataset (Milletari et al.,
2016). For each predicted tile $j$, the loss function is defined as:

$\text{loss}_{j}=1-\frac{2\sum_{\forall(m_{j},n_{j})}|Y_{GT_{i}}(m_{j},n_{j})Y_{s_{i}}(m_{j},n_{j})|+\epsilon}{\sum_{\forall(m_{j},n_{j})}(Y_{GT_{i}}(m_{j},n_{j}))^{2}+\sum_{\forall(m_{j},n_{j})}(Y_{s_{i}}(m_{j},n_{j}))^{2}+\epsilon}$

where $\epsilon>0$ is a small number for numerical stability. Subsequentially,
the cost function is represented as:

$\text{cost}_{\text{batch}}=\sum_{j\in\text{batch}}\text{loss}_{j}$

The best classifier architectures, among the three that were pre-tested, was
Arch_1, as shown in Table 1. Arch_1 yield the best accuracy, precision, F1
score, and recall on average; it also presented a considerable fast
computational time for each training epoch.

Fig. 7(a) shows three ground truth examples of the manual annotated brain
slices $s_{i}$. The predicted outcomes of the same brain slices $s_{i}$ with
the best mJ-Net after post-processing steps are shown in Fig. 7(b). At the
bottom of Fig. 4 (marked with 5), an example of the generated brain slice
$s_{i}$ with Arch_1, after the union of the various tiles; the ground truth
image is displayed on the left part of Fig. 7(a).

Table 3. Comparison of the mJ-Net results with literature methods. The input
column shows specific parametric maps with the corresponding threshold used as
input to extrapolate the infarcted areas, except our approach which uses in
input the entire 4D CTP images.
*Kasasbeh’ results were generetad not with our dataset, however a comparison is still applicable since the NIHSS score of both dataset is almost identical.

Method | Patients |  Infarcted Area |  Dice coef. | Sens. | Spec. | Prec. | Acc. | AUC | Input
---|---|---|---|---|---|---|---|---|---
mJ-Net (SGD) | 10 | Penumbra | 0.78 | 0.86 | 0.93 | 0.72 | 0.95 | 0.97 |  4D CTP
Core | 0.53 | 0.75 | 0.98 | 0.41 | 0.99 | 0.94
Wintermark et al. (Wintermark et al., 2006) | 10 | Core | 0.14 | 0.72 | 0.90 | 0.08 | 0.90 | 0.86 | CBV$<$33%
Campbell et al. (Campbell et al., 2012) | 10 | Penumbra | 0.40 | 0.30 | 0.96 | 0.62 | 0.85 | 0.71 |  $T_{\text{Max}}>$6s
Core | 0.05 | 0.35 | 0.86 | 0.03 | 0.85 | 0.30 | CBF$<$31% & TTP$>$4s
Cereda et al. (Cereda et al., 2016) | 10 | Penumbra | 0.34 | 0.32 | 0.88 | 0.35 | 0.79 | 0.69 |  $T_{\text{Max}}>$4s
Core | 0.07 | 0.85 | 0.74 | 0.04 | 0.74 | 0.87 | CBF$<$38%
Ma et al. (Ma et al., 2019), Lin et al. (Lin et al., 2014) | 10 | Penumbra | 0.37 | 0.26 | 0.97 | 0.65 | 0.85 | 0.69 |  $T_{\text{Max}}>$6s
Core | 0.09 | 0.80 | 0.81 | 0.05 | 0.80 | 0.87 | CBF$<$30%
Kasasbeh (Kasasbeh et al., 2019) * | 128 | Core | 0.48 | 0.91 | 0.65 | N.A. | N.A. | 0.87 |  CBF, CBV, TTP, $T_{\text{Max}}$ & clinical data

We excluded “background” predictions from the calculation of the statistical
results to balance the results only with the classes inside the skull. Since
the pixel values are in a [0, 255] domain, we decided to categorize them into
four classes. Based on the histogram of one output image (Fig. 7(c)) we choose
for the “brain” class the domain $[0,60)$, $[60,135)$ for “penumbra”,
$[135,234)$ for “core” and the remaining values are “background”. We noticed
from the histogram that most output values are centered around the target
values $\pm 15$, so the threshold is not very sensitive to small variations
outside these limits: future work is needed to find optimal thresholds on the
training dataset.

The statistical results for the best mJ-Net architectures are compared with
other reported methods in the literature, as showed in Table 3. All the
thresholding methods (Campbell et al., 2012; Cereda et al., 2016; Ma et al.,
2019; Lin et al., 2014; Wintermark et al., 2006) were tested on our dataset,
excluding the “background” from the calculation, to have comparable results.
We also compare our results with the reported results of the recent CNN-based
method by Kasasbeh et al. (Kasasbeh et al., 2019), however, we could not test
it on our dataset because there are no publicly available weights for the
proposed neural network. Nevertheless, since the median NIHSS score was 15 for
the 128 patients analyzed in (Kasasbeh et al., 2019), and this is
approximately the same as the NIHSS score of our dataset, a comparison is
still applicable.

The optimal mJ-Net, based on the SGD optimizer function, produces a Dice
coefficient score of 0.78 and 0.53, for penumbra and core respectively. The
area under the receiver operating characteristic curve (AUC), calculated with
an expanding upper and lower bound starting from the target value of each
class, is 0.97 and 0.94, corresponding to penumbra and core. A comparison with
the AUC of the various reference methods is given in Fig. 8(a) for the
penumbra area and in Fig. 8(b) for the core region. Our proposed method, with
4D CTP input, preprocessing and mJ-Net for multiclass segmentation, achieved
the best results both for penumbra and core.

##  5\. Discussion

(a)AUC for penumbra region

(b)AUC for core region

Figure 8. Comparison of AUC both for (a) penumbra and (b) core regions between
the literature methods and the mJ-Net model.

This study aimed to investigate the possibility of detecting the infarcted
regions of an ischemic stroke based solely on the 4D CTP data using the
proposed neural network. The performed experiments proved the validity of this
network compared with both the referred methods that use thresholding
approaches (Wintermark et al., 2006; Campbell et al., 2012; Cereda et al.,
2016; Ma et al., 2019; Lin et al., 2014) and the semi-automatic method that
uses the parametric maps as input (Kasasbeh et al., 2019). The results are
promising for both infarcted core and penumbra regions; this solution provides
an alternative way for detecting these regions without using thresholding
values and it could be used by medical doctors as a support instrument
together with the parametric maps.

To our knowledge, this is the first study that tries to segment penumbra and
core using the 4D set of CTP scans. The results achieved by the mJ-Net provide
the foundations for this objective but further research is necessary to
validate the results of the network.

The bottom part of Fig. 4 (marked with (5)) presents one of the outputs of the
Arch_1: the areas inside the brain are classified in a pixelated way; this is
not helpful for medical decisions during the treatment of a patient. For this
reason, mJ-Net was implemented to improve Arch_1 to overcome the approximation
in the regions’ segmentation. Arch_1 was the model that yields the best
results among the three classifiers; it also had the best ratio between the
predictions and computational time. The other two architectures presented
worse outcomes in comparison with Arch_1; moreover, they have some
disadvantages. Arch_2 did not convolve the entire third dimension (time) with
its kernel, thus its results were not as accurate as of the other two
classifiers. The number of parameters involved in Arch_3 is almost four times
larger than the one used in Arch_2: this increases significantly the
computational time of Arch_3 with improvements in the results of some patients
but not significantly close to Arch_1 (see Table 1).

Our proposed best network displays better results in the Dice coefficient,
precision, accuracy, and AUC (Table 3) both for penumbra and core compared
with the methods analyzed. Also, the mJ-Net shows high outcomes for
sensitivity (“penumbra”) and specificity (“core”) even if these are not the
best results for both the infarcted regions, since Kasasbeh et al. (Kasasbeh
et al., 2019) presented a better outcome in the sensitivity for the “core”
region, while the thresholding values both proposed by Ma et al. (Ma et al.,
2019) and Lin et al. (Lin et al., 2014) achieved a better specificity for the
“penumbra” area.

The dataset presents some limitations. First, due to the small number of
patients involved; however, the leave-one-patient-out cross-validation
performed over the networks helped to achieve consistent results. Second, the
manual annotations of the ischemic regions are not necessarily providing a
perfect ground truth. The annotations are enclosing the most important regions
of the infarcted areas but might leave out small regions of core spread in the
penumbra, etc. However, the annotations are made by a medical expert in the
field and reflect how the clinicians are working in a true setting.

##  6\. Conclusion & Future Work

To our knowledge, this is the first study to use the entire 4D set of CTP
brain slices over the injection period to train a 3D CNN to segment the
infarcted regions, both core and penumbra, in patients affected by an acute
ischemic stroke in a fully automated method.

The present study serves as a proof-of-concept, showing the possibilities of
this technology, and it displays promising results with the cross-validated
result: average Dice coefficient score of 0.78 and 0.53, and an AUC of 0.97
and 0.94 for penumbra and core respectively. It shows substantially improved
results compared to the different methods combining thresholding of the
different parametric maps. Further research with a larger dataset is required
to assert the validity of the proposed architecture and the achieved results.
Another possible future work might be a focus on the evaluation of a more
precise ground truth based on all the parametric maps generated during the CTP
in combination with medical expert input. Furthermore, another interesting
study for the future should be the test of the proposed method with a dataset
composed of patients affected by small or large vessel occlusions to
understand how the architecture can segment the infarct regions in a more
general dataset.

## References

  * (1)
  * Barman et al. (2019) Arko Barman, Mehmet E Inam, Songmi Lee, Sean Savitz, Sunil Sheth, and Luca Giancardo. 2019\.  Determining Ischemic Stroke From CT-Angiography Imaging Using Symmetry-Sensitive Convolutional Networks. In _2019 IEEE 16th International Symposium on Biomedical Imaging (ISBI 2019)_. IEEE, 1873–1877.
  * Campbell et al. (2012) Bruce CV Campbell, Søren Christensen, Christopher R Levi, Patricia M Desmond, Geoffrey A Donnan, Stephen M Davis, and Mark W Parsons. 2012\.  Comparison of computed tomography perfusion and magnetic resonance imaging perfusion-diffusion mismatch in ischemic stroke.  _Stroke_ 43, 10 (2012), 2648–2653.
  * Campbell et al. (2013) Bruce CV Campbell, Louise Weir, Patricia M Desmond, Hans TH Tu, Peter J Hand, Bernard Yan, Geoffrey A Donnan, Mark W Parsons, and Stephen M Davis. 2013.  CT perfusion improves diagnostic accuracy and confidence in acute ischaemic stroke.  _J Neurol Neurosurg Psychiatry_ 84, 6 (2013), 613–618.
  * Cereda et al. (2016) Carlo W Cereda, Søren Christensen, Bruce CV Campbell, Nishant K Mishra, Michael Mlynash, Christopher Levi, Matus Straka, Max Wintermark, Roland Bammer, Gregory W Albers, et al. 2016\.  A benchmarking tool to evaluate computer tomography perfusion infarct core predictions against a DWI standard.  _Journal of Cerebral Blood Flow & Metabolism_ 36, 10 (2016), 1780–1789.
  * Çiçek et al. (2016) Özgün Çiçek, Ahmed Abdulkadir, Soeren S Lienkamp, Thomas Brox, and Olaf Ronneberger. 2016\.  3D U-Net: learning dense volumetric segmentation from sparse annotation. In _International conference on medical image computing and computer-assisted intervention_. Springer, 424–432.
  * Committee et al. (2008) European Stroke Organisation (ESO) Executive Committee, ESO Writing Committee, et al. 2008\.  Guidelines for management of ischaemic stroke and transient ischaemic attack 2008.  _Cerebrovascular diseases_ 25, 5 (2008), 457–507.
  * Goshtasby (2005) Arthur Ardeshir Goshtasby. 2005\.  _2-D and 3-D image registration: for medical, remote sensing, and industrial applications_.  John Wiley & Sons.
  * Hacke et al. (2008) Werner Hacke, Markku Kaste, Erich Bluhmki, Miroslav Brozman, Antoni Dávalos, Donata Guidetti, Vincent Larrue, Kennedy R Lees, Zakaria Medeghri, Thomas Machnig, et al. 2008\.  Thrombolysis with alteplase 3 to 4.5 hours after acute ischemic stroke.  _New England journal of medicine_ 359, 13 (2008), 1317–1329.
  * Kamnitsas et al. (2017) Konstantinos Kamnitsas, Christian Ledig, Virginia FJ Newcombe, Joanna P Simpson, Andrew D Kane, David K Menon, Daniel Rueckert, and Ben Glocker. 2017.  Efficient multi-scale 3D CNN with fully connected CRF for accurate brain lesion segmentation.  _Medical image analysis_ 36 (2017), 61–78.
  * Kasasbeh et al. (2019) Aimen S Kasasbeh, Søren Christensen, Mark W Parsons, Bruce Campbell, Gregory W Albers, and Maarten G Lansberg. 2019.  Artificial Neural Network Computer Tomography Perfusion Prediction of Ischemic Core.  _Stroke_ 50, 6 (2019), 1578–1581.
  * Kim (1997) Yeong-Taeg Kim. 1997\.  Contrast enhancement using brightness preserving bi-histogram equalization.  _IEEE transactions on Consumer Electronics_ 43, 1 (1997), 1–8.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba. 2014.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * Kurz et al. (2016) KD Kurz, G Ringstad, A Odland, R Advani, E Farbu, and MW Kurz. 2016.  Radiological imaging in acute ischaemic stroke.  _European journal of neurology_ 23 (2016), 8–17.
  * Lin et al. (2014) Longting Lin, Andrew Bivard, Christopher R Levi, and Mark W Parsons. 2014.  Comparison of computed tomographic and magnetic resonance perfusion measurements in acute ischemic stroke: back-to-back quantitative analysis.  _Stroke_ 45, 6 (2014), 1727–1732.
  * Lucas et al. (2018) Christian Lucas, André Kemmling, Amir Madany Mamlouk, and Mattias P Heinrich. 2018.  Multi-scale neural network for automatic segmentation of ischemic strokes on acute perfusion images. In _2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018)_. IEEE, 1118–1121.
  * Ma et al. (2019) Henry Ma, Bruce CV Campbell, Mark W Parsons, Leonid Churilov, Christopher R Levi, Chung Hsu, Timothy J Kleinig, Tissa Wijeratne, Sami Curtze, Helen M Dewey, et al. 2019\.  Thrombolysis guided by perfusion imaging up to 9 hours after onset of stroke.  _New England Journal of Medicine_ 380, 19 (2019), 1795–1803.
  * Milletari et al. (2016) Fausto Milletari, Nassir Navab, and Seyed-Ahmad Ahmadi. 2016\.  V-net: Fully convolutional neural networks for volumetric medical image segmentation. In _2016 Fourth International Conference on 3D Vision (3DV)_. IEEE, 565–571.
  * Nesterov (2013) Yurii Nesterov. 2013\.  _Introductory lectures on convex optimization: A basic course_. Vol. 87.  Springer Science & Business Media.
  * Öman et al. (2019) Olli Öman, Teemu Mäkelä, Eero Salli, Sauli Savolainen, and Marko Kangasniemi. 2019\.  3D convolutional neural networks applied to CT angiography in the detection of acute ischemic stroke.  _European radiology experimental_ 3, 1 (2019), 8\.
  * Ronneberger et al. (2015) Olaf Ronneberger, Philipp Fischer, and Thomas Brox. 2015\.  U-net: Convolutional networks for biomedical image segmentation. In _International Conference on Medical image computing and computer-assisted intervention_. Springer, 234–241.
  * Saver (2006) Jeffrey L Saver. 2006\.  Time is brain-quantified.  _Stroke_ 37, 1 (2006), 263–266.
  * Sheth et al. (2019) Sunil A Sheth, Victor Lopez-Rivera, Arko Barman, James C Grotta, Albert J Yoo, Songmi Lee, Mehmet E Inam, Sean I Savitz, and Luca Giancardo. 2019.  Machine Learning–Enabled Automated Determination of Acute Ischemic Core From Computed Tomography Angiography.  _Stroke_ (2019), STROKEAHA–119.
  * Vargas et al. (2019) Jan Vargas, Alejandro Spiotta, and Arindram Rano Chatterjee. 2019\.  Initial Experiences with Artificial Neural Networks in the Detection of Computed Tomography Perfusion Deficits.  _World neurosurgery_ 124 (2019), e10–e16.
  * Vincent and Soille (1991) Luc Vincent and Pierre Soille. 1991.  Watersheds in digital spaces: an efficient algorithm based on immersion simulations.  _IEEE Transactions on Pattern Analysis & Machine Intelligence_ 6 (1991), 583–598.
  * Wang et al. (2016) Haidong Wang, Mohsen Naghavi, Christine Allen, Ryan M Barber, Zulfiqar A Bhutta, Austin Carter, Daniel C Casey, Fiona J Charlson, Alan Zian Chen, Matthew M Coates, et al. 2016\.  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 (2016), 1459–1544.
  * Wintermark et al. (2006) Max Wintermark, Adam E Flanders, Birgitta Velthuis, Reto Meuli, Maarten Van Leeuwen, Dorit Goldsher, Carissa Pineda, Joaquin Serena, Irene van der Schaaf, Annet Waaijer, et al. 2006\.  Perfusion-CT assessment of infarct core and penumbra: receiver operating characteristic curve analysis in 130 patients suspected of acute hemispheric stroke.  _Stroke_ 37, 4 (2006), 979–985.