A review of heart chamber segmentation for structural and functional analysis using cardiac magnetic resonance imaging

Cardiovascular magnetic resonance (CMR) has become a key imaging modality in clinical cardiology practice due to its unique capabilities for non-invasive imaging of the cardiac chambers and great vessels. A wide range of CMR sequences have been developed to assess various aspects of cardiac structure and function, and significant advances have also been made in terms of imaging quality and acquisition times. A lot of research has been dedicated to the development of global and regional quantitative CMR indices that help the distinction between health and pathology. The goal of this review paper is to discuss the structural and functional CMR indices that have been proposed thus far for clinical assessment of the cardiac chambers. We include indices definitions, the requirements for the calculations, exemplar applications in cardiovascular diseases, and the corresponding normal ranges. Furthermore, we review the most recent state-of-the art techniques for the automatic segmentation of the cardiac boundaries, which are necessary for the calculation of the CMR indices. Finally, we provide a detailed discussion of the existing literature and of the future challenges that need to be addressed to enable a more robust and comprehensive assessment of the cardiac chambers in clinical practice.


Introduction
Cardiovascular diseases (CVDs) consistently rank among the top major causes of morbidity and mortality. In 2008, 17.3 million people died due to CVDs worldwide, accounting for 30 % of total deaths [1]. Of these cases, about 7.3 million were due to coronary heart disease, and 6.2 million were due to stroke [2]. However, partially due to the aging population, this number keeps increasing. It is predicted that by the year 2030, a population of 23.3 million people will be killed by CVDs all over the world [1,3]. Consequently, major developments continue to be made in cardiovascular research and practice for improved early diagnosis of cardiac diseases.
In particular, magnetic resonance imaging (MRI), otherwise known as CMR (cardiovascular magnetic resonance), has become a key image modality in clinical practice due to its unique capabilities for non-invasive imaging of the cardiac chambers and great vessels [4]. A wide range of CMR sequences and protocols have been developed to assess various aspects of cardiac function, and significant advances have also been made in terms of imaging quality and acquisition times [5]. Furthermore, a lot of research has been dedicated to the development of global and regional quantitative CMR indices, which help distinguish between pathology and health.
The goal of this review paper is threefold. Firstly, we will review the various functional indices that have been proposed thus far in the literature. These will be presented in detail for each cardiac structure, including their definitions, calculation requirements, the exemplar Abstract Cardiovascular magnetic resonance (CMR) has become a key imaging modality in clinical cardiology practice due to its unique capabilities for non-invasive imaging of the cardiac chambers and great vessels. A wide range of CMR sequences have been developed to assess various aspects of cardiac structure and function, and significant advances have also been made in terms of imaging quality and acquisition times. A lot of research has been dedicated to the development of global and regional quantitative CMR indices that help the distinction between health and pathology. The goal of this review paper is to discuss the structural and functional CMR indices that have been proposed thus far for clinical assessment of the cardiac chambers. We include indices definitions, the requirements for the calculations, exemplar applications in cardiovascular diseases, and the corresponding normal ranges. Furthermore, we review the most recent state-of-the art techniques for the automatic segmentation of the cardiac boundaries, which are necessary for the calculation of the CMR indices. Finally, we provide a detailed discussion of the existing literature and of the future challenges that need to be addressed to enable a more robust and comprehensive assessment of the cardiac chambers in clinical practice. applications for cardiovascular diseases for each index, and the corresponding normal ranges. Subsequently, and because the calculation of these indices requires delineation of the cardiac boundaries, we will review the most recent state-of-the art techniques for the automatic segmentation of the various cardiac structures. These techniques have the advantage of delineating these boundaries of the heart more rapidly and objectively than clinical experts with manual contouring. We will focus on presenting the more practical properties of these segmentation techniques, such as in relation to their specific roles, and on the imaging materials these techniques need to accomplish the segmentation (e.g. long-axis vs. shortaxis images). Finally, we provide a detailed discussion of the existing literature and of the future challenges that need to be addressed to enable a more robust and comprehensive assessment of the cardiac chambers in clinical practice.
In comparison to existing reviews [6][7][8][9][10][11][12][13], our survey is more comprehensive in many aspects, including: • We review both the segmentation techniques and the structural and functional indices that use the extracted boundaries in their calculations. This will provide a better understanding of the role of cardiac segmentation to the final clinical assessment. • We also list the quantitative evaluations for the accuracy of both the segmentation and functional analysis to provide the readers with an overview of the level of performance of the existing techniques so far. • We provide a more comprehensive review of segmentation techniques for all cardiac chambers, including the left ventricle (LV), the right ventricle (RV), the left atrium (LA), and the whole heart. • We also include the use of long-axis images in cardiac segmentation as they play an important role in clinical use of CMR. • Finally, we include newly emerged concepts in machine learning-based cardiac image analysis, such as direct estimation of cardiac function [14][15][16][17].
This paper is organised as follows. In this introductory section, we will describe the anatomy of the heart, followed by a presentation of the CMR protocols. Subsequently, in section and Table 1, we will present in detail the existing indices of cardiac structure and function, which we will organize per cardiac chamber. In section and Tables 2, 3, 4, 5, and 6, we will then describe the most recent segmentation techniques that can be used to extract automatically the boundaries to be used to compute the functional indices of interest. Finally, we will conclude with a discussion of current cardiac segmentation challenges and future perspectives.

The anatomy of the heart
In this section, we briefly describe the anatomy of the heart to help readers establish a better association between the outcomes of various functional analysis methods and the actual structure of the heart (see Fig. 1). Essentially, the heart provides the blood circulation system with indispensable pressure. By contracting and relaxing in turns, it transports blood to different parts of the body through the vessels. The septum separates the heart into two halves that consist of an atrium and a ventricle. The left atrium (LA) and left ventricle (LV) are partitioned by the mitral valve, while the right atrium (RA) and the right ventricle (RV) are partitioned by the tricuspid valve. The semilunar valves are located between the pulmonary artery or the aorta and the ventricle. The RA recycles the low-oxygen blood while the RV delivers it to the lung. After it is oxygenated, the blood flows into the LA, while the LV pumps it to the rest of the body. The myocardium, the muscular tissue of the heart has an inner and outer border: the endocardium and the epicardium, respectively.

