Abstract:
A method of digital imaging includes receiving image data and fitting a curve to boundaries within the image data. The curve is fit to the boundaries within the image data by extracting a region of interest from the image data and computing a signed distance transform in a narrow band within the region of interest. Finite difference equations including various variables are solved to determine a rate at which the distance transform changes. The distance transform is then diffused at that rate. The technique is based on region-based diffusion propagation, pixel classification, and mathematical morphology. The method is implemented to run in the narrow band of the region of interest specified by the user and the computations are implemented using a fast marching method in the narrow band. While idealized for distinguishing segments of white matter, gray matter, and cerebral spinal fluid in the brain, the algorithm can applied to find contours in any digital image.

Description:
BACKGROUND OF THE INVENTION 
     The present application relates to diagnostic medical imaging. The invention finds particular application in segmenting pixel groups within a medical image for display and use in clinical diagnostics, real time image guided surgery, therapy planning, functional MRI, and the like. It finds particular application with computations which are implemented using a “fast marching method” in the “narrow band.” It is to be appreciated however, that the present invention finds further application in segmenting or defining borders in any digitized image. 
     Many different types of medical imaging equipment exist today. The uses of, and the analysis upon many of these images continue to improve. For example, medical sciences are in the process of searching for locations within the human brain for traits like spoken language, reading, and moral reasoning. Currently, imaging techniques are the least intrusive and most favorable techniques available to study different regions within the brain. For example, the recent growth of “fMRI” is revolutionizing the research in the behavior of the brain while engaged in an activity. This branch of MRI is highly dependent upon the classification of different regions in the brain. 
     Another field of brain imaging is magnetoencephalography (MEG) and electroencephalography (EEG). These techniques have enabled researchers to understand brain activity better than ever before. While all of these brain imaging techniques provide a valuable tool for studying the functions of the brain, these techniques typically rely on the inconsistent application of human hands to localize particular regions or areas within the particular medical image. Moreover, frequently there is a large time lag between image acquisition and image segmentation which can delay evaluation, diagnosis, and/or research. 
     The present invention contemplates an improved method and apparatus which overcomes the above-referenced problems and others. 
     SUMMARY OF THE INVENTION 
     In one embodiment of the present invention, a method of digital image presentation includes receiving image data and fitting a curve to boundaries within the image data. At least the curve and the image data are registered and processed for human readable display. 
     In accordance with another aspect of the present invention, a method of segmenting a medical image includes determining a regional interest on the medical image and computing a propagation speed indicative of a rate at which contour changes. The method further includes computing an altered contour within the region of interest based on a previous contour and the propagation speed. A final contour is eventually extracted from the region of interest and displayed to a user. 
     In accordance with another aspect of the present invention, the region of interest includes a first set of pixels distinguishable from other sets of pixels in the medical image. The extracting a final contour step includes repetitively adjusting the altered contour until the altered contour substantially circumscribes the first set of pixels. 
     In accordance with another aspect of the present invention, particularly for cerebral images, the first set of pixels includes one of the set of white matter, gray matter, and cerebral spinal fluid. 
     In accordance with another aspect of the present invention, the first set of pixels includes pixels having a defines similarity to each other. 
     In accordance with another aspect of the present invention, the propagation includes curvature speed relating to curvature of the contour. The speed of the curve or contour propagation is controlled by the regional constant which one can change depending on the size of the medical organ or object to be segmented. If the object is large and if the capture range is large, which implies a large distance to cover, then the waiting factor is automatically adjusted. 
     In accordance with another aspect of the present invention, the propagation speed includes regional speed relating to the determined region of interest. This regional speed is computed using fuzzy characteristics of the regions. These fuzzy characteristics are the membership functions which tell the contribution of each pixel in each of the identified classes. The number of classes are user defined and can thus be changed to improve the accuracy of the segmentation process. 
     In accordance with another aspect of the present invention, the propagation speed includes gradient speed relating to gradient information of the medical image. This information is computed from the pixel classification process. 
     In accordance with another aspect of the present invention, the propagation speed includes fuzzy gradient shape, so called shape-based speed. The shape-based speed is computed using gradient methods from the pixel classified image or membership images. This shape-based speed is a combination of gradient and fuzzy characteristics. 
     In accordance with another aspect of the present invention, the method further includes computing a signed distance transform of the previous contour using a curve layering method in a band surrounding the contour. 
     In accordance with another aspect of the present invention, the computing assigned distance transform step includes determining an accepted set of pixels. A trial set and a far set of pixels is then tagged and distances of the trial set from the accepted set, and of the far set from the accepted set are calculated. The curve layering or fast marching of the pixels is accomplished by testing 32 variant combinations and solving Eikonal equations. 
     In accordance with another aspect of the present invention, the medical image is registered with the final contour and displayed. Additionally, is an ability to register images from multiple sources for the same organ or object of interest. A segmented contour can then be computed for both images and displayed. 
     One advantage of the present invention resides in the increased capture range of a contour within a medical image. 
     Another advantage of the present invention resides in the derivation of the propagation of the curve from the parametric contour deformable model by incorporation of fuzzy characteristics. Of the two classes of deformable models, parametric class and level set class, the present invention derives the level set class from the parametric deformable class which yields all the properties of the parametric class. Implementation using the level set class offers advantages of both the parametric and level set classes. 
     Another advantage of the present invention resides in the ability to handle cavities, concavities, convolutedness, splitting or merging of the contours as they grow. 
     Another advantage of the present invention resides in the ability to prevent shock formations such as first order, second order, third order and fourth order. 
     Another advantage of the present invention lies in the controlled speed with which the curve evolves depending on features within the image itself. 
     Yet another advantage of the present invention resides in ability to duplicate image contours consistently or apply similar processing to variant images. 
     Still further advantages will become apparent to those of ordinary skill in the art upon reading and understanding the following detailed description. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     The invention may take form in various components and arrangements of components and in various steps and arrangements of steps. The figures are only for purposes of illustrating the preferred embodiments and are not to be construed as limiting the invention. 
     FIG. 1 is a diagrammatic illustration of an imaging apparatus employing the present invention; 
     FIG. 2 is a process-object flowchart overview of the segmentation processor as seen in FIG. 1; 
     FIG. 3 is a process-object flowchart illustrating the speed control system of FIG. 2; 
     FIG. 4 illustrates a process-object flowchart illustrating a regional speed determination as seen in FIG. 3; 
     FIG. 5 illustrates a process-object flowchart illustrating a gradient speed determination as seen in FIG. 3; 
     FIG. 6 illustrates a process-object flowchart illustrating a curvature speed determination as seen in FIG. 3; 
     FIG. 7 illustrates a process-object flowchart illustrating a shape-based speed determination as seen in FIG. 3; 
     FIG. 8 illustrates a process-object chart illustrating diffusion propagation process as illustrated in FIG. 2; 
     FIGS. 9A-9F illustrate successive examples of the evolution of a contour of an MRI image of a human brain; 
     FIG. 10 is an object-process flowchart illustrating the fast marching method algorithm; 
     FIG. 11 is a process-object flowchart illustrating the computation of the zero level curve or contour; 
     FIG. 12 is a process-object illustration of a pixel classification system; 
     FIG. 13 illustrates a process-object flowchart illustrating pixel membership for estimation using an iterative least squares method; 
     FIG. 14 is a process-object flowchart illustrating the quantification of various regions within a digitized medical image; 
     FIG. 15 is a process-object flowchart illustrating the computation of the accepted set as seen in FIG. 10; 
     FIG. 16 is a process-object flowchart illustrating diffusion of various images for visualization or display; 
     FIG. 17 is a process-object flowchart illustrating curve interpolation for fusion and overlay image generation; and, 
     FIG. 18 illustrates a process-object flowchart detailing subsystems of the graphical user interface according to the present invention. 
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
     Referring now to FIG. 1, a diagrammatic illustration of an apparatus to segment medical images is seen. An image generating device  12 , is provided such as a CT scanner, nuclear camera, magnetic resonance imaging, ultra sound imager, or the like. The illustrated CT scanner includes a region of interest  14  through which x-rays  16  are passed. The x-rays  16  are generated by a conventional x-ray tube  20  and are collimated into a thin fan beam  22  or a plurality of generally parallel fans which pass through the region of interest  14 . After passing through the region of interest  14 , the attenuated x-rays are received on a plurality of radiation detectors  30  located opposite of the x-ray tube  20 . Images from the detectors  30  are reconstructed by a reconstruction processor  32  into a volumetric image representation, a series of slice image representations, or a single slice image representation and stored in image memory  34 . Of course, the digital image representations can come from other types of imaging devices or even combinations of imaging devices. 
     The digitized image data is provided to a user interface  40  for interaction from a user (not shown) such as a surgeon. The user manually places an initial contour  42  or zero level curve onto a representation of the digitized image data  44 . The digitized image  44  and the zero level curve  42  are forwarded to a segmentation processor  50  which conforms the curve  42 ′ to distinguishable boundaries within the digitized image data  44 . The final curve  42 ′ is forwarded to a video processor  60  for combination or registration with processed digitized image data for output on display  62 . It should be appreciated that the final curve  42  may be continuous or split into two or more curves as the process progresses. 
     Referring now to FIG. 2, the segmentation processor  50  receives a two-dimensional gray scale representation of an organ cross-section  44 , for example a gray matter and white matter image (GM/WM) or a cerebral spinal fluid emphasized image (CSF). A region of interest (ROI) is identified in the image  44  using mathematical morphology tools  70 . This results in an estimate of the ROI  72 . 
     The initial field distribution or signed distance transform is completed as shown in process  74 . The process  74  receives two input tools; fast marching method tool  76  and the narrow band width  78 . The output is shown in  80 . 
     The next stage is the surface evolution or diffusion propagation system shown in process  84 . It employs a speed control system  86  more fully discussed below, the fast marching method  76  and narrow band method  78 . From the field distribution image we extract the contour using contour extraction program  88  which yields the segmented boundary  42 ′. 
     Referring now to FIG. 3, an overview of the speed control system  86  is detailed which drags the initial contour  42  using four exemplary propagation speed components. From the image  44 , four speed functions are determined called regional speed  90 , gradient or edge-based speed  92 , curvature speed  94 , and shape-based speed  96 . More detail on the individual components is provided in the following Figures. A process  100  integrates all the propagation speed control functions  90 ,  92 ,  94 ,  96  and yields the net speed  106 . This is the net force acting on the evolving curve. 
     Specifically, the propagation speeds are derived mathematically from classical parametric contouring modeling theory. The classical parametric contour model based on internal and external energy is given as: 
     
       
         γ∂ X/∂t=∂/∂s (α∂ X/∂s )−∂ 2   X/∂s (β∂ 2   X/∂s   2 )+ F   ext ( X ) 
       
     
     
       
         
               
               
               
               
               
             
               
               
             
           
               
                   
                   
               
             
             
               
                   
                 first term 
                 ± 
                 second term 
                 external-energy 
               
             
          
           
               
                   
                 internal energy terms 
               
               
                   
                   
               
             
          
         
       
     
     where, the first and second terms are the internal energy terms and the last term is the external energy terms. Also note the definition of the first two terms are the internal energy terms while the last term is the external energy term. 
     Since the second term of the internal energy term does not significantly affect the performance of the geometric deformable models, we can neglect this term and replace it with a force term which is a function of the normal to the curve N(X). 
     Thus, F press (X)=w p (X) N(X). 
     Using the relationship ε=α/γ, V P =w p /γ N(X) and V ext =F ext  (X)/γ and changing the transformation in terms of level frame work ∂φ/∂t=V(κ) N, where N=∇φ/¦∇φ¦, and ∂/∂s(α∂X/∂s)·N=ακ, we given the final curve/surface evolution as: 
     
       
         ∂φ/∂ t =(εκ+ V   p )¦∇φ¦− V   ext ·∇φ. 
       
     
     Generalizing this equation to incorporate other speed terms: 
     
       
         ∂φ/∂ t =(εκ+ V   p   +V   s )¦∇φ¦− V   ext ·∇φ, 
       
     
     where V s  is the speed due to shape. 
     So, we see the above equation has 4 speed terms which needs to computed. 
     The numerical implementation of the partial differential equation is given as: 
     
       
         φ n+1 ( x,y )=φ n ( x,y )−Δ t{V   reg ( x,y )+ V   grad ( x,y )+ V   shape ( x,y )− V   cur ( x,y )} 
       
     
     where, φ n+1  (x,y) and φ n  (x,y) are the level set functions at pixel locations (x,y) at times n and n+1. Δt is the time difference. 
     Regional speed is computed using partial differential equations, gradients of level sets and forward and backward differences as under: 
     
       
           V   reg ( x,y )=max{ V   p ( x,y ),0}∇ + +min{ V   p ( x,y ),0}∇ −   
       
     
     
       
           V   p ( x,y )= w   R *{1−2 u ( x,y )} −1   
       
     
     
       
         ∇ + ={∇ x   +   + ∇ y   + } 1/2   
       
     
     
       
         ∇ − ={∇ x   −   + ∇ y   − } 1/2   
       
     
     
       
         ∇ x   + =max{ D   −x ( x,y ),0} 2 +min{ D   +x ( x,y ),0} 2   
       
     
     
       
         ∇ y   + =max{ D   −y ( x,y ),0} 2 +min{ D   +y ( x,y ),0} 2   
       
     
     The difference operators are defined in terms of level set functions as: 
     
       
           D   −x ( x,y )={φ( x,y )−φ( x −1 ,y )}{Δ x } −1   
       
     
     
       
           D   +x ( x,y )={φ( x +1 ,y )−φ( x,y )}{Δ x } −1   
       
     
     
       
           D   −y ( x,y )={φ( x,y )−φ( x,y −1)}{Δ x } −1   
       
     
     
       
           D   +y ( x,y )={φ( x,y +1)−φ( x,y )}{Δ x } −1   
       
     
     Note w R  is the weighting factor which controls the regional speed and convergence speed of the deformable model. u(x,y) is the membership function computed from the fuzzy mean clustering algorithm or pixel classification algorithm, given the number of classes of tissues in the image. For example, for the brain image, the number of classes are 4, WM, GM, CSF, and background. For a CT image, the number of classes could be less while in pathology images the number of classes could be more. 
     Gradient speed is computed using partial differential equations, gradients of level sets and forward and backward differences as under: 
     
       
           V   grad ( x,y )= V   gradx ( x,y )+ V   grady ( x,y ) 
       
     
     
       
           V   gradx ( x,y )=max{ p   n ( x,y ),0 }D   −x ( x,y )}+min{ q   n ( x,y ),0 }D   +x ( x,y )} 
       
     
     
       
           V   grady ( x,y )=max{ q   n ( x,y ),0 }D   −y ( x,y )}+min{ q   n ( x,y ),0 }D   +y ( x,y )} 
       
     
     
       
           p   n ( x,y )=∇ x   {w   e ∇( Gσ*I )} 
       
     
     
       
           q   n ( x,y )=∇ y   {w   e ∇( G   σ   *I )} 
       
     
     where, I is the original image, G σ  is the Guassian operator with known standard deviation σ, and D −x (x,y), D +x (x,y), D −y (x,y), D +y (x,y) are the difference operators given as: 
     
       
           D   −x ( x,y )={(φ( x,y )−φ( x −1 ,y )}{Δ x}   −1   
       
     
     
       
           D   +x ( x,y )={(φ( x 1 ,y )−φ( x,y )}{Δ x}   −1   
       
     
     
       
           D   −y ( x,y )={(φ( x,y )−φ( x,y −1)}{Δ x}   −1   
       
     
     
       
           D   +y ( x,y )={(φ( x,y +1)−φ( x,y )}{Δ x}   −1   
       
     
     Note ∇ is the gradient operator and gradient is computed after smoothing the original image I with the Gaussian operator G σ . The output of the process ∇ (G σ *I) is an edge image controlled by w e , the edge weight factor and brings the robustness to the system. The output of p n (x,y)=∇ x {w e ∇(G σ *I)} and q n (x,y)=∇ y {w e ∇(G σ *I)} are the x and y components of the edge image for each pixel location. This is one method of computing the edge image. We can also incorporate any edge detection scheme such as likelihood method for computing the edges. 
     Shape speed is computed using partial differential equations, gradients of level sets and forward and backward differences as under: 
     
       
           V   shapex ( x,y )= V   shapex ( x,y )+ V   shapey ( x,y ) 
       
     
     The x and y components of the shape speed is computed as: 
     
       
           V   shapex ( x,y )=max{ p   n ( x,y ),0 }D   −x ( x,y )}+min{ q   n ( x,y )0 }D   +x ( x,y )} 
       
     
     
       
           V   shapey ( x,y )=max{ q   n ( x,y ),0 }D   −x ( x,y )}+min{ q   n ( x,y )0 }D   +x ( x,y )} 
       
     
     
       
           p   n ( x,y )=∇ x   {w   s ∇( G   σ   *U )} 
       
     
     
       
           q   n ( x,y )=∇ y   {w   s ∇( G   σ   *U )} 
       
     
     where U is the fuzzy membership image computed using the Fuzzy Mean Clustering algorithm, given the original image I and D −x (x,y), D +x (x,y), D −y (x,y), D +y (x,y) are the difference operators given as: 
     
       
           D   −x ( x,y )={φ( x,y )−φ( x −1 ,y )}{Δ x } −1   
       
     
     
       
           D   +x ( x,y )={φ( x +1 ,y )−φ( x,y )}{Δ x } −1   
       
     
     
       
           D   −y ( x,y )={φ( x,y )−φ( x,y −1)}{Δ x } −1   
       
     
     
       
           D   +y ( x,y )={φ( x,y +1)−φ( x,y )}{Δ x } −1   
       
     
     Note, ∇ (G σ *U) is again the edge detection process over the classified image. This brings the system very robust in clamping the deforming curves to the goal position. w s  ∇ (G σ *U) controls the weight of the edge computed from the membership function of the fuzzy clustering or pixel classification method. ∇ x {w s ∇(G σ *U)} and ∇ y {w s ∇(G σ *U)} are the x and y components of the shape speed terms. Note that ∇ x  and ∇ y  are the x-gradient and y-gradient operators run over the edge image. 
     Curvature speed is computed using partial differential equations, gradients of level sets and forward and backward differences as under: 
     
       
           V   cur ( x,y )=εκ n ( x,y ){( D   0x ( x,y )) 2 +( D   0y ( x,y )) 2 } 1/2   
       
     
     Where, κ n (x,y), D 0x (x,y), and D 0y (x,y) are given as: 
     
       
         κ n ( x,y )={φ 2   xx φ 2   y −φ 2   x φ 2   y φ 2   xy +φ 2   yy φ 2   x }{φ 2   x +φ 2   y } −3/2   
       
     
     where, the terms D −0x (x,y) and D −0x (x,y) are defined as: 
     
       
           D   0x ( x,y ){φ( x +1 ,y )−φ( x −1 ,y )}{2 ∇x}   −1   
       
     
     
       
           D   0y ( x,y ){φ( x,y +1)−φ( x,y −1)}{2 ∇y}   −1   
       
     
     Referring now to FIG. 4, a sub-process is used to compute the regional speed term  90  from a gray scale image  44 . Here the system computes a membership function for each pixel location. This is called partial volume compilation. Since each pixel or voxel is composed of a mixture of several tissue classes, a process  110  computes the percentage contribution of each class (tissue type) at a voxel. Such an algorithm is called pixel classification or voxel classification. The output of the process  110  includes an image where each pixel or voxel has been classified by tissue class  112 . If there are “N” classes in the images  44 , the process generates “N” different images, for example three images, one for each of GM, WM and CSF. These resulting classified images are combined with a region weight constant  114  in a regional evaluation processor  116  to give regional values for every pixel point (x, y). A regional force compilation processor  120  inputs calculations from a finite difference processor  122  and a signed distance transform processor  124 . This outputs the regional force term shown in process  126 . 
     Referring now to FIG. 5, the image/edge speed computation  92  is illustrated. This object-process diagram computes the image force due to edge velocity. The input is again the gray scale image  44 . An image gradient is computed  130  based on an edge weight constant  132 . The image gradient is then scaled between 0 and 1, in a process  136 . From the scaled image, edge strength is computed at every pixel location (x and y) in a process  138  and outputs U, V  140  become components in an image force computation process  142  which also employs the finite difference  122  and signed distance transform  124 . The edge component of speed results  144 . 
     Referring now to FIG. 6, the curvature speed term or curvature force determination  94  (FIG. 3) is illustrated. Again, the gray scale image  44  initiates the process. On the gray scale image  44  a signed distance transform process  160  operates using a curve layering method such as the fast-marching method (FMM)  76  and the narrow band method (NBM)  78  where the narrow band surrounds the contour. The output is a signed distance transform image or field phi  162 . A curvature force process  164  combines the curvature contour  166  and the signed distance transform image  162 . Curvature force or curvature velocity is output  168 . 
     Referring now to FIG. 7, the shape-based speed component is also computed from the gray scale image  44 . The shape-based speed is determined by first computing a pixel membership from a specified numbers of classes  170 . The membership values are between the range of 0 and 1. 
     Next, a gradient of the membership image is computed, process  172 , resulting in the gradient image. Components x and y of the shape-based gradient image are calculated, process  174 . These are called the U-V components at each pixel component. The force is then computed in a process  176  from the signed distance transform  124  and the finite difference tool  122 . The output of this sub-system is the shape velocity component  178 . 
     Referring now to FIG. 8, diffusion propagation in the initial field using adaptive narrow band and fast marching method is illustrated. This is the algorithm used for computing the zero-level-curve (ZLC), the final estimated boundary. A user initiates the process by entering the raw contour  42 . The initial field or initial signed distance function is computed in a process  182  resulting in the initial field in the narrow band. Now a new field is computed  184  based on the speed control functions  86  (FIG.  3 ). The output contains the new signed distance transform  186  in the narrow band. 
     Next, the new field is checked to determine whether “land mines” have been hit, decision block  188 . The “land mines” are the border points. If a land mine is not hit, then the loop repeats. If a land mine is hit, the loop is exited indicating that an output contour has been reached, i.e. a point on the interface between two tissue types has been identified. 
     The process of tube reconstruction is repeated, decision block  190 , until all the tubes have been processed. When this occurs, the system exits with a zero level curve (ZLC) estimate  192 . 
     Referring now to FIGS. 9A-9F, a progression illustrates the growth or evolution of the zero level curve from the raw contour entered by the user  42 FIG. 9A, to the final contour  42  illustrated in FIG.  9 F. The narrow band width in the illustrated example was twenty-five pixels on either side of the ZLC, with land mines being five pixels wide. Referring back generally to FIG. 8, the recursive nature of the algorithm disclosed results after a first pass in FIG. 9B evolving from FIG. 9A against a medical image (not illustrated) having defined pixel classifications. It is now apparent that the third iteration results in the illustration of FIG. 9C, and so on until the final contour  42  is reached as illustrated in FIG.  9 F. 
     Referring now to FIG. 10, the fast-marching algorithm ( 76 , FIG. 8) for signed distance transform computation in the narrow band using a neighborhood search method is illustrated. The input to this algorithm is a raw contour  42  (FIG. 1) as specified by the surgeon or user during image guided procedures. A flood-fill algorithm fills the region, process  204 . Now with the narrow band process  78 , the points or pixels which belong to an Accepted Set (AS) are computed, step  206 . Further detail on this process will be provided below in connection with FIG.  14 . 
     Selected pixel points in A Trial Set and a Far Set are tagged, step  208 , and output. The distances in Trial Set and Far Set points are computed, step  210 . In this way, the signed distance transform (SDT)  212  can be calculated from a raw contour. 
     Referring now to FIG. 11, a method of zero-level-curve computation or contour extraction from the field image is illustrated. This sub-system computes a zero-level-curve (ZLC) from a field image in a process called iso-contour extraction. For every pixel-point, the field signs of 4 neighbors (E, W, N, &amp; S) are checked, process  220 . Now the algorithm checks if there is a change in sign when we go from a central pixel to its neighboring pixel, decision block  222 . 
     If the product is −1, a sign has changed and the algorithm proceeds. If there is no change in sign, no changes are implemented, object  224 . If changes have occurred the intersection of the curve with the background grid is determined, process  226 . The output is the x,y location of the curve-grid-intersection  228 . After all points in the field have been checked, the ZLC or final contour  230  results. 
     Referring now to FIG. 12, an exemplary pixel classification system is illustrated. Previous figures have revealed the desirability of a pixel classification methodology to compute membership functions for each pixel location (see e.g. step  110  (FIG.  4 ), (FIG.  7 )). A vector is framed from a gray scale image process  236 . Now from an initial number of classes  238 , the initial centroid for each cluster is computed, process  240 . The output is the initial estimate of the centroids. A least squares algorithm is applied in process  242  to compute the labeled image and membership functions  244 . 
     Referring now to FIG. 13, an exemplary membership computation algorithm for regional force is illustrated. A new membership function is computed from the initial centroid, process  250 . A new centroid is computed and normalized, process  252 . If an error threshold  256  has not been reached, the process is repeated. If the centroid error is less than the threshold  256 , then the membership computation function exits. Upon exit, the final membership function for each pixel location are determined, process  258 . 
     Referring now to FIG. 14, an exemplary quantification of white matter, gray matter, and CSF is illustrated. The number of regions in a segmented boundary are counted, process  266 . The area of each of the regions (R 1  . . . R n ) is computed in process  268 , resulting in the output (quantified) regions having areas A 1, A   2  . . . A n  corresponding to “N” regions. 
     Referring now to FIG. 15, an exemplary method to compile the Accepted Set using nearest neighbor search is illustrated. This is a sub-system to compute the distances and tags for the Accepted Set (FIG. 9,  206 ). Process  270  checks the sign of the field image (Φ), given the initial field flow (Φ). Next all signs are checked  272  for positive and negative signs (32 combinations) relative to a central location. Next, fractional distances and tags are computed, process  276 . If all the points in the narrow band are finished, the Accepted Set is complete, if not, the method proceeds to the next point and cycles back to step  270  to check the next set of signs. 
     Referring now to Table 1, an exemplary set of results are illustrated for the case where the central pixel has a positive sign. This chart shows sixteen cases when the central pixel is positive and the neighboring pixels are negative The total number of neighboring pixels which are negative can be 1, 2, 3, or 4. 
     Case 1 to case 4 are when the neighboring pixel is negative. 
     Case 5 to case 10 are when 2 of the neighboring pixels are negative. 
     Case 11 to case 14 are when 3 of the pixels are negative. 
     Case 15 is when all the 4 neighboring pixels are negative. 
     Case 16, when none are negative. 
     
       
         
               
             
               
               
               
               
               
               
             
               
               
               
               
               
               
             
           
               
                 TABLE 1 
               
             
             
               
                   
               
               
                 32 Combination Cases 
               
               
                 If Central Pixel (C) is Positive Sign 
               
             
          
           
               
                   
                 E 
                 W 
                 N 
                 S 
                   
               
               
                   
                   
               
             
          
           
               
                  1. 
                 − 
                 + 
                 + 
                 + 
                 When 
               
               
                  2. 
                 + 
                 − 
                 + 
                 + 
                 1 
               
               
                  3. 
                 + 
                 + 
                 − 
                 + 
                 is 
               
               
                  4. 
                 + 
                 + 
                 + 
                 − 
                 Negative 
               
               
                  5. 
                 − 
                 + 
                 − 
                 + 
                 When 
               
               
                  6. 
                 − 
                 + 
                 + 
                 − 
                 2 
               
               
                  7. 
                 + 
                 − 
                 − 
                 + 
                 are 
               
               
                  8. 
                 + 
                 − 
                 − 
                 + 
                 negative 
               
               
                  9. 
                 − 
                 − 
                 + 
                 + 
               
               
                 10. 
                 + 
                 + 
                 − 
                 − 
                 When 
               
               
                 11. 
                 + 
                 − 
                 − 
                 − 
                 3 
               
               
                 12. 
                 − 
                 + 
                 − 
                 − 
                 are 
               
               
                 13. 
                 − 
                 − 
                 + 
                 − 
                 Negative 
               
               
                 14. 
                 − 
                 − 
                 − 
                 + 
                   
               
               
                 15. 
                 − 
                 − 
                 − 
                 − 
                 When 4 
               
               
                   
                   
                   
                   
                   
                 are negative 
               
               
                 16 
                 0 
                 0 
                 0 
                 0 
                 When none 
               
               
                   
                   
                   
                   
                   
                 are negative 
               
               
                   
               
             
          
         
       
     
     Referring now to Table 2, an exemplary set of results are illustrated for the case where the central pixel has a negative sign. 
     Case 17 to case 20 shows when 1 neighboring pixel is positive. 
     Case 21 to case 26 shows when 2 neighboring pixels are positive. 
     Case 27 to case 30 shows when 3 neighboring pixels are positive. 
     Case 31 shows when all 4 neighboring pixels are positive. 
     Case 32 shows when none of the 4 neighboring pixels are positive. 
     
       
         
               
             
               
               
               
               
               
               
             
               
               
               
               
               
               
             
           
               
                 TABLE 2 
               
             
             
               
                   
               
               
                 If Central Pixel (C) is Negative Sign 
               
             
          
           
               
                   
                 E 
                 W 
                 N 
                 S 
                   
               
               
                   
                   
               
             
          
           
               
                 17. 
                 + 
                 − 
                 − 
                 − 
                 When 
               
               
                 18. 
                 − 
                 + 
                 − 
                 − 
                 1 
               
               
                 19. 
                 − 
                 − 
                 + 
                 − 
                 is 
               
               
                 20. 
                 − 
                 − 
                 − 
                 + 
                 Positive 
               
               
                 21. 
                 + 
                 − 
                 + 
                 − 
               
               
                 22. 
                 + 
                 − 
                 − 
                 + 
                 When 
               
               
                 23. 
                 − 
                 + 
                 − 
                 + 
                 2 
               
               
                 24. 
                 − 
                 + 
                 + 
                 − 
                 are 
               
               
                 25. 
                 + 
                 + 
                 − 
                 − 
                 Positive 
               
               
                 26. 
                 − 
                 − 
                 + 
                 + 
                 When 
               
               
                 27. 
                 − 
                 + 
                 + 
                 + 
                 3 
               
               
                 28. 
                 + 
                 − 
                 + 
                 + 
                 are 
               
               
                 29. 
                 + 
                 + 
                 − 
                 + 
                 Positive 
               
               
                 30. 
                 + 
                 + 
                 + 
                 − 
               
               
                 31. 
                 + 
                 + 
                 + 
                 + 
                 When 4 
               
               
                   
                   
                   
                   
                   
                 are Positive 
               
               
                 32. 
                 0 
                 0 
                 0 
                 0 
                 When none 
               
               
                   
                   
                   
                   
                   
                 are Positive 
               
               
                   
               
             
          
         
       
     
     Referring now to FIG. 16, an exemplary method is illustrated to show a fusion of 2 curves over the gray scale image. The two curves could be the raw initial curve (as drawn by the surgeon or user) and the second curve could be the estimated boundary of either WM, GM, or CSF for example. 
     A number of points P 1  which constitute the raw curve are defined, step  300 . The points P 1  are interpolated into a modified curve P 2 , step  302 . Curve P 2  is converted into an image, in process  304 , such as a raw contour image. The gray scale image is also provided to the segmentation system  50  (FIG. 1) which yields the estimated boundary image. The boundary image is fused and/or registered with the raw image, process  308 . 
     This fused output is fused again with the raw contour from step  304  to yield 2 contours (raw and estimated) fused with the background gray scale image in process  310 . 
     Referring now to FIG. 17, an exemplary method of curve interpolation for fusion and over-image generation based on sampling is illustrated. This is a sub-system which generates an interpolated curve image given the discrete contour (say P 1  number of points). An arc length or partial perimeter is computed, process  312  from which an associated arc interval is computed, process  314 . The x-coordinate is interpolated (with P 2  number of points) in process  316 . In parallel, the y-coordinate is interpolated (with P 2  number of points) in process  318 . Each x and y are joined as one curve, process  320 , and converted to an image. This is called interpolated curve image. 
     Referring now to FIG. 18, a graphical user-interface of the segmentation engine is illustrated. This is an exemplary sub-system for creating the graphical user-interface (GUI), preferably using tcl/tk. 
     The main script is invoked, step  350 , having access to the gray scale image  44  and imaging package  354 . In the illustrated embodiment, the graphical interface has 3 buttons, namely an exit button  360 , raw contour button  362  and segmentation button  364 . On invoking the script  350  the image  44  is displayed, step  370 . The surgeon or user  372  manipulates a mouse  374  to draw the initial contour or points  42  on the image, step  378 . The points are registered and plotted over the image, step  382 . 
     Upon selection of the segmentation button  364 , the segmented boundary is created as seen in step  384 . Upon selection of the exit button  360 , the system exits as seen in step  390 . 
     The invention has been described with reference to the preferred embodiment. Obviously, modifications and alterations will occur to others upon reading and understanding the preceding detailed description. It is intended that the invention be construed as including all such modifications and alterations insofar as they come within the scope of the appended claims or the equivalents thereof.