# A Neutrino Detector Design for Safeguarding Small Modular Reactors

Emma Houston<EMAIL_ADDRESS>Sandra Bogetic  University of Tennessee,
Knoxville.  Oluwatomi Akindele  Marc Bergevin  Adam Bernstein  Steven Dazley
Lawrence Livermore National Laboratory

###### Abstract

Nuclear reactors have long been a favored source for antineutrino measurements
for estimates of power and burnup. With appropriate detector parameters and
background rejection, an estimate of the reactor power can be derived from the
measured antineutrino event rate. Antineutrino detectors are potentially
attractive as a safeguards technology that can monitor reactor operations and
thermal power from a distance. Advanced reactors have diverse features that
may present challenges for current safeguards methods. By comparison, neutrino
detectors offer complementary features, including a remote, continuous,
unattended, and near-real-time monitoring capability, that may make them
useful for safeguarding certain classes of advanced reactors. This study
investigates the minimum depth and size of an antineutrino detector for a SMR
to meet safeguards needs for advanced reactors. Extrapolating performance from
several prior reactor antineutrino experiments, this study uses an analytical
approach to develop a possible design for a remote antineutrino-based
monitoring device.

††preprint: APS/123-QED

##  I Introduction

The highly penetrating nature of neutrinos and their high rate of production
in reactors, makes them attractive as a potential reactor safeguards
monitoring technology. Several applications have been identified as potential
avenues for further exploration to determine if neutrino detectors can be
deployed for safeguard measurements [1]. A working group between the
International Atomic Energy Agency (IAEA) and antineutrinos experts determined
that antineutrino detectors may offer a complementary solution to the current
safeguards structure [1]. The implementation of neutrino detectors for
safeguarding advanced reactors was also a use case finding in the NuTools
Study. This assessment was based on expert opinions, including advanced
reactor developers, on applications of neutrino detection [2]. Several
advanced reactor types, including Small Modular Reactors (SMRs) and Molten
Salt Reactors (MSRs) were specifically identified as an opportunity for
potential neutrino detector applications. Previous studies on antineutrino
detection sensitivities to diversion of thorium-fueled MSRs have been
performed, while investigations into SMRs are less abundant [3]. When choosing
a detector size, economics and construction feasibility should be considered.
With this consideration, the ultimate goal of this study is to identify the
smallest and shallowest possible detector to meet IAEA detection goals.

##  II Background

###  II.1 Neutrino Detectors at Power Reactors

Nuclear beta decay is the primary driver for neutrino emissions in nuclear
reactors. The fission of 235U and 239Pu produces neutron-rich fission products
that are likely to beta decay. On average, these isotopes undergo $\sim 6$
beta decays per fission before reaching stability [4]. The neutrino detection
rate is proportional to the reactor’s thermal output. With a sufficient number
of events, the energy spectrum and number of emitted neutrinos can be used to
estimate the fissile content of reactor fuel [4]. However, this study will
solely investigate reactor power and operation as a function of neutrino rate,
due to the relatively modest number of events that can be acquired for
relevant standoff distances and dwell times.

Several previous experiments have used these neutrino rates from power
reactors to determine neutrino properties [5, 6, 7, 8]. These neutrino
experiments rely on inverse beta decay (IBD) interactions in gadolinium-loaded
liquid scintillator to detect neutrinos.

${\bar{\nu}_{e}}+p\longrightarrow e^{+}+n$

In this paper we will rely on previous gadolinium-loaded liquid scintillator
reactor experiments to analytically determine the feasibility of a neutrino
power monitoring detector for a generic Light Water SMR, similar to NuScale.

###  II.2 Challenges in SMR Safeguards

Nuclear reactors provide approximately 20% of the United States electrical
output [9]. Moving to carbon-free sources of electricity requires nuclear
energy but the last decade has shown the challenges behind our current
generation of nuclear reactors in terms of cost and timeliness of
construction. Additionally, the current generation of large scale nuclear
plants are well-suited for specific markets but not necessarily adaptable for
all facets of the electric grid. Several novel designs have been proposed that
differ in fuel type, coolant type, and size. The SMR design has adaptable
electrical output, dependent on the number of cores, which can be tailored to
the needs of the grid at a particular location [10]. This offers the
availability to put these reactors in more remote locations.

SMRs present both advantages and challenges to present day safeguard
techniques. On one hand, SMRs generally have smaller cores making them a less
attractive target for material theft or diversion. However, the anticipated
advantage of remote deployment is a challenge to performing inspections [11].
The flexibility in power generation make SMRs ideal for remote locations or
developing infrastructure, but also the number of cores increases the
frequency of inspections and the burden on IAEA inspectors[12].

The introduction of neutrino detectors as a remote monitoring tool for SMRs
may help address some of these challenges. The IAEA has used remote monitoring
technologies for the past few decades for other containment and surveillance
(C$\&$S) equipment [13]. As neutrino detectors advance through technology
readiness levels with ongoing experiments and are capable of remote operation,
neutrino detectors show promise as an additional monitoring tool in the
current C$\&$S toolbox [14].

##  III Diversion and Misuse Scenarios

To design a neutrino detector with the goal of advanced reactor safeguards, it
is useful to formulate diversion and misuse scenarios where its use could
provide potential indicators of diversion. Two hypothetical diversion and
misuse scenarios are defined below. The requirement to reduce total site
thermal power is a common element in both scenarios. The relevant neutrino
detector signature will be a signal rate decrease due to a reduction in power.

###  III.1 Scenario 1: Recovery from Unexpected Gaps in C$\&$S Measurements

The increased number of cores at SMR sites, and increased construction of SMRs
due to power demands, may increase the burden on the IAEA in-terms of
inspector-days needed for inventory and information verification, and/or re-
establishing Continuity of Knowledge following an anomaly, These new
circumstances may increase the opportunities for the failure of standards
safeguards measures.

Current IAEA safeguards at reactors rely primarily on C$\&$S in the form of
optical surveillance or seals. These elements can fail in a variety of ways,
including loss of power to surveillance cameras, degradation due to
environmental or radiative effects, and other potentially purposeful attacks.
Additionally, some C$\&$S technologies rely on in-person inspection of seals
that retrospectively indicate an item has been tampered with, and thus don’t
provide real-time monitoring of diversion. Further, many current C$\&$S
safeguards methods must be performed for the inspection of fresh fuel in each
core at each refueling. This increases the number of inspections required to
safeguard a multi-core SMR facility.

In our first concept of operations, the neutrino detector would monitor the
aggregate site reactor power, flagging anomalies such as the temporary
shutdown of a single core, or power reductions across several cores. As the
IAEA considers designs for SMR-specific safeguards concepts, a neutrino-based
component offers a complementary capability for continuous online monitoring
of aggregate reactor operations, and may help alleviate some of the above
problems. Currently, inspections are chosen at random, or based on the other
C$\&$S mechanisms in place. The indications and irregularities in the neutrino
signal can potentially aid in the decision to pursue unplanned inspections or
even routine inspection. Increased deployment has been identified as a
challenge behind SMR safeguards [12]. Remote neutrino-based power monitoring
may also help inform decisions about which cores to inspect and prioritize
within a larger number of operating reactors.

###  III.2 Scenario 2: Reducing the Reactor Power To Improve Plutonium Quality

At the beginning of cycle, the quality of the 239Pu to other Pu isotopes ratio
is higher than at the end of cycle. This is because the beginning of cycle
produces 239Pu through neutron capture on 238U. However, as 239 Pu is created,
it is converted to 240Pu through neutron capture. Operating at a lower power
can slow the production of the higher Pu isotopes such as 240Pu. Higher ratios
of 239Pu to other Pu isotopes would be desired by a proliferator because the
other Pu isotopes can be challenging to handle to design a device. 240Pu has a
large amount of both decay heat and spontaneous fission neutrons, 241Pu decays
to 241Am, which is very radioactive. Additionally, 242Pu has a high rate of
spontaneous neutrons and a large critical mass [15]. For these reasons, the
ratio of 239Pu to other isotopes is a key factor in determine the
attractiveness of the material to a proliferator. Most current safeguard
approaches for reactors use a combination of C$\&$S and inspections. As a
consequence, the production of weapons-grade plutonium through the cycle can
remain undetected between inspections. In addition to providing indications of
potential misuse, neutrino detectors may be used to provide redundancy to
current C$\&$S and thereby help to deter misuse in these instances.

Figure 1:  shows the decrease in 239Pu compared to other fission fractions
within the cycle using the approximation for fission fractions shown in [16]
for a 3$wt\%$ LEU fueled PWR. The decrease in 239Pu to total Pu ratio is
caused by neutron capture and the conversion of 239Pu to 240Pu. While this
graph does not include the relative amounts of 240Pu in the cycle, the ratio
of 239Pu decreases compared to the other Pu isotopes. Further comparisons of
the simplified fission fraction approximations to reactor data are shown in
the following [17].

Neutrino detectors provide a method to monitor reactor power levels from a
standoff distance outside of the site perimeter. Neutrinos are impossible to
shield and thus adversaries cannot hide the signals arising from reactor power
operation. These neutrino measurements may provide valuable information about
the period over which power has been decreased, which may give insight into
the quality of plutonium that could be produced from a reactor.

In both scenarios, the goal of the neutrino detector is to detect the
reduction in a timely manner. This reduction in neutrino signal is indicative
of the reduction in power of the aggregate site power.

##  IV Signal and Background Determination

###  IV.1 Reactor Signal

The detector mass and its distance from the reactor core affect the signal
rate. Here we vary these detector parameters to identify a decrease in reactor
power consistent with a single reactor shutdown from a six core generic SMR
configuration. Considering the goal of this study is to analyze a detector’s
safeguard capabilities from outside of the site perimeter, detector parameters
to meet current IAEA safeguards standards are investigated. The Exclusion Area
Boundary (EAB) and Low Population Zone (LPZ) are found using dose estimates,
among other factors, to require limiting distances from the reactor to the
EABs and LPZs [18]. There is significant cost reduction for SMRs due to their
smaller EABs [19]. The smaller EABs in turn makes SMRs attractive sites for
neutrino detectors, which can be much closer to the reactor and remain outside
the perimeter compared to traditional Light Water Reactors (LWRs). We define
the distance from the core to the EAB as the minimum distance of the detector
from the reactor to be considered non-intrusive.

The antineutrino signal can be found using Eq. (1). It is assumed the reactor
power is constant thus the number of fissions over time does not change. A
more robust approximation would include calculations with regard to burnup
over time. However, this study is intended to reflect the averaged detection
rates over a whole cycle.

$N_{det}(E_{\nu},t)=\frac{\varepsilon(E_{\nu})}{4\pi
L^{2}}\phi_{\nu}(E_{\nu},t)\sigma(E_{\nu})N_{T}P_{surv}(E_{\nu},L)$ (1)

In Eq. 1, $N_{det}$ is the number of neutrinos detected per second [20]. The
estimated detection efficiency, $\varepsilon$ is $80\%$, based on recent
experiments [6], [21], and was treated as the average over all neutrino
energies. The estimated distance from the reactor to the detector is L. The
flux averaged cross section, $\sigma$, is found by [22]. This study assumed an
operating thermal power of 160 MWt per core [23] with 6 cores operating during
normal operation. The oscillation dependent survival probability, $P_{surv}$,
is given by Eq. 65.9 in the Review of Particle Physics [24]. The energy-
dependent neutrino flux is given a fifth order polynomial with coefficients
for 235U, 238U, 239Pu, and 241Pu in the following references [25], [26]. The
fission fractions and energy per fission from the following source was used
[27]. $N_{T}$ is the number of free protons within the detector and is
directly proportional to the detector mass. Similar to previous antineutrino
detector experiments, we will assume the scintillator will be Gd-loaded linear
alkyl benzene (LAB) and thus has a free proton per molecule of 30 and a
density of 0.86 g/ml [28]. This is used to calculate the number of free
protons per ton.

