 Research article
 Open Access
 Published:
3D multimodal cardiac data reconstruction using angiography and computerized tomographic angiography registration
Journal of Cardiothoracic Surgery volume 10, Article number: 58 (2015)
Abstract
Background
Computerized tomographic angiography (3D data representing the coronary arteries) and Xray angiography (2D Xray image sequences providing information about coronary arteries and their stenosis) are standard and popular assessment tools utilized for medical diagnosis of coronary artery diseases. At present, the results of both modalities are individually analyzed by specialists and it is difficult for them to mentally connect the details of these two techniques. The aim of this work is to assist medical diagnosis by providing specialists with the relationship between computerized tomographic angiography and Xray angiography.
Methods
In this study, coronary arteries from two modalities are registered in order to create a 3D reconstruction of the stenosis position. The proposed method starts with coronary artery segmentation and labeling for both modalities. Then, stenosis and relevant labeled artery in Xray angiography image are marked by a specialist. Proper control points for the marked artery in both modalities are automatically detected and normalized. Then, a geometrical transformation function is computed using these control points. Finally, this function is utilized to register the marked artery from the Xray angiography image on the computerized tomographic angiography and get the 3D position of the stenosis lesion.
Results
The result is a 3D informative model consisting of stenosis and coronary arteries’ information from the Xray angiography and computerized tomographic angiography modalities. The results of the proposed method for coronary artery segmentation, labeling and 3D reconstruction are evaluated and validated on the dataset containing both modalities.
Conclusions
The advantage of this method is to aid specialists to determine a visual relationship between the correspondent coronary arteries from two modalities and also set up a connection between stenosis points from an Xray angiography along with their 3D positions on the coronary arteries from computerized tomographic angiography. Moreover, another benefit of this work is that the medical acquisition standards remain unchanged, which means that no calibration in the acquisition devices is required. It can be applied on most computerized tomographic angiography and angiography devices.
Background
Preliminary consideration
Among the cardiovascular system diseases, Coronary Artery Disease (CAD) is an important issue, which usually stands behind the loss of life around the world today. In fact, CAD is associated with blockage as well as narrowing the left or perhaps right coronary artery vessels. Therefore, a precise method to visualize coronary arteries is highly needed. There are several medical imaging techniques, which can be used for diagnosing heart diseases; Xray Angiography, Cardiac Computerized Tomography Angiography (CTA), Magnetic Resonance Angiography (MRA), Cardiac Positron Emission Tomography (PET), SinglePhoton Emission Computed Tomography (SPECT), Echocardiography [1] and so forth. Among these modalities, Cardiac Computerized Tomography Angiography (CT Angiography or briefly CTA), and also Xray Angiography (Xray arteriography or briefly angiography) are the best ways for visualizing coronary arteries.
Generally, the CT scan is a form of Xray, which utilizes a computer system to generate crosssectional images of the human body. The CTA is a type of medical exam that mixes a CT scan with the injection of a specific dye, referred to as a contrast material, to generate images of vessels in specific parts of the body. For this aim, the contrast material is often injected into a vein started in the hand or arm. When CTA is done, a series of images will be created, which can be observed as an axial view of cardiac components, such as cardiac chambers, aorta, heart’s muscle and coronary arteries as well. These images are referred to as CTA slices in this paper. One can get the 3D reconstruction of coronary arteries and also find out whether a plaque buildup has narrowed patient’s vessels or not. Second, modality is Angiography, which can be well used to visualize the coronary arteries and diagnose the blockage and stenosis in realtime. Therefore, most physicians prefer to use this modality instead of others to diagnose and treat cardiac coronary artery diseases. Coronary artery angiography is performed by injecting the radioopaque contrast agent into the coronary arteries and imaging using Xray based techniques such as fluoroscopy. A series of blood vessels radiographs is called angiograms (or angiographs). A sample of an angiogram and a CTA slice are shown in Figure 1.
Medical context
Cardiovascular system anatomy
Coronary arteries are a vital part of the cardiovascular system, which mostly relies on the surface of the heart and transports the blood to the heart muscle. As shown in Figure 2, coronary artery vessels are categorized into two principal parts: Right Coronary Artery (RCA) and Left Coronary Artery (LCA). RCA stems from the right aortic sinus and LCA stems from the left aortic sinus. The first part of the left coronary artery is referred to as the left main coronary artery (LM). This blood vessel can be about 5 mm wide and less than 30 mm long. LM branches directly into a pair of arteries: Left circumflex coronary artery (LCX) and left anterior descending coronary artery (LAD). LCX circles around the left side of the heart, which is embedded throughout the surface of the rear side of the heart. LAD is embedded throughout the surface of the front side of the heart. Each LCX and LAD artery bifurcates into smaller sub arteries; three septal arteries (S1 ~ S3) and three diagonal arteries (D1 ~ D3), which originate from LAD, and also two marginal arteries (OM1 and OM2), which come from LCX. The end of RCA bifurcates into a pair of smaller arteries: right posterior lateral branch (RPLB) and right posterior descending artery (RPDA).
Medical problem
Different modalities hold different information in assisting the physicians in making decision. Previously, they viewed these images inseparate or adjacent windows, and mentally fused these images together. Even though medical image registration and image fusion have been successfully implemented in other organs such as the brain and lungs, medical image registration and fusion for the heart present different challenges. First, since the heart is beating, fusion requires synchronization with the rhythm of the heart for different phases. Second, the heart is a nonsolid organ, thus acquired images give a vague impression.
Precise computerassisted coronary artery analysis is consequential in later diagnosis and treatment with cardiologists and cardiac surgeons. As explained in Section 1.1, angiography and CTA are the best modalities to visualize coronary arteries and CAD diagnosis. Individually, these modalities provide valuable information, but do not represent the complete information about coronary arteries. Thus, it is important for specialists (cardiologist or cardiac surgeon) to combine informative data in both modalities. Indeed, specialists believe that information about stenosis points in coronary arteries from angiogram, such as the 3D position, and their relative localization with respect to the corresponding segmented coronary artery from CTA is an important aid in CAD diagnosis.
There are few hybrid/combined devices like SPECT/CT, PET/CT, PET/MRI and MRI/PET, which give combined informative data. To the best of our knowledge, there no hybrid device is available for combining the results of angiography and CTA devices. Also, these devices are very expensive and none gives back the 3D position of stenosis lesion of the coronary arteries. In a successful manner, these two modalities are available in most hospitals nowadays. Therefore, one way to aid the specialist is to set up the connection between stenosis points from an angiography along with their 3D positions on the coronary arteries from CTA, so that they have a 3D informative model consisting of stenosis and coronary arteries’ information from both modalities.
Previous work
This work introduces a new method for 3D reconstruction of stenosis point on coronary arteries through the registration of both CTA and angiography modalities. Three main phases were defined for this method: (1) Coronary artery labeling and segmentation in angiography. (2) Coronary artery segmentation and labeling in CTA, and (3) Registration of CTA and angiography images. It is worth noting that for the registration phase in this work, it is highly required for the coronary arteries to be segmented, labeled and also, their control points are extracted precisely from both modalities. The previous works of each above phases are explained as follows.
Coronary artery segmentation and labeling are crucial steps, because further processes such as 3D reconstruction, fusion with other modalities, stenosis measurement, and blood flow analysis use the result of these steps. There are some difficulties in coronary artery segmentation from angiograms; such as poor signal to noise ratio and artefacts, which caused by organs such as the backbone and ribs. Several methods have been proposed to segment coronary arteries in angiograms [2,3], and [4]. However, they only segmented coronary arteries without labelling them. Xu et al. in [5] proposed an algorithm for coronary artery centreline tracking in angiograms using matched filter on the eigenvalues. In another study, HernándezVela et al. [6] introduced an accurate coronary centreline extraction in angiograms. Several semiautomatic algorithms were proposed relating to the coronary artery segmentation. One of the limitations of the proposed algorithms is that they involve users in defining seed points to locate the coronary arteries. For instance, Wang et al. [7] proposed a method for coronary artery segmentation in angiograms, but it requires a seed point to start. Meanwhile, other algorithms suffer from high computational complexity. For instance, Zhou et al. [8] proposed an automatic approach for segmenting coronary arteries in angiogram. One of the limitations of their work is the time complexity, and also nonprecise arteries extraction. Some other algorithms have proposed segmentation for general blood vessels. However, adopting these algorithms on the coronary artery in angiograms may lead to the appearance of some artefacts in the result. For instance, Li et al. [9] proposed a regionbased active contour model for vessel segmentation. Running their algorithms on angiograms would lead to inefficient artery segmentation, whereby parts of the background might appear, while parts of the arteries disappear. Also, none of the above method labels coronary arteries based on the angiogram segmentation result. The limitations of previous algorithm motivated us to propose a new method for coronary artery segmentation and labelling, which visualizes and labels all main coronary arteries from angiography images.
Coronary artery segmentation using CTA and their 3D reconstruction should be done precisely because registration using images from angiography modality is done based on this result. However, there are plenty of problems in coronary artery segmentation and labelling from CTA slices. Firstly, coronary arteries are shown as small parts, semicircular or tubular shape in each slice. Therefore, tracking via slices is not a straightforward process. Secondly, artefacts from other bodily organs in CTA slices, such as the backbone, ribs, cardiac chambers and other components should be removed from the final segmented image. Several algorithms have been proposed regarding coronary artery segmentation and 3D reconstruction from CTA slices [10]. Most of them used the 3D Frangi’s algorithm for coronary artery segmentation from CTA [11,12]. However, these algorithms suffer from artefacts, such as false step edge responses in some parts of the coronary arteries, especially near the right atrium. Yang et al. [13] solved this problem by proposing an improved 3D Frangi’s method. They improved the 3D Frangi’s vesselness filter by adding local geometrical features. Another drawback of the current methods is their noncapability of labelling, especially at the same time with the segmentation phase. The only study on coronary artery labelling from CTA was proposed in [14], where arteries are labeled after segmentation in a different phase. They identified the main branches using pointset registration method as proposed in [15]. However, the registration phase of their algorithm suffers from high computational complexity and also, the control points required for registration with angiogram in our work are not extracted. These mentioned limitations motivated us to propose a new method for coronary artery segmentation, labelling and also 3D reconstruction of CTA slices.
As explained before, the aim of this work is 3D reconstruction of stenosis point through CTA and angiography image registration. Therefore, registering both CTA and angiography modalities is one of the important phases of this work. Several works have been done for coronary artery registration in these two modalities. Some of them were done for 2D/3D coronary artery registration in both modalities [1618], and some proposed 2D/2D registration to guide endovascular stent grafting [19,20], and [21]. Also, another work on coronary artery registration was done in [22], but specific angiography devices (biplane) were needed in their work. Furthermore, all of the above studies proposed registration algorithms to align coronary arteries in both modalities CTA and angiography. The benefit of these works is the usage of CTA result (preinterventional/preoperative) as the image guidance in angiography (interventional/intraoperative) procedure for percutaneous coronary intervention (PCI). However, these methods cannot be applied in our work because they only aligned some correspondent arteries together, and therefore, it is not possible to use them for registering the start, bifurcations and end points of the correspondent coronary arteries from two modalities. The benefit of having these points is it enables the estimation of stenosis point between the nearest two points from one artery in angiogram to the correspondent artery in CTA. The above mentioned limitations motivated us to propose a new method for featurebased coronary artery registration using proper control points.
Methods
The keyideas of the proposed method are: (1) No calibration process required for the input images. (2) No need for 3D reconstruction images from 2D angiography. (3) No need for multiple views of angiography. (4) No need for specific angiography (like biplane) devices. Figure 3 provides a diagram outlining all steps of the proposed method.
Some steps in the proposed method can be described as follows.

