Keywords

1 Introduction

Magnetic Resonance Angiography (MRA) methods, like Time-of-flight (TOF) allow visualization of intracranial arteries without radiation dose or contrast agents. Beyond the clinical use of MRA for stenosis identification, a whole brain vasculature map can be generated through digital reconstruction of intracranial arteries, such as artery tracing and segmentation. The vascular features extracted from the vasculature map include artery length, volume and tortuosity, which provides quantitative measurements for pathological or blood flow conditions [1, 2]. However, due to the complex network structure of human intracranial arteries, and the low signal intensity and low image contrast for distal branches from TOF MRA, it is challenging for automated quantifications of the vasculature map, especially for patients with compromised cerebral blood flow due to cardiovascular diseases.

Artery segmentation and tracing from MRA are two methods used for digital reconstruction of intracranial arteries in order to measure the vascular features.

Artery segmentation which classifies every voxel into artery or non-artery, allows better visualization of artery structures and facilitates identification of stenoses and aneurysms [3]. Existing automated MRA segmentation approaches include region growing [4], active contours [3], and convolutional neural network [5]. However, artery segmentation alone cannot determine the layout and inter-connected relation of intracranial arteries, which limits the information in the vasculature map.

One solution is 3D artery tracing, which converts artery voxels into interconnected tree structures with radius varying cylinders. The criteria for connecting neighboring points is critical, for which local Hessian-based estimation [6], and Kalman filtering [7] have been attempted. But the performance of the methods by tracing directly from the original image may suffer due to the low artery signal or contrast, which is commonly seen in TOF MRAs. To ensure robust tracing, refinement of the artery centerline by re-centering is usually followed by tracing. Adjusting centerline positions by applying intensity features [6] or segmentation results [8] from re-sliced 2D cross-sectional planes is usually used, but the neighboring slice information is not considered. The multiplanar reformation (MPR) of arteries, considered as straightening the artery along its centerline (example in Fig. 4(c)), has been used in clinical reading and reported to be beneficial for coronary and renal artery stenoses detection [9, 10]. MPR view incorporates neighboring information and therefore is better suited for a global centerline refinement on the whole artery and correction of errors in artery tracing.

This paper introduces a novel method by simultaneously performing artery tracing and segmentation with the help of robust artery refinement in MPR views. Tracing and segmentation use results from each other to improve their individual performances. Instead of tracing directly from original images, artery segmentation from a deep neural network model is used to enhance the contrast of small arteries and constrain artery tracing in a restricted region. The artery refinement from MPR ensures trace smoothness, radius fitness, and avoids centerline deviation. The refined traces can then be used to further improve the artery segmentation results, so that both centerline features (length, tortuosity) and voxel features (volume) of arteries are accurately extracted.

2 Method

The method has four steps: artery rudimentary segmentation, tracing from segmentation, MPR refinement and segmentation from tracing. Flow chart is shown in Fig. 1.

Fig. 1.
figure 1

Flow chart for the proposed method

Full size image

2.1 Rudimentary Segmentation

Original images are normalized using the Nyul method [11], and rudimentarily segmented using a neural network model [5]. Trained from the semi-automatedly segmented arteries from more than 163 subjects [2], the model can segment small arteries with low contrast to the background. Centers of each 3D connected region in the segmentation image are considered as seed points to initiate tracing one at a time.

2.2 Tracing from Segmentation

From each position, a cross-sectional plane can be generated from each of the evenly distributed positions around the half sphere towards the tracing direction with 30° increments in each axis (in total 37 possible directions). Arteries are traced by iteratively finding the best matching cross-sectional planes based on the matching criteria. Direction selection during one iteration as an example is shown in Fig. 2.

Fig. 2.
figure 2

Direction selection in one iteration of tracing. For illustration, the red and yellow directions are selected from the 37 directions on the half sphere in the tracing direction. Three metrics are calculated from the cross-sectional slices for deciding the best tracing direction. (Color figure online)

Full size image

Matching Criteria.

Three metrics are used to select the best matching cross-sectional planes: circular similarity, neighboring similarity and signal change rate.

Based on the assumption that the cross-sectional plane along the centerline of an artery is a circle [12], the circular similarity is calculated as the Dice similarity coefficient (DSC) [13] of the segmented region in the cross-sectional plane with a perfect circle of the same area. The circular similarity is 1 for perfect circle.

A smooth trace should be continuous between neighboring cross-sectional slices both in geometry and intensity. The DSC for neighboring segmented regions, and the relative signal difference in their center pixels from the original image are used for evaluating neighboring relations.

Matching score is defined as weighted sum of circular and neighboring similarities (weights of 0.8 and 0.2 used in this study) minus signal change rate. The highest score is used for selecting the best match. As an example, in Fig. 2, the yellow direction has higher score than the red one, so it is selected as the stretching direction.

Tracing Procedure.