###  IV.2 Background Interpolation

This study considers three background types: 9Li-8He production, fast-neutron
production, and accidentals.

####  IV.2.1 Muon-Induced Backgrounds

Both cosmogenic radionuclide production (9Li and 8He), and fast-neutron
production are induced by atmospheric muons. When muons interact with the 12C
inside the scintillator, 9Li and 8He can be produced. Subsequent $\beta-n$
decays of these isotopes mimic antineutrino-induced IBD interactions [29, 30].
Fast neutrons can be produced when muons interact with nuclei in the rock
surrounding the detector. Some of these fast neutrons then produce recoil
protons in the detection medium, creating a prompt signal, then thermalize
within the detector and capture to create a delayed signal. The combination of
these two correlated signals and their temporal structure again mimics the IBD
interaction[30]. To reduce these backgrounds, the detector will require some
overburden to reduce the muon flux. However, as with detector size, detector
depth has a strong impact on detector construction and operational costs.

To approximate muon-induced backgrounds, several reactor neutrino experiments
were benchmarked and an analytical model was used to estimate background
rates. Table1 displays the data used to approximate the background from 9Li,
8He and fast-neutrons.

Table 1:  shows Li-He, fast-neutron, and accidental backgrounds from previous neutrino detector experiments. These values are used to estimate the production of muogenic and accidental backgrounds.  Experiment | Depth mwe | Li-He (events/day) | Fast Neutrons (events/day) | Accidentals (events/day)
---|---|---|---|---
RENO | 120 |  $12.45\pm 5.93$ [31] |  $5.0\pm 0.13$ [31] |  $4.3\pm 0.06$ [32]
Double Chooz | 120 |  $5.01\pm 1.43$[33] |  $3.42\pm 0.23$ [33] |  $0.344\pm 0.002$ [33]
Daya Bay | 250 |  $3.1\pm 1.6$ [31] |  $0.84\pm 0.28$ [31] |  $9.85\pm 0.06$ [31]
Daya Bay | 265 |  $1.8\pm 1.1$ [31] |  $0.74\pm 0.44$ [31] |  $7.67\pm 0.05$ [31]
Double Chooz | 300 |  $0.97_{-0.16}^{+0.41}$ [33] |  $0.586\pm 0.061$ [33] |  $0.106\pm 0.002$ [33]
RENO | 450 |  $2.59\pm 0.75$ [31] |  $0.97\pm 0.06$ [31] |  $0.68\pm 0.03$ [32]
Daya Bay | 860 |  $0.16\pm 0.11$ [31] |  $0.04\pm 0.04$ [31] |  $3.25\pm 0.03$ [31]

The background rates in the theoretical detector are then assumed to scale
with the volume of the detector in the case of 9Li, 8He and surface area in
the case of fast neutron backgrounds. The muon track lengths in the detector
are proportional to the detector volume,increasing the number of number of
potential 12C atoms to interact with. The rates for fast neutrons are scaled
to the surface area of the detector as the surrounding rock scales directly
with this type of background. The detector dimensions used to determine the
surface area and volume are shown in Table 2. Daya Bay has multiple detectors
at the same depth, each with their own background measurements. In these cases
the average background was used, rather than treating all detectors as
individual data points. For the Li-He background fit, the Daya Bay backgrounds
for 250 mwe and 265 mwe were averaged to avoid over fitting. (This was not
needed for fast-neutron background fitting as the fast neutrons backgrounds at
these two depths were very similar.) Lastly, RENO has more recent lower
backgrounds measurements that are the result of more stringent multiplicity
measurements [34]. However, the larger backgrounds measurements were chosen in
this approximation for conservatism.

Table 2: The detector and reactor characteristics of each experiment are listed by increasing depth. The radi, height, baselines, and signal are averaged amongst detectors from the same experiment in the same hall.  Experiment | Depth | Mass | Radius | Height | Site Power | Baseline | Efficiency | Signal
---|---|---|---|---|---|---|---|---
| mwe | tons | cm | cm | GWt | m | $\%$ | Events/Day
RENO | 120 | 16 [31] | 137.5 [35] | 315.0 [35] | 16.52 [31] | 408.56 [31] |  $64.70\%$ [32] |  $779.05$ [32]
Double Chooz | 120 | 9 [36] | 115.0 [37] | 245.8 [37] | 8.5 [33] | 400 [38] | - | 293.4 [33]
Daya Bay | 250 | 20 [21] | 155 [39] | 308.2 [39] | 17.4 [40] | 512 [41] |  $78.8\%$ [31] | 716.0 [31]
Daya Bay | 265 | 20 [21] | 155 [39] | 308.4 [39] | 17.4 [40] | 561 [41] |  $78.8\%$ [31] | 532.29 [31]
Double Chooz | 300 | 9 [36] | 115.0 [37] | 245.8 [37] | 8.5 [33] | 1050 [38] |  $78.0\%$ [42]111Double Chooz’s selection efficiency and Gd-capture fraction are listed separately in the following reference [42]. It is assumed these can be treated as follows to obtain the total efficiency: (91.5$\%\times 85.3\%=78.0\%$). |  $39.0$ [33]
RENO | 450 | 16 [31] | 137.5 [35] | 315.0 [35] | 16.52 [31] | 1443.99 [31] |  $74.50\%$ [32] |  $72.78$ [32]
Daya Bay | 860 | 20 [21] | 155 [39] | 308.5 [39] | 17.4 [40] | 1579 [41] |  $78.8\%$ [31] | 70.7 [31]

Daya Bay, Double Chooz, and RENO muon-induced background data were each fit to
an exponential function shown by Eq. 2. Each of these detectors have unique
selection criteria to optimize the signal to background rates. These
exponential fits were averaged for both types of correlated backgrounds.

$y=me^{-tx}$ (2)

Figure 2: The equations of best fit for these background measurements are:
$y=7.89\times 10^{-7}e^{-5.90\times 10^{-3}x}\frac{events}{day-cm^{3}}$ for
the Li-He induced background, where x is depth in mwe . Backgrounds for each
set up detectors are fit to Eq. 2 and then averaged. Fig. 2 show this
functions that are then sampled from to approximate background. Figure 3: The
equations of best fit for these background measurements are $y=2.10\times
10^{-5}e^{-6.50\times 10^{-3}x}\frac{events}{day-cm^{2}}$ for the fast-neutron
induced background, where x is depth in mwe . Backgrounds for each set up
detectors are fit to Eq. 2 and then averaged. Fig. 3 show this function that
are then sampled from to approximate background.

There are several considerations when extrapolating this data to approximate
background rates at shallower depths. This method for background
approximations does not account for veto capabilities of the varying
detectors. Additionally, this method assumes that these backgrounds are
directly related to the attenuated muon flux through the rock overburden, and
thus Eq. 2 can be applied. Higher energy muons are found at deeper rock
overburdens because lower energy muons are more likely to interact and stop
entirely. At shallow depths, both the muon energy and angular distribution
have a significant impact on muon intensity however this study approximates
muon flux using solely depth [43].

Lastly, RENO reported a contamination of 252Cf in their detector, with a more
significant impact on the far RENO detector than near. A multiplicity
requirement helps to mitigate the effects of the 252Cf source, causing the
contamination of the near detector to be negligible and the contamination of
the far detector to amount to $0.095\pm 0.018$ events/day [44]. Shallower
depth detector experiments, such as Palo Verde at 32 mwe, may offer some
insight into muon-induced backgrounds at shallower depths [5, 45]. However,
given the significant differences in detector designs, solely monolithic
detector experiments were chosen for extrapolation.

####  IV.2.2 Accidentals

Aside from the muon-induced backgrounds, accidentals, or uncorrelated
backgrounds, are another important source of IBD-like events in a detector.
This background is not depth dependent, and can be caused by a variety of
factors including detector contamination, selection criteria, shield
thickness, and detector size. For these reasons, accidentals can not be
correlated to a single parameter, such as depth in the previous case, and must
be treated separately.

The causes for uncertainty in accidentals can be the uncertainties
radioimpurities of the detector materials or surrounding rock overburden, or
any other natural radioactivity that can be attributed to coincidence. Ways to
minimize these accidental backgrounds include shielding thickness, stringent
multiplicity measurements in space, time, and energy, and general detector
cleanliness. Due to the unpredictability of accidental backgrounds, a scaled
prediction was not created and rather the largest accidentals rate from the
detector measurements was used. Like fast neutrons, accidentals are also
proportional to the surface area surrounding a detector. The largest
accidentals rate per surface area from the detectors using Tables 1 and 2 is
from the far RENO detector at 450 mwe. The maximum number of accidentals
divided by surface area is $5.57\times 10^{-5}\frac{events}{day-cm^{2}}$. As
this analysis will further show, a surface area is found for each of the
detector masses and multiplied by this value to find a conservative estimate
for accidental backgrounds in our theoretical detector.

###  IV.3 Uncertainty Estimates

The Double Chooz detector signal absolute uncertainty, including uncertainties
in cross section, fractional fission rates, thermal power, mean energy per
fission, and distance, is propagated to be 1.7 % [6]. Neglecting the effect of
burnup has the consequence of degrading the precision of the thermal power
estimate. Therefore the thermal power component of the uncertainty was
increased from 0.5% to 5%, while keeping all of the other uncertainties
constant. The updated signal uncertainty then propagated to 5.25%.

In all events, the largest systematic uncertainty percentage from each
respective data set, cosmogenic radionuclides (9Li and 8He), fast-neutrons and
accidentals, were used when generating assumed systematic uncertainties from
background. The largest cosmogenic radionuclide systematic uncertainty
percentage was 50% from Daya Bay at 250 mwe and the largest fast-neutron
uncertainty percentage was 100 % from Daya Bay at 860 mwe shown in Table 1.
The largest uncertainty percentage was $4.8\%$ from Palo Verde. These maximum
uncertainties using previous experiment data in Table 1 will be used as a
conservative approximation of background uncertainties.

Using these systematic uncertainties from previous detector experiments, 4
trials were conducted using combinations of $\pm$ reactor signal systematic
uncertainty and $\pm$ background systematic uncertainty. Using the
combinations of Gaussian distributed total signal rates with these
systematically shifted means, there was no significant different in detection
times, false alarm rates, and missed alarm rates. Given this, consideration of
systematic uncertainties was not included in the final detection times and
alarm rates.

##  V Detection Metrics

###  V.1 IAEA Requirements for Containment and Surveillance

With established signal and background rates that are detector depth, size,
and distance dependent, we can determine how quickly we can detect a reduction
in power as discussed in the diversion scenarios. The time to detect the
reduction in power from 6 to 5 operating reactors, and the Type I and II error
rates [46] are the metrics considered for the success of the detector. The
Type I error is the frequency of which the null hypothesis is declared to be
false when it is true. Type II error is the frequency of which the null
hypothesis is declared true when it is false [46]. The Type I and II error
rates will be further described as the false and missed alarm rates,
respectively. The IAEA timeliness detection goals depend on the quality and
quantity of fissile material [46]. One objective is to reduce the time
required to detect the diversion of a significant quantity of fissile
material. Although the IAEA timeliness component of the IAEA inspections goals
may vary by reactor [46], the desired goal to observe a reduction in power was
determined to be 30 days in this study for conservatism. This should not be
used as a sole indicator of diversion, and thus 30 days would allow ample time
for further inspections. These goals govern the technical requirements for
detector depth and size. The false alarm probability, $\alpha$, is the
frequency that there is an indication that diversion has occurred when there
has been no diversion. The goal for IAEA false alarms (Type I error) is
generally small [46], and therefore $\alpha<=0.05$ is the goal used in this
study. The Type II error goal, or missed alarm rate, was set to be 10 $\%$
Therefore the targets for this study are $\alpha=0.05$ and $\beta=0.1$.

