1

A Deep Structured Learning Approach Towards
Automating Connectome Reconstruction from
3D Electron Micrographs

arXiv:1709.02974v3 [cs.CV] 24 Sep 2017

Jan Funke?, Fabian David Tschopp?, William Grisaitis, Arlo Sheridan,
Chandan Singh, Stephan Saalfeld, Srinivas C. Turaga

F

Abstract—We present a deep structured learning method for neuron
segmentation from 3D electron microscopy (EM) which improves significantly upon the state of the art in terms of accuracy and scalability. Our
method consists of a 3D U-N ET architecture, trained to predict affinity
graphs on voxels, followed by a simple and efficient iterative region
agglomeration. We train the U-N ET using a new structured loss function
based on M ALIS that encourages topological correctness. Our M ALIS extension consists of two parts: First, we present an O(n log(n)) method to
compute the loss gradient, which improves over the originally proposed
O(n2 ) algorithm. Second, we compute the gradient in two separate
passes to avoid spurious gradient contributions in early training stages.
Our affinity predictions are accurate enough that simple learning-free
percentile-based agglomeration outperforms more involved methods
used earlier on inferior predictions. We present results on three EM
datasets (C REMI, F IB -25, and S EG E M) of different imaging techniques
and animals where we achieve relative improvements over previous
results of 27%, 15%, and 250%, respectively. Our findings suggest
that a single 3D segmentation strategy can be applied to both nearly
isotropic block-face EM data and anisotropic serial sectioned EM data.
The runtime of our method scales with O(n) in the size of the volume and
achieves a throughput of about 2.6 seconds per megavoxel, qualifying
our method for the processing of very large datasets.

1

I NTRODUCTION

Precise reconstruction of neural connectivity is of great
importance to understand the function of biological nervous
systems. 3D electron microscopy (EM) is the only available
imaging method with the resolution necessary to visualize
and reconstruct dense neural morphology without ambiguity. At this resolution, however, even moderately small
neural circuits yield image volumes that are too large for
manual reconstruction. Therefore, automated methods for
neuron tracing are needed to aid human analysis.
We present a structured deep learning based image segmentation method for reconstructing neurons from 3D electron microscopy which improves significantly upon state of
the art in terms of accuracy and scalability. For an overview,
see Fig. 1, top row. The main components of our method
are: (1) Prediction of 3D affinity graphs using a 3D U-N ET
architecture [1], (2) a structured loss based on M ALIS [2] to
train the U-N ET to minimize topological errors, and (3) an
?

these authors contributed equally

efficient O(n) agglomeration scheme based on quantiles of
predicted affinities.
The choice of using a 3D U-N ET architecture to predict
voxel affinities is motivated by two considerations: First,
U-N ETs have already shown superior performance on the
segmentation of 2D [3] and 3D [1] biomedical image data.
One of their favourable properties is the multi-scale architecture which enables computational and statistical efficiency.
Second, U-N ETs efficiently predict large regions. This is
of particular interest in combination with training on the
M ALIS structured loss, for which we need affinity predictions in a region.
We train our 3D U-N ET to predict affinities using an
extension of the M ALIS loss function [2]. Like the original
M ALIS loss, we minimize a topological error on hypothetical thresholding and connected component analysis on the
predicted affinities. We extended the original formulation to
derive the gradient with respect to all predicted affinities (as
opposed to sparsely sampling them), leading to denser and
faster gradient computation. Furthermore, we compute the
M ALIS loss in two passes: In the positive pass, we constrain
all predicted affinities between and outside of ground-truth
regions to be 0, and in the negative pass, we constrain affinities inside regions to be 1 which avoids spurious gradients
in early training stages.
Although the training is performed assuming subsequent thresholding, we found iterative agglomeration of
fragments (or “supervoxels”) to be more robust to small
errors in the affinity predictions. To this end, we extract
fragments running a watershed algorithm on the predicted
affinities. The fragments are then represented in a region
adjacency graph (RAG), where edges are scored to reflect
the predicted affinities between adjacent fragments: edges
with small scores will be merged before edges with high
scores. We discretize edge scores into k evenly distributed
bins, which allows us to use a bucket priority queue for
sorting. This way, the agglomeration can be carried out with
a worst-case linear runtime.
The resulting method (prediction of affinities, watershed,
and agglomeration) scales favourably with O(n) in the size
n of the volume, a crucial property for neuron segmentation
from EM volumes, where volumes easily reach several

2

3D U-N ET
+ M ALIS

A
(a) 3D U-N ET.

A
B
A
C
50nm

y

final segmentation

oversegmentation

affinities

D
C

B

D

A

C

segmentation

z
x

(b) Inter-voxel affinities.

b

A

C

[b,e]

C

b,e → b

threshold

raw volume

percentile agglomeration

seeded
watershed

b

A

C

c

D

[d,c,b,e]
b,d → b

