Probability mapping of scarred myocardium using texture and intensity features in CMR images

Earlier reported work

To our knowledge, visualization of the heterogeneous nature of the myocardium is mostly based on thresholding methods described in the literature where thresholds are defined at intensity levels corresponding to some percentage of the maximum intensity level inside the scar area [1–3]. Recently, Shokrollahi et al. was able to divide the scarred myocardium of animals into core and gray zones using fuzzy clustering algorithm and validated the results using in-vivo electro-anatomical voltage mapping [6]. Some previous publications reported automatic segmentation of scarred myocardium in late enhanced magnetic resonance images [7–9]. The segmentation of scar itself is an important application as the scar size contains information useful for finding the patients with high and low risk of arrhythmia [3]. Also, scar size is largely responsible in ventricular remodeling [10]. Dicki et al. use support vector machine (SVM) to classify the scar and the myocardium in delayed enhancement magnetic resonance images in [7] using three features. The three features used in [7] depend on the intensity of scar, myocardium and the whole CMR image.

From our group’s previous work [11], it is shown that patients with high and low risk of getting arrhythmia can be classified using textural measures from gray-level co-occurrence matrices. In our recent work [12], we presented some preliminary results where dictionary learning and sparse representations were used for representing the texture in the myocardial segments followed by a classification.

The proposed probability mapping

All the above discussed previous work have sought to identify methods to segment the hyper-enhanced area of the myocardium. However, by only delineating the hyper-enhanced myocardium, important information concerning the tissue within the infarct territory is not provided. Hence, our focus is on developing a method that provides probability values as a means to distinguish the areas of the myocardium rather than the crisp definitions offered by a scar vs no scar segmentation.

In this work, we propose a technique for transforming the scarred myocardium in the CMR image to a probability map to indicate varying degrees of its functioning capabilities: where high and low probabilities will indicate the scarred tissue and normal tissue, respectively. Section ‘Calculation of probability maps using Bayes rule’ describes the general principles for calculation of the probability maps, where each pixel is the posterior probability of that pixel being scar compared to healthy myocardium. The posterior probabilities are calculated using Bayes rule with the help of features from the myocardium. Features are calculated pixel by pixel based on a local neighborhood and is explained in Section ‘Extraction of features’. We have based our experiments on two features; a local mean (direct current, DC) that will correspond well to how the doctors manually segment the scar, and a textural feature using sparse representations, that might reveal information regarding gray zone areas.

Section ‘Experiments and results’ shows the experimental framework and examples of probability maps for both features. Also in section ‘Experiments and results’, cardiac segments example is demonstrated to show the diagnostic importance of texture feature probability maps. Finally, the paper is concluded in section ‘Conclusion’.

Training and testing

The process of generating probability maps is depicted in Figure 3. The whole CMR data set is divided into training and testing images. In the training process, training features are extracted along with some training parameters required in the testing phase. In the next step, the calculated training feature vectors are used for estimating the PDFs of scar and myocardium region using ML. The training process and its outputs are indicated by the dashed lines. The probability maps are generated on test images, different from the training images, and this process is indicated by solid lines. The extraction of features for training and testing phases are explained in section ‘Extraction of features’.

DC Feature

Historically, in electronics field, the mean is commonly referred to as DC (direct current) value [15]. Because of the intensity differences between the healthy and scarred myocardium tissues in the LG enhanced CMR images, DC value dc(i,j) is one of the features explored for developing the probability maps. This is probably the feature that would mimic the cardiologists segmentation best since they use the intensity values of scarred myocardium in CMR image to segment the scarred areas manually. We define the feature $\mathit{\text{dc}}(i,j)=\mathit{\text{mean}}\left(I_{\sqrt{N}\times \sqrt{N}}\right(i,j\left)\right)$ (where $I_{\sqrt{N}\times \sqrt{N}}$ is the neighborhood around the pixel I(i,j)), such that a new image I _dc , with dc(i,j) as values at pixel position (i,j) is made with a sliding window averaging over the image.

Dictionary-based textural feature

The pseudo-norm ∥.∥₀ is the number of nonzero elements and, s denotes the sparseness. A vector selection algorithm called Order Recursive Matching Pursuit (ORMP) algorithm is used for sparse approximation [16].

This is a very hard optimization problem. Recursive Least Squares Dictionary Learning Algorithm (RLS-DLA), [17] is used for dictionary learning in all our experiments. A brief explanation of the algorithm is given in Appendix.

