Continuous extraction of coronary artery centerline from cardiac CTA images using a regression-based method

1.   Introduction

Coronary artery disease (CAD) is one of the major leading causes of death all over the world ^[1]. Currently, the gold standard for evaluating CAD is coronary angiography, while it is limited to invasiveness, operational complexity, and high cost ^[2]. The three-dimensional (3D) model reconstructed from cardiac computed tomography angiography (CTA) can provide intuitionistic and accurate information on the geometry of the coronary artery, which crucially contributes to the non-invasive diagnosis of CAD ^[3,4]. There are already some visualized vessel structure techniques in clinical practice, such as curved planar Reformation (CPR) ^[5]. CPR can obtain straightened coronary artery images, helping planning intervention procedures and choosing stent size ^[6,7]. The 3D reconstruction of the coronary artery typically involves coronary centerline extraction and vessel lumen segmentation. Coronary centerline extraction as the first step can facilitate branch detection and lumen radius identification, which partly resists interference from nearby veins and myocardium ^[8]. This was in accordance with Shen et al.'s observation in which the model with centerline as prior knowledge could obtain a higher accuracy of vessel segmentation ^[9]. Whereas, manual extraction of coronary centerlines is a challenging task because the coronary arteries structure is complex and the annotation process is time-consuming.

In recent years, deep learning, especially U-Net and its variants, has been widely used in medical imaging segmentation tasks ^[10]. The U-Net and its variants have demonstrated the ability of blood vessel segmentation, which achieved promising results in medical image analysis ^[11,12]. Convolutional neural network (CNN) showed advantages in extracting high-level features of vascular structures ^[13]. Dorobanțiu et al. ^[14] adopted 3D U-Net architecture to segment the coronary artery centerline with voxel precision, estimating whether each voxel was the centerline point. Rjiba et al. ^[15] employed a local vessel filter based on a binary classification network, successfully extracting centerlines by determining whether the center point of the current patch was located on the centerline. However, these coronary artery structure extraction tasks still face numerous challenges. To start, the full CTA images are usually broken into small patches for training, which carry a great deal of non-vascular information. Accordingly, the process is time-consuming and tedious. Besides, these methods require complex post-processing techniques that can deal with the problem of extracted vessel discontinuities, preserving the completeness of the coronary centerline. Moreover, the coordinate of centerline point obtained by the binary classification network was imprecise, as there exist differences between the center point of patch and the centerline point.

Consequently, the path tracing method, a continuous process, has introduced great precision and efficiency to coronary artery segmentation. Wolterink et al. ^[16] proposed a CNN-based orientation classifier to determine the artery direction and subsequently extracted the centerline. Specifically, there were various points corresponding to different directions distributed on the sphere, and then the direction classification detected the possible point for direction. Combined with their research, Mostafa et al. ^[17] optimally developed a fully automated workflow. They employed a regression U-Net to identify the coronary ostia automatically and applied the improved CNNTracker to track it until extracting its endpoint. This method added a weight term to the loss function, which could calculate the angle between the predicted vector and each sampling vector, significantly improving the performance of the direction classifier. Yang et al. ^[18] presented a discriminative coronary artery tracking method by alternately training a tracker and discriminator. Their discriminator could provide a more robust learning-based stop criterion, effectively distinguishing coronary arteries and other cardiac tissues. However, the classification method used in these studies has relatively limited accuracy. The proposed models must increase the number of sampling points to achieve higher accuracy, undoubtedly causing large memory consumption in calculating training parameters. Besides, the results of multiple branches and distal parts extraction are not satisfying.

In this work, we proposed a regression method that learned to track centerlines continuously from coronary ostia. The regression method can calculate the deviation between the actual direction and the predicted direction, learning the transformation of the network output to the centerline direction. First, the branch classifier in this model could detect the endpoint, bifurcation point, and single-vessel point. Besides, the direction predictor was performed to predict the coronary artery centerline direction and lumen radius. Then, we developed a new loss function that associated the coordinates of the centerline points with the lumen radius, which achieved high accuracy extraction.

