Coronary artery segmentation based on transfer learning and UNet Architecture on Computed Tomography Coronary Angiography Images

Received 6 June 2023, accepted 26 June 2023, date of publication 7 July 2023, date of current version 26 July 2023. Digital Object Identifier 10.1109/ACCESS.2023.3293090 Coronary Artery Segmentation Based on Transfer Learning and UNet Architecture on Computed Tomography Coronary Angiography Images BELÉN SERRANO-ANTÓN 1,2,3, ALBERTO OTERO-CACHO 1,2,3, DIEGO LÓPEZ-OTERO4,5, BRAIS DÍAZ-FERNÁNDEZ4,5, MARÍA BASTOS-FERNÁNDEZ 4,5, VICENTE PÉREZ-MUÑUZURI 3,6, JOSÉ RAMÓN GONZÁLEZ-JUANATEY4,5, AND ALBERTO P. MUÑUZURI 2,3 1FlowReserve Labs S.L., 15782 Santiago de Compostela, Spain 2CITMAga, 15782 Santiago de Compostela, Spain 3Group of Nonlinear Physics, University of Santiago de Compostela, 15782 Santiago de Compostela, Spain 4Cardiology and Intensive Cardiac Care Department, Hospital Clínico Universitario de Santiago de Compostela, 15706 Santiago de Compostela, Spain 5Centro de Investigación Biomédica en Red de Enfermedades Cardiovasculares (CIBERCV), 28029 Madrid, Spain 6Institute CRETUS, Group of Nonlinear Physics, University of Santiago de Compostela, 15705 Santiago de Compostela, Spain Corresponding author: Alberto Otero-Cacho ([email protected]) This work was supported in part by Spanish Ministerio de Economía y Competitividad and European Regional Development Fund under Contract RTI2018-097063-B-I00 AEI/FEDER, UE; in part by Xunta de Galicia under Grant 2021-PG036; in part by Spanish Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 through the Industrial Doctorates Grant and European Union NextGenerationEU/PRTR Research under Grant DIN2020-011068; and in part by FEDER (UE). This work involved human subjects or animals in its research. Approval of all ethical and experimental procedures and protocols was granted by the Research Ethics Committee Santiago-lugo. ABSTRACT Coronary artery segmentation from CT scans is a helpful tool for coronary artery diseases diagnosis, which is frequently characterised by a vessel narrowing (stenosis). This is a highly demanded and high time-consuming process, thus automated procedures are becoming increasingly necessary. In this work, we propose an extremely light computationally 2D UNet that uses transfer learning for the first time in CT images. We compare the results, using different architectures and backbones, of a 2D UNet and a 3D UNet trained from scratch (i.e. weights are randomly initialised) and a 2D EfficientUNet. Both the amount of input data, with a total of 88 patients, and the extension of the structure to be recognised, the aorta and the coronary arteries (A+C.A), as well as the coronary arteries only (C.A) are analysed. Network outputs in clinically identified stenotic lesion areas are also assessed. The results show the advantage of using transfer learning when data is scarce, improving the F1score by up to 0.6 points for the 2D UNet. On the other hand, when data is sufficient, F1score values are close to 0.9 for all the networks. Besides, the results reveal that the 2D UNet distinguishes the thinnest and most distal vessels, although in the presence of a lesion, there is a clear tendency to overestimate it. The network with the best accuracy is the 3D UNet, with values above 95% and 75% in A+C.Aand C.A, respectively. Moreover, the proposed methods show dependence on the amount of training data and dataset structure (A+C.Aor C.A). INDEX TERMS Artery, convolutional neural network, coronary, CT, segmentation, UNet. I. INTRODUCTION Coronary artery disease (CAD) is one of the leading causes of death worldwide [1], [2]. Invasive tests, such as coronary The associate editor coordinating the review of this manuscript and approving it for publication was Juan A. Lara . angiography, with high spatial resolution and the possibility of intervening in the same procedure, can be used to diagnose this disease. However, due to the associated risks, noninvasive techniques such as computed tomography (CT) are increasingly being employed [3]. Coronary computed tomography angiography (CCTA) allows the extraction of 75484 This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 License. For more information, see https://creativecommons.org/licenses/by-nc-nd/4.0/ VOLUME 11, 2023 B. Serrano-Antón et al.: Coronary Artery Segmentation images of the heart in slices. If a contrast agent (dye) is used, these images show whether the arteries are narrowed or not. Additionally, each scan can be integrated to create a 3D visualization that greatly aids diagnosis. After the images are taken, lesions are visually assessed by an expert. Nonetheless, this is a tedious, error-prone, subjective and time-consuming process. One of the solutions to speed up the process and obtain an accurate diagnosis is to obtain a precise 3D geometry of the region of interest (aorta and coronary arteries). This enables not only the visualisation of the structure but also the calculation of other clinically-relevant diagnostic parameters such as Fractional Flow Reserve (FFR), using computational fluid dynamics or machine learning techniques [4], [5], [6], [7]. Semi-automatic and automatic segmentation (pixel-bypixel identification of structures) techniques are used to obtain these geometries. Although semi-automated techniques involve many tools, such as threshold structure filtering on Hounsfield units, centerline extraction, or 3D editing, they still require time and expertise. In automatic segmentation techniques, artificial intelligence (AI) is increasingly being applied, for example Convolutional Neural Networks (CNN) [8], [9]. The main advantage of these methods is that they are not user-dependent and can therefore be used by both expert personnel and personnel from areas outside the field of medical imaging. This is helpful for research and reproducibility of results, as well as saving a lot of time. Within the architectures described by CNN, UNet has played a relevant role in biological image segmentation [10]. This type of network consists of two parts: an encoder that reduces the spatial dimension of each layer and increases the number of channels, and a decoder that performs the reverse process. Some modifications of this architecture can be found in the literature in the field of coronary structure segmentation, such as 3D versions of UNet itself (3D UNet), VNet or UNet++ [11], [12], [13], [14], [15], [16]. Cheung et al. [12] propose a modified and computationally cheap 2D UNet that obtains a dice similarity coefficient (DSC) [17] of 91.20% and 88.8% when trained with the entire coronary tree (aorta and coronary arteries) and with the coronary arteries, respectively. Their dataset consists of 71 CT scans. Other works such as [13], propose a computationally intensive model based on 3D UNet with a dataset of 474 CT scans. In this case, the DSC in the coronary arteries increased to 96.91%. Huang et al. [14] use a 3D UNet with two settings, depending on whether the data include vessel centerlines (18 CCTA) or not (34 CCTA). Duan et al. [18] suggests a 3D UNet with prior spatial knowledge of coronary geometries, which allows to reduce the complexity of the model. The DSC received here is 71.46%. On the other hand, Gu and Cai [15] combine the advantages of 2D UNet (wide field of view) and 3D (its ability to favor structural continuity) in a 2D-3D UNet network. This implementation achieves a DSC of 97.54% in aorta and 86.62% in coronary arteries. This method achieves a DSC of 79.5% in the coronary arteries. An additional benefit of CNNs is their simplicity to transfer learning. This is particularly intriguing when the quantity of labeled data, time and assets are restricted [19]. It is not business as usual, consequently, that transfer learning procedure has been applied to clinical imaging and, specifically, to coronary artery and atherosclerotic plaque segmentation [20]. For example, [21] utilizes a pretrained (VGG16, ResNet50, and Inception-v3) CNN for stenosis region location in coronary angiography through transfer learning and fine tuning of a synthetic dataset of 10,000 pictures tested with various setups. With the best design, they got a F1score of 0.98. Candemir et al. [22] additionally features the benefits of utilising transfer learning in architectures, for example, VGG-16 [23], particularly when the dataset is small, as their weights have been optimised through more than one million pictures from ImageNet [24]. In any case, they additionally bring up the absence of pre-trained 3D networks. All the articles mentioned above, including this one, face some problems with respect to coronary artery segmentation. 1) Data scarcity. The results obtained by neural networks depend, to a large extent, on the input data. In the case of coronary artery segmentation, images are difficult to obtain (as they are clinical data) and require a previous manual segmentation. 2) The datasets obtained are class imbalanced, as the coronary vessels constitute a very small region within each slice. 3) Some approaches such as 3D networks are computationally very intensive, either in terms of memory or in computation time. The main contributions of the present study are: 1) Study of the dependence of the number of patients for obtaining an accurate coronary tree with UNet neural networks. 2) Comparison of 2D (with different backbones) and 3D networks. As a novelty to previous studies, transfer learning in 2D is included. 3) Training with a very detailed and accurate dataset. This article is structured as follows. Section II describes the methods, including image acquisition, dataset processing and technical specifications of the algorithms. Section III details the performance of the algorithms and shows examples of the volumes obtained. Finally, section IV contains the discussion and conclusions. II. METHODS A. CLINICAL DATA AND PRE-PROCESSING The dataset used for this study contains CT images from 88 patients that have been obtained at the University Hospital in Santiago de Compostela (Spain). CCTA is obtained for each patient consisting of 256 images with a 512x512 pixel resolution. Patients were chosen based on a clear image criterion and the lack of calcium-related lesions. If a patient develops some calcification (less than 5% of the patients), the calcium is washed out of the vessel, while maintaining VOLUME 11, 2023 75485 B. Serrano-Antón et al.: Coronary Artery Segmentation the blood flow region. CCTA images were obtained using a Revolution CT [25], with a one-beat 16cm wide coverage and 0.23mm spatial resolution. Acquisition parameters and patient premedication were chosen following the Society of Cardiovascular Computed Tomography recommendations. The images were cropped to a size of 400x400 pixels in order to reduce their size and non-essential information. The amount of the images per patient was also decreased from 256 to 200. Note that the quality of the images was never decreased and neither the aorta nor any other coronary structures of importance were eliminated. The Hounsfield Units used to measure the pixel intensity have been rescaled to yield values between 0 and 255. B. ETHICS STATEMENT The development of the project was carried out respecting the Declaration of Helsinki of the World Medical Association 1964 and ratifications of the following assemblies (Tokyo 75, Venice 83, Hong Kong 89, Somerset West 96, Scotland 00, Seoul 08 and Fortaleza 13) on ethical principles for medical research on human beings, RD 1090/2015, of December 24, on clinical trials, specifically the provisions of article 38 on good clinical practices, and the Convention on human rights and biomedicine), made in Oviedo on April 4, 1997 and successive updates. The researchers participating in this study agree that all clinical data collected from the study subjects will be separated from personal identification data ensuring the anonymity of the patient; respecting the Personal Data Protection Law (Organic Law 15/1999, of December 13), RD 1720/2007 of December 21, which approves the Regulations for the development of Organic Law 15/1999, Law 41/2002, of November 14 (basic regulation of patient autonomy and rights and obligations in terms of information and clinical documentation), as well as Law 3/2001, of May 28, (regulator of informed consent and the clinical history of patients), Law 3/2005, of March 7, modifying Law 3/2001 and Decree 29/2009 of February 5, which regulates access to history electronic clinic. The clinical data of the patients will be collected by the investigator in the Case Report Form (CRF) specific to the study. Each CRF will be encrypted, protecting the identity of the patient. Only the research team and the health authorities, who have a duty to maintain confidentiality, will have access to all the data collected for the study. Only information that cannot be identified may be transmitted to third parties. Once the study is finished, the data will be destroyed. The treatment, communication and transfer of data will be done in accordance with the provisions of the General Data Protection Regulation (Regulation (EU) 2016/679 of the European Parliament and of the Council, of April 27, 2016). The data collected will only be used for the purposes of the research study described in the protocol and kept for the time necessary to achieve the objectives of the study and in accordance with applicable legislation. As this is a retrospective study of medical records and archived samples that does not deviate from routine clinical practice, the Ethics Committee consider FIGURE 1. Example of patient segmentation. A) Slices of the CT image with the segmentation line of the aorta and the coronary arteries in red. B) 3D geometry of the coronary tree obtained from the slices. that patient informed consent and fully anonymization of the data before being access are sufficient requirements to carry out the study. C. GROUND TRUTH SEGMENTATION The process to obtain ground truth labels was carried out semi-automatically with 3D Slicer software (version 4.11.20210226) [26], [27]. The task was designed in three steps: 1) thresholding, 2) cleaning and 3) refinement. 1) Thresholding: A threshold is established to extract coronary structures with good contrast. During this process not only the aorta and coronary arteries are segmented, but also other parts of the muscle or ventricles. 2) Cleaning: From the previous step, structures other than the aorta or coronary arteries are removed. 3) Refinement: Narrow and distal vessels are segmented by thresholds lower than in step 1). The refinement step distinguishes our study from previous works as it performs a thorough segmentation of the entire coronary tree (see Fig. 1). Once the aorta and coronary arteries have been segmented, the aorta is removed, leaving only the coronary arteries to generate the two study datasets. One with the complete structure and the other set with only the coronary arteries. Finally, the volume is downsized to 200x400x400 and a binary mask is applied to obtain the ground truth labels (GT). D. METHOD 1) MODEL The UNet architecture underlies the models we compare. The first model is a 2D UNet. The peculiarity of this network is that it has an encoder that has already been trained through more than one million pictures from ImageNet [24] and whose weights cannot be changed, as well as a decoder that can be trained. The encoder is based on MobileNetV 2 [28] and is already implemented in tensorflow.keras.applications [29] (see model2Dpre.json file in Supplementary Material). This network was chosen because of its low computational cost and high performance. It is 32 times smaller and 27 times less computationally 75486 VOLUME 11, 2023 B. Serrano-Antón et al.: Coronary Artery Segmentation FIGURE 2. Segmentation ground truth (GT) of the 10 test patients. intensive than VGG16 while maintaining competitive accuracy [28]. Furthermore, by using a pre-trained encoder and keeping these weights constant, the number of variable parameters can be significantly reduced. In this case, the total number of parameters is 6,502,786, of which 4,658,882 are trainable and 1,843,904 are non-trainable. Because the input size in the first implementation of this network is fixed at 512x512, the 400x400 images of our dataset are filled with a black frame until the required size is obtained. This allows us to use the same dataset for both 2D and 3D algorithms. A 3D UNet is the second model. The input and output of this network are volumes made up of 16 images (CCTA slides). Following the implementation of [13], each 16-image patch contains information from the previous and next patch. This network has a total of 22,575,329 trainable parameters. The complete architecture of both models is in the Supplementary Material in Fig. S1 and Fig. S2. In addition, a diagram of the 3D UNet is shown in Fig. S3. 2) COMPARISON WITH OTHER NETWORKS To strengthen the study, the results of the pre-trained 2D UNet and the 3D UNet are compared with a 2D UNet, whose weights are randomly initialised (in what follows we will refer to it as ‘‘from scratch’’), and a pre-trained 2D UNet but with a different backbone, in this case an EfficientUNet [30]. Like the 2D pre-trained UNet, this network has also been pre-trained with the ImageNet dataset [24] and only the decoder is trained. The random initialised 2D network consists of 7,760,322 parameters and details of its architecture can be found in Supplementary Material, Fig. S4. Furthermore, the architecture and implementation of the EfficientUNet is based on a efficientnetb5 whose original code and information about its implementation can be found in [31]. The training datasets are exactly the same as for the 2D pre-trained UNet and 3D UNet, so an ‘‘apple to apple’’ comparison can be made. 3) TRAINING One goal of this work is to determine the relationship between the amount of data (in this case, patients) and the segmentation results. Patients were divided into training, validation, and test sets to achieve this goal, following train-test split strategy. The test set consists of 10 patients whose segmentation (ground truth (GT )) can be seen in Fig. 2. The number of training patients, N, increases in steps of 10 in the range [15,65], and the number of patients in the validation set corresponds to 20% of N. That is, if N=55, there are 55 ∗0.2=11 patients in the validation set. The same is truth for both A+C.Aand C.A. Such studies are not commonly found in the literature, however, they allow to gain intuition about the number of data needed to obtain quality coronary geometries. This is especially important when dealing with hard-to-obtain data, such as clinical data. In the case of 2D UNet from scratch and 2D EfficientUNet N=15 and N=65 was considered, since its implementation is intended only for comparison purposes. 4) IMPLEMENTATION The implementation of the models is developed in Python (v 3.7) [32], [33] and Tensorflow Keras API (v 2.6-tf) [29], [34]. The loss function we used for the study is binary cross entropy [34], [35], [36]. The problem of class imbalance has been tackled by giving five times more weight to the vessel class than to the background class, weighted cross entropy (focal loss). For the optimisation of the models, the Adam algorithm [34], [37] has been used. Convolutional kernel size is set to 3. A MaxPooling layer follows the convolutional layer with a stride of 2 and its followed by a ReLU activation function for VOLUME 11, 2023 75487 B. Serrano-Antón et al.: Coronary Artery Segmentation FIGURE 3. Example grow-shrink algorithm. A) Result segmentation from UNet. B) Grow step. C) Shrink step. both 2D and 3D networks. All information about the networks structure can be found in Section S.I in Supplementary Material. Each training consists of 50 epochs. In order to avoid overfitting, Early Stopping is used with patience =5 and monitor =val_loss (loss function evaluated on the validation set). E. POST-PROCESSING Two algorithms were used to generate cleaner and more connected structures. The first of them is based on the removal of components not connected to the main structure (islands) with values less than 50 voxels. This method will be referred to as Iin the following. Grow-shrink algorithm, on the other hand, is applied with a width of 2mm (see Fig. 3). Throughout the text, this method will be referred to as G. These algorithms were applied sequentially in both I)-G) and G)-I) orders. F. EVALUATION METRICS The fit goodness of the networks considered was evaluated, and, for that, 10 test patients were chosen (see Fig. 2). It should be noted that their manual segmentation is extremely detailed (very thin and distal vessels are included). Furthermore, due to variations in image contrast and the shape of the coronary artery tree, this test set is heterogeneous in terms of geometric complexity. The parameters chosen for the evaluation are well known in the field of segmentation. These include true positives (TP), false positives (FP), true negatives (TN) and false negatives (FN). The following indicators were defined based on these [17], [38]: •Precision: TP/(TP +FP). •Recall: TP/(TP +FN). •F1score or dice similarity coefficient (DSC): is the harmonic mean of precision and recall (TP/(TP +0.5∗ (FP +FN))). It can also be expressed as 2 ∗(Y∩ b Y)/(|Y|+| b Y|). Where Yand b Yrepresent ground truth and prediction, respectively, and Y, b Y∈ {0,1}. •Fa 1score (F1score attached): since we found false positives isolated from the coronary structure in the results, we define this Fa 1score in the same way as F1but only taking into account those FPs that are in contact with a TP. This quantifies how difficult it is to clean the final result of the network if we want to obtain a coronary tree without external structures. Moreover, it is an indicator of oversizing, which is critical if the purpose of coronary geometry is to calculate clinical parameters such as FFR [5]. Another requirement for clinical use of the networkpredicted geometries is that the entire geometry is composed by only one connected component (CC). As a consequence, the number of connected components in the coronary regions of interest is also calculated. To complete the statistical study and visualise the performance of the networks, Bland-Altman plots are included for the area of the coronary geometry in each axial slice of the volume (see Fig. S20 to Fig. S23 in Supplementary Material). G. LESION EVALUATION This work includes an evaluation of the network prediction in two stenosis diagnosed as severe by clinicians. This is critical because the network’s accuracy should improve in the lesionaffected regions. This assessment compares the segmentation performed by an expert (ground truth (GT)) to that predicted by the network from a geometrical standpoint. Although other papers, such as [39], assess lesions using the F1score, a geometric and visual assessment was preferred in this case because it is more informative than the F1score since over/undersized areas can be seen and measured. The volume at the intersection and the difference in the lesion region are provided, as well as their graphical representation. III. RESULTS A. GEOMETRY PREDICTION The network segmentation for test patient T001 is shown in Fig. 4. It depicts how the outcomes change as the number of training patients grows. In the instance of the 2D UNet, it is clear that vessel detection improves while false positives grow. The 3D UNet, on the other hand, starts with sliced and basic segmentations and ends up recognising all of the major vessels while keeping a cleaner segmentation. 75488 VOLUME 11, 2023 B. Serrano-Antón et al.: Coronary Artery Segmentation FIGURE 4. Example of predicted segmentations carried out by 2D UNet pre-trained and 3D UNet for test patient T001 on aorta and coronary arteries (A+C.A) and coronary arteries only (C.A). Datasets were trained with different number of patients (N). The validation set consisted of 20% of N. The last row shows the results obtained after post-processing the image using the IG algorithm (first the small islands are removed and then grow-shrink is applied). VOLUME 11, 2023 75489 B. Serrano-Antón et al.: Coronary Artery Segmentation FIGURE 5. Parameter values of 2D UNet pre-trained and 3D UNet results on aorta and coronary arteries (A+C.A) and coronary arteries only (C.A) datasets. The X-axis shows the number of patients used for the training set (N). On the Y-axis, the value of the corresponding parameter. A) F1score. B) Fa 1score. C) Recall. D) Recall in coronary arteries. E) Precision. F) False positive to background class number of pixels ratio. G) Number of connected components in coronary arteries. 75490 VOLUME 11, 2023 B. Serrano-Antón et al.: Coronary Artery Segmentation FIGURE 6. Segmentations of the coronary tree of Lesion 1. In red the manual segmentation performed by an expert. (A) In blue the prediction performed by the 2D pretrained UNet, N=65. (B) In green the prediction of the 3D UNet, N=65. The zoomed-in part shows the region of the lesion. FIGURE 7. Segmentations of the coronary tree of Lesion 2. In red the manual segmentation performed by an expert. (A) In blue the prediction performed by the 2D pretrained UNet, N=65. (B) In green the prediction of the 3D UNet, N=65. The zoomed-in part shows the lesion region. In addition, the segmentation for test patient T003 can be seen in Supplementary Material (Fig. S5). Furthermore, the last row shows the effect of post-processing. Structures become cleaner and more connected while maintaining the quality of the previously segmented vessels. Fig. 5shows the parameter values achieved for each type of network and training. Errorbars with mean and standard deviation represent the outcomes of the 10 test patients from 3 trainings with different data (all with the same number of patients, N), for a total of 30 values (see csv file in Supplementary Material). The first parameter to be addressed is the F1score, which is presented in Fig. 5.A. This is an indicator, in the range [0,1], of the similarity between two binary masks, where 0 implies no match between the two masks and 1 a full match. In a binary case, such as ours, it is equivalent to the DSC [17]. One item that strikes out are the numbers for A+C.A training that are around 0.9 and fluctuate slightly with N. This is because the aorta takes up a considerable portion of the volume and is a structure that both networks recognise. In particular, we always get numbers in the range [0.89,0.95]. In the instance of C.A, there is no discernible pattern for the pre-trained 2D UNet, which gives values around 0.7 for all N. This is in contrast to the 3D UNet’s behavior, which exhibits increasing values of the F1 score (from 0.53 to 0.75) as N increases. In Fig. 5.B, the F1score values are shown when only the FPs that are attached to the vessel are taken into account. These results do not show a different pattern with respect to Fig. 5.A, but they do show an increase of up to 0.2. The Recall and Precision metrics are given to provide a more precise picture of the network’s performance. Recall VOLUME 11, 2023 75491 B. Serrano-Antón et al.: Coronary Artery Segmentation FIGURE 8. Segmentation prediction for test patient T002 with networks trained with aorta and coronary arteries (A+C.A). FIGURE 9. Segmentation prediction for test patient T002 with networks trained with coronary arteries only (C.A). tells us how much of the cardiac tree has been recognised by the network. Fig. 5.C displays the outcomes. The same pattern is found for both networks, an increased trend with Nin general and a superior outcome with 2D pretrained UNet. We also notice the ‘‘aorta effect’’ (this refers to the fact that the aorta is a large structure easily recognisable by the network), having recall values greater than 0.8 when training A+C.A. To avoid the ‘‘aorta effect’’ Fig. 5.D displays the Recall values in the coronary arteries regardless of whether the trained structure is A+C.Aor C.A. The positive dependence on Ncan be demonstrated for all eventualities in this case. Mean values for the pre-trained 2D UNet with A+C.Astart at 0.57 when N=15 and rise sharply when N>35, reaching 0.93 when N=65. When N≤35, this identical network trained for C.Aexhibits values that are roughly 10% higher than for A+C.A; for bigger values of N, it follows an upward trend, reaching a mean value of 0.82 when N=65. Furthermore, when N>35, training A+C.Aproduces higher outcomes than C.A. The 3D UNet pattern is similar, albeit the Recall value is lower than in the pre-trained 2D pattern. When training with 75492 VOLUME 11, 2023