MRI protocols
Since pathological changes are related to abnormal structural and physiological indices, experts are seeking for a more accurate diagnosis or risk stratification of CVDs based on quantitative anatomical or functional information. Various imaging techniques for clinicians have been developed. Unlike radioisotopes, computed tomography, and angiography, CMR is a non-invasive imaging technique that is capable of generating images in decent resolution without ionising radiation. Compared to the traditional echocardiography, CMR does not suffer from speckle artifacts and produces good contrast between the different soft tissues. Images can be obtained in any orientation allowing for images to be acquired in specific anatomical planes. Owing to these properties, scientists have been developing diverse protocols providing varying information. Among them, cine CMR, flow CMR, tagged CMR, late gadolinium enhancement (LGE), and perfusion CMR are the mainstream applications.
Cine CMR aims at providing fine spatiotemporal resolution with high contrast between the tissues. One sample normally contains 20-30 consecutive frames, corresponding to 20-30 time points in the cardiac cycle. Each frame has multiple slices from base to apex (Fig. 2)-typically between 10 and 15. Generally, the images are captured along two axes: the long axis and the short axis views (Fig. 3). The long axis (LAX) goes across the LV from base to apex. The short-axis (SAX) slices are perpendicular to the LAX. Because the frame sequence loop reflects the dynamic process of a complete cardiac cycle during a  [19], cine CMR is widely employed in calculating global functional indices such as stroke volume and ejection fraction. Flow CMR is a velocity-encoded protocol, based on the principle that the pulse phase shifts of moving protons are proportional to their velocity along the magnetic field gradient direction [20]. Therefore, the motion of a tissue will generate an MRI signal variation. Flow CMR commences with a reference MRI scan, which uses stationary spins. Afterwards, a number of scans are produced to encode the velocity information by adjusting the direction of the gradient from +180° to −180°. In consequence, moving protons show different intensities from the initial scanning: the brighter areas on phase contrast images are drawn by the protons moving along a certain direction; the darker areas have protons going towards the opposite way; the regions where the stationary protons rest appear to be grey. This property gives flow CMR the advantage in measuring the cardiovascular flow and strain rate.
Tagged CMR builds a spatial line or grid pattern on the myocardium, which is then followed over the cardiac cycle to estimate cardiac motion. This is based on the received signal from myocardium by modulating saturated magnetisation inside the ventricular wall [21][22][23]. The dark pattern, which stands at a fixed position on the myocardial tissue, is usually added at end-diastole using radio-frequency excitation and gradient impulses before image acquisition. During the contractile cycle, the dark patterns will move with the tagged tissue, as shown in Fig. 4. By tracking the displacement and distortion of those saturated patterns marked on the tissue, researchers can compute the precise myocardial deformations or reconstruct the wall motion easily. Therefore, tagged CMR is efficient in regional assessments such as for the estimation of myocardial strain and torsion. The limitation of this promising protocol is that the markers always fade inevitably before the whole cycle ends. Also, the existence of the grids brings difficulties to automatic cardiac border identification. The progress and challenges of MRI tagging have been summarised in [23][24][25].
Displacement Encoding with Stimulated Echoes (DENSE), which combines the merits of flow CMR and tagged CMR, intends to map the myocardial displacement in high spatial resolution over long periods of cardiac cycles [26] without having serious fading. Different from flow CMR, DENSE uses stimulated echo to modulate the phase, which aims at capturing the emerging displacement between the second and third radio-frequency pulses. This technique can be applied to abnormal contraction diagnosis, myocardium deformation, and motion analysis. However, the imaging is usually time-consuming.
and HD are point-to-point, point-to-curve, point-to-surface, curve-to-curve, and Hausdorff distance, respectively; 60 (27)  Strain-Encoded (SENC) CMR is designed to obtain longitudinal strain straightforwardly, without dealing with displacement or velocity [27]. The dense estimation of longitudinal strain is achieved by processing the tag information extracted from two short-axis images, whose planes are orthogonal to the strain imaging orientation. The tags express the local strain as intensity and their surfaces are set to be parallel to the short-axis images. The short-axis images are generated with two phase-encodings based on slice selection. It has been shown that SENC is a reliable tool to quantify regional myocardial systolic and diastolic function [28].
Perfusion CMR produces contrast-enhanced images by injecting contrast agent (typically gadolinium-based chelates) [29]. The contrast agent travels through the vessels or lymphatic system as the blood flows past and finally reaches the target tissue, which leads to a variation in signal intensity of the agent. A fast scanner with high temporal resolution is responsible for monitoring this signal fluctuation and then sketching sequential images. Perfusion CMR is used for diagnosing ischemic heart disease, for which the myocardium is associated with less blood movement (see Fig. 5). However, perfusion CMR suffers from quantitative analysis degradation introduced by artifacts, ranging from surface coil inhomogeneity, dark rim to motion artifacts. Many researchers have proposed solutions to these inherent weaknesses [10].
Late Gadolinium Enhancement (LGE) CMR is an important technique for the estimation of scar tissue in the myocardium [31]. This technique acquires images (6 mm SAX slice thickness with 4 mm gap for contrast-enhancement match), followed by an injection of 0.10-0.15 mmol/kg intravenous gadolinium [32]. After a delay of 10-20 min, by using the inversion-recovery fast gradient echo (IR-FGE) pulse sequence, LGE images are collected from the same position with a decent spatial resolution. Normally, the contrast agent cannot enter the myocardial cells. In abnormal cases, gadolinium may gather in extracellular space or even break into the cells due to cell membrane rupture. As a consequence, healthy tissues stay dark while infarcted parts appear brighter on the image (see Fig. 5). Therefore, LGE can be very useful in examining injured tissue with infarction or scars.

Indices of cardiac function
The existing indices of morphology and function can be divided into two categories: global and regional. Global indices include chamber volumes, stroke volume, ejection fraction, cardiac output, and myocardial mass. Regional or local indices cover myocardial wall thickness and thickening. Strain analysis can be either global or local. We have listed in Table 1 many indices, which are frequently used in CMR research, for cardiac structural and functional analysis. We provide the parameter definitions, requirements for calculations, use cases for cardiovascular diseases, and published normal ranges.