2.   Methods

The overall workflow of our proposed network is illustrated in Figure 1. First, the acquired medical images were cropped into computed tomography (CT) patches for the input of the network. The CT patch was centered on ${C}_{i}$, where located on the centerline. Afterward, the CNN module extracted image features which were fed to branch classifier and path predictor. The branch classifier can predict the likelihood of single-vessel point, endpoint, and bifurcation point. Meanwhile, the output of path predictor consisted of direction vector $\left(x, y, z\right)$, the vessel radius r, and step. Subsequently, the coordinate of ${C}_{i+1}$ on the centerline can be calculated by Eq (1). If ${C}_{i+1}$ was found, it will be queued as the new current centerline point until tracking coronary artery centerline endpoint.

Figure 1.  Architecture of the proposed network.

DownLoad: Full-Size Img PowerPoint

where ${C}_{i+1}$ is the next centerline point; ${C}_{i}$ is the current centerline point; $\left(x, y, z\right)$ is the direction vector; step is the spatial distance between ${C}_{i}$ and ${C}_{i+1}$.

2.1.   Feature extraction network

The feature extraction network was composed of two convolution blocks, namely DimConv and DownConv. DimConv first implemented convolution, batch normalization, and ReLU activation function in sequence. The batch normalization was set before the ReLU activation function and limited the passing value between layers within a certain range, which significantly mitigated vanishing or exploding gradient problem during backpropagation. Afterwards, dilated convolution, batch normalization, and ReLU activation function were placed twice. The operation of dilated convolution instead of basic convolution contributed to obtain a broader receptive field without increasing the number of parameters. Besides, the patch size was about three times of the maximum coronary lumen diameter, which contained more neighborhood information to help locate the spatial location of the current vascular segment. DownConv was made up of convolution, batch normalization, and ReLU activation function. It achieved down sampling operation by using convolution with a step size of 2, so that the output image of DownConv module could be reduced by half to obtain the condensed feature expression.

The feature extraction network repeated DimConv and DownConv operations three times. The feature map was the output of the previous layer, by superimposing the channels of feature maps, multi-level features extracted from shallow network and deep network can be fused. In this case, the low-level feature information will not easily lose in the process of down sampling. DimConv can maintain the size of the input and output images while only changing the number of channels. Therefore, the output of DimConv had the same size as the output of the previous DownConv, whose channels can be superimposed and then input into DownConv for further extraction.

2.2.   Branch classifier

A branch classifier adopting spatial attention was developed, which can aggregate feature information and provide a terminal criterion for a continuous tracking process. In this part, the dimension of the input feature map was (64, 4, 4, 4), where 64 represented the number of channels and (4, 4, 4) was the size of 3D feature map. First, we calculated the mean and maximum values of per-pixel across all channels in feature maps respectively, generating two channels. Then, the convolution operation was performed on these two channels, effectively realizing multi-level feature fusion. The dimension of the output of convolution operation was (1, 4, 4, 4) and finally a fully connected layer was added to reach the final classification. The advantage of spatial attention was that dimensionality reduction could be rapidly achieved when integrating global features, which helps the number of input neurons of the fully connected layer more accurately close to the number of output classification. Each output channel of the branch classifier was considered as a feature classification, purely indicating the probability of the endpoint, bifurcation point, and single vessel point. Branch 0 and 2 indicated the endpoint and bifurcation point, respectively. When the output was branch 1, it represented that the current point is on the path of a single vessel, with only one possible direction.

2.3.   Path predictor

For path prediction, the feature maps were first fed into the global average pooling layer for information integration and the dimension of the output was (64, 1, 1, 1). Then, we input it into the ReLU activation layer and obtained nonlinearity. The results were sent to convolution with 1 × 1 × 1 kernel, which could customize dimensionality to satisfy the calculation requirement of the regression loss with different number of parameters. The number of output channels was $2\times 5$, where the maximum number of branches was 2 and each branch contained 5 parameters including the direction vector $\left(x, y, z\right)$, the vessel radius r, and step. The path predictor can predict three branch conditions including no branch, one branch, and two branches. It was worth noting that only when the judgment of branch classifier was valid can the regression loss operation be performed next.

