Vessel geometry can be represented using implicit (voxel-based) or explicit (mesh-based) data structures. In the implicit representation, the vascular network is segmented from surrounding tissue using a three-dimensional voxel grid similar to the original image, with individual voxels labeled as inside or outside vessels. Voxels are convenient for calculating volumes or performing voxel-wise statistical analyses.

Explicit representations, such as polygonal meshes, define vessel surfaces through interconnected vertices and edges. Meshes facilitate calculations of surface area and enable simulations requiring surface geometry. In this survey, all segmentation algorithms produce voxel-based results. We tested two skeletonization algorithms (tagliasacchi and antiga) that require triangular meshes as input, which we produce from segmentation results using the marching cubes algorithm.

This implicit representation can optionally be converted to a surface mesh using algorithms such as marching cubes (Lorensen and Cline, 1998). Measurements are taken across either structure as convenient. For example, surface area can be measured by integrating across a surface mesh, while volume can be measured by adding up the pixels within the network.

The vascular skeleton is almost always represented explicitly for analysis using connected curves. This explicit representation is fundamentally a connectivity graph, where each node represents a bifurcation and each edge represents a single non-branching vessel segment (Figure 1). The vessel centerlines are curves that can be integrated to calculate features such as length and tortuosity. Algorithms such as depth- and breadth-first searches can also be applied to calculate path lengths and branching characteristics.

In this paper, we evaluate the performance of algorithms for extracting the geometry and connectivity of a microvascular network in large images on imaging methods applicable to whole organs. Most algorithms first binarize the original image using semantic segmentation, using the result as the foundation for medial axis transforms that provide the skeleton (Figure 3).

Meeting the criteria for resolution and volume coverage is challenging because microscopes are diffraction-limited and tissue samples are highly scattering. However, recent advances are starting to enable sufficient resolution and volume coverage. This paper considers the following broadly-accessible imaging methods:

X-ray microtomography (µ-CT) is nondestructive and measures the absorbance of X-rays incident on a sample to create three-dimensional images (Ritman, 2004). While µ-CT enables large-volume imaging of whole mouse brains (Hong et al., 2020), its low contrast limits spatial resolution. Recent advances in vascular perfusion compounds such as Vascupaint 2 (MediLumine, Montreal, Quebec, Canada) improve µ-CT resolution to ≈ 20 ± 4.0 µm, whereas previous contrast agents limited features to ≈ 92 ± 25 µm (Margolis et al., 2024).

Light sheet fluorescence microscopy (LSFM) is characterized by separating illumination and detection. A thin sheet of light illuminates the sample (Hsu et al., 2022), and an orthogonally oriented objective (Figure 4) collects the emitted two-dimensional image using a CCD or CMOS camera. While the penetration depth is traditionally limited by tissue scattering, recent developments in clearing protocols (ex. CUBIC or iDISCO+) (Susaki et al., 2015; Renier et al., 2014) enable large-scale imaging of whole rodent brains. Recent cleared-tissue implementations span a broad range of resolutions. For example, hybrid open-top light sheet systems can achieve lateral resolution of 0.45 µm and axial resolution of 2.9 µm across millimeter-scale volumes (Glaser et al., 2022). When combined with 4 × expansion microscopy, LSFM has reached effective resolutions of 375 nm laterally and 750 nm axially in centimeter-scale samples (Glaser et al., 2025). More recently, axially swept LSFM designs have demonstrated nearly isotropic resolution of approximately 300 nm in fixed and cleared tissues (Lin et al., 2025). These gains in resolution typically trade off against imaging volume coverage or acquisition speed, so whole-organ datasets still exhibit significant anisotropy.

Milling microscopy removes layers of a sample during the imaging process to expose deeper tissue volumes. Initial experiments used on two-photon microscopy followed by photo-ablation (Tsai et al., 2003), demonstrating that volume constraints could be eliminated by systematically removing tissue. More recent techniques separate tissue sections using physical cutting. Knife-edge scanning microscopy (KESM) (Mayerich et al., 2008; Mayerich et al., 2011) separates a tissue slice from the rest of the block during imaging, while milling with ultraviolet surface excitation (MUSE) performs block-face imaging followed by ablation (Guo et al., 2019).

We evaluated our binarization and skeletonization methods on three datasets representing the modalities described in Section 3.1.

