Segmentation method and apparatus for medical images using diffusion propagation, pixel classification, and mathematical morphology6785409Abstract 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. Claims Having thus described the preferred embodiment, the invention is now claimed to be: Description BACKGROUND OF THE INVENTION
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.sub.press (X)=w.sub.p (X) N(X). Using the relationship .epsilon.=.alpha./.gamma., V.sub.P =w.sub.p /.gamma. N(X) and V.sub.ext =F.sub.ext (X)/.gamma. and changing the transformation in terms of level frame work .differential..phi./.differential.t=V(.kappa.) N, where N=.gradient..phi./.linevert split..gradient..phi..linevert split., and .differential./.differential.s(.alpha..differential.X/ .differential.s).multidot.N=.alpha..kappa., we given the final curve/surface evolution as: .differential..phi./.differential.t=(.epsilon..kappa.+V.sub.p).linevert split..gradient..phi..linevert split.-V.sub.ext.multidot..gradient..phi.. Generalizing this equation to incorporate other speed terms: .differential..phi./.differential.t=(.epsilon..kappa.+V.sub.p +V.sub.s).linevert split..gradient..phi..linevert split.-V.sub.ext.multidot..gradient..phi., where V.sub.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: .phi..sup.n+1 (x,y)=.phi..sup.n (x,y)-.DELTA.t{V.sub.reg (x,y)+V.sub.grad (x,y)+V.sub.shape (x,y)-V.sub.cur (x,y)} where, .phi..sup.n+1 (x,y) and .phi..sup.n (x,y) are the level set functions at pixel locations (x,y) at times n and n+1. .DELTA.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.sub.reg (x,y)=max{V.sub.p (x,y),0}.gradient..sup.+ +min{V.sub.p (x,y),0}.gradient..sup.- V.sub.p (x,y)=w.sub.R *{1-2u(x,y)}.sup.-1 .gradient..sup.+ ={.gradient..sub.x.sup.+.sub.+.gradient..sub.y.sup.+ }.sup.1/2 .gradient..sup.- ={.gradient..sub.x.sup.-.sub.+.gradient..sub.y.sup.- }.sup.1/2 .gradient..sub.x.sup.+ =max{D.sup.-x (x,y),0}.sup.2 +min{D.sup.+x (x,y),0}.sup.2 .gradient..sub.y.sup.+ =max{D.sup.-y (x,y),0}.sup.2 +min{D.sup.+y (x,y),0}.sup.2 The difference operators are defined in terms of level set functions as: D.sup.-x (x,y)={.phi.(x,y)-.phi.(x-1,y)}{.DELTA.x}.sup.-1 D.sup.+x (x,y)={.phi.(x+1,y)-.phi.(x,y)}{.DELTA.x}.sup.-1 D.sup.-y (x,y)={.phi.(x,y)-.phi.(x,y-1)}{.DELTA.x}.sup.-1 D.sup.+y (x,y)={.phi.(x,y+1)-.phi.(x,y)}{.DELTA.x}.sup.-1 Note w.sub.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.sub.grad (x,y)=V.sub.gradx (x,y)+V.sub.grady (x,y) V.sub.gradx (x,y)=max{p.sup.n (x,y),0}D.sup.-x (x,y)}+min{q.sup.n (x,y),0}D.sup.+x (x,y)} V.sub.grady (x,y)=max{q.sup.n (x,y),0}D.sup.-y (x,y)}+min{q.sup.n (x,y),0}D.sup.+y (x,y)} p.sup.n (x,y)=.gradient..sub.x {w.sub.e.gradient.(G.sigma.*I)} q.sup.n (x,y)=.gradient..sub.y {w.sub.e.gradient.(G.sub..sigma. *I)} where, I is the original image, G.sub..sigma. is the Guassian operator with known standard deviation .sigma., and D.sup.-x (x,y), D.sup.+x (x,y), D.sup.-y (x,y), D.sup.+y (x,y) are the difference operators given as: D.sup.-x (x,y)={(.phi.(x,y)-.phi.(x-1,y)}{.DELTA.x}.sup.-1 D.sup.+x (x,y)={(.phi.(x1,y)-.phi.(x,y)}{.DELTA.x}.sup.-1 D.sup.-y (x,y)={(.phi.(x,y)-.phi.(x,y-1)}{.DELTA.x}.sup.-1 D.sup.+y (x,y)={(.phi.(x,y+1)-.phi.(x,y)}{.DELTA.x}.sup.-1 Note .gradient. is the gradient operator and gradient is computed after smoothing the original image I with the Gaussian operator G.sub..sigma.. The output of the process .gradient. (G.sub..sigma. *I) is an edge image controlled by w.sub.e, the edge weight factor and brings the robustness to the system. The output of p.sup.n (x,y)=.gradient..sub.x {w.sub.e.gradient.(G.sub..sigma. *I)} and q.sup.n (x,y)=.gradient..sub.y {w.sub.e.gradient.(G.sub..sigma. *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.sub.shapex (x,y)=V.sub.shapex (x,y)+V.sub.shapey (x,y) The x and y components of the shape speed is computed as: V.sub.shapex (x,y)=max{p.sup.n (x,y),0}D.sup.-x (x,y)}+min{q.sup.n (x,y)0}D.sup.+x (x,y)} V.sub.shapey (x,y)=max{q.sup.n (x,y),0}D.sup.-x (x,y)}+min{q.sup.n (x,y)0}D.sup.+x (x,y)} p.sup.n (x,y)=.gradient..sub.x {w.sub.s.gradient.(G.sub..sigma. *U)} q.sup.n (x,y)=.gradient..sub.y {w.sub.s.gradient.(G.sub..sigma. *U)} where U is the fuzzy membership image computed using the Fuzzy Mean Clustering algorithm, given the original image I and D.sup.-x (x,y), D.sup.+x (x,y), D.sup.-y (x,y), D.sup.+y (x,y) are the difference operators given as: D.sup.-x (x,y)={.phi.(x,y)-.phi.(x-1,y)}{.DELTA.x}.sup.-1 D.sup.+x (x,y)={.phi.(x+1,y)-.phi.(x,y)}{.DELTA.x}.sup.-1 D.sup.-y (x,y)={.phi.(x,y)-.phi.(x,y-1)}{.DELTA.x}.sup.-1 D.sup.+y (x,y)={.phi.(x,y+1)-.phi.(x,y)}{.DELTA.x}.sup.-1 Note, .gradient. (G.sub..sigma. *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.sub.s .gradient. (G.sub..sigma. *U) controls the weight of the edge computed from the membership function of the fuzzy clustering or pixel classification method. .gradient..sub.x {w.sub.s.gradient.(G.sub..sigma. *U)} and .gradient..sub.y {w.sub.s.gradient.(G.sub..sigma. *U)} are the x and y components of the shape speed terms. Note that .gradient..sub.x and .gradient..sub.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.sub.cur (x,y)=.epsilon..kappa..sup.n (x,y){(D.sup.0x (x,y)).sup.2 +(D.sup.0y (x,y)).sup.2 }.sup.1/2 Where, .kappa..sup.n (x,y), D.sup.0x (x,y), and D.sup.0y (x,y) are given as: .kappa..sup.n (x,y)={.phi..sup.2.sub.xx.phi..sup.2.sub.y -.phi..sup.2.sub.x.phi..sup.2.sub.y.phi..sup.2.sub.xy +.phi..sup.2.sub.yy.phi..sup.2.sub.x }{.phi..sup.2.sub.x +.phi..sup.2.sub.y }.sup.-3/2 where, the terms D.sup.-0x (x,y) and D.sup.-0x (x,y) are defined as: D.sup.0x (x,y){.phi.(x+1,y)-.phi.(x-1,y)}{2.gradient.x}.sup.-1 D.sup.0y (x,y){.phi.(x,y+1)-.phi.(x,y-1)}{2.gradient.y}.sup.-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 42FIG. 9A, to the final contour 42 illustrated in FIG. 9F. 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. 9F. 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, & 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.sub.1 . . . R.sub.n) is computed in process 268, resulting in the output (quantified) regions having areas A.sub.1, A.sub.2 . . . A.sub.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 (.PHI.), given the initial field flow (.PHI.). 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.sub.1 which constitute the raw curve are defined, step 300. The points P.sub.1 are interpolated into a modified curve P.sub.2, step 302. Curve P.sub.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.sub.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.sub.2 number of points) in process 316. In parallel, the y-coordinate is interpolated (with P.sub.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.
|
Same subclass Same class Consider this |
||||||||||