2.4.   Loss function

The classification loss of branch classifier used cross entropy (CE) function, which was defined as follows:

where M is the number of classes; ${p}_{ic}$ represents the probability that sample $i$ belongs to class c; When the actual class of sample $i$ is equal to c, the value of ${y}_{ic}$ is 1, otherwise is 0.

To accurately segment coronary artery in CTA images, numerous regression methods have been proposed, including mean square error (MSE), mean absolute error (MAE), and other optimizations such as smooth L1 Loss piecewise function ^[19]. The core idea of vascular path regression is precisely learning the spatial coordinates of unit direction vector $\left(x, y, z\right)$ in the Cartesian coordinate system. Regression loss can be defined as the transformation of the predicted direction vector to the ground truth. However, when separately calculating the losses of parameters (e.g., x, y, z, r, and step), it is ultimately easy to ignore the correlation between these parameters (e.g., x, y, z, r, and step). This phenomenon is particularly serious in dealing with the centerline extraction problem of the coronary artery with extreme tortuosity, which increases the difficulty of model convergence.

Currently, the maximum endovascular sphere model has made significant achievements in the field of lumen centerline extraction ^[20,21]. This approach is based on the assumption that a sphere, whose diameter is adaptively varying to enclose the diameter of the vessel, moves in the direction of blood flow. When the sphere encounters the bifurcation, it will split into two spheres and continues tracking each branch parallelly. After tracking the current vessel end, the coronary centerline is formed by connecting the center of the sphere. Based on the idea of endovascular sphere model, we have designed a new loss function method, referred to as the spherical cross volume method, for vascular path tracking. While learning the deviation of the output parameters (e.g., x, y, z, r, and step), this suggested method can ensure the correlation among these parameters (e.g., x, y, z, r, and step).

Figure 2 shows two spheres centered on $O$ and ${O}_{t}$, respectively. The volume of the intersecting part of the two spheres is defined by the following formula:

Figure 2.  Illustration of sphere dice loss.

DownLoad: Full-Size Img PowerPoint

where $r$ and ${r}_{t}$ are the radiuses; $h$ and ${h}_{t}$ are obtained by taking radius $r$ and ${r}_{t}$ minus the distance between the center of the sphere to the intersecting plane, respectively; $cos\alpha$ and $cos\beta$ can be obtained as the following:

where $d$ is the distance between $O$ and ${O}_{t}$. Combining Newton-Leibniz formula and Eqs (3)–(5), ${V}_{\cap }$ can be simplified as follows:

We formulated the sphere dice loss function similar to the Dice coefficient and extended it to 3D space. Dice Loss is a typical tool using in image segmentation tasks, which can count the number of pixel points in the segmentation area and target area. Then the similarity of the segmentation area and target area can be described by the ratio of the number of points in the intersecting area to the total number of pixel points. We applied the concept of Diss Loss to geometry. Correspondingly, the similarity of two spheres was illustrated as the ratio of the volume of the intersecting parts to the total volume of two spheres. The value of similarity ranged from 0 to 1, while 0 and 1 denoted the two spheres non-intersection and perfectly coincidence, respectively. When two spheres partially intersected, their deviation can be obtained by subtracting the actual similarity from 1. Specifically, the regression loss can be calculated by Eq (7):

It is acknowledged that a sphere can be constructed based on the coordinates (x, y, z) and radius r in 3D space. Thereby, the maximal inscribed sphere located at the current coronary lumen can be characterized by the annotated coordinates of coronary centerline points and the actual vessel radius. As described above, we constructed two sphere, one centered on the current centerline point and another sphere centered on the predicted centerline point, relying on the output parameters (e.g., x, y, z, r, and step) by the path predictor. Subsequently, the sphere dice loss was calculated according to the location relationship between the two spheres, which achieved outstanding performance in evaluating the deviation between the predicted direction and the actual direction.

