Cardiovascular disease (CVD) is the leading cause of mortality and morbidity in the United States.1 An important factor in the pathophysiology of CVD is the composition and remodeling of the myocardium. Myocardial tissue includes muscle, adipose tissue, collagen fibers, and fibrotic myocardium, and the relative percentage of each varies by chamber and with the progression of the disease. Myofiber disarray can impair electrical conduction and result in arrhythmias or hypertrophic cardiomyopathy.2 The presence of adipose within the myocardium provides a high indication of arrhythmogenic cardiomyopathy3 and thickened layers of collagen fibers imply severe myocardial scar.4 The diffusion of myocardial fibrosis, a fundamental process in the remodeling of cardiomyopathy, is postulated to cause increased cardiac stiffness and poor clinical outcomes.5 Due to the importance of myocardial tissue composition on normal heart function and CVD, characterization of myocardial tissue can facilitate the evaluation of tissue remodeling, identification of arrhythmogenic substrates, and diagnosis of CVD.
In the past decades, medical imaging modalities including ultrasound,67.–8 multidetector computed tomography,9,10 and magnetic resonance imaging have been used to characterize cardiac tissue compositions such as collagen region during myocardial infarction,6,9,10 adipose tissue,8,11 or organization of myofibers within myocardium.7,12 However, the abovementioned modalities suffer from either a low image resolution or a long data acquisition time. Optical coherence tomography (OCT) has been demonstrated to have the ability to image biological tissue at a fast rate with a high resolution () with a 2 mm imaging range13,14 in the axial direction. Previous research efforts demonstrated that OCT can image important features within the heart15 such as the purkinjie network,16 atrial ventricular nodes,17,18 sinoatrial nodes,19 and myofiber organization.2021.22.–23 Given that the wall thickness in the human atria ranges from 2 to 5 mm,24 OCT has the ability to visualize a large percentage of the human atrial wall. There is a great potential to classify tissue compositions within human atria via OCT imaging. However, manual interpretation of OCT images is time consuming and not applicable for analysis on large three-dimensional (3-D) volumetric datasets. Therefore, automated identification of tissue composition in human atria from OCT images is greatly needed.
In this study, we present an image-processing algorithm to automate the classification of tissue compositions within atrial OCT images. Layer structures of multiple tissue compositions are automatically segmented using graph searching without any prior information. Within each layer, optical properties, statistical measurements, and texture features are extracted. Features are subsequently used to build a statistical classification model to distinguish tissue compositions of dense collagen, loose collagen, adipose tissue, normal myocardium, and fibrotic myocardium. Our work enables, to the best of our knowledge, the first automated classification of myocardial tissue compositions from human atrial OCT images.
Human hearts () were obtained under two approved protocols from the National Disease Research Interchange (NDRI). The inclusion criteria for the first NDRI protocol are based on the following diagnosis: end stage heart failure, cardiomyopathy, coronary heart disease, or myocardial infarction. The second protocol requests normal hearts. Fresh tissue samples were shipped submerged in ice-cold phosphate-buffered saline and received within 48 h of donor death. Detailed characteristics of the donor hearts within this study are listed in Table 1.
Clinical characteristics of heart donors in dataset.
|Age in yrs, median (interquartile range)||66.0 (62.25 to 69.75)|
|Male, n (%)||5 (33.3)|
|BMI, median (interquartile range)||29.25 (24.3 to 34.2)|
|Medical history, n (%)|
|Heart failure||2 (13.3)|
|Cause of death, n(%)|
|Cardiac arrest||3 (20.0)|
|Cardiopulmonary arrest||3 (20.0)|
|Respiratory arrest||1 (6.7)|
|Respiratory failure||2 (13.3)|
|Chronic obstructive pulmonary disease||4 (26.7)|
|Congestive heart failure||2 (13.3)|
|Complete characteristics were not available for all donors|
All samples were imaged ex vivo, using a spectral domain OCT system, Telesto I (Thorlabs GmbH, Germany). The system is an InGaAs-based system with its source centered at 1325 nm and a bandwidth of 150 nm. The axial and lateral resolutions are 6.5 and in air, respectively. All datasets were acquired at 28 kHz. In our experiments, each volume consists of , corresponding to a tissue volume of (in air). To extract the raw OCT data, the postprocessing algorithm, including to space interpolation, windowing, and Fourier transform, was implemented using MATLAB 2014b (Mathworks, Inc., Massachusetts).
Sections of tissue from the imaging field of view were processed for histopathology. Samples were sectioned parallel to the direction of the B-scans. Sample pieces were cut corresponding to the size of the OCT volume, fixed in formalin for , and then placed in ethanol (20%) for . After fixation, samples were stained with Masson Trichrome. For 33.3% of the samples, histology was taken every 2 mm to ensure multiple matches between histology and the OCT image set. For validation, two investigators, blind to the automated results, segmented and classified the images based on the histology. One-way analysis of variance with Tukey multiple comparison test tissue were performed to detect differences between the tissue compositions for each of the extracted features. A -value of 0.05 was considered statistically significant. All statistical analysis was conducted with the software package Prism 6.03 2013 (GraphPad Software, Inc, California).
To identify tissue compositions within OCT images, we present a region-based classification method. A schematic of the workflow for the analysis of two-dimensional (2-D) images is shown in Fig. 1. The algorithm consists of three steps: layer segmentation, feature extraction, and tissue classification. In each B-scan, OCT images were first segmented into multiple regions. Within the segmented region, features such as optical properties, texture analysis, and high order statistical moments were extracted. The features were inputs to a tissue classifier whose output was the tissue type for the region.
For the first step, we divided OCT images into multiple layers through segmentation for future feature extraction and classification. Compared with existing segmentation methods, segmenting OCT images of atrial tissue is challenging. Within prior work, layer boundaries were automatically determined by minimizing a cost function25 or building a minimum weight graph.26 The weighting scheme, and searching order are determined by prior knowledge of the layer structure,27,28 such as empirical thickness measurements and knowledge of bright-to-dark transition patterns between two layers. Unfortunately, neither empirical thickness measurement or transition patterns is consistent for atrial tissue. In the atrium, layer thicknesses and tissue composition vary within a normal heart depending on the region that is imaged. Furthermore, the layer thicknesses and tissue composition change with the progression of the disease. Therefore, the first step of our algorithm is layer segmentation, which includes preprocessing, layer information estimation, and boundary searching. A detailed flowchart of the layer segmentation steps is shown in Fig. 2(a). The preprocessing step improves the image quality through denoising and flattens the image to reduce the boundary searching range. The layer information estimation step determines the number of tissue compositions and identifies starting points for boundary searching. The image is segmented after boundary searching.
The preprocessing procedure includes image denoising and flattening. Given that OCT images are generally corrupted by speckle noise,29 we used a block matching 3-D (BM3-D)30,31 method to denoise OCT images and enhance the boundaries. Briefly, for the BM3-D algorithm we divide the original OCT image into multiple blocks and denoise similar blocks. The BM3D method exploits the sparsity of structural information and is thus considered to be a good tool to denoise speckle noise and enhance boundary information. To reduce the searching range and maintain a smooth searching shape, we flattened and shifted the filtered image based on the tissue surface. To flatten/shift images in a fast manner, we undersampled the original image. In cardiac tissue, the most hyper reflective surface is in the endocardium. Within a downsampled A-line, we thus estimated the location of the maximum pixel value as the axial location of the surface. Then the image was shifted based on an interpolation of the axial location of the maximum value-pixels within the downsampled image.
Layer information estimation
After image denoising and flattening, we estimated the layer information within each OCT image. The layer information consisted of the number of layers and the initial points for boundary searching between layers. To count the number of layers, we analyzed averaged A-lines in the OCT image. Due to the positioning of the sample and uneven sample surface, A-Lines around the center of the B-scan had better image quality than the regions toward the left and right edges of the image. To ensure an accurate estimation, five A-Lines were selected around the center of the image. The A-Lines were away from each other. For each A-Line, 20 A-lines were averaged. In each averaged A-line, the intensity curve was linearly fitted using a sliding window. The algorithm flow of linear fitting is shown in Fig. 2(b). We first set the location of the maximum pixel value as the first anchor. Within a window, the intensity was linearly fitted. We calculated the root mean square deviation (RMSD) error as following:
The boundaries were searched from the center of the image and progressed outward to the left and right. The boundary search algorithm minimizes the following function:Figure 2(c) shows a schematic of the boundary searching algorithm, starting from . The cost of pixels was examined and the pixel with the highest weight was considered to be the boundary of the layer within . We then estimated the boundary point for the next column. The searching algorithm was run in parallel for multilayers within an image.
Within each segmented region, we extracted features from the OCT images to study different patterns of tissue compositions. The extracted features can be divided into three categories: measured optical properties, statistical moments, and texture analysis.
Measured optical properties
Optical property parameters that we studied were attenuation coefficients () and penetration depth (mm). Attenuation coefficient was measured based on the method mentioned in Ref. 32. Penetration depth was defined as the depth at which the intensity drops to of its original intensity33 when light first reaches the layer. Additionally, we calculated the distance between the centers of layers to the tissue surface.
We performed histogram equalization and median filtering on raw OCT data. Then we calculated the statistics of high order moments (skewness and kurtosis), on the intensities of the denoised image within the whole layer to analyze the distribution of intensity for various tissue types.
We encoded OCT images with texture on equalized and filtered OCT images. Texture feature number (TCN)34 is assigned to each pixel. In TCN, the local feature of each pixel is represented by the intensity change of its eight neighboring pixels. We analyzed the statistics of the TCN number within each layer. In particular, we calculated the coarseness and homogeneity from the histogram of the TCN. We also quantified the mean and standard deviation from the default texture analysis tool in MATLAB, such as range filter and std filter, within each region. We also quantified entropy within the region. In addition, we constructed gray level cooccurrence matrices (GLCM)35 to extract more additional texture features. Specifically, contrast, energy, and correlation were computed by setting the number of levels to 16.
Representative parametric images obtained from left and right atrial samples are shown in Fig. 3, where we presented typical pixel-based, A-Line-based, and layer-based features. From pixel-based parametric images, such as attenuation coefficients and entropy in Figs. 3(c)–3(d) and Figs. 3(i)–3(j), large variations can be observed within a single layer. In A-line-based features, such as distance to the surface, there are smaller variations within each layer and the differences between tissue compositions can be well observed. For layer-based features, we performed texture analysis on pixels within the whole layer. Representatives are shown in Figs. 3(f) and 3(i). To simplify our model, we averaged the pixel-based and A-Line-based feature values within each layer. This resulted in a vector of features for each layer. The number of entries in the vector was the number of features. In this study, we extracted 16 features, listed in Table 2.
List of features used in the classifier.
|1||Attenuation coefficients (mean)||Mean value of attenuation coefficient|
|2||Attenuation coefficients (std)||Standard deviation of attenuation coefficient|
|3||Penetration depth (mean)||Mean value of penetration depth|
|4||Penetration depth (std)||Standard deviation of penetration depth|
|5||Std filter value (mean)||Mean value of standard deviation filtering results|
|6||Std filter value (std)||Standard deviation value of standard deviation filtering results|
|7||Range filter (mean)||Mean value of range filtering results|
|8||Range filter (std)||Standard deviation value of range filtering results|
|9||Entropy||Entropy of the pixel values|
|10||Coarseness (TCN)||Coarseness analysis of texture code number|
|11||Homogeneity (TCN)||Homogeneity feature of texture code number|
|12||Contrast (GLCM)||Contrast feature from GLCM|
|13||Energy (GLCM)||Energy feature from GLCM|
|14||Distance to surface||The distance between the center of the layer and the surface|
|15||Skewness||Skewness within the whole layer|
|16||Kurtosis||Kurtosis within the whole layer|
Relevance Vector Machine Classifier
The relevance vector machine (RVM)36 was used to classify tissue compositions. For each feature vector , the probability of the vector belonging to a specific tissue composition was determined by the following equation:
New values for the weight vector were estimated by calculating the derivative of the expectation of the weights and the Hessian matrix of the weights H. Using a Newton update method, the new weights were estimated as
RVM is a Bayesian framework of the support vector machine (SVM), which is widely used in classification37,38 and segmentation.39 Compared with the SVM model, RVM obtains sparser solutions for weight vector . This is done by adopting a nonGaussian prior for multiple hyperparameters , which only requires a limited number of weights to be “active”. Once the values for the hyperparameter are optimized, most of the hyperparameters tend to move toward infinity. This results in most weights getting closer to zero, and becoming “irrelevant” in establishing a decision boundary. Only relevant weights are retained, which produced a significantly lower number of relevance vectors compared to SVM.
A leave-one-out experiment was conducted on the whole dataset. We used OCT images from 14 hearts as training data and used the images from the remaining one heart as testing data. The experiment was repeated such that images from any heart would be the testing dataset once.
Given boundary information from the B-scans, we reconstructed the volumetric classification for myocardial tissue. Upon estimating the boundary from each B-scan, each layer can be roughly estimated. The detected boundary was arranged along the direction perpendicular to the B-scan. We further smoothed the estimated surface using a median filter and reconstructed the 3-D surface based on the smoothed plane. The estimated layer boundaries for each B-scan were modified accordingly. We then performed the tissue classification algorithm on each fine-tuned region in each B-scan. After performing the classification algorithm on each B-scan, the 3-D classification results were realigned. We overlaid the tissue composition with OCT volumetric data in an hue saturation value (HSV) scheme. In this study, tissue composition was encoded as hue; saturation and value are encoded as intensity. All 3-D results were visualized using the software package Amira 5.4.3 2012 (Zuse Institute Berlin, Germany).
Figure 4 shows two typical segmentation results of human atrial OCT volumes. In Figs. 4(a)–4(c), three layers are dense collagen (dark blue in histology), loose collagen (light blue in histology), and normal myocardium (red in histology) while in Fig. 4(d)–4(f), three layers are dense collagen (dark blue in histology), adipose tissue (white in histology), and normal myocardium (red in histology). In general, segmented boundaries matched the visible boundaries in the trichrome histology images. Moreover, we conducted quantitative comparisons between automated segmentations and manual segmentations from two observers for images from all 15 human hearts with corresponding histology. The results are listed in Table 3. The difference between automated segmentation and manual segmentation was , which is comparable to results provided by the two investigators, . To visualize boundaries in 3-D, multiple consecutive B-scans were segmented based on the filtered surface using the method in Sec. 2.8. The 3-D segmentation results are shown in Figs. 4(c) and 4(f).
Comparison between automated segmentation and manual measurements from two observers.
|Automated versus observer 1||Automated versus observer 2||Observer 1 versus Observer 2|
|Mean (μm)||Std (μm)||Mean (μm)||Std (μm)||Mean (μm)||Std (μm)|
We performed statistical comparison of features for five tissue compositions: normal myocardium, loose collagen, adipose tissue, fibrotic myocardium, and dense collagen. Measurements from optical properties, statistics, and texture analysis were shown in Fig. 5. In general, we found that features from OCT images of endomyocardium had a strong correlation with tissue composition. We found that normal myocardial tissue was only significantly different from the all four tissue compositions in homogeneity (). Loose collagen was not significantly different from fibrotic in homogeneity and energy, but showed statistical differences in statistical moments (kurtosis, ) and optical properties (attenuation coefficients ). Within texture analysis, dense collagen was not significantly different from adipose tissue in homogeneity, but had significantly different values in energy (). Similar observations were found in the rest of features, thus all 16 features were used within the classification model.
We performed classification experiments on 60 B-scans from 15 human hearts. Representative classification results are shown in Fig. 6. Similar to segmentation results, we performed comparison between automated classifications and trichrome histology and presented the comparison in Fig. 6(a)–6(f), of which 6(a) and 6(d) are raw OCT images, 6(b) and 6(e) are color-coded classifications, and 6(c) and 6(f) are histology images. The tissue compositions were color coded in HSV where hue encoded the tissue composition and saturation and value encoded intensity. Two layers are dense collagen (dark blue in histology) and adipose tissue (white in histology) in Figs. 6(a)–6(c) while dense collagen (dark blue in histology) and fibrotic myocardium (purple in histology) are shown in Figs. 6(d)–6(f). Both classification results agree with the trichrome histopathology.
A leave-one-out experiment was conducted on the whole dataset. We used OCT images from 14 hearts as training data and used the images from the remaining one heart as testing data. The experiment was repeated such that images from any heart will be the testing dataset once. To validate the accuracy, we compare the automated classification result with the tissue types shown in histology on a layer-wise basis. We obtained the confusion matrix to assess overall accuracy, Fig. 7. The final tissue classification estimation for the region was the class with the highest probability. Using this identification rule, we achieved an average accuracy of 80.41% for classifying the five tissue compositions.
Representative 3-D classification results are shown in Fig. 8. The classified tissue compositions are overlaid on the original 3-D dataset using the HSV scheme. Gold, yellow, red, and blue hues represent dense collagen, loose collagen, normal myocardium, and adipose tissue, respectively. In Fig. 8(c), three tissue compositions of dense collagen, loose collagen, and normal myocardium are well classified and tissue compositions of dense collagen and adipose tissue are well specified in Fig. 8(f).
In this paper, we presented an automated algorithm to characterize tissue compositions in human atrial tissue. To the best of our knowledge, it is the first work to automatically classify myocardial tissue using OCT. Our algorithm builds a probabilistic model to differentiate tissue compositions in OCT statistically, and achieved a high accuracy to identify tissue composition within human atrial samples ex vivo. Our method does not take into account prior information, and thus can have potential applications to other tissues including OCT images breast,40 oral,41 skin,42,43 and vascular tissues.44,45
Our classification algorithm had an 18% false positive rate of normal myocardium being misclassified as fibrotic myocardium. The tissue regions with small diffusion of fibrosis in myocardium contributed toward this high false positive rate. Moreover, the detection rate in adipose tissue was comparatively low because the adipose tissue is located at different depths, where the texture pattern varies. In Fig. 6(a), the adipose tissue is in a honeycomb structure when it is imaged in focus. When out of focus, the adipose tissue appears as an area of isolated dots.
In the preprocessing step, we detected the surface to shift image. The image shifting, widely used in the existing segmentation method,25,26 is necessary to determine the layer boundaries. We found that the difference between estimated layer boundaries and manual segmentation results were in LA and in RA if we process our algorithm without the flattening step. In the feature extraction step, our texture analysis is based on segmented atrial layers. We thus analyze the coarseness and homogeneity features in TCN. The index-wise TCN feature could be considered if we extend our algorithm to nonsegmented images.
One limitation of our study is that the time between a donor’s death and heart imaging is variable among all samples. We found that the viability of imaging is degraded when we compared OCT images from early shipped hearts with the late shipped hearts. We will consider the deliver time as a factor to normalize the OCT image in the future; moreover, the process of fixation during histology may result in shrinkage of the cardiac sample. It can impair the accuracy of correlation between histology and OCT image. In this study, we manually matched a small number of histology and the OCT images. Serial sections of histology can ensure that a variety of tissue features are observed. However, care will need to be taken to account for the fact that multiple samples are used in a training set from a single patient/chamber. Furthermore, imaging depth has an influence on image quality of OCT B-scans. Thus, to minimize the variations of imaging depth, we (i) maintain the distance between the tissue surface and the objective to be the same range in all our experimental scenarios and (ii) adjust image contrast based on the same thresholding and histogram equalization. With the translation towards in vivo use with a forward viewing catheter probe,46,47 contact imaging will be necessary and the sample location within the image will be constant.
In the future, we will need to implement a robust algorithm that take into account arbitrary shapes in which regions of tissue may present themselves. Although most of the endomyocardial tissue shows a layered structure, we notice there is a possibility that certain tissue compositions such as adipose tissue or a mixture of myocardium and blood vessels appear in a circular shape. It is possible that such a contour would not be not identified in the column-to-column search scheme. To overcome this drawback and make our algorithm more generalized, we will have a circular identification process and a refilling process. In particular, the Hough transform can be used to find circles and the snake algorithm to delineate the circular contour. Intensity values inside the circle will be refilled with the mean value of its neighboring pixels outside the circle. Then we will process layer segmentation and the tissue classification method as proposed in this study.
An important application of this work is for the development of improved cardiac models. We will incorporate atrial tissue characterization results into sample specific electrical physiology models of the human atria. A 3-D finite element model48 will be built based on the geometry of human atria. Tissue composition is important for understanding conduction properties and arrhythmia substrates.18 A second application is towards the goal of performing an “optical biopsy.” Most biopsy samples are around 1 to 2 mm3,49 similar to the size of volumetric data we usually acquired in the OCT system. We will extend our tissue classification algorithm to a model of ventricular myocardial characterization and data acquired with high resolution OCT systems.50 Our classification algorithm has the potential to guide the procedure of endomyocardial biopsy by avoiding areas with increased scar or/and performing optical biopsies where conventional biopsy is unsafe. Lastly, our classification method can be used to aid the treatment of atrial fibrillation by radiofrequency ablation (RFA).24,51,52 In particular, the identification of tissue composition can provide guidance for ablation operation and facilitate the evaluation of ablation performance.
We have developed an image processing tool to classify tissue compositions. We proposed an automated algorithm to segment layer structures within endomyocardium. Features including optical properties, high moment statistics, and texture analysis were extracted and compared. Based on extracted features, a probabilistic model was used to identify the tissue composition of segmentation. Segmentation results from human cardiac images agreed with histology. We achieved an accuracy of 80.41% in classification. Tissue composition information provided by this method can be used for a range of applications, to further understand the role of tissue composition on the electrical function of the heart, and translational applications for the monitoring and guidance of diagnostic and therapeutic interventions.
The authors would like to thank Christopher Hermawi, Christine Fung, Theresa Lye, Nathan Lin, Xinwen Yao for their technical assistance. The study was funded in part by National Institute of Health (NIH) 1DP2HL127776- 01 (CPH), National Science Foundation funding, NSF EEC-1342273 (CPH), NSF Career Award 1454365 (CPH), SPIE educational Scholarship (YG), and Wei foundation fellowship (YG).
Yu Gan received his MS degree in communications and information systems and electrical engineering, Chinese Academy of Sciences and Stevens Institute of Technology, respectively, and his BS degree in electronic and information engineering from Nanjing University of Science and Technology. He is a doctoral candidate in electrical engineering at Columbia University.
David Tsay received his MD, PhD, and BS degrees from Columbia University. Currently, he is a cardiac electrophysiology fellow and chief innovation fellow at New York-Presbyterian Hospital.
Syed Bin Amir received his MS degree in electrical engineering from Columbia University. Currently, he is a doctoral candidate in biomedical engineering at the University of Connecticut.
Charles C. Marboe is a professor of pathology and cell biology at Columbia University Medical Center. He has 34 years of experience in cardiovascular pathology.
Christine P. Hendon received her PhD in biomedical engineering from Case Western Reserve University and her bachelor's in electrical engineering and computer science from Massachusetts Institute of Technology. Her research interests are in developing biomedical optical systems and image analysis for applications for guidance of therapy. Currently, she is an assistant professor of electrical engineering at Columbia University.