Scalable Robust Graph and Feature Extraction for Arbitrary Vessel Networks in Large Volumetric Datasets
Dominik Drees, Aaron Scherzinger, René Hägerling, Friedemann Kiefer, Xiaoyi Jiang
11 Scalable Robust Graph and Feature Extraction forArbitrary Vessel Networks in Large VolumetricDatasets
Dominik Drees, Aaron Scherzinger, Ren´e H¨agerling, Friedemann Kiefer, Xiaoyi Jiang,
Senior Member, IEEE
Abstract —Recent advances in 3D imaging technologies providenovel insights to researchers and reveal finer and more detailof examined specimen, especially in the biomedical domain, butalso impose huge challenges regarding scalability for automatedanalysis algorithms due to rapidly increasing dataset sizes. Inparticular, existing research towards automated vessel networkanalysis does not consider memory requirements of proposedalgorithms and often generates a large number of spuriousbranches for structures consisting of many voxels. Additionally,very often these algorithms have further restrictions such asthe limitation to tree topologies or relying on the properties ofspecific image modalities. We present a scalable pipeline (in termsof computational cost, required main memory and robustness)that extracts an annotated abstract graph representation fromthe foreground segmentation of vessel networks of arbitrarytopology and vessel shape. Only a single, dimensionless, a-priorideterminable parameter is required. By careful engineering ofindividual pipeline stages and a novel iterative refinement schemewe are, for the first time, able to analyze the topology of volumesof roughly 1TB on commodity hardware. An implementation ofthe presented pipeline is publicly available in version 5.1 of thevolume rendering and processing engine Voreen . Index Terms —Graph Extraction, Centerline Extraction, Fea-ture Extraction, Vessel Network, Scalability
I. I
NTRODUCTION T HE study of vascular networks using modern 3D imagingtechniques has become an increasingly popular topic ofinterest in biomedical research [1], [2], [3], [4], [5], [6], as 2Dslice analysis is limited to a small subset of the available dataand cannot capture the 3D structure of vessel networks [7].At the same time, analysis of 3D datasets by visual inspectionis error prone [8]. Thus, for quantifiable results or novelinsights into data properties invisible to the human observer,automatic image processing techniques are required. In light ofdifferent imaging techniques, modalities and problem domains,the generation of a voxel-wise foreground segmentation, canbe seen as a sensible intermediate step for the automaticprocessing of vascular network images (Figure 1). However,the next step of calculating topological, morphological orgeometric features – a key requirement for the application inbiomedical research and beyond – has not received sufficient
D. Drees and X. Jiang are with the Faculty of Mathematics and ComputerScience, University of M¨unster, Germany.A. Scherzinger is with Ubisoft, Mainz, Germany.R. H¨agerling is with Charit´e – Universit¨atsmedizin Berlin, Germany.F. Kiefer is with The European Institute for Molecular Imaging, M¨unster,Germany.Corresponding author: Xiaoyi Jiang ([email protected]) attention from the research community [9], which is mostlyfocused on the segmentation domain [10], [11], [12].In recent years, improvements in imaging technology andprocedures are providing images of increasing resolution (interms of the number of voxels per unit volume) and thus offernew chances for finer grained analysis, but also pose newchallenges for image analysis algorithms. The analysis of finestructures down to the level of capillary networks promisesgreat insights. At the same time it is desirable to analyzelarge networks at once in order to obtain as much topologicalinformation as possible and avoid artifacts at the boundary ofthe volume. Modern microscopy hardware generates a singlevolumetric image file of hundreds of GBs or even TBs (thesame order of magnitude as hard drive capacity) per samplewhich existing vessel analysis approaches are not able toprocess. Due to rapid advances in imaging technology anddata sizes, this problem extends from commodity hardware tospecialized workstations. Consequently, for a vessel networkanalysis pipeline to be useful now and in the future, it has to be scalable both in terms of memory and runtime, and it shouldbe invariant with regard to increases in image resolution. Forapplication in a wide range of biomedical domains, a pipelineshould further not rely on specific imaging conditions, dataset-dependent parameters, a tubular vessel shape, isotropic imageresolution or a specific network topology.In this paper we present a pipeline engineered to fulfillthese requirements while extracting the topology, centerlinesand edge associated features from binary volumetric images.This is achieved by our main contribution , a novel iterativerefinement approach and careful design of all pipeline stages.In the remainder of this paper we will first discuss relatedwork and will then describe the proposed pipeline conceptu- CommunityFocus
ForegroundSegmentation FeatureExtraction
OurContribution
Fig. 1: Typical pipeline for processing of 3D images ofvascular networks. While the research community has mostlyfocused on the foreground segmentation task (e.g., [11]),this paper presents a widely applicable analysis based on(potentially very large) binary volumes as input. a r X i v : . [ c s . C V ] F e b ally. Next, we will elaborate the aforementioned requirementsbefore describing individual pipeline stages in more detailwith special attention to these requirements. Finally, we willquantitatively evaluate different aspects of proposed pipelinewith regard to the requirements and conclude our work.II. R ELATED W ORK
Chen et al. [13] use the neighborhood-based voxel classifi-cation definitions of [14] for analysis of hepatic vasculature.They also discuss other alternative voxel-based skeleton ex-traction methods, but note that they are not suitable: Distancetransformation-based methods [15] do not guarantee conn-tectedness of the extracted skeleton, while Voronoi-skeletonmethods [16] are comparatively time consuming. Furthermore,Chen et al. [13] perform extensive analysis to denote a singlevoxel of the skeleton as a branch point. They extract a graphby following skeleton voxels in a tree search, breaking upcycles in the graph. They do not perform any operations toremove erroneous branches – likely because for the processedlow resolution volumes the effect of noise is negligible.Drechsler and Laura [17] extend [13] by decomposing thegenerated skeleton into segments, computing a number ofproperties for each segment and generating a labeled graphstructure from the segments. The authors observe that theiralgorithm is very sensitive to noise in the foreground segmen-tation, but do not propose any means to reduce this effect.In [18] Chen et al. present an alternative thinning-basedskeletonization method for analysis of hepatic vasculature aspart of the previously published pipeline [13]. They use thevoxel classification methods of [14], but present a differentalgorithm. As a preprocessing step they suggest filling-holetechniques and morphological operations to remove cavitiesand irregularities on the surface of the object.Chothani et al. [19] propose a pipeline for tracing ofneurites from light microscopy image stacks. It comprisesforeground segmentation, skeletonization using a voxel codingalgorithm [20], tree construction using the scalar field of theskeletonization step and refinement including pruning basedon length as well as separation of branches that appear to beconnected due to limited image resolution. The separation ofbranches in 3D should not be necessary using higher resolutionimages. The authors note that in this case due to the increasedvolume size, advances in processing capabilities are required.Almasi et al. [21] present a method of extracting graphs ofmicrovascular networks from fluorescence microscopy imagestacks. In contrast to other methods operating on binary imagesthey use imaging method specific information to improve thedetection of vascular branches with disturbed image signal.Cheng et al. [22] analyze connected 3D structures such asvessel systems and metal foams in binary volumes. They usetopological thinning [23] and merge resulting branching pointsusing an algorithm based on regional maximal spheres andmaximal inscribed spheres with regard to the object surface.While their algorithm is linear in complexity, even for volumesof voxels (a gigabyte assuming one byte per voxel) itrequires hours of computation time.In [24], the authors use a non-topology preserving distancefield guided voxel thinning algorithm. They construct a graph Skeletonization Refinement
Done?
Topology Extraction
YesNo
Feature Annotation
Fig. 2: A schematic overview of the proposed pipeline.from an intermediate representation based on voxel neighbor-hood, but do not describe an efficient implementation. Thethinning step creates holes in irregularly shaped vessels whichis used to detect arteriovenous malformations, but makes isunsuitable for irregularly shaped (e.g., lymphatic) vessels.
Rayburst sampling [25] extracts precise measurements of ra-dius, volume and surface of spatially complex histopathologicstructures (given a centerline representation of a vessel). Themethod operates on grey value images and detects the surfaceby sampling along rays originating from the measurementpoint. For each point a large, but configurable number of raysis emitted to estimate the local morphology.Finally, as many vessel network analysis pipelines includea skeletonization step, the work of Bates et al. [26] is worthof note. The authors demonstrate binary vessel segmentationusing convolutional recurrent neural networks, including avariant centerline extraction. However, when compared totraditional skeletonization algorithms, the complex decisionprocess of ANNs has disadvantages: No hard guarantees ofconnectedness or the thickness of skeletal branches are avail-able (without further processing). Additionally, this approachmay have problems with high resolution datasets where thevessels are wider than the field of view of the network, thusmaking it impossible to determine the centerline position.All of the above works either do not mention processing ofvery large data samples (and thus likely do not support it) orexplicitly mention that large datasets pose a problem for theirmethod. In this work we want to fill this gap by presentinga widely applicable pipeline that is engineered to handle verylarge input volumes and presented in the next section.III. P
IPELINE
The proposed pipeline comprises four stages which are eval-uated iteratively (see Figure 2). In the first stage (
Skeletoniza-tion ) the foreground of the binary input volume is reduced toa voxel skeleton using a topological thinning-based algorithm,similar to [17]. We use a modified version of [14] whichcan be implemented efficiently both in terms of computationalcomplexity and disk access to the out-of-core dataset. Next,we build a graph representation from the voxel skeleton. Incontrast to [17] this
Topology Extraction stage is implementedas a single sweep over volume, engineered for memory localityresulting in fewer disc accesses. In the
Feature Annotation stage, we compute a set of morphological and geometricalproperties for all edges of the graph. This stage is subdividedinto two steps. First, the Voxel-Branch-Assignment determinesfor each foreground voxel which branch it is associated to.The Feature Extraction step uses this mapping to efficientlycalculate (geometrical [17] and additional morphological) fea-tures in single sweep over the volume. Then, in the
Refinement stage, the graph is pruned using previously calculated features,using the scale invariant, dimensionless property bulge size.Since pruning invalidates the Voxel-Branch-Assignment, thefirst three stages of the pipeline have to be reevaluated. Asshown in evaluation, this step is essential for application onhigh resolution images. This Extraction-Refinement cycle isrepeated until a fixed point is reached. We want to stress thatthe novel iterative refinement approach is essential to obtainmeaningful results from very large datasets, as demonstratedin Section VI.IV. R
EQUIREMENTS FOR B ROAD B IOMEDICAL A PPLICATION
As discussed before, the proposed pipeline, in contrast toprevious work, is engineered to fulfill a set of requirements thatare essential for a wide range of biomedical applications. Forthe discussion in the remainder of this paper, in this sectionwe will structure, motivate and expand on the requirementsbriefly mentioned in the introduction.
Primary Requirements:
The method should be scalablewith regard to runtime (P1). As all voxels of a volume can beexpected to be taken into account by the method, the runtimecan be expected to be at least linear in the number of voxels n . In order to be applicable to increasingly large datasets,a runtime of O ( n ) would be unacceptable, but quasilinearalgorithms (e.g., in O ( n log n ) ) can still be executed. Beyondcomputational complexity, special care has to be taken forexample in terms of memory access locality to avoid anincrease in computation time by a large constant factor.The method should be scalable with regard to memoryrequirements (P2). We can expect the input datasets to fit onto(commodity) non-volatile storage (disk), but not necessarilyin main memory. While the price per capacity of both disksand main memory fall exponentially, the rate of (exponential)decrease is higher for disks than for memory [27], [28].Consequently, the ratio of average disk size to average RAMsize grows over time. The exact relationship is not easy toestablish as it changes over time [28]. Here we assume thatdisk and volatile memory size roughly grow in a relationshipof m to m for a disk size of m . Therefore the required mainmemory size of the method should not exceed O ( m ) .The method should exhibit invariance with regard toimage resolution (P3). For a fixed size specimen, an increasein image resolution results in a dramatic increase in surfacenoise related artifacts for simple topological thinning-basedmethods which is unacceptable for analysis of large datasets. Secondary Requirements:
From the desire to apply themethod in practice we derive a set of further requirementsthat are not directly related to scalability:The method should be unbiased (S1), i.e., not dependent ona set of parameters that are required to be selected carefullyand correctly depending on the input image.The method should be robust in terms of deviations fromthe tube shape (S2) of the analyzed network structure. Whileblood vessels are tubular, lymphatic vasculature as an exampleis highly irregular, but should still be processed correctly.The method should handle volumes of anisotropic res-olution (S3). Volumetric imaging techniques often create
Skeletonization
Feature Annotation
Feature ExtractionVoxel-Branch-AssignmentRefinement
VesselgraphBinary VolumeID VolumeProto-Vesselgraph
Topology Extraction
Fig. 3: A schematic overview of the data flow in the proposedpipeline. Dashed lines signify an optional dependency on thedata for the particular stage. For skeletonization topologicalinformation is not available in the first iteration of the pipeline.datasets with voxel spacing that differs between coordinateaxes. As an example, in light sheet microscopy datasets areconstructed from a series of (2D) slices for which the spacingis independent from the x-y-resolution. Methods operating onthe voxel grid of a volume have to consider this.The method should be able to fully analyze networks of arbitrary topology (S4). Existing methods often assume atree structure and either fail to operate on other topologiesor drop information [13]. However, for lymphatic vessels andcapillaries, or even larger structures ( cerebral arterial circle )the assumption of a tree topology is invalid.The method should be able to analyze images indepen-dently of the imaging conditions (S5). Concrete imagingconditions (i.e., distribution of gray values or fluorescentstaining techniques in microscopy imaging) vary betweendomains. In order to be widely applicable, the method shouldnot be dependent on concrete models of imaging conditions.Our method achieves this by operating solely on binary inputdata and leveraging the broad range of other research on vesselimage segmentation [10], [12].V. M
ETHOD D ETAILS
In this section, we will present the algorithmic details of theproposed pipeline with special attention to the requirements forbroad biomedical application. Figure 3 provides an overviewof the data flow between individual stages of the pipeline.
A. Skeletonization
Like other vessel network analysis pipelines [13], [17], weuse a modified version of the established topological thinningalgorithm by Lee et al. [14] which has the advantage thatany voxel can be evaluated in terms of these properties byonly considering its 26-Neighborhood, which is advantageouswhen operating on very large volumes. Furthermore, we arenot aware that other approaches solve the inherent problemsof analysis of large datasets. In the original formulation, thealgorithm iteratively removes voxels that lie on the surface andwhose deletion does not change the topology of the object in6 (static) subiterations until no more voxels can be removed.For an efficient implementation, further considerations arenecessary. We model the surface of the object explicitly asa sequence of voxel positions, which is initialized by scan-ning the volume before the first iteration step. In subsequentsubiterations, we build the active surface (voxels that can potentially be deleted in the next iteration) by retaining voxelsfrom the previous active surface that are not considered duringthe current subiteration and adding all foreground voxels inthe 26-Neighborhood of a voxel after its deletion. If a voxelwas considered, but not deleted during the current subiteration,it will be removed from the active surface even althoughit is still part of the surface of the object. If one of itsneighbors is deleted, it will be re-added to the active surfaceand reconsidered for deletion in the following iteration. Thisimplementation has a runtime of O ( n ) in the number of voxelsin the volume (P1): In each subiteration, only voxels in theactive surface are considered. As a voxel is either deletedentirely or removed from the active surface after each iteration(i.e., 6 subiterations) and only added again once one of its26 neighbors is deleted, it will be considered for deletiona constant amount of times. In order to fulfill requirementP2, the volume is stored on disk and dynamically mapped tomemory using operating system capabilities. Constant runtimeimprovements could be observed using a compressed (2 bitsper voxel) representation and storing the volume as 32x32x32(8 KB) blocks in linear memory, thus reducing disk accessby exploiting the spatial locality of surface voxels. Activesurfaces are stored as on disk (linear) voxel positions insorted sequences. During a subiteration the previous activesurface is read from front to back. Simultaneously, the newactive surface is built slice-by-slice in memory by collectingpositions, and sorting and removing duplicates before writingto disk, requiring O ( m ) main memory (P2).To account for volumes with anisotropic resolution, we keeptrack of the real-world depth of voxel layers deleted for eachdirection and select the direction with the lowest total depth asthe direction for the next subiteration. This way, the speed ofvoxel deletion is equalized on average for highly anisotropicvolumes (S3). If the last subiterations of all 6 directions (whichmay have happened out of order) did not delete any voxels, theactive surface is empty and the skeletonization is completed. B. Topology Extraction
Previous work [13], [24] does not describe an efficientimplementation of topology extraction for large input datasets.For example, Chen et al. [13] perform a tree search followingskeleton voxels through the volume. This is unsuitable forlarge volumes as it either requires keeping the completevolume in memory or involves very frequent random access tohard disk with high latency. Additionally, their method breaksloops in the graph and always creates a tree topology.In constrast, the topology extraction in the proposed pipelineextracts the complete centerline graph in a single pass over thevolume using a modified version of the streaming connectedcomponent finding algorithm by Isenburg and Shewchuk [29].Instead of all foreground components, we consider the threeskeleton voxels classes (regular (2 neighbors), end ( < neigh-bors), and branch ( > neighbors) [13]) separately and extractcomponents for each. Individual connected components of endor branch voxels make up nodes in the graph. Connectedcomponents of regular voxels form lines within the volumeand have (except for the separately handled case of closed (a) (b) (c) (d) Fig. 4: If the centerline of a small vessel (blue) is closer tothe surface of a larger vessel (yellow) than its own centerline,some voxels of the larger vessel will be assigned to thesmaller vessel (a). After the initial Voronoi Mapping, con-nected foreground components of the same ID are remappedusing IDs from seed points of branches (b). Unlabeled regionsare identified (c) and flood-filled from labeled voxels (d).loops) exactly two end points. Unlike [13], we do not force theposition of nodes to be defined by a single voxel and insteaduse the barycenter of the connected region of voxels. Usingthe two end points of a segment of regular voxels we can findthe nodes that are connected by the path of regular voxels andthus construct an edge. To perform this node-edge-matchingprocess efficiently for large graphs, for each node voxel (i.e.,an end or branch voxel) we insert a reference to the originalnode into a (static, on-disk) k-d tree. Thus, for each run ofregular voxels we can efficiently query the position of nodevoxels that are 26-neighbors of its tips, and thus link node andedge data structures. Extracted nodes and edges (includingcenterlines) form the
Proto-Vesselgraph , are each stored infiles on the disk and are dynamically mapped to memory.Furthermore, the algorithm of [29] (and this modification)can be shown to only require O ( m ) main memory (P2).Moreover, we extract the whole graph in a single pass overthe volume in O ( n log n ) time (P1) and do not make anyassumptions about the topology of the extracted network (S4).Since the individual branch voxels necessarily lie on voxelpositions of the original volume, they resemble a raggedcenterline. This both artificially increases the line length andproduces a larger number of ambiguous cases in the Voxel-Branch-Assignment. We therefore smooth all centerlines of theProto-Vesselgraph using local bezier curves. C. Voxel-Branch Assignment
Drechsler and Laura [17] calculate the volume property ofedges by assigning voxels of the foreground segmentation tothe nearest centerline point. As illustrated in Figure 4a, thiscreates incorrect results in some cases. While these errors arepotentially not as severe for the calculated volume , morpho-logical properties based on the radius of the vessel can beheavily affected. In order to correctly calculate edge-associatedproperties for the previously created Proto-Vesselgraph struc-ture, we therefore first create a mapping between voxels ofthe foreground segmentation and the edges of the Proto-Vesselgraph. This process is performed in 4 steps that aredescribed below and illustrated in Figure 4.
1) Centerline Voronoi Mapping:
As a first approximationof the edge ID map, we find the nearest centerline point foreach voxel of the foreground segmentation. We accelerate thesearch using a static on-disk k-d tree of all centerline points of all edges of the Proto-Vesselgraph, each annotated with thecorresponding edge ID. As we iterate over all n foregroundvoxels and perform a log n search in the k-d tree for each, thetotal runtime is within O ( n log n ) , as is the construction of thek-d tree (P1). We iterate over the whole volume in 32x32x32blocks of voxels for spatial locality of lookups, thus reducingthe need to randomly reload parts of the k-d tree from disk.
2) Connected Component Remapping:
As illustrated inFigure 4a, matching voxels via minimal euclidean distancedoes not yield a proper voxel wise vessel segment map in allcases. When two vessels with highly dissimilar radii lie closeto each other, some voxels of the larger vessel will be ascribedto the smaller vessel. These cut off regions have to be identifiedand remapped. We perform a modified connected componentanalysis [29] on the generated edge ID volume in which weconsider two voxels mergable if they have the same ID. Thenwe map the connected component IDs back to the edge IDsusing a table constructed by sampling at one centerline pointof each edge in the result of the connected component analysis.
3) Cut-off Region Identification:
Cut-off regions do nothave an associated edge ID and can therefore be identifiedusing another pass of a connected component analysis [29]where only voxels with undefined IDs are considered. In asingle pass over the whole volume, we collect the axis alignedbounding boxes for each cut-off region (see Figure 4c).
4) Cut-off Region Flooding:
Previously identified cut-offregions are relabeled by flooding all voxels without edge IDsstarting from adjacent voxels with valid edge IDs. For thesimple case of a region cut-off from a single vessel section,this fills the region with the section’s edge ID. In a more com-plicated case of a cut-off region adjacent to multiple differentvessel sections, the resulting labeling approximates the L -Voronoi-cells spanned from the surface of adjacent regions. Aside effect of this procedure is that foreground regions that donot resemble (complete) vessels at the boundary of the volume(and thus do not contain any centerline voxels), are not floodedwith valid IDs in this step and ignored in subsequent steps.For efficiency, all identified disconnected regions are pro-cessed separately as cuboidal subvolumes and collected duringa single pass over the ID-volume. Currently valid boundingboxes of cut-off regions are organized in interval trees. For thez-axis, a single tree references all bounding boxes. For eachslice, the active regions can be queried to construct an intervaltree for the y-axis. For each line in the slice, using the y-axistree an interval tree for the x-axis can be constructed. For eachvoxel in the line, all active regions can be queried from thex-axis tree so that the current voxel value can be written tothe corresponding subvolumes. After this collection step, eachsubvolume is flooded separately. Similar to the skeletonizationstep, we use the concepts of an active surface, slice-wiseprocessing and memory mapping of files on disk to guaranteelinear runtime (P1) and O ( m ) memory requirement (P2).The results are written back to the ID volume in a processsimilar to the collection of the subvolumes. The number ofdisconnected regions is bounded by n . Consequently, theconstruction of the range trees and the querying for each voxelare within O ( n log n ) (P1). All subvolumes and range trees can be stored on disk and mapped dynamically to memoryand therefore do not require additional memory (P2). D. Feature Extraction
Using the Proto-Vesselgraph and the generated edge IDvolume, we can compute the same edge properties as Drechslerand Laura [17]. Some geometric features can be computeddirectly from the Proto-Vesselgraph. These include length (arc-length of the centerline), distance (euclidean distance ofconnected nodes) and straightness ( distancelength ). Drechsler andLaura [17] propose curveness ( lengthdistance ), but as the distance of any vessel segment is guaranteed to be smaller than its arc- length , we argue that straightness (with a guaranteed domainof [0 , ) is superior to curveness ( [1 , ∞ ) ).For other features we collect information from the edgeID volume and the Proto-Vesselgraph. For each foregroundvoxel in the edge ID volume, we find the corresponding edgeand query the closest of its centerline point(s). For efficientlookups, the centerline points are organized in a (static, on-disk) k-d tree for each edge. We equally distribute and storethe volume occupied by the foreground voxels alongside thecorresponding centerline point(s). Surface voxels are assignedto the closest centerline point, for which thus the minimum,maximum and average surface distance is computed.Values collected for individual centerline points can beaggregated to compute more per-edge features: For each ofthe surface distance-based values ( minimum , maximum , average , roundness := minimummaximum ) we compute the meanand standard deviation, totalling in 8 additional morphologicalfeatures not discussed in [17]. Compared to [25] this provides arelatively coarse understanding of the local vessel morphology,but is a good approximation that can be computed efficiently.Finally, the per-voxel volume can be accumulated to computethe total volume occupied by each vessel segment in the vesselgraph. The average cross-section of a vessel is obtained bydividing the volume by the length [17].In total, the number of centerline points is obviouslybounded by n so that the cost of building and querying thek-d trees is within O ( n log n ) . Additionally, the label volumeis accessed in 32x32x32 chunks improving memory localityof searches in the k-d trees (P1). All components of the Proto-Vesselgraph, the k-d trees and the created graph are stored ondisk so that main memory requirements are met (P2). E. Refinement
Skeletonization algorithms tend to produce a number ofsmall spurious branches, especially for vessels with irregularsurfaces. We therefore propose to refine the generated vesselgraph by pruning spurious branches. A simple, but scale-dependent way of defining deletability is to set a globalminimum length [19]. However, this is problematic if vesselsof different scales are present in a single dataset. Instead,we propose bulge size as a scale-independent dimensionlessmeasure (S1). Intuitively, the bulge size measures how far abump, bulge or branch has to extend from a parent vessel inorder to be considered a separate vessel. This size is expressedrelative to the radius of its parent vessel and itself, making it (a) bulge size ≈ (b) bulge size ≈ (c) bulge size ≈ avgRadiusMeanlengthinner lengthtip radius (d) legend Fig. 5: Schematic depiction of the calculation of the bulge size feature of an edge for three examples.scale-independent. More formally, the bulge size is an edgefeature, that is computed during the feature extraction V-D,and is only defined for bulging edges , i.e., edges that connecta leaf node (degree 1) and a branching point (degree > ).For all centerline points of the edge, we decide whetherthey are within a branching region ( inner points ) or not ( outerpoints ) during the feature extraction. If a point has associatedsurface points, neither of which are adjacent to surface voxelsof another edge, it is condidered an outer point. An inner pointis characterized by either not having any associated surfacevoxels or by having a surface voxel that is a 26-neighbor ofa point that belongs to another edge. The inner length of anode is the arc length defined by a run of centerline pointsstarting from that node and ending at the last inner voxel.This corresponds to the length of the piece of centerline thatlies within a branching area. The tip radius is defined as theminimum distance to the surface measured from the centerlinepoint closest to the node. Without loss of generality, let n b and n e be the branching point node and leaf node, respectively. bulge size ( e = ( n b , n e ))= length ( e ) − inner length ( n b ) + tip radius ( n e ) avgRadiusM ean ( e ) This definition provides a dimensionless measure of shape thatperforms well even for corner cases and can be computedefficiently during the feature extraction stage of the pipeline.The calculation and applicability is illustrated in Figure 5.The pruning of spurious branches is a node-oriented, itera-tive approach outlined in Algorithm 1. A pruning step iteratesover all nodes in the graph and collects deletable branches.All deletable edges and nodes that were only connected tothose edges are then removed from the graph. Additionally,two edges that share a common node of degree 2 are mergedby connecting their centerlines and recomputing propertiesfrom the attributes of points in the combined centerlines. Thepruning step is repeated until a fixed point is reached.While the refined graph does not retain any signs ofdeleted edges, the simultaneously computed centerlines andedge properties are not preserved by the refinement procedure.This happens because the volume regions associated with nowdeleted branches do not contribute to the properties of the
Algorithm 1
Graph refinement algorithm function REFINE ( G, t ) repeat D = {} for node n ∈ G do E = edges connected to n P = {} for e ∈ E do if bulging ( e ) ∧ bulge size ( e ) < t then P = P ∪ { e } end if end for while | E | − | P | < do (cid:46) Retain two edges P = P − { argmax e ∈ P bulge size ( e ) } end while D = D ∪ P end for Delete edges in D and remove orphaned nodes Remove nodes with degree 2 and merge edges until D = {} (cid:46) Iterate until reaching fixed point return G end function remaining branches until the ID map and features are recom-puted. While this affects the bulge size, the pruning schemewill not erroneously delete branches, since the missing regions(lowered radii) cause the bulge size to be overrestimated . F. Iterative Scheme
In order to ensure that centerlines and edge properties matchthe refined graph, we employ an iterative scheme where thepreviously extracted and refined graph is fed back into the firststage of the pipeline to improve the generated results.We modify the skeletonization algorithm to force it togenerate a voxel skeleton with a topology that matches therefined graph of the previous iteration. All voxels of nodeswith degree 1 in the graph are set to be fixed foregroundduring the skeletonization of the input volume. This meansthat they are never considered for deletion. These voxels markthe beginning/end points of lines to be extracted during theskeletonization. The skeletonization algorithm is modified tonot preserve voxels at the end of voxel lines (non- line voxels [14]). The resulting voxel skeleton connects all previouslyextracted endvoxels with medially positioned lines.The iteration can be stopped if it reaches a fixed point. Thisis the case if two consecutive iterations generate the samegraph. As the number of edges never increases, the checkcan be reduced to a simple integer comparison. Optionally,an upper limit for the number of iterations can be specified.VI. E
VALUATION AND D ISCUSSION
In this section we evaluate and discuss the proposed pipelinein terms of the primary and secondary requirements forbiomedical application (Section IV). All experiments wereconducted on a consumer grade PC (AMD Ryzen 7 2700X (a) Original Image (b) Resample Strategy (c) Mirror Strategy
Fig. 6: An exemplary demonstration of the resample (b) mirror (c) strategies for increasing volume size on a (2D) dataset (a)to scale 2, i.e., doubling the volume size in each dimension.For both strategies, the foreground (grey pixels) still form aplausible foreground segmentation of a vessel network.(3,7 GHz), 32GB RAM, and 1 TB hard disk (Samsung NVMeSSD 960 EVO)).
A. Runtime Scalability
The demonstration of scalability requires some way ofmodifying a scale parameter for a given dataset withoutchanging other characteristics. One possibility would be touse an already large dataset, scaling it down or clipping it to asmaller region. However, besides the problem of availablilityof a very large volume (representing the results of future microscopes), this could cause artifacts and loss of detail dueto downsampling, making the comparison of generated graphsvery difficult. Instead we use a small (real world) volume andartificially increase its size using two strategies (see Figure 6):The resample strategy: Resampling the volume with largerresolution, using nearest-neighbor sampling to avoid changingtopology near very thin connections in the original volume.The mirror strategy: Repeating the volume in each dimen-sion until the desired (integer) scale is achieved. In order togenerate a mostly connected network, in each dimension the i + 1 ’th volumes are mirrored.In a real world scenario, improved acquisition technologywould result in a combination of both increasing (voxel) size ofpreviously visible vessels and revealing previously undetectedfiner networks of larger complexity. When not specified other-wise, the dataset Lymphatic 1 (Figure 17a, × × vox-els) was used for evaluation, as it has a non-tree-like topology,anisotropic voxel resolution ( µ m × µ m × µ m ) and airregular vessel shape. In the following, scale denotes thefactor that each dimension of the dataset was multiplied withusing one of the above strategies.Figure 7 shows the average runtime of a single refinementiteration of the pipeline depending on the number of voxelsfor the mirror and resample strategies in a log-log plot.As shown, the runtime is only slightly worse than linear,as should be expected in the light of the derived runtimecomplexity of O ( n log n ) . Furthermore, it should be noted thatthe increase in slope around coincides with the exhaustionof main memory and thus a larger number of disk accesses(cf. Figure 11).Figure 8 once again shows the average iteration runtime,but using a linear scale and with details of individual stepsof the pipeline. It can be observed that the runtime of allsteps increases roughly linearly with respect to the number Number of Voxels10 A v e r a g e I t e r a t i o n R un t i m e [ s ] Scale Mirror c*nlog(n)c*n
Fig. 7: Average iteration runtime of the pipeline for the resample and mirror strategies in a log-log plot. The functions c · n log n and c · n are shown as visual guides. c was chosen sothat both guides match the mirror plot at the first data point. A v e r a g e I t e r a t i o n R un t i m e [ s ] SkeletonizationTopology ExtractionVoxel-Branch-AssignmentFeature Extraction (a) Resample Strategy A v e r a g e I t e r a t i o n R un t i m e [ s ] SkeletonizationTopology ExtractionVoxel-Branch-AssignmentFeature Extraction (b) Mirror Strategy
Fig. 8: Iteration runtime of individual stages of the pipelinefor the resample (a) and mirror (b) strategies. Note that theexecution time of the Refinement step is not shown as it isvery low compared to all other stages.of voxels. Here, it becomes apparent that for the resample strategy, an iteration appears to take slightly more time thanwith the mirror strategy. This appears to be mostly due tomore time spent in the Voxel-Branch-Assignment and FeatureExtraction steps – likely because larger vessels result in moretime spent in searching the k-d-trees in both steps.The total runtime of the pipeline, however, is not onlydependent on a single refinement iteration, but also on thenumber of iterations required to compute the final result. Thenumber of iterations required to reach a fixed point for anincreasing number of voxels is depicted in Figure 9 for bothscaling strategies. While the total number of iterations usingthe resample strategy can expected to be constant for scale 4and higher (since scale 8 corresponds to 8 copies of the datasetat scale 4), the number of iterations remains at 4 independentof the scale in this experiment. While we do not derive N u m . o f R e f i n e m e n t I t e r . Scale Mirror
Fig. 9: The total number of iterations required to reach a fixedpoint for the resample and mirror strategies. N u m b e r o f N o d e s i n G r a p h Scale181624 32404856
Fig. 10: The total number of nodes in the graph after differentnumbers of iterations for the various scales.a maximum for the number of iterations for the resample strategy, we can observe the required iterations to increaseonly moderately with the number of voxels in the experiment.While we acknowledge that the (total) processing time oflonger than one week for a 1 Terabyte dataset is certainlyfar from interactive use, this does not pose a problem inbiomedical research since the time required for the preparation,image aquisition and postprocessing of a sample is of the samemagnitude, and the presented pipeline is to the best of ourknowledge the first to make analysis of these samples possible.Furthermore, as our pipeline can be considered unbiased (seeSection VI-D), interactive parameter tuning is not required.Figure 10 shows that the total number of nodes in thegraph rapidly decreases even for higher scales. The finalnumber of nodes for all scales is very close to scale 1,but increases slightly for higher iterations. Section VI-C andFigure 12 demonstrate that the resulting graphs are actuallysimilar. Thus, if constant runtime (for different datasets of agiven scale) is required, a relatively low maximum numberof iterations can be specified to obtain a good approximationof running until a fixed point is reached. At the same time,Figure 10 shows that a single refinement step is not sufficientto obtain meaningful results for very large datasets, as alarge number of spurious branches remain. This underlinesthe necessity of the presented iterative refinement approach.
B. Main Memory Scalability
In the description of the pipeline we have argued for allsteps of the pipeline to allocate O ( m ) of memory. For thatreason and because total allocated heap memory is difficult tomeasure efficiently, we demonstrate how the implementationof the presented pipeline behaves with regard to the resident R e s i d e n t S e t S i z e [ G B ] Scale Mirror
Fig. 11: Maximum resident set size (RSS) of the graphextraction process for different problem sizes (number ofvoxels in the volume) for the resample and mirror strategies. set size (RSS) , which specifies the portion of the memorymap of a process that is actually held in main memory. Thisexcludes allocated memory that has been swapped out todisk, but includes (sections of) files that have been copied tomemory. Figure 11 shows the maximum RSS of the processperforming the graph extraction for different problem sizes(number of voxels in the volume). It can be observed that theRSS increases for higher problem sizes, but does not exceed alimit near the total available main memory of the test machine(32GiB). This demonstrates the aforementioned property of thepipeline to be scalable in terms of memory, but to use all ofthe available memory resources. Notably, the mirror strategyreaches the 32GiB limit earlier than the resample strategy. Thismay be due to the fact that by repeating the graph m -timesin each dimension, the total number of centerline points inthe graph (and thus the required memory to store it) increasesroughly by a factor of m , while for the resample strategy, thenumber of centerline points is (roughly) multiplied by m . C. Image Resolution Invariance
As a qualitative evaluation, Figure 12 demonstrates theeffect of increased image resolution on current skeletoniza-tion/graph extraction approaches and how the iterative refine-ment approach solves this problem. For a 16-fold resolutionincrease in each of the coordinate axes (using the resamplestrategy), the number of branches and nodes increases sharply(Figure 12b) compared to the the final graph extracted fromthe original volume (60-fold increase for the depicted example,Figure 12a). After 5 refinement iterations, the number ofbranches is reduced (Figure 12c) to a number comparableto the simple case (scale 1, Figure 12a). As the increase inresolution allows for finer grained skeletonization, the graphsstill differ in some details, but most of the structure is thesame. This is confirmed by biomedical domain experts whoconsider the refined version appropriate for further analysis,but reject the graph without refinement.Figure 13 shows how an increase in resolution (using the resample strategy) affects the extracted graph. Since no groundtruth graph is available for the dataset, we compare the allavailable graphs with the graph extracted using the lowest and highest scale, respectively. When using the graph ofscale 1 as a template, the edge match ratio drops relativelysharply for increasing, but smaller scales. For higher scales thedifference in similarity is not as pronounced. An explanationcould be that for low scales (and thus low image resolutions), (a) Scale 1, iterative refinement (b) Scale 16, no refinement (c) Scale 16, iterative refinement
Fig. 12: A demonstration of the effect of increased image resolution on the intermediate result and how the refinement iterationsmitigate this effect, using dataset Lymphatic 1: An increased image resolution highly increases the number of erroneous branches(in this example the total number of branches is increased roughly 60-fold). After 5 refinement iterations the resulting graphis very similar to the graph extracted from the original volume without artificially increased resolution. E dg e M a t c h R a t i o Template: Scale 1 Template: Scale 56
Fig. 13: A similarity comparison between graphs of differentscales using the edge match ratio [30] and the graphs extractedfrom the lowest and highest scale volume as template graph.skeletonization artifacts more heavily influence the final result.Small bulges with a bulge size near the specified thresholdcannot be treated with as much precision as necessary and arethus deleted although they should not be and vice versa. Thishypothesis is supported by the fact that when choosing thegraph extracted from the highest scale volume as the templatefor comparison, the graphs for most scales exhibit a muchhigher similarity.As an additional way of evaluating the image resolutioninvariance (P3), we want to focus on one of the primaryproblem of increased image resolution: The increase in surfacenoise on the individual voxel level. For this we generate10 datasets using Vascusynth [31], a tool for generatingvolumetric images of vascular trees as well as the corre-sponding ground truth segmentations and tree hierarchy usedhere. Before graph extraction we perturb the surface of thebinary volumes by iteratively selecting a random foregroundor background surface voxel and flipping its value if this doesnot change the topology of the object. The surface noise levelis defined as . An example of a volumewith applied surface noise is presented in Figure 14b. Foreach of the volumes, four variants with different random noiseseeds were evaluated. As shown in Figure 15, the edge matchratio [30], i.e., the ratio of edges that can successfully bematched between the extracted and the ground truth graph,rapidly decreases when not using iterative refinement whilethe refinement approach presented in this paper manages toreproduce the expected graph even for high noise levels with (a) Synthetic 1 (no noise) (b) Synthetic 1 (noise level 0.2)
Fig. 14: The foreground of one of the synthetic vessel datasetsused for evaluation without (a) and with maximal noise (b)applied. In (a) a very small bump which is categorized as aseparate branch according to the ground truth is highlighted. E dg e M a t c h R a t i o no refinement iterative refinement Fig. 15: The edge match ratio [30] between the groundtruth graph and the graph extracted from the correspondingvolume [31] with and without iterative refinement. Prior tothe graph extraction salt-and-pepper noise is added to thesurface (x-axis). For each noise level, minimum, maximumand average of 10 datasets for 4 surface noise seeds are shown.only slight decreases in similarity. Note that we chose a bulgesize of 1.5 for the iterative refinement, which is unusually lowfor blood vasculature. However, Vascusynth generates somebranches that are very short compared to the radius parentvessel, and in some cases are entirely enclosed, which areonly visible as very shallow bumps (or not visible at all) in thegenerated volume (see highlighted region in Figure 17a as an (a) bulge size=0.2 (b) bulge size=0.3(c) bulge size=1.0 (d) bulge size=3.2 Fig. 16: A demonstration on how different parameter valuesaffect the result of the graph extraction pipeline. In theexample, very low values allow even very minor bumps (suchas the cow’s teats in the model) to be considered separatebranches. Increasing values gradually remove further branchessuch as the belly, tail, and features of the head until finally fora very large value the graph consists of a single edge.example). For real world applications this is likely negligible.Instead, higher robustness is likely preferable, necessitatinghigher bulge size, as discussed in the following section.
D. Influence of the bulge size
Parameter
As explained in Section V-E, the parameter of the proposedmethod can be selected a-priori based on the application andthe desired outcome of the pipeline. The quality of the resultdoes therefore not directly depend on the parameter. Hence, weconsider our method unbiased (S1). Nevertheless, in Figure 16we present an overview of the influence of different bulgesize values on the result of the pipeline using a test datasetwith low-profile, but variable-depth bulges. As can be seen,a very low bulge size results in separate branches even forvery small surface features. Increasing the parameter graduallyreduces the complexity of the extracted topology. For realworld applications, the bulge size should be set a-priori basedon knowledge of the dataset. For example, lymphatic vessels(especially in the pathological case, illustrated in Figure 17b,but even in healthy individuals, see Figure 17a) often havecomparatively short branches, but also suffer from irregularsurface features, often near branching points. We thereforesuggest a bulge size of 1.5, where the edge case correspondsto nearly spherical, but slightly elongated bulges. On the otherhand, (healthy) blood vasculature often exhibits a smoothsurface, round diameter, and clearly defined branches, so thata bulge size of 3.0 or higher can be chosen confidently.
E. Anisotropic Resolution
The third secondary requirement states that a pipelineshould be able to operate on volumes with anisotropic res-olution (S3). In order to evaluate our proposed pipeline in (a) Lymphatic 1 (b) Lymphatic 2
Fig. 17: Some exemplary real world datasets used for evalua-tion rendered as the surface of the foreground. It is apparentthat the surface and topology of the real world, lymphaticvessel datasets are much more complex than that of thesynthetic blood vessel datasets (see Figure 14).terms of that property, we use a volume with known groundtruth and isotropic voxel spacing. We then artificially createdatasets with different, anisotropic resolutions by resamplingthe dataset in specific dimensions. As an example, confocalmicroscopes due to their anisotropic point spread function gen-erate volumes with a comparatively low z-resolution, which iswhy in this experiment we chose to leave the z-axis resolutionintact and increase the resolution in x- and y-direction (akin tothe resample strategy). In total, we evaluated 10 datasets withresolution scale differences of 1, 2, 4, 8, and 16. We comparedthe extracted topology using the graph matching method andedge match ratio similarity measure proposed in [30] and thequality of the extracted centerline geometry using the NetMetsframework [32]. In the case of NetMets the average radius ofall edges multiplied by 10 was chosen as the parameter σ .The results are shown in Figure 18. Although the groundtruth is not matched perfectly and some drift between scales isnoticeable (see discussion about ground truth quality above),in general no trend towards positive or negative impact onthe evaluated score can be observed based on the level ofanisotropy. Furthermore, it should be noted that the dynamicaxis selection discussed in Section V-A indeed improves theresult for larger anisotropy, although for smaller levels thereis next to no difference between both variants. We thereforeconclude that the proposed pipeline is in fact able to processdata of even highly anisotropic resolution. F. Robustness against shape deviations
For the evaluation of robustness of the pipeline, especiallywith regard to deviations from the tube (S2) shape of analyzedvasculature, we use the robustness test introduced in [30].Additionally, the FNR and FPR values (false negative/positiveratio) of [32] were used to measure the geometric error usingin the framework of [30]. The (real world) lymphatic andsynthetic blood vessel datasets were rotated and resampled(with doubled resolution in each axis) in ° steps around the z -axis, processed using the proposed pipeline, compared to thetemplate graph. The minimum similarity of all steps denotesthe robustness index. For the synthetic datasets the groundtruth graph provided by the tool was used as a template and E dg e M a t c h R a t i o dynamic axis selection no dynamic axis selection N e t M e t s F N R dynamic axis selection no dynamic axis selection1 2 4 8 16Spacing Scale Difference0.000.050.10 N e t M e t s F P R dynamic axis selection no dynamic axis selection Fig. 18: A demonstration of how anisotropy affects the topo-logical (Edge Match Ratio [30]) and geometrical (Center-line [32]) structure of the extracted graph. We compare theground truth graph of a volume generated by VascuSynth [31]to the graph extracted using the proposed pipeline after in-creasing the x- and y-resolution by a specified factor (but leav-ing z-axis resolution unchanged) and resampling the volume.Minimum, average and maximum for 10 datasets are shown.As can be seen, anisotropic resolution does not strongly affectthe pipeline and the dynamic axis selection is advantageous.thus are also measuring the accuracy of the method. For realworld datasets, generation of an annotated ground truth graphis virtually impossible [30] and thus only the robustness canbe evaluated. We compare the results without refinement tothe graphs extracted using a bulge size of 1.5 and iterationuntil reaching a fixed point.The aggregated results in Table I show that for all prop-erties the refinement procedure results in significantly higherrobustness for all datasets. As should be expected, the totalrobustness of the much more complex real world lymphaticvessel datasets is overall lower than for the (simpler) syntheticdatasets. However, it is noteworthy that even when applied tocomplex real world datasets, the proposed pipeline achieves ahigher score than the simple, but commonly applied approachwithout refinement on simple, synthetic datasets.
G. Remaining Secondary Requirements
The fourth secondary requirement (the ability to processvolumes of arbitrary topology) directly follows from how thetopology is extracted from the binary skeletonized volume(Section V-B) and, in particular, is not restricted to a treetopology, in contrast to other methods [17]. Furthermore, thepipeline is independent of the imaging conditions (S5) byoperating on an existing foreground segmentation.
H. Limitations
The proposed pipeline – like all topological thinning-basedmethods – by definition suffers from distortions in the binaryinput image that causes changes to the topology, i.e., cavitiesin foreground objects and loops on the boundary. While thesefeatures are not expected in the actual vessel structure, imagingartifacts, noise or problems with the segmentation procedurecan still produce such artifacts in practice. However, carefulpost-processing of the segmentation can mitigate these effects:Cavities in foreground objects can reliably be removed using amodified variant of [29] operating on the background withoutlabeling, but removing objects smaller than a specified size.The threshold can typically be chosen very liberally, e.g., asa percentage of the volume diagonal. Removal of loops onthe surface is more difficult, but in our experience a (binary)median filter performs well. The filter size should be chosento not disturb the smallest vessels in the image, but still beable to cover surface loops in the image.VII. C
ONCLUSION
Analyzing ultra large datasets by scalable algorithms is stilla huge challenge. We have presented a pipeline for extractingthe topology and various features of vessel networks fromlarge three-dimensional images. We were able to show thatour method fulfills all requirements identified in Section IV.Our main contribution, a novel iterative refinement approachand careful engineering of all pipeline stages, allowed usto demonstrate the scalability and thus applicability to verylarge datasets, e.g., generated by modern microscopes (primaryrequirements). At the same time, our pipeline can be appliedto a wide range of problem domains due to its robustness,unbiased nature, and lack of assumptions about the topologyand morphology of the analyzed vessel systems (secondaryrequirements).We will now apply the proposed pipeline and continueprevious work [4] using previously unachieveable level ofdetail, which hopefully leads to new biomedical discoveries.In the future, we would like to explore how even morespecific image features (c.f. [25]) can be integrated into thepipeline without compromising its scalability. Additionally,the current version of the pipeline is entirely single-threaded.While it is possible to more efficiently use the availableresources of modern hardware by simultaneously processingmultiple datasets in separate processes, it would be desirableto accelerate a single run using multiple processor coresor even GPUs. However, this is a challenging task: Whileapproaches to parallel skeletonization algorithms exist, theypose problems with regard to guarantees about the medialityof the centerline and reproducibility in general. Furthermore,at least in the current formulation, the connected-component-analysis [29], which is used in variations in several placesin this paper is inherently sequential. Offloading work to theGPU requires even more attention to detail with regard tomemory management. Furthermore, we would like to exploreautomatic segmentation approaches for vessel structures inlarge volumetric datasets, with special focus on (irregular)lymphatic vessel systems, to facilitate the usability of thepresented pipeline even further. Property G length G roundnessMean G straightness N FNR N FPR
Dataset RefinementLymphatic 1 Iterative Refinement
No Refinement 0.450 0.500 0.578 0.0515 0.0505Lymphatic 2 Iterative Refinement
No Refinement 0.486 0.536 0.614 0.0550 0.0539Lymphatic 3 Iterative Refinement
No Refinement 0.487 0.516 0.612 0.0653 0.0632Synthetic 1 Iterative Refinement
TABLE I: Results of the robustness test [30] using 36 rotations around the z-axis for each respective dataset. G property denotesthe GERoMe index for the given property . N F NR and N F P R use the same perturbation procedure, but measure the geometricerror [32] of the calculated centerlines and therefore show the aggregated maximum value. Refined graphs were extracted usinga bulge size of 1.5 and iteration until reaching a fixed point. Relative score differences above 10% are highlighted as bold text.The presented pipeline is part of version 5.1 of Voreen [33]a widely-used ([34], [35]) open-source volume processingand rendering framework. All scripts used in the evaluationare publicly available online to facilitate reproducing thepresented results. A CKNOWLEDGMENT
This work was funded by the Deutsche Forschungsgemein-schaft (DFG) – CRC 1450 – 431460824.R
EFERENCES[1] R. H¨agerling et al. , “A novel multistep mechanism for initial lymphan-giogenesis in mouse embryos based on ultramicroscopy,”
The EMBOJournal , vol. 32, no. 5, pp. 629–644, 2013.[2] P. Blinder et al. , “The cortical angiome: an interconnected vascularnetwork with noncolumnar patterns of blood flow,”
Nature Neuroscience ,vol. 16, no. 7, p. 889, 2013.[3] K. Bentley et al. , “The role of differential VE-cadherin dynamics in cellrearrangement during angiogenesis,”
Nature Cell Biology , vol. 16, no. 4,p. 309, 2014.[4] R. H¨agerling et al. , “VIPAR, a quantitative approach to 3D histopathol-ogy applied to lymphatic malformations,”
JCI Insight , vol. 2, no. 16,2017.[5] A. Nazir et al. , “OFF-eNET: An optimally fused fully end-to-endnetwork for automatic dense volumetric 3D intracranial blood vesselssegmentation,”
IEEE Transactions on Image Processing , vol. 29, pp.7192–7202, 2020.[6] Y. Ma et al. , “ROSE: A retinal OCT-angiography vessel segmentationdataset and new model,”
IEEE Transactions on Medical Imaging , 2021(accepted).[7] P. Campinho et al. , “Three-dimensional microscopy and image analysismethodology for mapping and quantification of nuclear positions intissues with approximate cylindrical geometry,”
Phil. Trans. RoyalSociety B , vol. 373, no. 20170332, 2018.[8] N. Hamilton, “Quantification and its applications in fluorescent mi-croscopy imaging,”
Traffic , vol. 10, no. 8, pp. 951–961, 2009.[9] Y. Zhao et al. , “Retinal vascular network topology reconstruction andartery/vein classification via dominant set clustering,”
IEEE Trans. Med.Imaging , 2019.[10] D. Lesage et al. , “A review of 3d vessel lumen segmentation techniques:Models, features and extraction schemes,”
Med. Image Analysis , vol. 13,no. 6, pp. 819–845, 2009.[11] C. Hu et al. , “Cerebral vessels segmentation for light-sheet microscopyimage using convolutional neural networks,” in
Medical Imaging 2017:Biomedical Applications in Molecular, Structural, and Functional Imag-ing, Orlando, Florida, United States, 11-16 February 2017 , 2017, p.101370K. https://zivgitlab.uni-muenster.de/d dree02/graph extraction evaluation [12] S. Moccia et al. , “Blood vessel segmentation algorithms - reviewof methods, datasets and evaluation metrics,” Computer Methods andPrograms in Biomedicine , vol. 158, pp. 71–91, 2018.[13] Y. Chen et al. , “Generation of a graph representation from three-dimensional skeletons of the liver vasculature,” in
Int. Conf. on BioMed-ical Engineering and Informatics, Proc. , 2009, pp. 1–5.[14] T. Lee et al. , “Building skeleton models via 3-D medial surface/axisthinning algorithms,”
Graphical Model and Image Processing , vol. 56,no. 6, pp. 462–478, 1994.[15] W. Choi et al. , “Extraction of the euclidean skeleton based on aconnectivity criterion,”
Pattern Recog. , vol. 36, no. 3, pp. 721–729, 2003.[16] M. Ilg and R. Ogniewicz, “The application of Voronoi skeletons toperceptual grouping in line images,” in
Conf. on Speech and SignalAnalysis, Proc.
IEEE, 1992, pp. 382–385.[17] K. Drechsler and C. O. Laura, “Hierarchical decomposition of vesselskeletons for graph creation and feature extraction,” in
Conf. on Bioin-formatics and Biomedicine, Proc.
IEEE, 2010, pp. 456–461.[18] Y. Chen et al. , “A thinning-based liver vessel skeletonization method,”in
Conf. on Internet Comp. and Inf. S., Proc.
IEEE, 2011, pp. 152–155.[19] P. Chothani et al. , “Automated tracing of neurites from light microscopystacks of images,”
Neuroinformatics , vol. 9, no. 2-3, pp. 263–278, 2011.[20] Y. Zhou and A. W. Toga, “Efficient skeletonization of volumetricobjects,”
IEEE Trans. Vis. & C.G. , vol. 5, no. 3, pp. 196–209, 1999.[21] S. Almasi et al. , “A novel method for identifying a graph-based repre-sentation of 3-D microvascular networks from fluorescence microscopyimage stacks,”
Med. Image Anal. , vol. 20, no. 1, pp. 208–223, 2015.[22] X. Cheng et al. , “Detecting branching nodes of multiply connected 3Dstructures,” in
Int. Sym. on Math. Morphology and Appl. to Signal andImage Processing, Proc.
Springer, 2019, pp. 441–455.[23] M. Couprie et al. , “Discrete bisector function and euclidean skeleton in2D and 3D,”
Img. Vis. Comput. , vol. 25, no. 10, pp. 1543–1556, 2007.[24] D. Babin et al. , “Skeletonization method for vessel delineation ofarteriovenous malformation,”
Comp. in Bio. and Med. , vol. 93, pp. 93–105, 2018.[25] A. Rodriguez et al. , “Rayburst sampling, an algorithm for automatedthree-dimensional shape analysis from laser scanning microscopy im-ages,”
Nature Protocols , vol. 1, no. 4, p. 2152, 2006.[26] R. Bates et al. , “Extracting 3D vascular structures from mi-croscopy images using convolutional recurrent networks,”
CoRR , vol.abs/1705.09597, 2017.[27] E. Grochowski, “Emerging trends in data storage on magnetic hard diskdrives,”
Datatech , pp. 11–16, 1998.[28] J. McCallum, “Historical cost of computer memory and storage,” https://jcmit.net/mem2015.htm, accessed: 2019-08-06.[29] M. Isenburg and J. Shewchuk, “Streaming connected component com-putation for trillion voxel images,” in
Works. on Mas. Data Alg. , 2009.[30] D. Drees et al. , “GERoMe – a method for evaluating stability of graphextraction algorithms without ground truth,”
IEEE Access , vol. 7, pp.21 744–21 755, 2019.[31] G. Hamarneh and P. Jassi, “Vascusynth: Simulating vascular trees forgenerating volumetric image data with ground truth segmentation andtree analysis,”
Computerized Med. Imaging and Graphics , vol. 34, no. 8,pp. 605–616, 2010. [32] D. Mayerich et al. , “NetMets: software for quantifying and visualiz-ing errors in biological network segmentation,” BMC Bioinformatics ,vol. 13, no. S-8, p. S7, 2012.[33] J. Meyer-Spradow et al. , “Voreen: A rapid-prototyping environment forray-casting-based volume visualizations,”
IEEE Computer Graphics andApplications , vol. 29, no. 6, pp. 6–13, 2009.[34] R. Landell et al. , “Material flow during friction hydro-pillar processing,”
Science and Technology of Welding and Joining , pp. 1–7, 2019.[35] F. Benz et al. , “Low wnt/ β -catenin signaling determines leaky vesselsin the subfornical organ and affects water homeostasis in mice,” eLifeeLife