X-ray microtomography (µ-CT) scans of mouse brain (Tg(Slco1c1-BAC-CreER); R26 ^-lsl-TdTom/+ ) vascular networks were acquired using a Skyscan 1276 (Bruker, Billerica, MA, United States) at an isotropic sampling rate of 10 µm per voxel. Mice were prepared based on previously published protocols (Suarez et al., 2024). Briefly, the vasculature was perfused with Vascupaint (MediLumine, Product number: MDL-121; Montreal, Quebec, Canada) to provide vascular and microvascular X-ray contrast (Figure 5a).

Light sheet fluorescence microscopy (LSFM) images of mouse brain microvasculature were acquired using a Cleared Tissue LightSheet (CTLS) Microscopy Workstation XL (3i, Intelligent Imaging Innovations, Denver, CO), equipped with a 30 fps Fusion BT sCMOS camera ( 2304 × 2304 resolution, 6.5 µm pixel size, 95% quantum efficiency, cooled to −50 °C; Hamamatsu Photonics, Japan). Adult C57BL/6N female mice (4.5 months old) were perfused with fluorescently labeled lectin (488 nm; Vector Laboratories, DL-1174-1) as previously described (Ahn et al., 2024), and brains were optically cleared using the iDISCO protocol (Hsu et al., 2022). Imaging was performed using a 1X (0.25 NA) objective and a 488 nm laser at 200 mW power with a 300 ms exposure time. Image stacks were acquired at a 6 µm step size in the Z-direction, with 15% right-side overlap and 50% center tile overlap, making the voxel size 2.0 × 2.0 × 6.0 µm. Images were stitched and reconstructed using SlideBook^™ software (Intelligent Imaging Innovations) with the LightSheet module for 3D multipoint acquisition. Raw image files were acquired in 16-bit TIFF format and rescaled to 8-bit using ImageJ for downstream processing. Final voxel spacing was adjusted to 1.0 × 1.0 × 0.75 µm using linear interpolation.

Milling microscopy images were acquired using a knife edge scanning microscope KESM at a voxel resolution of 0.6 × 0.7 × 1 µm, covering a 0.6 × 0.6 × 2 mm volume. The tissue was acquired from a normal mouse (C57BL/6J) perfused with India ink (Mayerich et al., 2011). The entire dataset is available using the KESM Mouse Brain Atlas (Chung et al., 2011).

Ground truth volumes for training and validation were manually annotated in Slicer (Kikinis et al., 2013). Binarization and skeletonization methods were evaluated on a 200 × 200 × 200 voxel dataset. The machine learning models were trained on six separate 128 × 128 × 128 voxel sub-volumes for each dataset (distinct from the 20 0 3 sample used for evaluation).

Each of these imaging modalities introduces sources of noise and other artifacts that challenge segmentation, including:

Resolution limitations are introduced due to both the signal strength and diffraction limit. Standard LSFM is typically anisotropic, with poorer axial than lateral resolution. Recent implementations (Glaser et al., 2022; Lin et al., 2025) can achieve sub-micron isotropic resolution, with trade-offs in speed or volume coverage. µ-CT is primarily limited by SNR to ≈ 20 µm in all three dimensions. Both of these constraints are larger than the diameter of the smallest microvessels (Figure 5a).

Staining and labeling used in LSFM relies on targeting cells in the vessel wall, producing hollow vessels when their diameter is significantly larger than the diffraction limit (Figure 5b). In practice, vascular perfusion with fluorescent compounds like dextran can overcome these artifacts while improving contrast (at the expense of molecular specificity). In that case we would expect images that are more comparable to KESM.

Blurred vessels can occur in LSFM images due to temporary misalignment during long imaging periods (Figure 5c).

Contrast variations are introduced by non-uniform illumination and/or staining in both LSFM and KESM. Contrast is also reduced in LSFM as a function of depth due to scattering.

Machining artifacts introduced by cutting tools in milling microscopy can cause lines (Figure 5e) and streaks (Figure 5f) in individual z-axis slices.

We evaluate a subset of available algorithms based on multiple criteria, prioritizing algorithms used in popular, domain-specific vessel analysis software. This includes Slicer (Kikinis et al., 2013), the Vascular Modeling Toolkit (VMTK) (Izzo et al., 2018), and VesselVio (Bumgarner and Nelson, 2022). Second, we prioritized algorithms with an established open-source implementation that was preferably provided by the authors. We selected representative algorithms from three classes of methods, including: (1) classical thresholding methods, (2) Hessian-based and gradient-based enhancement filters, and (3) deep learning-based semantic segmentation models. We selected skeletonization methods for a range of input data types, including: (1) binarized images, (2) meshes, and (3) point clouds. Our goal is to reflect diverse methodological approaches used for microvascular modeling. This selection process ensures that our evaluation adequately represents current best practices and popular trends in the field, providing a robust benchmark for future developments. Table 1 provides the dates and citations of the proposed algorithms.

