# The Silicon Electron Multiplier Sensor11footnotemark: 1

Marius Mæhlum Halvorsen  Victor Coco  Evangelos Leonidas Gkougkousis  Paula
Collins  Olivier Girard  CERN EP-LBD, Esplanade des Particules 1, Meyrin,
1211, Switzerland  University of Oslo, Department of Physics, Sem Sælands vei
24, Oslo, 0315, Norway  CEA-Leti, 17 avenue de Martyrs, Grenoble, 38054,
France

###### Abstract

The Silicon Electron Multiplier (SiEM) is a novel sensor concept for minimum
ionizing particle (MIP) detection which uses internal gain and fine pitch to
achieve excellent temporal and spatial resolution. In contrast to sensors
where the gain region is induced by doping (LGADs, APDs), amplification in the
SiEM is achieved by applying an electric potential difference in a composite
electrode structure embedded within the silicon bulk using MEMS fabrication
techniques. Since no gain-layer deactivation is expected with radiation
damage, such a structure is expected to withstand fluences of up to
1016neq/cm2. Various geometries and biasing configurations are studied, within
the boundaries imposed by the fabrication process being considered. The
effective gain, the field in the sensor, the leakage current and breakdown
conditions are studied for cell sizes in the range of 6 –
$15\text{\,}\mathrm{\SIUnitSymbolMicro m}$. Simulations show that gains in
excess of 10 can be achieved, and studies of the time structure of the induced
signals from a charge cloud deposited in the middle of the sensor show that
time resolutions similar to other sensors with internal gain can be expected.
Plans for the manufacture of a proof-of-concept sensor and for its subsequent
characterisation are discussed.

††journal: Nuclear Instrumentation Method A

[figure]style=plain,subcapbesideposition=top

##  1 Introduction

Charged particle detection in the innermost region of high energy physics
collider experiments combines the challenges of high fluence and high
occupancy with the need for measurements with excellent spatial and time
resolution. While demand for spatial resolutions below
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ is now commonplace, an ever
increasing number of applications now require in addition a time resolution
below 50 ps [1, 2, 3, 4, 5]. Existing solid state detector technologies suffer
various limitations. While thin planar technology allows fluences of several
times 1016cm2s-1 to be reached, and can in principle reach time resolution as
low as 30 ps [6], the limited signal amplitude makes this difficult to
exploit. High fluences can also be reached with 3D technologies, and thanks to
the decoupling between the drift direction and charge deposition direction,
the Landau fluctuation contribution to the time resolution becomes sub-
leading. Hence, time resolution below the 20-30 ps range have been measured
[7, 8]. Nonetheless, 3D technologies suffer high capacitance, fill factor
reduction and it will be difficult to reach pitches lower than
$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$ with present technologies. Another
approach consists of amplifying the signal within the sensors, such as done in
the LGADs [9] family. Gains of around 10 allow time resolutions in the range
20 – 30 ps to be reached. Recent variations of the LGAD design address the
limitation in segmentation (iLGAD [10], trenched LGAD, AC-LGADs, …). However
LGADs suffer gain loss for fluences above 2 1015neq/cm2 [11] due to acceptor
removal in the p+-doped gain layer.

We propose here a novel sensor design, the silicon electron multiplier (SiEM),
which provides intrinsic gain together with a pixel pitch better than
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$, without reliance on shallow
doping. The SiEM is expected to be more radaition hard than LGADs while
providing similar time resolution.

##  2 Sensor description

Figure 1: Schematic description of the SiEM detector. It consists of an
O($50\text{\,}\mathrm{\SIUnitSymbolMicro m}$) thick depleted silicon region
\raisebox{-.9pt} {1}⃝ adjacent to a region \raisebox{-.9pt} {2}⃝ made of
silicon pillars. In between the pillars two metallic electrodes, separated by
a dielectric are buried. The electrodes produce a high electric field region
in \raisebox{-.9pt} {2}⃝. As a particle passes through \raisebox{-.9pt} {1}⃝,
primary electrons are produced by ionisation and drift toward \raisebox{-.9pt}
{2}⃝ where they are multiplied and induce a signal while drifting towards the
readout electrode.

###  2.1 Detection principle

As shown in Figure 1 the SiEM can be delineated into two regions. In a
conversion and drift layer, primary charge carriers are created by the passage
of minimum ionising particles (MIPs) through the sensor which then drift under
the influence of an electric field towards an amplification and induction
layer. In this layer a region of high electric field is generated in the
material to achieve primary charge multiplication by impact ionisation. The
high electric field is achieved by applying electric potentials to a composite
conductive structure buried in the bulk of the material. The drift of the
secondary charges in this region induces a signal on the readout electrode.
The exact geometry of the amplification and induction regions depends on the
bulk material as well as the fabrication process that can be achieved. For the
design proposed in this paper, as illustrated in Figure 1, the bulk of the
sensor consists of a reversed biased n-in-p silicon junction. The readout side
is etched, leaving cylindrical Si-pillars which are terminated with the n+
implant and metal contact to the readout electrodes. A composite structure
consisting of two conductive grid electrodes separated by a dielectric is
deposited in the etched region. A difference of potential is applied to the
two grid electrodes creating a high field region in the silicon pillars. The
amplification is achieved in the pillars and charges drifting in this region
induces a signal on the readout electrodes. The amplification electrodes can
be biased either via bonding in a dedicated etched area, or through
metallisation running along the pillar wall in one of the etched region.
Deep reactive ion etching process (DRIE) followed by successive deposition and
patterning of dielectric and metal electrode is a possible approach to the
making of this device. Constraints associated to this fabrication technique
are considered in the following, in particular the interplay between the depth
of the etch, the thickness of the dielectric layers that will be deposited on
the wall and the distance needed between the pillar wall and the metal
electrodes to ensure no connection is left between the top and bottom of the
pillars. Alternative fabrication approaches may lift such limitations. Single
buried electrode structures may also be envisaged and will be discussed in
section 4

###  2.2 Geometry description

The text in this section will refer to Figure 2, shows a cross-section of the
SiEM design, identifying the most important dimensions. The thickness of the
drift region (D) and the pitch (p) between pillars can be optimised depending
on the application. Large thickness associated with small pitch readout could
be used to reconstruct the direction of the impinging particle while a thin
drift region is preferred for fast timing application, in order to minimise
the time distribution centroid of the signal and thus improving on time
resolution [6].

In the double buried electrodes geometry the electrode closest to the readout
electrode is numbered 1 while the one closest to the drift region is numbered
2.

Figure 2: Definition of the key geometry parameters for the etching based
fabrication process. The left schematic represents one ”pillar slice”.

The pillar width (w), the pillar height (h), and the pillar pitch (p), as well
as the electrodes geometry, their widths (w1, w2) and the inter-electrode
distance (d) are the principal drivers which define the electric field shape
and its intensity. A few other geometrical parameters are shown in Figure 2
and are related to the expected process limitation. The distances g1 and g2
are the guard distances between the electrodes and the pillar wall while td
designate the thickness of the dielectric between the bottom of the bulk and
the electrode 2 which is the first electrode to be deposited in the considered
fabrication process.
In the following sections the behaviour of the device and the impact of these
dimensions on the behaviour of the device is investigated using simulations.

##  3 Sensor behaviour

###  3.1 Setup

The TCAD device simulator Synopsis Sentaurus [12, 13] version S-2021.06 has
been used to perform the simulation of the proposed structure. The 2D-mesh has
been made by the Sentaurus Structure Editor with mesh cell sizes ranging from
200 to $10\text{\,}\mathrm{nm}$. The PARDISO solver [14, 15] along with the
Canali mobility model [16], Slotboom band gap model [17, 18, 19], and the
vanOverstreaten - de Man impact ionization model [20] have been applied for
quasi-static simulations. In transient simulations, the HeavyIon charge
deposition is used.

Unless otherwise specified, the values for the geometrical dimensions defined
in 2.2 are the following: the silicon parameters are defined by
D=$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$,
p=$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$,
h=$6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
w=$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ while the electrode geometry is
given by d=$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
td=$0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$. A slice of the sensor
corresponding to one single pillar, as shown on the left of Figure 2, is
referred to a ”cell”.

The geometry description in simulation reflects the expected constraints
coming from the expected production process. In particular it is expected that
when a layer of silicon oxide is deposited in the trenches, 60% of the oxide
thickness is also deposited along the wall. It is also expected that to
prevent metal remaining on the pillar walls, a guard of at least
$0.25\text{\,}\mathrm{\SIUnitSymbolMicro m}$ is introduced between the edge of
the electrodes and the pillar wall. The minimal value of g2 depends on td,
while the minimal value of g1 depends both on d and td.

In the default configuration the electrode 1 is as large as possible, i.e.
w1=$5.7\text{\,}\mathrm{\SIUnitSymbolMicro m}$ while the electrode 2 is
retracted from its maximal size in order to have w2=w1 and thus provide a
simple to interpret field geometry to start with.
The default voltage configuration used for the electrode biasing is the
following: while the readout electrode is at ground level, the potential of
electrode 1, V1, is set at $-5\text{\,}\mathrm{V}$ and V2, the potential of
electrode 2 can vary to produce the electric field needed for amplification
with V2=V1-$\Delta$V. The potential of the backside electrode
V${}_{\text{bias}}$ is constantly $-30\text{\,}\mathrm{V}$ below V2 to ensure
the depletion of the structure.

###  3.2 Electric field and leakage current

In Sentaurus TCAD quasi-stationary simulations, the voltages V1, V2 and Vbias
are ramp up simultaneously and linearly. During the voltage sweep, the leakage
current at each step is evaluated, leading to a characteristic curve similar
to the I-V of a planar sensor, even though the device is not bipolar. When the
target voltages are reached, the electric field and other figures of merit
discussed in 3.4 are evaluated. Figure 3 shows the I-V characteristic of a
cell for a $\Delta$V ramp up to a potential of $235\text{\,}\mathrm{V}$.

Figure 3: I-V characteristic of the device with default geometry and biasing
scheme. The current is estimated for a cell which correspond to a volume of
$512\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{3}$

At low voltages, the I-V curve follows a reversed biased diode-like shape,
until it reaches the saturation leakage current. At this point both the pillar
and the drift region are depleted and the leakage current is low enough to
allow primary charge particle detection ($\sim 0.1$pA per cell). As $\Delta$V
increases, the field eventually gets large enough for the creation of
secondary charges through impact ionisation. The leakage current increases as
the thermal electrons get multiplied in the pillar. In this regime high field
is concentrated at the top of the pillar as illustrated in Figure 4. Values of
the electric field between 200 and
$300\text{\,}\mathrm{k}\mathrm{V}\mathrm{/}\mathrm{c}\mathrm{m}$ can be
achieved in the silicon while keeping the field in the silicon oxide below
$3\text{\,}\mathrm{M}\mathrm{V}\mathrm{/}\mathrm{c}\mathrm{m}$.

(a)

(b)

Figure 4:  Maximum electric field amplitude reached just before the sensor
enters breakdown ($\Delta$V =$180\text{\,}\mathrm{V}$) (a) and the current
density just after breakdown ($\Delta$V =$195\text{\,}\mathrm{V}$) (b). On
figure (a) the black lines represent the electron velocity streamlines.

For the default pillar width of $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$,
above $\Delta$V of $200\text{\,}\mathrm{V}$ a sustained avalanche develops
along the pillar wall, where the field is the highest. This avalanche, which
leads to the breakdown of the device, is illustrated in Figure 4b.

Figure 5: Leakage current and maximal electric field at the center of the
pillar as function of $\Delta$V for different pillar width. Saturation of the
electric field amplitude is due to polarisation effects that occurs when
current flows in breakdown conditions.

