Robust adaptive 3-D segmentation of vessel laminae from fluorescence confocal microscope images and parallel GPU implementation

Abstract

This paper presents robust 3-D algorithms to segment vasculature that is imaged by labeling laminae, rather than the lumenal volume. The signal is weak, sparse, noisy, nonuniform, low-contrast, and exhibits gaps and spectral artifacts, so adaptive thresholding and Hessian filtering based methods are not effective. The structure deviates from a tubular geometry, so tracing algorithms are not effective. We propose a four step approach. The first step detects candidate voxels using a robust hypothesis test based on a model that assumes Poisson noise and locally planar geometry. The second step performs an adaptive region growth to extract weakly labeled and fine vessels while rejecting spectral artifacts. To enable interactive visualization and estimation of features such as statistical confidence, local curvature, local thickness, and local normal, we perform the third step. In the third step, we construct an accurate mesh representation using marching tetrahedra, volume-preserving smoothing, and adaptive decimation algorithms. To enable topological analysis and efficient validation, we describe a method to estimate vessel centerlines using a ray casting and vote accumulation algorithm which forms the final step of our algorithm. Our algorithm lends itself to parallel processing, and yielded an 8 x speedup on a graphics processor (GPU). On synthetic data, our meshes had average error per face (EPF) values of (0.1-1.6) voxels per mesh face for peak signal-to-noise ratios from (110-28 dB). Separately, the error from decimating the mesh to less than 1% of its original size, the EPF was less than 1 voxel/face. When validated on real datasets, the average recall and precision values were found to be 94.66% and 94.84%, respectively.

Fig. 1

Gallery of sample images illustrating some of the challenges to surface segmentation. (a) An open vessel end with large variations in vessel diameters. (b) Closely situated vessel segments. (c) Abrupt variations in vessel signal. (d) Irregular broken vessel segments. (e), (g) Large signal variations and gaps. (f) Shows vessel with strong depth-dependent attenuation. (h) Low-contrast and blurry vessel signal with complex branching. (i) Contrast stretched image showing significant signal overlap with other fluorescent imaging channels (nuclear channel, in this case). Our segmentation algorithm has been designed to robustly detect the vasculature overcoming these issues.

Fig. 2

Illustrating the performance of robust hypothesis testing based initial foreground detection compared to other approaches. (a) In the alternate hypothesis, there is a locally planar patch of the vessel surface passing through x within the rectangular neighborhood Γ(x). In the null case, all of the voxels in the neighborhood belong to the background. (b) maximum-intensity projection of the original 3-D confocal dataset that is contrast stretched to better show the imaging noise and overlap from the nuclear channel (small circular objects). (c) Adaptive Otsu thresholding. (d) Hessian filtering followed by adaptive Otsu thresholding. (e) Direct tube tracing results [4]. (f) Basic surface extraction algorithm described in Section II-B. (g) Results of adaptive region growth (Section II-D). These two panels have the same sensitivity (true positive rate). (h) The variation of the false alarm rate (with median estimate) for various background intensity values (5×5×5 neighborhood), for the case when λ_T = 3.

Fig. 3

Illustrating the automated vessel segmentation results for a confocal image from the cortical region of the rat brain (field size: 368.55 μm × 368.55 μm × 93.00 μm with field resolution 0.36 μm × 0.36 μm × 1.5 μm giving rise to 1024 × 1024 × 63 voxels). (a) Maximum-intensity projection of the grayscale image data displayed in shades of red. (b) Automated segmentation results rendered as a surface (white) viewed axially. The grayscale data is overlaid on the rendering to enable simultaneous viewing of the raw data and the segmentation. (c) Display of the segmentation confidence estimates over the 3-D field using the indicated color map. (d) Display of the local curvature estimates using the indicated color map. (e) Result of automatic vessel centerline extraction. (f) Illustrating manual editing of centerlines by an expert observer for edit-based validation.

Fig. 4

Illustrating the voting based centerline extraction method. (a) maximum-intensity projection of a confocal dataset. (b) Small vessel segment from panel A shown enlarged and volume rendered. (c) Voting results overlaid on the surface rendering, using the same color map as in Fig. 3. Small sliced region shows the cross-section. (d) Extracted centerline shown in red along with the surface rendering.

Fig. 5

Illustrating the parallel implementation in GPU. (a) Table of runtimes of the CPU code and the GPU code. (b) Illustration of the multipass bitonic sorting used in the GPU implementation. (c) Flowchart of the surface segmentation algorithm with adaptive region growing step. We implemented five different kernel functions to achieve a parallel implementation in the GPU.

Fig. 6

Showing the 3-D vessel phantom (512 × 512 × 100 × 8-bit) used for our experiments. (a) Maximum-intensity projection displayed using the colormap (grayscale units) shown on the left. The phantom was designed to simulate varying diameters, intensity values, and the presence of spectral unmixing artifacts (faint blobs in the background). (b) 3-D surface rendering of the phantom’s vessel structure. (c) Slice 19 of the phantom displayed using the colorscale on the left to show the hollow nature of the vessels. (d) Slice 19 data shown with added noise (Gaussian, σ = 90).

Fig. 7

Results on phantom image data. (a) Shows a plot of the average error per face (EPF) measure as a function of the peak signal-to-noise ratio (PSNR) for the case of Poisson and additive white Gaussian noise respectively for three algorithms. (b) Compares the best performance of basic algorithm with/without region growth in the presence of Gaussian noise for various values of the noise standard deviation (σ). The adaptive region growing step clearly picks up more vessel content without picking up the noise. (c) Compares the performance of the basic algorithm over the Hessian filtering based algorithm. (d) Shows one single slice of the noisy image data for Gaussian noise with σ = 90 along with the segmentation results from the Hessian filtering, basic algorithm, and adaptive region growth, respectively.

Fig. 8

Segmentation results on real data. (a)–(c) Volume rendering of the image data. (d)–(f) 3-D rendering of the mesh generated from the segmentation. (g)–(i) Confidence map of the segmentation is mapped over the segmentation (color scale shown in Fig. 6). (j)–(l) Results of edit-based validation (green lines are computed centerlines, blue lines were added by the human expert, and red lines were deleted).