BACKGROUND OF THE INVENTION

1. Field of the Invention

This invention relates broadly to cardiothoracic imaging. More particularly, this invention relates to segmenting cardiac chambers in a fourdimensional magnetic resonance image.

2. State of the Art

Accurate and reproducible assessment of the left and right ventricular volumes and ejection fraction is essential for the prognosis of patients with heart disease and for evaluating therapeutic response. Magnetic resonance (MR) imaging has proven to be an accurate and reproducible imaging technique for the quantitative analysis of left and right ventricular functions. Ventricular volumes and mass as well as regional functional parameters such as wall motion and wall thickening, especially for the left ventricle, can be obtained in the true shortaxis plane of the ventricle. These shortaxis views have various advantages. Since wall motion and wall thickening takes place mainly perpendicularly to the viewing direction, they can be followed best in this direction. Furthermore, partial volume effects are minimized. A major advantage of MR imaging is that it does not use ionizing radiation and is noninvasive since no contrast agent is used. Through using multiple imaged slices a combination of spatial and time coordinates is considered.

Chris A. Cocosco, in “Automatic cardiac regionofinterest computation in cine 3D structural MRI”, CARS and Elsevier, 2004, describes a method for extracting a region of interest around the heart in shortaxis MRI cardiac scans. The method produces a region within fourdimensional (4D) datasets wherein the left and right ventricles are localized. No attempt is made to differentiate between the left and the right ventricle. Neither are the base and apex locations of the ventricle(s) defined. The method may in particular fail whenever high local variances are present in the 4D dataset that may originate from main blood vessels (lung/aorta) or from specifically superimposed text such as trigger time. The method may also fail when the 4D image contains an inhomogeneity. This inhomogeneity can be caused by lower sensitivity at a larger distance between the heart and the imaging coils of the MR apparatus which then leads to lower pixel intensity in the image.
SUMMARY OF THE INVENTION

The present invention presents a more systematic approach, in that it does not rely on anatomy assumptions. First, it locates the left and right ventricles. Then the base short axis slice and the apex of the left ventricle shortaxis slice are located. Next, the start conditions for the segmentation process are determined. The locations of the left ventricle, the left ventricle epicard and finally the right ventricle endocard are detected approximately. Based on these start conditions the segmentation process starts. Note that various standard acquisitions of MR cardiac imaging are used, such as the short axis (SA) and long axis (LA) orientations. Finally the papillary muscles of the left ventricle are detected. Advantageously, the methodology of the present invention is fully automatic without user interactions based on intermediate results, which can yield 100% reproducibility.

Two procedures can be followed. In the first procedure, a 3D fuzzy extraction is followed by a fine tuning step to accurately locate the left ventricle endocard and epicard. In the second procedure, a new contour propagation algorithm locates the left ventricle endocard and epicard. Neither of the two above procedures relies on any anatomy assumption, which makes the segmenting also suitable to use with diseased ventricles that may be deformed to a certain extent.

The segmentation process is improved over the prior art segmentation techniques, which range from being purely data driven (using edge detection and region growing) on the one hand, to being model driven on the other hand.

One such prior art segmentation technique, the socalled “snake” model, presents a framework that is the basis of boundary driven image segmentation operations. The “snake” model is described by G. Kass in “Snakes: Active Contour Models”, Int. Journal of Computer Vision, 1, pp. 321332, 1988. Briefly, the snake model refers to an energy minimization technique that seeks for the lowest potential of a curvebased function. This function is a compromise between boundary terms and image terms. However, this segmentation method is quite sensitive to noise and physically artifacted image data. Various modifications thereto have been proposed, but the method remains sensitive to userdefined start conditions, which endanger the reproducibility of the segmentation.

Another example of a prior art segmentation technique is described in U.S. Pat. No. 5,239,591, “Contour extraction in multiphase, multislice cardiac MRI studies by propagation of seed contours between images”, which is incorporated by reference herein in its entirety.

The segmentation process of the present invention corrects for data acquisition and/or subjectcaused movement in the acquired MR image data set by using system coordinates from an MR scanner apparatus. In this way, accuracy is improved without disturbing the sequence of image pickup operations. Moreover, the segmentation of the endocardial and epicardial features starts from an operatordefined seed contour and is based on a radial minimum cost with dynamic sign definition and tissue classification. Apart from such elementary user input, all further operations are machinecontrolled, which yields excellent reproducibility.

In addition, the methodology of the present invention automatically detects Left Ventricle Right Ventricle (LVRV) connection position. Clinicians are recognizing this information as highly helpful.

Moreover, the methodology of the present invention estimates volume(s) of the heart which is/are derived from the spatial location of the heart's apex position on a curve that is obtained from the long axis cine sequences. In addition, such volume estimates preferably are derived from a mitral valve position obtained from long axis cine views. Such processing improves the accuracy of such volume estimate(s).

Additional objects and advantages of the invention will become apparent to those skilled in the art upon reference to the detailed description taken in conjunction with the provided figures.
BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flow chart illustrating various processing steps of the methodology of the present invention;

FIG. 2 is a schematic representation of inplane movement which is identified and corrected by the methodology of the present invention;

FIG. 3 is a threedimensional visualization of breathing movement of a subject person;

FIG. 4 is a flow chart illustrating the operations in locating the heart's left and right ventricle in accordance with the present invention;

FIG. 5 is an exemplary pixel intensity histogram of the heart's left ventricle;

FIG. 6 illustrates a mapping procedure from longaxes (LA) slices to shortaxis (SA) slices or inversely;

FIG. 7 is a flow chart illustrating segmentation operations that start from usersupplied longaxis (LA) input;

FIG. 8 is a highlevel schematic diagram of an exemplary MRI apparatus.
DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

Turning now to FIG. 1, there is shown a flow chart of various processing steps in accordance with the present invention. First, in block 10 (localize Apex/Base of left ventricle, estimate epicard, find end of diastole, end of systole), the 4D MRI dataset is processed to localize the left and right ventricles therein. This dataset contains a sequence of short axis cine images that cover the heart, combined with twochamber and fourchamber cine sequences.

As an alternative to the above, it is possible to work with shortaxis images only. Within this operation, the base and apex short axis slices are located as well as the end systole and end diastole heart phase in time. In this respect, FIG. 4 illustrates a flow diagram of the locating process for the heart's left and right ventricles, as a substitute for block 10 in FIG. 1. The remainder of FIG. 1 will be discussed later.