The specialist chooses an angiogram such that the stenosis lesion can be observed clearly.

2D main coronary artery segmentation and labeling for the chosen angiogram are done in 2D (u^{'}, v^{'}) plane.

3D main coronary artery segmentation and labeling are done for CTA modality in (x,y,z) Cartesian coordinate system.

Candidate control points are detected automatically for both modalities.

Stenosis lesion is marked in angiogram by the specialist. In addition, the coronary artery with the stenosis is chosen for registration.

CTA is rotated to each direction (x,y,z) to get an almost similar view of the artery as displayed in the angiogram.

Correspondent control points are marked by the specialist in both angiogram and CTA modalities.

Rotated artery and their control points in CTA are projected onto 2D (u, v) plane.

The selected control points are normalized for both modalities, individually.

Affine transformation function is computed for the normalized control points.

The correspondent coronary arteries from both modalities in (u^{'}, v^{'}) and (u, v) planes are registered.

Stenosis lesion from the angiogram is detected on the correspondent artery from CTA.

Stenosis point is back projected in the 3D Cartesian coordinate system (x,y,z).
The crucial steps of the proposed method are: (1) Coronary artery segmentation and labelling in angiogram. (2) Coronary artery segmentation and labelling in CTA. (3) Featurebased registration of coronary arteries in CTA and angiogram. These steps are respectively detailed as follows.
Coronary artery segmentation and labelling in angiogram
In this section, a new method for coronary artery labelling and segmentation from angiogram is proposed. The methodology of this work includes five phases, as described in the following subsections. In the first phase, after removing noise from the raw angiograms with Discrete Wavelet Transform (DWT), all arteries are sharpened with Starlet Wavelet Transform (SWT). In the next step, the main coronary arteries are segmented by applying the modified SWT. After that, coronary artery centrelines are segmented, and also detached from each other. Then, all centrelines are labeled, and finally, coronary arteries are labelled by constructing a proper mask.
Angiogram preprocessing
Each type of Xray angiography has a moving DICOM format result. Dealing with this type of imaging is too difficult. First, the angiograms have to be converted into n 2D bitmap frames {f_{1}, f_{2}, … f_{ n }} on each DICOM movie angle, such that n varies from 40 to 70 frames. One optimum frame (f_{ optimum }) must be chosen from these frames. The optimum frame is defined as a well Xray injected frame that the whole coronary arterial tree is contrasted with that, and also the stenosis can be detected by specialist. After choosing an appropriate frame, it is converted from RGB format to grayscale, which we call as the original input image \( {I}_A^o \) for next steps.
Since angiograms suffer from a large amount of noise, the first step is noise removal. The aim of the noise removal process is to eliminate all noise while preserving the quality of the images. Here, quality hints to retain the arteries in the angiograms. When we used the traditional algorithms for removing noise in angiograms, such as the smoothing method, some arteries disappeared. Hence, another technique for removing noise should be defined by converting images into a transformation domain, such as wavelet, and then compared the transformation’s coefficients to a proper threshold value. In this way, the arteries’ structure is kept. Therefore, DWT with wavelet type Haar and level 3 was used for removing noise from angiograms in this work. The details of using DWT for noise removal are well defined in [3,23]. This step is shown in Eq (1).
The next step in preprocessing is coronary artery sharpening. In this step, we intended to sharpen the edges of arteries and erased the surrounding background using SWT, which has been defined in [3]. Because this work concentrated on the main coronary artery labelling, this method was modified. SWT, or Isotropic Undecimated Wavelet Transform (IUWT), is well known in the field of biology [24], astronomy [25] and nowadays in medical applications [3] and [26]. SWT decomposes an image I into w_{ j }, as a wavelet coefficient and c_{ j } as a scale coefficient in each iteration j [27]. A preliminary algorithm for SWT [28] is presented in Figure 4.
In this section, the SWT algorithm is modified to obtain sharper coronary arteries in angiograms, and in the next step, it is changed again to segment the main coronary arteries from angiograms. For segmenting the main coronary arteries from angiograms, different filters were applied with various wavelet levels (l) and finally, we found that using h^{0} = [1,3,3,1]/8 on the input image with selective wavelet levels l = {1,2,3,4,5} were best for sharpening. The modified algorithm for sharpening main coronary arteries in angiogram is illustrated in Figure 5.
Main coronary artery segmentation
SWT is a type of wavelet transform, which is fast in computation and can be used for object segmentation [3]. In this phase, a new method for main coronary artery segmentation is proposed based on SWT application. Since intensity values of objects must be higher than background when segmenting with SWT, \( {I}_A^s \) should be inverted first. The output is called \( {I}_A^i \). Then, the modified SWT is applied on \( {I}_A^i \) for coronary artery segmentation. Generally, applying various filters and also several wavelet levels (l) in SWT returns different results. In the case of angiogram, thinner arteries are shown in smaller l values and thicker ones in higher l values. In addition, using various filters returns different qualities of coronary arteries. Empirically, it was found that using wavelet levels l = {2,3,4,5} applying filter h^{0} = [1,3,3,1]/8 on \( {I}_A^i \) was best for main coronary artery segmentation. The proposed algorithm is shown in Figure 6.
After the main coronary artery segmentation, a postprocessing step was considered, which includes thresholding, followed by length refinement and filling holes, as discussed in [3]. The final result is shown in Figure 7.
Centreline extraction and detachment
In this phase, the centrelines of coronary arteries are extracted using morphological operation on the results of the previous phase \( {I}_A^e \). This process removes pixels, so that all arteries are thinned to a minimally connected stroke. Then, arteries are detached. To do this, the start, the branch and the end points of arteries are detected by counting pixels’ neighbours. This process is done by convolving the centreline of the coronary arteries with a 3 by 3 kernel of ones and the results are saved in matrix \( {I}_A^b \). Then, centreline detachment is achieved after removing the branch points. In addition, short segmented arteries are removed by counting pixels in each detached centreline arteries and by using a proper threshold. The result of centreline extraction is saved in \( {I}_A^c \) and all starting, branching and ending points of arteries are saved separately in matrix \( {I}_A^p \). These steps are shown in Figure 8 for one selected angiogram.
Coronary artery centreline labeling
Before centreline labeling, all branch points are cleared. This process is done by computing the Euclidean distance to remove some centreline pixels in \( {I}_A^c \), which are closer to the branch points than to the nonvessel points. Then, centrelines are labelled by saving any connected centreline pixels with a different number in a new matrix \( {I}_A^{cl} \), which has the same size as the original angiogram. The result of centreline labelling is illustrated in Figure 9.
Coronary artery labeling
As explained before, all centrelines of the coronary arteries are detached and saved with unique numbers in matrix \( {I}_A^{cl} \). The aim of this phase is to label the coronary arteries from the labeled centrelines. First, the artery is selected by choosing the correspondent numbers in the labeled centrelines image \( {I}_A^{cl} \). The output is called \( {I}_A^{cs} \). Then, a proper mask is constructed to segment the correspondent artery in \( {I}_A^e \). Mask construction is done by computing the Euclidean distance transforms between all pixels on the selected labeled centrelines \( {I}_A^{cs} \) and the boundary of the correspondent artery in \( {I}_A^e \). Then, the maximum value of these distances is chosen as the mask radius r. After that, all pixels on the selected centrelines in \( {I}_A^{cs} \) are inflated with radius r. to achieve the proper mask. This mask is applied on \( {I}_A^e \) to construct the labelled coronary artery and the result is saved in I_{ A }. The start, the end and the branch points on the selected artery are obtained using Eq (2).
P_{ A } is the candidate control points in the angiogram and is used for the registration part of this work. The result of the main coronary artery labelling for one sample artery is shown in Figure 10. The candidate control points P_{ A } are also displayed on the start, the branch and the end points of the labeled artery.
Coronary artery segmentation, labelling and 3D reconstruction from CTA
In this part, a new method for coronary artery segmentation, labelling and 3D reconstruction from CTA slices proposed. Figure 11 displays this method, which consists of three main phases: 1) mask construction and aorta segmentation, 2) coronary artery enhancement, 3) coronary artery segmentation and labelling.
The details of every phase are described in the next subsections.
Mask construction and aorta segmentation
In the first phase, CTA slices are classified anatomically. Normally, every CTA slice is decomposed into these main components: pulmonary tissues, pulmonary vessels, aorta, diaphragm, pericardium, myocardium, aorta, bones, ventricles, hepatica tissues and (contrasted/uncontrasted) coronary arteries [29]. To distinguish these components, Hounsfield Unit (HU) is used, which is computed by voxel intensity. It is given by:
where Slope and Intercept parameters are obtained from the DICOM info of every CTA device. In this work, the histograms of 12 CTAs were assessed and it was found that most of them had four main regions through “1024” to “+1000” (HU) values, which specify different types of heart’s components. As result of this experiment, it is shown HU range for every component is different. For some such as contrasted coronary arteries, bone, aorta and ventricles, it is intrinsically high (between “+391” and “+1000” in region 4), while for others such as pulmonary tissue, it is markedly low because of the air (between “1024” to “225” in region 1). For particular soft tissues such as the diaphragm and pericardium tissues, the HU range is regarded as “226” to “+30” as considered in region 2, and for uncontrasted coronary arteries, myocardium and hepatica tissues, it is between “+31” and “+390”, which fits in region 3. The result of the categorization is illustrated in Figure 12.
Since the most important component in this work is the coronary artery, first, we categorized each CTA slice as it has the most probably coronary arteries through various other heart components. Therefore, we focused on regions 3 and 4 because the coronary arteries (both contrasted and uncontrasted) are shown in these two regions. Due to the fact that pulmonary vessels and tissues are located near to the coronary arteries in some slices, especially bordering the heart, and could be considered as coronary arteries incorrectly in some slices, we constructed a mask from region 1 to remove them in all slices. First, region 1 is extracted using HU. values for every CTA slice and then, the specific threshold is considered to create binary images. After that, the filling hole method is applied to remove the holes coming from the pulmonary vessels. Finally, the mask is constructed by reversing the black and white colors. These steps are shown in Figure 13 for a selected CTA slice.
Prior to coronary artery enhancement, aorta segmentation is done using Hough circle transform in initial axial slices. First, region 4 is selected based on HU measurement, small objects such as noise are removed. Then, the Hough circle transform is applied to detect circular objects and also some postprocessing methods are employed to remove residual noise and fill small holes. It is shown for a selected slice in Figure 14. Finally, the aorta is saved in a 3D matrix named M_{ Aorta } for all initial slices by considering the correspondent slice numbers as the rows for the 3D matrix.
Coronary artery enhancement
In the second phase, Frangi multiscale filter is applied to measure the vesselness components based on the eigenvalues in the Hessian matrix [30]. This method is used to find and enhance tubular components by simply computing the second order derivatives in the Gaussian kernel at various scales and giving a value between 0 and 1 for each pixel x at certain scale σ. This vesselness measure function is formulated in Eq (4).
where \( {\mathrm{\mathcal{R}}}_B=\frac{\left{\lambda}_1\right}{\left{\lambda}_2\right} \)., \( s = \sqrt{{\lambda_1}^2+{\lambda_2}^2} \). and λ_{1}, λ_{2} are eigenvalues in a 2D Hessian matrix. φ_{1}. and φ_{2} determine the level of sensitivity in the filter to the amounts ℛ_{ B } and S, respectively.
The 2D Hessian matrix for a given pixel x and scale σ is also given by:
where I_{ αβ }(x) indicates the second order derivative of the input image at pixel x obtained by convolving the image using the 2D Gaussian kernel G(,s) at scale σ. These definitions are formulate in Eq (6) and (7).
After applying the mentioned filter, the coronary arteries are enhanced as shown in Figure 15.
Coronary artery segmentation and labeling from CTA
As shown in Figure 15, there are many components that should be removed to obtain only the coronary arteries. For this, the Intersection Tracking method is proposed for tracking coronary arteries through 2D slices from ostium to the end. Before discussing the proposed method, the “Cardiovascular tree model” is defined based on Figure 2. As shown in Figure 16, LCA and RCA have a treelike structure and can be categorized as the left coronary arterial tree and right coronary arterial tree. We called this structure as the “Cardiovascular tree model”, which is used in this work as a prior knowledge for coronary artery labelling.
Before explaining the proposed algorithm, the Intersection Tracking method is explained, which is proposed for coronary artery segmentation and labelling from CTA. This method works based on the fact that each coronary artery has continuous pixels through slices from the start to the end point. Therefore, each part of the artery in every slice has intersection with the prior and the next parts in the previous and the next slices, respectively. See Figure 17 for LAD artery in three sequential slices.
As shown in the figure, the selected part of LAD has intersecting pixels with another part of LAD in the previous and next slices. Therefore, if it is possible to segment any part of the artery in specific slice i, other parts can be tracked and accessed through the previous and next slices. Based on this fact, and also the “Cardiovascular tree model” in Figure 16, the Intersection Tracking method is defined as follows:

Step 1: Start from a seed point preserved in slice S_{ i }; i is the number of slices and can be defined based on the artery considered for segmentation. We will discuss about the seed point for each artery later in this section.

Step 2: Segment the correspondent region and preserve it at the same position in a 2D temporary matrix called H. The size of matrix H is considered to be the same as slice S_{ i }.

Step 3: Construct a 3D matrix called M, and allocate the ith row from top to bottom by the 2D matrix H. To preserve the start, the bifurcation and the end points of the artery, at the first centroid point c_{ i } of the segmented region in matrix H should be computed. Then, a new 3D matrix called M^{p} is constructed, and the value “1” is allocated to the ith row from top to bottom and 2D position c_{ i } of this matrix. This point is later used as the start point of the artery in the registration phase. Indeed, the 3D matrix M^{p} preserves the candidate control points of the coronary arteries (such as, start, bifurcation and end points)

Step 4: Proceed to the next axial slice S_{i + 1} by increasing the counter i. Then, intersect it with the previous segmented region, which is currently saved in matrix H, to find the new segmented region. Plenty of segmented components are in slice S_{i + 1}. Based on the Intersection Tracking method, the component that has intersection with matrix H is selected. If there is no intersection between matrix H and slice S_{i + 1}, go to step 9.

Step 5: Replace matrix H with the new segmented region in slice S_{i + 1}, and allocate it in the next (i + 1)^{th} row in the 3D matrix M.