e
A

a

B

b

C

c

D

e
d
RAG and merge hierarchy

[a,d,c,b,e]
queue

(c) Percentile agglomeration.

Figure 1: Overview of our method (top row). Using a 3D U-N ET (a), trained with the proposed constrained M ALIS loss,
we directly predict inter-voxel affinities from volumes of raw data. Affinities provide advantages especially in the case of
low-resolution data (b). In the example shown here, the voxels cannot be labeled correctly as foreground/background: If A
were labeled as foreground, it would necessarily merge with the regions in the previous and next section. If it were labeled
as background, it would introduce a split. The labeling of affinities on edges allows B and C to separate A from adjacent
sections, while maintaining connectivity inside the region. From the predicted affinities, we obtain an over-segmentation
that is then merged into the final segmentation using a percentile-based agglomeration algorithm (c).
hundreds of terabytes. This is a major advantage over
current state-of-the-art methods that all follow a similar
pattern. First, voxel-wise predictions are made using a
deep neural network. Subsequently, fragments are obtained
from these predictions which are then merged using either
greedy (C ELIS [4], G ALA [5]) or globally optimal objectives
(M ULTI C UT [6] and lifted M ULTI C UT [7], [8]). Current efforts focus mostly on the merging of fragments: Both C ELIS
and G ALA train a classifier to predict scores for hierarchical
agglomeration which increases the computational complexity of agglomeration during inference. Similarly, the M UL TI C UT variants train a classifier to predict the connectivity
of fragments that are then clustered by solving a computationally expensive combinatorial optimization problem.
Our proposed fragment agglomeration method drastically
reduces the computation complexity compared to previous
merge methods and does not require a separate training
step.
We demonstrate the efficacy of our method on three diverse datasets of EM volumes, imaged by three different 3D
electron microscopy techniques: C REMI (ssTEM, Drosophila),
F IB -25 (FIBSEM, Drosophila), and S EG E M (SBEM, mouse
cortex). Our method significantly improves over the current

state of the art in each of these datasets, outperforming in
particular computationally more expensive methods without favorable worst-case runtime guarantees.
We made the source code for training1 and agglomeration2 publicly available, together with usage example scripts
to reproduce our C REMI results3 .

2

M ETHOD

2.1 Deep multi-scale convolutional network for predicting 3D voxel affinities
We use a 3D U-N ET architecture [1] to predict voxel affinities on 3D volumes. We use the same architecture for all
investigated datasets which we illustrate in Fig. 1a. In particular, our 3D U-N ET consists of four levels of different
resolutions. In each level, we perform at least one convolution pass (shown as blue arrows in Fig. 1a) consisting of
two convolutions (kernel size 3×3×3) followed by rectified
linear units. Between the levels, we perform max pooling
on variable kernel sizes depending on the dataset resolution
1. https://github.com/naibaf7/caffe
2. https://github.com/funkey/waterz
3. http://cremi.org/static/data/20170312 mala v2.tar.gz

3

for the downsampling pass (yellow arrows), as well as
transposed convolution of the same size for upsampling
(brown arrows). The results of the upsampling pass are
further concatenated with copies of the feature maps of
the same level in the downsampling pass (red arrows),
cropped to account for context loss in the lower levels.
Details of the individual passes are shown in Fig. 6. A more
detailed description of the U-N ET architectures for each of
the investigated datasets can be found in Fig. 5.
We chose to predict voxel affinities on edges between
voxels instead of labeling voxels as foreground/background
to allow our method to handle low spatial resolutions. As
we illustrate in Fig. 1b, a low z resolution (common for serial
section EM) renders a foreground/background labeling of
voxels impossible. Affinities, on the other hand, effectively
increase the expressiveness of our model and allow to
obtain a correct segmentation. Furthermore, affinities easily
generalize to arbitrary neighborhoods and might thus allow
the prediction of longer range connectivity.

2.2

Training using constrained M ALIS

We train our network using an extension of the M ALIS
loss [2]. This loss, that we term constrained M ALIS, is designed to minimize topological errors in a segmentation
obtained by thresholding and connected component analysis. Although thresholding alone will unlikely produce
accurate results, it serves as a valuable proxy for training:
If the loss can be minimized for thresholding, it will in
particular be minimized for agglomeration. To this end,
in each training iteration, a complete affinity prediction of
a 3D region is considered. Between every pair of voxels,
we determine the maximin affinity edge, i.e., the highest
minimal edge over all paths connecting the pair. This edge
is crucial as it determines the threshold under which the
two voxels in question will be merged. Naturally, for voxels
that are supposed to belong to the same region, we want
the maximin edge affinity to be as high as possible, and for
voxels of different regions as low as possible.
Our extension consists of two parts: First, we improve
the computational complexity of the M ALIS loss by presenting an O(n log(n)) method for the computation of the
gradient (thus improving over previous O(n2 )). Second, we
compute the gradient in two separate passes, once for affinities inside ground-truth objects (positive pass), and once for
affinities between and outside of ground-truth objects.
A maximin edge between two voxels u and v is an edge
mm(u, v) with lowest affinity on the overall highest affinity
∗ connecting u and v, i.e.,
path Pu,v
∗
Pu,v
= arg max min Ai, j
P ∈ Pu, v