The role played by w, the pillar width, was investigated and is illustrated in
Figure 5. In these simulations V1 is kept constant, while $\Delta$V is varied
for different pillar widths. The I-V curves display similar behaviour to the
one of Figure 3, but a higher $\Delta$V is required to reach charge
multiplication for thicker pillars. This is due to the fact that for the same
potential difference applied to the amplification electrodes, a lower electric
field amplitude is achieved in the center of the pillar. Furthermore, the
electric field just before breakdown voltage is reached is lower for thicker
pillars since the breakdown is initiated at the pillar wall where the field is
the highest. This reduces the expected maximum amplification capacity of
thicker pillars. With the DRIE based process, widths of down to
w$=$$0.6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ can most likely be achieved. A
pillar size of $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ is used in the rest
of the studies reported here. Note that, larger pillars are also of interest
since, as will be discussed in 3.4, other parameters than the maximum electric
field in the center of the pillar are relevant to the amplification mechanism.

###  3.3 Signal amplification

In the previous paragraph it was shown that fields above
$200\text{\,}\mathrm{k}\mathrm{V}\mathrm{/}\mathrm{c}\mathrm{m}$ can be
created in the pillar without breakdown of the device. Primary charge should
thus be multiplied by impact ionisation in the pillar, providing internal gain
to the sensor. Transient simulations are used to probe the multiplication
mechanism. A charge cloud of 80 electron–hole pairs is deposited in the drift
region center, and propagated in the device until all charges are collected at
the electrodes. During the charge transport, the induced currents at the
electrodes are recorded and further processed in order to extract the gain and
the signal shape. The gain is defined as the ratio of integrated current on
the readout electrode by the input charge. The integrated charge at the
readout electrode is leakage current subtracted.

(a)

(b)

Figure 6: Current induced on the readout electrode (a) and the backside
electrode (b) for $\Delta$V between $10\text{\,}\mathrm{V}$ and
$160\text{\,}\mathrm{V}$. $\delta$t indicates the time needed for the
electrons to reach region \raisebox{-.9pt} {2}⃝ for a $\Delta$V of
$160\text{\,}\mathrm{V}$.

In Figure 6a, the currents from transient simulations using different
$\Delta$V are displayed. Although the signal induced at the readout electrode
in Figure 6a is the one of interest for MIP detection, the signal induced on
backside electrode is interesting to understand the signal formation and is
show in Figure 6b.
At $\Delta$V =$10\text{\,}\mathrm{V}$ the device exhibits no gain. On the
backside electrode the primary holes induce a signal from t=0 ns until all
holes have reached the electrode. This signal mainly induced by the holes
drifting towards the backside electrode looks very much like what would be
observed on a pad sensor, as suggested by the distribution of the weighting
potential of the backside electrode shown on Figure 7b. On the read-out
electrode however, only the charge moving in the pillar induces signal. The
signal only start when the electrons reach the buried electrodes. Charges
moving in the drift region are shadowed by the amplification electrode as
expected from the low weighting potential variation seen in Figure 7a.
At $\Delta$V $>$$100\text{\,}\mathrm{V}$, when the electrons reach the top of
the pillar, where the electric field is the highest, they generate secondary
charges. Secondary holes drift back to the readout electrode inducing the
second wave of signal starting at $\delta$t. The amplitude of this second wave
depends on the amount of charge multiplication. The gain can be clearly seen
from looking at the amplitude of the signal induced on the readout electrode,
where the total secondary electron signal is added to the one from the primary
electrons. In this configuration, the gain vary from $\sim$2 for $\Delta$V
=$140\text{\,}\mathrm{V}$ to slightly above 10 for $\Delta$V
=$180\text{\,}\mathrm{V}$. It can be visualised on the Figure 6a by comparing
the area under the curves at $\Delta$V =$140\text{\,}\mathrm{V}$ and $\Delta$V
=$180\text{\,}\mathrm{V}$ with the area under the $\Delta$V
=$10\text{\,}\mathrm{V}$ in which the field is too low for multiplication to
occur.

(a)

(b)

Figure 7: Weighting potentials of the readout electrode (a) and backside
electrode (b). The representation is rotated by 90∘ with respect to the other
figures.

(a)

(b)

Figure 8: Electric field in pillar for different biasing configurations

###  3.4 Optimisation of the device

(a)

(b)

(c)

(d)

Figure 9: The maximum field in silicon and silicon oxide,
E${}_{\mathrm{Si~{}max}}$ (a), and E${}_{\mathrm{SiO_{2}~{}max}}$ (b)
evolution in the (V1,$\Delta$V) phase-space is presented, together with the
evolution of the field barycenter, $<\mathrm{Z}>_{\mathrm{E}}$/h (c) and the
average value of the field in the pillar, $<\mathrm{E}>$ (d). The gain is
overlaid in each figure. Simulation were performed on about 60 points of the
(V1,$\Delta$V) phase space and cubic 2D-interpolation is used to obtained a
continuous representation.

In the previous paragraphs, the configuration, both in terms of electrode
geometry and in terms of biasing scheme is close to that of a typical GEM
detector [21], with the high field region located in the interstice of the
multiplication electrodes structure as illustrated in Figure 8a. In the
following, configurations where the high field region goes further into the
pillar are investigated. In particular, high field can be generated by a
difference of potential between electrode 1 and the readout electrode. This
configuration could be achieved with a single electrode and, as illustrated in
Figure 8b, the highest field is shifted towards the readout electrode. A well
chosen balance between those two configurations would, on one hand, prevent
generation of a localised high field which could initiate breakdown either in
the silicon or in the silicon oxide and, on the other hand, increase the
fraction of the pillar volume where charge multiplication can occur. This can
be achieved by changing the electrode geometry or by modifying the biasing
scheme, as discussed in the following paragraphs.

#### Figure of merit

In general, gains of above a factor 10 are targeted. It is nonetheless useful
to define additional variables in order to understand the advantages or
limitations of each configuration. The maximum electric field in the silicon,
E${}_{\mathrm{Si~{}max}}$, and in the silicon oxide,
E${}_{\mathrm{SiO_{2}~{}max}}$, indicate if a configuration is close to a
breakdown conditions. High values of average electric field in the pillar,
$<\mathrm{E}>$ reflect higher field achieved over a larger region of the
pillar, especially when E${}_{\mathrm{Si~{}max}}$ is not large. Finally the
barycenter position of the electric field along the pillar
axis,$<\mathrm{Z}>_{\mathrm{E}}$, give information on the position of the high
field and how localised the field is. The closer it is to one electrode the
more peaked it is, the closer to the pillar center, the more spread along the
pillar it is. For this variable not to depend too much on the geometry, it is
normalised by the height of the pillar, h. Ideal configurations have high
$<\mathrm{E}>$, a medium $<\mathrm{Z}>_{\mathrm{E}}$ and as low
E${}_{\mathrm{Si~{}max}}$, and E${}_{\mathrm{SiO_{2}~{}max}}$ as possible.

#### Interplay between V1 and $\Delta$V

The impact of the biasing scheme on the aforementioned figures of merit is
studied in the following. $\Delta$V is scanned in steps of
$5\text{\,}\mathrm{V}$ until breakdown and for different values of V1. The
resulting gain and figures of merit are reported in Figure 9. As expected it
is possible to achieve gain above a factor 10 with different balances of
V1-$\Delta$V, nonetheless the way amplification is achieved in the structure
is very different. If $\Delta$V is large while V1 is small, high electric
field amplitudes are reached both in the silicon pillar, as shown in Figures
9a, and in the dielectric, as shown in Figure 9b. Depending on the quality of
the oxide this could lead to breakdown in the dielectric. When the gain is
obtained by high difference of potential between the readout electrode and the
electrode 1, ie. V1 is large and $\Delta$V small, the field in the silicon
reaches very high values close to the readout electrode and becomes the most
likely cause of breakdown of the structure. In contrast, for V1 values between
$50\text{\,}\mathrm{V}$ and $75\text{\,}\mathrm{V}$, for which gains of above
a factor 10 can be achieved with settings of $\Delta$V above
110-$60\text{\,}\mathrm{V}$. These settings limit the field in the silicon
below 250-$300\text{\,}\mathrm{k}\mathrm{V}\mathrm{/}\mathrm{c}\mathrm{m}$ and
in the silicon oxide below
3-$1.7\text{\,}\mathrm{M}\mathrm{V}\mathrm{/}\mathrm{c}\mathrm{m}$) which
would provide a more robust operation point further away from breakdown. The
field barycenter position along the pillar height in Figure 9c is, as
expected, the closest to the readout electrode for large V1, close to the
amplification electrodes in the default biasing scheme, but more central with
the ”robust operation” scheme which indicates that the pillar is better filled
with electric field. This is further supported by the fact that at this point
the average electric field shown in Figure 9d is high.

#### Electrode geometry

A similar distribution of the regions of high electric field amplitude all
along the pillar can be achieved by modifying the electrode geometry. When
retracting the electrode 1, ie. w${}_{1}<$ w2, the lower part of the pillar
receives a contribution from the difference of potential between the electrode
1 and the readout electrode. In Figure 10, the gain as function of $\Delta$V
is shown for various geometries of the amplification electrodes, from the
default geometry where w1=w2 to geometries with ever smaller dimension of w1
($<w_{2}$) down to the complete removal of electrode 1. The more retracted
electrode 1 is, the lowest $\Delta$V needs to be in order to achieve the same
gain. This is due to the fact that the smaller w1 is the more of the pillar is
filled with high electric field, even though the maximum field is lower.
This approach can be combined with modification to the (V1,$\Delta$V) balance.
The barycenter $<\mathrm{Z}>_{\mathrm{E}}$ and gain variation as function of
the biasing configuration are shown in Figure 11. As expected lower $\Delta$V
allow to reach the high gain region, and overall larger gains can be achieved.

Figure 10: Gain as function of $\Delta$V for different widths of electrode 1.
The electrode 1 width is downsized from the default
w1=w2=$5.7\text{\,}\mathrm{\SIUnitSymbolMicro m}$ down to its complete
removal. The points in light color correspond to values of $\Delta$V where the
structure is in breakdown. Figure 11: The evolution of the field barycenter,
$<\mathrm{Z}>_{\mathrm{E}}$/h and the gain as function of the biasing
configuration (V1,$\Delta$V) for a geometry with
w1=w2-$3\text{\,}\mathrm{\SIUnitSymbolMicro m}$. Simulations were performed on
60 points of the (V1,$\Delta$V) phase space and cubic 2D-interpolation is used
to obtain a continuous representation.

#### Inter-electrode spacing

In the default configuration, an inter-electrode spacing of
d=$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ is considered. In the DRIE based
process, after the patterning of the electrode 2 metal layer, silicon oxide is
deposited before depositing electrode 1. During this phase, about 60-70% of
the silicon oxide thickness d will also be deposited along the pillar wall.
Thus, for larger values of d, g1 is also larger. The variation of the gain
with increasing d is studied. In order not to mix effects, and although it was
shown in the previous paragraph that w${}_{1}<$w2 can be beneficial, we keep
w1=w2 while varying d, i.e. w1 is determined by g1 and the electrode 2 is
shortened to have the same length than electrode 1. The height of the pillar
has been modified for each inter-electrode spacing in order to keep a constant
distance, of $7\text{\,}\mathrm{\SIUnitSymbolMicro m}$, between the readout
electrode and electrode 1. This was done to limit the dependency of the study
on factors other than the spacing of the amplification electrodes. The gain
evolution with d is shown as function of $\Delta$V in Figure 12. Since the
inter-electrode distance is larger, a larger $\Delta$V is need to achieve the
same gain. On the other hand the breakdown voltage is higher and the highest
gain that can be achieved is larger with larger d as a bigger volume of the
pillar is filled with high electric field amplitude. Gains of up to a factor
40 are observed for d=$3\text{\,}\mathrm{\SIUnitSymbolMicro m}$ at $\Delta$V
=$320\text{\,}\mathrm{V}$. Together with optimised retraction of the electrode
and/or different V1/$\Delta$V balance, higher gains are likely to be
achievable.