LV quantification
The LV is the most investigated chamber in cardiac segmentation and structural and functional analysis due to its central role in blood circulation. It has relatively thick myocardial tissues that give blood circulation enough pressure. LV parameters can be abnormal in many CVDs, such as in hypertension or after myocardial infarction.
Left Ventricle End-diastolic and End-systolic Volumes (LVEDV and LVESV) are measurements of the amount of blood in the chamber, encompassed by the myocardial tissue, when the heart muscle is relaxed (LVEDV) or contracted (LVESV). The contour on the basal slice from the images stack is drawn on the aortic valve cusps level, resulting in an inclusion of the outflow tract as part of the LV volumes. There is no consensus as to whether to include or exclude papillary muscles from the LV blood pool [33,34,[41][42][43][44][45][46][47].
Left Ventricle Stroke Volume (LVSV) is the amount of blood ejected from the heart during each contraction. LVSV is the difference between the LVEDV and LVESV.
Left Ventricle Mass (LVM) measures the myocardial tissue. The volume of the myocardium can be obtained by subtracting the endocardial volume from the volume of within the epicardial border. Subsequently, the mass is the product of myocardial volume and the muscle density. LV mass is prognostic in hypertension [48,49]. Left Ventricle Ejection Fraction (LVEF) quantifies the quantity of blood pumped out of the heart in each beat as a percentage. It divides the LVSV by the LVEDV. Normal ranges for LVEF are gender-and age-dependent and also dependent on the analysis approach chosen (e.g. include or exclude papillary muscle from LV volumes). Reduced LVEF is a common final pathway in many CVDs (e.g. dilated cardiomyopathy, remodelling after myocardial infarctions). Hyperdynamic LV systolic function as seen by high LVEF can often be seen in LV hypertrophy (e.g. hypertrophic cardiomyopathy) [50].
Cardiac Output (LVCO) refers to the amount of systemic flow per minute. It can be estimated by multiplying the LVSV with Heart Rate (HR), which denotes the heartbeat frequency (beats per minute). LVCO is often normalised by the Body Surface Area (BSA), and then referred to as Left Ventricle Cardiac Index (LVCI). In patients with congestive heart failure, the LVCO and LVCI are reduced [51].
Left Ventricle Wall Thickness is the thickness of the myocardium typically measured on end-diastolic images in SAX view. The papillary muscles and trabecular tissues are usually excluded. First, both epicardial and endocardial boundaries   [24] are identified. Afterwards, a centre point or a centreline with reference points is specified to help compute the mean distance between the epicardial and endocardial contours [52], as displayed in Fig. 6. For regional analysis, researchers are encouraged to use the 17-segment model [53] (see Fig. 7). Wall thickness may be globally increased (and hence LVM is typically also increased) in conditions with increased afterload, such as hypertension. Some conditions lead to regional increased wall thickness (with or without increased LVM) typically referred to as showing asymmetric hypertrophy, such as seen in hypertrophic cardiomyopathy. In contrast, myocardial infarctions can lead to regional thinning in the area of infarct as a consequence of cardiac remodelling.
Left Ventricle Wall Thickening reflects the change of myocardial wall thickness during systole expressed as a percentage. Wall thickness may be employed to quantify regional dysfunction , such as those seen in myocardial ischemia or after myocardial infarction [54].
Left Ventricle Strain (LVS) indicates the degree of deformation of the ventricles, while Left Ventricle Strain Rate (LVSR) is the deformation rate. The required parameters can be received by echocardiography, such as tissue Doppler. Some MRI techniques, for example, SENC, DENSE, and tagging, can be complementary. LVS may play an important role in evaluating myocardial infarction, ischemia, and ventricular dyssynchrony [55].

RV quantification
The RV consists of the apical body, the inflow tract and the outflow tract. The existing RV functional indices quantify the amount of blood being transported to the lung in different forms. As can be seen from Table 1, the definitions of most RV indices, including Right Ventricle End-diastolic and End-systolic Volume (RVEDV and RVESV), Right Ventricle Stroke Volume (RVSV), Right Ventricle Ejection Fraction (RVEF), and Right Ventricle Cardiac Output (RVCO), are fundamentally similar to their LV counterparts. The papillary muscles and trabecular tissues are neglected in endocardial contour depiction. RV volumes may be increased in a number of conditions, including cardiac shunts, certain valve diseases, or pulmonary hypertension [56]. RVEF may also be decreased after myocardial infarctions including parts of the RV. Although efforts have been made to extract the boundaries of endocardium and epicardium simultaneously, the mass is still not regularly evaluated because the myocardial wall of RV is 3-6 times thinner than the wall of LV.

LA quantification
The LA plays an important role for the modulation of LV blood filling. The LA has a relatively complex geometric structure, surrounded by the aorta, pulmonary veins, and arteries. A dilated LA has prognostic value for cardiovascular death [57], stroke [58], congestive heart failure, and atrial fibrillation [59].
Left Atrium Volume (LAV) assessed when largest during ventricular systole just before mitral valve opening has been demonstrated as a reliable predictor of cardiovascular outcomes [60], including LV diastolic dysfunction [61], incident atrial fibrillation [62], ischemic stroke [63], hypertrophic cardiomyopathy [64], and lone atrial fibrillation [65]. Similar to the calculation of LV volumes, computational methods can be based on either contiguous summation or geometry assumption. In automatic segmentation, the confluence of pulmonary veins and the LA appendage (area under the mitral valve annulus) are abandoned. Three volumetric parameters, the maximum LA volume (LAV max ), the minimum LA volume (LAV min ), and the pre-atrial The reservoir describes the filling of LA in ventricular systole. It is modulated by the LV contraction, RV systolic pressure, LA relaxation, and LA chamber stiffness. When blood flows to the LA from the pulmonary veins, the mitral valve is closed and LAV increases to maximum. Total Emptying Volume (LAEV) and Total Emptying Fraction (LAEF) are used to quantify the total amount of blood the LA can pump into the LV.
The conduit function involves LV relaxation and LA afterloads. In early ventricular diastole, LAV grows and the atrial blood is suctioned by the LV. The LA acts like a passive conduit and three indices have been proposed to assess its function: the Passive Emptying Volume (LAPEV), the Passive Emptying Fraction (LAPEF), and the Conduit Volume (LACV), which indicate the amount of blood travelling from the pulmonary veins to LV.
The booster pump, also called contractile function, quantifies the amount of blood being pumped into the LV during LA contraction. It is modulated by the LV compliance, LA afterload, LA preload, and intrinsic LA contractility. In late ventricular diastole, the LA pumps all the remaining blood to the LV actively, and thus the LAV decreases to minimum. Corresponding measurements include Active Emptying Volume (LAAEV) and Active Emptying Fraction (LAAEF).

RA quantification
RA is not routinely assessed in CMR. However, its enlargement may indicate heart failure, as well as valvular and congenital diseases [40]. The filling of the RA is related to the functions of the RV. RA volume (RAV) indexed to BSA can predict pulmonary hypertension [66,67] and chronic systolic heart failure [68]. Because direct measurement of RAV can be time-consuming, scientists primarily do the estimation based on single or bi-plane area-length methods through inspecting RA areas on two-chamber and fourchamber LAX viewing (see Table 1).

Cardiac segmentation
In order to calculate the CMR structural and functional indices listed in this section, the boundaries of the heart chambers are necessary. However, delineating the heart manually on multiple slices and frames requires lots of time. Furthermore, this is subject to well-established intraand inter-subject variability. This has motivated engineers to develop automated cardiac segmentation techniques that can rapidly, objectively, and accurately extract the chamber boundaries from CMR in clinical practice. Although MRI provides decent soft tissue contrast among different protocols, accurate cardiac segmentation remains a great challenge for the researchers due to inevitable imaging inhomogeneity and high anatomical variability, as well to the inherent geometric and dynamic complexity of the heart. In this section, we will describe the existing segmentation methods published in popular journals and conferences from the year 2000, by focusing on their principles, functions, advantages, and limitations.

Segmentation methodologies
Generally, semi-automatic or fully automatic segmentation techniques fall into two categories: (1) image-driven approaches without or with weak prior models and (2) model-driven approaches based on strong prior knowledge. Training data are examples with ground-truth. Imagedriven methods identify the pixels or voxels belonging to the blood pool, myocardium, or appendage by visiting their intensity differences. Typical image-driven techniques consist of thresholding, region-growing, clustering, pixel or voxel classification, and active contour or surface. Strong prior knowledge, including cardiac atlases and statistical shape models (SSMs), make use of statistical information extracted from manually annotated training data that describe for example averages and modes of variations of the cardiac chambers. For the rest of this section, we briefly discuss each of these techniques.
Thresholding can be used to localise the region of interest (ROI), such as the blood pool or myocardium, based on analysing the intensity histogram. The latter is usually constructed as a discrete distribution of pixel intensities (counts vs. values). Then a threshold value, which corresponds to a specific intensity, is to divide the histogram into sub-intervals containing distinctive modes. The pixels having intensities in a same interval may belong to a certain type of tissue. This method is only effective when significant intensity diversity exists between the target and background areas. However, in some cases, the intensity of different tissue types overlap. Therefore, thresholding is often used as a pre-processing step and further combined with other segmentation techniques.
Region-growing starts with choosing one or multiple seed points in MR images in a selected region such as the myocardium. Afterwards, the initial region begins to grow by searching similar pixels nearby or inside a neighbourhood. If a pixel (x, y) meets the designed criterion, it will be allocated to region R i in the ith step: R i+1 = R i ∪ (x, y). When none of the surrounding pixels qualify, the region stops growing as it may have reached the boundary of the tissue. Merge behaves alike, but instead of judging single pixels, it combines similar small regions. While split performs in the opposite way, it shatters the region or suspends the membership if a sub area differs significantly from the rest of the area. Because of the continuity of growth, region-growing, or split and merge techniques, can often lead into over segmented target tissues, leaking into fragments of irrelevant parts. For instance, the aorta and cavity may have close intensities on basal SAX slices and cannot be distinguished using only thresholding. Watershed [69] combines thresholding and region merging by calculating the image gradient map and setting a threshold on the magnitude of the gradients. If a pixel and its adjacent neighbours all have similar magnitudes below the specified threshold, they are merged. Watershed is known to result in over-segmentation and poor performance in noisy regions due to the reliance on image gradient.
Pixel or voxel classification groups pixels in 2D or voxels in 3D in feature space. Patch-based features contain pixel intensity or textural appearance information. Unsupervised clustering, which is non-parametric, does not require manually labelled training data. Typical methods include K-means clustering and expectation-maximization (EM) [70]. K-means clustering randomly chooses K features as the initial centroids and classifies all other features according to their distances to the centroids, then calculates new centroids of those categories. These steps are repeated until centroids are converged and no longer change. EM finds the maximum likelihood (ML) or maximum a posteriori (MAP) estimates of parameters of a statistical model. For cardiac segmentation, a common model is the Gaussian Mixture Model (GMM), in which each tissue histogram follows a Gaussian distribution. Every pixel is classified to the region that maximises its corresponding class conditional probability. Supervised classifiers, such as K-nearest neighbour (KNN), random forest, and neural network, need manually labelled training data. In these methods, the training data and their associated labels are regarded as examples from which the parameters of the classifiers are learned by minimising a risk function that pertains to misclassification of the training labels. Each test pixel or voxel can be accordingly classified afterwards using the learned classifiers. However, annotating training data involves user interaction and the performance of these classifiers often depends on the quality of training samples. If the training and testing datasets statistically deviate by a large extent, the classification performance declines significantly. Moreover, classification-based segmentation methods often ignore the spatial dependencies of the local features. The advantage of supervised techniques is that they are trainable to segment more accurately, provided that the expert knowledge is properly employed in the classifier.
Active contour-based methods or snakes [71] search for chamber walls, instead of directly classifying the regions. A curve parametrised C(s) = (x(s), y(s)) where s denotes a free parameter, is morphed locally towards target boundaries by minimising a predefined energy. In order to achieve a better result, many researchers have designed different energy functions. Generally, the energy E can be written as E = ∫ E in (C(s)) + E ex (C(s)) + E c (C(s))ds, where E in indicates the internal force that aims at retaining the topology and smoothness of the curve, E ex is the external force pushing the curve to target boundary and E c stands for additional constraints. The last aims at improving convergence or penalising unwanted shape irregularities. Segmentation based on active contour may need user interaction, for example, roughly drawing or placing a contour for initialisation. An improvement over the traditional form of active contour is achieved using level-set formulation, in which the curve implicitly defined as the zero level-set of a higher dimension function [72][73][74]. The level-set can handle larger shape updates, when the morphology of the curve has to be evolved significantly. To segment the hearts in a multi-phase fashion, the converged segmentation can be propagated into images in the subsequent time points for a better initialisation, removing user interaction.
Direct estimation has been recently proposed as a means to estimate functional indices such as chamber volumes without segmenting MRI slices [14][15][16][17]. This approach has proved to be an efficient tool in myocardial abnormalities detection [14]. It uses regression-like models trained with discriminative image representations to estimate the ventricle volumes from image information. Different from the pixel or voxel classification problem, the whole or parts of the image act as global input features to establish similarities to reference samples with known functional indices.
Atlas-based segmentation methods rely on the spatial probability patterns of various tissue types of a typical heart. To segment a test case, the image is registered to the atlas, which serves as the prior information for the pixel labels given their locations. It has been demonstrated that multi-atlas-based segmentation methods outperform singleatlas approaches remarkably in terms of accuracy in other applications [75][76][77][78][79].
Statistical shape modelling introduced by Cootes et al. [80], is a powerful tool for cardiac quantitative assessment. Given a population of corresponding points or vertices from myocardial surface meshes, a mean shape is extracted and a set of variation modes can be built using principle component analysis (PCA). Then any novel shape from an individual can be represented as the mean shape varied by a linear weighted combination form of the PCA modes. This representation is called the point distribution model (PDM). For segmentation based on ASM [80] or AAM [81], the linear model is matched to the test image by matching landmark through global transform and finding an optimal 1 3 loading vector of PCA modes. This is usually achieved iteratively, updating one set of parameters at a time.
Despite popularity of ASM/AAM-based segmentation approaches, their training demands identifying a dense set of corresponding landmarks over the training population. Furthermore, the fitting procedure can be computationally slow and prone to local minimums. Models constructed from healthy populations many not fit pathological hearts well, as the model becomes too restrictive. Such large interclass variation provokes a need for training more generalizable models, as well as more sophisticated fitting processes. Nevertheless, model-based methods are still considered to be promising routes to accurate MRI segmentation as they are capable of preserving anatomical spatial knowledge while segmenting the heart. In the following sections, we detail the specific roles and efficacies of these techniques in quantitative analysis of LV, RV, LA, bi-ventricle, and whole heart segmentation.

LV
Among various compartments of the heart, LV has been studied the most extensively, as it pumps the blood into other parts of the body. A relatively thick myocardial wall leads to the popularity in the research of regional assessment, such as apical, middle, or basal wall thickness, local deformation, and myocardial strain. Only those approaches that process multiple phases (including end-diastole and end-systole) can measure LVEF, LVSV, and LVCO, since LVEDV and LVESV are known.
Methods able to extract both the endocardium and epicardium can be used to calculate the LVM and wall thickness. For wall thickening, the centreline method performs better, as its radial opponent often overestimates the distances between the contours of epicardium and endocardium. This is caused by the initial hypothesis of the radial method, which assumes the shape of myocardium as a circle [82].
LVS can be analysed by tracking myocardial motion, since regional muscular displacement and temporal information are both required. The global strain analysis in 3D begins with creating Cartesian coordinates. The extent of deformation, described as the change of length from an initial or reference status, can be calculated using Lagrangian or Eulerian formulae [83]. Because the heart deforms along different directions in Cartesian coordinates in 3D simultaneously, a matrix called a tensor is created to describe the process. For regional analysis, the local coordinates are with three mutually perpendicular axes: the radial (perpendicular to the epicardium and towards the outside), the longitudinal (tangent to the epicardium and towards the base), and the circumferential (according to the right-hand rule, from radius to longitude) axes. Therefore, the spatial orientation of three axes varies with the voxel position in the myocardium.
Theoretically, LV volumes can be estimated with any 2D segmentation outcome on SAX or LAX slices, not necessarily 3D, by making use of provided volumetric calculation methods in Table 1. However, when only SAX slices are in use, the segmentation must be completed on a stack of multiple slices from the base to the apex. We list all the LV segmentation techniques described in this section in Table 2.

Thresholding and region-growing
Thresholding is often integrated with region-growing. Lee et al. [84] and Codella et al. [85] use region-growing to find the full-blood LV region. They automatically identify a seed point by taking the pixel with the lowest energy in a window during slice propagation. Then in order to prevent the segmented LV region from diffusing to epicardial fat, fluids, and RV, they use an iterative thresholding mechanism that discovers a lower bound of myocardial intensity. Huang et al. [86] employ thresholding to distinguish the blood pool from the myocardium, followed by radial region-growing and extraction of convex hulling to identify the endocardial and epicardial boundaries. Lu et al. [87] apply thresholding to convert a ROI to a binary image for LV localisation and endocardial contour detection (Fig. 8), followed by region-growing to segment the LV epicardium. Ammar et al. [88] take the binary image produced by thresholding as the initial mask for a level-set segmentation method to extract the endocardium. Queiros et al. [89] perform class decomposition following thresholding step to search for the LV centroid. The method sets two thresholds for myocardium and cavity histogram in an EM algorithm to extract the endocardial contour. Kurkure et al. [90] localise LV in the thresholded image by finding a binary component that is closest to the intersection cross-hair generated by LAX vertical and four-chamber view projection in ED phase on a SAX slice. They have also proposed a novel fuzzy connectedness region-growing method taking the spatial adjacency, intensity homogeneity, and multiclass features into consideration. The myocardial boundaries are extracted by dynamic programming, which is an optimal path finding solution of overcoming obstacles such as papillary muscles or trabeculae carneae extrusion, and low liver-to-myocardium tissue contrast. Cousty et al. [91] extract epicardium by developing a spatial-temporal gradient computation for watershed cuts. It is noteworthy that approaches in [84][85][86][87]89] start their segmentation from mid-ventricular SAX slices, which might involve user interaction, and then propagate their initial results to other slices as prior knowledge. Also, the test image is usually mapped to the polar coordinate since LV roughly has a circular shape [84-87, 89, 90]. Furthermore, by making use of thresholding on the LV blood pool, the papillary muscles and the trabeculations can be easily outlined [84][85][86][87] due to their intensity diversity with surroundings. Thus, cardiac functional analysis such as LVEDV, LVESV, or LVM estimation be can varied by including or excluding papillary muscles and trabeculations or not, depending on the index definition and requirement of clinicians.

