Curvature based system for the segmentation and analysis of cardiac magnetic resonance images

[0002] The present invention relates generally to the field of imaging, and image processing and analysis techniques, such as those used in medical diagnostics and other fields. More particularly, the invention relates to a technique for rapidly segmenting and identifying features of interest in discrete pixel images and processing data extracted from the images following identification of the features.

[0003] A great number of applications have arisen in recent years employing discrete pixel images of various types. In general, such images are comprised of an array of picture elements organized in rows and columns of an image matrix. Each element in the matrix may be characterized by various parameters, such as intensity, color, shade, and so forth. These characteristics are most often associated with or defined by information expressed in digital terms and stored in data files defining the image. By organizing and displaying the matrix of picture elements with reference to the parameters, a recognizable and useful reconstructed image may be obtained. The resulting image is typically displayed on a computer monitor or printed on a tangible medium such as paper or photographic film.

[0004] A number of modalities are known for acquiring discrete pixel images, particularly in the medical diagnostics field. For example, magnetic resonance imaging (MRI) systems have been developed which permit high-quality images of anatomies of interest in a subject. Such systems make use of controllable magnetic fields which affect precession of gyromagnetic material within the subject. By exciting the gyromagnetic material in the presence of specific fields, the materials can be made to emit signals which are captured and processed. The signals can be transformed to identify the particular intensity of the emissions at specific locations within the subject. By further processing the digitized signals corresponding to these intensities, image data is obtained which can be analyzed, stored, transmitted, and further processed to create a reconstructed image. Other modalities in the medical diagnostics field include computed tomography (CT) systems, digital x-ray systems, positron emission tomography (PET) systems, ultrasound systems, and so forth. In modem systems of all of these types, data sets may be obtained directly, or indirectly, which define features of interest within a subject patient for image reconstruction and processing.

[0005] Many applications exist for analysis of identifiable structures within discrete pixel images. In conventional images, the features are typically recognizable by contrasts, outlines, expected shapes, and other somewhat intuitive boundaries or characteristics recognizable to the experienced viewer. Depending upon the image quality, a qualified viewer may identify these features and even take measurements of the structures. However, conventional techniques for manually doing so are quite cumbersome, time-consuming, and imprecise. In the medical diagnostics field, for example, measurements of internal anatomies such as structures of the heart, brain, and other organs and tissues, have conventionally been made by hand or through software requiring substantial user interaction. Due to the drawbacks of these techniques, however, attempts have been made to automate the process.

[0006] Techniques which have been advanced for automating feature identification and measurement are not entirely satisfactory. A key aspect of this is the “segmentation problem” which refers to the identification of homogeneous regions in a given image, defined by a set of edges or boundaries. Processes have been developed that use mathematical constructs to deform candidate boundaries to correspond to the limits of a feature of interest. Such deformable boundaries are sometimes referred to as “snakes.” However, the methods employed for segmentation based upon these techniques are particularly susceptible to error due to noise present in the image data. Such noise, which may be visible as bright or dark spots or regions on the reconstructed image, may cause the candidate boundary to stop its expansion or contraction before reaching the limits of the feature, or may cause other anomalies in the result. While such noise may be eliminated or reduced by smoothing and other filtering techniques, these also tend to reduce the resolution of the feature of interest, thereby reducing the likelihood that the snake will converge accurately on the feature.

[0007] Another serious drawback of heretofore known techniques for automatic segmentation via snakes is the extremely demanding computational requirements involved. In particular, known techniques adopt algorithms which require sequential calculation of a very large number of parameters as the boundary expands or contracts towards the feature of interest. As a result, to be carried out in a reasonable amount of time, very sophisticated and powerful computers are required, particularly for noisy images, and those including a large number of picture elements.

[0008] Further drawbacks in existing techniques for automatic segmentation are rooted in the particular algorithms used to generate the candidate boundary and to deform the boundary toward the feature of interest. In particular, the algorithms typically include mechanisms for converging and stopping the evolution of the candidate boundary. Depending upon the techniques employed, the boundary may not fully conform to the feature of interest, particularly where the feature includes concavities and similar contours. Also, where adjacent features are present, the boundaries may not recognize distinctions between the adjacent features and, again, fail to conform accurately to the feature of interest. In either case, the resulting analysis is unreliable due to the failure to accurately recognize the boundary (and hence the surface area, volume or other parameter) of the feature of interest.

[0009] There is a need, therefore, for an improved technique for performing segmentation on discrete pixel images. There is, at present, an increasing need for a technique which will permit such segmentation on a series of images, such as images acquired over a three-dimensional space or over time, or both. Such techniques would be particularly useful in analyzing moving tissues, such as those of the heart. The technique would also find useful applications in still structures such as the brain, and so forth.