Figure 12: Gain as a function of $\Delta$V for different inter-electrode
spacing d. w1 varies together with d due to the expected process constraints
and are reported in for completion. For the
d=$3\text{\,}\mathrm{\SIUnitSymbolMicro m}$ geometry, the gain obtained from
alternative Garfield++ simulation can be compared to the one obtained from
Sentaurus TCAD.

###  3.5 Optimisation for time resolution

While full simulations of the deposition, drift and amplifications of primary
charge carrier from MIPs, are needed to estimate the true time resolution of
the sensor, the impact of some of the sensor parameters on the time structure
of the signal can be assessed already. The time distribution of the signal
induced on the readout electrode by the secondary electrons is driven by the
time of arrival of the primary electrons at the amplification region.
Following the arguments for sensor with internal gain in [6] and assuming the
time for the secondary electrons to drift through the pillar is small222the
thickness of region \raisebox{-.9pt} {2}⃝ is small with respect to the
thickness of region \raisebox{-.9pt} {1}⃝ and the electron velocity in
\raisebox{-.9pt} {2}⃝ is expected to be close to saturation. , the time
resolution of SiEM will be dominated by the time it takes for electrons to
drift through region \raisebox{-.9pt} {1}⃝. Following [6], the centroid time
standard deviation will thus be given by
$\Delta_{t}=\omega\frac{\mathrm{T_{e}}}{\sqrt{12}}$ where $\omega$ is a factor
that depends on the charge deposition fluctuation and the sensor thickness
($\omega\sim$0.3 in this configuration) and $\mathrm{T_{e}}$ is the time it
takes electrons to drift through \raisebox{-.9pt} {1}⃝. For Vbias
=V2-$30\text{\,}\mathrm{V}$, TCAD simulation gives $\mathrm{T_{e}}\sim$
$500\text{\,}\mathrm{ps}$ thus $\Delta_{t}\sim$ $40\text{\,}\mathrm{ps}$. But
$\mathrm{T_{e}}$ can be further reduced by reducing region \raisebox{-.9pt}
{1}⃝ thickness or by operating at larger electric field in \raisebox{-.9pt}
{1}⃝ to increase the electron velocity. With gains of a factor 10 and above, a
sensor thickness of as low as $25\text{\,}\mathrm{\SIUnitSymbolMicro m}$
should provide sufficient signal for detection. For
$25\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thickness and Vbias
=V2-$30\text{\,}\mathrm{V}$, the field amplitude increase by a factor two,
resulting in a 50% increase in the electron velocity. Such a modification
would lead to $\mathrm{T_{e}}~{}\sim$ $170\text{\,}\mathrm{ps}$ thus
$\Delta_{t}\sim$ $15\text{\,}\mathrm{ps}$.
The time resolution of the sensor will also be affected by the fact that the
drift distance towards the induction region \raisebox{-.9pt} {2}⃝ will vary as
a function of the lateral position of the primary charge carrier in
\raisebox{-.9pt} {1}⃝. Indeed, as the primary charge carriers will drift along
the electron velocity streamlines visible in Figure 4a, their total drift
distance will vary as a function of their initial transverse position. The in-
homogeneity in the drift path depends on the ratio of the pitch to the pillar
width. The lower it is, the smaller the in-homogeneity. The standard deviation
of the drift path as function of the transverse position of the primary
electrons is illustrated in Figure 13. For the default
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ width and Vbias
=V2-$30\text{\,}\mathrm{V}$, it varies from
$0.3\text{\,}\mathrm{\SIUnitSymbolMicro m}$ to
$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ for pitch varying from
$6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ to
$15\text{\,}\mathrm{\SIUnitSymbolMicro m}$. The standard deviation of the
drifting time for uniform transverse position distribution then varies from
$5\text{\,}\mathrm{ps}$ to $23\text{\,}\mathrm{ps}$.
This contribution is sub-leading with respect to the contribution from the
centroid time standard deviation which may range from 15 to
$40\text{\,}\mathrm{ps}$ depending on the thickness of \raisebox{-.9pt} {1}⃝
and the electric field amplitude in this region.

Figure 13: The standard deviation of the drift distance,
$\sigma_{\mathrm{d}}$, (in blue) and of the drift time, $\sigma_{\mathrm{t}}$,
(in black) between charges spread uniformly in the transverse direction in
\raisebox{-.9pt} {1}⃝ and left to drift until \raisebox{-.9pt} {2}⃝ is shown
as function of the pitch p.

##  4 Discussion and future work

#### Impact ionisation model

In the results presented so far, the van Overstreaten - de Man impact
ionisation model [20] has been used. Several models are available and the I-V
characteristic curve shown in Figure 3 have been recalculated using the
University of Bologna [22], and the ”New” University of Bologna [23] models,
as well as the Okuto - Crowell[24] and the Lackner [25] models for impact
ionisation. The results are compared in Figure 14. The I-V characteristics and
all the curves show similar increase of the leakage current, corresponding to
charge multiplication, over about a $50\text{\,}\mathrm{V}$ range before
reaching breakdown conditions. The main difference between the various models
is the values of $\Delta$V at which charge multiplication appears, with the
van Overstreaten - de Man model predicting charge multiplication about
$20\text{\,}\mathrm{V}$ earlier than the other models.

Figure 14: IV-curves of quasi-stationary simulations for different impact
ionisation models.

#### Response to minimum ionising particles

A complete estimation of the time resolution that can be achieved by such
sensor requires on one hand to simulate the interaction of MIP with the
sensor, and on the other hand to study the impact of the front-end electronics
on the signal. The MIP simulation is performed with Garfield++ [26]. A full
study of the impact of the sensor geometry and field configuration will be the
topic of a future publication, but preliminary studies have been performed to
understand the compatibility of the gain estimated in TCAD from a charge cloud
deposited in the center of the bulk, as described in section 3.3, and the gain
estimated from the amount of charge induced on the readout electrode by a MIP
going through the detector with and without considering impact ionisation. In
Figure 12, the result of the gain evaluated with Garfield++ is compared to the
gain evaluated with the transient simulation of Sentaurus TCAD, using an
identical setup, ie. the deposit of a charge cloud in the center of region
\raisebox{-.9pt} {2}⃝. For the configuration which have been studied,
d=$3\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
w1=$3\text{\,}\mathrm{\SIUnitSymbolMicro m}$, the gain estimation with the two
methods agree to within five percent.

#### Number of electrodes

In order to provide further amplification and a better control of the field in
the pillar, and providing the guard distance between the electrode and the
pillar wall can be respected, geometries can be considered with more than two
stacked amplification electrodes. Conversely, simpler geometries with a single
electrode are also interesting, as they could be made with simpler fabrication
flow and no more overlay constraint between electrode 1 and electrode 2. As
shown in the discussion on electrode geometry, namely on Figure 10, a single
electrode geometry can indeed achieve gain. In those configurations, the
smaller the guard between the electrode and the pillar wall, and the smaller
the width over height pillar ratio, the better. It prevents highly peaked
electric field on the edge of the pillar just above the readout electrode.

#### Alternative materials

Any semi-conductor could be used for the bulk provided the band-gap is low
enough to produce significant amount of charge carriers. Whereas depleted
silicon offers more charges per MIP for low leakage current, silicon carbide
and diamond are intrinsically more radiation resistant and while less electron
holes pairs from ionisation by MIPs are expected, the amplification structure
could compensate for it.

#### Alternative processes

The default process, based on DRIE, will be studied in the coming year. In
parallel, an alternative process based on metal assisted chemical etching in
liquid or gaseous phase[27, 28], followed by plating of the etching catalyst
will be studied. It could achieve high aspect ratio pillars without
significant guard distance between the pillar wall and the electrode. This
process could be ideal for single electrode geometries. In the case of diamond
bulk, graphitisation of the diamond by laser could be used to produce the
embedded electrode as is the case in 3D sensor applications [29].

#### Possible adverse effects

It is possible that DRIE etching generates low quality silicon surface
terminations, leading to higher leakage currents. Annealing of the surface and
slow growth of a thermal oxide should allow this effect to be minimise.
Another concern is that filling the etched area with silicon oxide could
produce adverse effect after irradiation due to ionization in the oxide. To
some extent the build up of charge at the silicon interface should be
prevented by the presence of the buried amplification electrodes who may
collect the slow moving holes in the oxide before they reach the interface.
This will be further studied in simulation. Indeed, if the dielectric is
strictly needed in between electrode 1 and electrode 2 in the double electrode
configuration, filling up the rest of the etched region is expected to ease
the under bump metallisation deposition, but alternative techniques can be
investigated,
The production of a demonstrator will be a key achievement to study the impact
of those effects and assess the real performance of the device.

##  5 Conclusion

A new approach to achieve internal gain in solid state detector has been
presented. The Silicon Electron Multiplier (SiEM) combines a conversion region
where the passage of a MIP creates charge carriers, and an amplification and
induction region where the charge multiplication occurs and the moving charges
induce signal on the readout electrode. Simulations are used to assess and
optimise the relevant geometry variables. The expected performance is
estimated through simulations as well for a geometry that could be implemented
in silicon with a process based on deep reactive ion etching.

## Acknowledgment

The authors would like to thank warmly Jakob Haimberger and Heinrich Schindler
for their help with the simulations as well as Marc-Olivier Bettler for his
careful proofreading This work was carried out in the context of the CERN
Strategic R&D Programme on Technologies for Future Experiment https://ep-
rnd.web.cern.ch/.

## References

  * [1] Technical Design Report: A High-Granularity Timing Detector for the ATLAS Phase-II Upgrade, Tech. rep., CERN, Geneva (Jun 2020).
URL https://cds.cern.ch/record/2719855

  * [2] A MIP Timing Detector for the CMS Phase-2 Upgrade, Tech. rep., CERN, Geneva (Mar 2019).
URL http://cds.cern.ch/record/2667167

  * [3] Expression of Interest for a Phase-II LHCb Upgrade: Opportunities in flavour physics, and beyond, in the HL-LHC era, Tech. rep., CERN, Geneva (Feb 2017).
URL https://cds.cern.ch/record/2244311

  * [4] C. Lazzeroni, NA62/KLEVER Prospects for future high-intensity $K^{+}$ and $K_{L}$ running, including beam dump, physics Beyond Colliders Annual Workshop (2021).
URL
https://indico.cern.ch/event/1002356/contributions/4229594/attachments/2199520/3721249/PBCMar2021_Lazzeroni12_final.pdf

  * [5] G. Wilkinson, TauFV \- $\tau\rightarrow\mu\mu\mu$ and more: a fixed-target flavour-physics experiment at the Beam Dump Facility (BDF), Physics Beyond Colliders Annual Workshop (2021).