The combined loss was the sum of the CE loss (Eq (2)) and the sphere dice loss (Eq (7)). Note that the diameter of the vessel end was small, which easily resulted in a larger loss even from a small position deviation of the sphere.

3.   Experiment and results

3.1.   Description of datasets

The dataset consisted of 18 CTA scans provided by The Second Affiliated Hospital of Zhejiang University School of Medicine, whose resolution and voxel size were described in Table 1. We randomly picked 12 CTA scans for training and 6 CTA scans for testing. The requirement of written informed consent for participation was waived by the local ethics committee due to the retrospective nature of this study. The coronary centerlines of all images were manually annotated using Mimics 21.0 software (Materialise N·V., Belgium), which were corrected and confirmed by an experienced radiologist. Then, the manually marked centerlines were smoothed and the centerline points were resampled equidistantly. In addition, a resampling operation was performed on the CTA images, with uniform voxel spacing of 0.5 mm. Afterward, the CTA images were cropped into several CTA patches with the size of 32 × 32 × 32, which were centered at the centerline points. To achieve data augmentation, we made random offset for the location of centerline point and over-sampled nearby the bifurcation and the vessel end. As a result, in the training set and the testing set, 96,834 and 18,072 patches were extracted respectively.

Table 1.  Information of CTA datase.

CTA Data	Size (Voxels)	Resolution (mm3)
Train01	(512,512,275)	(0.422, 0.422, 0.500)
Train02	(512,512,275)	(0.391, 0.391, 0.500)
Train03	(512,512,275)	(0.422, 0.422, 0.500)
Train04	(512,512,297)	(0.447, 0.447, 0.500)
Train05	(512,512,275)	(0.367, 0.367, 0.500)
Train06	(512,512,275)	(0.422, 0.422, 0.500)
Train07	(512,512,275)	(0.422, 0.422, 0.500)
Train08	(512,512,275)	(0.422, 0.422, 0.500)
Train09	(512,512,275)	(0.422, 0.422, 0.500)
Train10	(512,512,256)	(0.488, 0.488, 0.625)
Train11	(512,512,256)	(0.488, 0.488, 0.625)
Train12	(512,512,256)	(0.488, 0.488, 0.625)
Test01	(512,512,275)	(0.422, 0.422, 0.500)
Test02	(512,512,275)	(0.422, 0.422, 0.500)
Test03	(512,512,256)	(0.488, 0.488, 0.625)
Test04	(512,512,256)	(0.488, 0.488, 0.625)
Test05	(512,512,275)	(0.367, 0.367, 0.500)
Test06	(512,512,275)	(0.391, 0.391, 0.500)

---

 | Show Table

DownLoad: CSV

3.2.   Experiments

We conducted the experiments on an NVIDIA GeForce GTX 1660 GPU. The algorithm was implemented in Python using PyTorch, an open-source deep neural network library. The continuous centerline extraction started from coronary ostia. The angle between the predicted direction vector and the previous direction vector was computed to prevent the following tracking path from overlapping with the already extracted centerline.

To verify the accuracy of the loss function, we respectively applied smooth L1 Loss function and sphere dice loss function to calculate regression loss, as shown in Figure 3. We can observe that the loss curve using smooth L1 Loss function oscillates obviously and is difficult to converge. This is probably attributed to the fact that the process for calculating deviations in the x, y, and z directions using the smooth L1 Loss function is separate. As the coronary artery structure is complex, the deviations in the x, y, and z directions are uncertain, resulting in poor performance. Additionally, the sensitivity of various sizes of coronary artery segments to deviation is different, which increases the difficulty of model convergence. Compared to the smooth L1 Loss function, the sphere dice loss function combines the deviations in three directions and relates deviations to the lumen radius, which can be effectively adapted to the complex coronary artery structure. As expected, the loss curve applying the sphere dice loss function is observed to converge quickly and smoothly.

Figure 3.  The loss curve in 200 epochs using different loss function. (a) the loss curve using smooth L1 Loss function. (b) the loss curve using sphere dice loss function.

DownLoad: Full-Size Img PowerPoint

3.3.   Evaluation measures