[0010] The present invention provides a technique for segmenting features of interest in discrete pixel images designed to respond to these needs. The technique makes use of a computationally efficient algorithm for initiating and expanding or contracting a candidate boundary for the feature. In accordance with certain aspects of the technique, a candidate boundary or snake is initialized and evolved by reference to mean curvature-weighted normals. Gradients of intensities within the image may be calculated before this evolution to greatly simplify the computational demands during the evolution. The candidate boundary thus converges extremely rapidly toward the feature and stops upon reaching the feature. The technique is particularly robust and tolerant of noise within the image. However, a diffused or smooth gradients may be used during the evolution to improve the tolerance to noise.

[0011] The technique may be used to segment both static and dynamic (or changing) structures. In both cases, the approach may be applied to a series of data sets representing individual images or slices of a subject. Over space, the technique may thus be used to reconstruct a three-dimensional representation of the feature. Over time, the technique may be used to analyze changes in the feature.

[0012] In a presently preferred implementation, the segmentation technique is used to measure and analyze features of the heart, such as the endocardium and the myocardium. By time analysis of this information, important characteristics of the heart may be evaluated, such as the ejection fraction. In other structures, both two-dimensional and three-dimensional, and time dependent or independent measurements may be made.

[0013] FIG. 1 is a diagrammatical representation of an imaging system for acquiring image data and for performing segmentation and analysis on acquired data in accordance with certain aspects of the present technique;

[0014] FIG. 2 is a representation of a heart and a series of images produced of the heart through use of a system such as that illustrated in FIG. 1;

[0015] FIG. 3 is a single discrete pixel image of the heart produced as illustrated in FIG. 2, illustrating the initialization of a pair of candidate boundaries for endocardium and myocardium structures;

[0016] FIG. 4 is an image of the heart as shown in FIG. 3 following partial expansion and contraction of the candidate boundaries initialized in FIG. 3;

[0017] FIG. 5 is an image of the heart as shown in FIG. 3 following full expansion and contraction of the candidate boundaries to the relevant features of interest;

[0018] FIG. 6 is a detailed representation of a portion of a candidate boundary expanding toward a feature of interest and converging on the feature of interest in a rapid manner by reference to mean curvature weighted normals;

[0019] FIG. 7 is an enlarged view of a portion of a discrete pixel image illustrating the portion of the analysis performed for expansion of a candidate boundary in accordance with aspects of the present technique;

[0020] FIG. 8 is a diagram illustrating steps in exemplary control logic for performing a segmentation analysis in accordance with the present technique on a discrete pixel image;

[0021] FIG. 9 is a time-based reconstruction of a series of measurements of an endocardium surface in accordance with the present technique for measurement and analysis of parameters of a patient's heart;

[0022] FIG. 10 is a perspective view of a reconstructed three-dimensional feature, such as the left ventricle of a heart, based upon analysis in accordance with the present technique;

[0023] FIG. 11 is a graphical representation of typical post-segmentation analysis which may be performed in accordance with the present technique, such as based upon the time-based reconstruction of FIG. 9;

[0024] FIG. 12 is a representation of a discrete pixel image of a brain obtained through an imaging system such as an MRI system; and,

[0025] FIG. 13 is a reconstructed three-dimensional image of the brain illustrated in FIG. 12, following segmentation analysis and reconstruction in accordance with aspects of the present technique.

[0026] Turning now to the drawings, and referring first to FIG. 1, an imaging and image analysis system 10 is illustrated diagrammatically as including a scanner 12, scanner control circuitry 14, and system control circuitry 16. It should be noted that while system 10 is based upon a magnetic resonance imaging (MRI) data acquisition modality in the illustrated and described embodiment, the segmentation technique provided herein is not limited to any particular acquisition modality. Accordingly, the technique may be applied to image data acquired through computed tomography systems, positron emission tomography systems, ultrasound systems, x-ray systems, photographic systems, and so forth. Moreover, while the segmentation and analysis functions may be performed in system control circuitry 16 in the illustrated embodiment, such analysis may be performed in remote or stand-alone computers, including work stations, conventional commercial computers, laptop computers, and so forth, either directly connected to the data acquisition circuitry or remote from the data acquisition circuitry. Accordingly, the segmentation and analysis functions may be performed based on archived image data stored or transferred via networks, and the like.

[0027] In system 10, as illustrated in FIG. 1, scanner 12 has a housing 18 in which a patient bore 20 is formed. Table 22 is positionable within patient bore 20 for supporting a patient 24 to be imaged. While specific systems may include larger or smaller bores, in the illustrated embodiment, scanner 12 is a full-body scanner in which a patient may be fully positioned.

[0028] Within housing 18, scanner 12 includes a primary coil 26 which generates a highly uniform magnetic field generally aligned with the patient. A series of gradient coils 28 are provided for generating controlled, pulsed magnetic fields in accordance with pre-established imaging protocols. A radio frequency coil 30 is provided for generating radio frequency pulses used to excite gyromagnetic material within the patient. In the illustrated embodiment, coil 30 also serves as a receiving coil. Thus, coil 30 may be coupled with driving and receiving circuitry in passive and active modes for receiving emissions from the gyromagnetic material and for outputting radio frequency excitation pulses, respectively. Alternatively, various configurations of receiving coils may be provided separate from coil 30, such as coils specifically adapted for target anatomies.