URL
https://indico.cern.ch/event/1002356/contributions/4229592/attachments/2200198/3721099/TauFV_BDF_2021.pdf

  * [6] W. Riegler, G. A. Rinella, Time resolution of silicon pixel sensors, Journal of Instrumentation 12 (11) (2017) P11017–P11017.  doi:10.1088/1748-0221/12/11/p11017.
  * [7] L. Anderlini, M. Aresti, A. Bizzeti, M. Boscardin, A. Cardini, G.-F. D. Betta, M. Ferrero, G. Forcolin, M. Garau, A. Lai, et al., Intrinsic time resolution of 3D-trench silicon pixels for charged particle detection, Journal of Instrumentation 15 (09) (2020) P09029–P09029.  doi:10.1088/1748-0221/15/09/p09029.
  * [8] G. Kramberger, et al., Timing performance of small cell 3D silicon detectors, Nucl. Instrum. Meth. A 934 (2019) 26–32.  arXiv:1901.02538, doi:10.1016/j.nima.2019.04.088.
  * [9] G. Pellegrini, P. Fernández-Martínez, M. Baselga, C. Fleta, D. Flores, V. Greco, S. Hidalgo, I. Mandić, G. Kramberger, D. Quirion, M. Ullan, Technology developments and first measurements of low gain avalanche detectors (LGAD) for high energy physics applications, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 765 (2014) 12–16.  doi:https://doi.org/10.1016/j.nima.2014.06.008.
  * [10] G. Pellegrini, M. Baselga, M. Carulla, V. Fadeyev, P. Fernández-Martínez, M. F. García, D. Flores, Z. Galloway, C. Gallrapp, S. Hidalgo, Z. Liang, A. Merlos, M. Moll, D. Quirion, H. Sadrozinski, M. Stricker, I. Vila, Recent technological developments on LGAD and iLGAD detectors for tracking and timing applications, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 831 (2016) 24–28.
  * [11] S. Mazza, E. Estrada, Z. Galloway, C. Gee, A. Goto, Z. Luce, F. McKinney-Martinez, R. Rodriguez, H.-W. Sadrozinski, A. Seiden, B. Smithers, Y. Zhao, V. Cindro, G. Kramberger, I. Mandić, M. Mikuž, M. Zavrtanik, R. Arcidiacono, N. Cartiglia, M. Ferrero, M. Mandurrino, V. Sola, A. Staiano, M. Boscardin, G. D. Betta, F. Ficorella, L. Pancheri, G. Paternoster, Properties of FBK UFSDs after neutron and proton irradiation up to $6\cdot 10^{15}\rm{n_{eq}/cm^{2}}$, Journal of Instrumentation 15 (04) (2020) T04008–T04008.  doi:10.1088/1748-0221/15/04/t04008.
  * [12] Sentaurus device, https://www.synopsys.com/silicon/tcad/device-simulation/sentaurus-device.html.
  * [13] Synopsys, SentaurusTM Device User Guide.
  * [14] O. V. Schenk, Scalable parallel sparse lu factorization methods on shared memory multiprocessors (2000).
  * [15] O. Schenk, K. Gärtner, W. Fichtner, Efficient sparse lu factorization with left-right looking strategy on shared memory multiprocessors, BIT 40 (2000) 158–176.  doi:10.1023/A:1022326604210.
  * [16] C. Canali, G. Majni, R. Minder, G. Ottaviani, Electron and hole drift velocity measurements in silicon and their empirical relation to electric field and temperature, IEEE Transactions on Electron Devices 22 (11) (1975) 1045–1047.  doi:10.1109/T-ED.1975.18267.
  * [17] J. Slotboom, H. de Graaff, Measurements of bandgap narrowing in Si bipolar transistors, Solid-State Electronics 19 (10) (1976) 857–862.  doi:10.1016/0038-1101(76)90043-5.
  * [18] J. Slotboom, H. de Graaff, Bandgap narrowing in silicon bipolar transistors, IEEE Transactions on Electron Devices 24 (8) (1977) 1123–1125.  doi:10.1109/T-ED.1977.18889.
  * [19] D. Klaassen, J. Slotboom, H. de Graaff, Unified apparent bandgap narrowing in n- and p-type silicon, Solid-State Electronics 35 (2) (1992) 125–129.  doi:https://doi.org/10.1016/0038-1101(92)90051-D.
URL https://www.sciencedirect.com/science/article/pii/003811019290051D

  * [20] R. Van Overstraeten, H. De Man, Measurement of the ionization rates in diffused silicon p-n junctions, Solid-State Electronics 13 (5) (1970) 583–608.  doi:https://doi.org/10.1016/0038-1101(70)90139-5.
  * [21] F. Sauli, GEM: A new concept for electron amplification in gas detectors, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 386 (2) (1997) 531–534.  doi:https://doi.org/10.1016/S0168-9002(96)01172-2.
URL https://www.sciencedirect.com/science/article/pii/S0168900296011722

  * [22] M. Valdinoci, Impact-ionization in silicon at large operating temperature, extraction of electron impact-ionization coefficient at large operating temperatures,, International conference on Simulation of Semiconductor Process and Devices (1999) 27–30.
  * [23] S. Reggiani, Experimental extraction of electron impact-ionization coefficient at large operating temperatures,, IEEE Transactions on ELectron Devices (2004) 407–410.
  * [24] Y. Okuto, C. Crowell, Threshold Energy Effect on Avalanche Breakdown Voltage in Semiconductor Junctions, Solid-State Electronics 18 (2) (1975) 161–168.
  * [25] T. Lackner, Avalanche multiplication in semiconductors: A modification of chynoweth’s law,, Solid-State Electronics 34 (2) (1991) 33–42.
  * [26] H. Schindler, Garfield++ User’s Guide, CERN.
URL https://garfieldpp.web.cern.ch/garfieldpp/

  * [27] L. Romano, J. Vila-Comamala, K. Jefimovs, M. Stampanoni, High-Aspect-Ratio Grating Microfabrication by Platinum-Assisted Chemical Etching and Gold Electroplating, Advanced Engineering Materials 22 (10) (2020) 2000258.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/adem.202000258, doi:https://doi.org/10.1002/adem.202000258.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/adem.202000258

  * [28] L. Romano, M. Kagias, J. Vila-Comamala, K. Jefimovs, L.-T. Tseng, V. A. Guzenko, M. Stampanoni, Metal assisted chemical etching of silicon in the gas phase: a nanofabrication platform for X-ray optics, Nanoscale Horiz. 5 (2020) 869–879.  doi:10.1039/C9NH00709A.
URL http://dx.doi.org/10.1039/C9NH00709A

  * [29] A. Oh, B. Caylar, M. Pomorski, T. Wengler, A novel detector with graphitic electrodes in CVD diamond, Industrial Diamond Review 38 (06 2013).  doi:10.1016/j.diamond.2013.06.003.


# Benefit from public unlabeled data: A Frangi filtering-based pretraining
network for 3D cerebrovascular segmentation

Gen Shi Hao Lu Hui Hui<EMAIL_ADDRESS>Jie Tian<EMAIL_ADDRESS>School of
Engineering Medicine and School of Biological Science and Medical Engineering,
Beihang University, Beijing, 100191, China, and also with the Key Laboratory
of Big DataBased Precision Medicine (Beihang University), Ministry of Industry
and Information Technology of China, Beijing, 100191, China  State Key
Laboratory for Management and Control of Complex Systems, Institute of
Automation, Chinese Academic of Science, Beijing 10086, China  CAS Key
Laboratory of Molecular Imaging, Institute of Automation, Chinese Academy of
Sciences, Beijing, 100190, China

###### Abstract

The precise cerebrovascular segmentation in time-of-flight magnetic resonance
angiography (TOF-MRA) data is crucial for clinically computer-aided diagnosis.
However, the sparse distribution of cerebrovascular structures in TOF-MRA
results in an exceedingly high cost for manual data labeling. The use of
unlabeled TOF-MRA data holds the potential to enhance model performance
significantly. In this study, we construct the largest preprocessed unlabeled
TOF-MRA datasets (1510 subjects) to date. We also provide three additional
labeled datasets totaling 113 subjects. Furthermore, we propose a simple yet
effective pertraining strategy based on Frangi filtering, known for enhancing
vessel-like structures, to fully leverage the unlabeled data for 3D
cerebrovascular segmentation. Specifically, we develop a Frangi filtering-
based preprocessing workflow to handle the large-scale unlabeled dataset, and
a multi-task pretraining strategy is proposed to effectively utilize the
preprocessed data. By employing this approach, we maximize the knowledge
gained from the unlabeled data. The pretrained model is evaluated on four
cerebrovascular segmentation datasets. The results have demonstrated the
superior performance of our model, with an improvement of approximately 3%
compared to state-of-the-art semi- and self-supervised methods. Furthermore,
the ablation studies also demonstrate the generalizability and effectiveness
of the pretraining method regarding the backbone structures. The code and data
have been open source at: https://github.com/shigen-StoneRoot/FFPN.

###### keywords:

###### MSC:

41A05, 41A10, 65D05, 65D17 Keyword1, Keyword2, Keyword3

††journal: Medical Image Analysis Fig. 1: The illustration of using unlabeled
MRA-TOF data through (a) semi-supervised learning, (b) self-supervised
learning and (c) our pretraining method.

##  1 Introduction

TOF-MRA is one non-invasive medical imaging technique to visualize cerebral
blood vessels. It provides detailed images of the blood vessels in the brain,
enabling early detection and treatment of potentially life-threatening
conditions [31, 14]. The accurate vessel segmentation is a crucial
preprocessing step in TOF-MRA image analysis, which provides surgeons with
essential information about the location, size, and connectivity of blood
vessels and facilitates surgical planning and intervention [28].

Compared to the segmentation of other biological tissues, brain vessel
segmentation based on TOF-MRA presents more challenges. Most organs typically
have a spherical shape and are distributed in a concentrated manner, while
blood vessels have a tubular structure and their distribution in TOF-MRA is
relatively sparse [4, 43], which increases the difficulty and cost in
obtaining manually labeled MRA-TOF data. On the contrary, the unlabeled TOF-
MRA data are easily accessible from the public resource (i.e., IXI and
OASIS3). Therefore, how to fully use the public unlabeled data is essential
for cerebrovascular segmentation with limited manually labeled data.

There are two fundamental deep learning-based methods that are capable of
using unlabeled data—-semi-supervised learning (SemiSL) [32, 45] and self-
supervised learning (SSL) [12, 19] methods (see Fig.1 a and b). SemiSL-based
methods typically utilize unlabeled data in conjunction with labeled data
during the training process, learning the consistency between them [9]. This
process requires repeated access to the unlabeled data and increases the
demand for computational resources, when different research institutes want to
use these data. Besides, the public unlabeled data are collected from multiple
resources and may be highly heterogeneous from the single-site dataset. The
issue of distribution shift, arising from this, has the potential to adversely
affect the performance of SemiSL methods [6]. Regarding the SSL methods, they
usually pretrain the model through the pretext task [47], and then the
pretrained model can be reused by different institutes without access to the
unlabeled data. This process decreases the computational cost when the
unlabeled data scale is large and protects patient privacy since different
institutes do not directly access the unlabeled data [2].

However, the performance of SSL methods highly depends on the choice of
pretext task. The pretext tasks, such as contrastive learning between
augmented data [7, 25] and masked imaging modelling [17], are significantly
different from the downstream task (e.g., cerebrovascular segmentation in this
study). This setting proves advantageous when dealing with multiple downstream
tasks, while a specific single downstream task might not gain much benefits
from SSL methods. Therefore, incorporating the prior knowledge of the specific
downstream task into the pretraining procedure may be extremely valuable.
Frangi filtering is an effective technique to enhance and extract tubular
structures in medical images [13]. The filtering process involves analyzing
the local intensity and Hessian matrix of the image to detect vessel-like
structures, making it useful for applications such as vessel segmentation and
analysis in medical imaging. Frangi filtering plays a fundamental role as a
preprocessing step in traditional vessel segmentation algorithms. It is worth
seriously considering incorporating this filtering technique into the
pretraining task.

