The proposed methodοlogy includes 6 stages: in the first stage the CTA images are pre-processed to remove the artifacts and detect the vessel silhouette. In the second stage the artery borders are roughly defined and used to extract the initial centerline of the vessel. In the third stage, the plaque parameters are adapted in the CTA's HU scale to optimize a 4-component Gaussian Mixture Model (GMM) [22] to the HU histogram. In the fourth stage, radial images are produced perpendicular by to the initial centerline which is further processed in the fifth stage where the lumen and the outer wall border of the vessel are detected and the plaque region is estimated. Finally, in the sixth stage the 3D surfaces of the lumen - outer vessel wall and CP plaques are constructed. The stages of the proposed methodology are shown schematically in Fig. 1.
Preprocessing
A preprocessing stage is applied on the acquired CTA images which includes:
-
contrast enhancement [23], which maps the intensity values of the image to new values such that, 1 % of the data is saturated at low and high intensities of the input data,
-
image thresholding [24], which allows identification of the vessel [25],
-
and the application of a Frangi Vesselness filter (using the default parameters suggested by Frangi et al.) [26] which permits detection of structures which correspond to coronary vessels.
Rough artery border detection
The following process is implemented to extract the initial lumen borders and define a first approximation of the luminal centerline:
-
The user defines a seed point in the main artery (the only user action in the whole methodology).
-
A region growing method is applied using 26-connected pixels neighbor connectivity for lumen border detection.
-
Edge Detection is applied to estimate the lumen borders.
-
The lumen border are used to define an initial vessel centerline using the center of gravity of each curve [27]. Then the centerline is smoothed (using moving average filter) in order to have a smooth first derivative of the centerline curve.
The initial centerline points are fitted using B-Splines. B-Spline is used to provide smooth first derivatives needed for producing a radial image perpendicular to a specific centerline point (Fig. 2).
Image HU modeling
The distribution of HU intensity over the artery is modeled by a Gaussian Mixture Model (GMM) [22]. A GMM model is a parametric probability density function which is represented as a weighted sum of M Gaussian component densities given as:
$$ p\left(x\Big|\lambda \right)={\displaystyle \sum_{i=1,\;..M}{w}_iN\left(x\Big|{\mu}_i,{\sigma}_i\right)}, $$
(1)
where x is a 1-dimensional continuous data (HU value), w
i
are the mixture weights, μ
i
and σ
i
are the mean vector and the standard deviation of the component Gaussian densities (N(x|μ
i
,σ
i
)), respectively and λ is the set:
$$ \lambda =\left\{{w}_i,\kern0.24em {\mu}_i,\kern0.24em {\sigma}_i\right\} $$
(2)
The maximum posterior probability of the sample x to belong to a specific component k, is given as:
$$ P\left(M=k\Big|x\right)=\frac{w_kN\left(x;{\mu}_k,{\sigma}_k\right)}{{\displaystyle \sum_{i=1..M}{w}_iN\left(x;{\mu}_i,{\sigma}_i\right)}}. $$
(3)
In our approach each component of the GMM models corresponds to one of the following classes: i) lumen, ii) calcified plaque, iii) non-calcified plaque and iv) background. In order to train the GMM model the user manually selects 4 points that correspond to the four classes ({X
L
,X
CP
,X
NCP
,XB}) and then each pixel is classified to one of the four classes based on the class/component with the maximum posterior probability. For each selected point, 5 random points are automatically detected in a circular neighborhood of radius 1.
Radial Image Projection (RIP)
In the initially extracted centerline which is sampled every 0.2 mm, a point is defined and a radial image is produced perpendicularly to the centerline point using the B-spline derivatives extracted on the specific point [27]. These images are called Radial Intensity Projection (RIP) images. For creating the RIP images we used the 3D derivative vector of the b-spline corresponding to the centerline. On this way the continuity of the RIP images is achieved. For an accurate RIP image the estimated centerline point is required to lie in the middle of the lumen. However, as described in Section II-B the centerline is smoothed in order to avoid surface intersections as it is depicted in Fig. 3. This smoothing results in an offset of the estimated centerline from lumen’s center. In order to overcome this problem an iterative correction process is applied. The iterative process includes the following steps:
-
The two bullets should be before Fig. 3. The radial image is produced for the given centerline and the current offset (initially the offset is set to zero).
-
The center of mass of the radial image is calculated and a new offset is computed.
The process is repeated five times (2, 3, 4 and 6 repetitions were also tested), number which is sufficient to minimize the offset. An example of the above procedure is depicted in Fig. 4 for a synthetic image with a binary circle with an initial large offset from the center (Fig. 4-(a)). The initial center of the circle should lie in the center of the image. It can be observed that after a few iterations the radial image is corrected (Fig. 4-(d)). After the correction step the final RIP image is used for the object extraction process. An example of the RIP image is given in Fig. 5.
Inner - outer vessel wall and plaque region detection
From the RIP of each slice we extract the Region of Interest (ROI) of the artery. The ROI is extracted based on the detection of the zero crosses (HU < −130) [25], (Fig. 5-b). The areas below the zero crossing are removed (Fig. 5-c). After defining the ROI, the HU profile of the RIP image is calculated for each radius interval (Fig. 5-d). Based on the profile, we classify the ROI using the GMM model (Fig. 5-e).
Based on the ROI classification we define profile regions. The lumen wall is set to the first non-lumen point. The outer vessel wall is set to the last non-background point. Therefore, two curves are created one for the lumen and one for the outer vessel wall. After smoothing, the curves are transformed back to the Cartesian coordinates.
Using the original radial image, the area between the inner and outer vessel wall defines the plaque region. The classified ROI (Fig. 6-(a)) is thresholded (the threshold is 0.5) according to the GMM probability estimations, keeping only pixels belonging to the specific plaque type (Fig. 6-(b)) resulting to binary objects. Then these binary objects are identified and the boundary of those objects is extracted. Smoothing is performed on the extracted boundary curves. Based on the physical size of each pixel the area of the plaque can be determined.
3D Surface construction
N points for each of the inner/outer vessel wall and plaque lesions are extracted clock wisely for each image of the centerline. Then a triangulation approach is implemented to construct the mesh surfaces. The border points (N) (inner/outer wall or plaque) of two sequential images, frame i and frame i + 1 are connected constructing a triangle mesh. Thus 2xN surfaces are constructed for each pair of sequential images.
Dataset
The data acquired from 20 patients, that underwent coronary angiography, IVUS imaging and CTA for clinical purposes, were used to validate the proposed methodology. Data collection were provided by the Invasive Cardiology Laboratory of the IFC-CNR, Pisa and by Heart Institute, University of Sao Paulo, São Paulo, Brazil. The study protocols were approved by the local Ethics Committees (IFC-CNR data: Comitato Etico, Azienda Ospedaliero-Universitaria Pisa, Pisa, Italy and Sao Paulo Data: Ethics Committee of the Clinics Hospital of the University of Sao Paulo Medical School - CAPPesq, Hospital das Clínicas da Faculdade de Medicina da U S P, Sao Paolo, Brazil) and written informed consent for participation in the study was obtained from participants.
IFC-CNR data
Coronary CTA was performed using a 64-slice scanner (Light Speed, General Electric Health Care). Eight patients with a significant (>50 %) luminal stenosis on CTA had further invasive investigations. Coronary angiography was performed with the use of a flat panel system (Innova 2100, GEHC). To further assess the severity of the detected lesions on X-ray angiography, IVUS imaging was performed with the use of 64-crystal electronic ultrasound probe (Eagle-Eye, Volcano corporation) that was pulled-back by an automated pull back device at a speed 0.5 mm/sec. The acquired IVUS and CTA data were stored in DICOM format and transferred to a workstation for off-line analysis.
Sao Paulo data
Twelve patients were enrolled in the study and had a clinical indication for CTA according to the guidelines and significant (>50 %) luminal stenosis was confirmed by CTA. The patients were further evaluated for coronary lesions using IVUS. The IVUS examination was performed using a 20 MHz electronic multi-array 2.9 F catheter (Eagle Eye®, Volcano Corporation Inc) connected to a dedicated console (InVision Gold®, Volcano Corporation Inc., San Diego, CA, USA). Acquisition was performed during motorized pullback at a constant speed of 0.5 mm/s (R-100® pullback device, Volcano Corporation Inc., San Diego, Ca, USA). Patients with a heart rate greater than 65 bpm received up to 15 mg of intravenous metoprolol before the acquisition of MDCT images, unless contra-indicated. Sublingual nitrate was given to all patients prior to the acquisition of MDCT images. All scans were performed utilizing a 64-slice MDCT scanner (Aquillion 64TM, Toshiba Medical Systems, Japan). After the intravenous injection of contrast media (Iopamidol, 370 mg iodine/ml, Bracco), the helical scan for CT coronary angiography was triggered once a threshold of 180 Hounsfield units (HU) was reached at the descending aorta.
All CTA (eight IFC-CNR and twelve Sao Paulo) data were processed using the proposed methodology and their corresponding IVUS data (eight IFC-CNR and twelve Sao Paulo) using a previously described plaque characterization method [21]. CP was detected in 20 IVUS examinations and were used to compare the 3D CP detected by the proposed CTA methodology.
CT-IVUS comparison
To compare the 3D reconstructed arteries and the detected 3D CP detected by the proposed methodology, IVUS plaque characterization [21] and 3D IVUS reconstruction [28] was performed in corresponding segments and the non-invasive-based (CTA) and invasive imaging (IVUS) based models were compared. Experts selected manually, segments of arteries which were assessed with both CTA and IVUS imaging. Side branches were used as fiducial points to define correspondence between CTA and IVUS.
Artery and CP reconstruction using IVUS
IVUS based reconstruction was performed using the methodology introduced in [28]. The method combines IVUS and X-ray angiography and places the lumen and media adventitia borders detected in IVUS frames [21] onto the 3D luminal centerline derived from two angiographic projections. In the present work instead of using the 3D centerline derived from angiography, we used the 3D center line derived from CTA. To detect the CP in the selected IVUS segments a recently developed plaque characterization method [21] was employed. Hence, the lumen-outer vessel wall morphology and the CP of the arteries are reconstructed by combining IVUS and CTA.
Comparison metrics
The inner and outer walls reconstructed by the two modalities were compared using their volume and surface area. The 3D CP objects detected by the two modalities were compared using the following metrics: volume, surface area (the total area of the object's faces), maximum length, inner angle and the overlapping volume (volume detected in both CTA and IVUS models) between two objects. As maximum length of an object we denote the distance between the first and the last voxel of the object along the centerline and as inner angle θ
In
, we denote the mean of the angles that define the circumferential extent of the object (Fig. 7). For an object having n 2D slices and angle θ angles, the θ
In
is defined as:
$$ {\theta}_{In}=\frac{{\displaystyle \sum_{i=1}^n{\theta}_i}}{n} $$
(4)