###  V.2 Simulating Signal

To identify the detection time for the reduction in power of a reactor, the
total signal from two scenarios is simulated. In the first, the expected
detector response is normally distributed around the average signal from 6
operating cores and the corresponding background at a given depth for 200
days. In the second instance, the detector response is normally distributed
around the signal from 6 operating cores and corresponding background for 100
days and then the detector response is simulated for 5 operating reactors and
the same background for the following 100 days. To create this data, the Numpy
random sampling feature was utilized [47]. Fig. 4 shows the Guassian-like
distribution of the number of signal events per day distributed data using the
Numpy package. This is for 6 reactors operating at 160 MWth each at a distance
of 244 m from the reactors. The following data and analysis is for a 5 ton
detector at 20 mwe, however this was repeated for various configurations.

When this data is sampled, especially in instances with smaller detectors,
there is a reasonable possibility that the signal randomly sampled from the
normal distribution will be negative. In these cases, the code re-samples from
the normal distribution until a realistically feasible value is chosen.
However, this is frequently not necessary as the vast majority of sampled
signals are positive values.

Figure 4:  The randomly Guassian-like distribution of signal rates for 100
days steady state reactor power and a reduced power, in a 5 ton detector at a
depth of 20 mwe. The decreased power is representative of either of the
diversion scenarios occurring. Figure 5:  shows the Gaussian distributed
signal over time for the diversion and non-diversion case. The vertical line
indicates when there is a power decrease, and marks the time at which
diversion has occurred.

Fig. 5 shows both the diversion and non-diversion simulated signal. Using
solely the neutrino signals on a per day basis, it can be difficult to
differentiate between signals in the non-diversion and diversion simulated
data. Alarming techniques are employed to be able to find how quickly after
diversion detection can occur. To lessen the inherent fluctuations in the
data, a rolling sum of the data is found for a specified ”window”. This window
can be optimized for each detector configurations, but in this case, a time
window of 10 days was used. The rolling sum of this data is then subtracted
from the expected sum (the full power signal expectation multiplied by the
window). This expected sum is found using a calibration window of 20 days,
where a full power signal is Gaussian-distributed and sampled. The mean of
this calibration period is then used as $s_{exp}$. This calibration period is
used to emulate the time before a detector is put into operation where the
expected signal is found. Eq. (3) shows this method for data analysis where
the time window is defined as $x$ in days, the expected full power signal is
$s_{exp}$, the actual signal is $s_{a}$, and $t$ is time in days.

$Diff(t,x)=x\times s_{exp}-\sum_{i=t-x}^{t}s_{a}(t)$ (3)

It is expected that in the first few weeks of operation, the expected signal
will be determined and used as a comparison to the remainder of operation. An
alarm indicates when the difference in this aggregate signal crosses 3 sigma
above the expected difference of aggregate signals.

If this alarm is greater than 30 days, the time to detection goal in this
study, then this is considered a missed alarm. If this occurs prior to the
reduction in power, then this is considered a false alarm. There can be
multiple false alarms for each trial as there can be multiple data points
above the 3-sigma threshold. There cannot be multiple missed alarms per trial
as the data can only miss the detection goal once. This is repeated for
$10^{5}$ trials and the time to detect the reduction in reactor power is
averaged, and the frequency of missed and false alarms are found. False alarms
are found on a per day, per trial basis whereas missed alarms are solely found
per trial. False alarms can be remediated by requiring sequential data to be
above the 3 sigma threshold, however this should optimized dependent on the
detector configuration.

Figure 6: shows the $Diff(t,x)$ as defined by Eq. (3) versus time in days. The
first 10 days, the specified time window in this case, are voided as there is
not sufficient previous simulated data to find the aggregate counts from the
previous days. When this difference in signal crosses the threshold for 3
sigma, an alarm is recorded.

In this example detector shown in Fig. 5, at 5 tons at a depth of 20 mwe, the
detection time is $\sim$ 10 days, $\alpha=0.006$ and $\beta=0.024$. These
values are well within the IAEA guidance for the detection of diversion. Of
course, this process would show an indication that diversion may have occurred
and should not be used as the sole tool for diversion. However, this should
encourage further inspections of a facility or be used in the prioritization
process of facility inspection. Additionally, this mechanism is complementary
to traditional material accountancy procedures to ensure fuel is not
immediately replaced following diversion.

This process was then repeated for detectors with fiducial volumes of 1 to 7
tons, and at depths of 20 mwe to 140 mwe in increments of 20. Table 3 shows
the the time to detect the reduction in power beyond the case described in
detail above. Table 4 shows the false alarm rates for these detectors and
Table 5 shows the missed alarm rates. The sample standard deviation shown in
Eq. 4 is used to find the standard deviation of the time to detection, where
$x_{i}$ is each time within a trial, $\bar{x}$ is the mean time of the trial,
and n is the number of trials ($10^{5}$).

$\sigma=\sqrt{\frac{\sum_{i=1}^{n}(x_{i}-\bar{x})^{2}}{n-1}}$ (4)

Table 3:  shows, in days, the average time over $10^{5}$ trials to detect a reduction in power from 6 to 5 operating cores 3 sigma confidence and the standard deviation found by Eq. 4 . Depth | 1 ton | 2 tons | 3 tons | 4 tons | 5 tons | 6 tons | 7 tons
---|---|---|---|---|---|---|---
20 mwe | 33.1 $\pm$ 24.9 | 23.4 $\pm$ 20.3 | 16.6 $\pm$ 15.1 | 12.4 $\pm$ 10.7 | 10.0 $\pm$ 7.5 | 8.5 $\pm$ 5.3 | 7.6 $\pm$ 4.0
40 mwe | 33.0 $\pm$ 25.0 | 23.1 $\pm$ 20.2 | 16.4 $\pm$ 15.0 | 12.3 $\pm$ 10.6 | 9.8 $\pm$ 7.2 | 8.4 $\pm$ 5.2 | 7.5 $\pm$ 3.9
60 mwe | 33.0 $\pm$ 25.0 | 23.0 $\pm$ 20.0 | 16.3 $\pm$ 14.8 | 12.2 $\pm$ 10.4 | 9.8 $\pm$ 7.1 | 8.4 $\pm$ 5.1 | 7.5 $\pm$ 3.9
80 mwe | 32.8 $\pm$ 25.0 | 22.9 $\pm$ 20.0 | 16.1 $\pm$ 14.6 | 12.1 $\pm$ 10.1 | 9.7 $\pm$ 7.1 | 8.3 $\pm$ 5.0 | 7.4 $\pm$ 3.7
100 mwe | 32.8 $\pm$ 24.9 | 22.9 $\pm$ 20.0 | 16.0 $\pm$ 14.4 | 12.0 $\pm$ 10.1 | 9.7 $\pm$ 7.0 | 8.3 $\pm$ 4.9 | 7.4 $\pm$ 3.7
120 mwe | 32.6 $\pm$ 24.9 | 22.6 $\pm$ 19.7 | 15.9 $\pm$ 14.5 | 11.9 $\pm$ 10.1 | 9.6 $\pm$ 6.9 | 8.2 $\pm$ 4.9 | 7.4 $\pm$ 3.7
140 mwe | 32.6 $\pm$ 24.9 | 22.5 $\pm$ 19.7 | 15.8 $\pm$ 14.3 | 11.8 $\pm$ 9.9 | 9.6 $\pm$ 6.8 | 8.2 $\pm$ 4.8 | 7.3 $\pm$ 3.7
Table 4:  The fractional false alarm rates show how frequently the aggregate signal crosses the alarm threshold without a simulated reduction in power. All of the following fractional false alarm rates are within the IAEA goals for C$\&$S metrics and thus not a factor when determining the minimum size and depth detector.  Depth | 1 ton | 2 tons | 3 tons | 4 tons | 5 tons | 6 tons | 7 tons
---|---|---|---|---|---|---|---
20 mwe | 0.005 | 0.006 | 0.006 | 0.006 | 0.006 | 0.006 | 0.006
40 mwe | 0.005 | 0.005 | 0.006 | 0.006 | 0.006 | 0.006 | 0.006
60 mwe | 0.005 | 0.005 | 0.006 | 0.006 | 0.006 | 0.006 | 0.006
80 mwe | 0.005 | 0.005 | 0.006 | 0.006 | 0.006 | 0.006 | 0.006
100 mwe | 0.005 | 0.006 | 0.006 | 0.006 | 0.006 | 0.006 | 0.006
120 mwe | 0.005 | 0.006 | 0.006 | 0.006 | 0.006 | 0.006 | 0.006
140 mwe | 0.005 | 0.006 | 0.006 | 0.006 | 0.006 | 0.006 | 0.006
Table 5:  The missed alarm rates are indicative of the fractional rate that the detection times are after the goal detection time (30 days), in which case the reduction in power was not detected soon enough.  Depth | 1 ton | 2 tons | 3 tons | 4 tons | 5 tons | 6 tons | 7 tons
---|---|---|---|---|---|---|---
20 mwe | 0.554 | 0.299 | 0.142 | 0.06 | 0.024 | 0.009 | 0.003
40 mwe | 0.551 | 0.296 | 0.138 | 0.059 | 0.023 | 0.008 | 0.003
60 mwe | 0.548 | 0.291 | 0.136 | 0.057 | 0.022 | 0.008 | 0.003
80 mwe | 0.546 | 0.288 | 0.133 | 0.054 | 0.021 | 0.008 | 0.003
100 mwe | 0.543 | 0.286 | 0.128 | 0.055 | 0.02 | 0.007 | 0.002
120 mwe | 0.541 | 0.281 | 0.128 | 0.053 | 0.02 | 0.007 | 0.002
140 mwe | 0.54 | 0.278 | 0.125 | 0.052 | 0.019 | 0.006 | 0.002

##  VI Conclusion

This exercise is entirely reactor power dependent to make determinations about
the steady state and reduced power signal. The diversifying advanced reactor
market requires new methods to safeguard nuclear material. As an example, the
current safeguards regime for LWRs is not designed for advanced fuel types
like molten salt reactors with liquid fuel and pebble bed reactors [48].
Neutrino detectors could be used as a way to determine the operational status
of these reactors, where traditional item accountancy may pose challenges.
This method can be applied to a range of advanced reactors. As discussed
previously, this paper determines the background and signal rate using
analytical methods and data from previous neutrino experiments to provide a
benchmark.

Depth likely has a minimal impact on the detection time and metrics as the
remaining 5 cores dominate the remaining signal rather than depth-dependent
muon-induced backgrounds.

The two diversion and misuse scenarios discussed in this paper are the
intentional shutdown of one core of six to divert fuel or the decrease in
power of all six cores to improve the ratio of Pu-239 in used fuel. For the
latter, it is assumed that all cores are reduced to 83 % of full operating
power making both of these scenarios appear identical with respect to neutrino
signals. Of course, reactors could operate at an even lower capacity to slow
other Pu isotopes production for longer. This would lead to even faster
detection times and better results with respect to the IAEA requirements.

