Method and system for image segmentation by evolving radial basis functions

[0001] This application claims the benefit of U.S. Provisional Application No. 60/865,688, filed Nov. 14, 2006, the disclosure of which is herein incorporated by reference.

[0002] The present invention relates to image segmentation, and more particularly, to image segmentation using radial basis functions (RBFs).

[0003] Since the introduction of active contour methods, they have been successfully used in computer vision applications. Such active contour methods have been used in many important applications, such as image segmentation, 3D scene reconstruction, and object tracking. One reason why active contours have been used so successfully is that they provide a structured approach, via energy minimization, to deform a contour or surface.

[0004] In the context of image segmentation, active contours deform based on various image-based and internal forces so that the contour's edges match object or region boundaries in the image, while maintaining smoothness. The smoothness, or regularization terms, provide robustness to noise while providing a measured approach to handling missing or low-confidence data. A typical application of an active contour will start with an initial contour, which is then iteratively deformed until it converges to a solution that balances the forces acting on the contour. Typically, these forces result from analytical expressions that are derived using variational calculus applied to an energy minimization problem. However, it is possible to define the forces directly without using an energy formulation.

[0005] In the convention active contour techniques, earlier methods represented the contour using a topologically fixed parametric representation, such as a polyline, spline, etc., specified by a fixed number of control points. Such representations are simple and efficient to implement, however, they lack straightforward mechanisms for topological control. Often in segmentation problems, the topology of the problem is unknown a priori, and the contour must break apart or merge during the evolution of the contour. Such topological changes cannot traditionally be performed by active contours using a topologically fixed parametric representation. Although various methods for providing topological changes have been proposed, the implementation of such methods is complicated and not natural to active contour techniques based on parametric representations.

[0006] More recent techniques use implicit active contours, or level-set methods, which represent the contour as a level set of a higher dimensional embedding function. The primary advantage of this representation is that topological changes occur naturally in the evolution of the contour. By manipulating the embedding function, the level set that represents the contour can innately split or merge without requiring any specialized implementation to handle topological changes. However, the embedding function must be updated on a dense set of points, which requires a significant amount of storage, even when using efficient narrowband techniques.

[0007] Recently introduced methods attempt to combine some of the advantages of explicit and implicit active contours through the use of radial basis functions (RBFs), or unstructured point clouds. Such methods model an implicit function defining an active contour as a superposition of RBFs. The RBFs define a set of points (each point representing an RBF), from which the embedding function can be calculated. The points (RBFs) are moved (and correspondingly, the embedding function is updated) in order to deform the active contour to solve a segmentation problem. These approaches have been demonstrated to provide the topological flexibility of level set methods with the low storage requirements of parametric representations, and provide flexibility in terms of RBF placement and interaction.

[0008] The present invention provides a variational method for evolution of radial basis functions (RBFs) for image segmentation. While conventional active contour methods using RBFs start with a force definition rather than deriving forces, embodiments of the present invention use forces derived through an energy formulation to evolve the RBFs. Based on this derivation, embodiments of the present invention can evolve multiple parameters of each RBF, such as position, weight, anisotropy (orientation, standard deviation along orthogonal directions) of each RBF, whereas conventional methods only evolved RBF positions.

[0009] In one embodiment of the present invention, a set of RBFs is initialized to define an implicit function, the zero level set of which defines a contour on an image. The contour is iteratively deformed by updating a weight and multiple parameters of each of the RBFs based on each pixel of the contour, until the contour converges. The parameters of each RBF can include a center point, standard deviation along orthogonal directions, and an orientation. The parameters and the weight of each RBF are updated using a coupled set of differential equations derived based from a region-based or boundary-based image segmentation energy formulation. The final contour at a convergence defines boundaries of a target object or region in the image.

[0010] According to another embodiment of the present invention, RBFs can be merged or added between iterations of updating the parameters of the RBFs. If two RBFs are sufficiently close together, they can be replaced by a single RBF, which is a sum of the two RBFs and is centered at a point halfway between center points of the two RBFs. If there is an area of the image that has a high gradient for the implicit function defined by the RBFs, and there is no existing RBF in that area, an RBF can be added to that area.

[0011] 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.

[0018] The present invention is directed to a method for image segmentation using radial basis functions (RBFs). Embodiments of the present invention are described herein to give a visual understanding of the segmentation method. A digital image is often composed of digital representations of one or more objects (or shapes). The digital representation of an object is often described herein in terms of identifying and manipulating the objects. Such manipulations are virtual manipulations accomplished in the memory or other circuitry/hardware of a computer system. Accordingly, is to be understood that embodiments of the present invention may be performed within a computer system using data stored within the computer system.