Initiated from a seed point in the artery region as a starting point \( \varvec{p}_{0} = \left( {x_{p} ,y_{p} ,z_{p} } \right) \), circular similarity is used to find the best matching trace direction \( \varvec{n}_{i} \) from the cross-sectional plane \( C_{i} \left( {u,v} \right) \). Neighboring similarity is not combined, as there are no neighbor slices available in the first iteration. The positions of \( p_{ - 1} = p - n_{i} \) and \( p_{1} = p + n_{i} \) are added to the trace with the radius \( r_{0,1, - 1} \) calculated from the masked region in \( C_{i} \left( {u,v} \right) \).

Starting from the second iteration, the trace is stretched at both sides in each iteration until each reaches the ending criteria. Taking positive tracing side at iteration \( j \) as an example, search space for next trace position is updated with \( \varvec{t}_{\varvec{j}} = \varvec{p}_{j - 1} - \varvec{p}_{j - 2} \) and the normal direction along the tracing direction \( \varvec{n'}_{i} \) yielding the highest metric from the rotated cross-sectional plane \( C'_{i} \left( {u,v} \right) \) is used for stretching the traces to the new position \( p_{j} = p_{j - 1} + n'_{i} r_{j - 1} \). Radius of the target position \( r_{j} \) is roughly estimated from the masked area and used as the stride for the next iteration. Traced region in the segmentation image is painted with zeros and the seeds inside the region are removed to avoid repeat tracing.

The iteration on each side of stretching is ended when \( \varvec{p}_{\varvec{j}} \) is out of image boundary, or the maximum combined similarity from all \( C'_{i} \left( {u,v} \right) \) is below a certain threshold (0.1 in this study).

2.3 Artery Refinement from Multiplanar Reformation

The position and radius of the centerline based on cross-sectional plane matching does not incorporate global information, therefore the refinement step applies MPR to overcome centerline deviation and abrupt radius change along the centerline. A three-stage optimization scheme is used for the artery refinement, i.e., trace position refinement, trace radius refinement, and trace deviation correction.

Trace Position Refinement.

The 3D position of points in the trace is refined using the optimization function considering losses for trace smoothness and their intensity.

$$ L_{1} \left( \varvec{p} \right) = \sum\nolimits_{i} {L_{1} \left( {\varvec{p}_{i} } \right) = \sum\nolimits_{i} {\left\{ {w_{dist} \left( {\left| {\left| {\varvec{d}_{i} } \right|} \right|^{\left( 1 \right)} + \gamma \left| {\left| {\varvec{d}_{i} } \right|} \right|^{\left( 2 \right)} + \gamma \left( {d_{x,i}^{\left( 2 \right)} + d_{y,i}^{\left( 2 \right)} + d_{z,i}^{\left( 2 \right)} } \right)} \right) - w_{int} \left[ {I_{n} \left( {\varvec{p}_{i} } \right) + I_{s} \left( {\varvec{p}_{i} } \right)} \right]} \right\}} } $$
(1)

where \( \varvec{d}_{i} = \left( {\varvec{p}_{i} - \varvec{p}_{i - 1} } \right) = \left( {d_{x,i} ,d_{y,i} ,d_{z,i} } \right) \), \( I_{n} \left( \varvec{p} \right) \) and \( I_{s} \left( \varvec{p} \right) \) are intensity values of normalized (\( M_{n} \)) and segmented (\( M_{s} \)) MPR images at position of \( \varvec{p}_{i} = \left( {x_{i} ,y_{i} ,z_{i} } \right) \). (1) and (2) represent 1st and 2nd order of derivative. \( \gamma \) is the parameter to control the first and second order weights. \( w_{dist} \) and \( w_{int} \) are weights for controlling the smoothness and intensity loss.

Trace Radius Refinement.

After the trace position refinement, centerline positions are fixed, and the radius of each point is refined using the following equation.

$$ \begin{aligned} L_{2} \left( \varvec{r} \right) = & \sum\nolimits_{i} {L_{2} \left( {u_{l,i} ,u_{r,i} ,v_{i} } \right) = } \sum\nolimits_{i} {L_{2} \left( {l\left( {v_{i} } \right),r\left( {v_{i} } \right),v_{i} } \right) = } \\ & \sum\nolimits_{i} {w_{\text{smooth}} } \left[ {l^{\left( 1 \right)} \left( {v_{i} } \right) + r^{\left( 1 \right)} \left( {v_{i} } \right) + \gamma l^{\left( 2 \right)} \left( {v_{i} } \right) + \gamma r^{\left( 2 \right)} \left( {v_{i} } \right)} \right] \\ & - w_{\text{grad}} \left[ {M_{u} \left( {l\left( {v_{i} } \right),v_{i} } \right) + M_{u} \left( {r\left( {v_{i} } \right),v_{i} } \right)} \right] \\ \end{aligned} $$
(2)

where \( l\left( v \right) \) and \( r\left( v \right) \) are the left and right boundary for artery radius in MPR image \( M_{n} \). \( M_{u} \) is the derivative of \( M_{n} \) in its horizontal direction.

Trace Deviation Correction.

Ideally, the mean location of the left and right radius boundaries \( \frac{{l\left( {u,v} \right) + r\left( {u,v} \right)}}{2} \) in the MPR image should always be in the vertical center of the MPR image \( \left( {v = v_{m} } \right) \). Any deviation away from the vertical center in \( u \) direction \( \varvec{o = }\frac{l\left( v \right) + r\left( v \right)}{2} - v_{m} \) needs to be re-centered.

Iterative Refinement from Different Angles.

MPR images \( M_{deg} \) are reconstructed using rotation angles from {0, 90, 45, 135} by repeating Rep times (3 in this study). Arteries are iteratively refined every 25 iterations of tracing and at the end of tracing.

Nelder-Mead algorithm [14] is used for optimization in this study.

2.4 Segmentation from Tracing

From the refined traces, regions inside the tubes are filled to be the refined segmentation results. As the trace is represented by a cylinder model, if a more detailed artery area information is needed, cross-sectional planes can be generated along the centerline and the artery region can be segmented based on refined radius boundary.

3 Evaluations

3.1 Accuracy for Quantification of Vascular Features

Four images each from a healthy community study [15] and a group of patients with intracranial atherosclerosis [16] are used for evaluation. The data collections followed local institutional review board guidelines. Three-dimensional TOF images were scanned on 3.0T MR scanners with: repetition time/echo time = 25/3.5 ms, flip angle = 20°, in-plane resolution = 0.35 mm × 0.35 mm, slice thickness = 1.4 mm.

Considering the unrealistic work load for manually labeling voxels for all regions of intracranial arteries in 3D images, a semi-automated tool [17] is used for tracing artery regions with manual corrections. The centerlines are considered as ground truth for the evaluation. As the completeness of artery detection is not our focus, only two major clinically important arteries (from the distal internal carotid artery (ICA), M1 segment of middle cerebral artery, until the most distally clearly visualized segment) per case are used for validation. Excessive traces are removed, and traces might be reconnected at bifurcations to allow same branches being compared with the ground truth.

Vascular features of artery length and tortuosity (length divided by the Euclidean distance between the first and last point) are used to evaluate the tracing performance, assessed by the mean percentage of absolute difference with the ground truth. DSC is used to evaluate the segmentation performance (artery volume difference).

From Table 1, the mean absolute length and tortuosity difference is within 4 percent, and the mean Dice similarity score of volume is more than 0.8, indicating excellent agreement [18]. Performance in the disease group is slightly worse than the healthy group. An example of 3D rendered artery traces with ground truth is shown in Fig. 3. All the 3D renderings of arteries are in the supplementary materials.

Table 1. Extracted vascular features compared with ground truth.
Full size table
Fig. 3.
figure 3

(a) MPR view for ICA (b) Tracing results generated from the proposed method (ICA_L: green and R: red). (c) Ground truth for the traced arteries. (d) Tracing results without artery refinement in the ablation test (Color figure online)

Full size image

The performance of this method is compared with other three artery segmentation and tracing methods (colliding fronts, fast marching, iso-surface) implemented in open-source Vascular Modeling Toolkit (VMTK, www.vmtk.org). The comparison results for left side arteries from healthy subjects are shown in Table 2. All three comparison methods suffer from weak signals in distal arteries, so that the tracing is stopped earlier. Right side arteries are not processed due to the processing time. VMTK methods fail for the diseased group. One visualization figure is in supplementary materials.

Table 2. Performance comparison with other segmentation and tracing methods.
Full size table

3.2 Continuity of Tracing

To evaluate the tracing continuity, a semi-automated snake based method [19] is used to trace the ICA. Due to the tortuous ICA structure and flow artifacts reducing luminal contrast, the snake method needs an average of 4.9 manually given seeds (including the seed used in our tracing method) to trace the whole artery segment, but all arteries traced by our method need only one seed, showing better performance.

The ablation test of artery tracing without refinement showed worse performance in Table 1, and the tracing iterations are aborted earlier for 4 of the arteries.

3.3 Improvement for Segmentation

Confined in regions within traces, artery segmentation is further improved, especially when multiple arteries are close to each other. An example is shown in Fig. 4.

Fig. 4.
figure 4

Improved segmentation from tracing when two arteries are close. Cross-sectional slice for normalized image (a), rudimentary segmentation (b) and improved segmentation (c)

Full size image

4 Limitations

Several limitations exist in this method. A selected set of major arteries from each subject is assessed. A small sample size (N = 8) is used for validation. Parameter tunings are also expected to further improve the performance.

5 Conclusion

In this paper, a simultaneous artery tracing and segmentation method with artery refinement from MPR view is proposed. The use of segmentation results to trace arteries allows tracing continuity and trace accuracy. The use of tracing allows centerline features quantified from a specific artery. The use of MPR view for artery refinement improves robustness by correcting mistakes in artery tracing.