Driven by the abovementioned analysis, we propose a Frangi filtering-based
pretraining network (called FFPN) to fully use the public unlabeled TOF-MRA
data. In particular, we have developed a preprocessing workflow to handle
large-scale unlabeled data. This preprocessing program utilizes Frangi
filtering to enhance vessel structures. Additionally, it incorporates
thresholding methods and connected component analysis to obtain a coarse
segmentation of the vessels. Then, the vessel-enhanced images and coarse
vessel segmentations are used for a multi-task pretraining procedure. Our
model is pretrained on the large-scale TOF-MRA dataset with 1510 volumetric
data and evaluated on four labeled datasets. The results demonstrate that it
notably outperforms the state-of-the-art SemiSL and SSL methods. Besides, the
ablation studies also confirm that our proposed pretraining task helps the
model achieve better performance with fewer labeled data. Moreover, this
pretraining task is not limited to specific model architectures and can
significantly improve the performance across various backbone structures
consistently.

The main contributions of this study can be summarized as follows:

  * 1.

We develop an automated preprocessing workflow to efficiently handle unlabeled
data. The preprocessing provides the vessel enhanced images and coarse vessel
segmentations that are used for the pretraining procedure.

  * 2.

We propose a simple yet effective pretraining strategy by using the
preprocessed TOF-MRA data. This pretraining task contains regression learning,
coarse label segmentation, and max intensity projection (MIP) consistency
learning. The results show that the proposed method shows great superiority
over existing models.

  * 3.

We constructed a large unlabeled dataset (1510 volumetric data) and three
manual-annotated datasets (a total of 113 volumetric data). These data are
expected to advance semi-supervised, weakly supervised, and self-supervised
learning methods in cerebrovascular segmentation.

Fig. 2: (a) The illustration of our proposed preprocessing workflow. (b) The
overall framework of our proposed pretraining procedure with multi-task
learning.

##  2 Related Work

###  2.1 Cerebrovascular Segmentation

Traditional methods for cerebrovascular segmentation typically rely on the
continuity of grayscale variations in images to perform threshold-based
segmentation [30, 13]. In addition, statistical modeling-based approaches
utilize the concept that different brain tissues exhibit distinct pixel
grayscale distributions during imaging [11, 1]. This has been successful in
the segmentation of cerebrovascular structures in TOF-MRA images. In recent
years, deep learning methods have shown remarkable performance in vessel
segmentation tasks [35, 38, 33]. For instance, Livne et al. utilized a UNet-
based model to segment vessels in patients with cerebrovascular disease [26].
Chen et al. proposed an attention-assisted generative adversarial network
(A-SegAN) for automatic cerebrovascular segmentation [8]. More related
research can be referred to this review [5].

###  2.2 SemiSL and SSL Methods for Medical Images

Numerous SemiSL and SSL methods have been developed to utilize unlabeled
medical data and improve model performance [41, 46, 36, 24]. For instance,
Chen et al. proposed a generative consistency for the semi-supervised (GCS)
model, which calculates the consistency of perturbed data to enhance
cerebrovascular segmentation performance [4]. Moreover, SSL-based methods have
also demonstrated the effectiveness of pretraining in medical imaging
research. To illustrate, Tang et al. introduced a multitask pretraining
strategy using a 3D Swin transformer, achieving significant success in
multiple organ segmentation tasks. For more detailed information on the topic,
comprehensive reviews [47, 21, 22] can be referred to.

##  3 Methods

In this section, we introduce the proposed FFPN model. The preprocessing
workflow and pretraining procedure are shown in Fig.2 a and b, respectively.
First, we preprocess the unlabeled TOF-MRA data to obtain the vessel enhanced
image (VEI) and coarse vessel segmentation (CVS). Subsequently, during the
pretraining procedure, the VEI and CVS are utilized as regression and
segmentation targets. Further details regarding this process will be
elaborated in the upcoming sections.

###  3.1 Preprocessing Method for Unlabeled Data

The pretraining process consists of three steps: image size cropping, vessel-
structure enhancement, and coarse segmentation.

####  3.1.1 Image size cropping:

The TOF-MRA images often contain a significant number of background pixels,
which have no relevance to the vessel segmentation task. Removing these
background pixels provides a key advantage by reducing computational costs,
especially when dealing with a large-scale unlabeled dataset. Previous
approaches have utilized non-zero masks to crop raw images [20]. However,
background pixels don’t always possess intensity values of zero. In this
study, we propose a novel method to calculate the cropping mask based on the
variation of pixel intensities along the Z-axis (slice direction).

For a given single image, denoted as $I_{i}$ with dimensions $(H,W,D)$, we
compute three projection images along the Z-axis: Average Intensity Projection
(AIP), Maximum Intensity Projection (MIP), and Intensity Variation Map (IVM).
The calculations are as follows:

$\displaystyle AIP_{i}=\mathrm{Average}(I_{i})$ (1) $\displaystyle
MIP_{i}=\mathrm{Maximum}(I_{i})$ $\displaystyle IVM_{i}=\mathrm{STD}(I_{i})$

To create the cropping mask, we employ threshold segmentation. Taking MIP as
an example, we generate a mask by retaining only the positions of the top
percentage (e.g., 35%) of pixels in the MIP image. The resulting masks, namely
the AIP mask, MIP mask, and IVM mask, are then merged using a bitwise AND
operation. Additionally, small regions containing fewer than 200 pixels are
removed. Finally, we obtain a square cropping mask that is applied to each
slice of $I_{i}$.

In summary, our proposed approach enables the generation of a cropping mask
based on the intensity variations along the Z-axis in TOF-MRA images. By
applying this mask, background pixels are effectively eliminated, resulting in
a significant reduction in computational costs.

####  3.1.2 Frangi fitering:

For 3D TOF-MRA images, the Frangi filtering involves several steps to enhance
the visibility of blood vessels. First, it computes the second-order mixed
partial derivatives for each pixel using the Hessian matrix. The Hessian
matrix is defined as:

$H=\begin{bmatrix}\frac{\partial^{2}I}{\partial
x^{2}}&\frac{\partial^{2}I}{\partial x\partial
y}&\frac{\partial^{2}I}{\partial x\partial z}\\\ \frac{\partial^{2}I}{\partial
y\partial x}&\frac{\partial^{2}I}{\partial
y^{2}}&\frac{\partial^{2}I}{\partial y\partial z}\\\
\frac{\partial^{2}I}{\partial z\partial x}&\frac{\partial^{2}I}{\partial
z\partial y}&\frac{\partial^{2}I}{\partial z^{2}}\\\ \end{bmatrix}$ (2)

where $\frac{\partial^{2}I}{\partial x^{2}}$ represents the second partial
derivative of the image intensity $I$ with respect to the $x$ coordinate, and
similarly for the other partial derivatives. Next, the eigenvalues of the
Hessian matrix, denoted as $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$,
are calculated. These eigenvalues provide information about the local
structure of the image. Specifically, the magnitude and orientation of the
eigenvalues help identify pixels that belong to blood vessels.

In general, larger positive eigenvalues correspond to the central region of
vessels, smaller positive eigenvalues represent the vessel edges, while
negative eigenvalues typically indicate irrelevant structures and backgrounds.
By analyzing the eigenvalues, pixels with vascular features are selected, and
the VEI is obtained at this step.

####  3.1.3 Coarse segmentation:

To obtain a coarse segmentation of the blood vessels, a thresholding method
and connected component analysis are applied. First, a vessel mask is created
by retaining only the positions of the top percentage (e.g., 5%) of pixels for
VEIs. This thresholding step helps separate the vessel pixels from the
background and other structures. Then, connected component analysis is
performed on the vessel mask. The largest $k$ connected regions, where $k$
could be empirically chosen between 3-5, are retained as the coarse
segmentation result of blood vessels. By combining the Frangi filtering,
thresholding, and connected component analysis, the CVS is obtained, which
provides an initial approximation of the blood vessel structure in the 3D TOF-
MRA images.

Datasets | TubeTK-42 | Brains | IXI-45 | EDEN
---|---|---|---|---
Method | Dice | clDice | HD95 | Dice | clDice | HD95 | Dice | clDice | HD95 | Dice | clDice | HD95
UNet (MICCAI’15) | 72.98 | 81.33 | 9.09 | 73.77 | 75.86 | 3.72 | 79.82 | 75.23 | 10.25 | 79.79 | 80.62 | 9.89
DTC (AAAI’21) | 69.31 | 76.35 | 25.60 | 64.36 | 64.86 | 34.72 | 70.44 | 63.48 | 63.31 | 71.65 | 65.65 | 91.93
GCS (TMI’22) | 18.11 | 33.40 | 22.50 | 10.37 | 9.28 | 24.92 | 25.05 | 21.32 | 55.10 | 61.36 | 67.85 | 28.42
GBDL (CVPR’22) | 36.40 | 55.67 | 16.16 | 39.79 | 52.84 | 63.71 | 12.46 | 20.37 | 42.43 | 51.78 | 64.48 | 11.41
SSNet (MICCAI’22) | 65.15 | 71.84 | 28.65 | 64.22 | 63.88 | 35.04 | 73.55 | 65.14 | 52.05 | 69.72 | 64.26 | 96.65
BCP (CVPR’23) | 73.36 | 81.95 | 8.33 | 77.25 | 80.31 | 2.88 | 81.37 | 80.64 | 6.17 | 86.44 | 89.41 | 2.71
PCRL (ICCV’21) | 72.90 | 81.92 | 7.82 | 75.94 | 79.09 | 2.74 | 83.16 | 83.69 | 4.05 | 86.15 | 88.71 | 2.91
MAE (CVPR’22) | 72.39 | 79.37 | 11.82 | 77.56 | 79.70 | 3.49 | 83.21 | 80.14 | 10.04 | 84.62 | 85.90 | 6.73
Swin UNETR (CVPR’22) | 74.66 | 82.18 | 9.71 | 77.76 | 80.40 | 3.30 | 80.64 | 79.01 | 6.37 | 81.23 | 81.87 | 12.90
UniMiSS (ECCV’22) | 70.54 | 82.59 | 9.03 | 77.63 | 79.62 | 3.03 | 85.45 | 83.49 | 3.83 | 88.02 | 88.51 | 2.91
GVSL (CVPR’23) | 74.27 | 82.57 | 8.30 | 77.87 | 79.68 | 3.16 | 85.92 | 84.82 | 3.35 | 85.15 | 89.39 | 2.61
FFPN (Ours) | 77.53 | 85.71 | 6.69 | 80.47 | 83.28 | 2.19 | 87.24 | 86.25 | 2.75 | 90.75 | 91.53 | 1.24
Preprocessing | 29.64 | 31.27 | 71.29 | 31.73 | 37.01 | 38.55 | 43.18 | 55.50 | 38.56 | 68.48 | 69.43 | 46.18
Table 1: The cerebrovascular segmentation results on TubeTK-42, Brains, IXI-45
and EDEN datasets

###  3.2 Pretraining task and Loss Function

Our proposed pretraining process contains three tasks: regression learning,
coarse vessel segmentation learning, and consistency learning in the MIP
projection. Let function $f(\cdot)$ be the backbone neural network, and
$V_{i}$ and $S_{i}$ are the VEI and CVS for subject $i$. $f(\cdot)$ is used to
encode the raw images, and one regression header $Rgn(\cdot)$and
$Seg(\cdot)$segmentation header are used to predict $V_{i}$ and $S_{i}$,
respectively. It can be formulated as:

$\displaystyle\hat{V}_{i}=Rgn(f(I_{i}))$ (3)
$\displaystyle\hat{S}_{i}=Seg(f(I_{i}))$

