DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

[0023] In the present system a method comprises quantifying the change in organ morphology using a semi-automatic segmentation method. A user approximately positions a model of the organ (or organs) in an image volume. After this, the process is fully automatic. Steps include intersecting the user-positioned model with the individual image planes of the volume to create a series of contours approximating the position of the underlying organ. Then, optimal 2-D segmentation algorithms are used to refine the boundary of the organ using the approximate contours as starting points. From these optimal contours a new model is applied. Volume and other morphometric measurements are then taken on the model. It is to be noted that in an alternative embodiment a model may be automatically positioned in an image volume.

[0024] For the case of recovering the left ventricle (LV) of the heart, for each stage of the cardiac cycle the following steps are taken: First, the approximate position of the LV is specified by the user. The approximate inner and outer boundaries are outlined with a series mouse clicks on both a short axis image as well as a long axis image, the model is automatically fit to these points. The resulting model is referred to as the initial model. Note that this model is not required to sit exactly on the boundary of the organ. It is enough that it is close to the boundary of the organ. Further, the initial model is intersected with each of the image planes in the image volume. The result is a series of contours (inner and outer) in each of the cross-sectional images of the volume. These individual contours form starting points for optimal 2-D segmentations of the LV. A new model (or the initial one) is automatically fit to these recovered contours. This procedure is then repeated for every temporal instance.

[0025] A cooperative framework includes a method for combining three dimensional (3D) models with two dimensional (2D) segmentations. In one example, a user initializes a parametric component of a hybrid model. A hybrid model is a parametric model with a spline-like mesh on its surface which is capable of deforming to describe fine detail in an image volume. Intersections of image slices with the parametric component serve as starting points for a segmentation of the object boundaries using active contours. The resulting segmentations are sampled to create a data set of 3D points to which the original hybrid model is fit. The hybrid model is free to deform both parametrically and locally in order to represent the object of interest. In this way, the user initializes the hybrid model, the hybrid model initializes the 2D active contours, and the active contours constrain the fit of the hybrid model. This fit is referred to as a physically motivated fit. The process may then be repeated for different time instances.

[0026] Referring to FIGS. 1 - 10 of the present invention and in particular to FIG. 1 which provides a flow chart of the method of the present invention. As described above, block 100 of FIG. 1 represents an initialization of the model. This includes taking a volumetric image of an object to be modeled. The intersection of the models with the imaging planes results in an approximation to the cross-section of the object in the image. This approximation serves as a starting point for active contour segmentations in the image, in block 105 . In block 110 , segmentation using active contours is performed on the object in the image. The segmented contours are used to create a data set of 3D points in block 120 . In block 130 , the model is fit to the data points of block 120 to represent the object. This process is repeated for subsequent time instances as shown in block 131 .

[0027] Referring to FIG. 2 of the drawings, a detailed flow chart is provided for block 100 initializing the model. This step includes identifying approximate positions on the surface from at least two views which bounds the extremities of the object. In the case of the left ventricle of the heart both inner and outer walls (endocardial and epicardial) must be delineated. The approximate positions of the points on the surface of the object are identified in blocks 102 and 106 . The orthogonal cross section is identified in block 104 . If necessary this step is repeated after the step in block 106 . A parametric component of the model is then fitted to the points as shown in block 108 .

[0028] Referring to FIG. 3 of the drawings, in block 111 , the initial model image and, the image planes intersect. FIG. 3 of the drawings illustrates, block 105 of FIG. 1 in a detailed flowchart, in a preferred embodiment a two dimensional image is represented as a graph in block 112 . These intersections are used as the starting points for 2D segmentation. The edge cost is found between graph nodes using image gradients in block 114 . A search space is defined around the approximate input contour in block 116 . A path of lowest energy between pixels/nodes is determined using Dijkstra's algorithm in block 118 . Dijkstra's algorithm is a standard graph algorithm and can be found in any algorithm's book such as C. G. Brassard and P. Bratley, Algorithmics: Theory and Practice, Prentice Hall, 1988. This step is repeated for a different search space if the resulting contour is not closed in block 122 . By determining the lowest energy path a contour (path) is defined in block 124 .

