A variational method for scar segmentation with myocardial contour correction in DE-CMR images

Most automatic scar segmentation methods for cardiac DE-CMR images rely on an existing myocardial segmentation (from CINE-CMR) that is registered to the DE-CMR volume, step where alignment errors are usually introduced. We present a variational method that, with the same inputs, identifies the healthy and scarred tissue and selectively corrects the endocardial and epicardial contours. For this, we tailor an existing multiphase segmentation method to provide different regularization costs for each region, and model the data fidelity energy term with a Bayesian approach that unifies the prior tissue probabilities and the myocardial labels. Experimental results show better overlapping for the ground truth and segmented myocardium, and the segmented scar compares favorably with respect to state of the art methods.


INTRODUCTION
Delayed Enhancement (DE) Cardiac Magnetic Resonance (CMR) allows the identification of scarred tissue in the myocardium. There is evidence that DE-CMR allows to predict adverse cardiovascular events in Hypertrophic Cardiomyopathy (HCM) patients [1]. Therefore, the development of methods that accurately detect these regions is of great interest. Most of the current scar segmentation methods assume there is an available segmentation of the endocardial and epicardial contours, either manually delineated or automatically generated from an anatomical modality (usually CINE-CMR) and aligned to the DE-CMR volume to be segmented. Then, the myocardial intensity values are classified into healthy or scarred tissue either by thresholding [2] or by means of a This work was partially supported by the Spanish Ministerio de Ciencia e Innovacion and the European Regional Development Fund (ERDF-FEDER) under Research Grant TEC2014-57428-R, the Spanish Ministerio de Ciencia e Innovacion under Research Grant TEC2013-44194-P and the Spanish Junta de Castilla y León under Grant VA136U13. Thanks are also due to the Group QDiagnóstica for steady financial support. probabilistic approach [3]. More elaborate proposals include the use of advanced segmentation methods [4] or postprocessing the outcome to remove false scars arising from errors in the location of the myocardial contours [5]. The latter, however, does not modify the location of the myocardial borders. Since the myocardium is a thin structure, these errors can be significant in the computation of the scar transmurality.
In our preliminary work [6], a segmentation method for DE-CMR that allows the displacement of the CINE myocardial borders, based on a Bayesian approach that takes into account both the tissue probability distributions and the probability of the anatomical labels provided by CINE was proposed. Here, we expand it by including that Bayesian approach as the data fidelity term of an existing variational framework [7], which we have modified so that it also depends on the CINE labels, and employs different weights for the regularization of the myocardial borders and the scar contours. This approach allows the displacement of the CINE myocardial borders where the edges are clearly defined, while borders in the surrounding of locations where there is uncertainty about the myocardial border placement are maintained.