Pixel or voxel classification
Classification-based methods in cardiac segmentation have also been thoroughly studied. Jolly [92] and Hu et al. [93] propose to classify regions by a 3-GMM with EM, based on the intensity histograms. Jolly [92] separates the muscle, air, blood, and fat, as presented in Fig. 9. Hu et al. [93] separate the muscle, the blood and the background. Pednekar et al. [94] fit a 5-GMM to the intensity histogram of the blood pool, the lung filled with air, myocardium, the region between the blood and myocardium, and the region between the air and myocardium. The EM algorithm is initialised through K-means clustering. Queiros et al. [89] use a 2-GMM. Some classifiers label the features from different tissues without making assumptions on intensity histogram distribution. To label the regions of the lung, the myocardium, and the blood pool, Stalidis et al. [95] make use of a neural network classifier, which is trained via a small number of representative tissue points. The input features of their classifier are the pixel position, pixel intensity, and slice location. Folkesson et al. [96] have presented that a trained KNN classifier is competent for the classification of the LV cavity, myocardium, and background, based on a feature selection scheme. The latter finds the most discriminative features for the pixel classifier or model fitting, aiming at increasing computational efficiency without degrading its accuracy. Bai et al. [97] have shown that the support vector machine (SVM) outperforms KNN in label fusion in a multi-atlas-based cardiac segmentation framework.