[0029] FIG. 4 of the drawings is a detailed flow chart of block 110 of FIG. 1 illustrating fitting the model of block 130 . This fitting of the model includes determining data and smoothing forces for the points of the data set in block 132 . The data and smoothing forces are scaled in block 134 . Shape functions that reduce the number of iterations needed to fit the model to the data are applied to finite element model (FEM) elements in block 136 . An (FEM) is a finite element model which is made up of a finite number of elements on the surface of the model and serve to describe the model. Data and smoothing forces are combined in block 138 . Applying the adjusted forces to the hybrid model to deform the hybrid model is performed in block 140 .

[0030] Segmentation is among the most widely researched topics in computer vision. A comparison of the cooperative framework described in this disclosure is compared with the most closely related work. However, the method as described is applicable to other techniques and systems.

[0031] A preferred embodiment of a hybrid 3-D model for the heart left ventricle is described in the following. It is understood that this invention is not restricted to this model only. The invention applies to any model of any object. The approach may make use of any 3D model form. In the present example the hybrid model is chosen. Referring to FIG. 5 , for the use of a hybrid model of the LV, the model formulation is hybrid in that it is made up of a parametric component 20 and a local, spline-like mesh component 22 over the parametric component. This local component may deform away from the parametric surface to describe fine detail. The parametric component of the model is based on two ellipsoids, one describing the inner wall of the LV and one describing the outer wall. The model has an underlying parametric component which resembles a double hulled ellipsoid (one for the inner wall of the left ventricle and one for the outer wall). The model is broken up into octants which are described by their own parametric formulations. The quadrants are called PosPosPos, NegPosPos, NegNegPos, PosNegPos, PosPosNeg, NegPosNeg, NegNegNeg, PosNegNeg with representations for the inner and outer walls. This parametric representation is swept out by two intrinsic spherical coordinates, u and v.

[0032] The extrinsic parameters of a1PosInner a1PosOuter rPosPosInner rPosPosOuter a2PosInner a2PosOuter rNegPosInner rNegPosOuter a1NegInner a1NegOuter rNegNegInner rNegNegOuter a2NegInner a2NegOuter rPosNegInner rPosNegOuter a3PosInner describe the general shape of the double hulled ellipsoid.

[0033] PosPosPosInner:

[0034] X=alpha

[0035] Cos[u] (a1PosInner Cos[v]+rPosPosInner Sin[2 v]);

[0036] Y=alpha Cos[u] Sin[v] (a2PosInner+rPosPosInner Sin[v]+rPosPosInner Sin[3 v]);

[0037] Z=a3PosInner alpha Sin[u];

[0038] NegPosPosInner:

[0039] X=alpha Cos[u] Cos[v] (a1NegInner−rNegPosInner Cos[v]+rNegPosInner Cos[3 v]);

[0040] Y=alpha Cos[u] Sin[v] (a2PosInner+rNegPosInner Sin[v]+rNegPosInner Sin[3 v]);

[0041] Z=a3PosInner alpha Sin[u];

[0042] NegNegPosInner:

[0043] X=alpha Cos[u] Cos[v] (a1NegInner−rNegNegInner Cos[v]+rNegNegInner Cos[3 v]);

[0044] Y=−(alpha Cos[u] Sin[v] (−a2NegInner+rNegNegInner Sin[v]+rNegNegInner Sin[3 v]));

[0045] Z=a3PosInner alpha Sin[u];

[0046] PosNegPosInner:

[0047] X=−(alpha Cos[u] Cos[v] (−a1PosInner−rPosNegInner Cos[v]+rPosNegInner Cos[3 v]));

[0048] Y=−(alpha Cos[u] Sin[v] (−a2NegInner+rPosNegInner Sin[v]+rPosNegInner Sin[3 v]));

[0049] Z=a3PosInner alpha Sin[u];

[0050] PosPosNegInner:

[0051] X=alpha Cos[u] (1+alpha tapxPosInner Sin[u]) (a1PosInner Cos[v]+rPosPosInner Sin[2 v]);

[0052] Y=alpha Cos[u] (1+alpha tapyPosInner Sin[u]) Sin[v] (a2PosInner+rPosPosInner Sin[v]+rPosPosInner Sin[3 v]);

[0053] Z=a3NegInner alpha Sin[u];

[0054] NegPosNegInner:

[0055] X=alpha Cos[u] Cos[v] (a1NegInner−rNegPosInner Cos[v]+rNegPosInner Cos[3 v]) (1+alpha tapxNegInner Sin[u]);