Sparse representations and learned dictionaries have been shown to work well for texture classification by Skretting and Husøy in [18] and by Mairal et al.[19]. Texture in a small image patch can be modeled as sparse linear combination of dictionary atoms in Frame Texture Classification Method (FTCM) presented by Skretting et al. in [18]. FTCM is developed by modeling a texture as a tiled floor where all the tiles are identical. The color or the gray level, at a given position in the floor is given by an underlying continuous periodic two-dimensional function. It is shown based on this model that a vector from spatial neighborhood is indeed a sparse combination of finite dictionary atoms.

where $w_l^s$ and$w_l^m$ are sparse coefficient vectors.

Pixel by pixel R _p  value can be interpreted as the scaling of residuals R _s  and R _m . Smaller values of R _p  means that the pixel is not likely to be scar (i.e. healthy myocardium), and larger values means that it is likely to be scar or border area. The texture feature, R _p (i,j) from the training set is used to estimate the PDFs for the probability maps computation. In the testing phase, texture feature vectors are collected from the test image set, in the same way as training feature vectors. In Figure 3, training and classification parameters are indicated. In fact, these training parameters are the dictionaries, D _s  and D _m . Using D _s  and D _m learned in the training phase, the residual images, R _s  and R _m  are calculated for test images. Texture is not pixel-by-pixel local except in the edge between two textures; thus some sort of smoothing is used on the test residual images before calculating the probability maps [18]. The test residual images R _s , and R _m  are smoothed using an a × a pixels separable Gaussian low pass filter with variance σ² before calculating R _p .

Probability maps obtained using DC feature, dc(i,j)

The neighborhood size 3×3 was used to form training vectors as explained in ‘DC Feature’. The same neighborhood size must be used while training and finding the DC images I _dc . The DC images obtained from the training images were used to form the training feature set. The parameters of the class specific PDFs: the mean and the standard deviation were found using ML estimation using the training feature set. DC-values were scaled to have zero mean and unit variance before finding the ML estimates. The scaling coefficients from the training were stored to scale the test vectors. The probability maps of the scarred myocardium were calculated using Bayes rule as explained in section ‘Calculation of probability maps using Bayes rule’. An example of a probability map obtained using the DC feature with color code is shown in Figure 2(b). Note the sharp transition from the segmented scar to the myocardium.

Probability maps obtained using texture feature, R _p (i,j)

The training vector and test vectors were generated in the same way as in the DC feature experiment using the same neighborhood size 3×3. Consider the pixels on the border zones, their neighborhood extends into other regions that are not under consideration. If we use training vectors from border regions, the dictionaries might learn the texture properties of other regions along with the texture properties they are intended to learn. So, the training vectors for the pixels whose neighborhood span other regions were not considered in our experiments. This is depicted in Figure 4.

The dictionaries were learned on the training vectors generated from scar and healthy myocardium using RLS-DLA. The dictionary size of 9 × 90 atoms was used in our experiments. The initial dictionaries were formed by randomly selecting 90 vectors of length 9 from the training sets. Dictionary size was selected based on previous experience, and it should be K > N < < L. The forgetting factor is initialized to 0.995 and slowly increased towards 1 according to the exponential method described in [17]. If the number of dictionary atoms used to represent the image patch increases, then the residual decreases, and a large number of atoms will provide a good approximation with any full rank dictionary. Therefore, higher sparsity lowers the difference between the residuals of healthy myocardium and scar which in turn decreases the differentiation between the classes. The sparsity s used in our work is two, i.e. we used two vectors from the dictionary to represent the original image patch.

After learning the dictionaries, they were used to obtain the residual images. For each pixel in the myocardium of training images, the scaled residual R _p (i,j) was found from residual images as explained in section ‘Dictionary-based textural feature’. The training feature set was then generated from these scaled residual values. The parameters of the PDFs to be used with the classifier was estimated from the labeled training vector set using ML estimation. The test vector set in this experiment was calculated after smoothing (using low pass Gaussian filter with σ = 5 and 9 × 9 window size) the test residual images. Using the ML estimates and Bayes rule, the probability maps for texture features of the scarred myocardium were obtained. An example of such a probability map with color code is shown in Figure 2(c). The color code used is explained under Figure 2. Here, the dictionaries were learned without the removal of the DC value in each image patch.