Active contours
Active contour, or deformable model, is one of the most widely applied techniques in heart segmentation. The breakthrough usually comes with the design of the energy function, such as using anatomical assumptions as the constraints on the level-set methods [98,99]. Paragios [98] propagates coupled endocardial and epicardial contours on SAX slices, where the edge, region, and anatomical  Gaussian distributed components representing the air, myocardial muscle, and blood/fat compartment; c the output image with classified pixels in different labels [92] 1 3 constraints are pre-defined. The edge constraint is used to push the curves to the myocardial walls. The region intensity criterion makes the model less sensitive to initial conditions. The GVF snake, a parametric active contour that overcomes the difficulty in evolving the curve to the boundary concavities, is introduced into LV segmentation in these works [100][101][102]. Wu et al. [103] claim that the gradient vector convolution (GVC) snake also conquers local minima such as artifacts and papillary muscles. Kaus et al. [104] integrate strong prior knowledge, in the form of PDM, into deformable contours by extending their internal energy, leading to an increase in the robustness of the model. They have considered the inter-spatial relationship of the inner and outer boundaries as well, which compensates for the error produced by incorrect feature detection. Lynch et al. [105] employ the probability density function as the prior information, which is created by a set of manually segmented boundaries on binary images. Furthermore, the evolutions of endocardial and epicardial curves are coupled by an extra level-set constraint. The deformable models proposed in [106,107] and the level-set approaches presented in [108,109] incorporate the spatial-temporal LV activation as prior knowledge and track the epicardium/endocardium boundaries on SAX slices in a complete cycle. A typical tracking result is shown in Fig. 10. The constraints in [106,107] are parameterised by Fourier descriptors. Moreover, an approach for the recognition of intra-ventricular dyssynchrony (IVD) is proposed in [106], where the non-uniform contraction of the ventricular walls brought by the activation delays can be discovered. Jolly [92] makes use of a deformable model to improve further the outcome of EM-based region segmentation. Chen et al. [110] propose to apply deformable models to LVS analysis on SAX slices in DENSE MRI. Their model is driven by minimising an energy function that consists of model intensity, edge attraction, shape prior, contours interaction, and smoothness. The shape prior can eliminate the concavities with negative curvatures in order to remove the papillary muscles from the ventricular walls. Huang et al. [111] have invented a novel deformable model called Metamorphs, whose energy functions are predefined on the distance maps of the object shape and its border. Metamorphs is not particularly designed for MRI cardiac segmentation while it outperforms the GVF snake as a result of better robustness to inferior initial conditions. Based on the motion trajectories in DENSE, representing the movement of myocardial wall between two consecutive phases, the initial manually drawn endocardial and epicardial contours can be propagated slice by slice to other frames [112,113]. The method proposed in [112] is applicable to both SAX and LAX images. Besides DENSE, Chen et al. [114] have also used tagged MRI to derive LVS. In their work, Gabor filters search the tag intersections. Through matching these intersections, the method is able to track the myocardial motion. The deformable model refines the tracking and displays a dense displacement map. Kermani et al. [115] draw a dense displacement map by fitting a 3D active surface model to an initial sparse displacement map, which is built by establishing point correspondence in cine images. Motion tracking makes LVS easier to be analysed, because the myocardial displacement and temporal scale are known at the meantime. The authors produce visualisations of LV strains in 3D (see Fig. 11). Khalifa et al. [116] measure the wall thickness and thickening with a stochastic speed function based level-set technique extracting inner and outer myocardial walls first. Subsequently, the points on the Fig. 10 LV epicardium (left) and endocardium (right) tracking: contours propagate through short-axis slices on all phases in a complete cardiac cycle [106] inner contour and the outer contour are paired. The Euclidean distance between each pair is the wall thickness. Wei et al. [117] use the myocardium contours from cine MRI as prior knowledge to guide the meshing of endocardium and epicardium, which are generated by contour registration to move towards the inner and outer edges in SAX and LAX slices in LGE. Grande et al. [118] model the image likelihood by sampling the intensity and gradient of pixels inside the myocardium or at the boundary of myocardium in different regions. After that, they create a Markov Random Field (MRF) to incorporate the prior and the likelihood models. The prior keeps the curve smooth and excludes the papillary muscles. The deformable model estimates the walls based on the MRF along the SAX radial direction.