Turning to FIG. 4, in block 40, the images are preprocessed to eliminate high local variance and inhomogeneity according to the expression (1) below. Therein, I_{p }is the 2D image array that is subjected to conventional Fourier transform methods to render the result more reliable.
$\begin{array}{cc}{I}_{p}={\mathcal{J}}^{1}\left\{\mathcal{J}\left(I\text{\hspace{1em}}\Theta \text{\hspace{1em}}b\right)\xb7H\right\}\xb7T\left({P}_{\mathrm{shift}}\right)\text{}\mathrm{where}\text{\hspace{1em}}H=1{e}^{\frac{{r}^{2}}{2\xb7{\partial}^{2}}}\text{\hspace{1em}}\mathrm{in}\text{\hspace{1em}}\mathrm{case}\text{\hspace{1em}}\mathrm{of}\text{\hspace{1em}}\mathrm{inhomogeneity}\text{}\text{\hspace{1em}}H=1\text{\hspace{1em}}\mathrm{otherwise}\text{}\text{\hspace{1em}}b\text{\hspace{1em}}\mathrm{is}\text{\hspace{1em}}a\text{\hspace{1em}}\mathrm{square}\text{\hspace{1em}}\mathrm{structured}\text{\hspace{1em}}\mathrm{element}\text{}\text{\hspace{1em}}{P}_{\mathrm{shift}}\text{\hspace{1em}}\mathrm{is}\text{\hspace{1em}}\mathrm{the}\text{\hspace{1em}}\mathrm{displacing}\text{\hspace{1em}}\mathrm{caused}\text{\hspace{1em}}\mathrm{by}\text{\hspace{1em}}\mathrm{acquisition}& \left(1\right)\end{array}$

In the preprocessing operations of block 40, grey level erosion or “thinning” is effected, after which a filter eliminates low frequencies caused by inhomogeneity. The inhomogeneity is detected by resampling along a scan line that runs from the centre of the two chamber and four chamber epipolar intersection lines, called l_{2ch }and l_{4ch }respectively, in the first frame of the middle shortaxis slice. This is followed by executing a first Gaussian derivative operation.

From the centre, the positions of the steepest absolute slopes l_{p }and r_{p }to the left and right are located. Between those positions a linear fit is performed based on the resampled image data. Whenever the absolute slope of this linear fit exceeds a certain threshold it is assumed that the short axis image data set contains an inhomogeneity.

Moreover, in the preprocessing operations of block 40, the displacement caused by the data acquisition is derived and corrected for. For acquiring the shortaxis data set typically 8 mm thick cine shortaxis slices with breathholding by the patient will encompass the entire left ventricle from apex to base. During this process the MR acquisition parameters for the acquired spatial slice location could have changed such that they do not fit perfectly in the acquired bounding box. By using the absolute system coordinates obtained from the MR system (through the “headers” of the conventional DICOM file format), it is possible to correct for these ‘inplane’ movements.

In this respect, FIG. 2 illustrates a spatial configuration of this inplane movement. In this FIG. 2, item 100 is the bounding box, two successive planes have been labelled 102 and 104, respectively, and the plane wherein the actual inplane movement occurs has been labelled 106. The line 108 should run continuously along the respective first elements that have been shown as perspective squares (110), of which the lowest one as shown has been shifted to the front of the bounding box.

Turning back to FIG. 4, the operations continue to block 41 to correct for such “inplane” movement. Such operations begin by identifying the starting plane (102) and then computes the threedimensional (3D) line l that runs perpendicularly through this starting plane (102) at the first image element (pixel/voxel) (110) of the image frame identified by starting plane (102). Next, for each one of a number of successor planes, the displacement between the first element of the successor plane under investigation and the 3D line/is computed. In pseudo code:
V _{1} =P _{1} +μR _{1} +γ·C _{1 }
l=P _{1}+(μ×·C _{1})·α (2)

Next, the intersection point P_{3D }of the line l with the successor planes is defined in the 3D world coordinate system:
P _{3D} =l∩V _{n }{n>1} (3)

Finally, the intersection point P_{3D }is transformed into the twodimensional (2D) coordinate system of the successor plane:
$\begin{array}{cc}{P}_{2D}=\left[\begin{array}{c}1\\ 1\\ 0\end{array}\right]\xb7S\left(1/{\mathrm{voxelsize}}_{n}\right)\xb7{R}_{n}\xb7T\left({P}_{n}\right)\xb7{P}_{3D}& \left(4\right)\end{array}$
Here, the rotation matrix R is defined as the cross product of the plane direction vectors:
$\begin{array}{cc}R=\left[\begin{array}{ccc}{R}_{1x}& {C}_{1x}& {N}_{1x}\\ {R}_{1y}& {C}_{1y}& {N}_{1y}\\ {R}_{1z}& {C}_{1z}& {N}_{1z}\end{array}\right]\text{}N=R\otimes C& \left(5\right)\end{array}$

For a particular plane (e.g., slice), the displacement is defined as P_{2D}. All frames of the current plane will be shifted according to this vector P_{2D}. The step is repeated for the remaining planes along the shortaxis stack.

In addition to such ‘inplane’ movement, the patient can move between the acquired shortaxis cine sequences, the twochamber and the fourchamber cine sequence. This moving is specific for each examination and its magnitude depends mainly on how well the patient reproduces each breathholding in an MR cine sequence. Typically, the motion ranges from zero to 10 mm. In block 42 a a method is presented giving the physician the feasibility to evaluate in a quick and easy to use environment the quality of the MR acquisition in respect to slice movement. Thus, the physician can decide if some cine sequences should be taken once more while the patient is still available for MR acquisition. This improves the acquired dataset for further therapeutic usage without requiring an expensive patient recall for another MR.

The method in block 42 a allows the physician to inspect the acquired MR dataset is by viewing the images from each acquired plane (used for treatment decisions) in a 3D environment conform the absolute system coordinates from the MR system. Generally, each image contains visible landmarks, such as organs and their details. Interpretation by the physician, e.g. the left ventricle, right ventricle, and viewing each image plane in a 3D environment spatially conform the 3D coordinates allows to quickly localize any movement visually.