[0029] The coils of scanner 12 are controlled by external circuitry for generating the desired magnetic fields and radio frequency pulses in accordance with the desired imaging protocol. Thus, a power source 32 is coupled to the scanner to supply power for generation of the desired fields and signals. In addition, one of more sensors may be provided for triggering the beginning or the end of a desired imaging sequence. For example, as illustrated in FIG. 1, a cardiac monitor 34 may be provided for sensing heart beats of the patient which are used to trigger the beginning of an acquisition sequence in which a series of images are acquired through at least one cardiac cycle.

[0030] The scanner control circuitry 14 includes functional components designed to interface with the scanner, to control the pulse sequences implemented for the data acquisition, and to receive and process image data. Thus, as illustrated in FIG. 1, scanner control circuitry 14 includes a control circuit 36 which commands operation of the scanner through the intermediary of an amplifier/control circuit 38 and a transmit/receive control circuit 40. Amplifier/control circuit 38 is coupled to gradient coils 28 to supply control signals which define pulsed field gradients produced by the coils during an examination sequence. Transmit/receive circuit 40 includes amplification circuitry for producing the radio frequency signals via coil 30, and for switching to passive modes for receiving detected emission signals from the gyromagnetic material. These signals are applied to control circuit 36 which may perform signal processing, filtration, scaling, and similar functions. Memory circuitry 42 is provided for storing configuration parameters, pulse sequence descriptions, and so forth. These are applied to or accessed by control circuit 36 during particular examination sequences. Interface circuitry 44 is coupled to control circuit 36 to permit the scanner control circuitry 14 to receive configuration and examination commands from system control circuitry 16 or directly from an operator interface 46, such as an operator station, work station, and so forth. Interface circuitry 44 also serves to transmit raw or processed image data from scanner control circuitry 14 to system control circuitry 16.

[0031] System control circuitry 16 includes various functional circuitry which may be based in a general purpose or application-specific computer or work station. As will be appreciated by those skilled in the art, such stations will include a central processing unit (CPU) or signal processor designed to execute programming code. Interface circuitry 48, which may be at least partially defined by such code, permits commands to be transmitted from the system control circuitry to scanner control circuitry 14, and to receive raw or processed image data back from the scanner. Among the functions executed by the system control circuitry 16, segmentation and analysis functions, as described below, are performed by a segmentation and analysis module 50. Module 50, which may also be at least partially defined by programming code executed by the control circuitry serves to analyze raw or processed image data to identify features of interest from the data and to process the data to evaluate the size, configuration, and other parameters of the particular features of interest. Parameters required for execution of the code implemented by control circuitry 16 are stored in a memory circuit 52. Memory circuit 52 may include various types of solid state or dynamic memory, such as random access memory, tape drives, optical storage devices, and so forth. An operator interface 54 is provided in circuitry 16 and is coupled to the segmentation and analysis module 50 for commanding and controlling the segmentation and analysis functions performed. A remote interface 56 may be provided to permit configuration and image data to be exchanged between circuitry 16 and remote stations, such as picture archive and communication systems as indicated at reference numeral 58. Moreover, systems 58 may include teleradiology stations to which image data is transferred. Each such remote station may include its own segmentation and analysis modules for performing the functions described below on image data.

[0032] As will be appreciated by those skilled in the art, when the present technique is employed in systems other than MRI diagnostic systems, the particular circuitry involved will be adapted accordingly. For example, in x-ray, computed tomography or ultrasound systems, data acquisition circuitry will be adapted to acquire image data sets in accordance with the particular physics of those systems. In general, however, the data sets will include a series of values descriptive of pixel characteristics in an image matrix. This image data is analyzed as described below to identify particular features of interest either in an automated or semi-automated fashion.

[0033] The present technique is particularly well suited to analysis of anatomical features of interest. These features may be included in both static and dynamic tissues of a patient, and may be defined in three dimensions as well as in time. FIG. 2 illustrates an exemplary sequence of images formed of a heart via an MRI system such as that illustrated in FIG. 1. As will be appreciated by those skilled in the art, various pulse sequences may be employed for generating selected slice images of such tissues. For example, in the embodiment illustrated in FIG. 2, image sets 60 are acquired of features within a heart 62 of a patient. In this embodiment, features of interest include the boundaries of the left ventricle which are used to generate measurements of particular interest such as the ejection fraction, volume at various stages of the cardiac cycle, and so forth.