Strong prior based techniques
Different from image-driven techniques, model-based approaches exploit strong prior knowledge such as by encoding the specific shape variability of the LV, instead of making simple assumptions on the boundaries. By taking advantage of the statistical shape information, segmentation becomes more robust to image noise by restricting the outcome to valid instances statistically. Mitchell et al. [119] introduced early application of 3D-AAM to LV segmentation in 2002. The method showed worthy results in quantifying LVV epi , LVV endo , and LVM on SAX volumes. Assen et al. [120] proposed a 3D-ASM segmentation method (SPASM) that can operate on sparse MR images scanned in arbitrary orientations. In most cases, the automated LV segmentation approaches require a stack of parallel SAX images. While SPASM can perform on a datasets of two orthogonal radial LAX slices, four radial LAX slices with a 45-degree angle between two neighbours, 11 equally spaced SAX slices, four SAX slices (one apical, one mid-cavity, two basal), or a combination of two LAX and two SAX slices. The processing pipeline of this method is implemented and shown in Fig. 12. Lekadir et al. [121] improve the 3D-ASM by incorporating an additional shape prior, which is invariant to transforms including translation, rotation, and scaling. This prior is used to detect and correct outliers, thus leading to more robust results. Andreopoulos and Tsotsos [122] use a hierarchical 2D-ASM that incorporates temporal constraints to enhance the fitting outcome of 3D-AAM. Assen et al. [123] replace the absolute intensity 3D-ASM with relative grey scales when ROI is being identified (fuzzy inference). Suinesiaputra et al. [124] propose to employ independent component analysis (ICA) [125] instead of PCA in SSM to extract myocardial contraction from SAX slices. Furthermore, due to the better performance on local description, ICA is used to design a classifier able to detect regional wall motion abnormalities. Fig. 11 Examples of detected LV myocardial strains visualised in 3D: a ED strain; b ES radial strain; c ES circumferential strain; d ES longitudinal strain [115] Fig. 12 A 3D-ASM (SPASM) LV segmentation technique [120] using GIMIAS platform: Step 1 user specifies three landmarks (the aorta, the mitral valve, and the apex) by three clicks on the cine MR volumes; Step 2 the platform automatically generates a model (a triangular surface mesh), which is pre-constructed in training stage, based on the three given landmarks; Step 3 the model fits to the target (feature point detected via fuzzy inference) through propagating the updates from the vertices close to the intersections between the surface and the image planes to distant regions on the earth Lekadir et al. [126] have also assessed myocardial motion through decomposing the global ventricular shape. They calculate the relationships between a series of spatiotemporal inter-landmarks. By tracking the epicardium and endocardium a dysfunction map is drawn to show abnormal contractions. O'Brien et al. [127] model shape, spatial, and temporal variation separately. They use a global contour optimisation instead of conventional ASM fitting. Roohi and Zoroofi [128] propose a kernel PCA (KPCA), in which the modes applied to represent a global ventricle shape are combined non-linearly. The distribution of landmarks is divided into intra-and inter-subspaces. A more recent work proposes to collect all the shapes learned from training data to build a dictionary [129]. The features of segmented frames from the test image are also added to the dictionary to create a patient-specific model dynamically. Each feature is classified into object (myocardial boundaries) or background (blood pool or muscles). A sparse shape model is then used to find the points on the ventricle walls based on their distances to the classified features. Unreliable points are abandoned and the complete LV shape is reconstructed according to the dictionary. Temporal constraints are not considered in this approach as current segmentation relies on the outcomes of the previous frame. Zhu et al. [130] developed a subject specific dynamical model that simultaneously handles inter-and intra-subject variabilities in a recursive Bayesian framework and a combined multi-linear PCA-ICA model. Starting from a manually segmented first frame, subsequent frames are segmented according to the current intensity pattern and a shape prior, predicted from the past frames.
Xenia et al. [131] proposed a framework for LV segmentation that is based on heuristic rules such as the brightness of the blood pool, sphericity of LV, and inter-slice smoothness of segmentation. A graph-cut algorithm was presented to infer the labels of the myocardium for robust optimization; however, the three-dimensional morphology of the heart was not fully exploited. A two-dimensional segmentation framework was proposed in [132] that effectively maps the edge patterns of each slice (centred at LV centroid) from polar into a Cartesian grid. Then a dynamic programming method is used to walk through the grid having the strongest edge values. Ayed et al. [133] propose another 2D segmentation framework that firstly learns intensity and the shape distributions for the blood pool from a manually segmented frame. Then, using a max-flow algorithm, it minimises a lower bound on the Battacharyya distance between the trained and subsequent distributions obtained from other frames. Similarly, Nambakhsh et al. [134] consider learning both intensity and shape constraints from a segmented first frame, and then minimise the distance between the test and trained distributions. However, rather than a graph-cut based method, a series of convex cost functions are solved for exact minimization. Yet, another training based LV segmentation method is proposed by Eslami et al. [135]; the prior information is implemented through kernel based approach, where test image is compared to the training data for the closest neighbour using a random walk paradigm. The method is shown to segment pathological hearts, whose data is usually overlooked using conventional PCA based statistical shape models.

Direct estimation
Afshin et al. [14] propose a direct estimation on heart abnormality detection. In their work, each subject comprises three SAX slices, whereas the apical slice is divided into four segments and the mid-ventricular and basal slices are divided into six segments. This user-provided segmentation is only needed in a single frame, which acts as a reference. All other subsequent frames in a complete cycle are automatically divided into 16 segments each, according to their distribution similarity. The local statistical descriptors, whose dimensions are reduced through linear discriminant analysis (LDA), are then constructed based on these segments. Because each segment can reflect the portion of blood filling, their statistics correlate well with the regional LV function. As a result, with a linear SVM trained by the ground-truth given by the radiologists and features from these local LDAs, regional abnormalities can be detected. They have classified 58 subjects (21 normal and 37 abnormal) with an accuracy of 86.09 %. Direct estimation for cardiac functional analysis can be a promising research direction and has achieved competitive performances in comparison with the state-of-the-art segmentation based assessment, with a significantly lower computational complexity.

RV
In the literature, the RV has received far less study compared to the LV, due to its more complex shape (in particular the variation of the complex crescent shape from the base to the apex), its thinner walls, and the similar intensity appearance with the trabeculations [11]. Because of these complexities, quantifying the regional myocardium wall thickness is not generally recommended for the MRIbased cardiac functional analysis of the RV. We discuss the available methodologies in the following and list them in Table 3.

Image-driven techniques
Maier et al. [136] segment the RV in MRI combining watershed filtering with graph-cut based region merging. They provide two initialisation options: the user either outlines the RV wall in 4-5 slices of the ED phase, or marks two points on the basal slice to register an atlas. Wang et al. [137] use a morphology-based algorithm, which considers the layout, shape, size, and relative locations to locate roughly the LV and RV first. The temporal discrepancy between two consecutive frames is then used to discover RV as the most active part. Ringenberg et al. [138] segment the endocardium by intersecting two ROI constrained binary images as follows. Firstly, an ROI window is selected and converted into binary with an optimal thresholding; next, the same ROI is convolved with a difference of Gaussian filter and thresholded at zero. The RV mask is roughly estimated as the intersection of these images. The window constraints label information from the previous slices and work as prior knowledge. The segmentation begins from the most basal slice in ED and ends at the apex. For ES, the prior is the union of labels from the previous slices at ES and the label of the current slice at ED. Punithakumar et al. [139] base their segmentation on registration and propagation. A 2D mesh delineating the endocardium or epicardium moves across all phases by establishing point-to-point correspondences. The manual segmentation of a single frame is required for initialisation. Mahapatra [140] uses a trained random forest classifier to give voxels two probability values, corresponding to the object and background. Based on this probabilistic map, a final segmentation is achieved by graph-cut. The image features they extracted for the discriminative description consist of intensity statistics, spatial context, textural, and curvature entropy. Nambakhsh et al. [141] propose a method based on the global shape and intensity similarity estimation. Based on the global distribution matching, the shape prior is intrinsically invariant with respect to translation and rotations. Centroid of LV and a small area of RV cavity have to be specified by the user. Compared to the learningbased approach in [140], this algorithm has the advantage of requiring only a single subject for training.

Model-driven techniques
Because of the geometrical complexity of the RV, robustness becomes a major concern in its automated segmentation. Amongst existing techniques, addressing this challenge, multi-atlas-based methods with label fusion [142][143][144] have received significant attention. In these frameworks, finding reliable correspondences between the patient and atlas spaces becomes critically important. For instance, Ou et al. [142] present a deformable registration algorithm that uses saliency of the matching for improved robustness versus variation of shape, intensity, and field of view. A "zoom-in" mechanism that uses the first round of RV segmentation to iteratively refine the registration and segmentation outcomes is employed. Alternatively, Grosgeorge et al. [145] employ a PCA-based SSM as the prior model, to guide the segmentation through a graph-cut method. The model is registered to the test case through a rigid transform, with two anatomical landmarks manually placed by the user on the ventricular septum. Oghli et al. [146] apply PCA on signed distance functions extracted from parametrised training contours as the shape prior for a deformable model. In addition, they use region and boundary based energies for improved fitting.