[0056] Y=alpha Cos[u] (1+alpha tapyNegInner Sin[u]) Sin[v] (a2PosInner+rNegPosInner Sin[v]+rNegPosInner Sin[3 v]);

[0057] Z=a3NegInner alpha Sin[u];

[0058] NegNegNegInner:

[0059] X=alpha Cos[u] Cos[v] (a1NegInner−rNegNegInner Cos[v]+rNegNegInner Cos[3 v]) (1+alpha tapxNegInner Sin[u]);

[0060] Y=−(alpha Cos[u] (1+alpha tapyNegInner Sin[u]) Sin[v] (−a2NegInner+rNegNegInner Sin[v]+rNegNegInner Sin[3 v]));

[0061] Z=a3NegInner alpha Sin[u];

[0062] PosNegNegInner:

[0063] X=−(alpha Cos[u] Cos[v] (−a1PosInner−rPosNegInner Cos[v]+rPosNegInner Cos[3 v]) (1+alpha tapxPosInner Sin[u]));

[0064] Y=−(alpha Cos[u] (1+alpha tapyPosinner Sin[u]) Sin[v] (−a2NegInner+rPosNegInner Sin[v]+rPosNegInner Sin[3 v]));

[0065] Z=a3NegInner alpha Sin[u];

[0066] PosPosPosOuter:

[0067] 2n X=alpha

[0068] Cos[u] (a1PosOuter Cos[v]+rPosPosOuter Sin[2 v]);

[0069] Y=alpha Cos[u] Sin[v](a2PosOuter+rPosPosOuter Sin[v]+rPosPosOuter Sin[3 v]);

[0070] Z=a3PosOuter alpha Sin[u];

[0071] NegPosPosOuter:

[0072] X=alpha Cos[u] Cos[v] (a1NegOuter−rNegPosOuter Cos[v]+rNegPosOuter Cos[3 v]);

[0073] Y=alpha Cos[u] Sin[v] (a2PosOuter+rNegPosOuter Sin[v]+rNegPosOuter Sin[3 v]);

[0074] Z=a3PosOuter alpha Sin[u];

[0075] NegNegPosOuter:

[0076] X=alpha Cos[u] Cos[v] (a1NegOuter−rNegNegOuter Cos[v]+rNegNegOuter Cos[3 v]);

[0077] Y=−(alpha Cos[u] Sin[v] (−a2NegOuter+rNegNegOuter Sin[v]+rNegNegOuter Sin[3 v]));

[0078] Z=a3PosOuter alpha Sin[u];

[0079] PosNegPosOuter:

[0080] X=−(alpha Cos[u] Cos[v] (−a1PosOuter−rPosNegOuter Cos[v]+rPosNegOuter Cos[3 v]));

[0081] Y=−(alpha Cos[u] Sin[v] (−a2NegOuter+rPosNegOuter Sin[v] +rPosNegOuter Sin[3 v]));

[0082] Z=a3PosOuter alpha Sin[u]; PosPosNegOuter:

[0083] X=alpha Cos[u] (1+alpha tapxPosOuter Sin[u]) (a1PosOuter Cos[v]+rPosPosOuter Sin[2 v] );

[0084] Y=alpha Cos[u] (1+alpha tapyPosOuter Sin[u]) Sin[v] (a2PosOuter+rPosPosOuter Sin[v]+rPosPosOuter Sin[3 v]);

[0085] Z=a3NegOuter alpha Sin[u];

[0086] NegPosNegOuter:

[0087] X=alpha Cos[u] Cos[v] (a1NegOuter−rNegPosOuter Cos[v]+rNegPosOuter Cos[3 v]) (1+alpha tapxNegOuter Sin[u]);

[0088] Y=alpha Cos[u] (1+alpha tapyNegOuter Sin[u]) Sin[v] (a2PosOuter+rNegPosOuter Sin[v]+rNegPosOuter Sin[3 v]);

[0089] Z=a3NegOuter alpha Sin[u];

[0090] NegNegNegOuter:

[0091] X=alpha Cos[u] Cos[v] (a1NegOuter−rNegNegOuter Cos[v]+rNegNegOuter Cos[3 v]) (1+alpha tapxNegOuter Sin[u]);