Variational Framework
the set of segmentation labels that a voxel may be classified as, A(x) : Ω ⊂ R D → A an a priori anatomical segmentation, estimated from CINE-CMR, and the set of labels that A(x) can adopt. For our specific case, we define the label sets L = {C, H, S, B} for DE-CMR and A = {C, M, B} for CINE-CMR, where C, M, H, S and B stand respectively for the blood cavity enclosed by the endocardium, the myocardium, the healthy tissue, the scar and the background.
The segmentation formulation extends the one proposed in [7], with a modified regularization term, which here contains independent cost functions for each label. The elements of the convex set can be used to express membership probability of an image location to the labels of L (soft indicator functions). A hard segmentationû(x) can be extracted from u(x) asû i (x) = 1 if i = arg max k=1,...,L (u k ) andû i (x) = 0 otherwise. The proposed primal variational problem to solve is: where, for each label, f l (x) models the fidelity of the segmentation to the imaged intensities and g l (x) weighs the segmentation total variation penalty |∇u l (x)|, which helps regularize the solution and improves boundary localization.
In DE-CMR scar identification, there is some prior knowledge about the region boundaries that may be employed to provide robustness to the segmentation. The epicardium and the endocardium without trabeculae and papillary muscles are smooth surfaces, while the scar/tissue interface can be more irregular. In addition, the endocardial and epicardial contours should not be far from their location in the CINE acquisition. For these reasons, we propose to employ the following tailored individual regularization functions g l (x): where γ 0 , γ c and γ t are scalar parameters that control the amount of regularization for the borders outside the mask, the myocardial contours, and the tissue interfaces, respectively; H(z, ε) = (1 + (2/π) arctan(z/ε))/2 is a smooth Heaviside function, where ε controls the smoothness; b(x) = |∇I(x)|; and σ b(x) is the standard deviation of b(x). The a i (x) maps are computed by smoothing each of the CINE ROI masks with a discrete truncated Gaussian filter G σ (x) and normalizing them so that with The expression for f Li (x) involves the maximum likelihood that an image location x belongs to a tissue L i subject to the given image intensity I(x) with its corresponding anatomical labels A(x) through the probabilities P (L i (x)|I(x)) (which are described in Section 2.2): since the problem is a minimization instead of a maximization. For the optimization of [7], an equivalent nonconvex dual model and its smoothed version is proposed where an iterative procedure similar to the expectation maximization algorithm is used, where at each iteration the primal and dual variables are updated by fixing their counterparts. For more details, see [7]. Here, the same methodology is employed.

Label Posterior Probability Formulation
We use the strategy of computing P (L i (x)|I(x)) presented in [6], which employs both I(x) and A(x) in its formulation and is summarized below. Applying the Bayes theorem and taking into account that both L and A are finite partitions of the label space, so that either allows for the law of total probability to be applied, we can express P (L i (x)|I(x)) as (the node location (x) has been removed for clarity): where P (A k (x)) is the probability that a node x belongs to the CINE label A k , P (I(x)|A k (x), L i (x)) is the likelihood of I(x) given L i and A k , and P (L i (x)|A k (x)) is the likelihood of L i (x) given A k (x). Next, the purpose and expression of each term is explained.
• P (I(x)|L i (x), A k (x)) is given by the intensity distributions of the tissues in L i , which are assumed to have invariant parameters with respect to the CINE anatomical segmentation A(x) and the location x; that is, P (I(x)|L i (x), A k (x)) = P (I|L i ). The Rician distribution is chosen to represent the blood and myocardial tissues, and the background is modeled by a non-parametric kernel distribution, since it is composed of several different tissues.
• P (L i (x)|A k (x)) is employed to locally weigh the influence of the CINE segmentation and the image likelihood, based on the assumption that segmentation label confusions may happen when ROI borders are not clearly defined. In that case, it is better to trust the atlas, i.e., P (A k (x)); if, on the contrary, the ROIs are easily told apart, the image likelihood should be given more consideration. Thus, three extreme situations are defined: (i) to fully trust A(x), (ii) to make the choice based only on P (I(x)|L i (x), A k (x)) around the endocardium, and (iii) the same around the epicardium. The final P (L i (x)|A k (x)) is a linear combination of these three cases, using a set of local weight maps w j (x) ≥ 0, j = 1, 2, 3, with w n (x) = (1−r(x))/2 if a n (x) > 0 and 0 otherwise for n = 1, 3, where r(x) is defined in Eq. (6); and w 2 (x) = 1 − w 1 (x) − w 3 (x). The full details are found in [6].