[0019] Embodiments of the present invention are directed to evolving a set of radial basis functions in order to segment a region or an object in an image. The set of RBFs define an implicit function, and the zero level set of the implicit function is the active contour. By evolving parameters of the RBFs, the active contour is deformed until it converges such that the contour's edges match object or region boundaries in the image. Embodiments of the present invention use forces derived through an energy formulation to evolve the RBFs. The derivation of these forces results in a set of differential equations that drive the evolution of the RBFs. These equations are derived herein for region-based segmentation problems and edge-based (or boundary-based) segmentation problems. Based on this derivation, it is possible to evolve all parameters of each RBF. Although certain RBFs are used in describing embodiments of the present invention, the present invention is not limited thereto, and can be used with the methodology described herein.

[0020] In region-based segmentation methods, the evolution of the active contour is based on attributes of entire regions in the image. As described here, the evolution of the contour is based on the attribute of intensity, but other attributes, such as intensity, color, texture, probabilities, etc., can also be used. Region-based segmentation methods model the image as being composed of distinct regions, each with its own statistics, which are used to deform the active contour towards the region boundaries. This is notably different from edge-based segmentation methods, where the evolution depends on local image gradients. In this sense, region-based segmentation methods are more global and more robust to noise than edge-based methods.

[0021] A well-known region-based model is the Mumford-Shah function, which models the image regions as piece-wise smooth functions. For example, in the case of two regions, region R inside the contour C and region R^C outside of the contour C,
E(ƒ,C)=∫_R(I(x)−ƒ_R(x))²dx+∫_RC(I(x)−f_RC(x))²dx+∫_Ω/C|∇ƒ|²dx+γ∫_Cds,  (1)
where E(ƒ,C) is the energy of the contour, I(x) is the image intensity at pixel x, f is a piecewise smooth function consisting of regions ƒ_R(x) inside C and ƒ_RC(x) outside C, and the last term is a regularization term weighted by a constant γ. While there are various active contour methods for solving the Mumford-Shah functional, the well-known piecewise constant model of Chan-Vese is a powerful approximation that replaces the piecewise smooth function ƒ with piecewise constant functions in each region. This simplified form can be expressed as:
E(c)=∫_R(I(x)−μ_in)²dx+∫_RC(I(x)−μ_out)²dx+γ∫_Cds,  (2)
where μ_in is the average intensity inside C and μ_out is the average intensity outside C. Other variations of the Chan-Vese function exist as well, for example, using various statistics other than or in addition to the mean.

[0022] In embodiments of the present invention, the contour CεΩ is represented by the zero level set of a Lipschitz function φ: Ω→, such that:
C={xεΩ:φ(x)=0},
R={xεΩ:φ(x)<0},
R^C={xεΩ:φ(x)>0}.

[0023] The embedding function φ(x) can be modeled using N RBFs ψ(x,g_i1 . . . g_iM) as:

[0024] $\begin{array}{cc}\varphi \left(x\right)=P\left(x\right)+\sum _{i=1}^N{w}_i{\psi }_i\left(x,g_{i1}\dots g_{\mathrm{iM}}\right),& \left(3\right)\end{array}$
where ψ(x, g_i1 . . . g_iM) is the ^th RBF parameterized by M variables g_i1 . . . g_iM, w_i is the weight of the ^th RBF, and P(x) is a polynomial term that spans the null space of the RBF.

[0025] Using the Heaviside function,