We evaluated the accuracy of centerline extraction by adopting the Rotterdam algorithm evaluation framework rules ^[22]. Four evaluation metrics including overlap (OV), overlap until first error (OF), overlap with clinically relevant vessel (OT), and average inside (AI) were measured, which can be computed by true positive method (TPM) and true positive reference (TPR). The TPR point refers to the point of the reference centerlines, whose defined radius is no less than the distance between this point and at least one connected point on the extracted centerlines. Then, The TPM point refers to the point of the extracted centerlines, whose defined radius is no less than the distance between this point and at least one connected point on the reference centerlines. OV stands for the completeness of the predicted coronary artery centerline. OF reveals the tracking ability of the model before making the first error. OT indicates the degree of overlaps between the extracted centerlines and the clinically relevant centerlines (radius   $\ge$ 0.75 mm). AI is an accuracy index, which measures the average distance between the extracted point on predicted centerline and the corresponding extracted point on reference centerline within the radius of the reference centerline.

Table 2 shows the average overlap and accuracy results for centerline extraction in all testing datasets. For each dataset, the OV, OF, OT, and AI results are listed. In terms of overlap, our method obtained an average OV of 89.19%, an average OF of 82.30%, and an average OT of 91.42%. In terms of accuracy, our method obtained an average AI of 0.32mm. Besides, we found that errors mainly occurred at bifurcations or sharply curved coronary arteries.

Table 2.  Evaluation results of the proposed method. Four metrics are reported: Overlap (OV), overlap until first error (OF), overlap with the clinically relevant part of the vessel (OT), and Average inside (AI).

Dataset	OV (%)	OF (%)	OT (%)	AI (mm)
Test01	91.87	87.41	90.82	0.32
Test02	94.53	90.47	95.18	0.26
Test03	88.64	89.58	94.55	0.28
Test04	82.67	63.45	86.84	0.40
Test05	85.04	73.72	87.93	0.36
Test06	92.36	89.18	93.21	0.29
Average	89.19	82.30	91.42	0.32

---

 | Show Table

DownLoad: CSV

To analyze the advantages of our proposed method, we conducted the ablation experiments. Three conditions were performed: (i) removing the Spatial Attention module and replacing it with basic convolutional layer; (ii) removing the global average pooling layer in path predictor; (iii) simultaneously removing the Spatial Attention module and the global average pooling layer in path predictor. The average OV, OF, OT, and AI results were calculated and illustrated in Table 3. Significant decreases in performance were observed in the ablation experiments compared with our proposed network, which proved that the specific feature integration method designed for the output task was effective.

Table 3.  The average OV, OF, OT, and AI results of ablation experiments.

Methods	Av.OV (%)	Av.OF (%)	Av.OT (%)	AI (mm)
proposed network	89.19	82.30	91.42	0.32
-Spatial Attention	85.83	72.67	90.55	0.36
-Avg Pool	88.26	75.41	89.97	0.39
-Spatial Attention & Avg Pool	84.45	70.86	89.23	0.38

---

 | Show Table

DownLoad: CSV

According to the guidelines of the American Heart Association (AHA) ^[23], the coronary arteries with a diameter $\ge$ 1.5 mm can be divided into 17 segments. Figure 4 shows the sensitivity of the proposed method for each of these 17 segments in the testing datasets. The testing sensitivity for right coronary artery (RCA), left anterior descending artery (LAD), left circumflex artery (LCX), right and left posterior descending arteries were 83.3–100%, while for the second obtuse marginal branch was 66.7%.

Figure 4.  Sensitivity of 17 segments of coronary arteries using the proposed network.

DownLoad: Full-Size Img PowerPoint

Table 4 shows the comparison results of our proposed method and the CNNTracker approach ^[16]. We reproduced the CNNTracker approach using our dataset, while the performance was not ideal. This may be due to the small size of our dataset. However, there was evidence that our method can improve the sensitivity of detecting distal segment of coronary arteries, especially for margin branches, left posterior descending branch, and other small vessels. Consequently, our proposed method has advantages in extracting more complete coronary artery structures.