[0092] Y=−(alpha Cos[u] (1+alpha tapyNegOuter Sin[u]) Sin[v] (−a2NegOuter+rNegNegOuter Sin[v]+rNegNegOuter Sin[3 v]));

[0093] Z=a3NegOuter alpha Sin[u];

[0094] PosNegNegOuter:

[0095] X=−(alpha Cos[u] Cos[v] (−a1PosOuter−rPosNegOuter Cos[v]+rPosNegOuter Cos[3 v]) (1+alpha tapxPosOuter Sin[u]));

[0096] Y=−(alpha Cos[u] (1+alpha tapyPosOuter Sin[u]) Sin[v] (−a2NegOuter+rPosNegOuter Sin[v]+rPosNegOuter Sin[3 v]));

[0097] Z=a3NegOuter alpha Sin[u];

Stage 1

User Initialization

[0098] In this stage, the user provides an approximation to the object of interest by approximately fitting the model to the data. The method for doing this depends on the topology of the object and the form of the parametric component. Here, a method for a left ventricle is described. Note that the approximation may be quite rough as the segmentation routine (Stage 2) will discover the exact object boundaries. In order for the segmentation routine to be effective, the initial model must be reasonably close, within 10 pixels for example.

[0099] The user defines points on the approximate border of the both the inner and outer (endocardial and epicardial) walls of the left ventricle for a cross-section. Then, in order to specify the general bounds of the shape, the user repeats this action on an orthogonal image. For the LV, these would be the short axis and long axis images. ( FIG. 6 ).

Stage 2

2D Optimal Active Contour Segmentation

[0100] The task of refining the contours is also known as segmentation. Here the segmentation module is described. It is given two approximate contours (endocardium and epicardium) from the previous step. These contours, however, do not exactly delineate the features in the image. The present invention seeks to refine the contours so that the endocardium exactly delineates the border between the myocardium and the blood pool, and the epicardium exactly delineates the border between the myocardium and the outside. See FIG. 7 .

[0101] FIGS. 8 and 9 also describe refining the contours following intersection of the model with the image from the approximate contour ( 150 ) to the recovered contour ( 160 ) as illustrated in FIG. 9 .

[0102] The way to do this is to use active contours (snakes). Snakes were first proposed by Kass et al. (M. Kass, A. Witkin, and D. Terzopoulos, “Snakes: Active contour models”, International Journal of Computer Vision, 2:321-331, 1988) are defined as contours that are pushed or pulled towards image features by constraining forces. An energy function based on the smoothness and curvature of the contours and the edge strength in the image along the contour. This energy function is to be minimized. Kass et al. Proposed to minimize the energy function using gradient descent. Amini et al. (A. A. Amini, T. E. Wemouth, and R. C. Jain, “Using dynamic programming for solving variational problems in vision”, IEEE Trans. Pattern Analysis and Machine Intelligence, 12(9):855-867, 1990) used dynamic programming.

[0103] The present invention uses an algorithm similar to Geiger et al. (D. Geiger, A. Gupta, L. A. Costa, and J. Vlontzos, “Dynamic programming for detecting, tracking, and matching deformable contours”, IEEE Trans. Pattern Analysis and Machine Intelligence, 17(3):294-302, 1995), but using Dijktra's graph search algorithm as in Mortensen and Barrett (E. N. Mortensen and W. A. Barrett, “Interactive segmentation with intelligent scissors”, Graphical Models and Image Processing, 60:349-384, 1998).

[0104] Dijkstra's algorithm is a standard graph theory algorithm and can be found in any algorithm's book such as (G. Brassard and P. Bratley, Algorithmics: Theory & Practice, Prentice Hall, 1988). The present system involves creating a graph by defining a node for every pixel in the image. Nodes that correspond to neighboring pixels are connected by a link whose weight is inversely proportional to the gradient in the image between the corresponding pixels. Then, given a source node, Dijkstra's algorithm is able to find the minimum path between the source node and all nodes in a set of sink nodes. Since the cost of a path is the sum of the cost of the links it goes through, the minimum path corresponds to the contour with maximum gradient. If one is not sure on the starting point, then one can define a set of source nodes and connect all of them to a “pseudo” source node with a zero cost link.