Bi-ventricle
Bi-ventricular segmentation uses slices covering both ventricles, from apex to the ventricular base, which is the valve plane, for full delineation of LV/RV myocardium. This research area has also been actively explored in the past 15 years (see Table 4). An overview of these methods is covered in this section.

Image-driven techniques
Sermesant et al. [147] have proposed a deformable biomechanical model based on tetrahedral geometric representation. The user specifies a proper mesh size to keep the data amount reasonable and retain a good mesh quality. The mesh is mapped to the test image using a non-rigid registration under influence of the internal and external forces, modelling elastic and imaging constraints, respectively. This method can be used for motion tracking, thus measuring local deformation and LVS is made feasible. Rougon et al. [148] employ a non-rigid registration method to assess myocardial contraction in both SAX and LAX slices. They use tagged MRI to infer the intra-myocardial motion and the cine MRI to extract the myocardial anatomy dynamically. Hautvast et al. [149] suggest a contour propagation scheme from ED to ES images. This method can be applied to SAX, two-chamber or four-chamber LAX slices, but requires manual segmentation on ED for initialisation. Cocosco et al. [150] convert a test image into a binary representation by optimally thresholding the intensity histogram of the ROI. The fat around ventricles is then removed by a thinning operation. All connected components are labelled and region-growing is performed on the SAX slices. Afterwards, they calculate the volume of each component at all frames and the maximum and minimum values are taken out. Two components having the most significant differences between their maximum and minimum volumes are selected as the LV and RV. The final delineation is obtained by merging the voxels classified in the first step along the LAX direction. Grosgeorge et al. [151] use the seminal model of active contours without edges [152] for bi-ventricular segmentation of a large dataset containing 1920 MR images, and obtained satisfactory results

3
comparable to the state-of-the-art. Mahapatra et al. [153] segment bi-ventricles using a graph-cut framework, guided by a shape prior based on the distribution of orientation angles from each pixel to the edge points, as extracted from a single manually annotated image. Wang et al. [154] adaptively use reinforcement learning to assimilate the knowledge provided by the user, such as edge point position correction, in LV/RV segmentation.

Model-driven techniques
Valdes et al. [155] propose to use a probabilistic atlas to guide the EM classification. The atlas provides a spatially and temporally varying probabilistic map for the LV, RV, myocardium, and background including the liver, stomach, lungs, and skin. The results of estimated volume of LV, RV, and myocardium demonstrate that the combination of the EM algorithm and a cardiac atlas improves segmentation accuracy. Bai et al. [156] fuse patched-based labels for a Bayesian formulation within a multi-atlas registration based segmentation framework. Furthermore, they refine the registration using intermediate label information. Figure 13 illustrates the procedure of multi-atlas label fusion and image registration refinement.
Ordas et al. [157] introduce a feature vector, which is invariant under Euclidean transforms in an ASM-based framework. Mitchell et al. [158] propose a hybrid AAM matching mechanism accomplished through three steps. Firstly, AAM alone is fitted to the image. Next, the hybrid AAM/ASM helps avoiding local minima by deploying the shape information. Finally, AAM is reapplied. Zhang et al. [159] also use a combined AAM-ASM model, which is based on novel spatial and temporal features, incorporating the motion. The combination of ASM and AAM yields better segmentation results and overcomes the drawbacks of using ASM or AAM individually. ASM requires good initialisation and can be trapped by incorrect nearby features, though it retains a fine global shape. While AAM performs well in tracking objects, but is easy to be trapped by local minima. Alba et al. [160] segment the LV and RV of highly abnormal hearts by using estimating a mapping between the abnormal image and the space of generic shape model built from a normal population, which can be thus used to segment any types of cardiac abnormality. Increased accuracy is demonstrated for both pulmonary hypertension and hypertrophic hearts.

Direct estimation
Wang et al. [15] estimate LV and RV cavity volumes on SAX slices without segmentation. This direct method relies on a likelihood function defined as the area correlation of the LV/RV cavities, and prior function specified by the product of the blobness, edgeness, and homogeneity. The framework consists of a training stage where the prior and likelihood probability functions are inferred. Given a test image, the posterior probability of observing a point in LV/ RV is derived using the Bayes rule. The mean cavity area of LV/RV is the expectation of a function of these posterior probabilities and the volumes are estimated using Simpson's method. However, as indicated by Zhen et al. [16], the limitations of [15] include a simple linear relationship assumption between LV and RV as well as an expensive computational requirement. Zhen et al. [16] make use of a three-layer convolutional deep network, which is learned from unlabelled images, to represent the input test case The label at a voxel (red dot) is given by the comparisons between the patch (yellow) on the target image and the patches (colourful boxes) on the warped atlases, weighed by the distance and similarity. Then the fusion of labels from all atlases assigns each voxel a final class. The segmentation result is used to refine the registration process [156] effectively in feature space. At the meantime, regression forests trained from manually labelled data, as discriminative learning, is responsible for estimating LV and RV volumes (see the flowchart in Fig. 14). They claim that their method significantly outperforms level-set and graph-cut methods. Another advantage of direct estimation is that the inconsistency of boundary and region intensity homogeneity has been excluded from the immediate influences on volumetric quantification.

LA
Segmentation of the LA is more challenging compared to the other structures in the heart. The shapes of LA may have different variations and its blood pool consists of other structures such as the auricular appendage and pulmonary veins; the surrounding pulmonary artery and the aorta have similar intensities to the atrium in MRI and the LA is typically much smaller than the ventricles, showing a relatively thin myocardium. The activity of mitral valve also makes the boundary between LA and LV invisible under some cases [12]. As a consequence, computer-aided LA segmentation has obtained much less progress. We list the methods in Table 5.

Image-driven techniques
John and Rahn [161] base their approach on thresholding and region merging. A thresholding roughly separates the blood pool voxels from the image. Afterwards, the Voronoi tessellation of the binarised mask is computed. The tessellated components are finally combined to segment the LA and other structures. Zhu et al. [162] also propose a region-growing framework, where the initial seed is found according to the anatomical knowledge from the middle SAX slice. A shape prior learned from training data is used to attract the growth to the statistically plausible region. This incorporation makes segmentation more robust to spatial variation and image quality.

Model-driven techniques
Karim et al. [163] construct a probabilistic atlas for atrium using 20 manually segmented training images. Given a test image, they apply an optimal thresholding to extract the blood pool and the vessel structures and obtain the Voronoi tessellation for the binarised image [161]. The narrow junctions, which are the connections between the atrium and its neighbouring structures, are then identified (as displayed in Fig. 15). Next, using the probabilistic atlas as prior, they present an MRF based cost function for segmenting cells that belong to the atrium. Additionally, a graph-cut method is applied for global optimisation. In order to deal with LA anatomical variations, Kutra et al. [164] have proposed a multi-component-based LA segmentation. The three most typical variations include the normal, common left trunk (CLT), and right middle pulmonary veins pattern (RMPV). Then a trained SVM is used to automatically select the model that fits the test image best. Eventually the model, which is a mesh of triangles, deforms towards the edge by the external and internal constraints.

Whole heart
The objective of whole heart segmentation includes delineation of LV, RV, LA, RA, and great vessels if required.  [16] Because of tissue diversity and indistinct boundaries between substructures; however, limited works have proved good efficacy in whole heart segmentation. The methods discussed in our review are summarised in Table 6.