(i, j)∈P

mm(u, v) = arg min Ai, j , (1)
∗
(i, j)∈Pu,
v

where Pu,v denotes the set of all paths between u and v, and
Ae denotes the predicted affinity of edge e. If we imagine a
simple thresholding on the affinity graph, such that edges
with affinities below a threshold θ are removed, then the
affinity of the maximin edge mm(u, v) is equal to the highest
threshold under which nodes u and v would still be part
of the same connected component. Taking advantage of the
importance of maximin edges, the M ALIS loss favors high

maximin affinities between nodes with the same label, and
low otherwise:
Õ
2
l(A) =
δid(u) id(v) − Amm(u,v) ,
(2)
u<v

where δ is the Kronecker delta and id(u) is the ground-truth
label of voxel u. It can easily be seen that maximin edges are
shared between voxel pairs. In fact, the union of all maximin
edges forms a maximal spanning tree (MST),
Ø
mm(u, v) = MST(A).
(3)
u,v

Consequently, we are able to identify the maximin edge and
compute its loss for each voxel pair in O(n log(n)) time. We
thus improve over a previous method [2] that required O(n2 )
and thus had to fall back to sparse sampling of voxel pairs.
Note that this only affects the training of the network, the
affinity prediction during test time scales linearly with the
volume size.
We further extend previous work by computing the
maximin edge losses in two passes: First for edges within
the same region (positive pass), second for edges between
different regions (negative pass). As shown in Fig. 2, in the
positive pass, we assume that all edges between regions
have been predicted correctly and set their affinities to zero.
Consequently, only maximin edges inside a region are found
and contribute to the loss. This obviates an inefficiency in
a previous formulation [2], where a spurious high-affinity
(i.e., false positive) path leaving and entering a region might
connect two voxels inside the same region. In this case, the
maximin edge could lie outside of the considered region,
resulting in an unwanted gradient contribution that would
reinforce the false positive. Analogously, in the negative
pass, all affinities inside the same region are set to one
to avoid reinforcement of false negatives inside regions.
Finally, the gradient contributions of both passes are added
together.
2.3

Hierarchical agglomeration

Our method for hierarchical agglomeration of segments
from the predicted affinities consists of two steps. First, we
use a heuristic to extract small fragments directly from the
predicted affinities. Second, we iteratively score and merge
adjacent fragments into larger objects until a predefined
threshold is reached.
2.3.1

Fragment extraction

The extraction of fragments is a crucial step for the subsequent agglomeration. Too many fragments slow down
the agglomeration unnecessarily and increase its memory
footprint. Too few fragments, on the other hand, are subject
to undersegmentation that cannot be corrected.
Empirically, we found a seeded watershed to deliver
the best trade-off between fragment size and segmentation
accuracy across all investigated datasets. For the seeded
watershed, we first average the predicted affinities for each
voxel to obtain a volume of boundary predictions. We subsequently threshold the boundary predictions at 0.5 and perform a distance transform on the resulting mask. Every local
maximum is taken as a seed, from which we grow basins

4

1
1
1
+
A
A
0
0.16599284152779
0.16033112623744
0.75213431443653
0.16599284152779
0.75213431443653
0.16033112623744
0.24786568556347
1
1
0
0.73812714085827
0.6049847870157
0.1163344565855
0.6908791841431 0.73812714085827
0.6049847870157
0.6908791841431
0.1163344565855
0.26187285914173
0.3950152129843
0.3091208158569
1
1
1
0
0
0.68449733643071
0.17299854181381
0.21019124170309
0.68449733643071
0.17299854181381
0.21019124170309 0.31550266356929
1
0
0
0
0
0.14125302771165
0.25307619644007
0.21293037781163
0.7037308046146
0.14125302771165
0.25307619644007
0.21293037781163
0.7037308046146 -0.25307619644007
-0.21293037781163
-0.7037308046146
C
C
0
0
0.65002569018399
0.16357205261643
0.74865410139256
0.74865410139256
0.65002569018399
0.16357205261643
-0.65002569018399
0.25134589860744
1
0.79635202996263
0
0.79635202996263
0
0.72731170483274
0.62714567856265
0.71476834561432
0.71476834561432
0.72731170483274
0.62714567856265
0.20364797003737
0.28523165438568
1
1 -0.72731170483274
0.27272447616454
0.62134905160468
0.64376995863568
0.27272447616454
0.62134905160468
0.64376995863568
0.72727552383546
0.37865094839532
1
1
1
B
B
0
0
0
-