[0105] Reference is made to the illustrative steps of the algorithm with an example (see FIG. 7 ). Block 200 shows an image of the left ventricle in a short axis cross section and an approximate contour for the endocardium.

[0106] In the first step the method of the present system defines a search space of a given width (say 10 pixels) around the approximate contour 202 ( a ). Then, a contiguous set of nodes is defined as the source nodes. Since the contour is closed, it does not really matter which nodes are defined as source nodes as long as they cut the search space like a sector. Their neighbors on one side are also defined as the sink nodes (see block 202 ).

[0107] The present system runs Dijkstra's algorithm to find the minimum cost path between the “pseudo” source node and the sink nodes. This generates a contour, but nothing in the algorithm guarantees that it is closed, and most of the time, it is not (see block 204 ).

[0108] So another set of source and sink points in the same search space is defined. The source node is uniquely defined as the mid-point of the previously recovered contour. The sink nodes are its neighbors on one side (see block 206 ) and Dijkstra's algorithm is called again. The recovered contour is guaranteed to be closed since the source and sink nodes are defined to be connected (see block 208 ).

[0109] In the case of the heart left ventricle, both contours are segmented independently. For the endocardium, the cost of a link between two nodes is defined as: 1 $e_1\left(p_1,p_2\right)=\{\begin{array}{cc}\frac{1}{{\nabla I\left(p_2\right)}^2+\varepsilon }& \mathrm{if}z\ge 0\\ 1/\varepsilon & \mathrm{otherwise}\end{array}\mathrm{where}\left(\begin{array}{c}0\\ 0\\ z\end{array}\right)=\left(\begin{array}{c}{x}_2-x_1\\ y_2-y_1\\ 0\end{array}\right)\times \left(\begin{array}{c}\cos \left(\overrightarrow{\nabla }I\left(p_2\right)\right)\\ \sin \left(\overrightarrow{\nabla }I\left(p_2\right)\right)\\ 0\end{array}\right)$

[0110] For the epicardium, it is defined as: 2 $e_2\left(p_1,p_2\right)=\frac{1}{{\nabla I\left(p_2\right)}^2+\varepsilon }$

[0111] Here, ∥∇I∥ is the gradient magnitude and {right arrow over (∇)}I is the gradient direction in the image, ε is a small constant to bound the energy function. In this case, the contour is built clockwise by the Dijkstra process and the image gradient points from bright to dark. By testing whether the sign of the z component of the cross product between the image gradient and the contour direction is positive, and setting the energy to a large number otherwise, the contour is forced to separate a bright region inside from a dark region outside.

Stage 3

3D Hybrid Model Fit

[0112] The initial approximation is fit to the object structure provided by the user in Stage 1 to the sampled segmented data garnered from Stage 2 ( FIG. 10, 170 ). The method employed is an iterative closest point method based on the physically-motivated recovery paradigm. In this paradigm, data points are thought of as exerting forces attracting the model towards them. A hybrid model is capable of rigid, parametric and local (spline-like) deformations. Local deformations may include a deformation penalty which serves to smooth the data. Finite element model (FEM) shape functions are used to distribute the data forces from element surfaces to element nodes. Both the data forces and smoothing forces are scaled, and the data forces are combined with the smoothing forces. The rigid and parametric fitting is discontinued once the local deformations are initiated. The result is shown in FIG. 10, 180 .

[0113] Typically, in hybrid model recovery rigid body characteristics (translation and rotation) are adjusted until they settle. Then, rigid body characteristics and the parametric component are fit. And, once they settle, the rigid body characteristics, the parametric component and the local component are fit. The present method is different from this standard scheme by freezing the rigid body characteristics and parametric component during the local model fit. The rationale is that if the parametric component is to provide a gross description of the data it should arrive at the same result irrespective of the presence of a local component. In addition, freezing the rigid body characteristics and parametric component during the local model fit, speeds the fitting process significantly.

Experimental Results

[0114] The technique was tested on a volunteer data set with two phases of the cardiac cycle, En-diastole (ED) and en-systole (ES). The blood pool volume for ED was 271 ml with a myocardial volume of 275 ml. For ES, the blood volume was 231 ml with a myocardial volume of 285 ml. The ejection fraction for the volunteer was 14%.

[0115] While certain embodiments have been shown and described, it is distinctly understood that the invention is not limited thereto but may be otherwise embodied within the scope of the appended claims.