Otsu’s method was evaluated as a baseline binarization algorithm (otsu3d), and used as the final binarization step as required for other algorithms. Vesselness filters, including both Frangi (frangi) and Beyond Frangi (bfrangi), were included due to their overwhelming popularity for vessel enhancement. Optimally oriented flux (oof) filters were included as a more recent innovation with the open-source implementations provided by the authors. U-Net architectures are extensively used for semantic segmentation in biomedical imaging, so we elected to test a baseline U-Net architecture (unet) while including recent work on self-configured U-Nets (nnunet).

We selected the thinning algorithms by Lee (lee) and Palágyi (palagyi) because of their popularity in skeletonization literature. While the level-set method by Kline (kline) and the confidence accumulation method by Kerautret (kerautret) are not used in existing modeling packages, the authors provide open-source implementations that were tested. We also selected the most recent skeletonization methods available for meshes (tagliasacchi). The selected methods and keywords are provided in Table 2.

Since microvasculature accounts for ≈ 4% of the total tissue volume, we focus on metrics that can accurately characterize unbalanced classifications. The Jaccard similarity index (J) and Dice similarity coefficient (D) are the most common for binarization. The Jaccard index is the normalized volume overlap between the binarization and ground truth: J A , B = ∑ A ∩ B ∑ A ∪ B The Jaccard index is in the range J ∈ [ 0,1 ] , where J = 0 indicates no overlap and J = 1 indicates that the binarization is identical to the ground truth.

We also calculate the precision, emphasizing the accuracy of the positive predictions, and recall, which quantifies the model’s ability to identify true positives. The precision (p) and recall (r) are often combined using the F-score (F), which is both equal to the Dice coefficient and a function of the Jaccard index in the case of binarization: F = 2 p r p + r = D = 2 J 1 + J

We adopt the definition of the skeleton as: a set of curves that have identical topology to the geometric surface and are equidistant to the boundary (Wei et al., 2018). We enforce these conditions by manually generating ground truth skeletons based on this definition.

Skeletonization is evaluated using NetMets (Mayerich et al., 2012), which calculates the similarity between two sets of curves A and B in 3D space: M A , B = 1 L ∑ a ∈ A ∫ e d a t , B 2 2 σ 2 d t where L is the total length of all curves in the set A , a ( t ) is a point on the i th curve parameterized by t , and σ is a sensitivity parameter. The function d ( x , B ) is the distance between a point x and the closest point in the set of curves B . Given two networks representing the set of curves in the ground truth G and test case T , the precision (positive predictive value) P and recall (true positive rate) R are calculated: P = M T , G R = M G , T The precision is the percentage of the test skeleton that correctly corresponds to the ground truth, and the recall is the percentage of the ground truth skeleton that is correctly detected in the test case. The sensitivity parameter σ = 1 µm is used for all images.

The most popular method for enhancing vessels is the Hessian-based approach described by Frangi et al. (1998), originally designed for 2D retinal fundus images. The Hessian matrix is calculated at each point in the image x ∈ R 3 using finite differences across S scales, creating a field H ( x , s ) ∈ R 3 × 3 . The filter response V ( x ) ∈ R is calculated using the Hessian matrix eigenvalues ( | λ 1 | ≤ | λ 2 | ≤ | λ 3 | ) (Supplementary Algorithm S1). In this study, the method introduced by Yang and Cheng (2014) was also used to accelerate the computations of the Hessian matrix by a factor of two.

The original “vesselness” algorithm (frangi) is outlined in Algorithm 1 and relies on four parameters: tuning parameters α , β , c that target the vessel shape, and a scale parameter γ ≈ 2 . The shape parameters are balanced to separate tube-like structures and bifurcations from background pixels (Figure 6). The scalar γ was proposed earlier (Lindeberg, 1998) to tune the derivatives used in the Hessian. Attempts to optimize this parameter yield a consistent value of γ = 2 in the literature. We found that this is due to the scale factor σ γ compensating for energy dissipation from the second-order Gaussian scale-space filter. A value of γ = 2 ensures that the most intense response comes from features near the scale-space parameter σ (Figure 7).