Table 4.  The sensitivity of detecting distal segment of coronary arteries.

	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17
CNNTracker	1	1	1	0.9	1	1	1	0.9	0.9	0.8	1	0.8	0.9	0.6	0.7	0.7	-
proposed method	1	1	1	0.8	1	1	1	0.8	1 ↑	0.8	1	1 ↑	1 ↑	0.6	0.8 ↑	1 ↑	0.8

---

 | Show Table

DownLoad: CSV

Except for the 17 reference segments, our method can also detect acute marginal branches and conus arteriae branches, clearly reflecting the strength of our method. As shown in Figure 5, the green and red coronary arteries both started from coronary ostium and extracted to the end of blood vessel. It was evident that the computed coronary artery centerlines were longer than that obtained merely by manual annotations.

Figure 5.  3D results of extracted coronary arteries by our method and manual annotations. The blue spheres indicate the location of the coronary ostia. The red coronary arteries were acquired by our method and the green coronary arteries were obtained by manual annotations.

DownLoad: Full-Size Img PowerPoint

Figure 6 provides the visual comparison between the results of manual annotations (Figure 6(a)) and our approach (Figure 6(b)), which verifies the effectiveness of our approach in predicting distal vessel branches without annotations. Additionally, our approach was able to successfully extract coronary artery centerlines in cases where intravascular stents (Figure 6(c)) or coronary artery calcifications (Figure 6(d)) were present.

Figure 6.  Straightened CPR (curved planar reconstruction) of coronary artery centerline. (a) the centerline extracted by manual annotations; (b) the centerline extracted by our method; (c) the centerline extracted by our method in the presence of intravascular stent; (d) the centerline extracted by our method in the presence of coronary calcification.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

In this paper, we have presented a new method for continuously tracking robust coronary arteries using a regression-based network. The proposed network first uses the CNN module to extract CTA image features, and then enters the branch classifier and the path predictor. There, the direction of the coronary centerline and the radius of lumen can be identified. Previous centerline extraction methods based on direction classification have certain limitations ^[16,17,18]. If the directions of sampling are too few, the direction deviation will accumulate with the continuously tracking process, and the tracking process will stop when a large deviation occurs at the distal coronary artery. This phenomenon results in the difficulty to the centerline extraction of small vessels. Whereas, dense sampling directions likely reflect the imbalance issue between positive and negative samples associated with direction classification, resulting in extraction failure. Besides, the direction classification method can only predict a single direction, which need to train additional seed points to deal with multi-branch problems. Due to the insufficient learning of effective information for bifurcation, some extracted centerlines of small coronary vessels cannot cross the bifurcation instead tracking back to the coronary artery ostia. While in our proposed network, the branch classifier can creatively handle multi-branch problems and provide a termination criterion for continuous extraction. Moreover, the path predictor can learn the transformation of the network output to the actual centerline direction based on the regression method, significantly improving the accuracy of direction prediction. Further, the proposed new loss function comprehensively considers the relationship between coordinates and lumen radius, which can adapt to the tortuous coronary structures, enhancing the ability to track small coronary vessels.

We conducted the tests using 6 CTA images to validate the reliability of the proposed network. The testing results demonstrated that the suggested method can extract the whole coronary artery tree structures, with a high performance of OV 89.19%, OF 82.30%, and OT 91.42%. While the manually annotated centerline cannot reach the most distal point of coronary artery due to the presence of nearby myocardium and adjacent coronary veins, our method allows more complete tracking of coronary artery structures. There is evidence that distal coronary arteries or small vessels also play an important role in coronary circulation, which can lead to myocardial ischemia ^[24]. Therefore, the extraction of distal coronary arteries and small vessels in our study could provide more comprehensive information for clinical diagnosis or operation. However, our suggested method still requires the starting point, whose coordinate should be set manually at coronary artery ostia. In future work, some interactive operations such as manually adding the starting points of the missing branches can be considered when the tracker drifts from the centerline or the extraction is incomplete.