Neutrino detectors can assist the IAEA with some of the anticipated challenges
behind SMR safeguards including increased deployment and remoteness of
reactors. The World Nuclear Association reports that around 30 countries are
beginning to develop nuclear programs [49]. This will inevitably create
increased costs and burden for the IAEA to safeguard these new reactors.
Neutrino detectors located outside of the reactor sites may offer the
capability for the IAEA to monitor these reactors remotely and offer an
additional layer of redundancy to the current C$\&$S safeguard regime.

Simulations using RAT(-PAC) can be used to compare these monte-carlo driven
predictions for signal and background to the analytical methods displayed in
this work [50]. Neutrino detectors are capable of measuring reactor power from
distances outside of a site boundary. This research investigated varying
detector sizes and depths to compare to IAEA metrics for success for C$\&$S.
Our results show that a 5 ton detector with an overburden of 20 mwe was able
to detect the reduction in power from 6 operating reactors to 5 in $\sim$ 10
days with a false and missed alarm rate of $\alpha=0.006$ and $\beta=0.024$,
respectively. Other methods of alarming systems may allow for smaller and
shallower detectors. This would reduce the cost of construction and operation.

## Acknowledgement

This material is based upon work supported by the Department of Energy
National Nuclear Security Administration under Award Number DE-NA0000979
through the Nuclear Science and Security Consortium. Lawrence Livermore
National Laboratory is operated by Lawrence Livermore National Security, LLC,
for the U.S. Department of Energy, National Nuclear Security Administration
under Contract DE-AC52-07NA27344. LLNL-JRNL-852145.

## References

  * [1] “Proceedings of the first meeting of the ad hoc working group on safeguards applications utilizing antineutrino detection and monitoring,” tech. rep., 2011\.
  * [2] T. Akindele, N. Bowden, R. Carr, A. Conant, M. Diwan, A. Erickson, M. Foxe, B. Goldblum, P. Huber, I. Jovanovic, J. Link, B. Littlejohn, H. Mumm, and J. Newby, “Nu tools: Exploring practical roles for neutrinos in nuclear energy and security,” tech. rep., aug 2021.
  * [3] O. A. Akindele, A. Bernstein, and E. B. Norman, “Antineutrino monitoring of thorium reactors,” Journal of Applied Physics, vol. 120, no. 12, p. 124902, 2016.
  * [4] A. Bernstein, N. Bowden, B. L. Goldblum, P. Huber, I. Jovanovic, and J. Mattingly, “Colloquium: Neutrino detectors as tools for nuclear security,” Reviews of Modern Physics, vol. 92, mar 2020.
  * [5] F. Boehm, J. Busenitz, B. Cook, G. Gratta, H. Henrikson, J. Kornis, D. Lawrence, K. B. Lee, K. McKinny, L. Miller, V. Novikov, A. Piepke, B. Ritchie, D. Tracy, P. Vogel, Y.-F. Wang, and J. Wolf, “Results from the palo verde neutrino oscillation experiment,” Phys. Rev. D, vol. 62, p. 072002, Sep 2000.
  * [6] R. Carr, Measurements of Electron Antineutrino Disappearance in the Double Chooz Experiment.  Columbia University, 2015.
  * [7] RENO Collaboration and J. K. Ahn, “Reno: An experiment for neutrino oscillation parameter $\theta_{13}$ using reactor neutrinos at yonggwang,” 2010\.
  * [8] J. Cao and K.-B. Luk, “An overview of the daya bay reactor neutrino experiment,” Nuclear Physics B, vol. 908, pp. 62–73, 2016.
  * [9] U. E. I. Administration, “What is u.s. electricity generation by energy source.” https://www.eia.gov/tools/faqs/faq.php?id=427&t=3 (11/08/2022).
  * [10] G. Locatelli, C. Bingham, and M. Mancini, “Small modular reactors: A comprehensive overview of their economics and strategic aspects,” Progress in Nuclear Energy, vol. 73, pp. 75–85, 2014.
  * [11] J. Whitlock and J. Sprinkle, “Proliferation resistance considerations for remote small modular reactors,” Nuclear Review, vol. 1, no. 2, pp. 9–12, 2014.
  * [12] S. Prasad, A. Abdulla, M. G. Morgan, and I. L. Azevedo, “Nonproliferation improvements and challenges presented by small modular reactors,” Progress in Nuclear Energy, vol. 80, pp. 102–109, 2015.
  * [13] J. K. Halbig, K. E. Gainer, and S. F. Klosterbuer, “An efficient network for interconnecting remote monitoring instruments and computers,” 8 1994.
  * [14] O. Akindele, A. Bernstein, M. Bergevin, S. Dazeley, F. Sutanto, A. Mullen, and J. Hecla, “Exclusion and verification of remote nuclear reactors with a 1-kiloton gd-doped water detector,” Physical Review Applied, vol. 19, no. 3, p. 034060, 2023.
  * [15] B. Pellaud, “Proliferation aspects of plutonium recycling,” Comptes Rendus Physique, vol. 3, no. 7-8, pp. 1067–1079, 2002.
  * [16] R. Mills, B. Slingsby, J. Coleman, R. Collins, G. Holt, C. Metelko, and Y. Schnellbach, “A simple method for estimating the major nuclide fractional fission rates within light water and advanced gas cooled reactors,” Nuclear Engineering and Technology, vol. 52, no. 9, pp. 2130–2137, 2020.
  * [17] S. Bogetic, R. Mills, A. Bernstein, J. Coleman, A. Morgan, and A. Petts, “Anti-neutrino flux from the edf hartlepool nuclear power plant,” arXiv preprint arXiv:2301.07153, 2023.
  * [18] N. R. C. (NRC), “Methods and assumptions for evaluating radiological consequences of design basis accidents at light-water nuclear power reactors,” tech. rep., 2003.
  * [19] W. L. Moe, “Opportunities in smr emergency planning,” 10 2014.
  * [20] A. Bernstein, N. Bowden, B. L. Goldblum, P. Huber, I. Jovanovic, and J. Mattingly, “Colloquium: Neutrino detectors as tools for nuclear security,” Reviews of Modern Physics, vol. 92, no. 1, p. 011003, 2020.
  * [21] F. An, J. Bai, A. Balantekin, H. Band, D. Beavis, W. Beriguete, M. Bishai, S. Blyth, K. Boddy, R. Brown, et al., “Observation of electron-antineutrino disappearance at daya bay,” Physical Review Letters, vol. 108, no. 17, p. 171803, 2012.
  * [22] Y. Abe, C. Aberle, J. Dos Anjos, J. Barriere, M. Bergevin, A. Bernstein, T. Bezerra, L. Bezrukhov, E. Blucher, N. Bowden, et al., “Reactor $\nu$ e disappearance in the double chooz experiment,” Physical Review D, vol. 86, no. 5, p. 052008, 2012.
  * [23] NuScale, “Nuscale standard plant design certification application.” https://www.nrc.gov/docs/ML2022/ML20224A482.pdf, 2020.
  * [24] P. D. Group, “Review of particle physics,” Phys. Rev. D, vol. 98, p. 030001, Aug 2018.
  * [25] T. A. Mueller, D. Lhuillier, M. Fallot, A. Letourneau, S. Cormon, M. Fechner, L. Giot, T. Lasserre, J. Martino, G. Mention, et al., “Improved predictions of reactor antineutrino spectra,” Physical Review C, vol. 83, no. 5, p. 054615, 2011.
  * [26] P. Huber, “Reactor antineutrino fluxes–status and challenges,” Nuclear Physics B, vol. 908, pp. 268–278, 2016.
  * [27] V. A. Li, S. A. Dazeley, M. Bergevin, and A. Bernstein, “Scalability of gadolinium-doped-water cherenkov detectors for nuclear nonproliferation,” Physical Review Applied, vol. 18, no. 3, p. 034059, 2022.
  * [28] M. Yeh, A. Garnov, and R. L. Hahn, “Gadolinium-loaded liquid scintillator for high-precision measurements of antineutrino oscillations and the mixing angle, $\theta$13,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 578, no. 1, pp. 329–339, 2007.
  * [29] J. Cao and K.-B. Luk, “An overview of the daya bay reactor neutrino experiment,” Nuclear Physics B, vol. 908, pp. 62–73, 2016.  Neutrino Oscillations: Celebrating the Nobel Prize in Physics 2015.
  * [30] C. N. Maesano, Lithium-9 Production from Stopped Muons in the Double Chooz Detector.  University of California, Davis, 2012.
  * [31] J. K. Ahn, S. Chebotaryov, J. Choi, S. Choi, W. Choi, Y. Choi, H. Jang, J. Jang, E. Jeon, I. Jeong, et al., “Observation of reactor electron antineutrinos disappearance in the reno experiment,” Physical Review Letters, vol. 108, no. 19, p. 191802, 2012.
  * [32] J. K. Ahn, S. Chebotaryov, J. H. Choi, S. Choi, W. Choi, Y. Choi, H. I. Jang, J. S. Jang, E. J. Jeon, I. S. Jeong, K. K. Joo, B. R. Kim, B. C. Kim, H. S. Kim, J. Y. Kim, S. B. Kim, S. H. Kim, S. Y. Kim, W. Kim, Y. D. Kim, J. Lee, J. K. Lee, I. T. Lim, K. J. Ma, M. Y. Pac, I. G. Park, J. S. Park, K. S. Park, J. W. Shin, K. Siyeon, B. S. Yang, I. S. Yeo, S. H. Yi, and I. Yu, “Observation of reactor electron antineutrinos disappearance in the reno experiment,” Phys. Rev. Lett., vol. 108, p. 191802, May 2012.
  * [33] S. Schoppmann, D.-C. Collaboration, et al., “Latest results of double chooz,” in XIII International Conference on Heavy Quarks and Leptons, vol. 274, p. 010, SISSA Medialab, 2017.
  * [34] M. Y. Pac, “Recent results from reno,” 2018.
  * [35] J. Park, J. Lee, I. Yeo, W. Choi, J. Ahn, J. Choi, S. Choi, Y. Choi, H. Jang, J. Jang, et al., “Production and optical properties of gd-loaded liquid scintillator for the reno neutrino detector,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 707, pp. 45–53, 2013\.
  * [36] R. Quéval, “Characterization, modelization and optimization of the double chooz acrylic vessels: Physics impact,” tech. rep., Paris VII Univ.(France), 2010\.
  * [37] F. Ardellier, I. Barabanov, J. Barriere, F. Beißel, S. Berridge, L. Bezrukov, A. Bernstein, T. Bolton, N. Bowden, C. Buck, et al., “Double chooz, a search for the neutrino mixing angle $\theta$13,” 2006.
  * [38] J. Crespo-Anadon, “Double chooz: latest results,” Nuclear and Particle Physics Proceedings, vol. 265, pp. 99–104, 2015.
  * [39] F. An, J. Bai, A. Balantekin, H. Band, D. Beavis, W. Beriguete, M. Bishai, S. Blyth, R. Brown, I. Butorov, et al., “The detector system of the daya bay reactor neutrino experiment,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 811, pp. 133–161, 2016.
  * [40] F. P. An, A. Balantekin, H. Band, M. Bishai, S. Blyth, D. Cao, G. Cao, J. Cao, W. Cen, Y. Chan, et al., “Measurement of electron antineutrino oscillation based on 1230 days of operation of the daya bay experiment,” Physical Review D, vol. 95, no. 7, p. 072006, 2017.
  * [41] F. An, A. Balantekin, H. Band, W. Beriguete, M. Bishai, S. Blyth, I. Butorov, G. Cao, J. Cao, Y. Chan, et al., “Search for a light sterile neutrino at daya bay,” Physical review letters, vol. 113, no. 14, p. 141802, 2014\.
  * [42] Y. Abe, J. Dos Anjos, J. Barriere, E. Baussan, I. Bekman, M. Bergevin, T. Bezerra, L. Bezrukov, E. Blucher, C. Buck, et al., “Improved measurements of the neutrino mixing angle $\theta$ 13 with the double chooz detector,” Journal of High Energy Physics, vol. 2014, no. 10, pp. 1–44, 2014.
  * [43] L. Bogdanova, M. Gavrilov, V. Kornoukhov, and A. Starostin, “Cosmic muon flux at shallow depths underground,” Physics of Atomic Nuclei, vol. 69, no. 8, pp. 1293–1298, 2006.
  * [44] C. D. Shin, Z. Atif, G. Bak, J. H. Choi, H. I. Jang, J. S. Jang, S. H. Jeon, K. K. Joo, K. Ju, D. E. Jung, J. G. Kim, J. Y. Kim, S. B. Kim, S. Y. Kim, W. Kim, E. Kwon, D. H. Lee, H. G. Lee, Y. C. Lee, I. T. Lim, D. H. Moon, M. Y. Pac, C. Rott, H. Seo, J. H. Seo, J. W. Seo, S. H. Seo, B. S. Yang, J. Y. Yang, J. Yoo, and I. Yu, “Observation of reactor antineutrino disappearance using delayed neutron capture on hydrogen at reno,” 2019.
  * [45] A. Piepke, P. V. Collaboration, et al., “Final results from the palo verde neutrino oscillation experiment,” Progress in Particle and Nuclear Physics, vol. 48, no. 1, pp. 113–121, 2002.
  * [46] IAEA Safeguards Glossary.  No. 3 (Rev. 1) in International Nuclear Verification Series, Vienna: INTERNATIONAL ATOMIC ENERGY AGENCY, 2022.
  * [47] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy,” Nature, vol. 585, pp. 357–362, Sept. 2020.
  * [48] R. Peel, G. Foster, and S. Aghara, “Nuclear security and safeguards considerations for novel advanced reactors,” King’s College London, London, UK, 2022.
  * [49] W. N. Association, “Emerging nuclear energy countries.” Available at https://world-nuclear.org/information-library/country-profiles/others/emerging-nuclear-energy-countries.aspx (2022/05).
  * [50] M. Askins, A. Bacon, Z. Bagdasarian, N. Barros, M. Bergevin, E. Callaghan, C. Grant, J. Gruszko, T. Kaptanoglu, J. Klein, et al., “Snowmass2021-letter of interest the rat (-pac) simulation and analysis code base,”