Dictionary learning was performed with and without the removal of the DC factor from the image patch of the training set. In order to illustrate some of the texture captured by the dictionary atoms, they were plotted as image patches in Figure 5. The difference between the scar and the healthy myocardium dictionaries was visible when the DC value was not removed, and this is according to the known intensity difference between scar and healthy myocardium. To investigate the variation in the dictionary atoms used as image patches, we had calculated a measure from the first derivative of dictionary atoms in both vertical and horizontal directions. The absolute mean of the first derivative of the dictionary atoms was used to compare the scar and healthy myocardium dictionaries learned with and without DC value. The calculated absolute means are shown under Figure 5. The difference between the scar and healthy dictionaries without DC removal was observed clearly, both from the depicted dictionaries and the absolute means. The difference between the scar and healthy myocardium dictionaries trained by removing the DC value from image patches could not be observed visually. A larger block for DC removal prior to textural capturing is probably more appropriate, more similar to suppression of difference in illuminance conditions in computer vision.

Figure 6 shows probability maps of some example CMR images and all of them were calculated using the texture feature, R _p (i,j). In the middle column, the DC value was not removed from each image patch before training dictionaries and calculating the probabilities. In the last column, the DC value was removed from every image patch both before learning dictionaries and before finding probabilities. The DC value is obviously important, so we do not expect this experiment to provide competitive results, but we see a tendency of higher values in the scar and border areas, but not necessarily in the core of the scar. We expect the scar core to be very bright and could be captured easily by some sort of intensity based measure (for example the DC feature). If the physical nature of the tissue is reflected in the image texture, we can expect the scar core to be quite homogeneous. Thus, we are not surprised that the regions that are the brightest in column II were actually quite dark in column III (core), but the intermediate regions were captured fairly. We have continued with the experiments in this work by letting the dictionaries incorporate the DC value as we did not combine the DC feature, dc and texture feature, R _p  into a 2D dimensional feature for computing the probability maps. Unless it is mentioned under figures of probability maps, all the probability maps of texture feature, R _p  were computed with dictionaries learned without removing DC from every image patch.

Examples of probability maps obtained from DC and texture feature

The second and the third column of Figure 7, shows more examples of probability maps of scarred myocardium using the DC, dc(i,j), and texture, R _p (i,j), features, respectively. Figure 8 shows examples of probability maps of scarred myocardium for both DC and texture feature for the patients without ICD. This shows that the PDFs estimated for the DC and texture features from the patient group implanted with ICD were also able to produce probability maps for low risk arrhythmia patients (i.e. patients without ICD). Both the features were able to identify the healthy myocardium without scars, and an example of this is shown in Figure 9. In order to compare the probability maps obtained by the two methods, they are plotted with the same color map. Experiments of probability mapping combining the DC feature and texture feature gave almost identical probability maps as the DC feature because of the dominance of the DC feature over the texture feature. Improved ways of combining both the features will be explored in our future work.

A possible application - cardiac segments

A cardiac segment, CS can be defined as a region of interest in the myocardium which might have diagnostic importance. We define a candidate cardiac segment by a lower and upper probability so that the segment consists of all pixels with probability values in the range between these. Different candidate cardiac segments might be evaluated, for example, calculating one or several features from the candidate segment and correlate result to a clinical meaningful hypothesis. We demonstrate the concept by making a comprehensive search for significant cardiac segments by comparing two patient groups and testing for statistical significance. The demonstration experiment is meant to illustrate how the probability mapping can be used to localize cardiac segments containing information discriminating between the high and low risk patient groups.

For each patient, cardiac segments were determined by a lower boundary L and an upper boundary U. These boundaries were varied in the range 0–1 in steps of 0.025 with the restriction that U - L ≥ 0.1. Figure 10(a) illustrates this by visualizing each of the resulting candidate cardiac segment as a black dot with coordinate (L,U). For example, the three yellow dots in coordinates (0.15,0.3), (0.3,0.7), and (0.7,0.825) correspond to the cardiac segments CS₁ = 0.15 - 0.3, CS₂ = 0.3 - 0.7, and CS₃ = 0.7 - 0.825, respectively. Figure 10(b) and (c) show one of the original image slices and the corresponding probability mapping (texture feature) for one of the high risk patients. The probability maps, which do not have 0 - 1 probability range, was extended using a sigmoid function for easy comparability. Figure 10(d–f) shows the cardiac segments corresponding to CS₁, CS₂ and CS₃ for the image slice shown in Figure 10(b).