where $Rgn(\cdot)$and $Seg(\cdot)$ are convolution operation with kernel size
1. $\ell_{1}$ is used to compute regression loss. For segmentation loss, Dice
loss and Binary Cross Entropy are used. These can be formulated as:

$\mathcal{L}_{rgn}=\|\hat{V}_{i}-V_{i}\|_{1}$ (4)
$\mathcal{L}_{seg}=\mathrm{Diceloss}(\hat{S}_{i}-S_{i})+\mathrm{BCEloss}(\hat{S}_{i}-S_{i})$
(5)

In this study, we do not conduct skull stripping for the preprocessed data. In
such case, the Frangi filtering will also enhance the skull pixels, while the
CVS will not contain the skull pixels after connected component analysis. In
order to focus more on learning the vascular structures, we incorporate an
extra consistency loss by calculating the $\ell_{1}$ loss between
$\hat{S}_{i}\odot\hat{V}_{i}$ and $\hat{S}_{i}\odot I_{i}$ after applying
maximum intensity projection along Z-axis. It can be formulated as:

$\displaystyle\mathrm{MIP}_{VS}=\mathrm{Maximum}(\hat{S}_{i}\odot\hat{V}_{i})$
(6)
$\displaystyle\mathrm{MIP}_{IS}=\mathrm{Maximum}(\hat{S}_{i}\odot\hat{I}_{i})$
$\displaystyle\mathcal{L}_{consistency}=\|\mathrm{MIP}_{VS}-\mathrm{MIP}_{IS}\|_{1}$

where $\hat{S}_{i}\odot\hat{V}_{i}$ and $\hat{S}_{i}\odot I_{i}$ have almost
removed the skull pixels. Meanwhile, the MIP images of TOF-MRA are commonly
used for the initial assessment of vascular morphology. The consistency loss
will help the model focus more on learning the vascular structures.

Finally, the total loss function $\mathcal{L}$ will be:

$\mathcal{L}=\gamma_{1}\mathcal{L}_{rgn}+\gamma_{2}\mathcal{L}_{seg}+\gamma_{3}\mathcal{L}_{consistency}$
(7)

where $\gamma_{1},\gamma_{2}$ and $\gamma_{3}$ are the loss weight. We
empirically set $\gamma_{1}=0.4,\gamma_{2}=0.4$ and $\gamma_{3}=0.2$ in this
study.

Fig. 3: Statistical validation on Brains dataset by using the T-test. ’*’
indicates $p<0.05$ and ’**’ indicates $p<0.001$.

###  3.3 Backbone

We select the Swin UNETR [15] as the backbone in this study. Different from
the vanilla model, successive convolutional layers are used in the initial
embedding layers [48]. This setting enhances the capture of spatial location
information, a feature that is potentially beneficial to the segmentation of
cerebrovascular structures with sparse distribution. Please note that we also
examined the impact of our pretraining strategy on other backbone structures
in the result section.

Fig. 4: The visualization comparison results of our method and the competing
baseline models from the 3D view.

##  4 Experimental Setup

###  4.1 Datasets

Unlabeled Datasets for Pretraining: We collect a large unlabeled TOF-MRA
consisting of 1510 subjects from 5 public datasets, including IXI (525
subjects), OASIS3 (525 subjects) [23], BrainAneurysm (280 subjects) [10], ADAM
(113 subjects) and TubeTK (67 subjects). We preprocessed these data by our
proposed workflow and the model is pretrained through the proposed tasks.

Labeled Datasets for Downstream Evaluation: We evaluate our proposed method in
four labeled datasets.

  * 1.

TubeTK-42: The full TubeTK dataset contains 109 subjects, with 67 unlabeled
subjects used for pretraining and 42 labeled subjects used for finetune
validation. The image size is (448, 448, 128) and the spacing is (0.5134,
0.5134, 0.8000).

  * 2.

Brains: This dataset contains 56 subjects. The size of the images is not
identical, with average size of (286, 320, 204). The spacing is (0.6188,
0.6188, 0.6200).

  * 3.

IXI-45: This dataset is provided by [8]. It contains 45 subjects from the
whole IXI dataset (570 subjects). The image size is (1024, 1024, 92). The
provided images have been processed and therefore do not contain spacing
information. The original spacing of the IXI dataset is (0.2637, 0.2637,
0.8000).

  * 4.

EDEN: This dataset contains 15 subjects. Unlike the previous dataset, we also
provide the blood vessels on the skull in this dataset. The image size is
(672, 672, 210) and the spacing is (0.2976, 0.2976, 0.4500).

Please note that the pretraining dataset and evaluation dataset are completely
isolated. For the labeled datasets, 80% subjects are used for training, and
the rest data are used for testing. Detailed information about the dataset can
be found in the Appendix.

###  4.2 Evaluation Metrics

Three metrics are used to assess the cerebrovascular segmentation performance,
including the Dice similarity coefficient (Dice), Hausdorff Distance 95%
(HD95) and centerline Dice coefficient (clDice) [37]. The Dice and HD95
metrics primarily focus on the overlap degree and shape differences in
segmentation results, clDice takes into account the topological structure of
the segmented regions, making it more suitable for evaluating tubular
structures.

###  4.3 Baseline Models and Implementation

We compare the performance of 3D cerebrovascular segmentation between the
proposed FFPN model and existing methods using unlabeled data, including the
SemiSL methods BCP [3], GBDL [40], GCS [4], SS-Net [42], DTC [27] and the SSL
methods MAE [17], PCRL [49], UniMiSS [44], Swin UNETR [39], GVSL [18].
Besides, we also include the basic UNet [34] model that does not utilize the
unlabeled data to provide a comparison.

For our proposed FFPN model, the encoder has 4 stages which comprise of [2, 2,
6, 2] transformer blocks at each stage. The initial patch size is 2 and the
feature size is 24. In pretraining task learning, the model is trained using
AdamW optimizer for 200k iterations. Warm-up cosine scheduler of 20k
iterations is used. The cropping patch size is (160, 160, 64) and the batch
size is 2. The weight decay is $3\times 10^{-5}$ and learning rate is $5\times
10^{-3}$. In finetune procedure, the model is trained using AdamW optimizer
for 5k iterations. The learning rate is $5\times 10^{-3}$ (half the learning
rate for the encoder) and the cropping patch size is (192, 192, 96).

##  5 Experimental Results

###  5.1 Comparison with state-of-the-art Methods

The cerebrovascular segmentation results for the four datasets are shown in
Table 1. The proposed method achieves the best performance compared with other
baseline models on the four datasets in terms of the three metrics. In
TubeTK-42 and Brains datasets, the clDice scores of our proposed method are
85.71 and 83.28, with an improvement of approximately 3% over the second
model. In IXI-45 and EDEN datasets, our proposed method also achieves an
improvement of around 2% compared with the second-best model.

To further demonstrate the superiority of FFPN, a statistical analysis is
shown in Fig. 3. Specifically, We conduct another 5-fold experiment on the
largest dataset (Brains) using FFPN and baseline models with competitive
performance (GVSL, Swin UNETR and PCRL). The results show that the improvement
of FFPN is statistically significant (T-test) among the three metrics (the
outcome for the metric clDice aligning closely with that of Dice).

We provide a 3D visualization result to support an intuitive evaluation (see
Fig. 4). Compared with other methods, the segmentation produced by FFPN
reveals marked superiority in achieving continuity for the vessel structure,
particularly noticeable at locations 1, 2, 3 and 4. Furthermore, it is worth
highlighting that the FFPN-enabled segmentation leads to a reduction in false
positive points (e.g., locations 5 and 6).

Fig. 5: Variation of the effect of the pretraining strategy with the labeled data percentage. Loss Function | Dice | clDice | HD95
---|---|---|---
Scratch | 72.08 | 80.16 | 11.46
$\mathcal{L}_{rgn}$ | 72.34 | 81.11 | 9.62
$\mathcal{L}_{seg}$ | 73.46 | 81.58 | 9.10
$\mathcal{L}_{consistency}$ | 72.15 | 80.78 | 10.13
$\mathcal{L}_{rgn}+\mathcal{L}_{seg}$ | 74.18 | 82.17 | 8.14
$\mathcal{L}_{rgn}+\mathcal{L}_{seg}+\mathcal{L}_{consistency}$ | 74.64 | 82.94 | 7.57
Table 2: The effectiveness of different loss functions.

###  5.2 Influence of the Percentage of Manual Labeled Data

We then conduct a study to evaluate the effectiveness of our pretraining
method using a reduced amount of manually labeled data. The results are
presented in Fig. 5.

It is observed that pretraining consistently improves the performance of the
model across different proportions of labeled data. When using all available
labeled data, our pretraining method achieves a 1.03% higher Dice score for
the IXI-45 dataset and 0.83% higher Dice score for the TubeTK-42 dataset
compared to a model trained from scratch. Additionally, with less than half of
the labeled data, the pretrained model demonstrates an improvement of 4.39%
for the IXI-45 dataset and 2.56% for the TubeTK-42 dataset. Furthermore, as
the availability of labeled images decreases, the performance gap between the
models trained with and without pretraining increases.

###  5.3 Effectiveness of Pretraining Loss Function

In this section, we evaluate the effectiveness of the proposed three loss
functions. Specifically, The model is pretrained using different loss
functions and then their effectiveness is assessed during the finetuning
procedure. We choose the TubeTK-42 dataset with 45.45% labeled data (i.e., 15
subjects) as an example, and the results are presented in Table 2.

For single-task pretraining, the segmentation learning achieves the best
performance with metrics Dice, clDice, and HD95 values of 73.46, 81.58, and
9.10, respectively. The regression learning also shows similar performance,
while the consistency learning does not significantly improve the model’s
performance. Moreover, when combining the regression and segmentation learning
loss functions, the model’s performance is further improved. Taking into
account all three learning tasks, the model attains the highest performance.

We also present a visualization result from the slice view in Fig. 6.
Incorporating pretraining with $\mathcal{L}_{seg}$ aids in reducing the
occurrence of false positive points in the generated output (e.g., locations 2
and 3). However, it is worth noting that due to the coarse segmentation
results potentially retaining some information from the skull, there may still
be instances where certain structures on the skull are annotated (e.g.,
location 4). Overall, by incorporating all pretraining tasks, the model
achieves improved segmentation, characterized by enhanced continuity of
vessels (see location 1) and a reduced presence of false positive points.

Fig. 6: Illustration of the effectiveness of pretraining loss functions from the slice view. The number indicates the clDice metric. Loaded component | Dice | clDice | HD95
---|---|---|---
Scratch | 80.79 | 78.08 | 8.41
Encoder | 81.71 | 79.90 | 12.08
Decoder | 81.14 | 76.22 | 11.19
Encoder + Decoder | 83.77 | 81.57 | 6.59
All | 85.19 | 83.09 | 4.43
Table 3: The effectiveness of different model component pretraining parameters
loaded.

###  5.4 Influence of Model Components for Pretraining

Different from previous SSL methods that typically only load the encoder
parameters for downstream tasks, our proposed approach loads both the encoder
and decoder parameters, as well as the segmentation header. In this section,
we explore the impact of loading different model components. We take the
IXI-45 dataset with 41.67% labeled data (i.e., 15 subjects) as an example, and
the results are presented in Table 3.

When only the encoder is loaded, there is a notable improvement in
performance, with the clDice score increasing to 79.90. Besides, when only the
decoder is loaded, there is a drop in performance compared to the encoder-
loaded model, as indicated by the clDice score of 76.22. Furthermore, when
both the encoder and decoder are loaded, there is a substantial boost in
performance. The clDice score increases to 81.57. When all components,
including the segmentation header, we observe the highest performance across
all metrics.