subsection section. subsection.

# Detection and Segmentation of Pancreas using Morphological Snakes and Deep
Convolutional Neural Networks

Agapi Davradou
<EMAIL_ADDRESS>

Instituto Superior Técnico, Lisboa, Portugal

(January 2021)

###### Abstract

Pancreatic cancer is one of the deadliest types of cancer, with 25% of the
diagnosed patients surviving for only one year and 6% of them for five.
Computed tomography (CT) screening trials have played a key role in improving
early detection of pancreatic cancer, which has shown significant improvement
in patient survival rates. However, advanced analysis of such images often
requires manual segmentation of the pancreas, which is a time-consuming task.
Moreover, pancreas presents high variability in shape, while occupying only a
very small area of the entire abdominal CT scans, which increases the
complexity of the problem. The rapid development of deep learning can
contribute to offering robust algorithms that provide inexpensive, accurate,
and user-independent segmentation results that can guide the domain experts.
This dissertation addresses this task by investigating a two-step approach for
pancreas segmentation, by assisting the task with a prior rough localization
or detection of pancreas. This rough localization of the pancreas is provided
by an estimated probability map and the detection task is achieved by using
the YOLOv4 deep learning algorithm. The segmentation task is tackled by a
modified U-Net model applied on cropped data, as well as by using a
morphological active contours algorithm. For comparison, the U-Net model was
also applied on the full CT images, which provide a coarse pancreas
segmentation to serve as reference. Experimental results of the detection
network on the National Institutes of Health (NIH) dataset and the pancreas
tumour task dataset within the Medical Segmentation Decathlon show 50.67% mean
Average Precision. The best segmentation network achieved good segmentation
results on the NIH dataset, reaching 67.67% Dice score.
Keywords: Pancreatic cancer, morphological snakes, deep learning,
segmentation, convolutional neural networks

##  1 Introduction

Cancer is the second leading cause of death worldwide, accounting for an
estimated 9.6 million deaths in 2018 [1]. Globally, about 1 in 6 deaths is
caused by cancer and pancreatic cancer is the seventh highest cause of death
from it [1]. In most cases, the symptoms are not visible until the disease has
reached an advanced stage. Due to its very poor prognosis, after diagnosis,
25% of people survive only one year and the 5-year survival rate is only 6%
[2].

This highlights the need for improved screening modalities and early
detection. In this context, radiomics, the process of applying machine
learning algorithms to extract meaningful information and automate image
analysis processes from computed tomography (CT) scans or magnetic resonance
imaging (MRI), has become a very active research area. This allows the
noninvasive characterization of lesions and the assessment of their
progression and possible response to therapies. However, the application of
such techniques, besides requiring high quality and reproducible data, often
require manual segmentation of the structures of interest which is a time-
consuming task. Moreover, manual segmentation is also user-dependent.
Consequently, it benefits from automatic segmentation techniques.

Automatic pancreas segmentation still remains a very challenging task, as the
target often occupies a very small fraction (e.g., $<$ 0.5% of an abdominal
CT) of the entire volume, has poorly-defined boundaries with respect to other
tissues and suffers from high variability in shape and size.

The goal of this dissertation is to develop a model to automate pancreas
segmentation in CT scans, in order to help the medical community to produce
faster and more consistent results.

The next section describes related work, while Section 3 presents the proposed
methods and implementation details. Section 4 details the obtained results.
Finally, Section 5 summarizes the conclusions and presents directions for
future work.

##  2 Related Work

###  2.1 Deep Learning for Pancreas Segmentation

Early work on pancreas segmentation from abdominal CT used statistical shape
models or multi-atlas techniques. In these approaches, the Dice similarity
coefficient (DSC) or Dice score in the public National Institutes of Health
(NIH) would not exceed 75%. Therefore, convolutional neural networks (CNNs)
have rapidly become the mainstream methodology for medical image segmentation.
Despite their good representational power, it was observed that such deep
segmentation networks are easily disrupted by the varying contents in the
background regions, when detecting small organs such as the pancreas, and as a
result produce less satisfying results. Taking that into consideration, a
coarse-to-fine approach is commonly adopted. These cascaded frameworks extract
regions of interest (RoIs) and make dense predictions on that particular RoIs.
More specifically, the state-of-the art methods primarily fall into two
categories.

The first category is based on segmentation networks originally designed for
2D images, such as the fully convolutional networks and the U-Net [33]. The
U-Net architecture utilizes up-convolutions to make use of low-level
convolutional feature maps by projecting them back to the original image size,
which delineates object boundaries with details. Many frameworks have used a
variant of the U-Net architecture to segment the pancreas [25], [27] and
others have made use of fully convolutional network (FCN) and U-Net in order
to build more complex models. For example, Zhou et al.[49] finds the rough
pancreas region and then a trained FCN-based fixed-point model refines the
pancreas region iteratively. Roth et al. [36] first segment pancreas regions
by holistically-nested networks and then refines them by the boundary maps
obtained by robust spatial aggregation using random forests. In addition, the
TernaryNet proposed by Heinrich et al. [17] applies sparse convolutions to the
CT pancreas segmentation problem, which reduce the computational complexity by
requiring a lower number of non-zero parameters. Combinations of convolutional
neural networks with recurrent neural networks [46], [44], [5] have also been
applied, as well as long short-term memory (LSTM) networks [24]. Recently,
Zheng et al. [48] proposed a model, which can involve uncertainties in the
process of segmentation iteratively, by utilizing the shadowed sets theory. In
some cases, in order to incorporate spatial 3D contextual information through
the 2D segmentation, slices along different views (axial, sagittal, and
coronal) are used, by fusing the results of all 2D networks, e.g. through
majority voting.

In the second category, the methods are based on 3D convolutional layers and
therefore operations such as 2D convolution, 2D max-pooling, and 2D up-
convolution are replaced by their 3D counterparts. Such networks are the 3D
U-Net [9] (which is a 3D extension of the U-Net), the Dense V-Net [15], ResDSN
[52], 3DFCN [37], and more [47]. In addition, OBELISK-Net proposed by Heinrich
et al. [18] applies sparse convolutions to 3D U-Net, Oktay et al. [31] applies
attention gates (AG) and Khosravan et al. [20] utilizes Projective Adversarial
Network to perform pancreas segmentation. Finally, Zhu et al. [51] applies
neural architecture search, to automatically find optimal network
architectures between the 2D, 3D and pseudo3D convolutions.

###  2.2 Active Contours and Morphological Snakes

Active contours, also known as ”snakes”, were first introduced in 1988 by Kass
et al. [19]. These are energy-minimizing methods, which are used extensively
in medical image processing, as a segmentation technique. The idea is to
initialize a position for the contour, and then define image forces that act
on the contour, making it change its position and adapt to the image’s
features. For instance, Kass et al. [19] define three energies: the internal
energy ($E_{int}$), the image energy ($E_{image}$) and the constraint energy
($E_{con}$), which represent the internal energy of the snake due to bending,
the image’s energy (which takes into consideration, for example, edges and
lines) and the energy of external constraint forces (which takes into
consideration constraints created by the user), respectively. In this work,
the authors propose different expressions for each energy, which are beyond
the scope of this work and will not be presented. After defining all the
energies, Kass et al. [19] present their active contour algorithm as the
solution of the following energy minimization problem:

$\displaystyle E_{snake}^{*}$
$\displaystyle:=\int_{0}^{1}E_{snake}\big{(}v(s)\big{)}\hskip 2.84544ptds$ (1)
$\displaystyle:=\int_{0}^{1}E_{int}\big{(}v(s)\big{)}$
$\displaystyle+E_{image}\big{(}v(s)\big{)}+E_{con}\big{(}v(s)\big{)}\hskip
2.84544ptds,$

where $\displaystyle\int_{C}f(r)\hskip
2.84544ptds:=\int_{a}^{b}f\big{(}r(t)\big{)}\hskip
2.84544pt\big{|}r^{\prime}(t)\big{|}\hskip 2.84544ptdt$ is the line integral
of $f$ along a piecewise smooth curve $C$, $r:[a,b]\rightarrow C$ an arbitrary
bijective parametrization of the curve $C$ such that $r(a)$ and $r(b)$ are the
endpoints of $C$, with $a<b$.

Since the minimization of energy leads to dynamic behavior in the segmentation
and because of the way the contours slither while minimizing energy, Kass et
al. [19] called them snakes. It should also be noted that this minimization
problem is generally not convex, which means that different snake
initializations may lead to different segmentations. In order to overcome the
problems of bad initialization and local minima, many variants of this method
have been proposed, such as using a balloon force to encourage the contour
expansion [11] or incorporating gradient flows [21]. Geometric models for
active contours were also introduced, by Caselles et al. [6], Yezzi et al.
[45], and more [41], [26], as well as a geodesic active contour (GAC) model by
Caselles et al. [7]. [12] et al. incorporated statistical shape knowledge in a
single energy functional, by modifying the Mumford-Shah functional [30] and
its cartoon limit and Chan and Vese [8] introduced a region-based method which
is not using an edge-function, but it minimizes an energy, which can be seen
as a particular case of the minimal partition problem. Additional region-based
methods were proposed, such as those by Li et al. [23], Zhu and Yuille [50],
and Tsai et al. [43].