(a) Predicted affinities.

(b) Positive pass.

(c) Negative pass.

(d) Gradient of loss.

Figure 2: Illustration of the constrained M ALIS loss. Given predicted affinities (blue low, red high) and a ground-truth
segmentation (a), losses on maximin edges are computed in two passes: In the positive pass, (b), affinities of edges between
ground-truth regions are set to zero (blue), in the negative pass (c), affinities within ground-truth regions are set to one
(red). In either case, a maximal spanning tree (shown as shadow) is constructed to identify maximin edges. Note that,
in this example, edge A is not a maximin edge in the positive pass since the incident voxels are already connected by a
high affinity path. In contrast, edge B is the maximin edge of the bottom left voxel to any other voxel in the same region
and thus contributes to the loss. Similarly, C is the maximin edge connecting voxels of different ground-truth regions and
contributes during the negative pass to the loss. The resulting gradients of the loss for each edge affinity is shown in (d)
(positive values in red, negative in blue).

averaged affinities

threshold at 0.5

distance transform

seeded watershed

Figure 3: Illustration of the seeded watershed heuristic.
using a standard watershed algorithm [9] on the boundary
predictions. For an example, see Fig. 3. As argued above,
voxel-wise predictions are not fit for anisotropic volumes
with low z-resolution (see Fig. 1b). To not re-introduce a flaw
that we aimed to avoid by predicting affinities instead of
voxel-wise labels in the first place, we perform the extraction
of fragments xy-section-wise for anisotropic volumes.
2.3.2 Fragment agglomeration
For the agglomeration, we consider the region adjacency
graph (RAG) of the extracted fragments. The RAG is an
annotated graph G = (V, E, f ), with V the set of fragments,
E ⊆ V × V edges between adjacent fragments, and f : E 7→ R
an edge scoring function. The edge scoring function is

designed to prioritize merge operations in the RAG, i.e., the
contraction of two adjacent nodes into one, such that edges
with lower scores are merged earlier. Given an annotated
RAG, a segmentation can be obtained by finding the edge
with the lowest score, merge it, recompute the scores of
edges affected by the merge, and iterate until the score
of the lowest edge hits a predefined threshold θ. In the
following, we will denote by Gi the RAG after i iterations
(and analogously by Vi , Ei , and fi its nodes, edges, and
scores), with G0 = G as introduced above. We will ”reuse”
nodes and edges, meaning Vi+1 ⊂ Vi and Ei+1 ⊂ Ei .
Given that the initial fragments are indeed an oversegmentation, it is up to the design of the scoring function
and the threshold θ to ensure a correct segmentation. The
design of the scoring function can be broken down into the
initialization of f0 (e) for e ∈ E0 (i.e., the initial scores) and the
update of fi (e) for e ∈ Ei ; i > 0 after a merge of two regions
a, b ∈ Vi−1 . For the update, three cases can be distinguished
(for an illustration see Fig. 4): (1) e was not affected by the
merge, (2) e is incident to a or b but represents the same
contact area between two regions as before, and (3) e results
from merging two edges of Ei−1 into one (the other edge will
be deleted). In the first two cases, the score does not change,
i.e., fi (e) = fi−1 (e), since the contact area between the nodes
linked by e remains the same. In the latter case, the contact
area is the union of the contact area of the merged edges, and
the score needs to be updated accordingly. Acknowledging
the merge hierarchy of edges (as opposed to nodes), we will
refer to the leaves under a merged edge e as initial edges,
denoted by E ∗ (e) ⊆ E0 .
In our experiments, we initialize the edge scores f (e) for
e ∈ E0 with one minus the maximum affinity between the
fragments linked by e and update them using a quantile
value of scores of the initial edges under e. This strategy has
been found empirically over a range of possible implementations of f (see Section 3).
Implemented naively, hierarchical agglomeration has a
worst-case runtime complexity of at least O(n log(n)), where

5

method
U-N ET M ALA
U-N ET
FlyEM [10]
CELIS [4]
CELIS+MC [4]

E
e
d
D

c

a

A

C
B

VOI split
0.891
1.205
1.490
1.426
1.037

VOI merge
0.180
0.316
0.462
0.208
0.229

VOI sum
1.071
1.520
1.952
1.634
1.266

(a) Results on F IB -25, evaluated on whole test volume.

b
method
U-N ET M ALA
U-N ET
FlyEM [10]
CELIS [4]
CELIS+MC [4]

ge
er

m

E
a

A

→
,B
A

e

D

c

A

b

C

VOI merge
0.198
0.471
0.251
0.166
0.216

VOI sum
2.151
2.914
3.411
3.568
2.570

(b) Results on F IB -25, evaluated on synaptic sites.
method
U-N ET M ALA
U-N ET
LMC [8]
CRunet
LFC