[0034] The particular features measured in this exemplary process are the endocardium and the myocardium. To permit these features to be imaged and measured through the technique described below, a series of images are acquired in the image set 60 in a short axis view of the heart. In these views, the magnetic fields generated by gradients coils 28 are physically oriented to produce a gradient system aligned along an axis 64 perpendicular to which a series of image slices 66 are oriented. The left ventricle 68 of the heart appears in each image generated along these slices. The endocardium 70 is visible in these images as is the myocardium 72, each having boundaries which vary over the sequence of images along axis 64, as well as over the cardiac cycle as the heart shape and configuration change. Thus, the image data includes data defining a series of images 74 for each point in time through the cardiac cycle. Each time sequence of images includes a series of images 74A, 74B, . . . , 74Y, and 74Z corresponding to the individual slices 66 in the acquired image set. At various points in the cardiac cycle, then, similar slice images are acquired to define subsequent image sequences 76, each comprising a series of slice images 76A, 76B, . . . , 76Y, and 76Z similar to those of set 74. In each of the images, the outline 78 of the heart is visible by contrast between pixels defining the image. Features 80 of interest, such as the left and right ventricles, are also visible in the image data and are subsequently identified as described below.

[0035] It should be noted that the particular orientation, composition, number and sequence of images may be adapted for any desired feature of interest. In particular, in a presently preferred configuration, images acquired for a cardiac cycle comprise 14 separate slices along axis 64 and 20 successive imaging times. The resulting image set 60 thus includes 280 images, each defined by a unique matrix of pixel data. Similarly, the matrix size may be adapted for any suitable dimensions or resolution. In the illustrated embodiment, images made up of 256×256 pixels are employed.

[0036] The present technique is based upon a mathematical construct in which one or more candidate boundaries are established in an image data set and evolved (i.e. expanded, contracted or conformed) toward a feature of interest. The feature may be defined by any one of a variety of characteristics, such as pixel intensity values. The characteristic used to define the feature provide a transition which sufficiently identifies the feature to permit the evolving boundary to rapidly converge on and stop at the feature. Subsequent data processing can then be performed based upon the configuration of the evolved boundary, such as by determining the number of pixels within the boundary, the boundary length, and so forth.

[0037] FIGS. 3, 4 and 5 illustrate the initiation and evolution of a set of candidate boundaries for identifying endocardium and myocardium boundaries in a cardiac image of the type illustrated in FIG. 2. As shown in FIG. 3, any one of the images in the image set may be employed for evolution of a candidate boundary to identify the features of interest. In the exemplary image 82 shown in FIG. 3, a series of tissue boundaries 84 appear, some or all of which define volumes of particular interest. As will be appreciated by those skilled in the art, the particular location of these boundaries is often less clear than illustrated in the Figures, represented by greater or lesser degrees of contrast between adjacent pixels in an image matrix. Moreover, various other features which are not of interest, some of which may include image artifacts, may be included in the image. The present technique permits the evolution of candidate boundaries to locate the boundaries of the features of interest with relatively high tolerance to such image artifacts and noise.

[0038] As illustrated in FIG. 3, the process, described in greater detail below, begins with initiation of one or more candidate boundaries or snakes as shown at reference numeral 86 for a candidate myocardium boundary and at reference numeral 88 for a candidate endocardium boundary. As shown in FIG. 4, through the present technique, the candidate boundaries are evolved in the direction of the features of interest as indicated by arrows 90. In particular, the candidate myocardium boundary is expanded toward the limits of the myocardium, both inwardly and outwardly, while the candidate endocardium boundary 88 is expanded towards the interior limits of the left ventricle. As shown in FIG. 5, following full evolution of the candidate boundaries, the limits of the features of interest are identified and further evolution of the candidate boundaries stops.

[0039] As will be appreciated by those skilled in the art, the present technique is based upon the grouping of parts of an image into units which are homogeneous with respect to one or more features. The candidate boundary, which may also be referred to as an active contour, is permitted to evolve, once initiated, under the influence of image-dependent factors, and certain parameters which may be set by the user. In the present technique, for example, rapid evolution is provided by computation of mean curvatures and mean curvature-weighted normals as described below, in the discrete domain. Moreover, the computations involved in the evolution may be greatly simplified by a single computation of gradients prior to the evolution of the candidate boundary. Where desired, additional tolerance to noise is provided by smoothing or diffusion of these gradients. However, the original image is not smoothed, but rather, the gradient direction along which the candidate boundary is evolved is diffused, preserving original image features and resolution.

[0040] The continuous domain implementation of the technique will now be described, followed by a description of the discrete domain implementation. The technique is based upon identification of a characteristic, denoted herein “φ,”which is an image-based function considered to be small for the feature of interest. For example, φ may represent intensity values or the inverse of intensity values for individual pixels, depending upon the nature of the image data and its intended presentation. The feature, which is identifiable in any one of a series of spatially and temporally defined images, as illustrated in FIG. 2, is identified by reference to the relationship (eq. 1): $\frac{\partial S}{\partial t}=\phi f\left(H\right)N+g\left(\phi ,\nabla \phi \right)$