Inspired by the active contour evolution, a new framework for image
segmentation was proposed in [3], which the authors called Morphological
Snakes. Instead of computing the snake that minimizes $E_{snake}^{*}$, Alvarez
et al. [3] proposed a new approach, which focuses on finding the solution of
certain set of partial differential equations (PDEs). This approach also
yields a snake-like curve and because it approximates the numerical solution
of the standard PDE snake model by the successive application of a set of
morphological operators (such as dilation or erosion) defined on a binary
level-set, it results in a much simpler, faster and stable, curve evolution.

The same authors [28] have introduced morphological versions of two of the
most popular curve evolution algorithms: morphological active contours without
edges (MorphACWE) [8] and morphological geodesic active contours (MorphGAC)
[7].

MorphACWE works well when pixel values of the inside and the outside regions
of the object to segment have different average. It does not require that the
contours of the object are well defined, and it can work over the original
image without any preprocessing.

MorphGAC is suitable for images with visible contours, even when these
contours are noisy, cluttered, or partially unclear. It requires, however,
that the image is preprocessed to highlight the contours and the quality of
the MorphGAC segmentation depends greatly on this preprocessing step.

Considering the application of this work, the MorphGAC algorithm was adopted

##  3 Methodology

###  3.1 Datasets

This work relied on two datasets, the pancreas tumour task dataset within the
Medical Segmentation Decathlon [40] and the NIH (National Institutes of
Health) dataset [34], [35], [10].

The Decathlon dataset is a recent multi-institutional effort to generate a
large, open-source collection of annotated medical image datasets of various
clinically relevant anatomies. The pancreas dataset was collected by the
Memorial Sloan Kettering Cancer Center (New York, NY, USA) and contains 420 3D
CT scans, 282 for training and 139 for testing. It consists of patients
undergoing resection of pancreatic masses (intraductal mucinous neoplasms,
pancreatic neuroendocrine tumours, or pancreatic ductal adenocarcinoma). The
CT scans have slice thickness of 2.5 mm and a resolution of 512$\times$512
pixels. An expert abdominal radiologist performed manual slice-by-slice
segmentation of the pancreatic parenchyma and pancreatic mass (cyst or
tumour), using the Scout application [13].

The NIH dataset contains 80 abdominal contrast-enhanced 3D CT scans from 53
male and 27 female subjects. It consists of healthy patients, since 17 of them
are healthy kidney donors scanned prior to nephrectomy and the remaining 65
subjects were selected by a radiologist from patients who neither had major
abdominal pathologies nor pancreatic cancer lesions. The CT scans have slice
thickness between 1.5 and 2.5 mm and a resolution of 512$\times$512 pixels,
with varying pixel sizes. The pancreas was manually segmented in each slice by
a medical student and the segmentations were verified/modified by an
experienced radiologist.

Due to computational power limit, only a certain amount of data was used. More
specifically, the detection model was trained and evaluated using both
datasets and the segmentation model was trained and evaluated using only the
NIH dataset. More details are given in Sections 3.2.2 and 3.3.2, respectively.
In addition, only the information from the transversal plane was taken into
consideration. In order to evaluate the variability in the data, a probability
map was recreated to define the most likely position of the pancreas.
Moreover, the average Hounsfield range and percentage of the volume occupied
by the pancreas was also evaluated. The results are presented in Section 4.1.

###  3.2 Detection Network Architecture

####  3.2.1 Pre-processing

Pre-processing is a crucial step in machine learning in order to improve the
performance of the models. For the detection task, the intensity of all data
was clipped between -200 to +300 HU, in order to capture the pancreas
intensity range and intensify its boundaries. More specifically, this range
was carefully chosen considering the data exploration results that will be
later presented in Section 4.1. The -125 to +225 intensity range was also
tested for clipping, as well as using 16-bit input images, in order to provide
the algorithm with more information. However, both proved to be inefficient
and were not implemented.

####  3.2.2 YOLOv4

In this work, the publicly available YOLOv4 [4] was deployed as a detection
network. The detection network was trained with both the Decathlon and the NIH
datasets. These datasets were shuffled and split into training and validation
sets. More specifically, 66 scans from NIH and 34 scans from Decathlon (100 in
total) were used for training with 20% validation. The hold-out method was
used to evaluate the performance of the model, by using 10% of the scans (5 of
each dataset) for testing.

Relying on studies regarding the effectiveness of transfer learning in deep
networks [32], [38], weights pre-trained on the ImageNet dataset [22] were
used. Afterwards, YOLOv4 was fine-tuned and re-trained with the pancreas
images. The model was trained for 20000 iterations with a batch size of 64 and
a learning rate of 0.0013.

####  3.2.3 Evaluation Metrics

The performance of all models was evaluated through the hold-out methodology,
using the corresponding test datasets. For the quality assessment of the
predictions, various metrics were used, i.e., precision, recall, IoU, and mean
average precision (mAP) for the detection task. All of the metrics are briefly
introduced in the following paragraphs.

Precision, also known as positive predictive value (PPV), measures how
accurate the model’s predictions are, i.e., the percentage of the predictions
that are correct. It is given as the ratio of true positives and the total
number of predicted positives, as seen in the following expression:

$PPV=\frac{TP}{TP+FP},$ (2)

where TP denotes the true positives (predicted correctly as positive) and FP
the false positives (predicted incorrectly as positive).

Similarly, recall, also known as sensitivity or true positive rate (TPR),
measures how well all the positives are predicted. It is given as the ratio of
the true positives and the total of ground truth positives, as seen in the
following equation:

$TPR=\frac{TP}{TP+FN},$ (3)

where FN are the false negatives (incorrectly predicted as negatives).

The IoU metric, also known as Jaccard index, is a popular similarity measure
for object detection problems using the predicted and ground-truth bounding
boxes. Evidently, as exemplified in Figure 1, a bigger overlap between the two
bounding boxes results in higher IoU score and therefore, detection accuracy.

Figure 1: Calculation of the IoU metric. With purple color the predicted
bounding box is depicted and with green the ground truth.

The average precision (AP) is a default evaluation metric in the PascalVOC
competition [14, p. 313], which derives from the area under the
precision/recall curve. The precision/recall curve is computed from a model’s
ranked output and the predictions are ranked with respect to the confidence
score of each bounding box. In this work, the detection model is set to keep
only predictions with a confidence score higher than 25%. The AP provides an
indication of the shape of the precision/recall curve, and is defined as the
interpolated average precision at a set of eleven equally spaced recall levels
[0, 0.1,…, 1], expressed as:

$AP=\frac{1}{11}\sum_{r\in\\{0,0.1,...,1\\}}p_{interp}(r).$ (4)

At each recall level $r$, the interpolated precision $p_{interp}$ is
calculated by taking the maximum precision measured for that $r$, given by the
following formula:

$p_{interp}(r)=\underset{\tilde{r}:\tilde{r}\geq r}{max}\ p(\tilde{r}).$ (5)

Since mAP is calculated by taking the average of the AP calculated for all the
classes, they will be used interchangeably in the current context.

###  3.3 Segmentation Network Architecture

####  3.3.1 Pre-processing

Intensity clipping was not adopted for the segmentation task, since it did not
show any improvement in the performance of the models. For the segmentation
network, curvature driven image denoising is applied on each slice, in order
to make the pixel distribution more uniform. Finally, contrast limited
adaptive histogram equalization (CLAHE) and normalization were also
investigated for the segmentation task, but they also proved to be
ineffective, in terms of Dice score.

####  3.3.2 U-Net model

An adaptation of the original U-Net architecture was chosen as a segmentation
model. More specifically, following the original implementation, 3$\times$3
convolutional kernels are used in the contracting path with a stride of 1,
each followed by a 2$\times$2 max-pooling operation with stride 2, and
2$\times$2 up-convolutional kernels with a stride of 2 are implemented in the
expansive path. Each convolution layer, except for the last one, is followed
by a ReLU, and dropout of 0.5 is selected. However, the number of levels in
both the descending and ascending paths of the network was increased from four
blocks to five blocks and the number of filters per layer was reduced compared
to the original implementation, as visualized in Figure 2. In addition,
residual connections were included in each convolutional block, padding was
added and the sigmoid function was used as an activation function for the last
layer. Due to memory limitations, all scans were resized from 512$\times$512
to 256$\times$256.

Figure 2: The 5-level U-Net segmentation network for an input image of size
256x256.

As mentioned in Section 3.1, the U-Net model was trained on 50 volumes from
the NIH dataset, 45 of them were used for training and 5 for validation.

Since the pancreas occupies only a very small region of a CT-Scan, the work
[29] is followed, and a DSC-loss layer is used to prevent the model from being
heavily biased towards the background class. In more detail, the DSC between
two voxel sets, $A$ and $B$, can be expressed as

$DSC(A,B)=\frac{2x|A\cap B|}{|A|+|B|},$ (6)

and this is slightly modified into a loss function between the ground-truth
mask Y and the predicted mask $\hat{Y}$, in the following way:

$L(\hat{Y},Y)=1-\frac{2x\sum_{i}\hat{y}_{i}y_{i}+\epsilon}{\sum_{i}\hat{y}_{i}+\sum_{i}y_{i}+\epsilon},$
(7)

where $\epsilon$ is an added term to avoid underflow.

The model is trained with batches of 128 instances and optimized using Adam
with a learning rate of 0.0001. Early-stopping is also used, in order to avoid
overfitting of the network, and the model with the lowest validation loss is
chosen.

Data augmentation is also implemented, by applying elastic deformation of
images as described in [39], as well as image shifting, rotation, zooming and
flipping, in order to expand the training dataset and make the network more
robust to such variations.

####  3.3.3 Cropped U-Net model

In order to reduce the irrelevant information in a CT-scan and train the U-Net
on scans with more useful information, the scans were cropped to a smaller
size. The rest of the training procedure described in Section 3.3.2 was kept
intact.

The image size of the cropped scans was carefully decided, as on the one hand
it should be as small as possible, but also contain the whole pancreas. Taking
into account the information retrieved from Figure 7, which shows the
probability of pancreas existing in a 512$\times$512 CT slice, as well as the
maximum size of pancreas in both datasets, a default cropping position was
decided. The default cropping position has its centroid at position
$(x=287,y=250)$ and image size of 224$\times$224, as visualized in Figure 3.

Figure 3: Default cropping position of a 512$\times$512 CT-scan into a
224$\times$224 scan.

####  3.3.4 YOLO + MorphGAC model

In order to overcome the problem of initialization, the MorphGAC algorithm was
combined with the YOLO detection model. The MorphGAC segmentation was applied
on cropped bounding box predictions using their centroid as an initialization
point. The segmentation results are then repositioned back to 512$\times$512
images, in order to form the final 3D segmented pancreas.

####  3.3.5 YOLO + U-Net model

An additional two-step approach with cropped CT-scans is proposed, which
combines the YOLO detection model to crop the images containing the pancreas,
as well as the default cropping position mentioned in Section 3.3.3, to crop
the images without pancreas. The U-Net architecture visualized in Figure 2 and
the training procedure described in 3.3.2 were used.

####  3.3.6 Post-processing

A post-processing step is adopted for all predicted segmentations, aiming for
optimization and a smoother result, by leveraging spatial information. Since
pancreas has a uniform and undivided shape, it is unlikely that a pixel in a
2D slice has a different value from both the previous and the next slice.
Taking that into consideration, all slices in a predicted segmentation are
compared in groups of three and the values of the middle ones are changed,
when different from the other two. In Figure 4, the segmentation improvement
on the middle slice is depicted, after post-processing it.