Figure 4: Illustration of the three different edge update cases
during a merge: (1) a is not affected by the merge operation,
(2) b and e link one of the merge nodes to a neighbor, which
is not adjacent to the other merge node, and (c) c and d link
one of the merge nodes to a node D, which is adjacent to
both A and B. c and d thus get merged into one edge and its
score needs to be updated.

n = |E0 | is the number of edges in the initial RAG. This is due
to the requirement of finding, in each iteration, the cheapest
edge to merge, which implies sorting of edges based on
their scores. Furthermore, the edge scoring function has to
be evaluated O(n) times, once for each affected edge of a
node merge (assuming nodes have a degree bounded by a
constant). For the merge function suggested above, a quantile of O(n) initial edge scores has to be found in the worst
case, resulting in a total worst-case runtime complexity of
O(n log(n) + n2 ).
To avoid this prohibitively high runtime complexity, we
propose to discretize the initial scores f0 into k bins, evenly
spaced in the interval [0, 1]. This simple modification has
two important consequences: First, a bucket priority queue
for sorting edge scores can be used, providing constant
time insert and pop operations. Second, the computation
of quantiles can be implemented in constant time and space
by using histograms of the k possible values. This way, we
obtain constant-time merge iterations (pop an edge, merge
nodes, update scores of affected edges), applied at most n
times, thus resulting in an overall worst-case complexity of
O(n). With k = 256 bins, we noticed no sacrifice of accuracy
in comparison to the non-discretized variant.
The analysis above holds only if we can ensure that the
update of the score of an edge e, and thus the update of
the priority queue, can be performed in constant time. In
particular, it is to be avoided to search for e in its respective
bucket. We note that for the quantile scoring function (and
many more), the new edge score fi (e) after merging an edge
f ∈ Ei−1 into e ∈ Ei−1 is always greater than or equal to its
previous score. We can therefore mark e as stale and f as
deleted and proceed merging without resorting the queue
or altering the graph. Whenever a stale edge is popped
from the priority queue, we compute its actual score and
insert it again into the queue. Not only does this ensure
constant time updates of edge scores and the priority queue,

VOI split
1.953
2.442
3.160
3.401
2.354

VOI split
0.425
0.979
0.597
1.081
1.085

VOI merge
0.181
0.546
0.272
0.389
0.140

VOI sum
0.606
1.524
0.868
1.470
1.225

CREMI score
0.289
0.793
0.398
0.566
0.616

(c) Results on C REMI.
method
U-N ET M ALA
U-N ET
SegEM [11]

IED split
6.259
6.903
2.121

IED merge
21.337
1.719
3.951

IED total
4.839
1.377
1.380

(d) Results on S EG E M.

Table 1: Qualitative results of our method (U-N ET M ALA)
compared to the respective state of the art on the testing volumes of each dataset and a baseline (U-N ET). Highlighted
in bold are the names of our method and the best value in
each column. Measures shown are variation of information
(VOI, lower is better), CREMI score (geometric mean of VOI
and adapted RAND error, lower is better), and inter-error
distance in µm (IED, higher is better) evaluated on traced
skeletons of the test volume. The IED has been computed
using the TED metric on skeletons [12] with a distance
threshold of 52 nm (corresponding to the thickness of two
z-sections). C REMI results are reported as average over all
testing samples, individual results can be found in Fig. 8.
it also avoids computing scores for edges that are never used
for a merge. This can happen if the threshold is hit before
considering the edge, or if the edge got marked as deleted
as a consequence of a nearby merge.

3

R ESULTS

Datasets We present results on three different and diverse
datasets: C REMI4 , F IB -25 [10], and S EG E M [11] (see Table 2
for an overview). These datasets sum up to almost 15
gigavoxels of testing data, with F IB -25 alone contributing
13.8 gigavoxels, thus challenging automatic segmentation
methods for their efficiency. In fact, only two methods have
so far been evaluated on F IB -25 [4], [10]. Another challenge
is posed by the C REMI dataset: Coming from serial section
EM, this dataset is highly anisotropic and contains artifacts
like support film folds, missing sections, and staining precipitations. Regardless of the differences in isotropy and
4. https://cremi.org

6

24

12

(84, 268, 268)

3
(56, 56, 56)

60

120

60

12

(60, 20, 20)

(80, 88, 88)

300

600

300

60

(64, 8, 8)

(76, 28, 28)

(72, 8, 8)

300

(68, 4, 4)

1500

Figure 5: Overview of the U-net architecture used for the C REMI dataset. The architectures for F IB -25 and S EG E M are
similar, with changes in the input and output sizes (in: (132, 132, 132), out: (44, 44, 44) for F IB -25 and in: (188, 188, 144), out:
(100, 100, 96) for S EG E M) and number of feature maps for F IB -25 (24 in the first layer, increased by a factor of 3 for lower
layers).
Name
C REMI