Step 6: Bifurcation detection is done in this step. As shown in the “Cardiovascular tree model”, each main artery such as LM, LCX, LAD and RCA has some subarteries. Therefore, it should be examined if the artery is bifurcated in the current slice or not. If it is, the subartery is removed to get only the main coronary artery. The thinning process is applied on the current segmented region, which is currently saved in matrix H, using morphological operation and convolution with the 3 by 3 kernel of ones to check the mentioned condition. If this region includes a bifurcation, the algorithm proceeds to the next step, and value “1” is allocated for the ith row from top to bottom of 3D the matrix M^{p} as a bifurcation position. Else, return to step 4 for the next slice.

Step 7: For discovering the continuing branch from the main artery and the subartery, every branched vessels are segmented and saved in two new temporary 2D matrices, N^{'} and N^{' '} (at the same position in the corresponding slice S_{i + 1}). Also, two 3D temporary empty matrices, M_{1} and M_{2} are constructed and 2D matrices, N^{'} and N^{' '} are added on (i + 1)^{th} row of each of them, respectively. Then next slice is considered by increasing the counter i. for each branched vessel regions in parallel, which is currently saved in two temporary 2D matrices N^{'} and N^{' '}, and they are tracked in the following slices using the mentioned Intersection Tracking method. This procedure continues and new segmented regions N^{'} and N^{' '} are added into the two 3D temporary matrices M_{1} and M_{2}, until one of conditions appear:

1)
One of vessels reaches the end: As shown in Figure 18, the ended vessel is a subartery and should be removed. The correspondent 3D matrix is removed and the remaining temporary 3D matrix (M_{1} or M_{2} based on which one is considered as continuous pixels of the main artery) is merged into a main 3D matrix M. Matrix H is also replaced with the last region in N^{'} or N^{' '} (based on the selected part, which is considered as the main artery).

2)
One of the vessels bifurcates before another one ends: If one of the arteries bifurcates again before another one ends, it means that most probably, this artery is as the main artery and the other one should be removed. But as shown in Figure 18, sometimes, the subarteries, such as septal or diagonal, make a subtree by bifurcating themselves. Therefore, whenever an artery is bifurcated, new 3D and 2D matrices as mentioned above are constructed for new branches. The same procedure is done for all branches in parallel to find the main artery and remove the other ones. The bifurcation points are reserved in the proper position in 3D matrix M^{p}, as explained before.

1)

Step 8: Go to step 4 for the next slice.

Step 9: The centroid point c_{ i }. of the segmented region in matrix H is computed. Then, the value “1” is allocated to the ith row from top to bottom of 3D matrix M^{p} and 2D position c_{ i }. This point is used as the end point of the coronary artery.