Figure 4: Example of the application of the optimization step on a slice. The
optimization is applied only on the middle slice of the three used.

####  3.3.7 Evaluation Metrics

The performance of all models was evaluated through the hold-out methodology,
using the corresponding test datasets. Since MorphGAC algorithm needs no
training, the YOLO+MorphGAC approach was evaluated on the detection’s model
test dataset. For the quality assessment of the segmentation task, the DSC or
Dice score was used. Because of the similarity between the IoU metric and the
Dice score, the former was not used for the evaluation of the segmentation
models.

The Dice score is the most common evaluation metric for the segmentation of
medical images and it is calculated as twice the area of overlap divided by
the total number of pixels in both images, as illustrated in Figure 5.

Figure 5: Calculation of the DSC metric.

Evidently, the DSC is equivalent to the IoU, using the following expression:

$DSC=\frac{2\mbox{IoU}}{1+\mbox{IoU}}.$ (8)

##  4 Results & discussion

###  4.1 Data Exploration

The average Hounsfield range and percentage of the volume occupied by the
pancreas were evaluated. In Figure 6 the pancreas intensity values for both
datasets are illustrated, where Decathlon has a mean +80.63 $\pm$ 57.91 HU and
NIH +86.71 $\pm$ 32.18 HU.

Figure 6: The Hounsfield unit values of the pancreas for the Decathlon (a) and
NIH (b) dataset, respectively. In Decathlon dataset the intensity range from
-998.0 HU to +3071.0, while in NIH from -875.0 HU to +458.0 HU.

In addition, the mean percentage of pancreas in the abdominal ct scans is
0.46% and 0.49% for Decathlon and NIH respectively, taking into consideration
only values $<$800 HU, in order to exclude the air. Moreover, a probability
map was recreated to define the most likely position of the pancreas. In
Figure 7 is the probability map of pancreas in a 2D slice for both datasets is
visualized. In the Decathlon dataset, the pancreas location ranges from pixel
150 to 434 in the x-axis and from pixel 139 to 348 in the y-axis. Similarly,
in NIH dataset, the pancreas location ranges from 167 to 405 in the x-axis and
from 143 to 360 in the y-axis.

Figure 7: Probability of pancreas in an 512x512 image for Decathlon (a) and
NIH (b) datasets, respectively.

###  4.2 Detection Task

The detection model was evaluated on the holdout test set containing volumes
from both the NIH and the Decathlon datasets and it can predict pancreas
successfully in both of them, as visualized in Figure 8.

Figure 8: Detection of pancreas on a NIH (a) and Decathlon (b) slice,
respectively. The green color represents the ground-truth and the purple the
YOLO prediction.

The results of the detection model evaluated by the PPV, TPR, IoU, and mAP
metrics described in Section 3.2.3 are presented in Table 1. The detection
model shows a mAP of 50.67% on the test set, with the IoU being 47.70%, the
precision 0.63, and the recall 0.52. However, the performance of the model
varies significantly between the two datasets, with the mAP being 71.43% for
the NIH dataset and 29.92% for the Decathlon dataset. For that reason, the
performance of the model in Table 1 is also evaluated separately for each
dataset. More specifically, the model shows a IoU of 57.43% with a precision
and recall of 0.63 and 0.52 respectively, on the NIH dataset. On the Decathlon
dataset, the model has 37.96% IoU, the precision is 0.5 and the recall 0.36.

Table 1: Evaluation results of the detection model.

Test dataset |  PPV |  TPR |  IoU |  mAP
---|---|---|---|---
NIH & Decathlon |  0.63 |  0.52 |  47.70% |  50.67%
NIH |  0.75 |  0.68 |  57.43% |  71.43%
Decathlon |  0.5 |  0.36 |  37.96% |  29.92%

The low detection accuracy of the YOLO model on the Decathlon dataset derives
from the multi-label nature of some ground-truth files. More specifically, it
was later realized the YOLO model was configured to be trained on one class
and therefore, since the pancreatic parenchyma and pancreatic mass (cyst or
tumour) are annotated separately as different classes in the Decathlon
dataset, the second class was ignored. as shown in Figure 9. In addition,
during the conversion from NIfTI to tiff format (which is required by YOLO),
all the individual globs of healthy pancreas surrounding the pancreatic mass,
are annotated as separate tiny bounding boxes. However, as visualized on this
figure, it is very interesting to notice that the model outputs a fairly good
prediction on the class that it was trained on (upper-right green bounding
box).

Figure 9: Detection of pancreas on a multilabel slice. The green color
represents the ground-truth and the purple the YOLO prediction. The two bigger
green bounding boxes represent the two classes.

Unfortunately, the problem was not discovered on time, since the YOLO model
was trained only on the NIH dataset at first with healthy pancreas.
Nevertheless, this could be overcome by considering a single class once the
goal of the work was to segment the healthy pancreas. Due to the
unavailability of the computer on which YOLO was trained afterwards, this
could not be realized in the current dissertation.

In Figures 10 and 11, the metrics used for the assessment of the detection
model during training are visualized.

Figure 10: The CIoU training loss with a minimum value of 0.18. Figure 11: The
validation mAP during training, reaching a maximum value of 90.9%.

###  4.3 Segmentation Task

The segmentation models were evaluated on the holdout set containing volumes
from the NIH dataset and the pancreas can be successfully segmented, as
visualized in Figure 13.

The segmentation performance of all models on the NIH dataset is presented on
Table 2. The YOLO+MorphGAC model has the best performance achieving a dice
score of 67.67 $\pm$ 8.62%, followed by the YOLO+U-Net model with a Dice score
of 64.87 $\pm$ 4.79%. The U-Net and cropped U-Net show a lower performance
with a Dice score of 59.91 $\pm$ 13.69% and 59.08 $\pm$ 12.76%, respectively.
The obtained DSC from the cropped U-Net when compared with the YOLO+U-Net was
lower and showed a higher standard deviation. The preliminary results suggest
that YOLO was able to improve segmentation when compared to using a
probability map. A similar behavior was observed for the U-Net model with
improved results when using the YOLO model. Regarding segmentation
performance, the highest average Dice score was obtained with the
YOLO+MorphGAC model and the lowest standard deviation with the YOLO+U-Net
model.

Table 2: Evaluation results of all proposed segmentation models. Test dataset |  Dice score
---|---
YOLO+MorphGAC |  67.67 $\pm$ 8.62 %
U-Net |  59.91 $\pm$ 13.69 %
Cropped U-Net |  59.08 $\pm$ 12.76 %
YOLO+U-Net |  64.87 $\pm$ 4.79 %

The loss function used for the assessment of the segmentation model during
training is visualized in Figure 12.

Figure 12: The metric of loss for the YOLO+U-Net during training with a
minimum value of 0.28.

In Figures 13 and 14 representative slices of the best and worst results of
the proposed segmentation models are presented, respectively.

Figure 13: Some of the best segmentations on the test set of the YOLO+MorphGAC
(a), U-Net (b), cropped U-Net (c) and YOLO+U-Net (d) model, respectively. The
blue line indicates the model’s predicted segmentation and the red line
represents the ground-truth. Figure 14: Some of the worst segmentations on the
test set of the YOLO+MorphGAC (a), U-Net (b), cropped U-Net (c) and YOLO+U-Net
(d) model, respectively. The blue line indicates the model’s predicted
segmentation and the red line represents the ground-truth.

####  4.3.1 Comparison with the state-of-the-art

In Tables 3 and 4 a comparison of the best two proposed models with the state-
of-the art deep learning models is presented with respect to the dice
coefficient and the training details, respectively. Only models which
implement a 2D approach were taken into consideration and trained on the NIH
dataset as well. The proposed models display a low performance compared to the
other networks. However, when examining the training methodology of all the
state-of-the-art methods, it is concluded that cross validation is the adopted
validation technique, as well as that they leverage the whole NIH dataset for
the training, which in the current work was not possible, due to computational
and time limitations. These two reasons could be the main causes for the low
segmentation performance of the proposed YOLO+U-Net model and following the
state-of-the-art training methods could result in a higher a Dice score.

Table 3: Comparison of the DSC with the state-of-the-art methods.

Method |  DSC
---|---
2D U-Net [17] |  71.07% $\pm$ 9.50
2D U-Net [27] |  86.93% $\pm$ 4.92
2D FCN + RNN [5] |  83.70% $\pm$ 5.10
Holistically Nested 2D FCN [36] |  81.27% $\pm$ 6.27
2D FCN [49] |  83.18% $\pm$ 4.81
YOLO+U-Net (proposed) |  64.87% $\pm$ 4.79
YOLO+MorphGAC (proposed) |  67.67% $\pm$ 8.62

Table 4: Comparison of the training details with the state-of-the-art methods.

Method |  Dataset |  Train/Test |  Validation
---|---|---|---
2D U-Net [17] |  NIH |  63/9 |  5-fold CV
2D U-Net [27] |  NIH |  62/20 |  4-fold CV
2D FCN + RNN [5] |  NIH |  62/20 |  4-fold CV
Holistically Nested 2D FCN [36] |  NIH |  62/20 |  4-fold CV
2D FCN [49] |  NIH |  62/20 |  4-fold CV
YOLO+U-Net (proposed) |  NIH |  50/5 |  holdout
YOLO+MorphGAC (proposed) |  - |  -/5 |  -

##  5 Conclusion and Future Work

###  5.1 Conclusions

Taking into account the state-of-the-art cascaded segmentation models, the
main goal of the current work was to investigate and present simpler two-step
approaches for the segmentation of the pancreas in CT. In the end, three
models were investigated and compared to the U-Net implementation. The two-
step approach was achieved by using a detection network for the pancreas
localization prior to the segmentation. This approach, which was combined with
both morphological snakes and a U-Net segmentation network, proved to be the
most efficient, showing that the organ detection benefits the segmentation
task. Another approach was to reduce the background information by cropping
the data using a pancreas probability map and then using a U-Net network to
segment this smaller area. However, this showed no improvement in segmentation
performance, when compared to U-Net segmentation results on whole CT-scans.
Due to the Covid-19 situation, many computational limitations existed, which
resulted in insufficient training of both the detection and the segmentation
network, which is the main reason of the low performance of the proposed
models. Since the segmentation performance lies below 70%, it is concluded
that further investigation and model improvement is needed, in order to tackle
the challenging nature of the pancreas segmentation problem efficiently.

###  5.2 Future Work

Regarding the proposed models, there are changes that could be investigated to
further improve the performance, focusing on YOLO+MorphGAC and YOLO+U-Net
architectures. First and foremost, since the detection network plays a vital
role on both of them, its improvement would also lead to overall improved
performance of the model. Considering that YOLOv4 is one of the most modern
and efficient state-of-the-art detection model, experimenting with different
ones would not be the main priority. A slight improvement could occur by
testing different backbone models for feature extraction, such as ResNet [16]
or EfficientNet [42]. More importantly, leveraging the whole Decathlon dataset
and training on a larger dataset, could also lead to a better detection
performance. Regarding the proposed U-Net architecture, similarly to the
detection task, the need to train on a larger dataset is undeniable and adopt
the training methodology of the state-of-the-art models. In addition, a
YOLO+FCN approach would also be interesting for investigation. Finally,
exploring 3D versions of the architectures could also improve the results. The
present work focuses on the segmentation of healthy pancreas, nevertheless,
once Decathlon dataset also presents tumor masks, tumor segmentation could
also be explored.

