The present invention relates generally to segmentation of an image, and more particularly to segmentation of a target structure of an image based on the solution of distinct partial differential equations (PDEs).
Content extraction from images typically relies on segmentation, i.e., extraction of the borders of target structures. Automated segmentation by computer algorithms has been a focus of research. In practice, the accuracy of a segmentation algorithm can be hampered by noise in the image acquisition and the complexity of the arrangement of target objects with respect to their surroundings within the image. In order to increase the accuracy of a segmentation algorithm, the algorithm typically has to become more complex and this often leads to an increase in computational cost.
Therefore, there remains a need to more efficiently and accurately segment a target structure of an image.
The inventors have invented a method and system for segmenting a target structure of an image more efficiently and accurately. In one embodiment, one or more inputs to designate (e.g., mark) a region of interest (ROI) of the image is received. The region of interest is the initial input from the user that points to the image domain that is going to be used by the segmentation algorithm. The segmentation aims to segment a target structure (e.g., a lymph node).
A first seed is then set outside of the target structure and a second seed is set inside the target structure. A first partial differential equation (PDE) associated with the first seed is solved (e.g., by a fast marching algorithm) and a second PDE associated with the second seed is also solved (e.g., by a fast marching algorithm) to segment the image. The first PDE and/or the second PDE can be either an Eikonal differential equation or a diffusion equation.
The solving of the PDEs enables the determination of an exterior distance function associated with a plurality of pixels outside the target structure and an interior distance function associated with a plurality of pixels inside the target structure. A competition criterion is then used between the interior distance function and the exterior distance function to determine and segment the target structure. The determining of the target structure can additionally include determining whether the interior distance function or the exterior distance function has a minimum distance from the first seed (or the second seed). The solution of the PDEs to produce the distance functions enables the increase in accuracy and efficiency in the segmentation of an image.
These and other advantages of the invention will be apparent to those of ordinary skill in the art by reference to the following detailed description and the accompanying drawings.
In accordance with an embodiment of the present invention, distance functions distinguishing between image locations on opposite sides of an edge or, alternatively, identifying locations of higher gradient magnitude (i.e., a greater change in light intensity between two adjacent pixels of an image) are computed. To obtain these distance functions, partial differential equations (PDEs), such as the Eikonal equation or diffusion equations, are solved. In one embodiment, the PDEs can be computed according to a function of the form O(N log N), where N is the number of image pixels. In the case where diffusion PDEs are used, edge information is propagated from the interior of the desired anatomical structure or from the boundary of the region of interest.
Once the distance functions are computed, a distance function competition occurs in order to determine the target structure of the ROI and to subsequently segment the target structure from the ROI.
Segmentation by Distance Function Competition
The following description describes the present invention in terms of the processing steps required to implement an embodiment of the invention. These steps may be performed by an appropriately programmed computer, the configuration of which is well known in the art. An appropriate computer may be implemented, for example, using well known computer processors, memory units, storage devices, computer software, and other components. A high level block diagram of such a computer is shown in
One skilled in the art will recognize that an implementation of an actual computer will contain other components as well, and that
In step 316, a distance function for each known image region is determined. In particular, a first distance function for the PDE associated with the first seed is determined and another second distance function for the PDE associated with the second seed is determined In one embodiment, the first distance function represents the distance to the nearest of a set of prespecified points (i.e., seed or seeds) interior to the desired structure The second distance function represents the distance to a set of prespecified points (i.e., seed or seeds) exterior to the structure.
A competition criterion between the two distance functions is then used to determine which image pixels belong to the interior region of the image (i.e., pixels within the target structure) and which pixels belong to the exterior region of the image (i.e., pixels outside the target structure) in step 320. This determination is then used to segment the target structure from the ROI of the image in step 324.
In one embodiment, the local travel cost for each distance function depends on the local image intensity variation. A travel cost is the cost associated by the step taken by the fictitious front that is propagating. The local travel cost at each image coordinate x is determined locally by the image characteristics in a neighborhood of x. Further, regions that are more likely to be edges are interpreted as regions that have higher local distance. Here, local distance is used interchangeably with the local travel cost. For example, the computer can first generate a binary map (i.e., a 2D array of ones and zeros, with ones representing the edges) after detecting edges of the image (e.g., using a Canny edge detector). Edges may be considered impassable obstacles and the distance function can be computed accordingly. Alternatively, the local distance can be defined as a function of the gradient magnitude of the image. In another embodiment, a set of local distances can be combined with different weights (i.e., coefficients).
In more detail, the following equation, which is referred to as the Eikonal equation, is a partial differential equation (PDE):
∥∇D∥=F,D=0 on G (1)
whose solution D:Ω⊂”→, where n is the dimensionality of the image, represents the arrival time of a moving front with spatially varying speed, 1/F, that starts at a given set of points, G⊂Ω, at time 0, and Ω is the domain on which the distance function D is defined. When the speed of the front is uniform within the domain Ω, the arrival time is proportional to the minimum distance to G, the set of starting points. As used herein, D denotes a distance function (i.e., the distance to the set G with locally varying travel cost, F).
Local distance or travel cost varies accordingly with the presence of an image edge or with local intensity variation. A “fast marching algorithm” can yield an efficient solution to the Eikonal equation on a uniform discrete grid. While the proposed steps apply for continuous image domains with differentiable images, discrete grid locations are referred to below and, as a result, finite difference approximations to the derivatives are applied.
The two distance functions Di and De are then computed in steps 512-532 by solving the two Eikonal PDE's (e.g., through fast marching). In particular, the pixels that are neighbors of the already Known points are labeled as Trial pixels in step 512. Other image pixels are labeled as Far points. It is then determined whether there is a Trial pixel present in step 516. If so, the Trial pixel with the lowest distance value, q, is selected, the pixel is labeled as a Known pixel, and it is verified that each neighbor pixel to q that is not Known is labeled as a Trial pixel while its value is updated according to the chosen distance function in step 520. The process then returns to step 516 and it is again determined whether there is a Trial pixel present.
If not, then the distance function to the interior set is computed with local travel cost, F, in step 524. The value of each pixel then corresponds to the weighted distance to the interior set (Di). In one embodiment, this step is initialized with interior points as Known set. The distance function to the exterior set (De) is then computed in step 528 by repeating steps 512-524 with respect to the exterior set (De). In one embodiment, this step is initialized with exterior points as Known set.
The interior region of the image is then determined in step 532. The interior region is the set of points whose interior distance value is less than the exterior distance value, i.e., the interior set is {(x,y):Di(x,y)<De(x,y)} in the case of a two-dimensional image. The local travel cost, F, of the distance functions are explained in more detail below.
In one embodiment, to solve an Eikonal PDE, a fast marching algorithm is used.
Fast Marching with Edge Map
In one embodiment, the distance function is computed in a way such that edge pixels represent points where the moving front cannot propagate. The Eikonal equation is then,
where E is the edge map that assumes the value of 1−ε where there are edges and the value of 0 at all other pixels. The computer performs this computation to determine the solution where ε>0 approaches 0 in order to represent locally infinite travel cost. The edge map can be derived from any edge detection algorithm that has binary output. In one embodiment, a Canny edge detector is used.
As one skilled in the art will recognize, Canny edge detection algorithms report an edge map that has no a priori known topology, i.e., it does not necessarily partition the region of interest into clear interior and exterior regions. The problem of obtaining a labeling for each pixel as an interior or exterior region is therefore not typically solved by edge detection alone. In accordance with the present invention, the competition algorithm is used in addition to the Canny edge detection algorithm.
The first column 612 depicts the two distance functions computed by starting from both exterior and interior seed points. The distance shown in
Fast Marching with Gradient
In one embodiment, regions of interest with high gradient magnitude are treated as having high local travel cost, and regions with low gradient magnitude are treated as having low local distance. The Eikonal equation then takes the form:
∥∇D∥=∥∇I∥, D=0 on G, (3)
In one embodiment, this method is more tolerant and resistant to errors in the edge map (e.g., from noise in the image) since the technique allows moderate levels of intensity variation to affect the local travel cost by a moderate amount instead of being classified as either edge or non-edge.
The above method can be further generalized by considering the Eikonal equation, ∥∇D∥=ƒ(∥∇I∥). In some embodiments of the function f, particularly those that resemble thresholding functions such as the sigmoid, show the relationship between the current method and the method described above since edge maps typically resemble such functions of the image gradient magnitude, ∥∇I∥. In one embodiment, the second column 616 of
Combined Method
The edge map and gradient methods described above can also be combined. In one embodiment, the computation of the edge map, E, results in a binary image assuming the value E=1 on edge pixels and E=0 elsewhere. This binary image is then directly added to the gradient image by a factor α. The Eikonal equation then takes the form:
∥∇D∥=(∥∇I∥+α*E). (4)
A specific choice of the parameter α depends on the level of trust that can be placed in the edge map, with higher trust corresponding to higher values of α. This will result in increased gradient effects where there are edges.
Note that this technique, which assigns additional travel cost to areas where there are definite image edges, considers the Eikonal equation ∥∇D∥=ƒ(∥∇I∥), i.e., with non-linear functions of the gradient image. For instance, choosing the function,
where t represents an image gradient magnitude threshold and can be chosen automatically based on the range values that the magnitude of the image gradient assumes. Many choices can be made for the function, f.
In accordance with an embodiment of the present invention, different PDEs can be used for each region of an image (i.e., the interior distance function can be obtained by solving a different PDE than the PDE solved to produce the exterior distance function). This flexibility may, for example, assist in the segmentation of interior regions that are textured.
Further, in one embodiment, an interior intensity based term can be added to the interior distance function. This added term can smooth the local gradient and decrease some texture or noise influence. Due to the nature of the main application (e.g., lymph node segmentation), the exterior distance function may not be smoothed in a similar fashion because exterior regions may include other structures that may interfere with the segmentation.
In one embodiment, in order to smooth the exterior distance function, the processor computes the mean intensity of a set of points adjacent to the foreground seed points as Î. The image at each pixel, p, then has a local weight of (I(p)−Î)2, which is added to the local travel cost in the Eikonal equation for the interior region with a weighting parameter, β, as follows:
∥∇Di∥=(∥∇I∥+αE+β(I−Î)2), (6)
where E is the binary edge map and Î is the mean intensity of points adjacent to the interior region seed points, Thus, the final combined method assumes the use of Eq. 4 for the computation of the exterior distance function and Eq. 6 for the computation of the interior distance function.
Diffusion Equation
The above embodiments have illustrated how Eikonal PDEs can be used. In another embodiment, a diffusion equation is used to propagate image information rather than the Eikonal equation. Although the nature of information propagation in diffusion equations is typically different than the information propagation in the Eikonal equation (e.g., diffusion equations often propagate information with infinite speed), the diffusion equation nonetheless propagates information in a similar manner to the previously mentioned distance function based techniques.
The linear heat equation on a function D(
where ΔD is the Laplacian operator. In one embodiment, the initial conditions are D(
hence diffusing edge information from the boundaries towards the non-boundary regions.
Similar to the Eikonal equation and fast marching techniques, where the information is propagated from the boundaries or the seeds of the image domain towards unlabeled points, diffusion equations can also be used for segmentation by creating two smooth distance functions, one for the interior seed points and one for the exterior seed points. For the interior distance function, Di, the boundary conditions are illustratively set to 1 at the interior seed points and the function is initialized illustratively to a value of 0 at all other points in the image domain. For the exterior distance function, De, the boundary conditions are illustratively set to 1 at the exterior seed points and the function is initialized illustratively to a value of 0 at all other points in the image domain.
To introduce image dependent terms to the diffusion equation, an anisotropic diffusion that depends on the local image variation is used, thereby allowing less diffusion in directions where the image derivative is higher and more diffusion where the image derivative is low, The definition of the four one-sided image derivatives around a pixel can be given by
Ix−(x,y)=I(x,y)−I(x−1,y),Ix+(x,y)=I(x+1,y)−I(x,y)
Iy−(x,y)=I(x,y)−I(x,y−1),Iy+(x,y)=I(x,y+1)−I(x,y)
An image-based discrete diffusion equation can be generated by introducing the image-driven weights to the discrete Laplacian equation as follows:
Note that γ represents a damping coefficient that affects the level of anisotropy inherent in this method. Higher values of γ allow for greater anisotropy, i.e., allow for the information diffusion to be more sensitive to differences in image intensity across the image. A reasonable range for this parameter may be, for example, 0.001<γ<0.01.
Using the set of seeds for the exterior region and the interior region as two distinct set of boundary conditions, the two distance functions, De and Di, corresponding to the exterior and interior are estimated after a set amount of diffusion time. Similar to the approach described above with respect to the Eikonal equation, the segmentation map is formed by considering the interior region (i.e., interior of the target structure) to be the set of points where the interior distance function is higher than the exterior distance function. The third column 620 in
In one embodiment, the diffusion is executed for a sufficiently long time for all pixels to be affected by the diffusion but for a short enough time for the diffusion to be practically useful, i.e., to avoid the constant solution D=1 along the entire domain. In practice, k=1000 forward Euler iterations often produces suitable results and the results are not particularly sensitive to moderate variations in k.
Two Region Segmentations
In one embodiment, the algorithm described above is not sensitive to the placement of the interior and exterior seed points. It may be possible to use any set of exterior or interior seed points as long as they are clearly outside or inside of the target structure, respectively. In one embodiment, a mouse drag operation is used on the image to set exterior seeds in the form of a 2D rectangular border, as shown in
Segmentation Results
The Canny edge detector typically propagates strong edges and discards the weak edges, and this can lead to errors, e.g., either “holes” in the edge map as in row one 728, or edge noise as in rows three 732, four 736, and five 740. This can directly influence the distance functions and, in turn, may influence the final segmentation.
These errors may be reduced by the second approach that uses image gradient in the Eikonal PDE. The distances found are more tolerant and resistant to errors in the edge map functions, and the segmentation matches the desired structure more closely,
Referring once again to the example of a lymph node, discussed herein above, in cases where a strong edge is situated near the edges of the lymph node, but yet external to the lymph node itself, the gradient method may be used, as shown in the bottom right of the image in row 1 728 and column 716 and row 6 744 and column 716. Note that such errors may be fixed by assuming more than two regions are present in the image. These errors can also be reduced via the diffusion method described herein above.
The diffusion method typically performs well when image edges are strong and robust to higher levels of image noise, especially to noise that occurs at a single pixel (i.e., “salt and pepper” noise). The reason for this is the ability for the diffusion equation to propagate information around a single pixel that, due to noise, has an abnormally high intensity value, This property does not typically hold for the image gradient based techniques because such noise affects the image gradient in a neighborhood around the pixel, thereby creating a larger region of high local travel cost. This method may be prone to error when the target objects are merged with other external structures containing high edge content (e.g., as shown in the images of the fifth row 740).
The results of the combined method discussed herein above are shown in column 724. This method specifically prohibits the propagation of information where an edge is detected (whereas such information often propagates in the image gradient method). However, when an edge is not detected, the combined method allows for varying degrees of information propagation depending on the magnitude of the image gradient.
Segmentation in 3D
Segmentation in 3D can also be achieved using the Eikonal PDEs and diffusion PDEs by extending the fast marching (in the case of the Eikonal PDE) and the gradient computations to three dimensions
Multi-Region Segmentation
In accordance with an embodiment of the present invention, the techniques described above can also be used for multi-region segmentation. Instead of choosing only interior and exterior seed regions as in the two region segmentation described above, the multiple region segmentation technique allows for any number of seed regions, where each seed region can correspond to a target structure of interest within the image.
For each seed region, the processor can compute a distance function from that seed region as described above. The segmentation labels are assigned by determining which distance function has the lowest value for a given pixel. That is, for each label index, I, the set of points corresponding to that region, Li is given by:
Li={(x,y):Di(x,y)≦Dj(x,y),∀j≠i}. (9)
Note that some pixels may be defined with multiple labels using such a definition. However, these pixels are typically pixels that are on the border between two or more regions and can be thus considered border pixels or can be assigned to one region or the other without any loss in the utility of the method.
The middle column 1212 shows the gradient of the smoothed ROI, and the right column 1216 shows the corresponding final segmentation result using the method of distance function locally weighted by the image gradient alone, i.e., ∥∇D∥=∥∇I∥. For each image, one seed region is the rectangular box surrounding the structure and two interior seed regions are chosen, one inside of each region of interest as shown in the first column.
The foregoing Detailed Description is to be understood as being in every respect illustrative and exemplary, but not restrictive, and the scope of the invention disclosed herein is not to be determined from the Detailed Description, but rather from the claims as interpreted according to the full breadth permitted by the patent laws. It is to be understood that the embodiments shown and described herein are only illustrative of the principles of the present invention and that various modifications may be implemented by those skilled in the art without departing from the scope and spirit of the invention. Those skilled in the art could implement various other feature combinations without departing from the scope and spirit of the invention.
This application claims the benefit of U.S. Provisional Application No. 60/724,374 filed Oct. 7, 2005, which is incorporated herein by reference.
Number | Date | Country | |
---|---|---|---|
60724374 | Oct 2005 | US |