Imaging
ssTEM

Tissue
Drosophila

Resolution
4×4×40 nm

F IB -25
S EG E M

FIBSEM
SBEM

Drosophila
mouse cortex

8×8×8 nm
11×11×26 nm

Training Data
3 volumes of
1250×1250×125 voxels
520×520×520 voxels
279 volumes of
100×100×100 voxels

Testing Data
3 volumes of
1250×1250×125 voxels
13.8 gigavoxels
400×400×350 voxels
(skeletons)

Table 2: Overview of used datasets.
presence of artifacts, we use the same method (3D U-N ET
training, prediction, and agglomeration) for all datasets.
The size of the receptive field of the U-N ET was set for
each dataset to be approximately one µm3 , i.e., 213×213×29
for C REMI, 89×89×89 for F IB -25, and 89×89×49 for S EG E M.
For the C REMI dataset, we also pre-aligned training and
testing data with an elastic alignment method [13], using
the padded volumes provided by the challenge.
Training We implemented and trained our network using
the CAFFE library on modestly augmented training data for
which we performed random rotations, transpositions and
flips, as well as elastic deformations. On the anisotropic
C REMI dataset, we further simulated missing sections by
setting intensity values to 0 (p = 0.05) and low contrast sections by multiplying the intensity variance by 0.5 (p = 0.05).
We used the Adam optimizer [14] with an initial learning
rate of α = 10−4 , β1 = 0.95, β2 = 0.99, and  = 10−8 .
Quantitative results On each of the investigated datasets,
we see a clear improvement in accuracy using our method,
compared to the current state of the art. We provide quantitative results for each of the datasets individually, where
we compare our method (labeled U-N ET M ALA) against
different other methods. We also include a baseline (labeled
U-N ET) in our analysis, which is our method, but trained
without the constrained M ALIS loss. In Table 1, we report
the segmentation obtained on the best threshold found in
the respective training datasets. In Fig. 7, we show the
split/merge curve for varying thresholds of our agglomeration scheme.
For S EG E M, we do not use the metric proposed by
Berning et al. [11], as we found it to be problematic: The
authors suggest an overlap threshold of 2 to compensate

for inaccuracies in the ground-truth, however this has the
unintended consequence of ignoring some neurons in the
ground-truth for poor segmentations. For the S EG E M segmentation (kindly provided by the authors), 195 out of 225
ground-truth skeletons are ignored because of insufficient
overlap with any segmentation label. On our segmentation,
only 70 skeletons would be ignored, thus the results are
not directly comparable. Therefore, we performed a new
IED evaluation using TED [12], a metric that allows slight
displacement of skeleton nodes (we chose 52 nm in this case)
in an attempt to minimize splits and merges. This metric
reveals that our segmentations (U-N ET M ALA) improve
over both split and merge errors, over all thresholds of
agglomeration, including the initial fragments (see Fig. 7c).
Qualitative results Renderings of 11 and 23 randomly
selected neurons, reconstructed using the proposed method,
are shown for the test regions of C REMI and F IB -25 in Fig. 9
and Fig. 10, respectively.
Dataset (an)isotropy Save for minor changes in the network
architectures and the generation of initial fragments, our
method works unchanged on both near-isotropic block-face
datasets (F IB -25, S EG E M) as well as on highly anisotropic
serial-section datasets (C REMI). These findings suggest that
there is no need for specialized constructions like dedicated features for anisotropic volumes or separate classifiers
trained for merging of fragments within or across sections.
Merge functions Our method for efficient agglomeration
allows using a range of different merge functions. In Table 3,
we show results for different choices of quantile merge
functions, mean affinity, and an agglomeration baseline
proposed in [15] on datasets C REMI and F IB -25. Even across
these very different datasets, we see best results for affinity

7

∗

conv( f 1, f 2 )

∗

(3, 3, 3)

w

(3, 3, 3)

w−2

w−4

max

kx

pool(k x, k y, k z )

w/k x

w

up(k x, k y, k z , f )

⊗

w

∗

+

(k x , k y, k z )

(1, 1, 1)

w·k

w·k

crop(cx, cy, cz )
w − cx

w

Figure 6: Details of the convolution (blue), max-pooling (yellow), upsampling (brown), and copy-crop operations (red). “∗”
denotes a convolution, “ ” a rectified linear unit, and “⊗” the Kronecker matrix product.

15%
25%
50%
75%
85%
mean
Zlateski [15]

VOI split
0.583
0.441
0.397
0.347
0.347
0.380
1.015

VOI merge
0.063
0.092
0.056
0.074
0.084
0.058
1.010

VOI sum
0.646
0.533
0.453
0.421
0.431
0.438
2.025

CREMI score
0.212
0.188
0.156
0.146
0.156
0.149
0.364

(a) C REMI (training data).

15%
25%
50%
75%
85%
mean
Zlateski [15]

