Fast Automatic Myocardial Segmentation in 4D cine
CMR datasets
Sandro Queirósa,b,1,∗, Daniel Barbosaa,c,d,1 , Brecht Heydea , Pedro Moraisa,b ,
João L. Vilaçab,d , Denis Fribouletc , Olivier Bernardc , Jan D’hoogea
a
Lab on Cardiovascular Imaging and Dynamics, KU Leuven, Belgium
ICVS/3B’s - PT Government Associate Laboratory, Braga/Guimarães, Portugal
c
Université de Lyon, CREATIS, CNRS UMR5220, INSERM U630, Université Lyon 1;
INSA-Lyon, France
d
DIGARC - Polytechnic Institute of Cávado and Ave (IPCA), Barcelos, Portugal
b
Abstract
A novel automatic 3D+time left ventricle (LV) segmentation framework is
proposed for cardiac magnetic resonance (CMR) datasets. The proposed
framework consists of three conceptual blocks to delineate both endo and
epicardial contours throughout the cardiac cycle: 1) an automatic 2D midventricular initialization and segmentation; 2) an automatic stack initialization followed by a 3D segmentation at the end-diastolic phase; and 3) a
tracking procedure. Hereto, we propose to adapt the recent B-spline Explicit Active Surfaces (BEAS) framework to the properties of CMR images
by integrating dedicated energy terms. Moreover, we extend the coupled
BEAS formalism towards its application in 3D MR data by adapting it to a
cylindrical space suited to deal with the topology of the image data. Furthermore, a fast stack initialization method is presented for efficient initialization
and to enforce consistent cylindrical topology. Finally, we make use of an
anatomically constrained optical flow method for temporal tracking of the
LV surface.
The proposed framework has been validated on 45 CMR datasets taken
from the 2009 MICCAI LV segmentation challenge. Results show the robustness, efficiency and competitiveness of the proposed method both in terms
Corresponding author
Email addresses: sandroqueiros@ecsaude.uminho.pt (Sandro Queirós),
dbarbosa@ipca.pt (Daniel Barbosa)
1
Joint first authorship.
∗
Preprint submitted to Medical Image Analysis
June 22, 2015
of accuracy and computational load.
Keywords: fast image segmentation, cardiac cine MRI, left ventricle
segmentation, automatic initialization
1. Introduction
Cardiovascular diseases (CVDs) are the leading cause of death in the
world. According to the World Health Organization (2012), an estimated
17.3 million people died from CVDs in 2008, which represents 30% of all
global deaths. This report and the known global population aging highlight
the increasing need for powerful and efficient techniques for diagnosis and
treatment follow-up of patients with CVDs. Among others, assessment of
the left ventricular morphology and global function using non-invasive cardiac
imaging is an interesting technique for such purpose.
Cardiac magnetic resonance (CMR) imaging is a successful and promising modality that has already demonstrated its excellent accuracy and reproducibility (Miller et al., 2012a). The quantification of left ventricular
volumes, mass (LVM) and ejection fraction (EF) is crucial for cardiac diagnosis, requiring the delineation of endocardial and epicardial contours of the
left ventricle (LV) from cine MR images. In clinical practice, manual segmentation is usually performed, being a tedious, time consuming and unpractical
task (Petitjean and Dacher, 2011), normally requiring between 6 to 20 minutes to segment all slices covering the whole LV in end-diastolic (ED) and
end-systolic (ES) cardiac phases only (Miller et al., 2012b; Hautvast et al.,
2012; Petitjean and Dacher, 2011). Moreover, manual delineation is prone to
intra- and inter-observer variability, mainly associated with the choice of ED
and ES phases, the choice of the most basal LV slice and the chosen endocardial edge selection approach (Anderson et al., 2007; Miller et al., 2012b). In
their recent study, Miller et al. (2012b) presented the “real-world” variability
associated with manual contouring, showing that significant differences in the
quantification of LV indices occur depending on the used analysis method. In
practice, such differences would require specific reference ranges depending
on the contouring methodology used by the physician.
To date, several methods have been presented for (semi-)automated LV
CMR segmentation. According to the recent review in Petitjean and Dacher
(2011), these methods can be divided in two main categories: weak or no
prior methods and strong prior methods. The first uses weak assumptions,
2
such as spatial or intensity-based relations and anatomical assumptions, and
include image-based techniques (threshold, dynamic programming, etc.) (Hu
et al., 2012; Liu et al., 2012), pixel classification methods (clustering, gaussian
mixture model fitting, etc.) (Ulén et al., 2013; Pednekar et al., 2006) and
deformable models (active contours, level-sets and variants) (Chen et al.,
2008; Pluempitiwiriyawej et al., 2005; Paragios, 2002; Lynch et al., 2008;
Constantinides et al., 2012; Cordero-Grande et al., 2011; Ben Ayed et al.,
2012). Strong prior methods include shape prior based deformable models,
active shape and appearance models and atlas based methods, which focus on
higher-level shape and intensity information and normally require a training
dataset with manual segmentations (Ulén et al., 2013; Eslami et al., 2013).
An extensive overview of the different available methodologies can be found
in Petitjean and Dacher (2011).
Although several approaches are already automatic, optimal boundaries
assessment from base to apex is still lacking, usually requiring the physician
to manually correct the contours (Hautvast et al., 2012; Petitjean and Dacher,
2011). Therefore, the LV segmentation remains an open problem with significant clinical value, easily noticed by the constant attention in the research
field (Cordero-Grande et al., 2011; Constantinides et al., 2012; Hu et al.,
2012; Liu et al., 2012; Ben Ayed et al., 2012; Ulén et al., 2013; Eslami et al.,
2013). In 2012, Hu et al. (2012) presented an image-based framework using
a gaussian-mixture model for endocardial delineation and region restricted
dynamic programming for epicardial contouring. Similarly, Liu et al. (2012)
also used region restricted dynamic programming for epicardial contouring,
but opted for topological stable-state thresholding to delineate the endocardium. Both methodologies are fast, automatic and proven to be robust
in a publicly available database. In the same year, Ben Ayed et al. (2012)
presented a fast LV segmentation framework based on max-flow optimization of discrete cost functions encompassing global intensity and geometry
constraints based on the Bhattacharyya similarity. Although sufficiently fast
in segmenting all cardiac phases, it still required a manual initialization in
the first frame. More recently, in 2013, Ulén et al. (2013) presented a multiregion segmentation framework based on graph cuts and Lagrangian duality.
They evaluated it for the segmentation of cardiac structures in CMR, namely
LV and right ventricle (RV) cavities, myocardium and trabeculae and papillary muscles (TPMs). Eslami et al. (2013) presented a framework based on
guided random walks for LV segmentation with high-level prior knowledge
but without the need for statistical shape models. Instead, they guide the
3
segmentation using the closest subject in the database. Both Ulén et al.
(2013) and Eslami et al. (2013) require some user interaction and most importantly the incorporation of training datasets, which are normally efficient
but intrinsically link their ability to handle clinical variability to the amount
of different morphological and physiological cardiac patterns included in the
training set.
The effort seen in the research field reinforce the difficulties in designing
a framework that could deal with the segmentation challenges presented in
CMR images (Petitjean and Dacher, 2011; Cordero-Grande et al., 2011; Chen
et al., 2008), such as:
• the poor contrast between myocardium and surrounding structures
(e.g., lungs and liver);
• the heterogeneities in the LV cavity brightness due to blood flow;
• the presence of TPMs with similar signal intensities to the myocardium;
• the partial volume effects (PVE) in apical slices due to the limited MRI
resolution;
• the inherent noise associated with cine CMR, which is caused by motion
artifacts, heart dynamics and intensity inhomogeneity;
• and the considerable shape and intensity variability of the LV across
patients, notably in pathological cases.
Thus, to be of interest in “real-world” clinical practice, a robust segmentation framework should be accurate, flexible and generic enough to be able
to segment a multitude of LV shapes with various myocardial contraction
patterns. Moreover, the time required to segment is a critical factor for
success when introducing a methodology in daily clinical practice.
Among the multitude of methodologies, active contour-based methods are
a flexible framework, allowing easy incorporation of prior knowledge about
the object’s shape and expected intensity distribution, as well as of temporal
constraints (Petitjean and Dacher, 2011). In these methods, the shape to
recover is captured by propagating an evolving interface while minimizing an
energy functional that reflects the object’s properties (Barbosa et al., 2012).
As image feature, the energy terms are usually classified as either edge-based
or region-based. Edge-based energy terms use image gradients to identify
4
object boundaries and attract the contours toward edges (Paragios, 2002;
Pluempitiwiriyawej et al., 2005; Chen et al., 2008), while region-based ones
use regional properties inside and outside the object to evolve the contour
(Paragios, 2002; Pluempitiwiriyawej et al., 2005; Chen et al., 2008; Lankton
and Tannenbaum, 2008; Barbosa et al., 2010). Some authors combine both
to improve accuracy, while adding smoothness constraints, contour interaction and/or shape priors (Paragios, 2002; Pluempitiwiriyawej et al., 2005;
Chen et al., 2008; Ben Ayed et al., 2009). Moreover, temporal integration
is sometimes included for 2D+time or 3D+time segmentation (Lynch et al.,
2008; Ben Ayed et al., 2012; Hautvast et al., 2006). Note that active contours
are usually dependent on their initialization, highlighting the importance of
this preliminary step (Petitjean and Dacher, 2011). Thus, active contours
usually require some user interaction, either in a minimal way (Paragios,
2002; Pluempitiwiriyawej et al., 2005; Chen et al., 2008) or in an advanced
way (e.g., requiring a first manual contour) (Ben Ayed et al., 2009, 2012;
Hautvast et al., 2006). Nevertheless, several authors proposed automatic
initialization procedures, focusing on obtaining a robust initialization while
avoiding the variability introduced by user input (Lu et al., 2009; Ciofolo
et al., 2008; Queirós et al., 2013).
The focus of the present work was to develop an automatic 3D+time
myocardial segmentation framework for CMR datasets with accurate shape
delineation and low computational burden, to be of interest in daily clinical
practice. To this end, we propose a novel framework corresponding to the
3D+time extension of the preliminary work previously presented in Queirós
et al. (2013). This work encompassed a fast automatic initialization approach
capable of detecting LV elliptical annular shapes and an adapted version of
the B-spline Explicit Active Surfaces (BEAS) framework, initially proposed
in Barbosa et al. (2012), to 2D CMR image segmentation (Queirós et al.,
2013). In the present paper, the key novelties introduced are three-fold:
• First, we extend the coupled BEAS framework introduced in Queirós
et al. (2013) towards its application in 3D MR data. To this end, we
adapt the mathematical formalism to a cylindrical space suited to deal
with the topology of the image data. Furthermore, we provide a rigorous mathematical derivation to support the variational formulation of
the coupled segmentation problem.
• Secondly, we develop a fast stack initialization method in order to provide an efficient and spatially consistent initialization for the 3D seg5
mentation algorithm. This step also enables to model the 3D LV surface
in a parametric cylindrical space, since it partially compensates for slice
misalignment. This is achieved through sequential threshold-based segmentation of adjacent slices, using the information gathered from the
automatic segmentation of the mid-ventricular short-axis (SAX) slice.
• Lastly, we introduce a hybrid energy term making use of both region
and edge-based image features. Although such strategy has been previously applied to level-set frameworks (Lin et al., 2003; Paragios, 2002),
it is the first time that edge-based terms are used in the BEAS formalism.
Additionally, the LV surface was tracked throughout the cardiac cycle
making use of an anatomical affine optical flow method, initially presented
for ultrasound (US) data (Barbosa et al., 2013b).
The accuracy of the proposed framework was assessed against manually
extracted references for 45 publicly available SAX datasets from the 2009
MICCAI LV segmentation challenge, including both healthy subjects and
patients. The results were compared to state-of-the-art methods that used
the same database and a very competitive performance was achieved by our
algorithm.
This paper is structured as follows. In Section 2, we focus on the detailed description of the proposed segmentation framework. In Section 3,
we discuss the implementation details of our method, particularly the choice
of parameters and their influence. Section 4 presents the validation experiments, whose results are presented in Section 5. In Section 6, we evaluate
and discuss the performance of each module of our framework and compare
with state-of-the-art methods. The main conclusions of this paper are given
in Section 7.
2. Methodology
2.1. General Overview
The ultimate goal of the proposed approach is to obtain an accurate
and fast 3D+time LV segmentation framework for CMR datasets in order
to assess both cardiac morphology and global function. We thus propose to
divide the approach in 3 conceptual blocks (Figure 1).
The first block corresponds to an automatic 2D mid-ventricular segmentation in the ED phase, divided in two modules. The first module focuses
6
Figure 1: General overview of the proposed fully automatic 3D+time LV segmentation
framework for CMR datasets.
on obtaining an initial contour for a mid-ventricular slice. Hereto, we use
the fast and robust automatic initialization procedure presented in a previous contribution (Queirós et al., 2013). The second module uses the BEAS
framework and a specifically designed energy functional for 2D segmentation
of the mid-ventricular slice (whose expression is given in Section 2.3).
The second block presents an extension to 3D segmentation of BEAS, with
a hybrid region-based and edge-based epicardial energy functional. However,
as it requires all slices in the ED phase to be initialized, an automatic stack
initialization is first performed using a simple threshold-based BEAS formulation using information from the 2D segmentation result (i.e., first block).
This step also allows to estimate possible slice misalignment to achieve spatial continuity during the 3D segmentation, as further explained in Sections
2.4 and 2.5.
At last, the third block focuses on tracking the segmented contours of the
ED phase over all cardiac phases, ultimately obtaining both endo and epicardial contours in the ES phase. Hereto, we make use of a global anatomical
affine optical flow method (Barbosa et al., 2013b), which allows performing
in-plane tracking of the myocardial motion during the cardiac cycle.
Using the resulting delineations, one is able to compute several clinical
cardiac indices, such as end-diastolic volume (EDV), end-systolic volume
(ESV), stroke volume (SV), LVM and EF, in an automated way.
2.2. Automatic Initialization using an Elliptical Annular Template Matching
Algorithm
In Queirós et al. (2013), we introduced a new automatic initialization
algorithm for a mid-ventricular slice based on a LV localization method and
an elliptical annular template matching algorithm.
7
In brief, the method detects the LV centroid by applying a multilevel
Otsu thresholding (Otsu, 1975) to divide the image in four classes (blood,
soft tissues, background and one related to bright artifacts), followed by a
class decomposition and a search procedure among the candidate objects.
To this end, the 4 classes are used to create 3 independent images to be
analyzed (Figure 2B), taking into account their possible distribution in the
presence or absence of artifacts: an image with the brighter class alone (since
the blood pool is usually the brightest object in the image); an image with
the two brightest classes combined (in the presence of artifacts or due to
heterogeneities related to blood flow, the cavity can be divided between these
classes); and an image with only the third brightest class (in case strong
artifacts occupy the two brightest classes). The LV object is chosen as the
most circular and close to the image’s center among the three analyzed images
(Figure 2C) (Queirós et al., 2013).
In a second step, in order to obtain two elliptical initial contours to be
used in the segmentation procedure, a template matching method designed
to search for a dark elliptical ring is applied (Queirós et al., 2013). The
strategy is to find the template (by varying the ring’s major axis, thickness,
orientation and eccentricity in a total of 240 possibilities - Figure 2D) that
best fits the subject-specific myocardium (Figure 2E). Note that to reduce the
search space, the object detected in the first step is used to fit an ellipse and
estimate the major axis and orientation, which defines the central value of the
search range for these parameters. For every kernel, a block matching search
is performed using normalized cross correlation and the position maximizing
the cross correlation value defines the optimal center for the kernel under
evaluation. Finally, the kernel having the maximum cross correlation value
among all kernels analyzed defines the optimal parameters, used to draw two
ellipses to be used as initial contours (Figure 2F-G) (Queirós et al., 2013).
2.3. Coupled Myocardial 2D Segmentation
2.3.1. B-Spline Explicit Active Surfaces
The key concept of the BEAS framework is to regard the boundary of
an object as an explicit function, where one of the coordinates of the points
within the surface, x = {x1 , · · · , xn }, is given explicitly as a function of the
remaining coordinates, i.e.x1 = ψ(x2 , · · · , xn ). In this framework (Barbosa
8
Figure 2: Automatic initialization algorithm based on LV localization and elliptical annular
template matching algorithm. a) Image ROI; b) class decomposition for LV object search
(from left to right and top to bottom: multilevel Otsu thresholding, brightest class, two
brightest classes combined, third brightest class); c) LV object and its centroid; d) several
constructed kernels for template matching with different orientations, radii, standard deviations and eccentricities; e) block matching strategy using normalized cross correlation,
f) original image overlay with optimal kernel; g) resulting elliptical initial contours.
et al., 2012), ψ is defined as a linear combination of B-spline basis functions:
x1 = ψ(x2 , · · · , xn ) = ψ(x∗ ) =
X
k∈Zn−1
c[k] β d (
x∗
− k),
h
(1)
where β d (·) is the uniform symmetric (n − 1)-dimensional B-spline of degree
d. The knots of the B-splines are located on a rectangular grid defined on the
chosen coordinate system, with a regular spacing given by h. The coefficients
of the B-spline representation are gathered in c[k].
2.3.2. Endocardial and Epicardial Contour Coupling
In typical active contour formulations, the endo and epicardial boundaries are defined separately and additional penalties guarantee their coupling
(Paragios, 2002; Chen et al., 2008). Inspired by the strategy used in the early
work of Barbosa et al. (2010) for the LV segmentation in US, we model the
9
endo and epicardial contours as a combination of two explicit functions representing the myocardial wall position (ψWP (θ)) and wall thickness (ψWT (θ)),
respectively. Thereby, each contour is defined as:
Γendo (θ) = ψWP (θ) − ψWT (θ),
(2)
Γepi (θ) = ψWP (θ) + ψWT (θ),
(3)
where ψWP (θ) gives the position of the myocardial wall as a function of θ and
ψWT (θ) encodes half of the myocardial wall thickness. Both Γendo and Γepi act
as interfaces between the blood pool (in), the myocardium (myo) and the
surrounding structures (out). The mathematical derivation that supports
the variational formulation of this coupling strategy, not included in the
early work of Barbosa et al. (2010), will be presented in Section 2.3.4 and
Appendix A. Moreover, the corresponding extension to a 3D formulation
will be proposed in Section 2.5.1.
The fundamental advantage of such formulation is that, by applying different scales (h in equation (1)) to ψWP and ψWT , one can implicitly smooth
the local variation of the myocardial thickness, in order to keep it smooth,
without affecting the overall ability to capture more complex myocardial
shapes. Moreover, it allows global shape constraints to be applied as shape
penalties in the energy functional and controlling the coupled segmentation
instead of each separated contour (Barbosa et al., 2010).
2.3.3. Proposed 2D Energy Functional
As commonly used in active contours methodologies (Lankton and Tannenbaum, 2008), our proposed energy functional is expressed as the sum of
a data attachment term, Ed , and a regularization term, Er , therefore being
expressed as:
E = Ed + Er
(4)
Since our model is defined based on the myocardial wall position and
thickness and given its B-spline formulation, it can be shown that the energy
functional is directly minimized through the optimization of the B-spline
coefficients c[k] (Barbosa et al., 2012). Therefore, the energy minimization
process can be achieved using the following evolution equations:
∂E
∂ Ed
∂ Er
=
+
∂ cWP [kj ]
∂ cWP [kj ] ∂ cWP [kj ]
10
(5)
∂E
∂ Ed
∂ Er
=
+
∂ cWT [kj ]
∂ cWT [kj ] ∂ cWT [kj ]
(6)
2.3.4. 2D Data Attachment term
Lankton and Tannenbaum (2008) have recently introduced the concept
of localized region-based energies, which overcome the limitations of classical edge-based and global region-based energy functionals. In the present
work, we propose to use two different localized region-based energies for the
endo and epicardial contours, which take the specific local appearance of
each boundary into account. Since the blood pool is brighter than the myocardial tissue, a signed version of the localized Yezzi (LY) energy (Lankton
and Tannenbaum, 2008; Barbosa et al., 2013a) is used for the endocardial
contour (henceforward referred as signed localized Yezzi - SLY). This energy
explicitly introduces the expected intensity relation between the in and myo
regions, forcing the contour to evolve towards an interface whose internal region is brighter than the external one. On the other hand, since the external
region (out) is relatively unpredictable due to the encountered variability in
structures (e.g., lungs, RV or soft tissues - Figure 2A), an unsigned constant
intensity model is chosen by using the localized Chan-Vese (LCV) energy
(Lankton and Tannenbaum, 2008).
Additionally, these localized energy functionals were modified by adding
weights to the inner and outer regions of the endocardial and epicardial interfaces, win and wout , respectively. This formulation allows an explicit control
over the equilibrium point in between the two regions, therefore pushing
the contour towards the myocardial region if a weight smaller than one is
considered (see Figure 3A). In practice, this strategy mimics the physician’s
knowledge when manually contouring, where they place the contours closer
to the myocardium to avoid TPMs to be included in this region, as illustrated
in Figure 3B-C.
Given the topology of the segmentation problem, the proposed energy
functional is defined as:
Z
Z
Ed =
δφendo (x) B(x, y) · Fendo (y) dy dx
(7)
Ω Z
Ω Z
+ δφepi (x) B(x, y) · Fepi (y) dy dx,
Ω
Ω
where
Fendo (y) = ux,myo − win · ux,in ,
11
(8)
Figure 3: a) Schematic representation of the weights influence in the proposed weighted
energy functionals. b) Example mid-ventricular CMR image. c) Two profiles of image B
with indication of reference contours points (indicated in red and green) and local means
of in, myo and out regions (indicated in black).
Fepi (y) = fmyo (y) · (Hφepi (y) − Hφendo (y))
+ fout (y) · (1 − Hφepi (y)),
(9)
fmyo (y) = (I(y) − ux,myo )2 ,
(10)
fout (y) = wout · (I(y) − ux,out )2 ,
(11)
and wi is the scalar weight applied to the region i; ux,in , ux,myo and ux,out
represent the local intensity means in the vicinity of x for the inner, myocardial and outer regions, respectively, as defined in Lankton and Tannenbaum
(2008); and I(y) is the image intensity at position y. Hφj is the Heaviside
operator applied to the level-set like function φj (x) = Γj (x∗ ) − x1 , representing the region inside the interface j. B(x, y) corresponds to a mask function
in which the local parameters that drive the evolution of the interfaces are estimated. In the present work, in order to maintain low computational costs,
we chose to use a squared region centered in x, whose size will be discussed
in Section 4.
Given the derivation available in the Appendix A, the minimization
of the data attachment term can be achieved using the following evolution
equations:
Z
∂ Ed
x∗
(12)
=
ḡendo (x∗ )β d ( − kj ) dx∗
∂ cWP [kj ]
h
Γendo
Z
x∗
ḡepi (x∗ )β d ( − kj ) dx∗
+
h
Γepi
12
∂ Ed
=−
∂ cWT [kj ]
x∗
− kj ) dx∗
h
Γendo
Z
x∗
ḡepi (x∗ )β d ( − kj ) dx∗
+
h
Γepi
Z
ḡendo (x∗ )β d (
(13)
where
¯ ∗
I(x ) − wmyo · ux,myo
ḡendo (x ) =
Aumyo
¯ ∗
I(x ) − win · ux,in
,
−
Auin
∗
(14)
ḡepi (x∗ ) = f¯myo (x∗ ) − f¯out (x∗ ).
(15)
and Aumyo and Auin are the areas of the myocardial and inner regions used
to estimate the local means ux,myo and ux,in , respectively (Lankton and Tannenbaum, 2008). Moreover, considering a generic function h(x) in Rn , we
note h̄ as the restriction of h over the interface Γ in R(n−1) (c.f. equation
¯ ∗ ) corresponds to the image value at position x = {x1 =
(A.10)). Thus, I(x
ψ(x∗ ), x2 , · · · , xn }.
2.3.5. 2D Regularization term
Since it is possible that both contours are attracted towards the same
image feature leading to their merging, a thickness penalty is included in the
evolution equations of the coupled segmentation problem. Such regularization term prevents the contours’ merging by locally penalizing a myocardial
thickness of less than 1 px. Let us note that in Chen et al. (2008), the
authors also introduced a thickness constraint in the energy, but they considered it as a global term (a deviation from the average thickness) modeled
as two distinct level-sets and not considering its local variation. Since the
term introduced in Chen et al. (2008) favors constant thickness, which limits
the flexibility of the model to deal with pathological cases, the formulation
introduced in Dietenbeck et al. (2011) was followed instead. However, the
term introduced in Dietenbeck et al. (2011) has an “on-off” behavior below
the minimum thickness, whereas a more gradual penalization of small thicknesses would be desirable. To this end, we have modified the term introduced
in Dietenbeck et al. (2011), thus expressing the thickness regularization term
as:
Z
1
ET =
φ(x + RT N)2 H(φ(x + RT N))δ(φ(x))dx,
(16)
2 Ω
13
where N corresponds to the inward normal of a point x belonging to the contour implicitly defined by the level-set function φ(x) and RT is the thickness
value below which the regularization term is activated. Following the deriva∂
tion presented in Dietenbeck et al. (2011), and given that ∂φ
φ(y)2 = 2φ(y),
the variation of this energy term wrt. φ can be expressed as:
∂ET
= φ(x + RT N)H(φ(x + RT N))δ(φ(x)).
∂φ
(17)
Thanks to the coupled formulation here presented, the term φ(x + RT N)
can be simplified since the local (half) thickness of the myocardium is explicitly defined by ψWT (x∗ ). Thus, the term φ(x+RT N) can be simply expressed
as (RT −2∗ψWT (x∗ )). Note that as φ(x+RT N) > 0 only if the local thickness
is smaller than RT , also the term (RT − 2 ∗ ψWT (x∗ )) is equally positive only
when half of the thickness, ψWT (x∗ ), is smaller than RT /2. Taking into account the details presented in the Appendix, namely (A.7), (A.8) and (A.10),
the thickness regularization term can be used to directly update the B-spline
coefficients controlling the myocardial thickness using:
Z
∂ ET
= (RT − 2 ∗ ψWT (x∗ ))
(18)
∂ cWT [ki ]
∗
∗
d x
H(RT − 2 ∗ ψWT (x )) β ( − ki ) dx∗ .
h
2.4. Automatic Stack Initialization
The previous sections focused on obtaining an accurate and automatic
segmentation for a mid-ventricular slice. However, the goal is to segment
all SAX slices from base to apex, which demands an initialization on these
slices. In this section, an automatic 3D stack initialization is proposed.
The simplest initialization would consist of naively using the segmentation result at the mid-ventricular level as initial contours of other slices in
the stack, plus a possible scaling to account for different sizes of the LV
SAX cuts from base to apex. However, this approach makes two unreliable
assumptions: it first considers that the relative LV size in different slices is
known or can be estimated a priori ; moreover, it assumes the cavity center is
always placed in the same position, which may not be true due to the known
misalignment of the slices in CMR datasets.
Thus, we propose a threshold-based BEAS formulation to initialize both
endo and epicardial contours in the full stack, while at the same time estimating possible slice misalignment.
14
The idea is to first estimate the expected intensity distribution of the
LV cavity and myocardial wall based on the previous 2D segmentation. We
propose to use an Expectation-Maximization (EM) algorithm (Moon, 1996)
assuming a mixture of two gaussians to describe the LV (cavity plus myocardium) histogram. After estimating the intensity distribution parameters
(mean and standard deviation), two thresholds are defined, namely the minimum cavity intensity and the minimum myocardium intensity (Figure 4C).
Due to the presence of TPMs and because they should be included in the
blood pool, the cavity’s distribution is discarded and only the myocardial one
is used to define both thresholds. The cavity intensity threshold is defined
as the 95 percentile of the myocardial distribution, which is a sufficiently
high value to capture the full blood pool plus TPMs, while simultaneously
avoiding the contour’s leakage to the myocardium. On the other hand, the
myocardium intensity threshold is defined as the 10 percentile of the myocardial distribution, allowing to capture the myocardial wall and avoid the
leakage to darker regions (such as the lungs). In Section 5, the choice of these
thresholds is addressed and the influence of their variation is investigated.
To obtain a coarse delineation of the endocardium, we start by placing an
initial circular contour with a small radius (2 px) centered in the LV cavity’s
center of the mid-ventricular slice and evolve it using the minimum cavity
threshold. In this sense, if the intensity of the contour point is higher than the
threshold, the contour is allowed to grow. In contrast, when it is lower than
the threshold, the contour is shrunk (Shi and Karl, 2008). Moreover, at each
iteration, an additional step is done focusing on re-initializing the center’s
position to be the contours’ center of gravity and consequently resample
the contour in the polar domain. The idea is to better estimate the cavity’s
center for each slice independently. Finally, the evolution stops when the area
enclosed by the contour does not significantly change for a few iterations.
We then initialize an epicardial contour with the thickness obtained in
the 2D mid-ventricular segmentation and evolve it using the minimum myocardium threshold. To improve the initialization robustness, the coupling
strategy described in Section 2.3.2 is used in this second stage. This allows
the epicardial contour delineation, while also enabling a better endocardial
definition by surpassing possible errors in the first stage. Again, the evolution stops when the area enclosed by both contours stays nearly constant for
a few iterations.
The proposed method is applied to each slice in the stack (Figure 4D-I),
progressively using the LV cavity’s center from previously initialized adjacent
15
Figure 4: Stack initialization procedure example. a) Mid-ventricular slice ROI overlay
with fully automatic 2D segmentation; b) LV mask obtained using epicardial contour; c)
LV histogram, corresponding intensity distributions’ parameters estimated based on EM
algorithm and threshold definition; d), f) and h) coarse delineation of the endocardium
using threshold-based BEAS; e) g) and i) final initialization for both endo and epicardial
contours using coupled threshold-based BEAS.
slices to start the endocardial contour. In the end, it provides a consistent
initialization for the subsequent segmentation method, while estimating the
slice misalignment.
2.5. Coupled Myocardial 3D Segmentation
2.5.1. BEAS 3D extension
The proposed automatic initialization procedure allows to obtain an initial contour plus an estimate for the LV cavity’s center for each slice. Using
such initialization, a segmentation could be independently applied in each
slice. However, it would not take into account the known spatial continuity
of the LV. Thus, a 3D segmentation enforcing a spatial continuity allows a
more robust segmentation, minimizing possible errors in individual slices.
In the BEAS framework, an important step towards this 3D modeling is
the choice of an appropriate coordinate system. As the CMR slices have a
relatively large thickness and the pixel intensities between slices cannot be
reliably estimated, considering a spherical coordinate system (which would
be the most straightforward BEAS extension) is unfeasible. Therefore, we
propose to use a cylindrical coordinate system, where the third dimension
has only a finite number of possible values corresponding to the number of
slices in the stack (Figure 5). At the same time, we improved the formulation
introduced in Barbosa et al. (2010) (Section 2.3.2) for an implicit coupled
16
propagation by extending it to a three-dimensional parameterization of the
myocardium, making use of the cylindrical topology of the CMR data. Thus,
the 3D extension of equations (2) and (3) yields:
Γendo (θ, z) = ψWP (θ, z) − ψWT (θ, z),
(19)
Γepi (θ, z) = ψWP (θ, z) + ψWT (θ, z),
(20)
where ψWP (θ, z) gives the position of the myocardial wall as a function of
θ and the slice level z, and ψWT (θ, z) encodes half of the myocardial wall
thickness. One should note that the proposed formulation only allows the
surface points to evolve radially, inside its own initial z level. Such constraint
is crucial as any inter-slice information is unknown.
Moreover, several points are important to ensure a proper surface evolution. On the one hand, extending the neighborhood to a cubic local volume
would lead to incorrect local statistics calculations due to the data anisotropy
and considerable slice thickness. As such, we continue to use an in-slice
squared region to update the local means around each interface. On the other
hand, the introduction of a new dimension requires an additional B-spline
scale to be defined in the formulation. However, it has to be a low scale to
avoid over-smoothing of the 3D surface in the third dimension, as a reduced
number of slices are acquired to cover the left ventricle.
Importantly, one should note the spatial continuity of the cylindrical parametric representation might be affected in case of slice misalignment due to
the smooth nature of our B-spline representation. Therefore, in order to
represent the LV surface with a consistent cylindrical topology and to allow
for spatial continuity, we propose to reduce the slice misalignment prior to
the 3D segmentation by considering the LV cavity’s centers estimated during
the automatic stack initialization (Section 2.4). To this end, a translation
is applied to each slice correcting the center’s position, resulting in a virtually straight LV axis, as illustrated in Figure 6. Such alignment allows
for a synergistic interaction between the segmented contours in neighboring
slices during surface evolution, therefore enabling a truly three-dimensional
analysis of the data without compromising the performance of the proposed
3D segmentation algorithm.
2.5.2. 3D Data Attachment term
Due to the contrast in MR images and subsequent well-defined edges,
several authors use edge-based energy terms to evolve the contours and segment the LV (Paragios, 2002; Pluempitiwiriyawej et al., 2005; Chen et al.,
17
Figure 5: Endocardial LV 3D-surface representation through an explicit function in the
cylindrical domain. a) Stack of 2D slices. b) LV 3D-surface in cartesian coordinates. c)
Correspondence between the cartesian and cylindrical domains; d) Explicit function in the
cylindrical domain.
2008). Although such edges are well defined at the endocardium, they usually
lead to wrong contouring when TPMs are present unless another strategy is
added. Regarding the epicardial delineation, even with low contrast between
myocardium and surrounding structures (e.g., lungs), edge-based terms can
substantially help if properly included. Taking this into account, we propose
to use a hybrid formulation for the epicardial energy functional by introducing an edge-based term. This term is only applied to the epicardium, because
the presence of TPMs in the blood pool may mislead the contours’ evolution
due to additional edges. Furthermore, as a result of higher contrast between
cavity and muscle, the segmentation result in the endocardium is expected
to be better than the epicardial one. The proposed edge-based term relies
on the computation of an edge map and its subsequent post-processing. Inspired by the method in Lee et al. (2010), for each slice in the ED phase an
edge map is computed following 4 steps:
1. The image is mapped to polar space using the LV cavity’s center estimated during the automatic stack initialization (Figure 7B);
2. A Canny edge detector (Canny, 1986) is applied to extract the edge
information (Figure 7C). Note that the high threshold was defined in
order to obtain 60% of the edge information, while the lower threshold
was 0.4 times the high one;
3. A post-processing step is applied to filter the edge information (Figure
7D). Specifically, the initial endocardial contour is used to filter out any
edge belonging to the cavity’s boundary or inner region. Moreover, the
edges in more remote regions are also cleaned as they are associated
with a myocardium thickness that is unlikely (over 20 mm);
18
Figure 6: Misalignment correction based on translation applied according to the estimated
LV cavity’s centers of each slice.
4. The edge energy map is computed by calculating the unsigned distance of every image’s pixel to the closest edge point (Figure 7E). This
distance function allows to improve the edge detection robustness by
filling small gaps in missed edges.
Using the resulting edge energy map, Υ(y), as a new term in the epicardial energy, the complete proposed energy functional for 3D segmentation is
defined as:
Z
Z
Ed =
δφendo (x) B(x, y) · Fendo (y) dy dx
(21)
Ω Z
Ω Z
+ δφepi (x) B(x, y) · (Fepi (y) + α · Υ(y)) dy dx,
Ω
Ω
where α is a positive hyper-parameter to balance the two terms in the epicardial energy.
19
Figure 7: a) Mid-ventricular slice ROI; b) Input image in polar domain; c) Edge information by Canny edge detector; d) Post-processed edge map; e) edge energy map overlay
with radial vector field.
The energy minimization process can be achieved using the following
evolution equations:
Z
∂ Ed
x∗
(22)
=
ḡendo (x∗ ) β d ( − ki ) dx∗
∂ cWP [ki ]
h
Γendo
Z
x∗
(ḡepi (x∗ ) + α · ∇r Ῡ(x∗ )) β d ( − ki ) dx∗ ,
+
h
Γepi
x∗
− ki ) dx∗
h
Γendo
x∗
(ḡepi (x∗ ) + α · ∇r Ῡ(x∗ )) β d ( − ki ) dx∗ ,
h
Γepi
∂ Ed
=−
∂ cWT [ki ]
Z
+
Z
ḡendo (x∗ ) β d (
(23)
where ∇r Ῡ is the edge energy map gradient along the radial direction.
2.5.3. 3D Regularization term
Similarly to the 2D case, the 3D regularization term comprehends a thickness term for contours’ merging prevention.
Moreover, even though the B-spline support of the coupled contours ensures a smooth interface, strong image features near the papillary muscles
can still wrongly attract the endocardial contour. To avoid this, a curvaturebased shape prior is also included as a regularization term in the endocardial
energy functional to avoid concave shapes. This can be easily achieved using a binary curvature-based switch that locally controls if the endocardial
contour should evolve according to the image features or if it should rather
20
be regularized to avoid concave shapes. The proposed shape regularization
strategy was simply implemented making use of a conventional mean curvature flow term, as the one used in Chan and Vese (2001), which is only applied
to regions of negative curvature. Thus, the endocardial feature function was
modified by:
ĝendo (x∗ ) = H(κendo (x∗ ))ḡendo (x∗ )
− H(−κendo (x∗ ))κendo (x∗ ),
(24)
where κendo (x∗ ) is the local curvature on the endocardial contour. Note that
the binary weighting term H(−κendo (x∗ )) constrains the curvature regularization to concave regions in the segmented endocardial interface, whereas
H(κendo (x∗ )) allows to normally evolve the endocardial contour according to
the segmentation functional introduced in (21). The relevance of this regularization term in the proposed 3D segmentation energy is discussed in Section
6.
2.6. 3D+time Tracking using Anatomical Affine Optical Flow
In the previous sections, a 3D segmentation for all slices covering the LV
was obtained for the ED phase. In order to compute the left ventricular EF,
one requires the segmentation to be performed in both ED and ES phases.
Several approaches can be used for this purpose, namely temporal segmentation or tracking procedures. The pure segmentation-based approaches focus
on identifying the relevant object on consecutive images, using specific energy
terms for temporal integration (Lynch et al., 2008). In contrast, trackingbased approaches normally focus on following the motion of specific patterns
along the sequence of images (Hautvast et al., 2006).
We propose to use a tracking algorithm to segment the LV in the whole
cardiac cycle, by estimating the motion between adjacent cardiac phases
(starting from the previously segmented contours in the ED phase). Because
CMR data is highly anisotropic and the motion cannot be reliably estimated
in the third dimension, each slice sequence should be independently analyzed.
Moreover, as both endo and epicardial contours need to be tracked through
all phases, the motion estimation should be performed independently for each
boundary as well.
In this sense, the motion estimation is performed using the global anatomical affine optical flow method proposed by Barbosa et al. (2013b), firstly
presented for three-dimensional US data. The principle is to estimate an
21
anatomically constrained global affine transform between two consecutive
frames by taking into account the estimated motion over the reference contour points in the previous frame. Such constrained estimation only considers
anatomical regions of interest and avoids the influence of adjacent tissues.
Moreover, by including an iterative displacement refinement scheme, the algorithm is still able to accurately capture large displacements (Barbosa et al.,
2013b).
Let I(x1 , x2 , t) denote the pixel intensity at location x = [x1 , x2 ] and time
t for a given slice temporal sequence. Based on the least square solution
of the optical flow equation, the 2D affine motion (for each slice and each
boundary) can be estimated by minimizing the following energy:
Z
W(x1 − c1 , x2 − c2 )(Ix1 u + Ix2 v + It )2 dx,
(25)
EM =
R2
where
u(x1 , x2 ) = u0 + ux1 (x1 − c1 ) + ux2 (x2 − c2 ),
(26)
v(x1 , x2 ) = v0 + vx1 (x1 − c1 ) + vx2 (x2 − c2 ),
(27)
encode the local motion field along x1 and x2 respectively. From such estimated motion field, the affine transformation to be applied to the analyzed
slice and contour can be written as the augmented matrix:
1 + ux 1
ux 2
u0
1 + v x2 v 0 ,
(28)
A = v x1
0
0
1
∇I = [Ix1 , Ix2 ] is the local image spatial gradient, whereas It corresponds to
the temporal derivative. In order to reduce the effect of noise in the assessment of the spatial gradient, ∇I(x) was computed using gaussian derivative
kernels along x1 and x2 , with σ = 1. It was simply estimated with finite
forward differences. W is a local window function centered in the position
c = [c1 , c2 ] and, according to the anatomical constraint introduced in Barbosa
et al. (2013b), is given by:
X
W(x) =
N (xk ),
(29)
xk ∈Γ
22
where N (x) is a neighborhood function defined as a 2D square center in x
and xk ∈Γ stands for the sampled points in the LV contour being tracked. N
was defined as a 7 × 7 square centered in the target point. Such region is
smaller than the one in US data (Barbosa et al., 2013b) due to the larger pixel
spacing in CMR images and reduced uncertainty associated with this data.
To decrease the method’s computational burden, we propose to subsample
the contours by a factor of 4, reducing the required local motion estimations.
Such subsampling is possible because superposition of local regions can be
avoided without losing robustness for the global estimation.
The proposed pure tracking-based energy is used to find the global affine
transformation between consecutive frames, independently applied to each
SAX cut and each boundary. Using such principle, the propagation is applied sequentially between adjacent phases using the previous LV contours,
until every cardiac phase has been tracked. Note that to minimize error accumulation, we ensure a minimal number of propagation steps by propagating
the contours in both directions for half the total number of cardiac phases.
3. Implementation Details
Concerning BEAS-related parameters in the 2D segmentation, the Bspline scale, h, was set to 22 and 23 for the wall position and thickness respectively, while for the 3D segmentation a value of 23 and 24 was used. The
increased spacing allows a higher degree of smoothness of the surface’s shape
and thickness, which helps overcoming possible errors due to the presence of
TPMs or in cases of low contrast and reduced spatial information. Regarding the smoothing across slices, a value of 20 was set. Similar to Queirós
et al. (2013), 64 points per slice were used on each boundary for the contour
discrete representation. Furthermore, a local region of size 15x15 and 11x11
was used to assess the local region means for the endo and epicardial contours, respectively. The local region is smaller for the epicardium because of
the unpredictability of the outer region, which may introduce errors in the
local mean estimation. The proposed region weights, wi , were set to 0.3 for
both in and out regions. As stated in Section 2.3.3, a value lower than one
pushes the contours towards the myocardial region, surpassing the presence
of TPMs. The chosen values allow a good compromise for contour delineation, without leakage into the myocardial region. At last, the parameter
α in the hybrid epicardial energy functional for 3D segmentation was empirically set to 0.05. For further general considerations and implementation
23
details regarding BEAS, the reader is referred to the work of Barbosa et al.
(2012). In Section 5, the choice of these parameters is addressed and the
influence of their variation in the results is investigated.
4. Validation experiments
To assess the performance of the proposed framework, we used the cine
steady state free precession (SSFP) MR SAX datasets available from the
MICCAI 2009 Cardiac MR Left Ventricle Segmentation Challenge (Radau
et al., 2009). The database includes 45 CMR datasets obtained with a 1.5T
GE Signa MRI, during 10- to 15-s breath holds with a temporal resolution
of 20 phases over the cardiac cycle. Six to twelve images were obtained
from the atrioventricular ring to the apex (thickness=8mm, FOV=320mm ×
320mm, matrix=256x256), with an isotropic pixel spacing ranging from 1.29
to 1.56mm. The data includes 12 patients with heart failure and myocardial
infarction (HF-I), 12 with heart failure but no myocardial infarction (HFNI), 12 with hypertrophic cardiomyopathy (HYP) and 9 normal cases (N).
Moreover, a ground truth, drawn by an experienced cardiologist, is available
in all slices for the epicardium and endocardium at the ED phase and for the
endocardium at the ES phase. All reported results consider manual contours
drawn to include TPMs in the LV cavity.
The segmentation performance was assessed against reference contours
using the evaluation code provided in Radau et al. (2009), which evaluates
the Dice coefficient, the average perpendicular distance (APD) and the percentage of “good contours” (contours for which the APD is less than 5 mm).
Moreover, the Hausdorff distance (defined as the maximum perpendicular
distance between contours) was also computed to understand the maximum
local error of the proposed framework. Each measure was computed slice by
slice and a mean value for all slices of a dataset was calculated. However,
the original evaluation code only considers “good contours” when computing
these means. Besides this limitation, it also removes the slices only present in
one of the cardiac phases (ED or ES), disregarding longitudinal contraction.
Taking into account the interest in future comparisons with methods that also
use this database, the results obtained with the original evaluation code are
reported. Nevertheless, to overcome both limitations and fully understand
the framework performance, the segmentation metrics were also computed
considering all the available contours using a customized evaluation code.
24
The above measures were used to validate both 3D and 3D+time segmentation results. Moreover, our automatic 3D+time framework was compared
to the current state-of-the-art methods that used the same database, using
the original evaluation code for metrics computation.
In addition to quantitative technical measures, 5 clinical parameters were
computed, namely EDV, ESV, SV, LVM and EF, to validate the performance of the proposed 3D+time framework in a clinical context (Petitjean and Dacher, 2011; Hu et al., 2012). Therefore, linear regression and
Bland-Altman analysis were performed to assess the correlation between reference clinical indices and the automatically computed ones. Again, the evaluation code presents the aforementioned limitations preventing the correct
assessment, but such results are still reported to allow future comparisons for
both EF and LVM, as stated in Radau et al. (2009). Nevertheless, the customized code is also used to report clinically correct values using all contours
initially considered by the expert during manual segmentation.
Finally, the influence of the parameters variation was assessed to study
the robustness and stability of the proposed framework. To this end, the
8 most relevant parameters were included in the analysis, namely the inner
and outer region weights (win and wout ), the sizes of the local squared regions
([n × n]endo and [n × n]epi ), the threshold values for automatic stack initialization (Ti,endo and Ti,epi ), the hyperparameter α and the number of points
for contour discrete representation (npoints ). The range where the parameter
influence is evaluated was defined as a ±50% variation around the aforementioned chosen parameters values. Each dataset was then segmented using
the new value of one parameter while keeping the remaining ones fixed to
their original value. Please note that due to the multiscale operation involved
in the BEAS framework, npoints should be a base-2 number. Therefore, its
range was defined as a -50% to +100% variation around the initial predefined
value (64 points). Moreover, in the case of the hyperparameter α, the zero
value was also considered to study the importance of the edge-based term
introduced in the 3D epicardial energy functional.
The reported mean computational times were obtained using MATLAB
code running on a Intel(R) Core(TM) i5 CPU at 2.4 GHz (shared memory
of 4 GB).
5. Results
Table 1 summarizes the performance for the 3D segmentation, computed
25
Table 1: 3D segmentation performance for endocardial and epicardial contours in ED
phase (# = 45, µ ± σ).
Dice
Endo
Epi
APD (mm)
Endo
Epi
Hausdorff (mm)
Endo
Epi
Good Contours (%)
Endo
Epi
OC 0.93 ± 0.03 0.94 ± 0.02 1.50 ± 0.47 1.80 ± 0.41 3.97 ± 1.24 4.65 ± 1.14
97.4 ± 8.2 95.4 ± 9.6
CC 0.92 ± 0.04 0.93 ± 0.03 1.70 ± 0.89 2.06 ± 0.84 4.34 ± 1.91 5.14 ± 1.79
Original code (OC) only considers “good contours” for metrics computation, while the customized
code (CC) uses all segmented slices.
Figure 8: Boxplots of computed a) Dice, b) APD and c) Hausdorff distance for endocardial
and epicardial contours using the proposed automatic 3D segmentation (# = 45). The
ends of the whiskers represent the lowest and highest data point still within 1.5 times
the inter-quartile range of the lower and upper quartile, respectively. Any data value
larger is considered an outlier and plotted as a dot. The results were computed using the
customized code.
using both original evaluation code (OC) and customized code (CC). To also
assess the overall performance distribution rather than only average results,
Figure 8 presents the boxplots for each metric for both endocardial and epicardial contours using the customized code. Moreover, the distribution from
basal to apical slices of both APD and Hausdorff distance metrics are illustrated in Figure 9. Note that the Dice metric distribution is not presented,
because it is influenced by the relative size of the object (which is smaller in
apical slices) and thus not suitable for comparison between slices. Figure 10
gives three segmentation examples with several slices from base to apex.
Regarding 3D+time LV segmentation, Table 2 presents the overall segmentation plus tracking performance (using CC) for all 45 datasets and also
categorized by pathology. Figure 11 gives an example tracking result over
the cardiac cycle.
The comparison of the proposed approach against state-of-the-art methods that use the same database is presented in Table 3, using the original
26
Figure 9: Boxplots for APD and Hausdorff metrics distribution from basal to apical slices
computed for the proposed automatic 3D segmentation results (# = 45). a) APD for
endocardial contour; b) APD for epicardial contour; c) Hausdorff distance for endocardial
contour; d) and Hausdorff distance for epicardial contour. The ends of the whiskers
represent the lowest and highest data point still within 1.5 times the inter-quartile range
of the lower and upper quartile, respectively. The crosses represent the mean values for
each metric. The results were computed using the customized code.
Table 2: 3D+time segmentation performance for endocardial and epicardial contours,
categorized by pathology and average result (# - number of datasets, µ ± σ).
#
HF-I
HF-NI
HYP
N
Total
12
12
12
9
45
Dice
Endo
a b
0.92±0.05 ,
c d
0.91±0.03 ,
a c
0.81±0.06 ,
b d
0.85±0.05 ,
0.87±0.07
Epi
0.94±0.03
0.92±0.04
0.91±0.03
0.92±0.03
0.93±0.03
APD (mm)
Endo
Epi
1.93±1.28
2.01±0.57
2.40±0.61
2.25±0.79
2.14±0.88
1.93±0.96
2.17±0.86
2.14±0.63
2.01±0.83
2.06±0.84
Hausdorff (mm)
Endo
Epi
4.66±2.38
5.42±2.01
4.95±0.93
5.22±1.58
5.05±1.85
5.09±1.91
5.37±2.09
5.18±1.36
4.87±1.65
5.14±1.79
Good Contours (%)
Endo
Epi
93.3±12.0 95.1±12.6
94.2±6.3 93.8±9.4
90.1±9.6 98.3±3.7
93.2±8.6 94.0±9.9
92.7±9.5 95.4±9.6
ap
< 0.05, unpaired t-test between HF-I and HYP; b p < 0.05, unpaired t-test between HF-I and N;
< 0.05, unpaired t-test between HF-NI and HYP; d p < 0.05, unpaired t-test between HF-NI and N.
Metrics computed using our customized code, which considers all available contours.
cp
evaluation code (Radau et al., 2009) for metrics computation. Moreover, the
average total computational time to segment a 4D dataset is also presented.
Regarding the clinical validation, the correlation coefficient obtained by
regression analysis and both bias and standard deviation obtained by BlandAltman analysis are presented in Table 4 for both versions of the evaluation
code. Figure 12 and Figure 13 present the linear regression and BlandAltman plots, respectively, obtained for EDV, ESV, LVM and EF (using
CC).
To better understand the computational burden of each module of the
proposed framework, the average computational time for each step is presented in Table 5 along with the average total computational time to segment
a dataset.
Finally, Figure 14 illustrates the influence of the variation of each parameter in the automatic 3D+time segmentation performance. Furthermore, in
27
Figure 10: Automatic 3D segmentation result for three example CMR datasets (red:
endocardium; green: epicardium; yellow: ground truth).
order to stress the importance of the hybrid energy over a pure region-based
one, the anatomical distribution of the epicardial APD errors with and without the proposed edge-based term is presented in Figure 15.
6. Discussion
The first module of the proposed framework was evaluated in the preliminary contribution (Queirós et al., 2013), in which the 2D automatic initialization algorithm showed a high feasibility (91.1%) and robustness. However, four wrongly initialized images were observed in the 45 datasets used
(Queirós et al., 2013). Nonetheless, we still included these datasets in the
present 3D and 3D+time analysis, as it was found that the proposed stack
initialization is able to recover from these bad initializations. Regarding the
28
Figure 11: Automatic 3D+time segmentation result for an example CMR dataset (magenta: automatic contour at ED phase; red: endocardium; green: epicardium; yellow:
ground truth).
second module, the proposed energy for 2D segmentation combined with the
automatic 2D initialization resulted in an accurate endo and epicardial delineations (data not shown). Hence, the first block of the proposed framework
showed to be efficient and robust to obtain a first mid-ventricular delineation
for further initialization.
Concerning the proposed automatic stack initialization, the ability to include the 4 wrongly initialized images may be linked to the robustness of the
EM algorithm and the simple but powerful principle of thresholding. In other
words, even if the initial 2D mid-ventricular segmentation has a sub-optimal
result, the applied mask and subsequent EM algorithm still allows a good
estimation of the intensity distribution of the LV cavity and myocardial wall.
After properly estimating the two gaussian distributions, both myocardial
and cavity thresholds are also appropriately defined. Consequently, the initial contours obtained for all slices from base to apex are fairly close to the
true boundaries, which permits the proposed 3D localized energy functional
to obtain an accurate result during segmentation. In fact, 2 of these 4 cases
29
Table 3: Comparison of segmentation performance between proposed approach and stateof-the-art methods (# - number of datasets evaluated, µ ± σ).
Dice
Endo
Epi
APD (mm)
Endo
Epi
Good Contours (%)
Endo
Epi
Method
#
Time (s)
Proposed
45
11.13±2.76
30
-
0.89±0.04 0.92±0.02 2.04±0.47 2.35±0.57
90.35
92.56
45
∼60
0.86±0.05 0.91±0.03 2.44±0.56 2.80±0.71
80±16
71±26
Constantinides
et al. (2009)
Constantinides
et al. (2012)
Hu et al.
(2012)
Huang et al.
(2011)
Jolly (2009)
Liu et al.
(2012)
Schaerer et al.
(2010)
Uzunbas et al.
(2012)
Wijnhout
et al. (2009)
0.90±0.05 0.94±0.02 1.76±0.45 1.80±0.41 92.7±9.5
95.4±9.6
45 153.86±32.13 0.89±0.03 0.94±0.02 2.24±0.40 2.19±0.49 91.06±9.4 91.21±8.5
45
-
0.89±0.04 0.93±0.02 2.16±0.46 2.22±0.43 79.2±19.0 83.9±16.8
30
∼60
0.88±0.04 0.93±0.02 2.26±0.59 1.97±0.48 95.62±8.8 97.29±5.8
45 129.45±30.39 0.88±0.03 0.94±0.02 2.36±0.39 2.19±0.49 91.17±8.5 90.78±10.7
45
-
0.87±0.04 0.92±0.02 2.97±0.38 3.14±0.33
-
-
15
∼45
0.82±0.06 0.91±0.03 2.98±0.88 1.78±0.35
-
-
15
∼60
0.89±0.03 0.93±0.01 2.29±0.57 2.28±0.39 86.4±11.0
94.2±7.0
Metrics computed using original evaluation code, which considers only “good contours” for metrics
computation.
Table 4: Linear regression and Bland-Altman analysis for clinical cardiac indexes (# = 45).
R
EDV
Bias
σ
R
ESV
Bias
σ
R
SV
Bias
σ
R
LVM
Bias
σ
R
EF
Bias
σ
a
OC
0.906 -0.01 18.24 0.977 2.45 4.83
a
a
CC 0.985 -2.46 15.03 0.988 -3.83 14.78 0.955 1.38 8.07 0.951 5.51 14.13 0.976 3.00 5.60
a
p < 0.05, paired t-test against zero.
Original code (OC) only considers “good contours” for metrics computation, while the customized
code (CC) uses all segmented slices.
R is the Pearson product-moment correlation coefficient.
where the 2D initialization failed ended up having a 3D segmentation with errors below average. Moreover, the proposed method is at the same time able
to adequately estimate and reduce occasional misalignment between slices.
Such feature is possible due to the center estimation during contour growing, which adapts itself to dataset-specific acquisition displacement artifacts
(Figure 6). Although such strategy might not be anatomically correct, as
it could overcompensate for ventricles presenting curvilinear long-axis (e.g.,
banana-shaped heart) and thus deviate from the physiologic LV axis, it allows
for a 3D spatial continuous delineation when using the proposed cylindrical
30
Figure 12: Linear regression for a) end-diastolic volume (EDV); b) end-systolic volume
(ESV); c) LV mass (LVM); d) and ejection fraction (EF).
Figure 13: Bland-Altman analysis for a) end-diastolic volume (EDV); b) end-systolic
volume (ESV); c) LV mass (LVM); d) and ejection fraction (EF).
Table 5: Average computational time for
each step of the proposed approach (# =
45, µ ± σ).
Method
Time (s)
LV Location
Template Matching
2D Segmentation
Stack Initialization
3D Segmentation
Anatomical Affine Optical Flow
0.10±0.02
1.16±0.37
0.22±0.05
0.55±0.10
6.21±2.41
2.89±0.49
Total
11.13±2.76
topology.
Please note that the proposed framework assumes that the stack of slices
was previously cropped to only include slices covering the LV in the stack
initialization and subsequent segmentation. In case this was not done, a
simple solution would be to have the user input two points at run-time,
corresponding to the choice of basal and apical slices. Although the presented
framework does not require any user interaction during the segmentation
itself, such pre-processing of the data is still required. In order to obtain a
fully automatic methodology, we intend to address this issue in future work.
Regarding Table 1, the proposed 3D segmentation method and associ31
Figure 14: Influence of the variation of 8 parameters in the 3D+time framework performance. (A) and (B) Inner and outer region weights ([0.15; 0.45]); (C) and (D) sizes of the
local region ([9×9, 21×21]) and [5×5, 17×17] for endo and epicardial region, respectively);
(E) and (F) threshold for stack initialization ([92.5%; 97.5%] and [5%; 15%] for endo and
epicardial thresholds, respectively); (G) parameter α for hybrid epicardial energy (0 and
[0.025; 0.075]); and (H) number of points for contours discrete representation ([32; 128]),
(red: endocardium; green: epicardium). * p < 0.05, paired t-test against baseline settings.
32
Figure 15: Boxplots of epicardial APD distribution from basal to apical slices computed
using a pure region-based epicardial energy and a hybrid energy functional. The ends of
the whiskers represent the lowest and highest data point still within 1.5 times the interquartile range of the lower and upper quartile, respectively. The crosses represent the
mean values for each metric. The results were computed using the customized code.
ated hybrid energy functional proved to give accurate results in all metrics
used. A high percentage of “good contours” was obtained, with an average
of 97.4% and 95.4% for endo and epicardial contours, respectively. The slices
considered bad (AP D ≥ 5mm) were mostly located in the apex (due to small
size and low contrast with respect to the background) or in basal slices in the
presence of the left ventricular outflow tract (LVOT), as illustrated in Figure
9. The latter problem (i.e., myocardium appearing as an incomplete annular shape) was still not directly addressed in the proposed framework, being
only avoided in most cases due to the intra- and inter-slice spatial continuity
of the LV surface (see Figure 10A, at the base). Nevertheless, it increases
the measured distances and is the main contributor to the reported average
Hausdorff distance. Note the increased average and standard deviation in
this metric when considering all slices (CC compared to OC), although only
a few percentage of contours are added. As expected, the reported values for
our customized code are lower due to the inclusion of slices considered poorly
segmented. Such results reinforce the need to look together to both metrics
and the percentage of “good contours” when using the original evaluation
code.
Figure 8 gives an overview of the proposed methodology performance
in segmenting the ED phase, with more than 75% of the studied datasets
having a Dice metric higher than 0.90 for both endo and epicardial contours.
At the same time, these datasets also present an APD under 2 mm for both
contours and a maximum error around 4-5 mm. Taking into account the
33
observed outliers, it is important to note that two of these are from the 4
cases having a wrong 2D initialization. The other outliers are associated with
datasets with considerable heterogeneity in blood image intensity between
slices, hampering the correct automatic initialization of apical and basal slices
further away from the initial mid-ventricular slice when using our thresholdbased strategy. Nonetheless, despite its low-level intensity assumptions, the
proposed method proved to be globally generic.
Regarding the error distribution from basal to apical slices in Figure 9,
higher errors are observed for the two ends of the left ventricle, mainly due
to the LVOT in basal slices and to PVE in apical ones (as stated above).
Notwithstanding, overall the segmentation is robust independently from the
slice level. Figure 10 confirms the robustness of the proposed 3D segmentation against different myocardial intensities, contrast and presence of TPMs
(note the third column of Figure 10B).
The same robustness is observed in the 3D+time approach (Table 2).
However, a decrease in the endocardium accuracy for the 3 measured metrics
is detected compared to the above 3D results (Table 1). This is partially
related with the greater difficulty in distinguishing the compacted TPMs from
the myocardial tissue in ES slices, hampering the correct identification of the
endocardial boundary. This, in combination with the LVOT in basal slices
and PVE in apical ones, as already discussed for the ED frame, is the main
cause for contours considered badly segmented in the ES frame. Nevertheless,
the percentage of “good contours” is still high for the endocardium, with an
average 92.3%. Moreover, the APD is kept around 2 mm when considering
all contours in both ED and ES phases. Note that epicardial results are
unchanged as they are not evaluated in ES phase due to lack of manual
contouring. When dividing the studied datasets by pathology, no significant
difference is found between average APD between these groups. Regarding
the Dice metric, significant differences were found for the endocardial contour
between some groups (Table 2). However, this can be related to the different
cavity size between groups of subjects, since the differences in average APD
are not significant. In summary, these results suggest that the proposed
algorithm is robust against different pathologies and ventricular anatomy
and function. At the same time, it corroborates the feasibility of using an
affine transform to describe the in-plane LV motion during the 3D+time
tracking stage of the proposed algorithm.
According to Table 3, the proposed framework proved to be competitive,
as compared to the state-of-the-art frameworks that used this database for
34
validation, providing leading results for both accuracy and computational
times. This observation was valid for all segmentation performance measures employed, with a smaller average APD, superior average Dice metric
and one of the highest percentage of “good contours”. Note that the original
evaluation code only considers “good contours” for computation of means,
implying that our results are more accurate even considering a higher percentage of slices. On the other hand, it is important to note that the results
for some of the state-of-the-art methods were obtained during the MICCAI
challenge, therefore being tested in one set and validated in another. In
our case, we tested and validated on the full database. Nevertheless, the
presented sensitivity analysis (Figure 14) proved the robustness of the algorithm, discrediting the possibility of over-tuning of the different parameters
of the proposed algorithm. Finally, besides presenting high accuracy, the
framework showed also to be very fast (Table 5), presenting an average total computational time of 11.16 ± 2.76s per dataset (around 0.06 − 0.09s
per image) in a Matlab implementation. Since a fully manual segmentation
takes between 6 to 20min for only ED and ES phases (Hautvast et al., 2012;
Petitjean and Dacher, 2011), the proposed framework seems to be fast and
accurate enough to be of interest in daily clinical practice.
Linear regression revealed a high correlation for all clinical cardiac indices,
ultimately showing the accuracy and interest of the proposed framework for
automatic assessment of LV morphology and global function (Table 4 and
Figure 12). Regarding Bland-Altman limits of agreement (Figure 13), the
reported values are slightly higher than other frameworks in the literature
(Eslami et al., 2013; Cordero-Grande et al., 2011). Nonetheless, these methods require some user interaction or were only validated on patients with
a specific pathology. In the present work, an automatic framework is being proposed and validated on a database with both healthy subjects and
patients.
Notably, the sensitivity analysis of the proposed framework (Figure 14)
reveals that the overall performance is not significantly impaired within a
fairly large range of parameter settings. The most critical choice for the
overall segmentation performance was found to be the number of points used
to represent the myocardial contours (Figure 14H). However, using 32 points
instead of the proposed 64 points corresponds to a quite severe downsampling of the surface representation, which strongly constrains the degrees of
freedom of the segmented LV shape. On the other hand, using 128 points to
represent the contours did not lead to any significant accuracy improvement,
35
thus reflecting that the proposed 64 points provides an optimal balance between overall segmentation accuracy and robustness, as well as reduced computational burden. The influence of the endocardial scalar weight win was
found to also have a statistically significant influence on the segmentation
performance for variations of ±30% (Figure 14A). This is intrinsically linked
to the fact that this term controls the equilibrium point between the region
forces from the blood pool and myocardium, thus being translated as the
amount of trabeculae included in the myocardial region. On the other hand,
the influence of the epicardial scalar weight wout did not have such variability
pattern, statistically showing significant influence only on the Dice metric for
a variation of +50% (Figure 14B). All the other parameters either did not
show any statistically significant influence on the segmentation performance
or only showed small but significant effects for variations of ±50%.
Additionally, in the case of a pure localized region-based energy for the
epicardial contour in the 3D segmentation (α = 0), there was a small yet
statistically significant decrease on the global performance, mostly perceived
in the apical slices (see Figure 15). Such result is primarily linked to the increased capture range of the edge-based term, which helps handling the blurry
appearance of the myocardium in apical slices (caused by partial volume effects). A similar effect is achieved in cases of slightly worst initializations,
where the edge-based term counteracts the increased sensitivity to the initial
position of the region-based term. Overall, this result reinforces the interest
of the proposed hybrid energy functional in the 3D segmentation.
A note should be addressed to the proposed regularization terms. Whereas
the minimum thickness penalty was found to have an important contribution by preventing the local overlap of the endocardial and epicardial interfaces, the proposed curvature regularization approach, included to deal with
concave endocardial regions, was found to have a small segmentation error
reduction (< 5% in all metrics). However, the segmentation results became
visually more appealing, while it also positively contributed in more challenging slices. Thus, this term was kept in the 3D segmentation step of the
proposed framework. Note that the curvature-based prior was only included
in this step, since the 2D case is simply used as an initialization step towards
the subsequent 3D segmentation. In fact, and as stated above, the automatic
stack initialization step is capable of recovering in cases of sub-optimal 2D
segmentation results, thus avoiding the need for such extra 2D regularization
term.
Overall, the present analysis showed the robustness of the proposed frame36
work, with only slight variations in the segmentation performance for most of
the parameters analyzed. This result is of utmost importance to demonstrate
its feasibility against variability, highlighting the robustness of the method.
7. Conclusion
In summary, a novel automatic 3D+time segmentation framework for
CMR datasets was presented. The proposed approach was shown to be
efficient, accurate, robust and fast compared to the other state-of-the-art
methods, presenting leading results in both accuracy and computational time
in the common database used for comparative purposes. The computational
time to analyze an entire 4D dataset took on average 11s, with an overall
correlation with reference contours of 0.99, 0.99, 0.96, 0.95 and 0.98 for EDV,
ESV, SV, LVM and EF, respectively. Such results show the appeal of the
proposed framework for automatic assessment of LV morphology and global
function in daily clinical practice.
In the future, the proposed algorithm should be prospectively evaluated
in a larger database comprehending “real-world” variability in image quality
and anatomical and functional characteristics.
Appendix A. Coupled BEAS energy derivation
We consider here the differentiation of the coupled localized region-based
energy functional, E (7), with respect to a given B-spline coefficient c[kj ].
Given (7), one can rewrite it as:
Z
2 Z
2
X
X
E=
δφi (x) B(x, y) · Fi (y) dy dx =
Ei ,
(A.1)
i=1
Ω
Ω
i=1
with i = 1 and i = 2 representing the endocardium and epicardium, respectively,
Now, let consider the derivation of Ei only:
Z
Z
∂ Ei
∂ Fi
∂ φi
=
δφi (x) B(x, y)
·
dy dx
(A.2)
∂ cWP [kj ]
∂ φi ∂ cWP [kj ]
Ω
Ω
The term ∂∂ Fφii depends only on the feature function used and can be
expressed as (Barbosa et al., 2012, 2013a):
∂ Fi (y)
= gi (y) · δφi (y)
∂ φi
37
(A.3)
Noting moreover that B(x, y)δφ (x) is different from zero only if x = y,
we have B(x, y)δφ (x) = δφ (x − y) and therefore,
Z
∂ φi (x)
∂ Ei
=
gi (x)
δφ (x) dx
(A.4)
∂ cWP [kj ]
∂ cWP [kj ] i
Ω
Moreover, we have:
φendo (x) = Γendo (x∗ ) − x1 = ψWP (x∗ ) − ψWT (x∗ ) − x1
(A.5)
φepi (x) = Γepi (x∗ ) − x1 = ψWP (x∗ ) + ψWT (x∗ ) − x1 ,
(A.6)
and
which, from (1), leads to:
and
x∗
∂ φendo (x)
= β d ( − kj )
∂ cWP [kj ]
h
(A.7)
∗
∂ φepi (x)
d x
= β ( − kj )
∂ cWP [kj ]
h
(A.8)
Equation (A.4) thus becomes:
Z
∂ Ei
x∗
=
gi (x)β d ( − kj )δφi (x) dx
∂ cWP [kj ]
h
Ω
(A.9)
Taking into account the following property (Barbosa et al., 2012):
Z
Z
f (x)δφ (x) dx = f¯(x∗ ) dx∗ ,
(A.10)
Ω
we get:
∂ Ei
=
∂ cWP [kj ]
Γ
Z
ḡi (x∗ )β d (
Γi
x∗
− kj ) dx∗
h
(A.11)
and, thus, yielding equation (12). Using the same reasoning, we can get
equation (13).
38
References
Anderson, J., Weaver, A., Horne, B., Jones, H., Jelaco, G., Cha, J., Busto, H.,
Hall, J., Walker, K., Blatter, D., 2007. Normal cardiac magnetic resonance
measurements and interobserver discrepancies in volumes and mass using
the papillary muscle inclusion method. The Open General and Internal
Medicine Journal 1, 6–12.
Barbosa, D., Bernard, O., Savu, O., Dietenbeck, T., Heyde, B., Claus, P.,
Friboulet, D., D’hooge, J., 2010. Coupled B-spline active geometric functions for myocardial segmentation: A localized region-based approach, in:
IEEE International Ultrasonics Symposium (IUS), pp. 1648–1651.
Barbosa, D., Dietenbeck, T., Heyde, B., Houle, H., Friboulet, D., D’hooge,
J., Bernard, O., 2013a. Fast and fully automatic 3-D echocardiographic
segmentation using b-spline explicit active surfaces: Feasibility study and
validation in a clinical setting. Ultrasound in medicine & biology 39, 89–
101.
Barbosa, D., Dietenbeck, T., Schaerer, J., D’hooge, J., Friboulet, D.,
Bernard, O., 2012. B-spline explicit active surfaces: An efficient framework for real-time 3-D region-based segmentation. IEEE Transactions on
Image Processing 21, 241–251.
Barbosa, D., Heyde, B., Dietenbeck, T., Friboulet, D., D’hooge, J., Bernard,
O., 2013b. Fast left ventricle tracking in 3d echocardiographic data using
anatomical affine optical flow, in: Proceedings of Functional Imaging and
Modeling of the Heart (FIMH2013), pp. 191–199.
Ben Ayed, I., Chen, H., Punithakumar, K., Ross, I., Li, S., 2012. Max-flow
segmentation of the left ventricle by recovering subject-specific distributions via a bound of the bhattacharyya measure. Medical Image Analysis
16, 87–100.
Ben Ayed, I., Li, S., Ross, I., 2009. Embedding overlap priors in variational
left ventricle tracking. IEEE Transactions on Medical Imaging 28, 1902–
1913.
Canny, J., 1986. A computational approach to edge detection. IEEE Transactions on Pattern Analysis and Machine Intelligence , 679–698.
39
Chan, T., Vese, L., 2001. Active contours without edges. IEEE Transactions
on Image Processing 10, 266–277.
Chen, T., Babb, J., Kellman, P., Axel, L., Kim, D., 2008. Semiautomated segmentation of myocardial contours for fast strain analysis in cine
displacement-encoded MRI. IEEE Transactions on Medical Imaging 27,
1084–1094.
Ciofolo, C., Fradkin, M., Mory, B., Hautvast, G., Breeuwer, M., 2008. Automatic myocardium segmentation in late-enhancement MRI, in: IEEE
International Symposium on Biomedical Imaging (ISBI), pp. 225–228.
Constantinides, C., Chenoune, Y., Kachenoura, N., Roullot, E., Mousseaux,
E., Herment, A., Frouin, F., 2009. Semi-automated cardiac segmentation
on cine magnetic resonance images using GVF-Snake deformable models.
The MIDAS Journal - Cardiac MR Left Ventricle Segmentation Challenge
.
Constantinides, C., Roullot, E., Lefort, M., Frouin, F., 2012. Fully automated
segmentation of the left ventricle applied to cine MR images: Description
and results on a database of 45 subjects, in: IEEE Engineering in Medicine
and Biology Society (EMBC), IEEE. pp. 3207–3210.
Cordero-Grande, L., Vegas-Sánchez-Ferrero, G., Casaseca-de-la Higuera, P.,
Alberto San-Román-Calvar, J., Revilla-Orodea, A., Martı́n-Fernández, M.,
Alberola-López, C., 2011. Unsupervised 4d myocardium segmentation with
a markov random field based deformable model. Medical Image Analysis
15, 283–301.
Dietenbeck, T., Alessandrini, M., Barbosa, D., Dhooge, J., Friboulet,
D., Bernard, O., 2011. Detection of the whole myocardium in 2Dechocardiography for multiple orientations using a geometrically constrained level-set. Medical Image Analysis .
Eslami, A., Karamalis, A., Katouzian, A., Navab, N., 2013. Segmentation
by retrieval with guided random walks: Application to left ventricle segmentation in mri. Medical image analysis .
Hautvast, G., Lobregt, S., Breeuwer, M., Gerritsen, F., 2006. Automatic
contour propagation in cine cardiac magnetic resonance images. IEEE
Transactions on Medical Imaging 25, 1472–1482.
40
Hautvast, G., Salton, C., Chuang, M., Breeuwer, M., O’Donnell, C., Manning, W., 2012. Accurate computer-aided quantification of left ventricular
parameters: Experience in 1555 cardiac magnetic resonance studies from
the Framingham Heart Study. Magnetic Resonance in Medicine .
Hu, H., Liu, H., Gao, Z., Huang, L., 2012. Hybrid segmentation of left ventricle in cardiac MRI using gaussian-mixture model and region restricted
dynamic programming. Magnetic Resonance Imaging .
Huang, S., Liu, J., Lee, L., Venkatesh, S., Teo, L., Au, C., Nowinski, W.,
2011. An image-based comprehensive approach for automatic segmentation
of left ventricle from cardiac short axis cine MR images. Journal of Digital
Imaging 24, 598–608.
Jolly, M., 2009. Fully automatic left ventricle segmentation in cardiac cine
MR images using registration and minimum surfaces. The MIDAS Journal
- Cardiac MR Left Ventricle Segmentation Challenge .
Lankton, S., Tannenbaum, A., 2008. Localizing region-based active contours.
IEEE Transactions on Image Processing 17, 2029–2039.
Lee, H., Codella, N., Cham, M., Weinsaft, J., Wang, Y., 2010. Automatic left
ventricle segmentation using iterative thresholding and an active contour
model with adaptation on short-axis cardiac MRI. IEEE Transactions on
Biomedical Engineering 57, 905–913.
Liu, H., Hu, H., Xu, X., Song, E., 2012. Automatic left ventricle segmentation in cardiac MRI using topological stable-state thresholding and region
restricted dynamic programming. Academic Radiology .
Lu, Y., Radau, P., Connelly, K., Dick, A., Wright, G., 2009. Segmentation
of left ventricle in cardiac cine MRI: An automatic image-driven method.
Functional Imaging and Modeling of the Heart , 339–347.
Lynch, M., Ghita, O., Whelan, P., 2008. Segmentation of the left ventricle of
the heart in 3D+ t MRI data using an optimized nonrigid temporal model.
IEEE Transactions on Medical Imaging 27, 195–203.
Miller, C., Pearce, K., Jordan, P., Argyle, R., Clark, D., Stout, M., Ray, S.,
Schmitt, M., 2012a. Comparison of real-time three-dimensional echocardiography with cardiovascular magnetic resonance for left ventricular vol41
umetric assessment in unselected patients. European Heart Journal 13,
187–195.
Miller, C.A., Jordan, P., Borg, A., Argyle, R., Clark, D., Pearce, K., Schmitt,
M., 2012b. Quantification of left ventricular indices from SSFP cine imaging: Impact of real-world variability in analysis methodology and utility
of geometric modeling. Journal of Magnetic Resonance Imaging .
Moon, T., 1996. The expectation-maximization algorithm. Signal Processing
Magazine, IEEE 13, 47–60.
Otsu, N., 1975. A threshold selection method from gray-level histograms.
Automatica 11, 23–27.
Paragios, N., 2002. A variational approach for the segmentation of the left
ventricle in cardiac image analysis. International Journal of Computer
Vision 50, 345–362.
Pednekar, A., Kurkure, U., Muthupillai, R., Flamm, S., Kakadiaris, I.A.,
2006. Automated left ventricular segmentation in cardiac mri. Biomedical
Engineering, IEEE Transactions on 53, 1425–1428.
Petitjean, C., Dacher, J., 2011. A review of segmentation methods in short
axis cardiac MR images. Medical Image Analysis 15, 169–184.
Pluempitiwiriyawej, C., Moura, J., Wu, Y., Ho, C., 2005. STACS: New active
contour scheme for cardiac MR image segmentation. IEEE Transactions
on Medical Imaging 24, 593–603.
Queirós, S., Barbosa, D., Heyde, B., Morais, P., Friboulet, D., Claus, P.,
Bernard, O., D’hooge, J., 2013. A fast fully automatic segmentation of the
myocardium in 2d mr images, in: Proceedings of Functional Imaging and
Modeling of the Heart (FIMH2013), pp. 71–79.
Radau, P., Lu, Y., Connelly, K., Paul, G., Dick, A., Wright, G., 2009. Evaluation framework for algorithms segmenting short axis cardiac MRI. The
MIDAS Journal - Cardiac MR Left Ventricle Segmentation Challenge .
Schaerer, J., Casta, C., Pousin, J., Clarysse, P., 2010. A dynamic elastic
model for segmentation and tracking of the heart in MR image sequences.
Medical Image Analysis 14, 738–749.
42
Shi, Y., Karl, W., 2008. A real-time algorithm for the approximation of
level-set-based curve evolution. IEEE Transactions on Image Processing
17, 645–656.
Ulén, J., Strandmark, P., Kahl, F., 2013. An efficient optimization framework for multi-region segmentation based on lagrangian duality. IEEE
Transactions on Medical Imaging .
Uzunbas, M., Zhang, S., Pohl, K., Metaxas, D., Axel, L., 2012. Segmentation
of myocardium using deformable regions and graph cuts, in: IEEE International Symposium on Biomedical Imaging (ISBI), IEEE. pp. 254–257.
Wijnhout, J., Hendriksen, D., Assen, H., Geest, R., 2009. LV challenge LKEB
contribution: Fully automated myocardial contour detection. The MIDAS
Journal - Cardiac MR Left Ventricle Segmentation Challenge , 683.
World Health Organization, 2012. Cardiovascular diseases. Fact sheet 317.
43
Denis Friboulet