Method for deformable registration of images

[0001] This application claims priority from U.S. Provisional application Ser. No. 60/728,224 filed on Oct. 19, 2005, which is incorporated herein by reference.

[0002] This invention relates generally to the registration of images and more particularly to the deformable registration of images.

[0003] As is known in the art, registration of images has a wide range of applications. One application is in medical imaging. Registration of pairs of images (2D or 3D) has been extensively studied for medical images, see for example Maintz, J. B. A., Veirgerver, A survey of medical image registration, Medical Image Analysis, 2(1),1-36, 1998.

[0004] For the most of the approaches, the main assumption is that warping functions should be continuous, smooth and invertible, so that every point in image one (fixed) maps to exactly one point in image two (moving), and vice-versa. Such smooth, invertible functions are known as diffeomorphisms. Diffeomorphism can be enforced (but not guaranteed) through regularization of the dense deformation field assuming of elastic (see Bajscy R., Lieberson R. and Reivich M., A computerized system for the elastic matching of deformed radiographic images to idealized atlas images, J. Comput. Assis. Tomogr., 7:618-625, 1983), viscous fluid (see Christensen G. E., Joshi S. C., and Miller M., volumetric transformation of brain anatomy, IEEE Trans. Medical Image, 16:864-877, 1997, or splines (see Rueckert D., Frangi A. and Schnabel J., Automatic construction of 3D statistical deformation models using non-rigid deformation,. In MICCAI, pages 77-84, 2001 properties. Diffeomorphic transformations maintain topology and guarantee that connected sub-regions remain connected. The main problem arises, when there are no correspondences available. When dealing with medical images, for example, image from abdominal area of a patient taken at two time points, it is foreseeable that there would not be explicit correspondences for all the points within the sets. This in turn causes strong violation of the intensity similarity and failure of methods which are assuming and enforcing these. There are numerous examples for this case. Let us take two sets of MR scans of same patient before and after surgery where a tumor is present in the first and not in the second image. Another example that happens very frequently is the CT images of male pelvic region for prostate cancer therapy. Rectum, bladder, and intestine content change drastically for one therapeutic session to another makes the process of establishing correspondences impossible and more importantly meaningless. For all these cases, any regularization on the erroneous deformation field caused from naive similarity metric enforcing a complete match, results in serious error on parts of the image even places, where the correspondences can be established. Regularizing the flow causes erroneous results to disperse to other parts as well.

[0005] In accordance with the present invention, a method is provided for registration of two image data sets. The method includes compartmentalizing a first one of the two image data sets into a plurality of regions with each one of such regions having a corresponding region in the other one of the two image data sets. For each one of the regions in the first one of the two image data sets and for the corresponding one of the regions in the other one of the two image data set, the method computes an energy function related to the degree such two regions match one another. The method minimizes the sum of the energy functions for each one of the regions in the first one of the two image data sets and for the corresponding one of the regions in the other one of the two image data set by deforming the image data set of such region in the other one of the image data set where the energy functions for each region is defined separately.

[0006] With such method, prior information regarding the parts (i.e., regions) of the images where correspondence and conventional similarity could be violated is obtained. The method uses prior spatial knowledge of such regions on one of the images. The method registers the two image data sets and at the same time propagates the specified region boundaries from one image data set to the other, while trying to preserve the diffeomorphic property of the field all over the image.

[0007] The method deals with this kind of scenarios through a framework that requires rough spatial knowledge of areas, where correspondences cannot be found. The method incorporates a constraint replacing image similarity on parts, where correspondence cannot be established. Using this, the method avoids having penalizing effect on the deformation field on those areas. The constraint could be realized, similar to a segmentation approach, through computation of the probability of the intensity belonging to a certain distribution.

[0008] The details of one or more embodiments of the invention are set forth in the accompanying drawings and the description below. Other features, objects, and advantages of the invention will be apparent from the description and drawings, and from the claims.

[0009] FIG. 1 is a flow diagram of a process used in the registration of deformable images in accordance with the invention;

[0010] FIGS. 2A and 2B are sketches of a baseline image (i.e., a fixed image data set) and a current image (i.e., a moving data set), the method of FIG. 1 being used to register different regions in the moving data set with corresponding regions in the fixed data set in accordance with the invention.

[0011] Like reference symbols in the various drawings indicate like elements.

[0012] In the deformable registration problem, we are given two intensity images I_f (i.e., fixed that is a prior or baseline image) and I_m (i.e., moving that is a current image which is to be deformed to match the fixed image) over the space Ω, where an unknown transformation T:Ω→Ω has to be recovered. In order to solve this problem, there should be a spatial metric available to measure intensity dissimilarity between the two images. Equation 1 shows an energy functional that has to be minimized with respect to the transformation T. $\begin{array}{cc}{E}_{\mathrm{image}}\left(T\right)={\int }_{\Omega }M\left(I^f\left(x\right),I^m\left(T\left(x\right)\right)\right)dx,& \left(1\right)\end{array}$
where x=[x y z]εΩ for three dimensional images, M ( . , . ) estimates the (dis) similarity. Since, the equation 1 is ill-conditioned, we need to consider another set of constraints, which regularize the transformation T. Equation 2 shows an energy functional for the diffusion regularization. $\begin{array}{cc}{E}_{\mathrm{reg}}\left(T\right)={\int }_{\Omega }\mathrm{trace}\left(\left(I-J\left(T\left(x\right)\right)\right){\left(I-J\left(T\left(x\right)\right)\right)}^{\prime }\right)dx& \left(2\right)\end{array}$
Finally, concurrent minimization of the two equations delivers the solution, as follows: $\begin{array}{cc}\hat{T}=\arg \underset{T}{\min }\left(E_{\mathrm{image}}\left(T\right)+\alpha E_{\mathrm{reg}}\left(T\right)\right).& \left(3\right)\end{array}$
where α is a parameter defining the degree of regularization applied over the deformation field and hat (ˆ) denotes the estimate for the variable. If we assume that images are mono-modal, where the brightness constancy constraint holds, we can use optical flow framework for computing the registration parameters. In this framework, since we are concerned with flow, we re-formulate the transformation T as T(x)=x+u, where u is the flow. Furthermore, replacing M with sum of square differences, transforms the equation 3, as follows: $\begin{array}{cc}\begin{array}{c}\hat{u}=\arg \min {\int }_{\Omega }\left({I^f\left(x\right)-I^m\left(x+u\right)}^2\right)dx+\\ \alpha {\int }_{\Omega }\left({\nabla u_x}^2+{\nabla u_y}^2+{\nabla u_z}^2\right)dx\end{array}& \left(4\right)\end{array}$
where u=[u_x, u_y, u_z]^T. Equation 4 is a non-linear one with respect to u. Solution must satisfy a set of Euler-Lagrange equations, which can be solved using an iterative scheme.

[0013] In some application, we need to deal with images, in which for some parts no correspondences can be established. Enforcing both geometrical and radiometrical correspondences, as in it is done in sum of square difference in equation 4 is too penalizing and cause serious errors. In these cases, we need to apply much softer constrains. One feasible constrain is to consider that the intensities of these corresponding parts are belonging to a known probability density function. This is a rather global constraint defined over the specific parts of the image, and has no specificity on the local flow field over that area. Let us assume that Φ_i for iε[0, n−1] are non-overlapping subsets of Ω, where the intensity probability density function, i.e. pdf is defined as p_i. We can then define an energy function that constrains the deformation field based on the either pre-defined pdf or estimated pdf from the pre-defined compartmentalization step, as follows: $\begin{array}{cc}{E}_{\mathrm{seg}}\left(T\right)=-\sum _i{\int }_{{\Phi }_i}\log (p_i\left(I^m\left(x+u\right)\right)dx.& \left(5\right)\end{array}$

[0014] In special case, where pi is a Gaussian distribution, the equation 5 is constraining the flow field in a way that the intensity of the area defined by Φ_i remain close to the mean of the distribution. In one scenario, the pdf for each region can be estimated using the intensity histograms of those regions specified on the I_f, in another scenario pdfs could be known as a priori.

[0015] Here, the method incorporates the spatial soft constraints on the parts of the image, as it is described above into the optical flow framework. The motivation is to compensate for the fact that the brightness constancy constraint and diffeomorphism do not hold on specific parts of the image. Penalizing the flow field to provide accurate correspondences as it is done in equation 4, does have an adverse effect and results in an erroneous mapping. Deformation field can be extracted using the following equation: $\begin{array}{cc}\begin{array}{c}\hat{u}=\arg \underset{u}{\min }{\int }_{\Omega -\bigcup _i{\Phi }_i}\left({I^f\left(x\right)-I^m\left(x+u\right)}^2\right)dx-\\ \beta {\int }_{\bigcup _i{\Phi }_i}\log \left(p_i\left(I_m\left(x+u\right)\right)\right)dx+\\ \alpha {\int }_{\Omega }\left({\nabla u_x}^2+{\nabla u_y}^2+{\nabla u_z}^2\right)dx\end{array}& \left(6\right)\end{array}$
where ∪_iΦ_i denotes the union of Φ_i's, and α, β>0. If p_i is a Gaussian distribution β is one, otherwise, we choose an experimental value close to one. Thus, the second term in equation $\begin{array}{cc}\left(i.e.,\beta {\int }_{\bigcup _i{\Phi }_i}\log \left(p_i\left(I_m\left(x+u\right)\right)\right)dx\right)& \left(6\right)\end{array}$
represents the sum of the energy functions for Ui each one of the regions in the first one of the two image data sets and for the corresponding one of the regions in the other one of the two image data set by deforming the image data set of such region in the other one of the image data set where the energy functions for each region is defined separately.

[0016] According the calculus of variation, the minimizer of the equation in 6 must fulfill the Euler-Lagrange equations: $\begin{array}{cc}\begin{array}{cc}\left(I^f\left(x\right)-I^m\left(x+u\right)\right)I_x^m\left(x+u\right)+\mathrm{\alpha \Delta }{u}_x& \mathrm{for}x\in \Omega -\bigcup _i{\Phi }_i\\ \frac{1}{p_i\left(I_m\left(x+u\right)\right)}{p}_i^{\prime }\left(I_m\left(x+u\right)I_m^x\left(x+u\right)\right)+\mathrm{\alpha \Delta }{u}_x& \mathrm{for}x\in {\Phi }_i\end{array}& \left(7\right)\end{array}$
where u_x, u_y, and u_z are the components of u and Δ denotes Laplacian operation. By changing of the subscripts of I_x^m(x+u) and ∇u_x to y and z, we get additional equations for the three dimensional case. The equation (7) can be solved using a fixed point iteration scheme that continually solves for updated of the deformation field in a multi-resolution setting.

[0017] Referring to FIG. 1, the process is as follows:

[0018] The method obtains the fixed data set and the moving data set, Step 100 and 102. Next, conventional pre-processing (such as smoothing, de-noising, etc) is performed on both the fixed data set and the moving data set, Steps 104 and 106. Next, the fixed data is compartmentalized (i.e., segmenting) into regions by, for example, any known manual or automatic contouring techniques, as for example using a digital pen to outline the regions of interest, Step 108. For example in FIG. 2A one region C₁ having a probability density of p_i(I) may for example be the bladder and another, C₂, probability density of p2 (I) the rectum. The regions C₀ are other things in the image. FIG. 2B shows the two regions (i.e., the rectum and the bladder) in the moving data. The images may be obtained with any conventional imaging equipment, such as CT, MRI, X-ray apparatus, not shown, and the process described herein performed by a computer program stored in a memory of the processor therein.

[0019] Next, multi-resolution pyramid of both fixed and moving data sets are set up, Steps 110 and 112. Next, the process starts from the lowest resolution i until the highest resolution is reached, Step 114. Thus, since the first data set is not at the highest resolution, the process increases the resolution by 1, i.e., i+1. Step 118. The process then determines whether there is convergence and if not, the process is aborted (i.e., stopped). Steps 120 and 130. If there is convergence the moving data is warped according to the k th iteration deformation u_kⁱ to I_mⁱ(x+u_kⁱ), Step 122. Any image transformation method maybe used such as that outlined in the following publications Digital Image Warping, George Wolberg, IEEE Computer Society Press, Los Alamitos, Calif., 1990.

[0020] Next, the process computes the image force for each compartment separately C_j is for the compartment with known intensity distribution of P_j and C₀ for elsewhere and, Step 124 $\begin{array}{cc}f\left(I_m^i,I_f^i,u_k^i\right)=\{\begin{array}{cc}-\left(I_f^i\left(x\right)-I_m^i\left(x+u_k^i\right)\right)\nabla I_m^i\left(x+u_k^i\right)& \mathrm{for}{C}_0\\ \frac{\partial p_j/\partial I_m^i\left(I_m^i\left(x+u_k^i\right)\right){\nabla }_m^i\left(x+u_k^i\right)}{p_j\left(I_m^i\left(x+u_k^i\right)\right)}& \mathrm{for}{C}_j\end{array}& \left(8\right)\end{array}$

[0021] Next, (Step 126) the process solves for deformation to balance force in accordance with:
Au_k+1ⁱ=αf(I_mⁱI_fⁱ,u_kⁱ)
u_k+1ⁱ=αA⁻¹f(I_mⁱ,I_fⁱu_kⁱ)  (9)
where A is a linear operator derived from spatial discretization of the regularization term in equation (7) and the A⁻¹ is the symbolic inverse of the operation. The equation (9) is the *iterative form that is performed within a loop till convergence. A multi-grid based successive over-relaxation method can be employed to compute the inverse operation [William L. Briggs, Van Emden Henson, and Steve F. MacCormick. A multigrid tutorial. SIAM, Society for Industrial and Applied Mathematics, 2. ed. edition, 2000].

[0022] The resulting deformation (initialized with zero) u_kⁱ, Step 128 and the process returns to Step 120 to determine whether there is convergence.

[0023] Thus, considering FIGS. 2A and 2B, for the rectum, for example, and energy function is determined that relates the degree such two regions match one another, i.e., relates the degree to which the image of the rectum in the. fixed image data set matches the rectum in the moving image data set as described above in equation (6).

[0024] The process then deforms the image data set of the rectum in the moving image data sets and by minimizing the energy function. This process is performed concurrently for the other regions, such as the bladder, in this example.

[0025] A number of embodiments of the invention have been described. Nevertheless, it will be understood that various modifications may be made without departing from the spirit and *scope of the invention. Accordingly, other embodiments are within the scope of the following claims.