VOI split
1.480
1.393
1.115
1.085
1.176
1.221
1.054

VOI merge
0.364
0.163
0.234
0.318
0.394
0.198
1.017

VOI sum
1.844
1.555
1.350
1.402
1.570
1.418
2.071

(b) F IB -25.

Table 3: Results for different merge functions of our method
compared with the agglomeration strategy proposed in [15].
We show the results at the threshold achieving the best
score in the respective dataset (CREMI score for C REMI,
VOI for F IB -25). Note that, for this analysis, we used the
available training datasets which explains deviations from
the numbers shown in Table 1.

dataset
C REMI
F IB -25
S EG E M

U-N ET
3.04
0.66
2.19

watershed
0.23
0.92
0.25

agglomeration
0.83
1.28
0.14

total
4.10
2.86
2.58

Table 4: Throughput of our method for each of the investigated datasets in seconds per megavoxel.
quantiles between 50% and 75%.
Throughput Table 4 shows the throughput of our method
for each dataset, broken down into affinity prediction (UN ET), fragment extraction (watershed), and fragment agglomeration (agglomeration). For C REMI and S EG E M, most
time is spent on the prediction of affinities. The faster
predictions in F IB -25 are due to less feature maps used in
the network for this dataset.
To empirically confirm the theoretical speedup of using
a bucket queue for agglomeration, we show in Fig. 7d a
speed comparison of the proposed linear-time agglomeration against a naive agglomeration scheme for volumes of
different sizes.

4

D ISCUSSION

A remarkable property of our method is that it requires
almost no tuning to operate on datasets of different characteristics, except for minor changes in the size of the receptive

8
2

2
U-N ET M ALA

U-N ET M ALA
U-N ET

U-N ET
1.5

1.5

LMC [8]

FlyEM [10]
CELIS [4]

VOI split

VOI split

CRunet
LFC

1

0.5

0

CELIS+MC [4]

1

0.5

0

0.2

0.6

0.4

0.8

0

1

0

0.5

1

2

1.5

2.5

VOI merge

VOI merge

(a) VOI split/merge on C REMI.

(b) VOI split/merge on F IB -25.

3

bucket queue O(n)

U-N ET

priority queue O(n log(n))

SegEM [11]
100.5

2
time in s

distance between split µ m

·104
U-N ET MALA

1

100
0
101

102

0

0.2

0.4

0.6

0.8

1

distance between merge µ m

size in voxels

(c) IED split/merge on S EG E M.

(d) Runtime analysis.

1.2

1.4
·109

Figure 7: (a-c) Split merge curves of our method (lines) for different thresholds on the C REMI, F IB -25, and S EG E M datasets,
compared against the best-ranking competing methods (dots). (d) Performance comparison of a naive agglomeration
scheme (priority queue, O(n log(n))) versus our linear-time agglomeration (bucket queue, O(n)).
field of the U-N ET and initial fragment generation. This suggests that there is no need for the development of dedicated
algorithms for different EM modalities. Across all datasets,
our results indicate that affinity predictions on voxels are
sufficiently accurate to render sophisticated post-precessing
obsolete. It remains an open question whether fundamentally different approaches, like the recently reported floodfilling network [16], also generalize in a similar way. At
the time of writing, neither code nor data were publicly
available for a direct comparison.
Furthermore, the U-N ET is the only part in our method
that requires training, so that all training data can (and
should) be used to correctly predict affinities. This is an
advantage over current state-of-the-art methods that require careful splitting of precious training data into nonoverlapping sets used to train voxel-wise predictions and
an agglomeration classifier (or accepting the disadvantages
of having the sets overlap).
Although linear in runtime and memory, parallelization
of our method is not trivial due to the hierarchical agglomeration. However, we can take advantage of the fact that scores
of merged edges only increase (as discussed in Section 2.3.2):
This allows to split a given volume into regions, where
within a region agglomeration can be performed independently, as long as the score does not exceed any boundary
score to a neighboring region.

R EFERENCES
[1]

Çiçek, Ö., Abdulkadir, A., Lienkamp, S. S., Brox, T. & Ronneberger,