Step 10: Save the 3D matrix M as M_{ LM }, M_{ RCA }, M_{ LAD } and M_{ LCX }, based on the coronary artery considered for segmentation in this phase. The candidate (start, bifurcation and end) points for each main coronary arteries are saved in one of the 3D matrices \( {M}_{LM}^p, \)\( {M}_{RCA}^p, \)\( {M}_{LAD}^p \) and \( {M}_{LCX}^p, \) as well.
As mentioned in step 1, the Intersection Tracking algorithm starts from a seed point for each LM, LAD, LCX and RCA arteries. Different conditions are considered for choosing them. As shown in the “Cardiovascular tree model” in Figure 16, LCA is started from LM artery. Axial ices are sought from top to bottom to find slice S_{ i } as LM’s seed point, where LCA starts from the aorta and i is the number of slice. Then, the above Intersection Tracking algorithm is applied on LM segmentation. When a bifurcation is detected in step 6, it means that the LM ends in that slice and both LAD and LCX would start from that point. Therefore, the bifurcation point is selected as a seed point for LAD and LCX segmentations, as shown in Figure 19.
For RCA, axial slices are sought from top to bottom to find slice S_{ i } as RCA’s seed point, where the RCA starts from the aorta. Then, the same Intersection Tracking method is used for RCA segmentation. As shown in the “Cardiovascular tree model”, RCA has one main branch with some subarteries and also two main branches, called RPDA and RPLB. The sub arteries should be removed, except for RPDA and RPLB because these two subarteries are considered as the main parts of RCA, and sometimes, stenosis can appear in these parts. A proper threshold is defined in step 7 of the above algorithm for preserving these two subarteries. The steps for LAD and LCX segmentation in few sequential slices are shown in Figure 20.
The pseudo code of the proposed method in this phase is described in the following.
Featurebased CTA and angiogram registration
As analyzed in other papers dealing with cardiac image registration [3133], there is no previous work on coronary artery registration from CTA and angiography modalities to create a 3D reconstruction of stenosis point. In this section, the defined control points and a proper transformation function are used to register the coronary arteries from both CTA and angiography modalities.
First, different angiograms from several angles are examined by a specialist (cardiologist or cardiac surgeon) to find the best view, in which stenosis can be visualized and detected. Then, the coronary artery (which has stenosis on it) is segmented and labeled from the selected angiogram using the proposed algorithm in section 2.1, and the candidate control points are detected. Then, the position of the stenosis point is marked by the specialist. After that, the correspondent artery and their candidate control points are segmented from CTA using the proposed algorithm in section 2.2. Since the CTA results are preserved as 3D points in one of the 3D matrices M_{ LM }, M_{ RCA }, M_{ LAD } or M_{ LCX } (based on which one is detected as having the stenosis), it is possible to rotate the artery so that it is in the same angle as the selected angiogram. Then, the rotated CTA is projected from (x, y, z) Cartesian coordinate (3D) to the (u, v) plane (2D). The same projection onto the (u, v) plane is done for the candidate control points in 3D matrices \( {M}_{LM}^p, \)\( {M}_{RCA}^p, \)\( {M}_{LAD}^p \) and \( {M}_{LCX}^p, \) based on the selected artery. The result of the candidate control points projection from CTA is called P_{ C }.
Let us consider I_{ C } and I_{ A } as two 2D images containing the segmented correspondent arteries from the projected CTA and the angiogram, respectively. Also, C and C^{'} denote the correspondent control points selected from P_{ C } and P_{ A } (defined in Eq (2)), respectively.
By considering i as a point number, C_{ i } and \( {C}_i^{\prime } \) denote the correspondent control points, which is mathematically formulated in Eq (9).
In the first step of the registration, C_{ i } and \( {C}_i^{\prime } \) are normalized to improve the accuracy of the result. The normalization comprises of translating and scaling the coordinates of I_{ C } and I_{ A } images, defined as the following steps:
1. First, the control points’ coordinates inside every image are separately translated in order to bring the centroid toward the origin, formulated as follows:
In which \( \overline{u}=\frac{{\displaystyle {\sum}_{i=1}^n}{u}_i}{n}, \)\( \overline{v}=\frac{{\displaystyle {\sum}_{i=1}^n}{v}_i}{n}, \)\( {\overline{u}}^{\hbox{'}}=\frac{{\displaystyle {\sum}_{i=1}^n}{u_i}^{\hbox{'}}}{n}, \)\( {\overline{v}}^{\hbox{'}}=\frac{{\displaystyle {\sum}_{i=1}^n}{v_i}^{\hbox{'}}}{n}. \)
And the centroid points for each modality are defined as: \( c=\left[\begin{array}{cc}\hfill \overline{u}\hfill & \hfill \overline{v}\hfill \end{array}\right] \) and \( {c}^{\hbox{'}}=\left[\begin{array}{cc}\hfill {\overline{u}}^{\hbox{'}}\hfill & \hfill {\overline{v}}^{\hbox{'}}\hfill \end{array}\right] \)
In the next step, C_{ T } and \( {C}_T^{\hbox{'}} \) are scaled so that the root mean squared (RMS) distance from the origin is equivalent to \( \sqrt{2} \). It is formulated for C_{ T } as:
The RMS, root mean squared, distance is given by Eq (12).
By using Eq (11), the scale factor S for the control points in I_{ C } can be calculated as:
In the same way, the scale factor S′ for the control points in angiogram I_{ A } can be calculated as:
a result, new two matrices containing normalized control points are computed using Eq (15) and (16).
And the new correspondent control points are shown as:
After normalization, the geometrical transformation function is computed using the normalized control points. We used the affine transformation matrix for this function, which has the following format:
When there are at least three normalized control points, the affine transformation function can be written as Eq (19).
where \( \tilde{C} \) and \( {\tilde{C}}^{\prime } \) are defined in Eq (15) and (16), respectively.
This equation can be simplified as:
where
Therefore, the affine transformation matrix of the normalized control points can be computed by Eq (22).
By calculating the above equation, all variables in \( {\tilde{T}}^{\prime } \) are computed and \( \tilde{T} \) are obtained.
Eventually, a proper correction should be done to bring the (normalized) transformation function \( \tilde{T} \) to the original coordinate system. This step is called denormalization and mathematically formulated in Eq (23).
with
By applying the transformation function T, the coronary artery from angiogram I_{ A } is registered on the correspondent one in I_{ C }. The registered version of angiogram \( {I}_A^{\hbox{'}} \) is formulated in Eq (25).
As both arteries have been registered, they can be considered to be in the same plane, called image I_{ R } in (u, v) plane, which contains both \( {I}_A^{\prime } \) and I_{ C }.
After registration, the stenosis point from angiogram should be located on the correspondent artery in CTA.ost probably, the stenosis point from angiogram intersects with the coronary artery from CTA in I_{ R }. Otherwise, the stenosis point is estimated as follows. First, a proper thinning algorithm is used to extract the centrelines of both arteries in the registered plane, called \( {I}_R^{\hbox{'}} \). Then, the points p_{ k } and p_{ l } are detected as the two nearest points to the stenosis point in \( {I}_R^{\hbox{'}}, \) such that the two arteries intersect at those points. Then, the length between the stenosis point and p_{ k }, and also the length between stenosis point and p_{ l }, for the artery from angiogram in \( {I}_R^{\hbox{'}} \) are computed, called \( {l}_1^{\hbox{'}} \) and \( {l}_2^{\hbox{'}}, \) respectively. In addition, the length of the centreline of the artery between two points p_{ k } and p_{ l } for the artery from CTA in \( {I}_R^{\hbox{'}} \) is computed, called L. By calculating the ratio between the above lengths, the stenosis point is localized in the correspondent artery from CTA. It is formulated in Eq (26).
After locating the stenosis point on I_{ C }, it is reconstructed in 3D. Since the 3D Cartesian coordinates (x, y, z) r each point in I_{ C } are saved in the original 3D matrix of that artery, the stenosis point can be obtained by back projecting that point to the original 3D matrix. Some of the above steps are illustrated in Figure 21.
Results and discussion
In this section, first the optimum values for parameters are reported and afterwards, the efficiency of each algorithm is evaluated. We obtained the dataset of angiogram and CT Angiography from the UiTM Medical Center. Five patients’ angiograms acquired with the PHILIPS Angiography device are used for the system configuration (training dataset) and 12 patients’ angiograms are used for the system analysis (test dataset), to evaluate the performance of the proposed method in this modality. The angiogram dataset includes 493angiograms of 12 different patients acquired in eight angles for each of patients. The format of them was 24bit grayscale BMP in 512 × 512 pixels size. We randomly selected 60 angiograms (30 LCA and 30 RCA) from this as our dataset for this work, and manually created the ground truth for each. Also, four patients’ CTAs acquired with the SIEMENS 3D CT Angiography device are used in best diastole for the system configuration (training dataset) and 12 patients’ CTAs are used for the system analysis (testing dataset), to evaluate the performance of the proposed method in CTA modality. The images has the size of 512 × 512 pixels in the horizontal plane and about 450 slices of 0.75 mm for the z axis. The corresponding ground truth was created from the slices, as well. It is worth noting that all images in ground truth were validated accurately by the expert cardiac surgeon, cardiologist and radiologist. In this work, we applied the qualitative, and also quantitative methods such as accuracy, precision, sensitivity, specificity and the normalized sum of false detections values, to determine the performance of this work [34,35].
where T_{ p }, T_{ n }, F_{ p }, F_{ n } and q are defined:
T_{ p }, (True positive): the number of coronary artery pixels detected correctly.
T_{ n }, (True negative): the number of noncoronary artery pixels detected correctly.
F_{ p } (False positive): the number of noncoronary artery pixels detected as coronary artery.
F_{ n } (False negative): the number of coronary artery pixels not detected.
q: the quantity of all pixels inside an image.
Setting the parameter values
Prior to accomplishing experiments regarding the overall performance assessment, suitable values in the proposed method’s parameters need to be determined. Therefore, we evaluated the effect of different values for each of the parameters and chose the best in each phase of this work.
Parameters in coronary artery segmentation and labeling in angiogram
For coronary artery segmentation and labelling in angiogram, parameters filters (h^{0}) and wavelet levels (l) parameters of the proposed algorithm were considered based on the accuracy of the proposed method for the main coronary artery segmentation. We examined these parameters on 20 different angiograms from the dataset, and finally calculated the average of each for optimization purposes. To obtain the best result in this method, we applied filters (h^{0}) [1,2,1]/4, [1,3,3,1]/8 and [1,4,6,4,1]/16 on the dataset, and then calculated the accuracy for each. Results are shown in Figure 22.
As shown here, the best result was related to h^{0} = [1,3,3,1]/8.
As explained before, thinner arteries are detected by smaller levels and thicker ones by higher levels. We applied the range of wavelet levels (l ⊆ {1,2,3,4,5}) for this research and examined their effects individually and as a whole. The obtained results are shown in Figure 23.
As demonstrated in Figure 23, the best result was obtained using wavelet levels l = {2,3,4,5}.
Parameters in coronary artery segmentation and labelling in CTA
For mask construction, the slope and intercept parameters in Eq. (3) were obtained from the CTA DICOM info. In our dataset, they typically had “1” and “1024” for slope and intercept, respectively. For coronary artery enhancement, two φ_{1} and φ_{2} parameters in Eq. (4) were obtained based on the dataset characteristics. We examined different value combinations on the vesselness measure function v. We considered the range of 0.25 to 0.75 for φ_{1} and the range of 15 to 35 for φ_{2} parameters. The effect of these combinations on the proposed Intersection Tracking method in a selected slice is shown in Figure 24, by considering a specific range of 1 ≤ σ ≤ 4 as the scale value. The value of φ_{1} controls the shape structure and should be less than 1, and φ_{2} establishes the particular effect associated with contrast robustness on artery enhancement. By choosing larger values of φ_{2}, lowcontrast items were usually disregarded in support of arteries along with considerable increase of contrast. Smaller values for this parameter, for instance, φ_{2} = 15, added more noise around the artery in the result in comparison with larger values. Nevertheless, a large value, for instance, φ_{2} = 35, caused considerable decline of the vesselness result, possibly for higher contrast arteries as shown in Figure 24.
The optimal values for φ_{1} and φ_{2} were determined by comparing the results of the proposed method on the slices of the dataset to the created ground truth. For this, the normalized sum of false detections (_{ F }) was calculated. Here, _{ F } is an objective discrepancy measure that quantifies the change of the artery enhanced images acquired by utilizing various values associated with parameters φ_{1} and φ_{2}, from the corresponding ground truth images. In this study, three values for φ_{1}, φ_{1} = {0.25,0.5,0.75} and three values for φ_{2}, φ_{2} = {15,25,35} were applied on the constructed dataset and the minimum discrepancy measure ε_{ F } was calculated for all combinations. Based on the result of this experiment, the values 0.75 and 25 were selected for parameters φ_{1} and φ_{2}, respectively.
Experimental evaluation
Validation is usually a difficult, but essential phase for any coronary artery segmentation, labeling and registration techniques. With the exception of specific cases, such as centreline extraction, the main way to validate the results of coronary artery segmentation, labeling and registration in both CTA and angiography is by visual observation by an expert cardiologist and radiologist, or by comparison with the manually segmented case as the ground truth. We evaluated the proposed method in every phase of this work, as follows.
Evaluation of main coronary artery segmentation and labeling in angiogram
In this section, we evaluated the capability of the main coronary artery segmentation of the proposed method on the dataset. First, we applied the proposed algorithm on 60 angiograms in the dataset, which included 30 LCA and 30 RCA. Then, we calculated the averages of accuracy, precision, sensitivity and specificity values for each.
As shown in Table 1, the average accuracy, precision, sensitivity and specificity values of the proposed method in LCA angiograms were 0.95153, 0.86812, 0.85897 and 0.97154, respectively. And as shown in Table 2, the average accuracy, precision sensitivity and specificity values of the proposed method in RCA angiograms were 0.96547, 0.92871, 0.89105 and 0.98154, respectively. With an emphasis on the capability of the proposed method, we compared our algorithm with a number of stateoftheart coronary artery segmentation methods on the same dataset. To this end, the proposed methods by Khaleel et al. [4], Li et al. [9], Frangi et al. [30] and Bankhead et al. [26] were used for comparison. The result was summarized in Tables 1 and 2. As shown in these tables, all performance metrics such as, accuracy, precision, sensitivity and specificity for our method were much higher than the others for both LCA and RCA angiograms. Since Bankhead et al. [26] applied SWT method for vessel extraction in retina images, we also used their method here. They used the original SWT algorithm using the filter [1, 4, 6, 4, 1]/16, but as mentioned before, we modified this algorithm, especially for the filters and wavelet levels for main coronary artery segmentation in angiograms. Also, using the original SWT algorithm as used in Bankhead et al. [26], demonstrated that segmented arteries are thicker than the actual size and also increased F_{ n } rate; because it could not detect the coronary artery pixels in some parts and as a result, it decreased the performance, especially in the sensitivity value as shown in Tables 1 and 2. This part of our proposed method was implemented in MATLAB R2014a. In the implementation, we used the original functions in the Image Processing Toolbox only, and because we aimed to decrease the running time, the additional compiled MEX code was not used. In Table 3, we compared the running time of the proposed method with some stateoftheart methods. In this table, we calculated the running time for the methods proposed by Khaleel et al. [4], Frangi et al. [30], Li et al. [9] and also our proposed method using the same PC and conditions, while the running time for the method proposed by Zhou et al. [8] was obtained from their paper. As shown in Table 3, the proposed method needs lower computational time, less than 1 second, and it is a positive point for using this method in realtime systems. Therefore, we can consider it as a fast method for main coronary artery segmentation in angiograms. The results of manually labelling were evaluated by an expert cardiologist for 80 angiograms in our dataset. If the segmentation is done accurately, the proposed method for labelling will be robust to label all main coronary arteries in the angiograms.
Evaluation of the intersection tracking method for main coronary artery segmentation in CTA
The first part of the Intersection Tracking method evaluation is the seed point detection. As mentioned before, ostiums for every LCA and RCA are detected manually by seeking axial slices from top to bottom. Some algorithms have proposed automatic coronary artery seed point detection in CTA by considering the starting point from the aorta [36]. However, in some patients, there is a third coronary artery, Conus artery, which arises independently from the aorta. In this case, other algorithms failed to achieve the target seed points. Therefore, it is better to choose these points manually. As mentioned in the Intersection Tracking method, other seed points for LAD, LCX, are detected automatically. Therefore, the proposed method can be categorized as requiring minimal userinteraction by finding only the initial artery ostiums, as seed points.
Another part of the proposed method is bifurcation detection for removing subarteries. This process was evaluated using our datasets. According to the results, the LAD and LCX’s bifurcations were detected in all images in 12 datasets, but in one of them, only one of the diagonal’s bifurcation and also one of the septal’s bifurcation were missed. In addition, in two patients, the spurious bifurcations were detected. The spurious bifurcations were ignored because they did not have any intersecting region in consecutive slices. Therefore, we obtained legible results for bifurcation detection for 95.8% of our dataset.
The last part of evaluating the Intersection Tracking method is the coronary artery segmentation. As mentioned in the proposed algorithm, every coronary artery can be segmented by tracking from the corresponding seed point and the following slices, until the artery ends. This method works based on the intersection between each slice and the next. This process was evaluated for each LM, LAD, LCX and RCA arteries on 12 datasets by comparing the results with the ground truth. Based on the comparison, the LM, LAD and LCX were segmented in most patients. But, part of the RCA failed at the end slices. The reason for this problem was the corresponding slices of the end of RCA, especially in parts of the RPLB branch, were shown in previous slices but are not shown in the current slice. As shown in Figure 25, for example, parts of RPLB for “Patient8” from the dataset was missed in slice no. 236. Even though they were visible in few previous slices, slices no. 231 to 235, but they did not have intersecting region with the current slice (slice no. 236). Detecting these parts in previous slices was enough to solve this problem. We added a backtracking step to the proposed algorithm to segment all parts of the artery, especially for those shown in previous slices.
The results of the comparison between the Intersection Tracking method on the dataset and the corresponding ground truth slices are shown in Table 4. In this table, the number of slices tracked using the proposed method was shown in comparison with the same slices in the ground truth.
As shown in the table, the coronary arteries were tracked in most slices, and only failed in some, especially in those which had low contrast. The result of overlapping evaluation is shown in Figure 26.
The next step of this phase is 3D coronary artery reconstruction and labelling. As discussed before, coronary arteries are segmented as cross sectional in each slice and kept in different 3D matrices using the proposed Intersection Tracking method. It means that for all LM, LAD, LCX and RCA, 3D matrices M_{ LM }, M_{ LAD }, M_{ LCX } and M_{ RCA } are constructed respectively, which preserve each pixel of the corresponding artery in the 3D coordinate R^{3} domain. For 3D reconstruction and visualization, the OPENGL library in VC++ program was applied. A typical example of the labeled arteries for “Patient12” from the dataset is shown in Figure 27, where the various branches were labelled separately.
As all arteries are preserved in different matrices, it is possible to show all of them together. This type of visualization has many advantages for cardiologist and cardiac surgeon. All matrices M_{ LM }, M_{ RCA }, M_{ LAD }, M_{ LCX } and M_{ Aorta } can be merged together, to create one 3D matrix for the whole arterial tree. The same method was used for 3D reconstruction from the merged matrix. The whole coronary arterial tree for “Patient12” is visualized in Figure 28, in four different views. As discussed before, x and y are the coordinates of every segmented slice and z is the slice number in CTA, from top to bottom.
Our method was compared with the one and only studied on coronary artery labeling in CTA [14], where arteries were labeled in separate phases after the segmentation, centreline extraction and registration steps. In [14], first the arteries were segmented and arteries’ centreline was extracted. Then, the pointset registration method in [15] was applied for identifying main branches. However, in our proposed Intersection Tracking method, the segmentation and labelling were done together, at the same time and there was no need for other processes, such as centreline extraction or prior image registration. Therefore, it is faster and also more accurate in comparison with [14]. The overall overlap of their labelling was 91.41%. In our study, the overlap amount is related to the overlap in segmentation, which means that, if coronary artery segmentation is done accurately, the labelling will be done completely without any error as well. Since the overall segmentation amount of the proposed work is 97.38%, it is obvious that the overall overlap for labelling has the same value.
Validation of CTA and angiography registration and evaluation of the 3D position of stenosis
In this experiment, the result of the proposed method for coronary artery registration from two modalities is validated. Then, the 3D position of the stenosis of the coronary artery is evaluated. For the registration, it is not possible to evaluate the accuracy of real data, such as the heart, using a mathematical criterion. However, the specialist can visually validate them. For this, first the results of applying conventional methods, which did not use control points, on two correspondent images A and B from the dataset, are illustrated in Figure 29.
As shown in this figure, the coronary arteries were not registered well without control points. The result of extracting {C_{1},C_{2},C_{3}} and \( \left\{{C}_1^{\prime },{C}_2^{\prime },{C}_3^{\prime}\right\} \) as the two sets of the correspondent control points in A and B, respectively, and also applying the proposed method for registration, are illustrated in Figure 30.
The specialist found that the result of the proposed method for registration was much better than the conventional method. As shown in this figure, the control points in the affine transformation are registered exactly and most parts of the correspondent arteries are registered. It is worth noting that if the control points in both modalities could be registered exactly and the correspondent arteries could be registered exactly, the 3D position of the stenosis lesion would be more accurate. In addition, if the distance between the stenosis point and its two nearest control points is small, the result would be more accurate.
To evaluate the 3D position of the stenosis point on the coronary artery, 10 patients were chosen from the database as the dataset for this part. The criterion for choosing these patients were based on the fact that their stenosis points could be manually detected by the specialist in both CTA and angiography modalities. In this manner, coronary arteries and their stenosis positions on the CTA modality were known before. Therefore, the mean 3D distance between the stenosis point on the result of the proposed method and the stenosis point detected by the specialist in CTA, was calculated for each patient. In fact, the mean distances should be ideally null. It is worth noting that from the dataset, three patients had stenosis lesion in the right coronary artery (RCA) and seven patients had stenosis lesion in the left coronary artery (LCA); including 3 patients in LAD, 2 patients in LCX and 2 patients in LM arteries. Table 5 shows the results of the evaluation using the above dataset.
The mean value and standard deviation were 7.19 mm and 3.07 mm, respectively. The main reasons for these distances were related to the nature of the heart, whereby the heart is a nonsolid organ and also is beating. Therefore, there are dissimilarities between the coronary arteries of the two different modalities.
Another problem in the conventional registration algorithms is that they are time consuming; the registration part only takes about 1 min. However, in the proposed method, the registration part using the control points took about 0.31 s. This specific computational time was based on a common personal computer set up with Intel core i5 (CPU 3.2 GHz) and 8 GB RAM, and the program being created in MATLAB R2014a. Furthermore, this computational time can be decreased further by utilizing better stateoftheart computer. It is worth noting that the proposed method provides further information, which can complement the angiogram and CTA. By using this method, it is possible to visualize both the correspondent coronary arteries from angiogram and CTA together. This kind of information assists the cardiac surgeon and cardiologist to make a decision regarding whether an artery needs to be dilated or not.
Conclusions
A new method is proposed for coronary artery registration in both CTA and angiography modalities to provide the 3D position of the stenosis lesion diagnosed in an angiogram. The result has benefits for CAD diagnosis and treatment. Tests using the dataset demonstrated that the proposed method aided the specialists to find the location of the stenosis lesion and also determine the visual relationship between the correspondent coronary arteries. To the best of our knowledge, there is no hybrid device for both CTA and angiography modalities, yet. In addition, based on the literature, no algorithm has been proposed for registering these two modalities in order to obtain the 3D position of the stenosis point. Therefore, the proposed method can be seen as a new contribution. The aim of the proposed algorithm is applicable and portable for common personal computers, also with respect to the standard medical acquisition methods. In the future, we intend to perform a comprehensive research on stenosis detection in angiography images and apply the proposed method to automatically mark the stenosis point in the angiogram.
Abbreviations
 CTA:

Computed tomography angiography
 CAD:

Coronary artery disease
 RCA:

Right coronary artery
 LCA:

Left coronary artery
 LCX:

Left circumflex coronary artery
 LAD:

Left anterior descending coronary artery
 DWT:

Discrete wavelet transform
 SWT:

Starlet wavelet transform
 HU:

Hounsfield unit
 RMS:

Root mean squared
References
 1.
Mazaheri s., Sulaiman P., Wirza R. and Moosavi Tayebi R., “Echocardiography Image Segmentation: A Survey”, 2nd International Conference on Advanced Computer Science Applications and Technologies – ACSAT2013, Sarawak, Malaysia, 2013.
 2.
Moosavi Tayebi R, Sulaiman P, Wirza R, Dimon Z, Kadiman S, Nurliyana A, et al. Coronary artery segmentation in angiograms with pattern recognition techniques – A survey. In: IEEE Proceeding of International Conference on Advanced Computer Science Applications and Technologies. 2013. p. 321–6.
 3.
Moosavi Tayebi R, Sulaiman P, Wirza R, Dimon Z, Kadiman S, Mazaheri S. A fast and accurate method for automatic coronary arterial tree extraction in angiograms. J Comput Sci. 2014;10(10):2060–76.
 4.
Khaleel H, Wirza R, Zamrin Dimon M, Ramlan M, Norwati M. Extraction of coronary artery trees in angiocardiography images. Sci Res Essays. 2012;7(47):4014–36.
 5.
Xu Y, Zhang H, Li H, Hu G. An improved algorithm for vessel centerline tracking in coronary angiograms. Comput Methods Programs Biomed. 2007;88(2):131–43.
 6.
HernándezVela A, Gatta C, Escalera S, Igual L, MartinYuste V, Sabate M, et al. Accurate coronary centerline extraction, Caliber estimation and catheter detection in angiographies. IEEE Trans Inf Technol Biomed. 2012;16(6):1332–40.
 7.
Wang Y, Shu H, Zhou Z, Toumoulin C, Coatrieux JL. Vessel extraction in coronary Xray angiography. In: 27th IEEE Annual International Conference of the Engineering in Medicine and Biology Society. 2005. p. 1584–7.
 8.
Zhou S, Yang J, Chen W, Wang Y. New approach to the automatic segmentation of coronary artery in Xray angiograms. Sci Chi S F Inf Sci. 2008;51(1):25–39.
 9.
Li C, Kao CY, Gore JC, Ding Z. Minimization of regionscalable fitting energy for image segmentation. IEEE Trans Image Process. 2008;17(10):1940–9.
 10.
Kirişli H, Schaap M, Metz C, Dharampal A, Meijboom WB, Papadopoulou S, et al. Standardized evaluation framework for evaluating coronary artery stenosis detection, stenosis quantification and lumen segmentation algorithms in computed tomography angiography. Med Image Anal. 2013;17(8):859–76.
 11.
Goldenberg R, Eilot D, Begelman G, Walach E, BenIshai E, Peled N. Computeraided simple triage (CAST) for coronary CT angiography (CCTA). Int J Comput Assist Radiol Surg. 2012;7(6):819–27.
 12.
Öksüz İ, Ünay D, Kadıpaşaoğlu K. A hybrid method for coronary artery stenoses detection and quantification in CTA images. In: MICCAI Workshop 3d Cardiovascular Imaging: A MICCAI Segmentation. 2012.
 13.
Yang G, Kitslaar P, Frenay M, Broersen A, Boogers MJ, Bax JJ, et al. Automatic centerline extraction of coronary arteries in coronary computed tomographic angiography. Int J Cardiovasc Imaging. 2012;28(4):921–33.
 14.
Yang G, Broersen A, Petr R, Kitslaar P, de Graaf MA, Bax JJ, et al. Automatic coronary artery tree labeling in coronary computed tomographic angiography datasets. Comput Cardiol. 2011;38:109–12.
 15.
Myronenko A, Song X. Point set registration: Coherent point drift. IEEE Trans Pattern Anal Mach Intell. 2010;32(12):2262–75.
 16.
Turgeon GA, Lehmann G, Guiraudon G, Drangova M, Holdsworth D, Peters T. 2D3D registration of coronary angiograms for cardiac procedure planning and guidance. Med Phys. 2005;32(12):3737–49.
 17.
Ruijters D, ter Haar Romeny BM, Suetens P. Vesselnessbased 2D–3D registration of the coronary arteries. Int J Comput Assist Radiol Surg. 2009;4(4):391–7.
 18.
Metz CT, Schaap M, Klein S, Baka N, Neefjes LA, Schultz J, et al. Registration of coronary CTA and monoplane Xray angiography. IEEE Trans Med Imaging. 2013;32(5):919–31.
 19.
Imamura H, Ida N, Sugimoto N, Eiho S, Urayama SI, Ueno K, et al. Registration of preoperative CTA and intraoperative fluoroscopic images for assisting aortic stent grafting. In: Springer Medical Image Computing and ComputerAssisted Intervention  MICCAI 2002. 2002. p. 477–84.
 20.
Chen X, Gilkeson RC, Fei B. Automatic 3Dto2D registration for CT and dualenergy digital radiography for calcification detection. Med Phys. 2007;34(12):4934–43.
 21.
Gatta C, Balocco S, MartinYuste V, Leta R, Radeva P. Nonrigid multimodal registration of coronary arteries using SIFTflow. In: Springer Pattern Recognition and Image Analysis. 2011. p. 159–66.
 22.
Serradell E, Romero A, Leta R, Gatta C, MorenoNoguer F. Simultaneous correspondence and nonrigid 3D reconstruction of the coronary tree from single Xray images. In: IEEE International Conference on Computer Vision (ICCV). 2011. p. 850–7.
 23.
Mazaheri S, Sulaiman P, Wirza R, Dimon Z, Khalid F, Moosavi Tayebi R. “Hybrid Pixelbased Method for Cardiac Ultrasound Fusion Based on Integration of PCA and DWT”, in Computational and Mathematical Methods in Medicine Journal, (MMMI14). 2014.
 24.
Genovesio A, Liedl T, Emiliani V, Parak WJ, CoppeyMoisan M, OlivoMarin JC. Multiple particle tracking in 3D + t microscopy: method and application to the tracking of endocytosed quantum dots. IEEE Trans Image Process. 2006;15(5):1062–70.
 25.
Starck JL, Murtagh F. Astronomical image and data analysis’, Springer Science & Business Media. 2007.
 26.
Bankhead P, Scholfield CN, McGeown JG, Curtis TM. Fast retinal vessel detection and measurement using wavelets and edge location refinement. PLoS One. 2012;7(3):1–12.
 27.
Starck JL, Fadili J, Murtagh F. The undecimated wavelet decomposition and its reconstruction. IEEE Trans Image Process. 2007;16(2):297–309.
 28.
Starck, JL, Murtagh F, Fadili JM. ‘Sparse image and signal processing: wavelets, curvelets, morphological diversity’, Cambridge University Press, 2010. https://scholar.google.com.my/scholar?hl=en&q=Sparse+image+and+signal+processing%3A+wavelets%2C+curvelets%2C+morphological+diversity&btnG=&as_sdt=1%2C5&as_sdtp=.
 29.
Gu S, Gupta R, Kyprianou I. Computational highresolution heart phantoms for medical imaging and dosimetry simulations. Phys Med Biol. 2011;56(18):5845–64.
 30.
Frangi AF, Niessen WJ, Vincken KL, Viergever MA. Multiscale vessel enhancement filtering, Medical Image Computing and ComputerAssisted Interventation  MICCAI’98, Springer, 1998: 130–37. http://link.springer.com/chapter/10.1007/BFb0056195#.
 31.
Makela T, Clarysse P, Sipila O, Pauna N, Pham QC, Katila T, et al. A review of cardiac image registration methods. IEEE Trans Med Imaging. 2002;21(9):1011–21.
 32.
Baka N, Metz C, Schultz C, van Geuns RJ, Niessen W, van Walsum T. Oriented Gaussian mixture models for nonRigid 2D/3D coronary artery registration. IEEE Trans Med Imaging. 2014;33(5):1023–34.
 33.
Baka N, Metz C, Schultz C, Neefjes L, van Geuns RJ, Lelieveldt BP, et al. Statistical coronary motion models for 2D+ t/3D registration of Xray coronary angiography and CTA. Med Image Anal. 2013;17(6):698–709.
 34.
Han J, Kamber M, Pei J. ‘Data mining: concepts and techniques’, Morgan Kaufmann Publisher, 2006. https://books.google.com.my/books?hl=en&lr=&id=AfL0tYzOrEC&oi=fnd&pg=PP2&dq=Data+mining:+concepts+and+techniques%E2%80%99,+Morgan+Kaufmann+Publisher&ots=UwUYvV6kx0&sig=wp_SNybYxFwmR4Tohu3u5cw7MQw#v=onepage&q=Data%20mining%3A%20concepts%20and%20techniques%E2%80%99%2C%20Morgan%20Kaufmann%20Publisher&f=false.
 35.
Monteiro FC, Campilho AC. Performance evaluation of image segmentation. In: Springer Image Analysis and Recognition. 2006. p. 248–59.
 36.
Bouraoui B, Ronse C, Baruthio J, Passat N, Germain P. 3D segmentation of coronary arteries based on advanced mathematical morphology techniques. Comput Med Imaging Graph. 2010;34(5):377–87.
Acknowledgment
This research was supported by Ministry of Science, Technology and Innovation, Malaysia, under Science Fund number UPM0007353.
Author information
Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
All authors have made substantial contributions to conception and design, acquisition of data, analysis and interpretation of the data. All authors read and approved the final manuscript.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Moosavi Tayebi, R., Wirza, R., Sulaiman, P.S.B. et al. 3D multimodal cardiac data reconstruction using angiography and computerized tomographic angiography registration. J Cardiothorac Surg 10, 58 (2015). https://doi.org/10.1186/s1301901502492
Received:
Accepted:
Published:
Keywords
 Angiography
 Computerized tomography angiography
 Segmentation
 Labeling
 Multimodal registration
 3D reconstruction