[0026] $\begin{array}{cc}H\left(z\right)=\{\begin{array}{cc}1,& z\ge 0\\ 0,& \mathrm{otherwise}\end{array}& \left(4\right)\end{array}$
and the Dirac delta function

[0027] $\delta \left(z\right)=\frac{d}{dz}H\left(z\right),$
the energy of Equation (2) can be expressed as

[0028] $\begin{array}{cc}E\left(\varphi \left(w,g\right)\right)={\int }_{\Omega }{\left(I\left(x\right)-{\mu }_{\mathrm{in}}\right)}^{2}H\left(\varphi \left(x\right)\right)dx+{\int }_{\Omega }{\left(l(x)-{\mu }_o\right)}^2\left(1-H\left(\varphi \left(x\right)\right)\right)dx+\gamma {\int }_{\Omega }\delta \left(\varphi \left(x\right)\right)\nabla \varphi \left(x\right)dx,& \left(5\right)\end{array}$
where w=[w₁ . . . w_N]^T, and g=[g_ij], where i=1 . . . N, j=1 . . . M.

[0029] Under this formulation of the problem, the variation of the energy with respect to the RBF parameter g_ij and the RBF weight w_i the following expression can be derived:

[0030] $\begin{array}{cc}\frac{\partial E}{\partial g_{\mathrm{ij}}}={\int }_{\Omega }{\left(I\left(x\right)-{\mu }_{\mathrm{in}}\right)}^2\delta \left(\varphi \left(x\right)\right)\frac{\partial \varphi }{\partial g_{\mathrm{ij}}}dx-{\int }_{\Omega }{\left(I\left(x\right)-{\mu }_{\mathrm{out}}\right)}^2\delta \left(\varphi \left(x\right)\right)\frac{\partial \varphi }{\partial g_{\mathrm{ij}}}dx+\gamma {\int }_{\Omega }\mathrm{div}\left(\frac{\nabla \varphi }{\nabla \varphi }\right)\delta \left(\varphi \left(x\right)\right)\frac{\partial \varphi }{\partial g_{\mathrm{ij}}}dx,& \left(6\right)\end{array}$
which simplifies to the contour integral,

[0031] $\begin{array}{cc}\begin{array}{c}\frac{\partial g_{\mathrm{ij}}}{\partial t}=\frac{\partial E}{\partial g_{\mathrm{ij}}}\\ ={\int }_C\left[{\left(I\left(x\right)-{\mu }_{\mathrm{in}}\right)}^2-{\left(I\left(x\right)-{\mu }_{\mathrm{out}}\right)}^2+\gamma \mathrm{div}\left(\frac{\nabla \varphi }{\nabla \varphi }\right)\right]\frac{\partial \varphi }{\partial g_{\mathrm{ij}}}dC,\end{array}& \left(7\right)\end{array}$
and similarly,

[0032] $\begin{array}{cc}\begin{array}{c}\frac{\partial w_i}{\partial t}=\frac{\partial E}{\partial w_{\mathrm{ij}}}\\ ={\int }_C\left[{\left(I(x)-{\mu }_{\mathrm{in}}\right)}^2-{\left(I\left(x\right)-{\mu }_{\mathrm{out}}\right)}^2+\gamma \mathrm{div}\left(\frac{\nabla \varphi }{\nabla \varphi }\right)\right]\frac{\partial \varphi }{\partial w_{\mathrm{ij}}}dC.\end{array}& \left(8\right)\end{array}$

[0033] These expressions take the form

[0034] ${\int }_C\frac{\partial E}{\partial \varphi }\frac{\partial \varphi }{\partial p}dC,$
where p is g_ij or w_i. This form results from the functional composition that resulted in application of the chain rule. Conceptually, these equations state that to determine the update of a parameter of an RBF, the contour C must be traversed, accumulating gradients at each point That is, each point on the zero level set contributes to the update of the RBF parameters. Accordingly, unlike conventional level set update equations, the integral in Equations (7) and (8) combines measurements from all points along the zero level set, providing increased robustness to noise. Based on Equations (7) and (8), it is possible to update all parameters of each RBF, including the weights, position, orientation, anisotropy, etc. In such a case, each iteration of evolving all RBFs is order O(LMN), where L is the number of pixels on the zero level set, and N and M are the number of RBFs and the number of parameters for each RBF, respectively, as described above.

[0035] As described above, the energy minimization results in the form

[0036] ${\int }_C\frac{\partial E}{\partial \varphi }\frac{\partial \varphi }{\partial p}dC.$
In general, this relationship can be utilized to transform any level set flow or curve evolution into an RBF evolution in this framework. For boundary-based image segmentation, a similar expression can be deduced for the geodesic flow,

[0037] $\begin{array}{cc}\frac{\partial g_{\mathrm{ij}}}{\partial t}=\frac{\partial E}{\partial g_{\mathrm{ij}}}={\int }_C\left[F\kappa \nabla \varphi -\nabla F·\nabla \varphi \right]\frac{\partial \varphi }{\partial g_{\mathrm{ij}}}dC& \left(9\right)\\ \frac{\partial w_{\mathrm{ii}}}{\partial t}=\frac{\partial E}{\partial w_i}={\int }_C\left[F\kappa \nabla \varphi -\nabla F·\nabla \varphi \right]\frac{\partial \varphi }{\partial w_i}dC,& \left(10\right)\end{array}$
where F(I(x)) is a function of an edge detector response, such as

[0038] $F=\frac{1}{1+{\nabla I}^2},$
and κ is the curvature of the active contour at point x. It can be noted that with boundary-based segmentations, it may be helpful to run a GVF diffusion on the vector field ∇F before evolving the active contour.

[0039] As described above, basic equations are derived for evolution of RBFs to achieve either region-based or boundary-based image segmentation. However, the complete derivation depends on the RBFs chosen and their derivatives. According to an embodiment of the present invention, anisotropic Gaussian RBFs can be used for region-based or boundary-based image segmentation. However, the present invention is not limited to anisotropic Gaussian RBFS, and any RBF with an analytic derivative can be similarly used applying the framework described herein. The use of anisotropic RBFs may be advantageous because they are effective for approximating sharp corners, and they are a more general form of isotropic RBFs. Although 2D RBFs are used in the embodiments of the present invention described herein, it is to be understood that 3D RBFs can be used as well.

[0040] The equation for a 2D anisotropic Gaussian centered at point (c_x,c_y), standard deviation σ_x,σ_y, and orientation angle θ has the function form:

[0041] $\begin{array}{cc}\begin{array}{c}\psi \left(x,y\right)=\psi \left(x,y\right)\exp [\frac{-1}{2o_x^2{\sigma }_y^2}({a_1\left(x-c_x\right)}^2-\hspace{1em}\hspace{1em}\\ 2a_2\left(x-\hspace{1em}\hspace{1em}{c}_x\right)\left(y-c_y\right)+{a_3\left(y-c_y\right)}^2)],\end{array}& \left(11\right)\\ \mathrm{where}& \\ a_1={\sigma }_x^2{\cos }^2\theta +{\sigma }_y^2{\sin }^2\theta & \left(12\right)\\ a_2=\left({\sigma }_y^2-{\sigma }_x^2\right)\cos \theta \sin \theta & \left(13\right)\\ a_3={\sigma }_x^2{\sin }^2\theta +{\sigma }_y^2{\cos }^2\theta .& \left(14\right)\end{array}$

[0042] This RBF is parameterized by the M=5 parameters, g_j=[c_x,c_y, σ_x,σ_y,θ]^T. To implement the RBF flow (evolution), the derivative of the RBF must be determined with respect to each of these parameters. Doing so yields the following equations:

[0043] $\begin{array}{cc}\frac{\partial \psi }{\partial c_x}=\psi ·\left(\frac{a_{1}X-a_{2}Y}{{\sigma }_x^2{\sigma }_y^2}\right)& \left(15\right)\\ \frac{\partial \psi }{\partial c_y}=\psi ·\left(\frac{-a_{2}X+a_{3}Y}{{\sigma }_x^2{\sigma }_y^2}\right)& \left(16\right)\\ \frac{\partial \psi }{\partial \theta }=\psi ·\left(-\frac{X^2\frac{\partial a_1}{\partial \theta }-2\mathrm{XY}\frac{\partial a_2}{\partial \theta }+Y^2\frac{\partial a_3}{\partial \theta }}{2{\sigma }_x^2{\sigma }_y^2}\right)& \left(17\right)\\ \frac{\partial \psi }{\partial {\sigma }_x}=-\frac{\psi }{{\sigma }_x}+\psi ·\left[\begin{array}{c}\frac{a_1{X}^2-2a_2\mathrm{XY}+a_3{Y}^2}{{\sigma }_x^3{\sigma }_y^2}-\\ \frac{1}{2{\sigma }_x^2{\sigma }_y^2}\left(X^2\frac{\partial a_1}{\partial {\sigma }_x}-2\mathrm{XY}\frac{\partial a_2}{\partial {\sigma }_x}+Y^2\frac{\partial a_3}{\partial {\sigma }_x}\right)\end{array}\right]& \left(18\right)\\ \frac{\partial \psi }{\partial {\sigma }_y}=-\frac{\psi }{{\sigma }_y}+\psi ·\left[\begin{array}{c}\frac{a_1{X}^2-2a_2\mathrm{XY}+a_3{Y}^2}{{\sigma }_x^2{\sigma }_y^3}-\\ \frac{1}{2{\sigma }_x^2{\sigma }_y^2}\left(X^2\frac{\partial a_1}{\partial {\sigma }_y}-2\mathrm{XY}\frac{\partial a_2}{\partial {\sigma }_y}+Y^2\frac{\partial a_3}{\partial {\sigma }_y}\right)\end{array}\right]& \left(19\right)\\ \mathrm{where}& \\ X=\left(x-c_x\right)& \left(20\right)\\ Y=\left(x-c_y\right)& \left(21\right)\\ \frac{\partial a_1}{\partial \theta }=-2\left({\sigma }_x^2-{\sigma }_y^2\right)\cos \mathrm{\theta sin}\theta & \left(22\right)\\ \frac{\partial a_2}{\partial \theta }=-2\left({\sigma }_x^2-{\sigma }_y^2\right)\left({\cos }^2\theta -{\sin }^2\theta \right)& \left(23\right)\\ \frac{\partial a_3}{\partial \theta }=-2\left({\sigma }_x^2-{\sigma }_y^2\right)\cos \mathrm{\theta sin}\theta & \left(24\right)\\ \frac{\partial a_1}{\partial {\sigma }_x}=2{\sigma }_x{\cos }^2\theta & \left(25\right)\\ \frac{\partial a_2}{\partial {\sigma }_x}=-2{\sigma }_x\cos \theta \sin \theta & \left(26\right)\\ \frac{\partial a_3}{\partial {\sigma }_x}=2{\sigma }_x{\sin }^2\theta & \left(27\right)\\ \frac{\partial a_1}{\partial {\sigma }_y}=2{\sigma }_y{\sin }^2\theta & \left(28\right)\\ \frac{\partial a_2}{\partial {\sigma }_y}=2{\sigma }_y\cos \mathrm{\theta sin}\theta & \left(29\right)\\ \frac{\partial a_3}{\partial {\sigma }_y}=2{\sigma }_y{\cos }^2\theta .& \left(30\right)\end{array}$

[0044] Using the above equations (15-30) an equation for each parameter (and weight) of an anisotropic Gaussian RBF can expressed in the form derived above (e.g., equations (7) and (8)). Accordingly, the following is a set of coupled differential equations that drive the evolution of the i^th RBF:

[0045] $\begin{array}{cc}\frac{d{c}_{\mathrm{ix}}}{dt}={\int }_C{\mathrm{Dw}}_i{\psi }_i·\left(\frac{a_{i1}{X}_i-a_{i2}{Y}_i}{{\sigma }_{\mathrm{ix}}^2{\sigma }_{\mathrm{iy}}^2}\right)dC& \left(31\right)\\ \frac{d{c}_{\mathrm{iy}}}{dt}={\int }_C{\mathrm{Dw}}_i{\psi }_i·\left(\frac{-a_{i2}{X}_i+a_{i3}{Y}_i}{{\sigma }_{\mathrm{ix}}^2{\sigma }_{\mathrm{iy}}^2}\right)dC& \left(32\right)\\ \frac{d{w}_i}{dt}={\int }_{C}D{\psi }_{i}dC& \left(33\right)\\ \frac{d{\theta }_i}{dt}={\int }_C{\mathrm{Dw}}_{\mathrm{Xi}}·{\psi }_i·\left(-\frac{X_i^2\frac{\partial a_{i1}}{\partial {\theta }_i}-2X_i{Y}_i\frac{\partial a_{i2}}{\partial {\theta }_i}+Y_i^2\frac{\partial a_{i3}}{\partial {\theta }_i}}{2{\sigma }_{\mathrm{ix}}^2{\sigma }_{\mathrm{iy}}^2}\right)dC& \left(34\right)\\ \frac{d{\sigma }_{\mathrm{ix}}}{dt}={\int }_C{\mathrm{Dw}}_i·\left\{-\frac{{\psi }_i}{{\sigma }_{\mathrm{ix}}}+{\psi }_i·\left[\frac{a_{i1}{X}_i-2a_{i2}{X}_i{Y}_i+a_{i3}{Y}_i^2}{{\sigma }_{\mathrm{ix}}^3{\sigma }_{\mathrm{iy}}^2}-\frac{1}{2{\sigma }_{\mathrm{ix}}^2{\sigma }_{\mathrm{iy}}^2}\left(X_i^2\frac{\partial a_{i1}}{\partial {\sigma }_{\mathrm{ix}}}-2X_i{Y}_i\frac{\partial a_{i2}}{\partial {\sigma }_{\mathrm{ix}}}+Y_i^2\frac{\partial a_{i3}}{\partial {\sigma }_{\mathrm{ix}}}\right)\right]\right\}dC& \left(35\right)\\ \frac{d{\sigma }_{\mathrm{iy}}}{dt}={\int }_C{\mathrm{Dw}}_i·\left\{-\frac{{\psi }_i}{{\sigma }_{\mathrm{iy}}}+{\psi }_i·\left[\frac{a_{i1}{X}_i^2-2a_{i2}{X}_i{Y}_i+a_{i3}{Y}_i^2}{{\sigma }_{\mathrm{ix}}^2{\sigma }_{\mathrm{iy}}^3}-\frac{1}{2{\sigma }_{\mathrm{ix}}^2{\sigma }_{\mathrm{iy}}^2}\left(X_i^2\frac{\partial a_{i1}}{\partial {\sigma }_{\mathrm{iy}}}-2X_i{Y}_i\frac{\partial a_{i2}}{\partial {\sigma }_{\mathrm{iy}}}+Y_i^2\frac{\partial a_{i3}}{\partial {\sigma }_{\mathrm{iy}}}\right)\right]\right\}dC& \left(36\right)\\ \mathrm{where}& \\ X_i=\left(x-c_{\mathrm{ix}}\right)& \left(37\right)\\ Y_i=\left(y-c_{\mathrm{iy}}\right)& \left(38\right)\end{array}$
and for region-based segmentation,

[0046] $\begin{array}{cc}D=\left[{\left(I\left(x\right)-{\mu }_{\mathrm{in}}\right)}^2-{\left(I\left(x\right)-\mu \mathrm{out}\right)}^2+\gamma \mathrm{div}\left(\frac{\nabla \varphi }{\nabla \varphi }\right)\right]& \left(39\right)\end{array}$
while for boundary-based segmentation,
D=[Fκ|∇φ|−∇F·∇φ].  (40)

[0047] In order to utilize the above described systems of equations in an image segmentation method, these equations can be readily implemented by one skilled in the art as computer executable instructions on a computer.

[0048] Although the derivation of the set of differential equations to drive an anisotropic Gaussian RBF is described above, a similar derivation can be performed for any RBF interpolant that has an analytic derivative, including for example, multi-order RBFs and Wedland's RBFs. It may be advantageous for image segmentation to use RBFs that decay from their center location.

[0049] FIG. 1 illustrates an example of image segmentation by evolving RBFs using the RBF evolution equations (31)-(36) derived above according to an embodiment of the present invention. As illustrated in FIG. 1, a white square having sharp corners is segmented from a dark background. Images (a)-(c) of FIG. 1 show region-based segmentation using the RBF evolution equations derived above, and images (e)-(g) show boundary-based images segmentation using the RBF evolution equations derived above. Images (a) and (e) show initialization of a contour 102, which is the zero level set of the implicit function φ(x, y). The implicit function φ(x, y) is defined as a summation of RBFs 104 and 106. In FIG. 1, reference numeral 104 refers to all negative weighted RBFs, and reference numeral 106 refers to all positive weighted RBFs. In images (a) and (e), negative weighted RBFs 104 are located in the center of the image, and positive weighted RBFs 106 are located around the edges of the image. The ellipse around each of the RBFs 104 and 106 is a visualization of the anisotropy of each of RBFs 104 and 106. Initially (in images (a) and (e)), the RBFs 104 and 106 are all isotropic and are rendered as circles. Images (b) and (f) show an intermediate stage of the segmentation methods, and images (c) and (g) show final segmentation results after the methods have converged. As shown in images (b), (c), (f), and (g), the RBFs 104 and 106 changed their positions, weights, and anisotropy in order to deform the contour 102 to segment the white square. These parameters of the RBFs are evolved using the RBF evolution equations (31)-(36) derived above. Images (d) and (h) respectively show results of conventional region-based and boundary-based image segmentation techniques which evolve only the weights and locations of isotropic RBFs. As shown in images (d) and (h), the conventional methods are do not accurately capture the sharp corners of the square.

[0050] FIG. 2 illustrates a method of segmenting an object in an image according to an embodiment of the present invention. This method utilizes the RBF evolution equations derived above to evolve multiple parameters of RBFs in order to deform a contour such that the contour matches boundaries of the target object. This method can be performed using region-based or boundary based image segmentation.

[0051] As illustrated in FIG. 2, at step 202, a set of RBFs is initialized. In order to initialize the RBFs, initial RBFs are distributed around the image. The initial RBFs define an initial contour for the image segmentation method. The RBFs may be distributed automatically, for example, by distributing positive weighted RBFs around the image border or in a large ring around the image center and negative weighted RBFs in the center of the image. This initial distribution is shown in images (a) and (e) of FIG. 1. It is also possible that the RBFs be distributed manually. The initialization is flexible, as the initial contour can be inside, outside, or both inside and outside of the object being segmented.

[0052] At step 204, an implicit function is calculated based on the set of RBFs. The implicit function is calculated as a summation of the set of RBFs. This is expressed in Equation (3). At step 206, the initial contour is determined by calculating the zero level set of the implicit function. The zero level set of the implicit function is all points (pixels) at which the implicit function is equal to zero. Accordingly, the initial contour is a contour consisting of pixels for which the summation of all of the RBFs is equal to zero.

[0053] At step 208, each parameter and weight of each RBF is updated based on all points on the contour. At this step, the method cycles through each RBF and updates all the parameters of each RBF using the RBF evolution equations derived above. For example, for anisotropic Gaussian RBFs, the center point, anisotropy, and weight of each RBF are updated for each RBF. The anisotropy for Gaussian RBFs can be defined by an orientation, and standard deviation in orthogonal directions. The center point is updated using Equations (31) and (32), the weight is updated using Equation (33), the orientation is updated using Equation (34) and the standard deviation is updated using Equations (35) and (36). As described above, these equations can be used for region-based segmentation (Equation (39)) or boundary-based segmentation (Equation (40)). The equations for updating each of the parameters and the weight of each RBF integrate over the contour, thus all points (pixels) on the contour are used to update each parameter of each RBF. The updates to each parameter of each RBF changes the position, eight, orientation, and anisotropy of each RBF, thus evolving each RBF.

[0054] At step 210, an updated implicit function is calculated based on the updated RBFs. The updated implicit function is the summation of the updated RBFs. At step 212, the contour is deformed by calculating the zero level set of the updated implicit function. Accordingly, the evolution of the RBFs causes the contour to deform towards the boundaries of the target object to be segmented.

[0055] At step 214, RBFs that are sufficiently close together are merged into a single RBF. For two RBFs, it is determined whether the RBFs are sufficiently close together, and if so, they are combined and replaced by a new RBF. In the case of a merge, let the i^th RBF be parameterized by └c_ix,c_iy,σ_ix,σ_iy,θ_i┘ and weight w_i, and the j^th RBF by └c_jx,c_jy,σ_jx,σ_jy,θ_j┘ and weight w_j, i≠j. If the distance between two RBFs, d_ij=√{square root over ((c_ix−c_jx)²+(c_iy−c_jy)²)}{square root over ((c_ix−c_jx)²+(c_iy−c_jy)²)} becomes less than a threshold T_M, then the RBFs are combined by deleting them and replacing them with a new RBF,

[0056] $\left[\frac{1}{2}\left(c_{\mathrm{ix}}+c_{\mathrm{jx}}\right),\frac{1}{2}\left(c_{\mathrm{iy}}+c_{\mathrm{jy}}\right),{\sigma }_{\mathrm{ix}}+{\sigma }_{\mathrm{jx}},{\sigma }_{\mathrm{iy}}+{\sigma }_{\mathrm{jy}},{\theta }_i+{\theta }_j\right]$
and weight w_i+w_j. This new RBF is the sum of the two RBFs being merged and is centered halfway between the two RBFs being merged.

[0057] At step 216, RBFs are added in regions of the image where the gradient of the implicit function is high and there are no existing RBFs. In order to avoid adding new RBFs too close to an existing RBF, the method utilizes a function

[0058] $A\left(x,y\right)=\frac{\nabla \varphi \left(x,y\right)}{\max \left(\varphi \left(x,y\right)\right)}·S\left(x,y\right),\mathrm{where}\frac{\nabla \varphi \left(x,y\right)}{\max \left(\varphi \left(x,y\right)\right)}$
is the normalized gradient of the implicit function and S(x, y) is a splat buffer formed as:

[0059] $\begin{array}{cc}S\left(x,y\right)=\prod _i\left(1-e^{-\lfloor {\left(c_{\mathrm{ix}}-x\right)}^2+{\left(c_{\mathrm{iy}}-y\right)}^2\rfloor /\left(2{\sigma }^2\right)}\right)& \left(41\right)\end{array}$
where σ is a standard deviation of a 2D Gaussian function. The splat buffer is close to one where there are no RBFs, and decreases to zero near RBF centers. Therefore, function A(x, y) is large where there are no RBFs and the gradient of the implicit function is high. The value of A(x, y) is calculated over the image, and an RBF is added at any position (x, y) where A(x, y) has a value above a threshold T_A. This generates a new RBF that is not too close to an existing RBF, yet is located at a point of high gradient in the implicit function. The threshold T_A prevents the method from adding too many constraints where they are not needed. The weight of the newly added isotropic RBF is initially zero, and subsequent iterations will update the weight and anisotropy of the new RBF.

[0060] At step 218 it is determined whether the contour has converged. According to an embodiment of the present invention, it is possible to determine convergence of the contour by monitoring the energy E(C) of the contour. For example, in region-based segmentation, the energy E(C) given by Equation (2). The energy E(C) for the contour can be calculated at each iteration and compared to the energy E(C) calculated at the previous iteration. If the difference is less than a threshold, it can be determined that the contour has converged. If the contour has not converged, the method returns to step 208 and updates the parameters of each RBF again, thus deforming the contour again. This is repeated until the contour converges. If the contour has converged, the method proceeds to step 220.

[0061] At step 220, the final contour is output as the boundaries of the target object. For example, the final contour can be displayed on the image to segment the target object, or stored in memory or storage of a computer system.

[0062] FIG. 3 illustrates exemplary image segmentation results in which RBFs are merged and added. As illustrated in FIG. 3, two white circles are segmented from the dark background using the method of FIG. 2. Image (a) shows initialization of the RBFs and the initial contour, image (b) shows an intermediate stage of the segmentation method, and image (c) shows a final contour at convergence. As shown in image (b), the location, orientation, anisotropy, and weights of the RBFs in image (a) have been updated. Furthermore, new RBFs 302 were added in image (b) that were not in image (a). The new RBFs 302 were added as described above in step 216. Image (c) shows that RBFs of image (b) have been merge into a single RBF 304 in image (c). The RBFs were merged as described above in step 214. Image (c) also illustrates a contour that underwent a topological change from the initial contour in image (a) to the final contour in image (c). The final contour split to separately segment both of circles.

[0063] FIG. 4 illustrates exemplary image segmentation results of a noisy image. As illustrated in FIG. 4, a light circle is segmented from a dark background using the method of FIG. 2 with a region-based active contour. In the images of FIG. 4, Gaussian white noise of unit standard deviation was added to the image to test the method's robustness to noise. Image (a) shows initialization of the RBFs and the initial contour, image (b) shows an intermediate stage of the segmentation method, and image (c) shows a final contour at convergence. As shown by images (a), (b), and (c), the square is successfully segmented despite the low signal to noise ratio.

[0064] FIG. 5 illustrates exemplary image segmentation results of an ultrasound fetal image. As illustrated in FIG. 5, the method of FIG. 2 can be applied to medical images, such as ultrasound images to segment anatomic structures in medical images The images of FIG. 5 show segmentation of a fetal structure in an ultrasound image. Image (a) shows initialization of the RBFs and the initial contour, image (b) shows an intermediate stage of the segmentation method, and image (c) shows a final contour that has converged to the boundaries of the fetal structure.

[0065] The above-described methods for image segmentation by evolving RBFs may be implemented on a computer using well-known computer processors, memory units, storage devices, computer software, and other components. A high level block diagram of such a computer is illustrated in FIG. 6. Computer 602 contains a processor 604 which controls the overall operation of the computer 602 by executing computer program instructions which define such operation. The computer program instructions may be stored in a storage device 612 (e.g., magnetic disk) and loaded into memory 610 when execution of the computer program instructions is desired. Thus, an image segmentation application may be defined by the computer program instructions stored in the memory 610 and/or storage 612 and controlled by the processor 604 executing the computer program instructions. The computer 602 also includes one or more network interfaces 606 for communicating with other devices via a network. The computer 602 also includes other input/output devices 608 that enable user interaction with the computer 602 (e.g., display, keyboard, mouse, speakers, buttons, etc.) One skilled in the art will recognize that an implementation of an actual computer could contain other components as well, and that FIG. 6 is a high level representation of some of the components of such a computer for illustrative purposes.

[0066] 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 Descriptions 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.