The hyper-parameters directly affect computational complexity and method performance, which should be modified to achieve effectiveness. During the hyper-parameter exploration, we trained the CNN models for 200 epochs with a 2e-5 initial learning rate and 0.8 decay rate per 50 epochs. The model was trained using Adam optimizer with a batch size of 32 in which the model can achieve better performance by observing accuracy metrics (e.g., OV, OF, OT) and visualization results. After experimenting with various batch size, the results reveal that the training curve would oscillate when the batch size was too small, which may be due to the individual coronary artery structures and variable vessel sizes. Whereas larger batch size will cause the gradient descent to be too slow, which may be attributed to the calculation of the average loss of large amounts of data affecting the learning efficiency. As a consequence, the batch with 32 patches was found suitable as it could reflect conditions such as different numbers of branches, various coronary artery sizes, and whether deviation existed or not.

Despite promising results in extracting coronary centerlines, our proposed framework suffers from bias, which is caused by the generation inconsistency between the training and inference process. To be specific, at training time, the input of the path prediction network is the CT patch, and the output is the parameters (e.g., x, y, z, r, and step) of the predicted centerline point. Overall, the training process is to learn features and optimize model from the image data. While at inference the trained model is constantly used to estimate outcomes until a terminal criterion is fulfilled. Hence, as the training grows, the errors accumulate. To address this problem, this work uses data augmentation which allows the tracker to make corrections to process off-centerline points. However, considering that the coronary artery has delicate tubular structures and complex tree shapes, the continuous path tracking method remains challenging as the contextual information of multiple patches is independent. In this case, some recent studies combined the centerline topological structure. For example, Gao et al. ^[25] constructed a hybrid model of CNNTracker and graph convolutional network (GCN), with the advantage of combining coronary centerline extraction and lumen segmentation. This approach has achieved state-of-the-art performance in capturing the geometry of the RCA, LAD, and LCX. Besides, Jeon et al. ^[26] presented a deep regression Bayesian network for extracting the main coronary artery centerline, creatively detecting branching sites by adopting a particle filter. The proposed method performed more effectively when extracting coronary arteries with sharp curvature. However, the extraction process combined with topological structure is complicated and needs large amounts of training data. In the future, based on more-powerful computers and large-scale data, the topological structure can be considered to improve the accuracy of coronary centerline extraction.

Furthermore, our proposed method still has some limitations. Firstly, all cases in this work had nonsignificant coronary stenosis, whose Fractional Flow Reserve (FFR) was greater than 0.80 ^[27,28]. It is difficult to extract coronary centerline in the presence of severe epicardial stenosis (diameter severity [DS] $\ge$ 70%) ^[29], hence, the multi-scale model would be an interesting exploration in the future. Secondly, a major limitation is the low number of cases included in this study. Generally, the CTA examination is more suitable for early screening of CAD and the prognosis of interventional treatment. A large number of nonsignificant coronary stenosis patients have not received sufficient attention and complete examination in clinical practice. Therefore, it is difficult to collect large amounts of data in a short period. Besides, to more accurately identify coronary small vessels, we used 256-slice spiral CT images in this study, whose high cost also increased difficulties in data collection. In future studies, large-scale studies can be considered to improve the accuracy of coronary centerline extraction and optimize our method for severe coronary artery stenosis. However, a large training set requires time-consuming manual annotations. In this case, semi-automatic supervision method is recently used in medical image analysis, which can effectively address the insufficient data problem and has achieved significant performance ^{[30,31,32,33]}. Accordingly, our approach can be further investigated by combining the semi-automatic supervision method.

5.   Conclusions

In summary, we proposed a regression-based tracking method for coronary artery centerline extraction. Our approach can effectively deal with multi-branch problems and realize further coronary artery extraction. Future work can combine multi-scale model or semi-automatic supervision method for improvement, powerfully assisting CAD diagnoses.

Acknowledgments

This work is supported by the Natural Science Foundation of China (NSFC) under grant number 62171408, and the Key Research and Development Program of Zhejiang Province (2020C03060, 2020C03016, 2022C03111).

Conflict of interest

The authors declare there is no conflict of interest.