## References

  * [1] Who cancer fact sheet 2018.  https://www.who.int/news-room/fact-sheets/detail/cancer.  Accessed: 2020-09-26.
  * [2] A. Acs.  Cancer facts and figures 2010.  American Cancer Society, National Home Office, Atlanta, pages 1–44, 2010.
  * [3] L. Álvarez, L. Baumela, P. Henríquez, and P. Márquez-Neila.  Morphological snakes.  In 2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, pages 2197–2202. IEEE, 2010.
  * [4] A. Bochkovskiy, C.-Y. Wang, and H.-Y. M. Liao.  Yolov4: Optimal speed and accuracy of object detection.  arXiv preprint arXiv:2004.10934, 2020.
  * [5] J. Cai, L. Lu, F. Xing, and L. Yang.  Pancreas segmentation in ct and mri images via domain specific network designing and recurrent neural contextual learning.  arXiv preprint arXiv:1803.11303, 2018.
  * [6] V. Caselles, F. Catté, T. Coll, and F. Dibos.  A geometric model for active contours in image processing.  Numerische mathematik, 66(1):1–31, 1993.
  * [7] V. Caselles, R. Kimmel, and G. Sapiro.  Geodesic active contours.  In Proceedings of IEEE international conference on computer vision, pages 694–699. IEEE, 1995.
  * [8] T. F. Chan and L. A. Vese.  Active contours without edges.  IEEE Transactions on image processing, 10(2):266–277, 2001.
  * [9] Ö. Çiçek, A. Abdulkadir, S. S. Lienkamp, T. Brox, and O. Ronneberger.  3d u-net: learning dense volumetric segmentation from sparse annotation.  In International conference on medical image computing and computer-assisted intervention, pages 424–432. Springer, 2016.
  * [10] K. Clark, B. Vendt, K. Smith, J. Freymann, J. Kirby, P. Koppel, S. Moore, S. Phillips, D. Maffitt, M. Pringle, et al.  The cancer imaging archive (tcia): maintaining and operating a public information repository.  Journal of digital imaging, 26(6):1045–1057, 2013.
  * [11] L. D. Cohen.  On active contour models and balloons.  CVGIP: Image understanding, 53(2):211–218, 1991.
  * [12] D. Cremers, C. Schnorr, and J. Weickert.  Diffusion-snakes: combining statistical shape knowledge and image information in a variational framework.  In Proceedings IEEE Workshop on Variational and Level Set Methods in Computer Vision, pages 137–144. IEEE, 2001.
  * [13] B. M. Dawant, R. Li, B. Lennon, and S. Li.  Semi-automatic segmentation of the liver and its evaluation on the miccai 2007 grand challenge data set.  3D Segmentation in The Clinic: A Grand Challenge, pages 215–221, 2007.
  * [14] M. Everingham, L. Van Gool, C. K. Williams, J. Winn, and A. Zisserman.  The pascal visual object classes (voc) challenge.  International journal of computer vision, 88(2):303–338, 2010.
  * [15] E. Gibson, F. Giganti, Y. Hu, E. Bonmati, S. Bandula, K. Gurusamy, B. Davidson, S. P. Pereira, M. J. Clarkson, and D. C. Barratt.  Automatic multi-organ segmentation on abdominal ct with dense v-networks.  IEEE transactions on medical imaging, 37(8):1822–1834, 2018.
  * [16] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [17] M. P. Heinrich, M. Blendowski, and O. Oktay.  Ternarynet: faster deep model inference without gpus for medical 3d segmentation using sparse and binary convolutions.  International journal of computer assisted radiology and surgery, 13(9):1311–1320, 2018.
  * [18] M. P. Heinrich, O. Oktay, and N. Bouteldja.  Obelisk-net: Fewer layers to solve 3d multi-organ segmentation with sparse deformable convolutions.  Medical image analysis, 54:1–9, 2019.
  * [19] M. Kass, A. Witkin, and D. Terzopoulos.  Snakes: Active contour models.  International journal of computer vision, 1(4):321–331, 1988.
  * [20] N. Khosravan, A. Mortazi, M. Wallace, and U. Bagci.  Pan: Projective adversarial network for medical image segmentation.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 68–76. Springer, 2019.
  * [21] S. Kichenassamy, A. Kumar, P. Olver, A. Tannenbaum, and A. Yezzi.  Gradient flows and geometric active contour models.  In Proceedings of IEEE International Conference on Computer Vision, pages 810–815. IEEE, 1995.
  * [22] A. Krizhevsky, I. Sutskever, and G. E. Hinton.  Imagenet classification with deep convolutional neural networks.  In Advances in neural information processing systems, pages 1097–1105, 2012.
  * [23] C. Li, C.-Y. Kao, J. C. Gore, and Z. Ding.  Minimization of region-scalable fitting energy for image segmentation.  IEEE transactions on image processing, 17(10):1940–1949, 2008.
  * [24] H. Li, J. Li, X. Lin, and X. Qian.  Pancreas segmentation via spatial context based u-net and bidirectional lstm.  arXiv preprint arXiv: 1903.00832, 2019.
  * [25] L. Lu, L. Jian, J. Luo, and B. Xiao.  Pancreatic segmentation via ringed residual u-net.  IEEE Access, 7:172871–172878, 2019.
  * [26] R. Malladi, J. A. Sethian, and B. C. Vemuri.  Shape modeling with front propagation: A level set approach.  IEEE transactions on pattern analysis and machine intelligence, 17(2):158–175, 1995.
  * [27] Y. Man, Y. Huang, J. Feng, X. Li, and F. Wu.  Deep q learning driven ct pancreas segmentation with geometry-aware u-net.  IEEE transactions on medical imaging, 38(8):1971–1980, 2019.
  * [28] P. Marquez-Neila, L. Baumela, and L. Alvarez.  A morphological approach to curvature-based evolution of curves and surfaces.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 36(1):2–17, 2013.
  * [29] F. Milletari, N. Navab, and S.-A. Ahmadi.  V-net: Fully convolutional neural networks for volumetric medical image segmentation.  In 2016 Fourth International Conference on 3D Vision (3DV), pages 565–571. IEEE, 2016.
  * [30] D. B. Mumford and J. Shah.  Optimal approximations by piecewise smooth functions and associated variational problems.  Communications on pure and applied mathematics, 1989.
  * [31] O. Oktay, J. Schlemper, L. L. Folgoc, M. Lee, M. Heinrich, K. Misawa, K. Mori, S. McDonagh, N. Y. Hammerla, B. Kainz, et al.  Attention u-net: Learning where to look for the pancreas.  arXiv preprint arXiv:1804.03999, 2018.
  * [32] S. J. Pan and Q. Yang.  A survey on transfer learning.  IEEE Transactions on knowledge and data engineering, 22(10):1345–1359, 2009.
  * [33] O. Ronneberger, P. Fischer, and T. Brox.  U-net: Convolutional networks for biomedical image segmentation.  In International Conference on Medical image computing and computer-assisted intervention, pages 234–241. Springer, 2015.
  * [34] H. R. Roth, A. Farag, E. B. Turkbey, L. Lu, J. Liu, and R. M. Summers.  Data from pancreas-ct., 2016.
  * [35] H. R. Roth, L. Lu, A. Farag, H.-C. Shin, J. Liu, E. B. Turkbey, and R. M. Summers.  Deeporgan: Multi-level deep convolutional networks for automated pancreas segmentation.  In International conference on medical image computing and computer-assisted intervention, pages 556–564. Springer, 2015.
  * [36] H. R. Roth, L. Lu, N. Lay, A. P. Harrison, A. Farag, A. Sohn, and R. M. Summers.  Spatial aggregation of holistically-nested convolutional neural networks for automated pancreas localization and segmentation.  Medical image analysis, 45:94–107, 2018.
  * [37] H. R. Roth, H. Oda, Y. Hayashi, M. Oda, N. Shimizu, M. Fujiwara, K. Misawa, and K. Mori.  Hierarchical 3d fully convolutional networks for multi-organ segmentation.  arXiv preprint arXiv:1704.06382, 2017.
  * [38] H.-C. Shin, H. R. Roth, M. Gao, L. Lu, Z. Xu, I. Nogues, J. Yao, D. Mollura, and R. M. Summers.  Deep convolutional neural networks for computer-aided detection: Cnn architectures, dataset characteristics and transfer learning.  IEEE transactions on medical imaging, 35(5):1285–1298, 2016.
  * [39] P. Y. Simard, D. Steinkraus, J. C. Platt, et al.  Best practices for convolutional neural networks applied to visual document analysis.  In Icdar, volume 3, 2003.
  * [40] A. L. Simpson, M. Antonelli, S. Bakas, M. Bilello, K. Farahani, B. Van Ginneken, A. Kopp-Schneider, B. A. Landman, G. Litjens, B. Menze, et al.  A large annotated medical image dataset for the development and evaluation of segmentation algorithms.  arXiv preprint arXiv:1902.09063, 2019.
  * [41] G. Sundaramoorthi, A. Yezzi, and A. C. Mennucci.  Sobolev active contours.  International Journal of Computer Vision, 73(3):345–366, 2007.
  * [42] M. Tan and Q. V. Le.  Efficientnet: Rethinking model scaling for convolutional neural networks.  arXiv preprint arXiv:1905.11946, 2019.
  * [43] A. Tsai, A. Yezzi, and A. S. Willsky.  Curve evolution implementation of the mumford-shah functional for image segmentation, denoising, interpolation, and magnification.  IEEE transactions on Image Processing, 10(8):1169–1186, 2001.
  * [44] Z. Yang, L. Zhang, M. Zhang, J. Feng, Z. Wu, F. Ren, and Y. Lv.  Pancreas segmentation in abdominal ct scans using inter-/intra-slice contextual information with a cascade neural network.  In 2019 41st Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pages 5937–5940. IEEE, 2019\.
  * [45] A. Yezzi, S. Kichenassamy, A. Kumar, P. Olver, and A. Tannenbaum.  A geometric snake model for segmentation of medical imagery.  IEEE Transactions on medical imaging, 16(2):199–209, 1997.
  * [46] Q. Yu, L. Xie, Y. Wang, Y. Zhou, E. K. Fishman, and A. L. Yuille.  Recurrent saliency transformation network: Incorporating multi-stage visual cues for small organ segmentation.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 8280–8289, 2018.
  * [47] N. Zhao, N. Tong, D. Ruan, and K. Sheng.  Fully automated pancreas segmentation with two-stage 3d convolutional neural networks.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 201–209. Springer, 2019.
  * [48] H. Zheng, Y. Chen, X. Yue, C. Ma, X. Liu, P. Yang, and J. Lu.  Deep pancreas segmentation with uncertain regions of shadowed sets.  Magnetic Resonance Imaging, 68:45–52, 2020.
  * [49] Y. Zhou, L. Xie, W. Shen, Y. Wang, E. K. Fishman, and A. L. Yuille.  A fixed-point model for pancreas segmentation in abdominal ct scans.  In International conference on medical image computing and computer-assisted intervention, pages 693–701. Springer, 2017.
  * [50] S. C. Zhu and A. Yuille.  Region competition: Unifying snakes, region growing, and bayes/mdl for multiband image segmentation.  IEEE transactions on pattern analysis and machine intelligence, 18(9):884–900, 1996.
  * [51] Z. Zhu, C. Liu, D. Yang, A. Yuille, and D. Xu.  V-nas: Neural architecture search for volumetric medical image segmentation.  In 2019 International Conference on 3D Vision (3DV), pages 240–248. IEEE, 2019.
  * [52] Z. Zhu, Y. Xia, W. Shen, E. Fishman, and A. Yuille.  A 3d coarse-to-fine framework for volumetric medical image segmentation.  In 2018 International Conference on 3D Vision (3DV), pages 682–690. IEEE, 2018.