Dataset | EDEN
---|---
Backbone | Scratch | Finetune | Improvement
UNet | 79.79 | 84.26 | +4.47%
AttentionUNet | 78.44 | 88.20 | +9.76%
UNETR | 83.32 | 88.06 | +4.74%
Swin UNETR | 88.87 | 90.75 | +1.88%
Dataset | IXI-45
Backbone | Scratch | Finetune | Improvement
UNet | 79.82 | 85.24 | +5.41%
AttentionUNet | 81.52 | 84.62 | +3.10%
UNETR | 80.01 | 83.85 | +3.84%
Swin UNETR | 86.21 | 87.24 | +1.03%
Table 4: The effectiveness of our proposed pretraining method on various
backbone structures. The number indicates the Dice scores.

###  5.5 The Effectiveness on Different Backbone Structures

In this section, we assess the effectiveness of our proposed pretraining
learning on different backbone structures. We focus on four fundamental
medical segmentation backbones, which encompass both convolutional neural
network (CNN)-based models (UNet [34] and AttentionUNet [29]) and Transformer-
based models (UNETR [16] and Swin UNETR [15]). The corresponding results are
summarized in Table 4.

For the EDEN dataset, we observe improvements when using our proposed
pretraining learning compared to training from scratch across all backbone
architectures. Specifically, for the UNet backbone, the Dice score increases
from 79.79 to 84.26, resulting in a significant improvement of 4.47%.
Similarly, the AttentionUNet and UNETR backbones show improvements of 9.76%
and 4.74%, respectively. Even for the advanced Swin UNETR backbone, there is
still a notable 1.88% improvement.

In the case of the IXI-45 dataset, we observe similar trends. Our proposed
pretraining learning consistently leads to improvements compared to training
from scratch. The UNet backbone shows an improvement of 5.41%, while the
AttentionUNet, UNETR, and Swin UNETR backbones show improvements of 3.10%,
3.84%, and 1.03%, respectively.

###  5.6 The Robustness of Our Pretraining Method Regarding the Unlabeled
Dataset Heterogeneity

We carried out this experiment to assess the robustness of our proposed method
in relation to multi-source heterogeneous pretraining data, using the IXI
dataset as an illustrative case.

In particular, we maintain the quantity of unlabeled data constant (525
subjects) while varying the proportion of the unlabeled IXI dataset used. Then
we test the model performance on the labeled IXI-45 dataset. The best-
performing semiSL method BCP is used for comparison and The outcomes are
depicted in Figure 7. From these results, it is observed that the performance
of the semi-supervised learning method, BCP, significantly diminishes as the
number of homogeneous data decreases. Conversely, our proposed method
consistently shows a stable and high segmentation performance. This result
demonstrates that the proposed pretraining strategy may maintain a high level
of robustness.

Fig. 7: The robustness of our proposed method. Horizontal coordinates indicate
the composition of unlabeled data..

##  6 Discussion

To unleash the potential of unlabeled TOF-MRA data for 3D cerebrovascular
segmentation, our approach incorporates Frangi filtering as a key element in
the pretraining phase. The proposed FFPN model demonstrates enhanced
robustness and superior segmentation performance, notably outshining other
SemiSL and SSL methodologies.

Our proposed cropping strategy can reduce the computation cost when processing
such a large 3D MRA dataset. We take IXI dataset as an example to show the
effectiveness of the proposed cropping strategy. We compute the cropping rate
(CR) for each subject as the following formula:

$\mathrm{CR}=\frac{H\cdot W\cdot D-H^{{}^{\prime}}\cdot W^{{}^{\prime}}\cdot
D^{{}^{\prime}}}{H\cdot W\cdot D}$ (8)

where $(H,W,D)$ is the original image size and $(H^{{}^{\prime}}\cdot
W^{{}^{\prime}}\cdot D^{{}^{\prime}})$ is the image size after cropping.
Higher CR represents more background pixels being removed. The CR of the IXI
dataset is shown in Fig. 8. The average CR is 42.99%, and 80% of the samples
subtracted nearly 40% (38.49%) of the background points. Besides, two
representative data are shown in Fig. 9. The strategy of non-zero cropping is
widely employed in the analysis of 3D medical images. Yet, the presence of
noise in MRI scans undermines the effectiveness of this approach, leading to
suboptimal image cropping. Our proposed methodology addresses this issue by
meticulously eliminating the maximum amount of background pixels. This
approach gains importance in scenarios involving extensive 3D datasets,
ensuring a more precise and efficient analysis.

Our proposed method is not limited by the backbone structure and has
relatively high generalizability. More importantly, the proposed pretraining
strategy may reduce the performance gap among the various backbone models. It
is observed from Table 4 that the performance gap between Swin UNETR and UNet
is around 6.5 before pretraining, while it becomes approximately 2 after using
the proposed pretraining strategy. This indicates that the lightweight model
(e.g., vanilla UNet) can also exhibit competitive performance after employing
the proposed pretraining method, which holds significant implications for the
use of lightweight networks in clinical settings.

One limitation of this study is that our method may not effectively use the
other modal data that can not fully reflect the vessel morphology (e.g., T1).
The proposed pretraining method relies on the preprocessing workflow for
cerebrovascular extraction, while other modal data might not include any
vessel-related information. Any modality of data, such as computed tomographic
angiography, digital subtraction angiography, and optical coherence
tomography, capable of capturing vessel morphology, can derive advantages from
the proposed method. Besides, our pretraining method can serve as an auxiliary
training approach for a SSL model trained using other modal data or general
large vision models, and help to construct a domain-specific model for
cerebrovascular segmentation.

##  7 Conclusion

In this study, we introduce a Frangi filtering-based pretraining network to
effectively leverage the unlabeled TOF-MRA data. Our method capitalizes on the
prior knowledge of tubular structures by incorporating it into the pretraining
tasks, and shows great superiority over other SemiSL and SSL methods on four
cerebrovascular segmentation datasets. The ablation studies also show that our
proposed method significantly improves model performance, particularly with a
limited number of labeled data. Besides, our method is not restricted by model
architecture and can enhance the performance of various backbone structures.

Fig. 8: The cropping rate in IXI dataset. Fig. 9: Three representative
subjects for the cropping strategy. The number indicates the image size.

## Acknowledgments

This work was supported in part by the National Natural Science Foundation of
China under Grant: 62027901, 81930053, 81227901; Beijing Natural Science
Foundation: JQ22023; CAS Youth Innovation Promotion Association under Grant
Y2022055. The authors would like to acknowledge the instrumental and technical
support of Multimodal Biomedical Imaging Experimental Platform, Institute of
Automation, Chinese Academy of Sciences.