Since the parameters in the original algorithm can be challenging to select, we provide two comparisons. First, we examine binarization results based on parameters used in the original paper (frangi), and after parameter optimization (ofrangi) for each dataset using training data. A parallel implementation was used to create sensitivity maps for each parameter (Figure 8). The relationship between parameters for frangi is shown in Figure 6a.

A modified implementation of the vesselness framework replaces the parameters with a single normalization term τ (Jerman et al., 2015). The Beyond Frangi (bfrangi) algorithm (Supplementary Algorithm S2) regularizes the largest eigenvalue at each scale, yielding a more uniform response across vessel cross-sections and intensities. The sensitivity of each dataset to this parameter is illustrated in Figure 9.

The second derivatives within the Hessian are numerically calculated using finite-difference methods, introducing artifacts when features are adjacent. An alternative approach proposes optimally oriented flux (OOF) (Law and Chung, 2010), which relies on first derivatives and integrating gradient information over a spherical region. This applies a hard cut-off that localizes the gradient information used.

Given a point x in the 3D image I ( x ) ∈ R and a direction vector ρ , the gradient flux aligned with ρ is integrated over a sphere Ω with radius σ s (Figure 10). F x , r , ρ = ∬ Ω ∇ I x ⋅ ρ ρ ⋅ n d A = ρ T Q x , r ρ where d A is the differential area of the sphere and n is the sphere normal at x . This integral is used to calculate the flux matrix Q . The optimal vessel orientation is the vector ρ that maximizes the outward-oriented flux: max ρ ρ T Q x , r ρ Our tests used the OOF method available in MATLAB (oof), where the only tunable parameters are the sphere radii Σ . In our tests, we use the same Σ values for OOF and scale-space arguments (frangi, ofrangi, bfrangi).

The U-Net (unet) architecture uses parallel encoder/decoder paths bridged by skip connections that merge high-resolution features from the encoder into the decoder. The network meta-parameters were determined through iterative trial and error based on validation performance. This U-Net uses 3 × 3 × 3 convolutional layers, ReLU activation, and 2 × 2 × 2 max-pooling layers (Figure 10). The data is divided into smaller volumes, selecting only those with intensities above the global mean to avoid empty or near-empty voxels that reduce training quality (Çiçek et al., 2016). The unet model tested here is trained and tested on these selected volumes using the Focal Tversky loss function (Abraham and Khan, 2019).

The nnU-Net (nnunet) architecture builds on U-Net by overcoming the need for manual meta-parameter tuning. nnU-Net maintains the encoder-decoder structure with skip connections while automatically adapting its architecture and hyperparameters. Like U-Net, it employs 3 × 3 × 3 convolutional layers, ReLU activation, and pooling layers but integrates additional strategies, such as ensemble modeling and deep learning-based boundary enhancement to improve segmentation accuracy (Schott et al., 2023). During training, nnU-Net uses standard loss functions and ensemble modeling, where multiple models are trained and combined. We used the source code provided by the authors ¹ for these experiments.

Our binarization results are shown in Table 3, providing the Jaccard Index, precision, and recall. Otsu’s method (otsu) was tested on raw data and used to select the final threshold for all preprocessing algorithms. The vesselness filter was tested using published parameters in the original paper (frangi) and after optimization of all parameters (ofrangi), while Beyond Frangi (bfrangi) was tested after optimization of τ . Optimally oriented flux (oof) does not require parameter tuning outside of the scale-space parameters used to select sphere radii. Precision-recall curves are provided for the raw output so that the reader can evaluate the results without Otsu’s method (Figure 11). Like most semantic segmentation networks, unet and nnunet integrate binarization into the model and do not require Otsu’s method for thresholding.

Raw image sections are provided from each dataset, along with the best corresponding binarization result for each algorithm (Figure 12). Further discussion is provided in Section 7.

Skeletonization is the method of determining the connected vessel centerlines, usually based on the extracted geometry. A binarized image is used as an input to thinning algorithms, while other approaches calculate the skeleton (or medial axis) from an explicit geometry.

Thinning uses morphological erosion to calculate an implicit representation of the skeleton from a binarized image. Voxels are iteratively removed until a thin single-voxel skeleton remains. The concept was originally proposed in Blum’s grassfire propagation (Blum, 1967), while Tsao and Fu (1981) were the first to introduce a three-dimensional algorithm based on formalized topological and geometric rules. The main concept in erosion-based methods is determining the “simple points” that can be removed without changing the model topology. Thinning algorithms usually differ in how they detect and remove simple points. The most popular method used in current software packages (Fedorov et al., 2012; Bumgarner and Nelson, 2022) is Lee et al. (1994) approach.