Image-driven techniques
Makowski et al. [165] have proposed an active contourbased procedure to segment the heart and vessels in 2D transversal slices. This shape-independent method uses a balloon force to place the segmenting contour roughly and then uses a snake model to refine the segmentation. However, due to the complex geometry of the whole heart, later published works tend to use shape priors to increase robustness.

Model-driven techniques
Lotjonen et al. [166] reconstruct the 3D geometry of atria and ventricles from both SAX and LAX views. The pulmonary artery, pulmonary veins, and vena cava are excluded for volumetric measurement. The shape variability is modelled using PDM, a novel landmark distribution model, and a probabilistic atlas. Then the mean shape model is nonrigidly registered to the test image, and the model deforms towards the boundaries based on the shape priors. Since the performance of SSM-based methods is related to the richness of training samples, Koikkalainen et al. [167] have shown the feasibility to improve the segmentation of four-chamber and major vessels by artificially enlarging the training sets. Wierzbicki et al. [168] build PCA-based models for LV, RV plus RA, LA plus aorta, and the entire heart separately, using high quality training data. Each model is then registered to the mid-diastole frame of a low quality sequence, and propagated to all other frames by animating motion dynamics. Peters et al. [169] developed a deformable model by proposing a novel and robust boundary identifying technique called simulated search, whose mesh matching functions are previously trained. For the prior information-based approaches, model registration is always a critical step. Zhuang et al. [170] and Zuluaga et al. [171] find the breakthrough herein. They present a locally affine registration mechanism, which is further refined by a free-form deformation registration. This atlas-propagationbased method has turned out to be robust against various pathologies. Examples of the segmented whole heart in different views and a visualisation of segmentation errors between the result and the ground-truth are shown in Fig. 16.

Direct estimation
Zhen et al. [17] explored the feasibility in applying direct estimation to four-chamber volume measurement as well, by representing the MR images in a compact and discriminative way. The image features are generated using a supervised descriptor learning algorithm. Then the volume estimation becomes a multi-output regression problem solved with random forest.

Discussion
Despite the advances in cardiac image segmentation listed in this review, there are plenty of challenges waiting to be addressed to allow a more comprehensive assessment of cardiac function in clinical practice and medical research with MRI.

Choice of segmentation techniques
From this review, it can be seen there is a wide range of techniques and approaches that can be used for cardiac MR image segmentation. The choice of a particular technique is thus not trivial. However, a number of recommendations can be made. Firstly, the choice of the technique to be implemented can be constrained by the specific protocol. For example, a model-based technique can be used to obtain the walls of the LV is combined with thresholding to eliminate the effect of the papillary muscles. Secondly, Fig. 15 An LA blood pool (left) subdivided to Voronoi cells (middle). The narrow junction is the smaller sphere (right) locating between two larger components [163] the choice of a particular approach can depend on the availability of large training datasets. In such situations, modelbased approaches can be very powerful tools to restrict the segmentation results to valid instances. In the contrary, when only small cohorts are available for training, modelbased techniques can be too restrictive and methods that do not use any prior are preferred.
Finally, the obvious criterion for the choice of techniques should be the segmentation accuracy. However, while we provided a detailed list of the evaluations of the existing techniques in Tables 2, 3, 4, 5, and 6 for an overview of their performance, their direct comparison is difficult as the error metrics between the segmentation outcomes and the ground-truth are defined differently in different articles (point-to-surface errors, point-to-point errors, Hausdorff distance, dice similarity, correlation and linear regression coefficients, etc.). Also, the datasets are not the same in terms of image sequences, their numbers (sample size), and the classes (healthy vs. abnormal cases). For this reason, the emergence of challenges in international conferences is a very important initiative that will be able to highlight more objectively the merits and limitations of the existing methods. We can list, for example, the Left Ventricle Segmentation Challenge 1 (MICCAI 09), Right Ventricle Segmentation Challenge 2 (MICCAI 12), as well as Left Atrial Segmentation Challenge 3 (MICCAI 13).

Segmentation of the whole heart
Among the four chambers, the LV has received the most attention in cardiac segmentation and MRI-based cardiac functional assessment. This is because it plays a key role on the process of the blood circulation, and thus its function/dysfunction is associated with most cardiac diseases. Furthermore, the LV has a relatively simple geometry with thick myocardial walls, making its automated segmentation more feasible. In contrast, as it can be seen through comparing the list of works reviewed in this paper (Tables 2, 3, 4, 5, 6), the RV and LA have received less attention from the cardiac image analysis community (Fig. 17). This is due to the more complex geometry of these chambers and their much thinner walls. Yet, these chambers are associated with many critical diseases, such as modelling in patients with pulmonary hypertension [160] or left atrial enlargement [57][58][59]. Further research is thus required to develop techniques capable of coping with the difficulties of segmenting complex and thin structures such as the RV and RA, and more generally to segment the whole heart to enable an assessment that takes into account the combined motion of all chambers. Fig. 16 An evaluation of segmentation accuracy using surface-to-surface (S2S) distance between the segmented result and the manually delineated ground-truth from two different views in 3D [118]

Number of Publications
Training Free Training Required Fig. 17 The amount of referred publications in each section

Segmentation of large-scale CMR datasets
Another future perspective is related to the segmentation of large-scale datasets. In the era of big data, there is a demand for computational techniques that are scalable for the processing of thousands of cases and for the extraction of novel clinical knowledge from existing databases. However, previous cardiac analysis methods in MRI have been developed and validated with at most a few dozen cases, often with well-controlled imaging protocols, and based on a homogeneous class of subjects (e.g. healthy). A major research topic in the future will consist of extending the existing techniques such that they can handle the large variability in anatomy and MRI image sequence that are typically found in large-scale databases. Furthermore, current cardiac MR segmentation methods are rarely fully automatic. User interaction often used for example to define manually the apical and valve points. However, this becomes impossible when dealing with large numbers of datasets, and thus fully automatic techniques will be required.

Segmentation of abnormal cases
One major research challenge in cardiac segmentation is the development of approaches that are robust to different groups of individual and classes of disorders. In the existing literature, however, most techniques have been mainly developed and validated with normal subjects, and in some exceptional cases with mildly abnormal hearts, i.e. mostly few regional septal defects such as hypertrophic cardiomyopathy (HCM) [124,135]. These techniques are developed in a generic form for both normal and abnormal cases and do not have a mechanism to handle explicitly large remodelling effects owing to cardiac diseases. Recently, Alba et al. [160] developed a technique specifically designed to segment severely abnormal hearts, with a promising validation to pulmonary hypertension patients with highly remodelled RV. Such techniques need to be further investigated using large cohorts and with multiple diseases to make the tools more robust for clinical use, where routine cardiac MRI quantification is concerned mostly with diseased subjects or subjects suspected to be diseased.

Clinical translation
Finally, significant effort is being dedicated, in parallel to the consolidation of the cardiac image analysis techniques, to the clinical translations of software tools that can be used robustly and routinely in clinical practice. Table 7 presents some of the existing software used in clinical practice or in cardiovascular research in alphabetical order. We recommend the readers to check the details of the available software on their websites as the functionalities tend to evolve continuously over time as the result of new advances in CMR research.

Conclusions
This review paper has summarised the most recent advances in cardiac image segmentation methods, which can be employed for the assessment of cardiac structure and function with CMR. These approaches range from image classification based techniques to statistical shape models. We have highlighted the properties of each of these approaches and their links to cardiac structure and functional assessment in MRI. After years of continuous developments, cardiac segmentation has become an interdisciplinary subject associating cardiology, medical imaging, and image processing. Further research is required to consolidate these advances with validation to larger cohorts, as well as to extend these approaches to the segmentation of all chambers and pathological hearts, ultimately allow for a more comprehensive application of the existing tools in clinical practice.