[0041] where S(x₁, . . . , x_n, t) is the evolving candidate boundary, ∂S/∂t represents the partial derivative of S with respect to t, φ is a stopping term defined below, H represents a mean curvature of the moving surface, N is the unit normal, and 𝑓 and g are functions which arise from the minimization of an energy functional derived from weighted length, area and volume. It should be noted that equation 1 is valid in any number of space and time dimensions. In the present embodiment, four dimensions may be employed, including three space dimensions and a time dimension. Alternatively, two-dimensional implementations may be employed in which only length and width spatial dimensions are used, and with or without an additional time dimension.

[0042] As a first step in the analysis, φ is defined to be a positive image-dependent function which is small near the boundary of the feature of interest. In this implementation, the function φ (x₁, . . . , x_n) is employed as a stopping term in the evolution of the candidate boundary. For example, the term φ (x₁, . . . , x_n) may be chosen to be small near an edge, and acts to stop the evolution of the candidate boundary when the contour approaches or overlies the edge. In this case, the definition of φ may be expressed (eq. 2): $\phi \text{:}=\frac{1}{1+{\Vert \nabla G_{\sigma }*I\Vert }^2}$

[0043] where I is the n-dimensional gray-scale image and G_σ is a Gaussian (smoothing) filter. In general for two vectors a=(a₁, . . . , a_n) and b=(b₁, . . . , b_n), we set <a, b>=a·b=a₁b₁+a₂b₂+ . . . +a_nb_n, and ∥a∥={overscore (<a, a>)}. This stopping term is designed to find intensity-based edges. Others may be formulated for textures and for various features in color imagery as well. The evolution of the candidate boundary is considered to progress through a flow in which the boundary of the feature of interest resides at the bottom of a potential well. The value of the boundary may be expressed (eq. 3): $\frac{\partial S}{\partial t}=\phi \left(H+ϵ\right)N+\left(1\text{/}nS·\nabla \phi -\nabla \phi ·N\right)N$

[0044] where S (x₁, . . . , x_n, t), H, N, and φ have the values defined above, and ε is a binary metric threshold constant, and ∇ is a gradient term of the value φ in the directions of the coordinate system of the image. Here, H is the mean curvature which is the sum of the principal curvatures, and N is the inward unit normal. The direction of evolution of the candidate boundary can be controlled by the parameter ε. As will be appreciated by those skilled in the art, this expression has the interpretation of balancing the volume and surface area minimizing forces, so that the edge in n-dimensions is regarded as a minimal surface/volume in an appropriate conformal geometry. To enhance the tolerance of the technique to image artifacts and other image noise, a diffused gradient term ∇φ may be used as described below.

[0045] In implementing the technique based upon equation 3, the equation is suitably discretized on a triangulated surface so as to move the vertices of the candidate boundary to decrease the energy of which equation 3 is Euler-Lagrange derivative at each step. This technique entails the computation of the mean curvature, the normal, and the Gaussian curvature (which informs the method of the concavity or convexity of the surface at any given point in its evolution).

[0046] The triangulated surface employed in the technique is defined by a series of flat surfaces, edges and points as illustrated in FIG. 6. These basic geometric elements are shown in FIG. 6 as adjacent surfaces 94, joined to one another by straight edges 96, the edges meeting in common vertices 98. In the technique, these elements are referred to by identifiers which also contain an orientation for each element. In particular, vertices 98 are points in space, edges 96 are straight line segments joining pairs of vertices, and surfaces 94 are faces of flat triangles bounded by three edges. The surface topology is given by cross-references among the elements. The identifiers may be defined by four byte words, with one of the bits being an orientation bit. In particular, each face-edge is defined in two double-linked lists, one of face-edges around the edge, and another of face-edges around the face. Each edge points to its two end points, and to one face-edge loop around it. Each face points to one face-edge and the face-edge loop around it, and to the two bounding faces on its two sides. As described below, normals 100 extend from each point or vertex 98. Evolution of the candidate boundary 88 follows the mean curvature weighted normals 100 to evolve rapidly towards the boundaries of the feature of interest indicated at reference numeral 84 in FIG. 6.

[0047] Assuming a three dimensional implementation of the technique (defining a three dimensional surface), a maximum curvature k₁and a minimal curvature k₂may be computed. The mean curvature H is, then, defined by the relationship (eq. 4):

[0048] H=k₁+k₂

[0049] If S:[0,1] X [0,1]→R³(R=real numbers) denotes a compact imbedded surface with local coordinates (u, v), the values S_uand S_vmay be set as (eq. 5A): $S_u=\frac{\partial S}{\partial u}$

[0050] and (eq. 5B) $S_v=\frac{\partial S}{\partial v}$

[0051] The infinitesimal area on S is then given by the relationship (eq. 6):

[0052] dS=(∥S_u∥²∥S_v∥²−<S_u,S_v>²)^½dudv