Preliminary experiments indicate that corresponding segments determined without using the probability mappings depends on the signal intensity values and seems to lack the ability to define contiguous areas in the same way as the probability mapping for some of the cardiac segments. For example, the cardiac segment determined directly from thresholding relative to the maximum 3D scar signal intensity is shown in Figure 10(g–i) for boundaries corresponding to CS₁CS₂ and CS₃, respectively. These plots are comparable to the corresponding cardiac segments based on the probability mapping in Figure 10(d–f), and show less contiguous cardiac segment, especially CS₁. Figure 11 shows similar plots of the cardiac segments corresponding to CS₁, CS₂ and CS₃ for one of the low risk patients.

For each candidate cardiac segment, its relative size (accumulated for all slices) to the total myocardium was calculated. The computed relative sizes for high and low risk patient groups were compared using a Mann-Whithney test at significance level 0.05. The LU-plots in Figure 12 show the p-values for the test of all the candidate cardiac segments derived using the texture based (a) and DC based (b) probability map and those derived without the use of probability mapping (c). The blue colored dots indicate cardiac segments for which the high and low risk groups showed significant differences of the relative size.

Interpreting the comparison results for the texture feature probability map based cardiac segments in Figure 12 (a), the significant CS with L = 0 and U > 0.5 correspond to the cardiac segment containing all of the healthy myocardium and gradually including the scar area. More interesting is the localization in the upper right corner which makes more specific localization in the myocardial tissue in the regions close to the scar. Note that CS₃ is identified as significant while CS₁ and CS₂ are not. Figure 12 (b) shows the LU-plot from the probability maps of DC feature and, the plot shows that cardiac segments does not give any localization. This is due to the narrow transition between high and low probabilities of the DC based probability maps. The DC probability maps will predominantly have extreme values close to 0 or 1 and few values in between 0–1. This provides little information to explore. These plots (Figure 12 (a) and Figure 12 (b)) of cardiac segments clearly shows that the probability maps obtained from the DC features reflect the way cardiologist perceive the scarred myocardium in CMR images and, the probability maps from texture features gives information about the characteristics of myocardial tissue and might emphasize information that is not easy for inspection with a human eye in the original CMR image. The difference between the narrow and broad transitions for the DC- and texture-based probability mappings might explain this. As the narrow transition gives little information besides the crisp core/scar discrimination corresponding to the cardiologist’s delineation based on the visual information. Thus, most of the candidate cardiac segments will be empty as there are virtually no intermediate probability values. This is not the case for the texture based probability mappings, which discloses the information present in this intermediate probability range.

In comparison the results for the cardiac segments derived directly without using the probability plots (Figure 12 (c)), the localization seems vaguer, not specifically identifying any coherent region in the myocardium. It is our belief that probability mapping might be used to identify cardiac segments that might be used in future studies.

The discriminatory power of the DC and texture features

The probability maps calculated for the two features, DC and textural feature, give different information. The discriminatory power of the two features were explored by comparing to the manual segmentation of scar and healthy myocardium to obtain Receiver Operating Characteristic (ROC) curves. The ROC curves calculated for the DC and texture features on the CMR images of the 18 test patients with high risk of getting arrhythmia are shown in Figure 13. The area under the average ROC curves is 0.9052 and 0.8428 for the DC and texture features, respectively. From ROC curves, it is observed that the discriminative power of the DC feature is comparatively higher than that of the texture features. The sensitivity, and specificity reported in [7] are plotted as a single point in Figure 13. This reported result is not calculated on our CMR data. This single point lies below within 95% confidence interval of DC feature. This reported segmentation result performed well compared to our average DC and texture ROC curves as Dikici et al.[7] does not include all the CMRI slices of a patient for training and testing where as in our work we include all the CMRI slices i.e. the scar volume.

Recursive least squares dictionary learning algorithm

The dictionary learning problem can be solved in a two steps algorithm. Step 1) Find sparse coefficient matrix, W, using vector selection algorithms, keeping dictionary, D, fixed. Step 2) Keeping W fixed, find D. The dictionary update step depends on the dictionary learning method we choose to use. D and C are initialized randomly.

where u = C_t-1w _t  and $\alpha =1/(1+w_t^{T}u)$, r _t  = x _t  - D_t-1w _t  is the approximation error.

where $u=\left({\lambda }_t^{-1}{C}_{t-1}\right)w_t$. The remaining equations remain the same. Due to introduction of the forgetting factor, RLS-DLA has good converging properties and is less dependent on the initial dictionary.