## References

  * Alderliesten et al. [2006] Alderliesten, T., Bosman, P.A., Niessen, W.J., 2006.  Towards a real-time minimally-invasive vascular intervention simulation system.  IEEE Transactions on Medical Imaging 26, 128–132.
  * Asadian et al. [2022] Asadian, A., Weidner, E., Jiang, L., 2022.  Self-supervised pretraining for differentially private learning.  arXiv preprint arXiv:2206.07125 .
  * Bai et al. [2023] Bai, Y., Chen, D., Li, Q., Shen, W., Wang, Y., 2023.  Bidirectional copy-paste for semi-supervised medical image segmentation, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 11514–11524.
  * Chen et al. [2022a] Chen, C., Zhou, K., Wang, Z., Xiao, R., 2022a.  Generative consistency for semi-supervised cerebrovascular segmentation from tof-mra.  IEEE Transactions on Medical Imaging 42, 346–353.
  * Chen et al. [2023] Chen, C., Zhou, K., Wang, Z., Zhang, Q., Xiao, R., 2023.  All answers are in the images: A review of deep learning for cerebrovascular segmentation.  Computerized Medical Imaging and Graphics , 102229.
  * Chen et al. [2019] Chen, K., Yao, L., Zhang, D., Chang, X., Long, G., Wang, S., 2019.  Distributionally robust semi-supervised learning for people-centric sensing, in: Proceedings of the AAAI Conference on Artificial Intelligence, pp. 3321–3328.
  * Chen et al. [2020] Chen, T., Kornblith, S., Norouzi, M., Hinton, G., 2020.  A simple framework for contrastive learning of visual representations, in: International conference on machine learning, PMLR. pp. 1597–1607.
  * Chen et al. [2022b] Chen, Y., Jin, D., Guo, B., Bai, X., 2022b.  Attention-assisted adversarial model for cerebrovascular segmentation in 3d tof-mra volumes.  IEEE Transactions on Medical Imaging 41, 3520–3532.
  * Cheplygina et al. [2019] Cheplygina, V., de Bruijne, M., Pluim, J.P., 2019.  Not-so-supervised: a survey of semi-supervised, multi-instance, and transfer learning in medical image analysis.  Medical image analysis 54, 280–296.
  * Di Noto et al. [2023] Di Noto, T., Marie, G., Tourbier, S., Alemán-Gómez, Y., Esteban, O., Saliou, G., Cuadra, M.B., Hagmann, P., Richiardi, J., 2023.  Towards automated brain aneurysm detection in tof-mra: Open data, weak labels, and anatomical knowledge.  Neuroinformatics 21, 21–34.
  * Duan et al. [2019] Duan, H.H., Su, G.Q., Huang, Y.C., Song, L.T., Nie, S.D., 2019.  Segmentation of pulmonary vascular tree by incorporating vessel enhancement filter and variational region-growing.  Journal of X-ray science and technology 27, 343–360.
  * Ericsson et al. [2022] Ericsson, L., Gouk, H., Loy, C.C., Hospedales, T.M., 2022.  Self-supervised representation learning: Introduction, advances, and challenges.  IEEE Signal Processing Magazine 39, 42–62.  doi:10.1109/MSP.2021.3134634.
  * Frangi et al. [1998] Frangi, A.F., Niessen, W.J., Vincken, K.L., Viergever, M.A., 1998.  Multiscale vessel enhancement filtering, in: Medical Image Computing and Computer-Assisted Intervention—MICCAI’98: First International Conference Cambridge, MA, USA, October 11–13, 1998 Proceedings 1, Springer. pp. 130–137.
  * Hassouna et al. [2006] Hassouna, M.S., Farag, A.A., Hushek, S., Moriarty, T., 2006.  Cerebrovascular segmentation from tof using stochastic models.  Medical image analysis 10, 2–18.
  * Hatamizadeh et al. [2021] Hatamizadeh, A., Nath, V., Tang, Y., Yang, D., Roth, H.R., Xu, D., 2021.  Swin unetr: Swin transformers for semantic segmentation of brain tumors in mri images, in: International MICCAI Brainlesion Workshop, Springer. pp. 272–284.
  * Hatamizadeh et al. [2022] Hatamizadeh, A., Tang, Y., Nath, V., Yang, D., Myronenko, A., Landman, B., Roth, H.R., Xu, D., 2022.  Unetr: Transformers for 3d medical image segmentation, in: Proceedings of the IEEE/CVF winter conference on applications of computer vision, pp. 574–584.
  * He et al. [2022] He, K., Chen, X., Xie, S., Li, Y., Dollár, P., Girshick, R., 2022.  Masked autoencoders are scalable vision learners, in: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 16000–16009.
  * He et al. [2023] He, Y., Yang, G., Ge, R., Chen, Y., Coatrieux, J.L., Wang, B., Li, S., 2023.  Geometric visual similarity learning in 3d medical image self-supervised pre-training, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 9538–9547.
  * Huang et al. [2023] Huang, S.C., Pareek, A., Jensen, M., Lungren, M.P., Yeung, S., Chaudhari, A.S., 2023.  Self-supervised learning for medical image classification: a systematic review and implementation guidelines.  NPJ Digital Medicine 6, 74.
  * Isensee et al. [2021] Isensee, F., Jaeger, P.F., Kohl, S.A., Petersen, J., Maier-Hein, K.H., 2021.  nnu-net: a self-configuring method for deep learning-based biomedical image segmentation.  Nature methods 18, 203–211.
  * Jiao et al. [2022] Jiao, R., Zhang, Y., Ding, L., Cai, R., Zhang, J., 2022.  Learning with limited annotations: a survey on deep semi-supervised learning for medical image segmentation.  arXiv preprint arXiv:2207.14191 .
  * Krishnan et al. [2022] Krishnan, R., Rajpurkar, P., Topol, E.J., 2022.  Self-supervised learning in medicine and healthcare.  Nature Biomedical Engineering 6, 1346–1352.
  * LaMontagne et al. [2019] LaMontagne, P.J., Benzinger, T.L., Morris, J.C., Keefe, S., Hornbeck, R., Xiong, C., Grant, E., Hassenstab, J., Moulder, K., Vlassenko, A.G., et al., 2019.  Oasis-3: longitudinal neuroimaging, clinical, and cognitive dataset for normal aging and alzheimer disease.  MedRxiv , 2019–12.
  * Liao et al. [2023] Liao, W., Li, X., Wang, Q., Xu, Y., Yin, Z., Xiong, H., 2023.  Cupre: Cross-domain unsupervised pre-training for few-shot cell segmentation.  arXiv preprint arXiv:2310.03981 .
  * Liao et al. [2022] Liao, W., Xiong, H., Wang, Q., Mo, Y., Li, X., Liu, Y., Chen, Z., Huang, S., Dou, D., 2022.  Muscle: Multi-task self-supervised continual learning to pre-train deep models for x-ray images of multiple body parts, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer. pp. 151–161.
  * Livne et al. [2019] Livne, M., Rieger, J., Aydin, O.U., Taha, A.A., Akay, E.M., Kossen, T., Sobesky, J., Kelleher, J.D., Hildebrand, K., Frey, D., et al., 2019.  A u-net deep learning framework for high performance vessel segmentation in patients with cerebrovascular disease.  Frontiers in neuroscience 13, 97.
  * Luo et al. [2021] Luo, X., Chen, J., Song, T., Wang, G., 2021.  Semi-supervised medical image segmentation through dual-task consistency, in: Proceedings of the AAAI Conference on Artificial Intelligence, pp. 8801–8809.
  * Ni et al. [2020] Ni, J., Wu, J., Wang, H., Tong, J., Chen, Z., Wong, K.K., Abbott, D., 2020.  Global channel attention networks for intracranial vessel segmentation.  Computers in biology and medicine 118, 103639.
  * Oktay et al. [2018] Oktay, O., Schlemper, J., Folgoc, L.L., Lee, M., Heinrich, M., Misawa, K., Mori, K., McDonagh, S., Hammerla, N.Y., Kainz, B., et al., 2018.  Attention u-net: Learning where to look for the pancreas.  arXiv preprint arXiv:1804.03999 .
  * Otsu [1979] Otsu, N., 1979.  A threshold selection method from gray-level histograms.  IEEE transactions on systems, man, and cybernetics 9, 62–66.
  * Özsarlak et al. [2004] Özsarlak, Ö., Van Goethem, J.W., Maes, M., Parizel, P.M., 2004.  Mr angiography of the intracranial vessels: technical aspects and clinical applications.  Neuroradiology 46, 955–972.
  * Qi and Luo [2020] Qi, G.J., Luo, J., 2020.  Small data challenges in big data era: A survey of recent progress on unsupervised and semi-supervised methods.  IEEE Transactions on Pattern Analysis and Machine Intelligence 44, 2168–2187.
  * Qi et al. [2023] Qi, Y., He, Y., Qi, X., Zhang, Y., Yang, G., 2023.  Dynamic snake convolution based on topological geometric constraints for tubular structure segmentation.  arXiv preprint arXiv:2307.08388 .
  * Ronneberger et al. [2015] Ronneberger, O., Fischer, P., Brox, T., 2015.  U-net: Convolutional networks for biomedical image segmentation, in: Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18, Springer. pp. 234–241.
  * Sanchesa et al. [2019] Sanchesa, P., Meyer, C., Vigon, V., Naegel, B., 2019.  Cerebrovascular network segmentation of mra images with deep learning, in: 2019 IEEE 16th international symposium on biomedical imaging (ISBI 2019), IEEE. pp. 768–771.
  * Shi et al. [2023] Shi, G., Yin, L., An, Y., Li, G., Zhang, L., Bian, Z., Chen, Z., Zhang, H., Hui, H., Tian, J., 2023.  Progressive pretraining network for 3d system matrix calibration in magnetic particle imaging.  IEEE Transactions on Medical Imaging 42, 3639–3650.  doi:10.1109/TMI.2023.3297173.
  * Shit et al. [2021] Shit, S., Paetzold, J.C., Sekuboyina, A., Ezhov, I., Unger, A., Zhylka, A., Pluim, J.P., Bauer, U., Menze, B.H., 2021.  cldice-a novel topology-preserving loss function for tubular structure segmentation, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 16560–16569.
  * Sichtermann et al. [2019] Sichtermann, T., Faron, A., Sijben, R., Teichert, N., Freiherr, J., Wiesmann, M., 2019.  Deep learning–based detection of intracranial aneurysms in 3d tof-mra.  American Journal of Neuroradiology 40, 25–32.
  * Tang et al. [2022] Tang, Y., Yang, D., Li, W., Roth, H.R., Landman, B., Xu, D., Nath, V., Hatamizadeh, A., 2022.  Self-supervised pre-training of swin transformers for 3d medical image analysis, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 20730–20740.
  * Wang and Lukasiewicz [2022] Wang, J., Lukasiewicz, T., 2022.  Rethinking bayesian deep learning methods for semi-supervised volumetric medical image segmentation, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 182–190.
  * Wang et al. [2023] Wang, X., Tang, F., Chen, H., Cheung, C.Y., Heng, P.A., 2023.  Deep semi-supervised multiple instance learning with self-correction for dme classification from oct images.  Medical Image Analysis 83, 102673.
  * Wu et al. [2022] Wu, Y., Wu, Z., Wu, Q., Ge, Z., Cai, J., 2022.  Exploring smoothness and class-separation for semi-supervised medical image segmentation, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer. pp. 34–43.
  * Xia et al. [2022] Xia, L., Zhang, H., Wu, Y., Song, R., Ma, Y., Mou, L., Liu, J., Xie, Y., Ma, M., Zhao, Y., 2022.  3d vessel-like structure segmentation in medical images by an edge-reinforced network.  Medical Image Analysis 82, 102581.
  * Xie et al. [2022] Xie, Y., Zhang, J., Xia, Y., Wu, Q., 2022.  Unimiss: Universal medical self-supervised learning via breaking dimensionality barrier, in: European Conference on Computer Vision, Springer. pp. 558–575.
  * Yang et al. [2022] Yang, X., Song, Z., King, I., Xu, Z., 2022.  A survey on deep semi-supervised learning.  IEEE Transactions on Knowledge and Data Engineering .
  * Yu et al. [2023] Yu, K., Sun, L., Chen, J., Reynolds, M., Chaudhary, T., Batmanghelich, K., 2023.  Drasclr: A self-supervised framework of learning disease-related and anatomy-specific representation for 3d lung ct images.  Medical Image Analysis , 103062URL: https://www.sciencedirect.com/science/article/pii/S1361841523003225, doi:https://doi.org/10.1016/j.media.2023.103062.
  * Zhang et al. [2023] Zhang, C., Zheng, H., Gu, Y., 2023.  Dive into the details of self-supervised learning for medical image analysis.  Medical Image Analysis , 102879.
  * Zhou et al. [2021a] Zhou, H.Y., Guo, J., Zhang, Y., Yu, L., Wang, L., Yu, Y., 2021a.  nnformer: Interleaved transformer for volumetric segmentation.  arXiv preprint arXiv:2109.03201 .
  * Zhou et al. [2021b] Zhou, H.Y., Lu, C., Yang, S., Han, X., Yu, Y., 2021b.  Preservational learning improves self-supervised medical image models by reconstructing diverse contexts, in: Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 3499–3509.

## Supplementary Material

###  7.1 Pretraining Datasets

The pretraining dataset consists of 5 publicly accessible datasets:

  * 1.

IXI111http://brain-development.org/ixi-dataset/: This dataset contains 570
subjects. 45 subjects with vessel masks Chen et al. [2022b] are used for
evaluation dataset. The rest 525 subjects are used for pretraining. The image
size of most data (497 subjects) is (512, 512, 100). The image size of 27
subjects is (1024, 1024, 92) and one subject is (1024, 1024, 91). The spacing
is (0.4688, 0.4688, 0.8000).

  * 2.

OASIS222https://www.oasis-brains.org: This dataset contains 525 subjects. The
spacing and image size are not identical among the subjects. The average image
size is (583, 765, 216). The spacing of most data (337 subjects) is (0.2995,
0.2995, 0.6000), and the average spacing is (0.2977, 0.2977, 0.5933).

  * 3.

TubeTK333https://public.kitware.com/Wiki/TubeTK/Data: This dataset contains
109 subjects. 42 subjects with vessel masks are used for the evaluation
dataset. The rest 67 subjects are used for pretraining. The image size is
(448, 448, 128) and the spacing is (0.5134, 0.5134, 0.8000).

  * 4.

ADAM444https://adam.isi.uu.nl/data/: This dataset contains 113 subjects. Each
subject contains one original MRA-TOF image and one bias field corrected MRA-
TOF image. Only the original data are used in this study. The average image
size is (556, 556, 131). The spacing is various and the average spacing is
(0.3524, 0.3524, 0.5447).

  * 5.

BrainAneurysm555https://openneuro.org/datasets/ds003949/versions/1.0.1: This
dataset contains 284 subjects, of which 127 are healthy controls and 157 are
patients with brain aneurysms. Four subjects (sub-200, sub-235, sub-315 and
sub-450) are removed due to the incomplete image data. The average spacing is
(0.4021, 0.4021, 0.6613), and the average image size is (466, 546, 147).

###  7.2 Evaluation Datasets

The raw image data of the three labeled datasets are accessible at the
following public websites:

  * 1.

TubeTK-42: https://public.kitware.com/Wiki/TubeTK/Data

  * 2.

Brains: https://www.nitrc.org/projects/icbmmra/

  * 3.

EDEN: https://zenodo.org/record/3994749

We provide the voxel-wise vessel mask for the three datasets in this study.

###  7.3 The implementations of baseline models

For the baseline models, we use the following publicly available
implementations.

  * 1.

GCS: https://github.com/MontaEllis/SSL-For-Medical-Segmentation

  * 2.

DTC: https://github.com/HiLab-git/DTC

  * 3.

GBDL: https://github.com/Jianf-Wang/GBDL

  * 4.

SSNet: https://github.com/ycwu1997/SS-Net

  * 5.

BCP: https://github.com/DeepMed-Lab-ECNU/BCP

  * 6.

MAE: https://github.com/facebookresearch/mae

  * 7.

Swin UNETR: https://github.com/Project-MONAI/research-
contributions/tree/main/SwinUNETRR

  * 8.

UniMiSS: https://github.com/YtongXie/UniMiSS-code

  * 9.

GVSL: https://github.com/YutingHe-list/GVSL

  * 10.

PCRL: https://github.com/Luchixiang/PCRL

The UNet is implemented in MONAI666https://monai.io/.