EXPERIMENTAL RESULTS
For this work, 11 CMR studies from HCM patients were used, each of which contained CINE sequences in short axis (SAx), two chamber (2C) and four chamber (4C) long axis (LAx), and a SAx DE-CMR sequence. Six of them had scarred tissue in the myocardium. All sequences were acquired with a 3T Philips Achieva MR scanner, and their main acquisition parameters are summarized in Table 1. Both in CINE and DE-CMR, the in-plane spacing was the same in both slice axes, and there were no gaps between slices. Expert delineations for the endocardium and the epicardium were manually drawn on the SAx DE-CMR volume and the SAx CINE in end-diastole. Along with the myocardial contours of the DE-CMR volume, the cardiologists also provided a validated threshold for the scarred tissue within the myocardium. In hyperenhanced regions considered as artifacts by the cardiologists, the SAx DE-CMR myocardial contours were adjusted in order to exclude them from the scar segmentation; these regions are not taken into account in the metrics, for any of the segmentation methods. The CINE and DE-CMR volumes were spatially aligned using the framework described in [8]. Summarizing, the SAx volumes are interpolated in the long axis direction in order to provide quasi-isotropic resolution. Then, the respiratory motion shift between slices in the SAx CINE is corrected using the 2C and 4C LAx CINE sequences, and after that, the corrected SAx CINE and the DE-CMR volumes are aligned by means of a registration using the Mattes mutual information as metric [9]. All manual delineations and the masks orig-inated from them were subjected to the same spatial transformations as their CMR counterparts. Then, the DE-CMR volumes were segmented using the threshold-based method with a feature analysis false positives and false negatives removal (TFA1) proposed in [5]; a similar method that also employs thresholding and feature analysis, proposed in [10] (TFA2); the method (PWD) described in [3], which performs a watershed segmentation using probability distribution models and distances to the endocardium; and the proposed variational method (PROP). Small scar islands are removed and scar holes are filled in the watershed output. For TFA1, the maximum size for a scar island rejected as a false positive was set to 24 voxels. For the proposed method, ε = 0.1, σ = 2, γ 0 = 3, γ c = 0.5 and γ t = 5 · 10 −4 were chosen.
The Dice Index (DI) is used to assess overlapping between regions of interest (ROIs) in the segmentations and the ground truth. The DI expression is: DI = 2|GT ∩ S|/(|GT | + |S|), where S and GT respectively stand for the segmented ROI and the ground truth against which it is compared. It must be considered, however, that the ground truth for the scar is taken from a thresholding, with the drawbacks it entails (possibility of false positives and false negatives). In Table 2, the DI values of the scar ROIs yielded by the considered methods are shown for all the studies where there was scarred tissue in the ground truth (1)(2)(3)(9)(10)(11). For studies 10 and 11, all DIs are null and not shown in Table 2. This is due to the fact that all methods failed in obtaining a good segmentation criterion: either because a too high threshold was computed in TFA1 and TFA2, or because the distribution estimation was not accurate, in the rest of the methods. The best DI for each study is in bold. The DI for study 9 and method PWD could not be computed because the distribution parameter estimation crashed. The proposed method achieved higher DI values than the methods without contour correction, except for TFA2 in study 9, where it was second best. In order to assess the contour correction effect of the PROP method over the aligned CINE, the DI between the ground truth myocardium and the myocardial masks yielded by the respective methods was computed, and then, a Bland-Altman plot was generated by plotting, for each study, the averaged PROP myocardium and aligned CINE DIs against their difference, which is shown in Figure 1. We observe that in general, the PROP method obtained a higher DI than the aligned CINE, and that the cases where the DI of the aligned CINE was lower, the increase in DI yielded by the contour correction is higher. Figure 2 shows an example of the aligned CINE contours and the endocardial and epicardial contours using the proposed method over a checkerboard of the DE-CMR and the CINE images, which also allows to assess the alignment between sequences, in a case where the probability maps show irregularities on the epicardial border. We observe that the myocardial contours provided by the proposed method have been smoothed by the regularization included in the variational framework.

CONCLUSIONS
We have described a segmentation method for cardiac DE-CMR, where, in addition to identifying the scar, the myocardial borders may also be modified. To do this, a variational method where the data fidelity uses a Bayesian approach that takes into account both the image intensity probability distributions and a registered myocardial segmentation coming from CINE. We have observed that the correction of the myocardial contours improves the scar identification, and that the correction is stronger in the volumes with lower CINE myocardial alignment.