US20070253609A1

US20070253609A1 - Method, Apparatus and Computer Program Product for Automatic Segmenting of Cardiac Chambers

Info

Publication number: US20070253609A1
Application number: US11/380,681
Authority: US
Inventors: Jean-Paul Aben
Original assignee: Pie Medical Imaging BV
Current assignee: Pie Medical Imaging BV
Priority date: 2006-04-28
Filing date: 2006-04-28
Publication date: 2007-11-01
Also published as: US7864997B2

Abstract

A method, apparatus and computer program product for machine-based segmentation of the cardiac chambers in a four-dimensional magnetic resonance image data that contains at least short-axis cine sequences. The segmenting is based on temporal variations in the data set and on features extracted from the image data. In particular, the method and apparatus process the image data to automatically segment an endocard based on fuzzy object extraction and to automatically segment an epicard through finding a radial minimum cost with dynamic sign definition and tissue classification. In addition, image data is preferably preprocessed to eliminate high local image variances and inhomogeneity. The image data is also preferably processed to identify positions of a base slice, an apex slice, an end-diastole, an end-systole and an LVRV connection.

Description

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 four-dimensional 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 short-axis plane of the ventricle. These short-axis 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 non-invasive 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 region-of-interest computation in cine 3D structural MRI”, CARS and Elsevier, 2004, describes a method for extracting a region of interest around the heart in short-axis MRI cardiac scans. The method produces a region within four-dimensional (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 short-axis 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 so-called “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. 321-332, 1988. Briefly, the snake model refers to an energy minimization technique that seeks for the lowest potential of a curve-based 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 user-defined 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 multi-phase, multi-slice 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 subject-caused 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 pick-up operations. Moreover, the segmentation of the endocardial and epicardial features starts from an operator-defined 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 machine-controlled, 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 in-plane movement which is identified and corrected by the methodology of the present invention;
FIG. 3 is a three-dimensional 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 long-axes (LA) slices to short-axis (SA) slices or inversely;
FIG. 7 is a flow chart illustrating segmentation operations that start from user-supplied long-axis (LA) input;
FIG. 8 is a high-level 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 two-chamber and four-chamber cine sequences.
As an alternative to the above, it is possible to work with short-axis 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_pis the 2D image array that is subjected to conventional Fourier transform methods to render the result more reliable. $\begin{matrix} I_{p} = 𝒥^{- 1} {𝒥 (I Θ b) \cdot H} \cdot T (- P_{shift}) where H = 1 - ⅇ^{\frac{- r^{2}}{2 \cdot \partial^{2}}} in case of inhomogeneity H = 1 otherwise b is a square structured element P_{shift} is the displacing caused by acquisition & (1) \end{matrix}$
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_2chand l_4chrespectively, in the first frame of the middle short-axis slice. This is followed by executing a first Gaussian derivative operation.
From the centre, the positions of the steepest absolute slopes l_pand r_pto 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 short-axis data set typically 8 mm thick cine short-axis slices with breath-holding 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 ‘in-plane’ movements.
In this respect, FIG. 2 illustrates a spatial configuration of this in-plane 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 in-plane 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 “in-plane” movement. Such operations begin by identifying the starting plane (102) and then computes the three-dimensional (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 ₁ =P ₁ +μR ₁ +γ·C ₁
l=P ₁+(μ×·C ₁)·α (2)
Next, the intersection point P_3Dof 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_3Dis transformed into the two-dimensional (2D) coordinate system of the successor plane: $\begin{matrix} P_{2 D} = [\begin{matrix} 1 \\ 1 \\ 0 \end{matrix}] \cdot S (1 / {voxelsize}_{n}) \cdot R_{n} \cdot T (- P_{n}) \cdot P_{3 D} & (4) \end{matrix}$
Here, the rotation matrix R is defined as the cross product of the plane direction vectors: $\begin{matrix} R = [\begin{matrix} R_{1 x} & C_{1 x} & N_{1 x} \\ R_{1 y} & C_{1 y} & N_{1 y} \\ R_{1 z} & C_{1 z} & N_{1 z} \end{matrix}] N = R \otimes C & (5) \end{matrix}$
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 short-axis stack.
In addition to such ‘in-plane’ movement, the patient can move between the acquired short-axis cine sequences, the two-chamber and the four-chamber cine sequence. This moving is specific for each examination and its magnitude depends mainly on how well the patient reproduces each breath-holding 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 three-dimensional 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 two-chamber 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 two-chamber image is created by resampling each short-axis image along the intersection line (epipolar line) of the two-chamber image, which gives an image which can be used for registration against the true two-chamber plane found during the acquisition. Breathing movement can now be computed by correlating in the Fourier domain each resampled line against the corresponding two-chamber intersection. This method can be extended with the four-chamber view.
In block 43, after preprocessing all short-axis frames, the total variance image of the complete short-axis dataset is generated. Because the short-axis images are perpendicular to the long axis of the ventricle (its z-axis 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{matrix} VAR = \frac{1}{n_frames} \sum_{i = 1}^{n_frames} {(I_{p_{i}} - {\overline{I}}_{p})}^{2} I_{variance} = MAX (VAR (short axis {slice}_{n})) + \frac{\sum_{n = 1}^{n_slices} (VAR (short axis {slice}_{n})}{n_slices} & (6) \end{matrix}$
The mean intersection points derived from the two-chamber and four-chamber epipolar lines projected on each short-axis slice are smoothed and added to the overall variance image: $\begin{matrix} \overline{p} = \frac{\sum_{1}^{n_slices} l_{2 ch} ⋂ l_{4 ch}}{n_slices} l_{2 ch} = V_{{SA}_{i}} ⋐ V_{2 ch} mapped on {SA}_{i} l_{4 ch} = V_{{SA}_{i}} ⋐ V_{4 ch} mapped on {SA}_{i} \partial_{p} = STDEV (p) & (7) \end{matrix}$
The mapping back on the SA slice is performed by: $\begin{matrix} P_{2 D_{SA}} = [\begin{matrix} 1 \\ 1 \\ 0 \end{matrix}] \cdot S (1 / {voxelsize}_{SA}) \cdot R_{SA} \cdot T (- P_{SA}) \cdot P_{3 D} & (8) \end{matrix}$
The variance image is now adjusted with this information according to: $\begin{matrix} I_{totalvariance} = I_{variance} \cdot \sum_{i = p_{y} - \partial_{p_{y}}}^{p_{y} + \partial_{p_{y}}} \sum_{p_{x} - \partial_{p_{x}}}^{p_{x} + \partial_{p_{x}}} I_{{variance}_{(i, j)}} ⅇ^{- (\frac{(x - {\overline{p}}_{x})}{2 \cdot \partial_{p_{x}}} + \frac{(y - {\overline{p}}_{x})}{2 \cdot \partial_{py}})} & (9) \end{matrix}$
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 grey-level histograms” IEEE Trans. Syst. Man Cybernet SMC-01 (1970) pp. 62-66, 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 short-axis dataset, defined as mask M. The short-axis slices are perpendicular to the long axis of the left ventricle, which allows to approximate the left ventricle by a circle in these short-axis 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 well-known 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 short-axis slice. Within this slice the most likely end-diastole 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 so-called Canny edge detector, which is well-known 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{matrix} \nabla I = \sqrt{\langle G_{x}^{2} + G_{Y}^{2} \rangle} ⋂ M θ = \tan^{- 1} (Gy / G_{x}) ⋂ M & (11) \end{matrix}$
The kernels used to derive G_xand G_yare: $\begin{matrix} G_{X} = [\begin{matrix} - 1 & - 2 & - 1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \end{matrix}], G_{y} = [\begin{matrix} - 1 & 0 & 1 \\ - 2 & 0 & 2 \\ - 1 & 0 & 1 \end{matrix}] & (12) \end{matrix}$
The gradient magnitude image ∇I is binarized to I_binusing 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_binis parsed and suggested to be on a circle boundary with radius r between r_minand 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{matrix} \nabla^{2} H = H \cdot (\frac{(x^{2} + y^{2}) - \partial^{2}}{σ^{4}}) \cdot ⅇ^{- (\frac{(x^{2} + y^{2})}{2 σ^{2}})} ⋂ M & (14) \end{matrix}$

Now the position with the highest value in ∇²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_pis 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_mit 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{matrix} \begin{matrix} \begin{matrix} \begin{matrix} Single Otsu on H_{p_{M}} \Rightarrow T_{M}; (μ_{M_{0}}, σ_{M_{0}}); (μ_{M_{1}}, σ_{M_{1}}) \\ Dual Otsu on H_{p_{M}} \Rightarrow T_{{Md}_{1}}; T_{{Md}_{2}}; (μ_{{Md}_{0}}, σ_{{Md}_{0}}); (μ_{{Md}_{1}}, σ_{{Md}_{1}}); \end{matrix} \\ (μ_{{Md}_{2}}, σ_{{Md}_{2})} \end{matrix} \\ Single Otsu on H_{p_{circle}} \Rightarrow T_{C}; (μ_{C_{0}}, σ_{C_{0}}); (μ_{C_{1}}, σ_{C_{1}}) \\ T = {\begin{matrix} μ_{C_{1}} - σ_{C_{1}} & if (T_{M} > μ_{C_{1}})^(T_{{Md}_{1}} > μ_{C_{1}}) \\ T_{{Md}_{1}} & if (T_{M} > μ_{C_{1}})^(T_{{Md}_{1}} \leq μ_{C_{1}}) \\ T_{M} & otherwise \end{matrix} \end{matrix} & (16) \end{matrix}$
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 four-chamber 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 end-diastole and end-systole 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_iis the area of the object extracted with the threshold T obtained through applying expression 16 on the frame i from the current short-axis slice. Again the binarization is performed within the mask M.
Finally, in block 47 the base short-axis slice and apex short-axis slice are defined at the end-diastolic phase as computed by expression 17. For each frame I_pwithin the end-diastolic 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_paccording to expression 16 that gives B_i.
- Intersect with mask M
- Fill holes located in B_iby using a morphology flood fill algorithm.
- Extract the object that obtains the centre position C_pThe extracted object gives a global estimation of the endocardial contour location B_endo _i
- Define Area_iand perimeter_iof extracted object B_endo _i
- Compute the circularity (C) of the extracted object by: $\begin{matrix} C = \frac{4 \cdot π \cdot {Area}_{i}}{{Perimeter}_{i}^{2}} & (18) \end{matrix}$
  - 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_pis adjusted according to: $\begin{matrix} C_{p} \Rightarrow {\begin{matrix} x = \frac{\sum_{pixels} x}{{Area}_{i}} \\ y = \frac{\sum_{pixels} y}{{Area}_{i}} \end{matrix} & (19) \end{matrix}$
- 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{matrix} \begin{matrix} \overline{A} rea 1 = \frac{\sum_{i = {seedframeb}_{1}}^{{Vb}_{1}} {Area}_{i}}{{Vb}_{1} - seedframe} \\ \overline{A} rea 2 = \frac{\sum_{i = {Vb}_{2}}^{seedframe} {Area}_{i}}{seedframe - {Vb}_{2}} \\ apex_slice = {\begin{matrix} {Vb}_{1} & if \overline{A} rea 1 < \overline{A} rea 2 \\ {Vb}_{2} & otherwise \end{matrix} \\ base_slice = {\begin{matrix} {Vb}_{1} & if \overline{A} rea 1 \geq \overline{A} rea 2 \\ {Vb}_{2} & otherwise \end{matrix} \end{matrix} & (20) \end{matrix}$
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_paccording 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 end-diastolic phase, end-systolic 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 disc-structured element is used with radius one. For the segmentation a fuzzy connectedness algorithm can be used. The algorithm is preferably based upon on the well-known 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, “Fuzzy-Connectedness and Object Definition: Theory, Algorithms and Applications in Image Segmentation” Graphic Models, Image Processing, vol. 58, 246-261, 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) 4-connectivity 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_Cand P_Dare 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₂. 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{matrix} \begin{matrix} F = S \cdot (1 - N_{1}) \\ S = {\begin{matrix} N_{2} & if P \leq μ_{2} \\ N_{2} & with σ = 3 σ_{2} if P > μ_{2} \end{matrix} \end{matrix} & (22) \end{matrix}$
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₂, but for intensities above the mean μ₂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_maxof these values for all 6 neighbours is derived then.
  - c. If the maximum value f_maxis 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 so-called “scene” giving at each position in the 3D image/stack the connectivity or measure of “hanging-togetherness” 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, σ_xand σ_yare the standard deviations in x and y direction, respectively, of the elements included in the thresholded scene, and σ_gis 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_gare derived, and min_nare 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{matrix} pos = \arg \max_{x} {\langle \frac{\partial}{\partial x} g (x) \rangle} & (24) \end{matrix}$
  - g represents the scanline
- Define sign by: $\begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} l_{mean} = \frac{\sum_{i = 0}^{pos} g (i)}{pos} \\ r_{mean} = \frac{\sum_{i = pos}^{resample_width} g (i)}{resample_width - pos} \end{matrix} \\ if (l_{mean} < r_{mean}) \Rightarrow sign = 1 \Leftrightarrow DARK_TO_BRIGHT \end{matrix} \\ if (l_{mean} \geq r_{mean}) \Rightarrow sign = - 1 \Leftrightarrow BRIGHT_TO_DARK \end{matrix} & (25) \end{matrix}$

The epicard contour is now described by: $\begin{matrix} \begin{matrix} C_{s} (i) = \min (- C_{s} (i, j) ❘ samples j on scanline) \\ C_{S} (i, j) = {\begin{matrix} (α \cdot (sign (j) \cdot (\frac{\partial}{\partial i} g_{j} (i))) + \\ (1 - α) (sign (j) \cdot (\frac{\partial^{2}}{\partial_{i}^{2}} g_{j} (i)))) + \frac{1}{2} n (i, j) \end{matrix}} \cdot grey_cos t (i, j) \\ n (i, j) = \langle \frac{\partial}{\partial t} g_{{previous}_{j}} (i) \rangle + \langle \frac{\partial}{\partial i} g_{{next}_{j}} (i) \rangle \\ grey_cos t (i, j) = {\begin{matrix} ⅇ^{- 0.5 (\frac{{(g_{j} (i) - (μ - σ))}^{2}}{σ^{2}})} & if sign = BRIGHT_TO_DARK \\ ⅇ^{- 0.5 (\frac{{(g_{j} (i) - (μ + σ))}^{2}}{σ^{2}})} & otherwise \end{matrix} \end{matrix} & (26) \end{matrix}$
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 short-axis 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 short-axis 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{matrix} \begin{matrix} C_{S} (i, j) = \begin{matrix} (\begin{matrix} α \cdot (sign (j) \cdot (\frac{\partial}{\partial i} g_{j} (i))) + \\ (1 - a) (sign (j) \cdot (\frac{\partial^{2}}{\partial_{i}^{2}} g_{j} (i))) \end{matrix}) \cdot \\ grey_cos t (i, j) \cdot prev_cos t (i, j) \end{matrix} \\ where prev_cos t (i, j) = ⅇ^{(\frac{{({contour_position}_{j} - i)}^{2}}{σ_{p}^{2}})} \end{matrix} & (27) \end{matrix}$
This approach is suitable for contour corrections in short-axis 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 two-chamber and four-chamber image set to define the 3D spatial left ventricle virtual apex position. Due to partial volume effects and the relative thick short-axis frame thickness used during the acquisition, the short-axis 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 signal-to-noise 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 end-diastolic phase and in the end-systolic phase as identified by expression 17. The derived endocardial contours in the short-axis frames in the end-diastolic 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_2chand p_4chare 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{matrix} P_{3 D} = T (P_{SA}) \cdot R_{SA} \cdot S ({voxelsize}_{SA}) \cdot {point}_{2 D} & (29) \\ {point}_{2 D_{LA}} = [\begin{matrix} 1 \\ 1 \\ 0 \end{matrix}] \cdot S (1 / {voxelsize}_{LA}) \cdot R_{LA} \cdot T (P_{LA}) \cdot P_{3 D} & (30) \end{matrix}$
The above steps repeat for all short axis frames until reaching the apex short-axis slice. In this respect, FIG. 6 illustrates a mapping procedure from long-axes (LA) slices to short-axis (SA) slices or inversely. Through the various nodes indicated by black dots for both 2CH and 4CH cine sequences a well-known Catmull-Row spline curve is fitted. Also the apex containing slice has been shown: $\begin{matrix} s (t) \Rightarrow [\begin{matrix} 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{matrix}] t \in [0, 1] & (31) \end{matrix}$
This model is used to segment the left ventricles endocardial in the long axis view of the end-diastolic 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{matrix} C_{s} (i) = \min (- C_{s} (i, j) | samples j on scanline) C_{s} (i, j) = \begin{matrix} (α \cdot \frac{\partial}{\partial_{i}} g_{j} (i) + (1 - α) \frac{\partial^{2}}{\partial_{i}^{2}} g_{j} (i)) \cdot \\ (I_{ES} - I_{ED}) \cdot ⅇ^{- 0.5 (\frac{{(g_{j} (i) - (μ + σ))}^{2}}{σ^{2}})} \end{matrix} & (32) \end{matrix}$
The apex position is now located (expression 33) in the two-chamber and four-chamber for both long axis views, the 3D spatial apex being defined by expression 34: $\begin{matrix} {\vec{p}}_{midbase} = \frac{{\vec{p}}_{1_{base}} + {\vec{p}}_{2_{base}}}{2} {\vec{p}}_{{apex}_{LA}} = \max dist ({endo}_{{LA}_{i}} - {\vec{p}}_{midbase}) & (33) \\ {\vec{p}}_{apex} = \frac{\begin{matrix} T (P_{2 CH}) \cdot R_{2 CH} \cdot s ({voxelsize}_{2 CH}) \cdot {\vec{p}}_{{apex}_{2 CH}} + \\ T (P_{4 CH}) \cdot R_{4 CH} \cdot s ({voxelsize}_{4 CH}) \cdot {\vec{p}}_{{apex}_{4 CH}} \end{matrix}}{2} & (34) \end{matrix}$
The same steps are repeated for the end-systolic phase. The apex position in the remaining phases is interpolated or extrapolated from the spatial 3D apex positions in both the end-diastolic and end-systolic 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 four-chamber 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 2-chamber and 4-chamber view both the end-diastole and end-systole 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, ‘in-plane’ 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 short-axis 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 _2chs(t) 2chamber spline
pt _4ch [ ]=s _4ch(t)∩l _4chs(t) 4chamber spline (35)
Convert the point—lists pt_2chand pt_4chinto 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., CD-ROM 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 machine-based 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.

Claims

1. A method for machine-based segmentation of cardiac chambers in a four-dimensional magnetic resonance image that contains at least short-axis cine sequences, the segmentation based on temporal variations in the image, the method comprising:

preprocessing image data to eliminate high local image variances and in homogeneity;

processing image data to identify positions of a base slice, an apex slice, an end-diastole phase, an end-systole phase and an LVRV connection, the LVRV connection being a posterior junction between a right ventricle and the interventricular septum;

followed by machine-based segmentation of an endocard based on fuzzy object extraction, and machine-based segmentation of an epicard through finding a radial minimum cost with dynamic sign definition and tissue classification.

2. A method as claimed in claim 1, wherein:

said cine sequences furthermore comprise two-chamber cine sequences and four-chamber cine sequences.

3. A method as claimed in claim 1, further comprising:

correcting for data acquisition movements.

4. A method as claimed in claim 1, further comprising:

automatic segmentation of the papillary muscles of the heart.

5. A method as claimed in claim 4, further comprising:

automatic detection of the position of the LVRV connection.

6. A method as claimed in claim 1, further comprising:

allowing for operator inspection of a subject-caused movement in the acquired MR image dataset by using system coordinates from an MR scanner apparatus.

7. A method as claimed in claim 6, further comprising:

automatic correction for subject-caused movement in the acquired MR image dataset by using system coordinates from an MR scanner apparatus.

8. A method as claimed in claim 1, wherein:

segmentation of the endocard and epicard starts from an operator-defined seed contour and being based on radial minimum cost with dynamic sign definition and tissue classification.

9. A method as claimed in claim 1, wherein:

volume of a ventricle is determined by a correction that is derived from the spatial location of the heart's apex position on a curve which is obtained from the long axis cine sequences.

10. A method as claimed in claim 9, wherein:

said correction takes into account a mitral valve position obtained from long axis cine views.

11. An apparatus for machine-based segmentation of cardiac chambers in a four-dimensional magnetic resonance image that contains at least short-axis cine sequences, the segmentation based on temporal variations in the data set and on features extracted from the image, the apparatus comprising:

preprocessing means for preprocessing image data to eliminate high local image variances and inhomogeneity;

means for processing image data to identify positions of a base slice, an apex slice, an end-diastole, an end-systole and an LVRV connection, the LVRV connection being a posterior junction between a right ventricle and an interventricular septum, and

means for machine-based segmentation of an endocard based on fuzzy object extraction, and means for machine-based segmentation of an epicard through finding a radial minimum cost with dynamic sign definition and tissue classification.

12. An apparatus as claimed in claim 11, wherein:

13. An apparatus as claimed in claim 11, wherein:

machine-based segmentation over time of endocardial, epicardial and papillary muscles is based on a seed contour.

14. A computer program product that is installable onto a computer processing platform, the computer program product readable by the computer processing platform, tangibly embodying a program of instructions executable by the computer processing platform to perform methods steps for machine-based segmentation of the cardiac chambers in four-dimensional magnetic resonance image that contains at least short-axis cine sequences, the segmentation based on temporal variations in the image and on features extracted from the image, the method steps comprising:

15. A computer program product as claimed in claim 14, the method steps further comprising:

correcting for data acquisition movements.

16. A computer program product as claimed in claim 14, the method steps further comprising:

automatic segmentation of the papillary muscles of the heart.

17. A computer program product as claimed in claim 14, the method steps further comprising:

18. A computer program product as claimed in claim 17, the method steps further comprising:

19. A method for machine-based segmentation of cardiac chambers in a four-dimensional magnetic resonance image that contains short-axis cine sequences, two-chamber cine sequences and four-chamber cine sequences, the segmentation based on temporal variations in the image and on features extracted from the image, the method comprising:

interacting with an operator to indicate the epicard and endocard in a long axis data set; and

automatic segmentation of the endocard based on fuzzy object extraction, and automatic segmentation of the epicard based on radial minimum cost with dynamic sign definition and tissue classification.

20. An apparatus for machine-based segmentation of cardiac chambers in a four-dimensional magnetic resonance image that contains short-axis cine sequences, two-chamber cine sequences and four-chamber cine sequences, the segmentation based on temporal variations in the image and on features extracted from the image, the apparatus comprising:

means for interacting with an operator to indicate the epicard and endocard in a long axis data set; and

means for automatic segmentation of the endocard based on fuzzy object extraction, and means for automatic segmentation of the epicard based on radial minimum cost with dynamic sign definition and tissue classification.