[0053] It is noted in the present technique that the gradient flow associated with the area functional for such surfaces (ie., the direction in which the area is changing most rapidly) is given by the relationship (eq. 7): $\frac{\partial S}{\partial t}=\mathrm{HN}$

[0054] Thus, the mean curvature vector HN on a triangulated surface may be determined by taking the direction in which the area is changing most rapidly, and, given the length of this vector, the normal direction can be computed. This technique is employed in the present procedure as a basis for computing the mean curvature on a triangulated surface. In particular, if p is a vertex as described above, and each vertex has a star of faces surrounding it of area A_p, the force due to the surface tension on the vertex is given by the relationship (eq. 8): $F_p=-\frac{\partial A_p}{\partial p}$

[0055] Because each face has three such vertices, the area associated with p is A_p/3. Thus, the (average) mean curvature at v is (eq. 9): $h_p=\frac{3F_p}{2A_p}$

[0056] Moreover, the Gaussian curvature K may be defined as (eq. 10):

[0057] K=k₁k₂

[0058] where the values of k₁and k₂are those described above. This value may be computed at each vertex on a triangulated surface by the relationship (eq. 11):

[0059] (Σ_iθ_i)−π

[0060] where the sum is taken over the adjacent angles θ to the given vertex. This provides an indication of whether the surface is convex or concave at a given point (convex being defined as a quantity <π, and concave >π).

[0061] As noted above, the foregoing technique may be employed with diffused gradients to enhance the tolerance of the technique to image artifacts and noise. While the tolerance to noise may be provided by pre-processing smoothing, this may have the undesired effect or removing salient features, including the definition of the feature of interest. In a presently preferred technique, enhanced tolerance to noise is provided by smoothing or diffusing a gradient term ∇φ, which represents the actual vector field. In particular, for the desired vector field (eq. 12):

[0062] v(x₁, . . . ,x_n)=(v₁(x₁, . . . ,x_n), . . . , v_n(x₁, . . . , x_n))

[0063] (representing the desired smooth vector field), v(x₁, . . . , x_n) may be derived as the minimal energy solution to the variational problem (eq. 13):

[0064] ∫∫2αΣ_i=0ⁿ∥∇v_i∥+∥∇φ∥²∥v−∇φ∥²dx₁, . . . dx_n

[0065] This minimization problem may be solved by using gradient descent associated Euler-Lagrange equations (eq. 14): $\alpha D\left(\frac{v_i}{\Vert v_i\Vert }\right)-\left(v_i-\nabla \phi \right){\Vert \nabla \phi \Vert }^2=0,i=1,\dots ,n$

[0066] noting that for a function f: Rⁿ→R, the value div f is set in accordance with the relationship defined by (eq. 15): $\mathrm{Df}=\sum _{i=0}^n\frac{\partial f}{\partial x_i}$

[0067] In the relationship defined above, the value α acts as a scale space parameter, the magnitude of which is directly related to the degree of smoothing effect on the gradient vector field. Moreover, the term (eq. 16): $D\left(\frac{v_i}{\Vert v_i\Vert }\right)$

[0068] defines the mean curvature on the isophotes of v_i, that is, on the sets (eq. 17):

[0069] v_i(x₁, . . . , x_n)=c

[0070] for a constant c.

[0071] The foregoing technique provides a robust and rapid mechanism for evolving a candidate boundary toward a boundary of a feature of interest. To facilitate implementation of the technique with discrete pixel image data, differences between values of image characteristics φ (eg. intensity values) for adjacent pixels (or voxels) are divided by pixel spacing to approximate the derivative of values described above. This technique may be employed in three spatial dimensions for three-dimensional implementation. The same technique may be employed in only two dimensions for two-dimensional implementations as illustrated in FIG. 7. As indicated above, the image is representative of a series of voxels 102 which, in a two-dimensional implementation corresponds to image pixels 104. Characteristics of particular interest in the pixels will define the feature of interest, as indicated at pixels 106 in FIG. 7. In a coordinate system defined by x, y and z axes, 108, 110 and 112, respectively in FIG. 7, the slice may be considered to extend in x and y directions for a constant location along the z direction. In the cardiac implementation described above, a series of such slices are defined by image slices 66, and correspond to each image in the series 74 and 76, as well as in intermediate series. Thus, the image matrix extends in the x direction as indicated at reference numeral 114 in FIG. 7, and in the y direction as indicated at reference numeral 116. Rows 18 and columns 120 of pixels thus define the image.

[0072] As mentioned above, intensity differences may be defined by comparing intensities of adjacent pixels in both the x direction 114 and in the y direction 116 in the discrete pixel image data. Similarly, intensity differences over time may be defined by comparing identical locations of image data in different time images (see image series 74 and 76 in FIG. 2), and an approximation of the time differential determined by dividing these differences by the time step between images. By sequential iteration of the computations described above, the candidate boundary automatically converges on and stops at the desired feature boundary or edge. Again, this boundary may be defined as the location where the term φ is small (i.e. minimal). There is generally no need to provide, a priori, a bound on the number of iterations.