Mesh-based medial axis transforms approximate the skeleton from a geometric surface. These include point-sampling methods that leverage ray casting Kerautret et al. (2016) to localize the skeleton. While these methods face challenges when differentiating the “inside” and “outside” of complex shapes (Wei et al., 2018), statistical methods can overcome these problems using vector accumulation maps to probabilistically approximate the skeleton (Kerautret et al., 2016). This can be overcome by directly evolving the geometry towards its centerline using an active contour model (Tagliasacchi et al., 2012).

Another subcategory of mesh-based transforms leverages a Voronoi diagram of the surface (Antiga et al., 2003). The boundaries of the Voronoi cells within vessels can be tracked to extract the medial axis. Despite their precise results, Voronoi-based techniques require careful handling of boundary conditions and specifying seed points to avoid generating spurious branches.

Thinning algorithms reduce an object to a single-voxel-thick centerline while maintaining topology (ex. connectivity, cavities, tunnels). This is accomplished by examining the local neighborhood of each voxel to identify simple points whose deletion preserves topology. These points are iteratively removed until a thinned skeleton remains.

The lee algorithm (Lee et al., 1994) identifies simple points using an octree to recursively subdivide the local neighborhood (Supplementary Algorithm S3). A later approach (Palágyi and Kuba, 1999) uses predefined templates to detect and flag simple points, which is easier to parallelize. Both the sequential lee and parallel palagyi algorithms are parameter-free and take a binarized image as input.

The Kerautret et al. (2015) algorithm extracts a centerline from a surface mesh by casting rays from each surface point along its normal. The rays are accumulated in a voxel grid where the maximal ridges represent centerlines. The most recent modification (Kerautret et al., 2016) determines an accumulation confidence value for each voxel, where points with higher confidence are more likely to be located on the centerline.

The only required parameter is the accumulation distance ( d acc ) or ray length, which is slightly larger than the maximal radius of the input shape. This value is straightforward to approximate for each volume using the radius of the larger vessels.

The kline algorithm (Kline et al., 2010) calculates a distance field using the fast marching method (FMM), and then internal ridges (maxima) are followed from user-specified seed points. This algorithm assumes that the geometry is a connected tree of tubular structures branching from a single root point. Since microvascular networks do not have a unique root, a large number of seed points are required to comprehensively trace the network.

The tagliasacchi algorithm (Tagliasacchi et al., 2012) calculates the skeleton from a closed mesh using mean curvature flow (MCF). The medial axis is extracted by moving the mesh along its surface normal proportional to the surface curvature. This iteratively collapses the mesh into a one-dimensional curve (Figure 13). We tested the tagliasacchi method using an implementation provided by StarLab ² , which requires a watertight manifold mesh and is highly sensitive to mesh boundaries. We forced a watertight mesh at all boundaries and applied the algorithm to each closed mesh component separately.

The antiga algorithm (Antiga et al., 2003) builds a Voronoi diagram from point samples of the vascular surface. The intersections of each Voronoi region lie equidistant to the sampled surface, such that medial axis curves can be traced between user-specified seed points. The Delaunay tessellation is calculated from a set of points P that densely samples the vascular surface, which provides the vertices V for the Voronoi diagram. The user specifies a set of seed points S that are connected by selecting points in V that maximize the distance from the vascular surface. This is done by solving an Eikonal equation using the fast sweeping method, ensuring a maximum distance from the boundary surface and providing stability in complex regions like bifurcations. This study uses the implementation of antiga provided in Slicer ³ . However, the method failed to process seed points or centerlines for messy surfaces created in some LSFM datasets (Table 4).

The precision and recall for the proposed skeletonization methods are shown in Table 4. The NetMets algorithm is used to quantify performance using σ = 1 µm. The results are provided using all binarization algorithms as input. The ground truth binarization is also provided as a baseline to demonstrate skeletonization from an “ideal” starting point (Table 5).

The lee and palagyi methods achieved the highest F-score across all images, capturing vessel continuity at the cost of adding spurious branches to larger vessels. The antiga method provided the highest precision, however, it identified fewer vessels. The method failed to extract skeletons from binarized results using Otsu’s method due to the excessive number of vertices and edges. Gradient-based methods (kerautret and kline) performed poorly across most images. While kline was designed for vessels, it is poorly suited to microvascular networks since they are not tree-like. The difference in the performance of the kerautret for KESM data also stands out, likely due to the higher resolution of KESM data. This provides (1) a more completely connected mesh and (2) larger distances between adjacent unconnected vessels.