FIG. 3 illustrates a threedimensional visualization of breathing movement of a subject. Rotating this 3D visualisation renders evaluating the quality regarding movements in the MR image dataset straightforward. In FIG. 3, item 130 presents a single short axis slice, which has the shape of a closed form. Item 132 shows a twochamber cine sequence that has the shape of a curve that is closed at one end, the apex, and open at the other. Combining these two without intermediate movement at item 134 shows the two curves to intersect. With breathing (or other patient caused movement) at item 136 show the two curves to be relatively displaced. This may lead to decide an image retake for correction in block 42 a. Even better visualisation is acquired by image enhancing on the image dataset of each plane.

Automatically correction for breathing movements is achieved in block 42 b by image registration techniques. In such techniques, a twochamber image is created by resampling each shortaxis image along the intersection line (epipolar line) of the twochamber image, which gives an image which can be used for registration against the true twochamber plane found during the acquisition. Breathing movement can now be computed by correlating in the Fourier domain each resampled line against the corresponding twochamber intersection. This method can be extended with the fourchamber view.

In block 43, after preprocessing all shortaxis frames, the total variance image of the complete shortaxis dataset is generated. Because the shortaxis images are perpendicular to the long axis of the ventricle (its zaxis corresponds to the voxel thickness) and the heart generally remains in motion, it is possible to effectively compress the 4D motion/variance into a single 2D image:
$\begin{array}{cc}\mathrm{VAR}=\frac{1}{\mathrm{n\_frames}}\sum _{i=1}^{\mathrm{n\_frames}}{\left({I}_{{p}_{i}}{\stackrel{\_}{I}}_{p}\right)}^{2}\text{}{I}_{\mathrm{variance}}=\mathrm{MAX}\left(\mathrm{VAR}\left(\mathrm{short}\text{\hspace{1em}}\mathrm{axis}\text{\hspace{1em}}{\mathrm{slice}}_{n}\right)\right)+\frac{\sum _{n=1}^{\mathrm{n\_slices}}(\mathrm{VAR}\left(\mathrm{short}\text{\hspace{1em}}\mathrm{axis}\text{\hspace{1em}}{\mathrm{slice}}_{n}\right)}{\mathrm{n\_slices}}& \left(6\right)\end{array}$
The mean intersection points derived from the twochamber and fourchamber epipolar lines projected on each shortaxis slice are smoothed and added to the overall variance image:
$\begin{array}{cc}\stackrel{\_}{p}=\frac{\sum _{1}^{\mathrm{n\_slices}}{l}_{2\mathrm{ch}}\bigcap {l}_{4\mathrm{ch}}}{\mathrm{n\_slices}}\text{}{l}_{2\mathrm{ch}}={V}_{{\mathrm{SA}}_{i}}\u22d0{V}_{2\mathrm{ch}}\text{\hspace{1em}}\mathrm{mapped}\text{\hspace{1em}}\mathrm{on}\text{\hspace{1em}}{\mathrm{SA}}_{i}\text{}{l}_{4\mathrm{ch}}={V}_{{\mathrm{SA}}_{i}}\u22d0{V}_{4\mathrm{ch}}\text{\hspace{1em}}\mathrm{mapped}\text{\hspace{1em}}\mathrm{on}\text{\hspace{1em}}{\mathrm{SA}}_{i}\text{}{\partial}_{p}=\mathrm{STDEV}\left(p\right)& \left(7\right)\end{array}$
The mapping back on the SA slice is performed by:
$\begin{array}{cc}{P}_{2{D}_{\mathrm{SA}}}=\left[\begin{array}{c}1\\ 1\\ 0\end{array}\right]\xb7S\left(1/{\mathrm{voxelsize}}_{\mathrm{SA}}\right)\xb7{R}_{\mathrm{SA}}\xb7T\left({P}_{\mathrm{SA}}\right)\xb7{P}_{3D}& \left(8\right)\end{array}$
The variance image is now adjusted with this information according to:
$\begin{array}{cc}{I}_{\mathrm{totalvariance}}={I}_{\mathrm{variance}}\xb7\sum _{i={p}_{y}{\partial}_{{p}_{y}}}^{{p}_{y}+{\partial}_{{p}_{y}}}\sum _{{p}_{x}{\partial}_{{p}_{x}}}^{{p}_{x}+{\partial}_{{p}_{x}}}{I}_{{\mathrm{variance}}_{\left(i,j\right)}}{e}^{\left(\frac{\left(x{\stackrel{\_}{p}}_{x}\right)}{2\xb7{\partial}_{{p}_{x}}}+\frac{\left(y{\stackrel{\_}{p}}_{x}\right)}{2\xb7{\partial}_{\mathrm{py}}}\right)}& \left(9\right)\end{array}$

Next in block 44 (threshold, dilate, convex hull), the generated variance image I_{total variance }is smoothed with a Gaussian. Also, a dual threshold is defined. Preferably, the dual threshold is defined according to the Otsu threshold method described in N. Otsu “A threshold selection method for greylevel histograms” IEEE Trans. Syst. Man Cybernet SMC01 (1970) pp. 6266, herein incorporated by reference in its entirety.

After threshold selection, the occurrence of the detected objects is defined. When the first object covers an area of at least 85% of the histogram (FIG. 5), the threshold is based on the single method. Otherwise the threshold is selected between second and third objects derived from the histogram of I_{total variance}.

After thresholding/binarization, all holes in the binary image are filled in as based on a morphology operation. In a labelling operation the object with the largest area is extracted. The grey level erosion is undone again by a binary dilation with the same kernel as the earlier grey level erosion. Then, a binary dilation has a kernel size of approximately the myocard wall thickness. Finally a convex hull is fitted which localizes the complete left and right ventricles in the 4D shortaxis dataset, defined as mask M. The shortaxis slices are perpendicular to the long axis of the left ventricle, which allows to approximate the left ventricle by a circle in these shortaxis slices.

In block 45, an adaptive Hough transform combined with the generated mask of the heart area, enables automatic location of the left ventricle. The Hough transform s wellknown to persons skilled in the art for detecting features that can be described as a parametric equation. For a circle with a and b the coordinates of the centre and radius r, this is:
x(t)=r·cos(t)+a
y(t)=r·sin(t)+b (10)

The localization of the left ventricle is done on the middle shortaxis slice. Within this slice the most likely enddiastole frame is found by selecting the frame which has the most elements within the mask M and above the threshold as defined during the localization phase of the heart.

Within this frame, the edges are located preferably by a socalled Canny edge detector, which is wellknown to persons skilled in the art. First, the frame is smoothed with a Gaussian filter to reduce noise and unwanted details. Then, magnitude and direction of the gradient are calculated.
$\begin{array}{cc}\nabla I=\sqrt{\uf603{G}_{x}^{2}+{G}_{Y}^{2}\uf604}\bigcap M\text{}\theta ={\mathrm{tan}}^{1}\left(\mathrm{Gy}/{G}_{x}\right)\bigcap M& \left(11\right)\end{array}$
The kernels used to derive G_{x }and G_{y }are:
$\begin{array}{cc}{G}_{X}=\left[\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right],{G}_{y}=\left[\begin{array}{ccc}1& 0& 1\\ 2& 0& 2\\ 1& 0& 1\end{array}\right]& \left(12\right)\end{array}$

The gradient magnitude image ∇I is binarized to I_{bin }using a predefined threshold. This gives a feature map that contains the edges and background. The edges are likely located on the perimeter of the circle that must be detected. Therefore, every detected edge point at (x,y) of I_{bin }is parsed and suggested to be on a circle boundary with radius r between r_{min }and r_{max}. The position of the centre (C_{x},C_{y}) is then predicted from the edge location (x,y), the gradient direction θ and its radius r.
C _{x} =x−r·cos(θ)
C _{y} =y−r·sin(θ) (13)

 Where rε[r_{min},r_{max}.]
All found elements (C_{x},C_{y}) are accumulated into the Hough image space H which is processed with a second order Gaussian derivative within the mask M.
$\begin{array}{cc}{\nabla}^{2}H=H\xb7\left(\frac{\left({x}^{2}+{y}^{2}\right){\partial}^{2}}{{\sigma}^{4}}\right)\xb7{e}^{\left(\frac{\left({x}^{2}+{y}^{2}\right)}{2{\sigma}^{2}}\right)}\bigcap M& \left(14\right)\end{array}$

Now the position with the highest value in ∇^{2}H is selected as the radius of the circle (C_{x},C_{y}), thus yielding the centre point of the left ventricle.

The frame I_{p }is binarized preferably with the Otsu threshold method within the mask M, so the threshold is based on the histogram:
H _{p} _{ m } =HIST(I _{p} ∩M) (15)

Because the threshold is based on the frame histogram Hp_{m }it could happen that next to blood tissue and myocardium tissue, also background were represented in the histogram Hp_{m}. Therefore, a second threshold is defined based on the histogram within the circle defined by its centre (C_{x},C_{y}) and radius r_{min}. The threshold preferably used during the binarizing process is:
$\begin{array}{cc}\begin{array}{c}\begin{array}{c}\begin{array}{c}\mathrm{Single}\text{\hspace{1em}}\mathrm{Otsu}\text{\hspace{1em}}\mathrm{on}\text{\hspace{1em}}{H}_{{p}_{M}}\Rightarrow {T}_{M};\left({\mu}_{{M}_{0}},{\sigma}_{{M}_{0}}\right);\left({\mu}_{{M}_{1}},{\sigma}_{{M}_{1}}\right)\\ \mathrm{Dual}\text{\hspace{1em}}\mathrm{Otsu}\text{\hspace{1em}}\mathrm{on}\text{\hspace{1em}}{H}_{{p}_{M}}\Rightarrow {T}_{{\mathrm{Md}}_{1}};{T}_{{\mathrm{Md}}_{2}};\left({\mu}_{{\mathrm{Md}}_{0}},{\sigma}_{{\mathrm{Md}}_{0}}\right);\left({\mu}_{{\mathrm{Md}}_{1}},{\sigma}_{{\mathrm{Md}}_{1}}\right);\end{array}\\ ({\mu}_{{\mathrm{Md}}_{2}},{\sigma}_{{\mathrm{Md}}_{2})}\end{array}\\ \mathrm{Single}\text{\hspace{1em}}\mathrm{Otsu}\text{\hspace{1em}}\mathrm{on}\text{\hspace{1em}}{H}_{{p}_{\mathrm{circle}}}\Rightarrow {T}_{C};\left({\mu}_{{C}_{0}},{\sigma}_{{C}_{0}}\right);\left({\mu}_{{C}_{1}},{\sigma}_{{C}_{1}}\right)\\ T=\{\begin{array}{cc}{\mu}_{{C}_{1}}{\sigma}_{{C}_{1}}& \mathrm{if}\text{\hspace{1em}}\left({T}_{M}>{\mu}_{{C}_{1}}\right)^\left({T}_{{\mathrm{Md}}_{1}}>{\mu}_{{C}_{1}}\right)\\ {T}_{{\mathrm{Md}}_{1}}& \mathrm{if}\text{\hspace{1em}}\left({T}_{M}>{\mu}_{{C}_{1}}\right)^\left({T}_{{\mathrm{Md}}_{1}}\le {\mu}_{{C}_{1}}\right)\\ {T}_{M}& \mathrm{otherwise}\end{array}\end{array}& \left(16\right)\end{array}$

Within this binarized image the object is extracted in which the centre (C_{x},C_{y}) and the mean position derived from the two and fourchamber intersection lines p are located. The left ventricle may contain papillary muscles, so the centre (C_{x},C_{y}) or p need not lie within the extracted object. In such case the process is repeated for the next slice, above and/or below.

After successfully locating the left ventricle according to block 45, the enddiastole and endsystole phases are automatically determined in block 46 (estimate ED/ES heart phase), by:
ED=arg max(Area_{i})
ES=arg min(Area_{i}) (17)
Herein, Area_{i }is the area of the object extracted with the threshold T obtained through applying expression 16 on the frame i from the current shortaxis slice. Again the binarization is performed within the mask M.

Finally, in block 47 the base shortaxis slice and apex shortaxis slice are defined at the enddiastolic phase as computed by expression 17. For each frame I_{p }within the enddiastolic phase the frame is binarized according to expression 16 which gives B_{i}. The process starts at the “seed” frame in which the left ventricle centre position has been determined previously, looking at the frame above and below.

Within each frame the following steps are performed, looking once at the frames spatially above and once looking at the frames spatially below.

 Perform thresholding on frame I_{p }according to expression 16 that gives B_{i}.
 Intersect with mask M
 Fill holes located in B_{i }by using a morphology flood fill algorithm.
 Extract the object that obtains the centre position C_{p }The extracted object gives a global estimation of the endocardial contour location B_{endo} _{ i }
 Define Area_{i }and perimeter_{i }of extracted object B_{endo} _{ i }
 Compute the circularity (C) of the extracted object by:
$\begin{array}{cc}C=\frac{4\xb7\pi \xb7{\mathrm{Area}}_{i}}{{\mathrm{Perimeter}}_{i}^{2}}& \left(18\right)\end{array}$
 If the circularity is below a predefined threshold, or if the computed area (Area_{i}) is three times above the area of the frame in which the left ventricle position is defined, the process stops and a ventricle boundary is defined as Ventricle boundary[ ]=i
 Centre position C_{p }is adjusted according to:
$\begin{array}{cc}{C}_{p}\Rightarrow \{\begin{array}{c}x=\frac{\sum _{\mathrm{pixels}}x}{{\mathrm{Area}}_{i}}\\ y=\frac{\sum _{\mathrm{pixels}}y}{{\mathrm{Area}}_{i}}\end{array}& \left(19\right)\end{array}$
 The global estimation of the epicardial border is defined B_{epi} _{ i }

The algorithm for segmenting the epicardial border is based on a radial minimum cost algorithm and uses as starting entity the estimated boundary of the endocardial B_{endo} _{ i }. A detailed description will be given herebelow when describing the segmentation of the left ventricle epicardial.

The mean area is found from the frames above and below the seed frame. Therewith the apex and base slice are defined as:
$\begin{array}{cc}\begin{array}{c}\stackrel{\_}{A}\mathrm{rea}\text{\hspace{1em}}1=\frac{\sum _{i={\mathrm{seedframeb}}_{1}}^{{\mathrm{Vb}}_{1}}{\mathrm{Area}}_{i}}{{\mathrm{Vb}}_{1}\mathrm{seedframe}}\\ \stackrel{\_}{A}\mathrm{rea}\text{\hspace{1em}}2=\frac{\sum _{i={\mathrm{Vb}}_{2}}^{\mathrm{seedframe}}{\mathrm{Area}}_{i}}{\mathrm{seedframe}{\mathrm{Vb}}_{2}}\\ \mathrm{apex\_slice}=\{\begin{array}{cc}{\mathrm{Vb}}_{1}& \mathrm{if}\text{\hspace{1em}}\stackrel{\_}{A}\mathrm{rea}\text{\hspace{1em}}1<\stackrel{\_}{A}\mathrm{rea}\text{\hspace{1em}}2\\ {\mathrm{Vb}}_{2}& \mathrm{otherwise}\end{array}\\ \mathrm{base\_slice}=\{\begin{array}{cc}{\mathrm{Vb}}_{1}& \mathrm{if}\text{\hspace{1em}}\stackrel{\_}{A}\mathrm{rea}\text{\hspace{1em}}1\ge \stackrel{\_}{A}\mathrm{rea}\text{\hspace{1em}}2\\ {\mathrm{Vb}}_{2}& \mathrm{otherwise}\end{array}\end{array}& \left(20\right)\end{array}$
The right ventricle is localize by subtracting the global estimation of the epicardial border in each frame (B_{epi} _{ i }) from the thresholded image I_{p }according to expression 16 masked with M, this yields to a global estimation of the right ventricles endocardial border, defined as RV_{endo} _{ i }. This ends the discussion of FIG. 4.

Now, turning back to FIG. 1, in block 10, the location of the left ventricle is identified. Thus, the indicating in the 4D stack of the left/right ventricle position, the enddiastolic phase, endsystolic phase, apex slice, base slice and the global boundary positions of the left ventricle epicardial have been done.

In block 11, the segmentation of the endocardial boundary is performed in a complete slice, which gives a spatiotemporal 3D image (i.e. 2D plus time) within the area eroded n times around the estimated epicardial boundary and masked with M when segmenting the left ventricle. During the segmentation of the right ventricle, the estimated endocardial boundary masked with M is used. During the erosion a discstructured element is used with radius one. For the segmentation a fuzzy connectedness algorithm can be used. The algorithm is preferably based upon on the wellknown concept of fuzzy objects in which image elements (pixels/voxels) exhibit a similarity or “hanging togetherness” both in geometry and in grey scale. See J. K. Udupa, S. Samarasekera, “FuzzyConnectedness and Object Definition: Theory, Algorithms and Applications in Image Segmentation” Graphic Models, Image Processing, vol. 58, 246261, 1996, herein incorporated by reference in its entirety.

Pixels/voxels relate to each other by their adjacency and similarity. For simplicity, adjacency is defined as two neighbouring pixels in a (3D equivalent of) 4connectivity sense. Therefore, only left/right, top/bottom and front/back neighbours are considered to have a relation. This relation is expressed by their similarity with respect to pixel grey values (also called affinity). It depends on the application which relation is used for finding similarity. Among various standard relations, in the preferred embodiment affinity μ is expressed as a weighted sum of a Gaussian and an additive Gaussian. The additive Gaussian is based on the intensity (grey value) P and the Gaussian on the intensity difference (gradient) G.
μ(P _{C} ,P _{D})=F(P)·N(G)
with
P=(P _{c} +P _{D})/2
G=P _{c} −P _{D} (21)

 in which P_{C }and P_{D }are the intensities of two neighbour pixels. The function N is a Gaussian or standard Normal function, specified by the mean μ, and standard deviation σ.

The function F is an adaptive Gaussian: the pixels inside the endocardium have an intensity distribution given by the right peak in the image histogram in FIG. 5. First, we could take this distribution modelled as a Gaussian N_{2}. However, additional knowledge is present: pixels outside the endocardium may not be included in the object to be extracted. Therefore, we have to penalise these pixels by adding a second term in F:
$\begin{array}{cc}\begin{array}{c}F=S\xb7\left(1{N}_{1}\right)\\ S=\{\begin{array}{cc}{N}_{2}& \mathrm{if}\text{\hspace{1em}}P\le {\mu}_{2}\\ {N}_{2}& \mathrm{with}\text{\hspace{1em}}\sigma =3{\sigma}_{2}\text{\hspace{1em}}\mathrm{if}\text{\hspace{1em}}P>{\mu}_{2}\end{array}\end{array}& \left(22\right)\end{array}$

In this respect, FIG. 5 illustrates an exemplary pixel intensity histogram of the heart's left ventricle. In particular, S generally equals a Gaussian distribution function N_{2}, but for intensities above the mean μ_{2 }it is broadened by taking a standard deviation that is 3 times the true value. This pseudo Gaussian will cause pixels with large intensities to also be extracted when they would otherwise be excluded. This takes the fact into account that one or two frames sometimes show hyper intensities within the endocardium occur that are caused by blood flow which leads to locally high concentrations of protons. The second term in expression (22) is the penalty term: grey values below the threshold (dashed line in FIG. 5 separating the myocardium from the endocardium) will have an extra low affinity.

To start extracting the fuzzy object, a seed point is needed which should be located inside the expected object. Taking the symbol f_{o}(c) to represent the scene value at position c, the pseudo code for the algorithm follows:

 1. All scene values are set to zero and the scene value at the seed point is set to 1.
 2. Push each neighbour of the seed point to the queue if its affinity strength is greater than 0.
 3. Enter a while loop:
 a. Remove an element from the queue. Its position is denoted by c. If queue is empty then the extraction is completed.
 b. Find maximum of weakest affinity strength around this element. For each neighbour d the affinity strength is computed and the minimum (weakest) of this strength and the existing scene value f_{o}(d) at the position of the neighbour is taken. The maximum f_{max }of these values for all 6 neighbours is derived then.
 c. If the maximum value f_{max }is larger than the existing scene value f_{o}(c) then
 1. set scene value at c equal to the found maximum value, so f_{o}(c)=f_{max}.
 2. push all elements around c with affinity strength >0 to queue.
 d. endif
 2. endwhile

The output of the algorithm is a socalled “scene” giving at each position in the 3D image/stack the connectivity or measure of “hangingtogetherness” regarding the seed point. Also here, a segmentation step must assign a scene value to the object or not. The scene is threshold conform:
T _{scene}=scene_{i}∃_{i}{σ_{x}>min_{x}ˆσ_{y}>min_{y}ˆσ_{g}>min_{g} ˆi≧min_{n}} (23)

Here, σ_{x }and σ_{y }are the standard deviations in x and y direction, respectively, of the elements included in the thresholded scene, and σ_{g }is the standard deviation of the grey values when implementing the threshold. From the seed frame a threshold is deduced from this segmented image: min_{x}, min_{y}, min_{g }are derived, and min_{n }are the numbers of elements found from this segmented frame.

After thresholding the scene with the value defined from expression 23, the endocardial contours for the left ventricle or right ventricle are extracted. The final step in segmenting the left ventricles endocardial contour is the fine tune step. Therein two methods are presented. In one method the endocardial contour is fine tuned by means of a convex hull operation. In the second method the endocardial contour is forced into a smooth circular behaviour, which is more realistic from a physiologic point of view. This applies especially when analysing regional functional parameters of the left ventricle, such as wall thickness, wall thickening and wall motion. It might mean that the endocardial volume is overestimated, for example in that the papillary muscle volume is added to the left ventricle blood volume. Deducing the papillary muscles and subtracting their volume from the endocardial volume will calculate the correct blood volume, see hereinafter. During the detection of the right ventricle endocardial contour the fine tune step is always the first method. The right ventricles epicardial boundary is not detected by the current procedure, caused by the limited spatial resolution of the MR images.

Block 12 presents the left ventricle epicardial segmentation. This will segment the epicardial contour as based on a radial minimum cost algorithm. The detected endocardial contour is used as seed data. From the endocardial contour a polar model is created with an origin wherefrom radial lines originate and in which the radials are separated by n degrees. The origin is the centre of gravity of the endocardial contour. This polar model is used to resample the original image.

The epicardial boundary also contains various transitions in grey levels, alternating from BRIGHT_TO_DARK and vice versa. Because we use derivatives for edge detection, the sign of the edge should be given. Thereto, the following approach with dynamic sign definition is introduced, as follows. For each polar scanline:

 Compute the location of the position whereas the absolute derivative is maximum.
$\begin{array}{cc}\mathrm{pos}=\mathrm{arg}\text{\hspace{1em}}{\mathrm{max}}_{x}\left\{\uf603\frac{\partial}{\partial x}g\left(x\right)\uf604\right\}& \left(24\right)\end{array}$
 g represents the scanline
 Define sign by:
$\begin{array}{cc}\begin{array}{c}\begin{array}{c}\begin{array}{c}{l}_{\mathrm{mean}}=\frac{\sum _{i=0}^{\mathrm{pos}}g\left(i\right)}{\mathrm{pos}}\\ {r}_{\mathrm{mean}}=\frac{\sum _{i=\mathrm{pos}}^{\mathrm{resample\_width}}g\left(i\right)}{\mathrm{resample\_width}\mathrm{pos}}\end{array}\\ \mathrm{if}\text{\hspace{1em}}\left({l}_{\mathrm{mean}}<{r}_{\mathrm{mean}}\right)\Rightarrow \mathrm{sign}=1\iff \mathrm{DARK\_TO}\mathrm{\_BRIGHT}\end{array}\\ \mathrm{if}\text{\hspace{1em}}\left({l}_{\mathrm{mean}}\ge {r}_{\mathrm{mean}}\right)\Rightarrow \mathrm{sign}=1\iff \mathrm{BRIGHT\_TO}\mathrm{\_DARK}\end{array}& \left(25\right)\end{array}$

The epicard contour is now described by:
$\begin{array}{cc}\begin{array}{c}{C}_{s}\left(i\right)=\mathrm{min}\left({C}_{s}\left(i,j\right)\u2758\mathrm{samples}\text{\hspace{1em}}j\text{\hspace{1em}}\mathrm{on}\text{\hspace{1em}}\mathrm{scanline}\right)\\ {C}_{S}\left(i,j\right)=\left\{\begin{array}{c}(\alpha \xb7\left(\mathrm{sign}\text{\hspace{1em}}\left(j\right)\xb7\left(\frac{\partial}{\partial i}{g}_{j}\left(i\right)\right)\right)+\\ \left(1\alpha \right)\left(\mathrm{sign}\text{\hspace{1em}}\left(j\right)\xb7\left(\frac{{\partial}^{2}}{{\partial}_{i}^{2}}{g}_{j}\left(i\right)\right)\right))+\frac{1}{2}n\left(i,j\right)\end{array}\right\}\xb7\mathrm{grey\_cos}\text{\hspace{1em}}t\left(i,j\right)\\ n\left(i,j\right)=\uf603\frac{\partial}{\partial t}{g}_{{\mathrm{previous}}_{j}}\left(i\right)\uf604+\uf603\frac{\partial}{\partial i}{g}_{{\mathrm{next}}_{j}}\left(i\right)\uf604\\ \mathrm{grey\_cos}\text{\hspace{1em}}t\left(i,j\right)=\{\begin{array}{cc}{e}^{0.5\left(\frac{{\left({g}_{j}\left(i\right)\left(\mu \sigma \right)\right)}^{2}}{{\sigma}^{2}}\right)}& \mathrm{if}\text{\hspace{1em}}\mathrm{sign}=\mathrm{BRIGHT\_TO}\mathrm{\_DARK}\\ {e}^{0.5\left(\frac{{\left({g}_{j}\left(i\right)\left(\mu +\sigma \right)\right)}^{2}}{{\sigma}^{2}}\right)}& \mathrm{otherwise}\end{array}\end{array}& \left(26\right)\end{array}$

Based on the sign determination described earlier, the posterior junction of the right ventricle and the interventricular septum (the LVRV connection) is defined in block 13. For each frame inside the shortaxis dataset the area is defined with the largest closed DARK_TO_BRIGHT sign by using expression 24 and 25. Then the first scan line sign in this frame going from dark to bright in a clockwise resample orientation of the above polar model is defined as the LVRV connection point.

An alternative segmentation for the left ventricle is presented in block 16. This method is based on expression 25. Therein, both the left ventricle endocardial and epicardial borders are segmented with a minimum cost approach. The main difference is that this method assumes an approximately correct location of a seed contour. From this assumption, segmentation propagates in time with each shortaxis frame. The previous location of the contours is introduced into the cost function. Also the sign definition and tissue statistics are updated for each frame. The cost function is given by:
$\begin{array}{cc}\begin{array}{c}{C}_{S}\left(i,j\right)=\begin{array}{c}\left(\begin{array}{c}\alpha \xb7\left(\mathrm{sign}\text{\hspace{1em}}\left(j\right)\xb7\left(\frac{\partial}{\partial i}{g}_{j}\left(i\right)\right)\right)+\\ \left(1a\right)\left(\mathrm{sign}\text{\hspace{1em}}\left(j\right)\xb7\left(\frac{{\partial}^{2}}{{\partial}_{i}^{2}}{g}_{j}\left(i\right)\right)\right)\end{array}\right)\xb7\\ \mathrm{grey\_cos}\text{\hspace{1em}}t\left(i,j\right)\xb7\mathrm{prev\_cos}\text{\hspace{1em}}t\left(i,j\right)\end{array}\\ \mathrm{where}\text{\hspace{1em}}\mathrm{prev\_cos}\text{\hspace{1em}}t\left(i,j\right)={e}^{\left(\frac{{\left({\mathrm{contour\_position}}_{j}i\right)}^{2}}{{\sigma}_{p}^{2}}\right)}\end{array}& \left(27\right)\end{array}$
This approach is suitable for contour corrections in shortaxis images.

The final step in segmenting the left ventricle is the identifying/segmenting of the papillary muscles in block 14, as follows. From knowing the boundary of the myocardium, the endocardial contour and the epicardial contour, the grey level distribution of the myocardium and the grey level associated to the blood pool can be derived accurately. This is done by computing the normal distribution of both tissues; then the associated threshold separates both distributions by using the Otsu threshold method in each frame I_{p}. In places where papillary muscles are located no blood can be present. The segmentation of the papillary is then reduced to thresholding the area given by the endocardial contour and extracting all segmented objects which exceed a predefined area based on binary morphology.

Block 15 presents a further improvement in the use of the twochamber and fourchamber image set to define the 3D spatial left ventricle virtual apex position. Due to partial volume effects and the relative thick shortaxis frame thickness used during the acquisition, the shortaxis apex is difficult to localize, even by an experienced clinician. Now, in current state of the MRI art cine acquisition slices are about 8 mm thick, to have sufficient signaltonoise ratio. Knowing the 3D spatial position of the left ventricle apex, it is possible to afterwards correct the volume of the left ventricle.

The spatial 3D apex position is computed in the enddiastolic phase and in the endsystolic phase as identified by expression 17. The derived endocardial contours in the shortaxis frames in the enddiastolic phase are mapped on the two chamber and four chamber frames starting from the base slice to the apex slice as found when localizing the left ventricle. The mapping is done by:
p _{2ch}[ ]endo∩l _{2ch }
p _{4ch}[ ]endo∩l _{4ch} (28)

The point lists p_{2ch }and p_{4ch }are mapped on to the long axis frames by transforming them first into the 3D world and then back onto the long axis plane as follows:
$\begin{array}{cc}{P}_{3D}=T\left({P}_{\mathrm{SA}}\right)\xb7{R}_{\mathrm{SA}}\xb7S\left({\mathrm{voxelsize}}_{\mathrm{SA}}\right)\xb7{\mathrm{point}}_{2\text{\hspace{1em}}D}& \left(29\right)\\ {\mathrm{point}}_{2\text{\hspace{1em}}{D}_{\mathrm{LA}}}=\left[\begin{array}{c}1\\ 1\\ 0\end{array}\right]\xb7S\left(1/{\mathrm{voxelsize}}_{\mathrm{LA}}\right)\xb7{R}_{\mathrm{LA}}\xb7T\left({P}_{\mathrm{LA}}\right)\xb7{P}_{3D}& \left(30\right)\end{array}$

The above steps repeat for all short axis frames until reaching the apex shortaxis slice. In this respect, FIG. 6 illustrates a mapping procedure from longaxes (LA) slices to shortaxis (SA) slices or inversely. Through the various nodes indicated by black dots for both 2CH and 4CH cine sequences a wellknown CatmullRow spline curve is fitted. Also the apex containing slice has been shown:
$\begin{array}{cc}s\left(t\right)\Rightarrow \left[\begin{array}{c}x={a}_{x}{t}^{3}+{b}_{x}{t}^{2}+{c}_{x}t+{d}_{x}\\ y={a}_{y}{t}^{3}+{b}_{y}{t}^{2}+{c}_{y}t+{d}_{y}\end{array}\right]\text{}t\in \left[0,1\right]& \left(31\right)\end{array}$

This model is used to segment the left ventricles endocardial in the long axis view of the enddiastolic phase. Again, a radial minimum cost algorithm localizes the endocardial border, wherein both the 2CH and 4CH cine sequences contribute two border points, in the lower Figure part. The following function (expression 32) defines the boundary, by using CS as a cost function. Herein, g represents the resample image by the polar model.
$\begin{array}{cc}{C}_{s}\left(i\right)=\mathrm{min}\left({C}_{s}\left(i,j\right)\mathrm{samples}\text{\hspace{1em}}j\text{\hspace{1em}}\mathrm{on}\text{\hspace{1em}}\mathrm{scanline}\right)\text{}{C}_{s}\left(i,j\right)=\begin{array}{c}\left(\alpha \xb7\frac{\partial}{{\partial}_{i}}{g}_{j}\left(i\right)+\left(1\alpha \right)\frac{{\partial}^{2}}{{\partial}_{i}^{2}}{g}_{j}\left(i\right)\right)\xb7\\ \left({I}_{\mathrm{ES}}{I}_{\mathrm{ED}}\right)\xb7{e}^{0.5\left(\frac{{\left({g}_{j}\left(i\right)\left(\mu +\sigma \right)\right)}^{2}}{{\sigma}^{2}}\right)}\end{array}& \left(32\right)\end{array}$

The apex position is now located (expression 33) in the twochamber and fourchamber for both long axis views, the 3D spatial apex being defined by expression 34:
$\begin{array}{cc}{\overrightarrow{p}}_{\mathrm{midbase}}=\frac{{\overrightarrow{p}}_{{1}_{\mathrm{base}}}+{\overrightarrow{p}}_{{2}_{\mathrm{base}}}}{2}\text{}{\overrightarrow{p}}_{{\mathrm{apex}}_{\mathrm{LA}}}=\mathrm{max}\text{\hspace{1em}}\mathrm{dist}\text{\hspace{1em}}\left({\mathrm{endo}}_{{\mathrm{LA}}_{i}}{\overrightarrow{p}}_{\mathrm{midbase}}\right)& \left(33\right)\\ {\overrightarrow{p}}_{\mathrm{apex}}=\frac{\begin{array}{c}T\left({P}_{2\mathrm{CH}}\right)\xb7{R}_{2\mathrm{CH}}\xb7s\left({\mathrm{voxelsize}}_{2\mathrm{CH}}\right)\xb7{\overrightarrow{p}}_{{\mathrm{apex}}_{2\mathrm{CH}}}+\\ T\left({P}_{4\mathrm{CH}}\right)\xb7{R}_{4\mathrm{CH}}\xb7s\left({\mathrm{voxelsize}}_{4\mathrm{CH}}\right)\xb7{\overrightarrow{p}}_{{\mathrm{apex}}_{4\mathrm{CH}}}\end{array}}{2}& \left(34\right)\end{array}$

The same steps are repeated for the endsystolic phase. The apex position in the remaining phases is interpolated or extrapolated from the spatial 3D apex positions in both the enddiastolic and endsystolic phases: this will present the correct volume contribution for the apex slice, where the slice may, or may not intersect the apex position proper. In fact, the left ventricle volume computed through Simpson's rule is now corrected by fitting an ellipsoid through the 3D apex position and the apex slice.

In this respect, FIG. 7 illustrates the starting of segmentation through presenting user LA input. Briefly the user identifies the start conditions in the two and fourchamber image dataset. The advantage of this method is that the volume can be corrected according to the spatial location of the mitral valve derived from the long axis views as well. Now first, in block 70, the endocardial and epicardial left ventricle border in the 2chamber and 4chamber view both the enddiastole and endsystole cardiac heart phase must be present, for instance by manually drawing a spline based contour based on expression 32. Although, this represents operator input, all of the further blocks in FIG. 7 are executed automatically. Now, in block 71, ‘inplane’ movement and breathing movement are corrected as described with reference to blocks 41 and 42 of FIG. 4. Next, in block 72 these long axis contours are mapped on the shortaxis dataset. In pseudo code:

Walk through all slices
l _{2ch} =SA(slice)∩2CH Transformed into the 2CH
l _{4ch} =SA(slice)∩4CH Transformed into the 4CH
pt _{2ch} [ ]=s _{2ch}(t)∩l _{2ch }s(t) 2chamber spline
pt _{4ch} [ ]=s _{4ch}(t)∩l _{4ch }s(t) 4chamber spline (35)

Convert the point—lists pt_{2ch }and pt_{4ch }into the SA slice

Define the slices which are part of the ventricle

The segmenting of the left ventricle endocardial (block 11), epicardial (block 12), LVRV connection point (block 13) and the papillary muscles (block 14) are identical to previously described procedures. Next, in block 73, the volume as computed by Simpson's rule, corrected according to the spatial left ventricle apex and mitral valve position while excluding the volume ‘outside’ the left ventricle. Thereupon, the procedure is finished again.

FIG. 8 illustrates a high level functional block diagram of an exemplary MRI apparatus. For clarity, only three subsystems have been shown. Block 112 represents the MRI apparatus proper that operates under commands from user interface module 116 and will provide data to data processing module 114. The latter may be represented by any computer processing platform, including a workstation, personal computer and the like; it executes the generating of intermediate data and final results to user interface module 116, the latter furthermore allowing a dialog with data processing module 114. User interface 116 represents all kinds of input and output devices; it may also be embodied in a program storage device (e.g., CDROM or DVD or in one or more downloadable files) that store a computer program that is loaded onto and executed by a computer processing platform to carry out the inventive operations described herein. Further steps and procedures have either been detailed supra, or have been common knowledge to persons skilled in the art.

There have been described and illustrated herein several embodiments of a machinebased methodology for segmenting the image of the chambers of the human heart. While particular embodiments of the invention have been described, it is not intended that the invention be limited thereto, as it is intended that the invention be as broad in scope as the art will allow and that the specification be read likewise. It will therefore be appreciated by those skilled in the art that yet other modifications could be made to the provided invention without deviating from its spirit and scope as claimed.