[0073] For smoothing the gradient, and thereby enhancing the tolerance to image noise, a “time” parameter is added to equation 14 above to define the relationship (eq. 18): $\frac{\partial v_i}{\partial t}=\alpha D\left(\frac{v_i}{\Vert v_i\Vert }\right)-\left(v_i-\nabla \phi \right){\Vert \nabla \phi \Vert }^{2}i=1,\dots ,n$

[0074] To simplify the nomenclature, in the discrete embodiment two dimensions in a two-dimensional implementation may be defined as x=x₁, and y=x₂. The two-dimensional image, such as that illustrated in FIG. 7, is discretized with the grid points placed at the pixel locations. The values Δx and Δy represent the step sizes in the two directions, which in a simplified implementation may be assumed to be equal, such that (eq. 19):

[0075] Δx=Δy=h

[0076] The discretization in a temporal direction is introduced and the time step Δt is used to compute the temporal intensity gradient change at each pixel. These discretizations are then employed to compute gradients at each point (iΔx, jΔy), as follows (eq's 20, 21):

[0077] φ_x(i,j)≈0.5(I(i+1,j)−I(i−1,j))/Δx

[0078] φ_y(i,j)≈0.5(I(i,j+1)−I(i,j−1))/Δy

[0079] In this implementation, the characteristic of interest in the pixels is considered to be a pixel intensity. As noted above, however, other characteristics such as color, and so forth may be employed. It should also be noted that in the foregoing relationships, the base relationship (eq. 22):

[0080] φ_x(i,j)=φ_x(iΔx,jΔy)

[0081] is employed, as is a similar base relationship in the y direction.

[0082] As a subsequent step in the discretization of the system, the value of v₁is defined as (eq. 23):

[0083] v_1ij,k=v₁(iΔx,jΔy,kΔt)

[0084] Based upon this relationship, the shock derivative may be defined as (eq. 24): $v_{1i,j,k+1}=v_{1i,j,k}+\frac{\Delta t}{h}\left[{\Delta }_{-}^x\left(\frac{{\Delta }_{+}^x{v}_{1i,j,k}}{{\left({\rho }_{v_1^x}\right)}_{i,j,k}}\right)+{\Delta }_{-}^y\left(\frac{{\Delta }_{+}^y{v}_{1i,j,k}}{{\left({\rho }_{v_1^y}\right)}_{i,j,k}}\right)\right]$

[0085] where (eq. 25):

[0086] (ρ_v^y)_i,j,k=[(Δ₊^xv_1i,j,k)²+(m(Δ₊^yv_1i,j,k,Δ₋^yv_1ij,k))²]^½

[0087] and (eq. 26):

[0088] (ρ_v^x)_i,j,k=[(m(Δ₊^xv_1ij,k,Δ₋^xv_1i,j,k))²+(Δ₊^yv_1i,j,k)²]^½

[0089] In these relationships, the following definitions are held (eq. 27):

[0090] Δ_∓^xv_1i,j,k=∓(v_1i∓i,j,k−v_1i,j,k)

[0091] and (eq. 28):

[0092] Δ_∓^yv_1i,j,k=∓(v_1i∓i,j,k−v_1i,j,k)

[0093] Moreover, the function m is defined as (eq. 29): $m\left(a,b\right)=\left(\frac{\mathrm{sgn}\left(a\right)+\mathrm{sgn}\left(b\right)}{2}\right)\min \left(\left|a\right|,\left|b\right|\right)$

[0094] In a preferred implementation, a current Courant-Friedrichs-Lewy (CFL) step size restriction is imposed which may take the following form (eq. 30): $\frac{\Delta t}{h^2}\le c$

[0095] The technique may be expressed as a sequence of logical steps as defined in FIG. 8. These steps are preferably defined by appropriate programming code executed by module 50 (see FIG. 1). In the steps of the sequence, indicated generally by reference numeral 130, image data is first collected as indicated at step 132. Again, this image data may be collected by various means, such as through MRI systems, CT systems, x-ray systems, ultrasound systems, and so forth. The values of interest, φ, as indicated at step 134, are then computed for the image data, for each pixel in the discrete implementation. Again, these values may comprise any suitable image characteristic, typically intensity. For analysis of a series of images, distributed over space and/or time, characteristic values are computed for each of the space and/or time images. As indicated at step 136, these values φ are then stored, such as in the system's digital memory.

[0096] At step 138, the gradients of the values φ are computed. As describe above, in the discrete implementation, these gradients are based upon differences between adjacent pixel values divided by the spacing between the pixels. For time gradients, identically located pixels are compared and differences between the pixels are divided by the sample time step. At step 140 the gradients are diffused as described above, and the diffused gradients are stored at step 142. It will be noted by those skilled in the art, that the present technique provides for calculation of the gradients only one time at the initial phase of analysis. Subsequently, evolution of the candidate boundary proceeds without the need to again calculate the gradients, significantly improving the performance and speed of the technique, greatly enhancing the computational efficiency, and reducing the need for sophisticated or expensive processing capabilities.

[0097] At step 144 the candidate boundary is initialized. In a presently preferred technique, a user may trace a candidate boundary on an image at a location within a feature of interest for evolution within the feature, outside a feature of interest for evolution from outside inwardly towards the feature, or between various features. Where multiple images are analyzed, a single initialization may be used for the series of images and the initial candidate boundary referenced in the analysis of the other images in the sequence. Alternatively, this boundary initialization may be fully automated as well.

[0098] At step 146, the mean curvatures H are computed, and at step 148 the mean curvature-weighted normals HN are computed. In the discrete implementation, steps 146 and 148 are performed inherently through the discretizations summarized in equations 20-30. At step 140 the candidate boundary is evolved through a series of iterative steps which converge on and automatically stop at the feature boundary. At step 152, any post feature identification processing may be performed, such as computation of areas, volumes, and so forth, based upon the final location of the fully evolved boundary.

[0099] Post feature identification processing may take a number of forms, including calculation of scalar and vector quantities, and reconstruction of two and three-dimensional images. Moreover, images may be reconstructed based upon several dimensions, including time. By way of example, FIG. 9 represents reconstruction of a time profile for a specific location or slice of an endocardium, comprising a series of filly evolved boundaries generated in accordance with the foregoing technique. The reconstructed profile, designated generally by the reference numeral 160, thus comprises a first boundary 162 followed by a subsequent boundary 164 which will typically correspond to the same spatial slice in the sequence of cardiac images, but subsequent in time. By reconstructing the profile based on a series of such evolved boundaries, a composite profile 166 may be obtained which provides a graphical representation of the change in configuration of the feature over time. As indicated below, by measuring parameters of the boundary, such as by determining the number of pixels within each evolved boundary, specific measurements and comparisons over time may be obtained. Similar time reconstructions may be generated based upon any feature of interest which can be identified and measured by the foregoing technique, including boundaries of the endocardium, the myocardium, as well as other static and moving boundaries.

[0100] FIG. 10 represents a spatial reconstruction 170 based upon a series of filly evolved candidate boundaries from a series of spatially related images. In particular, FIG. 10 represents a three-dimensional profile of a left ventricle defined by a sequence of boundaries 172, 174, and so forth, reconstructed to provide a fill three-dimensional profile 176. As will be appreciated by those skilled in the art, the profile permits a composite volume of the feature to be measured by summation of the individual volumes of the evolved boundaries. Such composite images may be presented as time reconstructions as well, to provide simulated two or three-dimensional motion reconstructions offering an intuitive illustration of changes in the feature over time.

[0101] Further processing may be provided to identify particular parameters of interest both in space, time, or both. FIG. 11 represents a graphical presentation of individual volumes of an endocardium as computed for specific times in a cardiac cycle in accordance with the foregoing technique. The resulting time profile 180 thus represents a change in volume over time, with each volume being calculated by summation of individual slices of the feature as described above. In the example illustrated in FIG. 11, a minimal volume 182 may be identified by analysis of the sequence of resulting volumes. Similarly, a maximum volume 184 may be identified. It should be noted that the volumes may be expressed or illustrated in terms of reconstructed images as shown in FIG. 10, as well as quantities representative of overall or actual volumes measured. Between the minimum and maximum volumes, a series of volume steps 186 are defined for each time frame in which the slice images are generated. By computing the minimum and maximum volumes, as represented at reference numerals 188 and 190, the ejection fraction may be determined as a ratio of these quantities.

[0102] As noted above, the foregoing technique may be employed for identifying and measuring features of interest in any of a variety of contacts. Moreover, the technique is equally applicable to static and dynamic structures or tissues. FIGS. 11 and 12 represent a single reconstructed image of a brain and a subsequent processed reconstruction of the brain as developed from a series of such views. The image 200 is based upon a magnetic resonance image produced in accordance with conventional techniques. The image includes a variety of structures as defined by varying intensity pixels in an image matrix, including brain tissue 202 surrounded by an image boundary 204 identifiable through the foregoing technique. Where such boundaries are identified for a series of related images, such as a series of slices produced of the entire brain of a patient, reconstructed images 206 may be produced as described above in the case of cardiac images. In particular, the reconstructed images are defined by a series of external boundaries 208 which, in a composite rendering provide surfaces 210 in three dimensions. Such images may be associated with programming code for rotation, enlargement, reproduction, measurement, and so forth.

[0103] While the invention may be susceptible to various modifications and alternative forms, specific embodiments have been shown in the drawings and have been described in detail herein by way of example only. However, it should be understood that the invention is not intended to be limited to the particular forms disclosed. Rather, the invention is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the invention as defined by the following appended claims.