O. 3D U-Net: learning dense volumetric segmentation from sparse
annotation. In MICCAI, 424–432 (Springer, 2016).
[2] Briggman, K., Denk, W., Seung, S., Helmstaedter, M. N. & Turaga,
S. C. Maximin affinity learning of image segmentation. In Bengio,
Y., Schuurmans, D., Lafferty, J. D., Williams, C. K. I. & Culotta, A.
(eds.) Advances in Neural Information Processing Systems 22, 1865–
1873 (Curran Associates, Inc., 2009).
[3] Ronneberger, O., Fischer, P. & Brox, T. U-Net: Convolutional
networks for biomedical image segmentation. In MICCAI, 234–
241 (Springer, 2015).
[4] Maitin-Shepard, J. B., Jain, V., Januszewski, M., Li, P. & Abbeel,
P. Combinatorial energy learning for image segmentation. In
Advances in Neural Information Processing Systems 29, 1966–1974
(Curran Associates, Inc., 2016).
[5] Nunez-Iglesias, J., Kennedy, R., Plaza, S. M., Chakraborty, A. &
Katz, W. T. Graph-based active learning of agglomeration (gala):
a python library to segment 2d and 3d neuroimages. Frontiers in
neuroinformatics 8, 34 (2014).
[6] Andres, B. et al. Globally optimal closed-surface segmentation for
connectomics. In European Conference on Computer Vision, 778–791
(Springer, 2012).
[7] Keuper, M. et al. Efficient decomposition of image and mesh
graphs by lifted multicuts. In The IEEE International Conference
on Computer Vision (ICCV) (2015).
[8] Beier, T. et al. Multicut brings automated neurite segmentation
closer to human performance. Nature Methods 14, 101–102 (2017).
[9] Coelho, L. P. Mahotas: Open source software for scriptable
computer vision. Journal of Open Research Software 1(1):e3 (2013).
[10] Takemura, S.-y. et al. Synaptic circuits and their variations within
different columns in the visual system of drosophila. Proceedings
of the National Academy of Sciences 112, 13711–13716 (2015).
[11] Berning, M., Boergens, K. M. & Helmstaedter, M. Segem: efficient
image analysis for high-resolution connectomics. Neuron 87, 1193–
1206 (2015).
[12] Funke, J., Klein, J., Moreno-Noguer, F., Cardona, A. & Cook, M.
Ted: A tolerant edit distance for segmentation evaluation. Methods
115, 119–127 (2017).

9
0

1

U-N ET M ALA

U-N ET M ALA
U-N ET M ALA
U-N ET

U-N ET
LMC

RAND merge

VOI merge

LMC
CRunet

0.6

LFC
0.4

CRunet

0.4

LFC

0

0.5

1

1

2

1.5

1

0.6

0.8

0.2

0

1.5

U-N ET M ALA

0.2

VO;I

0.8

0.5

1

0.8

0.6

0.4

0.2

0

U-N ET M ALA

LMC

CRunet

LFC

LMC

CRunet

LFC

LMC

CRunet

LFC

RAND split

VOI split

(a) C REMI sample A+
0

1
U-N ET M ALA

U-N ET M ALA

U-N ET M ALA

U-N ET M ALA

0.2

U-N ET

LMC
RAND merge

VOI merge

LMC
CRunet

0.6

LFC
0.4

0.2

0

0

0.5

1

1.2

CRunet

0.4

LFC
0.6

1

0.8

0.8

1

2

1.5

1.4

U-N ET

VO;I

0.8

1

0.8

0.6

VOI split

0.4

0.2

0

U-N ET M ALA

RAND split

(b) C REMI sample B+
0

1

U-N ET M ALA

U-N ET M ALA
U-N ET M ALA

U-N ET
LMC

RAND merge

VOI merge

LMC
CRunet

0.6

LFC
0.4

CRunet

0.4

LFC

0

0.5

1

1.5

2

1.5

0.6
1

0.8

0.2

0

2

U-N ET M ALA

0.2

U-N ET

VO;I

0.8

1

1

0.8

0.6

0.4

0.2

0

U-N ET M ALA

RAND split

VOI split

(c) C REMI sample C+

Figure 8: Comparison of the proposed method against competing methods on the C REMI testing datasets A+, B+, and C+.
Shown are (from left to right) variation of information (VOI, split and merge contribution), Rand index (RAND, split and
merge contribution), and the total VOI (split and merge combined). Baseline U-N ET is our method, but without M ALIS
training (i.e., only minimizing the Euclidean distance to the ground-truth affinities during training) For U-N ET M ALA, the
red dot indicates the best threshold found on the training data.
[13] Saalfeld, S., Fetter, R., Cardona, A. & Tomancak, P. Elastic volume reconstruction from series of ultra-thin microscopy sections.
Nature methods 9, 717–720 (2012).
[14] Kingma, D. & Ba, J. Adam: A method for stochastic optimization.
arXiv preprint arXiv:1412.6980 (2014).
[15] Zlateski, A. & Seung, H. S. Image segmentation by size-dependent
single linkage clustering of a watershed basin graph. CoRR
abs/1505.00249 (2015). URL http://arxiv.org/abs/1505.00249.
[16] Januszewski, M. et al. Flood-filling networks. CoRR abs/1611.00421
(2016). URL http://arxiv.org/abs/1611.00421.

10

Figure 9: Reconstructions of 11 randomly selected neurons of the 100 largest found in the C REMI test volume C+.

11

Figure 10: Reconstructions of 23 randomly selected neurons of the 500 largest found in the F IB -25 test volume.

