Method and apparatus for image registration using large deformation diffeomorphisms on a sphere

[0002] This is a continuation-in-part of prior U.S. patent application Ser. No. 09/186,359 filed Nov. 5, 1998, now U.S. Pat. No. 6,226,418. This continuation-in-part claims priority to U.S. patent application Ser. No. 09/186,359 and incorporates the same by reference herein. This continuation-in-part also claims priority to U.S. provisional application 60/155,141 filed Sep. 22, 1999 (now expired) and incorporates the same by reference herein.

[0004] The present invention relates to image processing systems and methods, and more particularly to image registration systems that combine two or more images into a composite image; in particular the fusion of anatomical manifold-based knowledge with volume imagery via large deformation mapping which supports both kinds of information simultaneously, as well as individually, and which can be implemented on a rapid convolution FFT based computer system.

[0005] Image registration involves combining two or more images, or selected points from the images, to produce a composite image containing data from each of the registered images. During registration, a transformation is computed that maps related points among the combined images so that points defining the same structure in each of the combined images are correlated in the composite image.

[0006] Currently, practitioners follow two different registration techniques. The first requires that an individual with expertise in the structure of the object represented in the images label a set of landmarks in each of the images that are to be registered. For example, when registering two MRI images of different axial slices of a human head, a physician may label points, or a contour surrounding these points, corresponding to the cerebellum in two images. The two images are then registered by relying on a known relationship among the landmarks in the two brain images. The mathematics underlying this registration process is known as small deformation multi-target registration.

[0007] In the previous example of two brain images being registered, using a purely operator-driven approach, a set of N landmarks identified by the physician, represented by x_i, where i=1 . . . N, are defined within the two brain coordinate systems. A mapping relationship, mapping the N points selected in one image to the corresponding N points in the other image, is defined by the equation u(x_i)=k_i, where i=1 . . . N. Each of the coefficients, k_i, is assumed known.

[0008] The mapping relationship u(x) is extended from the set of N landmark points to the continuum using a linear quadratic form regularization optimization of the equation: $\begin{array}{cc}u=\arg \underset{u}{\min }\int {\Vert \mathrm{Lu}\Vert }^2& \left(1\right)\end{array}$

[0009] subject to the boundary constraints u(x_i)=k_i′. The operator L is a linear differential operator. This linear optimization problem has a closed form solution. Selecting L=α∇²+β∇(∇·) gives rise to small deformation elasticity.

[0010] For a description of small deformation elasticity see S. Timoshenko, Theory of Elasticity,McGraw-Hill, 1934 (hereinafter referred to as Timoshenko) and R. L. Bisplinghoff, J. W. Marr, and T. H. H. Pian, Statistics of Deformable Solids,Dover Publications, Inc., 1965 (hereinafter referred to as Bisplinghoff). Selecting L=∇²gives rise to a membrane or Laplacian model. Others have used this operator in their work, see e.g., Amit, U. Grenander, and M. Piccioni, “Structural image restoration through deformable templates,” J. American Statistical Association.86(414):376-387, June 1991, (hereinafter referred to as Amit) and R. Szeliski, Bayesian Modeling of Uncertainty in Low-Level Vision,Kluwer Academic Publisher, Boston, 1989 (hereinafter referred to as Szeliski) (also describing a bi-harmonic approach). Selecting L=∇⁴gives a spline or biharmonic registration method. For examples of applications using this operator see Grace Wahba, “Spline Models for Observational Data,” Regional Conference Series in Applied Mathematics. SIAM, 1990, (hereinafter referred to as Whaba) and F. L. Bookstein, The Measurement of Biological Shape and Shape Change,volume 24, Springer-Verlag: Lecture Notes in Biomathematics, New York, 1978 (hereinafter referred to as Bookstein).

[0011] The second currently-practiced technique for image registration uses the mathematics of small deformation multi-target registration and is purely image data driven. Here, volume based imagery is generated of the two targets from which a coordinate system transformation is constructed. Using this approach, a distance measure, represented by the expression D(u), represents the distance between a template T(x) and a target image S(x). The optimization equation guiding the registration of the two images using a distance measure is: $\begin{array}{cc}u=\arg \underset{u}{\min }\int {\Vert \mathrm{Lu}\Vert }^2+D\left(u\right)& \left(2\right)\end{array}$

[0012] The distance measure D(u) measuring the disparity between imagery has various forms, e.g., the Gaussian squared error distance ∫|T(h(x))−S(x)|²dx, a correlation distance, or a Kullback Liebler distance. Registration of the two images requires finding a mapping that minimizes this distance.

[0013] Other fusion approaches involve small deformation mapping coordinates x ε Ω of one set of imagery to a second set of imagery. Other techniques include the mapping of predefined landmarks and imagery, both taken separately such as in the work of Bookstein, or fused as covered via the approach developed by Miller-Joshi-Christensen-Grenander the '212 patent described in U.S. Pat. No. 6,009,212 (hereinafter “the '212 patent”) herein incorporated by reference mapping coordinates x ε Ω of one target to a second target. The existing state of the art for small deformation matching can be stated as follows:

[0014] Small Deformation Matching: Construct h(x)=x−u(x) according to the minimization of the distance D(u) between the template and target imagery subject to the smoothness penalty defined by the linear differential operator L: $\begin{array}{cc}h(·)=·-u(·)\mathrm{where}\hat{u}=\arg \underset{u}{\min }\int {\Vert \mathrm{Lu}\Vert }^2+\sum _{i=1}^{N}D\left(u\right).& \left(3\right)\end{array}$

[0015] The distance measure changes depending upon whether landmarks or imagery are being matched.

[0016] 1. Landmarks alone. The first approach is purely operator driven in which a set of point landmarks x_i, i=1, . . . N are defined within the two brain coordinate systems by, for example, an anatomical expert, or automated system from which the mapping is assumed known: u(x_i)=k_i, i=1, . . . , N. The field u(x) specifying the mapping h is extended from the set of points {x_i} identified in the target to the points {y_i} measured with Gaussian error co-variances Σ_i: $\begin{array}{cc}\hat{u}=\arg \underset{u}{\min }\int {\Vert \mathrm{Lu}\Vert }^2+\sum _{i=1}^N{\left(y_i-x_i-u\left(x_i\right)\right)}^t{{\Sigma }^{-1}\left(y_i-x_i-u\left(x_i\right)\right)}^T& \left(4\right)\end{array}$

[0017] Here (·)^τdenotes transpose of the 3×1 real vectors, and L is a linear differential operator giving rise to small deformation elasticity (see Timoshenko and Bisplinghoff), the membrane of Laplacian model (see Amit and Szeliski), bi-harmonic (see Szeliski), and many of the spline methods (see Wahba and Bookstein). This is a linear optimization problem with closed form solution.

[0018] 2. The second approach is purely volume image data driven, in which the volume based imagery is generated of the two targets from which the coordinate system transformation is constructed. A distance measure between the two images being registered I₀, I₁is defined as D(u)=∫|I₀(x−u(x))−I₁(x)|²dx. The corresponding optimization is: $\begin{array}{cc}h(·)=·-u(·)\mathrm{where}u=\arg \underset{u}{\min }\int {\Vert \mathrm{Lu}\Vert }^2+\int {\left|I_0\left(x-u\left(x\right)\right)-I_1\left(x\right)\right|}^{2}dx.& \left(5\right)\end{array}$

[0019] The data function D(u) measures the disparity between imagery and has various forms. Other distances are used besides the Gaussian squared error distance, including correlation distance, Kullback Liebler distance, and others.

[0020] 3. The algorithm for the transformation of imagery I₀into imagery I₁has landmark and volume imagery fused in the small deformation setting as in the '212 patent. Both sources of information are combined into the small deformation registration: $\begin{array}{cc}u=\arg \underset{u}{\min }\int {\Vert \mathrm{Lu}\Vert }^2+D\left(u\right)+\sum _{i=1}^N\frac{{\Vert y_i-x_i-u\left(x_i\right)\Vert }^2}{\sigma \underset{1}{2}}& \left(6\right)\end{array}$

[0021] Although small deformation methods provide geometrically meaningful deformations under conditions where the imagery being matched are small, linear, or affine changes from one image to the other. Small deformation mapping does not allow the automatic calculation of tangents, curvature, surface areas, and geometric properties of the imagery. To illustrate the mapping problem FIG. 9 shows an oval template image with several landmarks highlighted. FIG. 10 shows a target image that is greatly deformed from the template image. The target image is a largely deformed oval that has been twisted. FIG. 11 shows the results of image matching when the four corners are fixed, using small deformation methods based on static quadratic form regularization. These figures illustrate the distortion which occurs with small deformation linear mapping when used with landmark points which define a motion corresponding to large deformation. As can be seen in FIG. 11, landmarks defined in the template image often map to more than one corresponding point in the target image.

[0022] Large deformation mapping produces maps for image registration in which the goal is to find the one-to-one, onto, invertible, differentiable maps h (henceforth termed diffeomorphisms) from the coordinates x ε Ω of one target to a second target under the mapping

[0023] h:x→h(x)=x−u(x),xεΩ  (7)

[0024] To accommodate very fine variations in anatomy the diffeomorphic transformations constructed are of high dimensions having, for example a dimension greater than 12 of the Affine transform up-to the order of the number of voxels in the volume. A transformation is diffeomorphic if the transformation from the template to the target is one-to-one, onto, invertible, and both the transformation and it's inverse are differentiable. A transformation is said to be one-to-one if no two distinct points in the template are mapped to the same point in the target. A transformation is said to be onto if every point in the target is mapped from a point in the template. The importance of generating diffeomorphisms is that tangents, curvature, surface areas, and geometric properties of the imagery can be calculated automatically. FIG. 12 illustrates the image mapping illustrated in FIG. 11 using diffeomorphic transformation.

[0025] The present invention overcomes the limitations of the conventional techniques of image registration by providing a methodology which combines, or fuses, some aspects of techniques where an individual with expertise in the structure of the object represented in the images labels a set of landmarks in each image that are to be registered and techniques that use mathematics of small deformation multi-target registration, which is purely image data driven.

[0026] Additional features and advantages of the invention will be set forth in the description which follows, and in part, will be apparent from the description, or may be learned by practicing the invention. The embodiments and other advantages of the invention will be realized and obtained by the method and apparatus particularly pointed out in the written description and the claims hereof as well as in the appended drawings.

[0027] To achieve these and other advantages and in accordance with the purpose of the invention, as embodied and broadly described, an embodiment consistent with the present invention registers a template and a target image by defining a manifold landmark point in the template image, identifying a point in said target image corresponding to the defined manifold landmark point, and selecting a coordinate frame suitable for spherical geometries. According to an embodiment of the present invention, a large deformation transform is computed using the selected coordinate frame, a manifold landmark transformation operator, and at least one manifold landmark transformation boundary value. The large deformation transform relates the manifold landmark point in the template image to the corresponding point in the target image. The template and target image are registered using the large deformation transform in the selected coordinate frame. Another embodiment consistent with the present invention registers a template image with a target image by selecting a coordinate frame suitable for spherical geometries and registering the target and template image using a large deformation transform in the selected coordinate frame.

[0028] Both the foregoing general description and the following detailed description are exemplary and explanatory and are intended to provide further explanation of the invention as claimed.

[0029] The accompanying drawings provide a further understanding of the invention. They illustrate embodiments of the invention and, together with the description, explain the principles of the invention.

[0030] FIG. 1 is a target and template image of an axial section of a human head with 0-dimensional manifolds;

[0031] FIG. 2 is schematic diagram illustrating an apparatus for registering images in accordance with the present invention;

[0032] FIG. 3 is a flow diagram illustrating the method of image registration according to the present invention;

[0033] FIG. 4 is a target and a template image with 1-dimensional manifolds;

[0034] FIG. 5 is a target and a template image with 2-dimensional manifolds;

[0035] FIG. 6 is a target and a template image with 3-dimensional manifolds;

[0036] FIG. 7 is sequence of images illustrating registration of a template and target image; and

[0037] FIG. 8 is a flow diagram illustrating the computation of a fusing transform;

[0038] FIG. 9 is an oval template image which has landmark points selected and highlighted;

[0039] FIG. 10 is a deformed and distorted oval target image with corresponding landmark points highlighted and selected;

[0040] FIG. 11 is an image matching of the oval target and template images; and

[0041] FIG. 12 is an image matching using diffeomorphism.

[0042] 1. Method for Image Registration Using Both Landmark Based knowledge and Image Data

[0043] A method and system is disclosed which registers images using both landmark based knowledge and image data. Reference will now be made in detail to the present preferred embodiment of the invention, examples of which are illustrated in the accompanying drawings.

[0044] To illustrate the principles of this invention, FIG. 1 shows two axial views of a human head. In this example, template image 100 contains points 102, 104, and 114 identifying structural points (0-dimensional landmark manifolds) of interest in the template image. Target image 120 contains points 108, 110, 116, corresponding respectively to template image points 102, 104, 114, via vectors 106, 112, 118, respectively.

[0045] FIG. 2 shows apparatus to carry out the preferred embodiment of this invention. A medical imaging scanner 214 obtains the images show in FIG. 1 and stores them on a computer memory 206 which is connected to a computer central processing unit (CPU) 204. One of ordinary skill in the art will recognize that a parallel computer platform having multiple CPUs is also a suitable hardware platform for the present invention, including, but not limited to, massively parallel machines and workstations with multiple processors. Computer memory 206 can be directly connected to CPU 204, or this memory can be remotely connected through a communications network.

[0046] Registering images 100, 120 according to the present invention, unifies registration based on landmark deformations and image data transformation using a coarse-to-fine approach. In this approach, the highest dimensional transformation required during registration is computed from the solution of a sequence of lower dimensional problems driven by successive refinements. The method is based on information either provided by an operator, stored as defaults, or determined automatically about the various substructures of the template and the target, and varying degrees of knowledge about these substructures derived from anatomical imagery, acquired from modalities like CT, MRI, functional MRI, PET, ultrasound, SPECT, MEG, EEG, or cryosection.

[0047] Following this hierarchical approach, an operator, using pointing device 208, moves cursor 210 to select points 102, 104, 114 in FIG. 1, which are then displayed on a computer monitor 202 along with images 100, 120. Selected image points 102, 104, and 114 are 0-dimensional manifold landmarks.

[0048] Once the operator selects manifold landmark points 102, 104, and 114 in template image 100, the operator identifies the corresponding template image points 108, 110, 116.

[0049] Once manifold landmark selection is complete, CPU 204 computes a first transform relating the manifold landmark points in template image 100 to their corresponding image points in target image 120. Next, a second CPU 204 transform is computed by fusing the first transform relating selected manifold landmark points with a distance measure relating all image points in both template image 100 and target image 120. The operator can select an equation for the distance measure several ways including, but not limited to, selecting an equation from a list using pointing device 208, entering into CPU 204 an equation using keyboard 212, or reading a default equation from memory 206. Registration is completed by CPU 204 applying the second computed transform to all points in the template image 100.

[0050] Although several of the registration steps are described as selections made by an operator, implementation of the present invention is not limited to manual selection. For example, the transforms, boundary values, region of interest, and distance measure can be defaults read from memory or determined automatically.

[0051] FIG. 3 illustrates the method of this invention in operation. First an operator defines a set of N manifold landmark points x_i,where i=1, . . . , N, represented by the variable M, in the template image (step 300). These points should correspond to points that are easy to identify in the target image.

[0052] Associated with each landmark point, x_i, in the template image, is a corresponding point y_iin the target image. The operator therefore next identifies the corresponding points, y_i, in the target image are identified (step 310). The nature of this process means that the corresponding points can only be identified within some degree of accuracy. This mapping between the template and target points can be specified with a resolution having a Gaussian error of variance σ².

[0053] If a transformation operator has not been designated, the operator can choose a manifold landmark transformation operator, L, for this transformation computation. In this embodiment, the Laplacian $\left(\nabla =\frac{{\partial }^2}{\partial x_1^2}+\frac{{\partial }^2}{\partial x_2^2}+\frac{{\partial }_x}{\partial x_3^2}\right)$

[0054] is used for the operator L. Similarly, the operator can also select boundary values for the calculation corresponding to assumed boundary conditions, if these values have not been automatically determined or stored as default values. Here, infinite boundary conditions are assumed, producing the following equation for K, where K(x,x_i) is the Green's function of a volume landmark transformation operator L²(assuming L is self-adjoint): $\begin{array}{cc}K\left(x,x_i\right)=\left[\begin{array}{ccc}\Vert x-x_i\Vert & 0& 0\\ 0& \Vert x-x_i\Vert & 0\\ 0& 0& \Vert x-x_i\Vert \end{array}\right]& \left(8\right)\end{array}$

[0055] Using circulant boundary conditions instead of infinite boundary conditions provides and embodiment suitable for rapid computation. One of ordinary skill in the art will recognize that other operators can be used in place of the Laplacian operator, such operators include, but are not limited to, the biharmonic operator, linear elasticity operator, and other powers of these operators.

[0056] In addition, the operator may select a region of interest in the target image. Restricting the computation to a relatively small region of interest reduces both computation and storage requirements because transformation is computed only over a subregion of interest. It is also possible that in some applications the entire image is the desired region of interest. In other applications, there may be default regions of interest that are automatically identified.

[0057] The number of computations required is proportional to the number of points in the region of interest, so the computational savings equals the ratio of the total number of points in the image to the number of points in the region of interest. The data storage savings for an image with N points with a region of interest having M points is a factor of N/M. For example, for a volume image of 256×256×256 points with a region of interest of 128×128×128 points, the computation time and the data storage are reduced by a factor of eight.

[0058] In addition, performing the computation only over the region of interest makes it necessary only to store a subregion, providing a data storage savings for the template image, the target image, and the transform values.

[0059] Following the identification of template manifold landmark points and corresponding points in the target image, as well as selection of the manifold transformation operator, the boundary values, and the region of interest, CPU 204 computes a transform that embodies the mapping relationship between these two sets of points (step 350). This transform can be estimated using Bayesian optimization, using the following equation: $\begin{array}{cc}\hat{u}=\arg \underset{u}{\min }{\int }_{\Omega }{\left|\mathrm{Lu}\right|}^2+\sum _{i=1}^N\frac{{\left|y_i-x_i+u\left(x_i\right)\right|}^2}{{\sigma }_i^2},& \left(9\right)\end{array}$

[0060] the minimizer, u, having the form $\begin{array}{cc}\hat{u}\left(x\right)=b+Ax+\sum _{i=1}^N{\beta }_{i}K\left(x,x_i\right)& \left(10\right)\end{array}$

[0061] where A is a 3×3 matrix, b=[b₁, b₂, b₃] is a 3×1 vector, [β_i1,β_i2,β_i3] is a 3×1 weighting vector.

[0062] The foregoing steps of the image registration method provide a coarse matching of a template and a target image. Fine matching of the images requires using the full image data and the landmark information and involves selecting a distance measure by solving a synthesis equation that simultaneously maps selected image landmarks in the template and target images and matches all image points within a region of interest. An example of this synthesis equation is: $\begin{array}{cc}\hat{u}=\arg \underset{u}{\min }\gamma {\int }_{\Omega }{\left|T\left(x-u\left(x\right)\right)-S\left(x\right)\right|}^{2}dx+{\int }_{\Omega }{\left|Lu\right|}^2+\sum _{i=1}^N\frac{{\left|y_i-x_i+u\left(x_i\right)\right|}^2}{{\sigma }_i^2}& \left(11\right)\end{array}$

[0063] here the displacement field u is constrained to have the form $\begin{array}{cc}u\left(x\right)=\sum _{k=0}^d{\mu }_k{\varphi }_k\left(x\right)+\sum _{i=1}^N{\beta }_{i}K\left(x,x_i\right)+Ax+b& \left(12\right)\end{array}$

[0064] with the variables β_i, A, and b, computed at step 350 in FIG. 3. The operator L in equation (11) may be the same operator used in equation (9), or alternatively, another operator may be used with a different set of boundary conditions. The basis functions φ are the eigen functions of operators such as the Laplacian Lu=∇²u, the bi-harmonic Lu=∇⁴u, linear elasticity Lu=α∇²u+(α+β)∇(∇.u), and powers of these operators L^pfor p≧1.

[0065] One of ordinary skill in the art will recognize that there are many possible forms of the synthesis equation. For example, in the synthesis equation presented above, the distance measure in the first term measures the relative position of points in the target image with respect to points in the template image. Although this synthesis equation uses a quadratic distance measure, one of ordinary skill in the art will recognize that there are other suitable distance measures.

[0066] CPU 204 then computes a second or fusing transformation (Step 370) using the synthesis equation relating all points within a region of interest in the target image to all corresponding points in the template image. The synthesis equation is defined so that the resulting transform incorporates, or fuses, the mapping of manifold landmarks to corresponding target image points determined when calculating the first transform.

[0067] The computation using the synthesis equation is accomplished by solving a sequence of optimization problems from coarse to fine scale via estimation of the basis coefficients μ_k. This is analogous to multi-grid methods, but here the notion of refinement from coarse to fine is accomplished by increasing the number of basis components d. As the number of basis functions increases, smaller and smaller variabilities between the template and target are accommodated. The basis coefficients are determined by gradient descent, i.e., $\begin{array}{cc}{\mu }_k^{\left(n+1\right)}={\mu }_k^{\left(n\right)}-\Delta \frac{\partial H\left(u^{\left(n\right)}|S\right)}{\partial {\mu }_k}\mathrm{where}& \left(13\right)\\ \frac{\partial H\left(u^{\left(n\right)}\right)}{\partial {\mu }_k}=-\gamma {\int }_{\Omega }\left(T\left(x-u^{\left(n\right)}\left(x\right)\right)-S\left(x\right)\right)\nabla T\left(x-u^{\left(n\right)}\left(x\right)\right)·{\varphi }_k\left(x\right)dx+{\lambda }_k^2{\mu }_k^{\left(n\right)}+\left(2\sum _{i=1}^N\frac{y_i-x_i+u^{\left(n\right)}\left(x_i\right)}{{\sigma }_i^2}\right)·{\varphi }_k\left(x_i\right)\mathrm{and}& \left(14\right)\\ u^{\left(n\right)}\left(x\right)=\sum _{k=0}^d{\mu }_k^{\left(n\right)}{\varphi }_k\left(x\right)+\sum _{i=1}^N{\beta }_{i}K\left(x,x_i\right)+Ax+b& \left(15\right)\end{array}$

[0068] also ∇ is a fixed step size and λ_kare the eigenvalues of the eigenvectors φ_k.

[0069] The computation of the fusion transformation (step 370) using the synthesis equation is presented in the flow chart of FIG. 8. Equation (12) is used to initialize the value of the displacement field u(x)=u⁽⁰⁾(x) (step 800). The basis coefficients μ_k=μ_k⁽⁰⁾are set equal to zero and the variables β_i, A, and b are set equal to the solution of equation (11) (step 802). Equation (13) is then used to estimate the new values of the basis coefficients μ_k⁽ⁿ⁺¹⁾given the current estimate of the displacement field u⁽ⁿ⁾(x) (step 804). Equation (15) is then used to compute the new estimate of the displacement field u⁽ⁿ⁾(x) given the current estimate of the basis coefficients μ_k⁽ⁿ⁾(step 806). The next part of the computation is to decide whether or not to increase the number d of basis functions φ_kused to represent the transformation (step 808). Increasing the number of basis functions allows more deformation. Normally, the algorithm is started with a small number of basis functions corresponding to low frequency eigen functions and then on defined iterations the number of frequencies is increased by one (step 810). This coarse-to-fine strategy matches larger structures before smaller structures. The preceding computations (steps 804-810) are repeated until the computation has converged or the maximum number of iterations is reached (step 812). The final displacement field is then used to transform the template image (step 814).

[0070] Once CPU 204 determines the transform from the synthesis equation fusing both landmark manifold information and image data, CPU 204 uses this transform to register the template image with the target image (step 380).

[0071] The spectrum of the second transformation, h, is highly concentrated around zero. This means that the spectrum mostly contains low frequency components. Using the sampling theorem, the transformation can be represented by a subsampled version provided that the sampling frequency is greater than the Nyquist frequency of the transformation. The computation may be accelerated by computing the transformation on a coarse grid and extending it to the full voxel lattice e.g., in the case of 3D images, by interpolation. The computational complexity of the algorithm is proportional to the dimension of the lattice on which the transformation is computed. Therefore, the computation acceleration equals the ratio of the full voxel lattice to the coarse computational lattice.

[0072] One of ordinary skill in the art will recognize that composing the landmark transformations followed by the elastic basis transformations, thereby applying the fusion methodology in sequence can provide an alternatively valid approach to hierarchial synthesis of landmark and image information in the segmentation.

[0073] Another way to increase the efficiency of the algorithm is to precompute the Green's functions and eigen functions of the operator L and store these precomputed values in a lookup table. These tables replace the computation of these functions at each iteration with a table lookup. This approach exploits the symmetry of Green's functions and eigen functions of the operator L so that very little computer memory is required. In the case of the Green's functions, the radial symmetry is exploited by precomputing the Green's function only along a radial direction.

[0074] The method described for fusing landmark information with the image data transformation can be extended from landmarks that are individual points (0-dimensional manifolds) to manifolds of dimensions 1, 2 and 3 corresponding to curves (1-dimensional), surfaces (2-dimensional) and subvolumes (3-dimensional).

[0075] For example, FIG. 4 shows a template image 400 of a section of a brain with 1-dimensional manifolds 402 and 404 corresponding to target image 406 1-dimensional manifolds 408 and 410 respectively. FIG. 5 shows a template image 500 of a section of a brain with 2-dimensional manifold 502 corresponding to target image 504 2-dimensional manifold 506. FIG. 6 shows a template image 600 of a section of a brain with 3-dimensional manifold 602 corresponding to target image 604 3-dimensional manifold 606.

[0076] As with the point landmarks, these higher dimensional manifolds condition the transformation, that is we assume that the vector field mapping the manifolds in the template to the data is given. Under this assumption the manually-assisted deformation (step 350, FIG. 3) becomes the equality-constrained Bayesian optimization problem: $\begin{array}{cc}u\left(x\right)=\arg \underset{u}{\min }{\int }_{\Omega }{\left|Lu\left(x\right)\right|}^{2}dx\mathrm{subject}\mathrm{to}& \left(16\right)\\ u\left(x\right)=k\left(x\right),x\in \bigcup _{i=0}^{3}M\left(i\right).& \left(17\right)\end{array}$

[0077] If M(i) is a smooth manifold for i=0, 1, 2, 3, the solution to this minimization is unique satisfying L^†Lû(x)=0, for all template points in the selected manifold. This implies that the solution can be written in the form of a Fredholm integral equation: $\begin{array}{cc}\hat{u}\left(x\right)={\int }_{\bigcup _{i=0}^{3}M\left(i\right)}K\left(x,y\right)\beta \left(y\right)dS\left(y\right),\mathrm{where}K=G{G}^{†}& \left(18\right)\end{array}$

[0078] and G the Green's function of L.

[0079] When the manifold is a sub-volume, M(3),dS is the Lebesgue measure on R³. For 2-dimensional surfaces, dS is the surface measure on M(2), For 1-dimensional manifolds (curves), dS is the line measure on M(1) and for point landmarks, M(0), dS is the atomic measure. For point landmarks, the Fredholm integral equation degenerates into a summation given by equation (10).

[0080] When the manifold of interest is a smooth, 2-dimensional surface, the solution satisfies the classical Dirichlet boundary value problem:

[0081] L^†Lû(x)=0,∀xεΩ\M  (19)

[0082] The Dirichlet problem is solved using the method of successive over relaxation as follows. If u^k(x) is the estimate of a deformation field at the k^thiteration, the estimate at the (k+1)^thiteration is given by the following update equation:

[0083] u^k+1(x)=u^k(x)+αL,^†Lu(x),xεΩ\M.

[0084] u^k+1(x)=k(x),xεM,  (20)

[0085] where α is the over relaxation factor.

[0086] It is also possible to compute the transform (step 370) with rapid convergence by solving a series of linear minimization problems where the solution to the series of linear problems converges to the solution of the nonlinear problem. This avoids needing to solve the nonlinear minimization problem directly. Using a conjugate gradient method, the computation converges faster than a direct solution of the synthesis equation because the basis coefficients μ_kare updated with optimal step sizes.

[0087] Using the conjugate gradient, the displacement field is assumed to have the form $\begin{array}{cc}u\left(x\right)=\sum _{k=0}^d{\mu }_k{\varphi }_k\left(x\right)+f\left(x\right)\mathrm{where}& \left(21\right)\\ f\left(x\right)=\sum _{i=1}^N{\beta }_{i}K\left(x,x_i\right)+Ax+b.& \left(22\right)\end{array}$

[0088] Begin by assuming that f(x) is fixed. This is generalized below. The eigen functions in the expansion are all real and follow the assumption that {φ_i(x)} are R³valued.

[0089] The minimization problem is solved by computing

[0090] μ_j^new=μ_j^old+Δ_jj=0 . . . d   (23)

[0091] to update the basis coefficients in equation (21) where μ_j=0,j=0 . . . d initially, Δ_jis computed using the equation $\begin{array}{cc}\Delta j={\left({\int }_{\Omega }\left(\sum _{k=0}^d{h}_k\left(x\right)\right)h_j\left(x\right)dx+{\lambda }_j^2+\frac{1}{{\sigma }^2}\sum _{i=1}^N{\theta }_{\mathrm{jj}}\left(x_i\right)\right)}^{-1}\left({\int }_{\Omega }\left(T\left(x-u\left(x\right)\right)-S\left(x\right)\right)h_j\left(x\right)dx+\sum _{i=1}^N{\beta }_i-{\varphi }_j\left(x_i\right)+\frac{1}{{\sigma }^2}\sum _{i=1}^N\left(y_i-x_i+u\left(x_i\right)\right)·{\varphi }_j\left(x_i\right)+\frac{1}{{\sigma }^2}\sum _{k=0}^{j-1}\Delta k\sum _{i=N}^N{\theta }_{\mathrm{kj}}\left(x_i\right)\right).& \left(24\right)\end{array}$

[0092] where h_i(x)=∇T|_x−u(x)·φ_i(x), and where θ_ij(x)=φ_j(x)·φ_j(x). The notation f·g is the inner-product, i.e., $f·g=\sum _{i=1}^3{f}_i{g}_i\mathrm{for}f,g\in R^3.$

[0093] Similarly, since u(x) is written in the series expansion given in equations (21) and (22), the identical formulation for updating β_iarises. Accordingly, the fusion achieved by the present invention results. Computation of equation (23) repeats until all Δ_jfall below a predetermined threshold solving for each Δ_jin sequence of increasing j, and Δ_jis computed using the values of Δ_kfor 0≦k<j.

[0094] A further improvement over prior art image registration methods is achieved by computing the required transforms using fast Fourier transforms (FFT). Implementing an FFT based computation for image registration using a synthesis equation, as required at step 370 of FIG. 3, provides computational efficiency. However, to exploit the known computational efficiencies of FFT transforms, the solution of the synthesis equation must be recast to transform the inner products required by iterative algorithms into shift invariant convolutions.

[0095] To make the inner-products required by the iterative algorithms into shift invariant convolution, differential and difference operators are defined on a periodic version of the unit cube and the discrete lattice cube. Thus, the operators are made cyclo-stationary, implying their eigen functions are always of the form of complex exponentials on these cubes having the value: $\begin{array}{cc}{\Psi }_k^{\left(r\right)}=\left[\begin{array}{c}{c}_{1k}^{\left(r\right)}\\ c_{2k}^{\left(r\right)}\\ c_{3k}^{\left(r\right)}\end{array}\right]e^{j⟨{\omega }_k,x⟩}.& \left(25\right)\end{array}$

[0096] r=1, 2, 3 with x=(x₁,x₂,x₃)ε[0, 1]³,

[0097] ω_ki=2πk_i, i=1, 2, 3, and the Fourier basis for periodic functions on [0, 1]³takes the form $\begin{array}{cc}{e}^{j⟨{\omega }_k,x⟩},{\omega }_{k}x={\omega }_{k_1}{x}_1+{\omega }_{k_2}{x}_2+{\omega }_{k_3}{x}_3.{\omega }_k=\left({\omega }_{k_1},{\omega }_{k_2},{\omega }_{k_3}\right)\mathrm{On}\mathrm{the}\mathrm{discrete}{N}^3=\mathrm{periodic}\mathrm{lattice},{\omega }_k=\left(\frac{2\pi k_1}{N},\frac{2\pi k_2}{N},\frac{2\pi k_3}{N}\right),x\in {\left\{0,1\dots N-1\right\}}^3.\mathrm{For}\mathrm{real}\mathrm{expansions},\mathrm{the}\mathrm{eigen}\mathrm{vectors}\mathrm{becomes}{\varphi }_k\left(x\right)={\Psi }_k\left(x\right)+{\Psi }_k^{*}\left(x\right)\mathrm{and}\mathrm{the}\mathrm{real}\mathrm{expansion}\mathrm{in}\mathrm{equation}\left(21\right)\mathrm{becomes}\text{:}& \\ u\left(x\right)=\sum _{k=0}^d{\mu }_k\left({\Psi }_k\left(x\right)+{\Psi }_k^{*}\left(x\right)\right){\mathrm{where}}^{*}\mathrm{means}\mathrm{complex}\mathrm{conjugate},\mathrm{and}0<d\le \frac{N}{2}.& \left(26\right)\end{array}$

[0098] This reformulation supports an efficient implementation of the image registration process using the FFT. Specifically, if step 370 of FIG. 3, computing the registration transform fusing landmark and image data, is implemented using the conjugate gradient method, the computation will involve a series of inner products. Using the FFT exploits the structure of the eigen functions and the computational efficiency of the FFT to compute these inner-products.

[0099] For example, one form of a synthesis equation for executing Step 370 of FIG. 3 will include the following three terms: $\text{Term 1:}$ ${\int }_{\Omega }\left(T\left(x-u\left(x\right)\right)-S\left(x\right)\right)h_j\left(x\right)dx$ $\text{Term 2:}$ ${\int }_{\Omega }\left(\sum _{k=0}^d{h}_k\left(x\right)\right)h_j\left(x\right)dx$ $\text{Term 3:}$ $u\left(x\right)=\sum _{k=0}^d{\mu }_k{\varphi }_k\left(x\right)$

[0100] Each of theses terms must be recast in a suitable form for FFT computation. One example of a proper reformulation for each of these terms is:

[0101] Term 1: $\begin{array}{cc}{\int }_{\Omega }\left(T\left(x-u\left(x\right)\right)-S\left(x\right)\right)\nabla T·\left({\Psi }_i^r\left(x\right)+{\Psi }_i^{\left(r\right)*}\left(x\right)\right)dx=2\mathrm{Re}\left({\int }_{\Omega }\left(T\left(x-u\left(x\right)\right)-S\left(x\right)\right)\left(\sum _{r=1}^3\nabla T·c_k^{\left(r\right)}\right)e^{j⟨{\omega }_k,x⟩}dx\right),& \left(27\right)\end{array}$

[0102] where c_k^(r)=[c_1k^(r), c_2k^(r), c_3k^(r)]^t. This equation is computed for all k by a Fourier transformation of the function. $\begin{array}{cc}\left(T\left(x-u\left(x\right)\right)-S\left(x\right)\right)\left(\sum _{r=1}^3\nabla T·c_k^{\left(r\right)}\right)& \left(28\right)\end{array}$

[0103] and hence can be computed efficiently using a 3-D FFT.

[0104] Term 2: $\begin{array}{cc}\sum _{r=1}^3\sum _{k=0}^d\int {\left({\Psi }_j^{\left(r\right)}+{\Psi }_j^{\left(r\right)*}\right)}^t\left(\nabla {T\left(\nabla T\right)}^t\right)\left({\Psi }_j^{\left(r\right)}+{\Psi }_j^{\left(r\right)*}\right)dx=2\mathrm{Re}\left(\sum _{r=1}^3\sum _{k=0}^d\int {\left(c_k^{\left(r\right)}\right)}^t\left({\int }_{\Omega }\left(\nabla {T\left(\nabla T\right)}^t\right)e^{⟨{\omega }_k·{\omega }_j,x⟩}dx\right)\right)& \left(29\right)\end{array}$

[0105] The integral in the above summation for all k can be computed by Fourier transforming the elements of the 3×3 matrix:

[0106] ∇T(∇T)^t  (30)

[0107] evaluated at ω_k+ω_j.Because this matrix has diagonal symmetry, the nine FFTs in this reformulation of term 2 can be computed efficiently using six three dimensional FFTs evaluated at ω_k+ω_j.

[0108] Term 3:

[0109] Using the exact form for the eigen functions we can rewrite the above equation as $\begin{array}{cc}u\left(x\right)=2\mathrm{Re}\left(\sum _{r=1}^3\sum _{k=0}^d{\mu }_k^{\left(r\right)}\left[\begin{array}{c}{c}_{1k}^{\left(r\right)}\\ c_{2k}^{\left(r\right)}\\ c_{3k}^{\left(r\right)}\end{array}\right]e^{j⟨{\omega }_k,x⟩}\right).& \left(31\right)\end{array}$

[0110] This summation is precisely the inverse Fourier transforms of the functions $\sum _{r=1}^3{\mu }_k^{\left(r\right)}{c}_{\mathrm{ik}}^{\left(r\right)}\mathrm{for}i=1,2,3$

[0111] and hence can be computed efficiently by using a 3-D FFT.

[0112] One of ordinary skill in the art will recognize that restructuring the computation of registration transforms using FFTs will improve the performance of any image registration method having terms similar to those resulting from a synthesis equation fusing landmark and image data. Improvements results from the fact that many computer platforms compute FFTs efficiently; accordingly, reformulating the registration process as an FFT computation, makes the required computations feasible.

[0113] A distance function used to measure the disparity between images is the Gaussian square error distance ∫|T(x−u(x))−S(x)|²dx. There are many other forms of an appropriate distance measure. More generally, distance functions, such as the correlation distance, or the Kullback Liebler distance, can be written in the form

[0114] ∫D(T(x−u(x)), S(x))dx.

[0115] An efficient convolution implementation can be derived using the FFT for arbitrary distance functions. Computing the fusing transform using the image data follows the equation: $\begin{array}{cc}\hat{u}=\arg \underset{u}{\min }\gamma {\int }_{\Omega }D\left(T\left(x-u\left(x\right)\right),S\left(x\right)\right)dx+{\int }_{\Omega }{\left|\mathrm{Lu}\right|}^2+\sum _{i=1}^N\frac{{\left|y_i-x_i+u\left(x_i\right)\right|}^2}{{\sigma }_i^2}& \left(32\right)\end{array}$

[0116] where D(.,.) is a distance function relating points in the template and target images. The displacement field is assumed to have the form: $\begin{array}{cc}u\left(x\right)=\sum _{k=0}^d{\mu }_k{\varphi }_k\left(x\right)+f\left(x\right)\mathrm{where}& \left(33\right)\\ f\left(x\right)=b+Ax+\sum _{i=1}^N{\beta }_{i}K\left(x,x_i\right)& \left(34\right)\end{array}$

[0117] is fixed. The basis coefficients {μ_k} are determined by gradient descent, i.e., $\begin{array}{cc}{\mu }_k^{\left(n+1\right)}={\mu }_k^{\left(n\right)}-\Delta \frac{\partial H\left(u^{\left(n\right)}|S\right)}{\partial {\mu }_k}& \left(35\right)\end{array}$

[0118] where the gradient is computed using the chain rule and is given by the equation $\begin{array}{cc}\frac{\partial H\left(u^{\left(n\right)}\right)}{\partial {\mu }_k}={\int }_{\Omega }{D}^{\prime }\left(T\left(x-u^{\left(n\right)}\left(x\right)\right),S\left(x\right)\right)\nabla T\left(x-u^{\left(n\right)}\left(x\right)\right)·{\varphi }_k\left(x\right)dx+{\lambda }_k^2{\mu }_k^{\left(n\right)}+\left(2\sum _{i=1}^N\frac{y_i-x_i+u^{\left(n\right)}\left(x_i\right)}{{\sigma }_i^2}\right)·{\varphi }_k\left(x\right)& \left(36\right)\end{array}$

[0119] where D′ (.,.) is the derivative with respect to the first argument. The most computationally intensive aspect of the algorithm is the computation of the term ${\int }_{\Omega }{D}^{\prime }\left(T\left(x-u^{\left(n\right)}\left(x\right)\right),S\left(x\right)\right)\nabla T\left(x-u^{\left(n\right)}\left(x\right)\right)·{\varphi }_k\left(x\right)dx$

[0120] Using the structure of the eigen functions and the computational efficiency of the FFT to compute these inner-products, the above term can be written as $2\mathrm{Re}{\int }_{\Omega }{D}^{\prime }\left(T\left(x-u^{\left(n\right)}\left(x\right)\right),S\left(x\right)\right)\left(\sum _{r=1}^3\nabla T·c_k^{\left(r\right)}\right)e^{j⟨{\omega }_k,x⟩}dx$

[0121] where

[0122] c_k^(r)=[c_1k^(r), c_2k^(r), c_3k^(r)]^t. This equation is computed for all k by a Fourier transformation of the function $D^{\prime }\left(T\left(x-u^{\left(n\right)}\left(x\right)\right),S\left(x\right)\right)\left(\sum _{r=1}^3\nabla T·c_k^{\left(r\right)}\right)$

[0123] and hence can be computed efficiently using a 3-D FFT.

[0124] The following example illustrates the computational efficiencies achieved using FFTs for image registration instead of direct computation of inner-products. Assuming that a target image is discretized on a lattice having N³points, each of the inner-products in the algorithm, if computed directly, would have a computational complexity of the order (N³)². Because the inner-products are computationally intensive, the overall complexity of image registration is also (N³)². In contrast, each of the FFTs proposed has a computational complexity on the order of N³log₂N³. The speed up is given by the ratio N⁶/(N³log₂N³)=N³/(3 log₂N). Thus the speed up is 64 times for a 16×16×16 volume and greater than 3.2×10⁴speed up for a 256×256×256 volume.

[0125] A further factor of two savings in computation time can be gained by exploiting the fact that all of the FFTs are real. Hence all of the FFTs can be computed with corresponding complex FFTs of half the number of points. For a development of the mathematics of FFTs see, A. V. Oppenheim and R. W. Schafer, Digital Signal Processing,Prentice-Hall, New Jersey, 1975 (hereinafter referred to as Oppenheim).

[0126] Alternative embodiments of the registration method described can be achieved by changing the boundary conditions of the operator. In the disclosed embodiment, the minimization problem is formulated with cyclic boundary conditions. One of ordinary skill in the art will recognize that alternative boundary conditions, such as the Dirichlet, Neumann, or mixed Dirichlet and Neumann boundary conditions are also suitable. The following equation is used in an embodiment of the present invention using one set of mixed Dirichlet and Neumann boundary conditions: $\begin{array}{cc}{\left(\frac{\partial u_t}{\partial x_t}|\right)}_{\left(x|x_i=k\right)}=u_i\left(x|x_j=k\right)=0\mathrm{for}i,j=1,2,3;i\ne j;k=0,1\text{:}& \left(37\right)\end{array}$

[0127] where the notation (x|x_i=k) means x is in the template image such that x_i=k. In this case, the eigen functions would be of the form: $\begin{array}{cc}{\varphi }_k^{\left(r\right)}=\left[\begin{array}{c}{C}_{1k}^{\left(r\right)}\cos {\omega }_{\mathrm{k1}}{x}_1\sin {\omega }_{\mathrm{k2}}\sin {\omega }_{\mathrm{k3}}{x}_3\\ C_{2k}^{\left(r\right)}\sin {\omega }_{\mathrm{k1}}{x}_1\cos {\omega }_{\mathrm{k2}}\sin {\omega }_{\mathrm{k3}}{x}_3\\ C_{3k}^{\left(r\right)}\sin {\omega }_{\mathrm{k1}}{x}_1\sin {\omega }_{\mathrm{k2}}\cos {\omega }_{\mathrm{k3}}{x}_3\end{array}\right]\mathrm{for}r=1,2,3.& \left(38\right)\end{array}$

[0128] Modifying boundary conditions requires modifying the butterflies of the FFT from complex exponentials to appropriate sines and cosines.

[0129] In FIG. 7, four images, template image 700, image 704, image 706, and target image 708, illustrate the sequence of registering a template image and a target image. Template image 700 has 0-dimensional landmark manifolds 702. Applying the landmark manifold transform computed at step 350 in FIG. 3 to image 700 produces image 704. Applying a second transform computed using the synthesis equation combining landmark manifolds and image data to image 700 produces image 706. Image 706 is the final result of registering template image 700 with target image 708. Landmark manifold 710 in image 708 corresponds to landmark manifold 702 in template image 700.

[0130] Turning now to techniques to register images which may possess large deformation characteristics using diffeomorphisms, large deformation transform functions and transformations capable of matching images where the changes from one image to the other are greater than small, linear, or affine. Use of large deformation transforms provide for the automatic calculation of tangents, curvature, surface areas, and geometric properties of the imagery.

[0131] 2. Methods for Large Deformation Landmark Based and Image Based Transformations

[0132] An embodiment consistent with the present invention maps sets of landmarks in imagery {x_i, i=1,2, . . . , N} ⊂ Ω into target landmarks {y_i, i=1, . . . , N}, and or imagery I₀into target I₁both with and without landmarks. For example when there is a well-defined distance function D(u(T)) expressing the distance between the landmarks and or imagery. The large deformation maps h: Ω→Ω are constructed by introducing the time variable,

[0133] h(x,t)=(x₁,x₂,x₃,t)εΩ×[0,T]→h(x,t)=(x−u₁(x,t),x₂−u₂(x,t),x₃−u₃(x,t))εΩ.

[0134] The large deformation maps are constrained to be the solution h(x, T)=x−u(x,T) where u(x, T) is generated as the solution of the ordinary differential equation $\begin{array}{cc}u\left(x,T\right)={\int }_o^T\left(I-\nabla u\left(x,t\right)\right)v\left(x,t\right)dt,\mathrm{where}v\left(x,t\right)=\sum _{k-0}^{\infty }{v}_k{\varphi }_k\left(x,t\right),\left(x,t\right)\in \Omega \times \left[0,T\right]& \left(39\right)\end{array}$

[0135] assuming that the {φ_k} forms a complete orthonormal base.

[0136] A differomorphism for the landmark and image matching problem is given by the mapping of the imagery given by ĥ(x,T)=x−û(x,T) satisfying the Ordinary Differential Equation (O.D.E.) $\begin{array}{cc}\hat{u}\left(x,T\right)={\int }_o^T\left(I-\nabla \hat{u}\left(x,t\right)\right)\hat{v}\left(x,t\right)dt& \left(40\right)\end{array}$ $\begin{array}{cc}\mathrm{where}\hat{v}\left(x,t\right)={\arg }_{\mathrm{vi}},\underset{k=1,\dots }{\min }{\int }_o^T{\int }_{\Omega }{\Vert Lv\left(x,t\right)\Vert }^{2}dxdt+D\left(u\left(T\right)\right).& \left(41\right)\end{array}$

[0137] L is preferably a linear differential operator with the φ_kforming a complete orthonormal base as the eigenfunctions Lφ_k=λ_kφ_k.

[0138] 2.1 Large Deformation Landmark Matching

[0139] In order to register images with large deformation using a landmark matching technique, N landmarks are identified in the two anatomies targert {x_i, y_i, i=1,2, . . . , N}. The landmarks are identified with varying degrees of accuracy {Σ_i, i=1, . . . , N}, the Σ_i, 3×3 covariance matrices. The distance between the target and template imagery landmarks is preferably defined as

[0140] $\begin{array}{cc}{D}_1\left(u\left(T\right)\right)\stackrel{.}{=}\sum _{i=1}^N\left(y_i-x_i-u\left(x_i,T\right)\right){\Sigma }_i^{-1}\left(y_i-x_i-u\left(x_i,T\right)\right).& \left(42\right)\end{array}$

[0141] A preferable differomorphism is the minimizer of Eqns. 40, 41 with D(u(T)) the landmark distance: $\begin{array}{cc}\hat{u}\left(x,T\right)={\int }_0^T\left(I-\nabla u\left(x,t\right)\right)\hat{v}\left(x,t\right)dt& \left(43\right)\end{array}$ $\begin{array}{cc}\mathrm{where}\hat{v}\left(x,t\right)=\arg \underset{v_k,k=1,\dots }{\min }\left(H\left(u\right)\left(T\right)\right)\stackrel{.}{=}{\int }_o^T{\int }_{\Omega }{\Vert Lv\left(x,t\right)\Vert }^{2}dxdt+D_1\left(u\left(T\right)\right)).& \left(44\right)\end{array}$

[0142] A method for registering images consistent with the present invention preferably includes the following steps:

[0143] STEP 0: Define a set of landmarks in the template which can be easily identified in the target {x_i:x_iεΩ, i=1,2, . . . , N} with varying degrees of accuracy in the target as {yi, i=1, . . . , N} with associated error covariance {Σ_i, i=1, . . . , N}, and initialize for n=0, ν.=0.

[0144] STEP 1: Calculate velocity and deformation fields: $\begin{array}{cc}{v}^{\left(n\right)}\left(x,t\right)=\sum _{k=0}^{\infty }{v}_k^{\left(n\right)}{\varphi }_k\left(x,t\right),u^{\left(n\right)}\left(x,T\right)={\int }_o^T\left(I-\nabla u^{\left(n\right)}\left(x,t\right)\right)v^{\left(n\right)}\left(x,t\right)dt.& \left(45\right)\end{array}$

[0145] STEP 2: Solve via a sequence of optimization problems from coarse to fine scale via estimation of the basis coefficients {ν_k}, analogous to multi-grid methods with the notion of refinement from coarse to fine accomplished by increasing the number of basis components. For each ν_k, $\begin{array}{cc}{v}_k^{\left(n+1\right)}=v_k^{\left(n\right)}-\Delta \frac{\partial H\left(v^{\left(n\right)}\right)}{\partial v_k}=v_k^{\left(n\right)}-\Delta \left({\lambda }_k^2{v}_k^{\left(n\right)}=\frac{\partial D\left(u^{\left(n\right)}\left(T\right)\right)}{\partial v_k}\right)\mathrm{where}\begin{array}{c}\left(1\right)\frac{\partial D\left(u\left(T\right)\right)}{\partial v_k}=-2\sum _{i=1}^N⟨y_i-x_i-u^{\left(n\right)}\left(x,T\right),\frac{\partial u^{\left(n\right)}\left(x,T\right)}{\partial {\gamma }_k}⟩\\ \begin{array}{c}\left(2\right)\frac{\partial u^{\left(n\right)}\left(x,T\right)}{\partial v_k}={\int }_o^T\left(I-\nabla u^{\left(n\right)}\left(x,t\right)\right){\varphi }_k\left(x,t\right)dt-\\ {\int }_o^T\nabla \frac{\partial u^{\left(n\right)}\left(x,t\right)}{\partial v_k}{v}^{\left(n\right)}\left(x,t\right)dt.\end{array}\end{array}& \left(46\right)\end{array}$

[0146] STEP 3: Set n←n+1, and return to Step 1.

[0147] 2.2 Large Deformation Image Matching

[0148] Turning now to the technique where a target and template image with large deformation are registered using the image matching technique. The distance between the target and template image is chosen in accordance with Eqn. 47. A diffeomorphic map is computed using an image transformation operator and image transformation boundary values relating the template image to the target image. Subsequently the template image is registered with the target image using the diffromorphic map.

[0149] Given two images I₀, I₁, choose the distance $\begin{array}{cc}{D}_2\left(u\left(T\right)\right)\stackrel{.}{=}\gamma {\int }_{{\left[0,1\right]}^3}{\left|I_0\left(x-u\left(x,T\right)\right)-I_1\left(x\right)\right|}^{2}dx.& \left(47\right)\end{array}$

[0150] The large deformation image distance driven map is constrained to be the solution ĥ(x,T)=x−û(x,T) where $\begin{array}{cc}\hat{u}\left(x,T\right)={\int }_o^T\left(I-\nabla \hat{u}\left(x,t\right)\right)\hat{v}\left(x,t\right)dt\mathrm{and}& \left(48\right)\end{array}$

[0151] A method for registering images consistent with the present invention preferably includes the following steps: $\begin{array}{cc}\hat{v}\left(x,t\right)=\arg \underset{\left\{V\text{:}v\left(x,t\right)=\sum _{k}vk\varphi k\right\}}{\min }\left(H\left(v\right)\stackrel{.}{=}{\int }_o^T{\int }_{\Omega }{\Vert Lv\left(x,t\right)\Vert }^{2}dxdt+D_2\left(u\left(T\right)\right)\right).& \left(49\right)\end{array}$

[0152] STEP 0: Measure two images I₀, I₁defined on the unit cube Ω≈[0,1]³⊂R⁸, initialize parameters for n=0, ν_k⁽ⁿ⁾=0, k=0,1, . . . , and define distance measure D(u(T)) Eqn. 47.

[0153] STEP 1: Calculate velocity and deformation fields from ν⁽ⁿ⁾: $\begin{array}{cc}{v}^{\left(n\right)}\left(x,t\right)=\sum _{k=0}^{\infty }{v}_k^{\left(n\right)}{\varphi }_k\left(x,t\right),u^{\left(n\right)}\left(x,T\right)={\int }_0^T\left(I-\nabla u^{\left(n\right)}\left(x,t\right)\right)v^{\left(n\right)}\left(x,t\right)dt.& \left(50\right)\end{array}$

[0154] STEP 2: Solve optimization via sequence of optimization problems from coarse to fine scale via re-estimation of the basis coefficients {ν_k}, analogous to multi-grid methods with the notion of refinement from coarse to fine accomplished by increasing the number of basis components. For each ν_k, $\begin{array}{cc}{v}_k^{\left(n+1\right)}\stackrel{.}{=}{v}_k^{\left(n\right)}-\Delta \frac{\partial H\left(v^{\left(n\right)}\right)}{\partial v_k}=v_k^{\left(n\right)}-\Delta \left({\lambda }_k^2{v}_k^{\left(n\right)}+\frac{\partial D\left(u^{\left(n\right)}\left(T\right)\right)}{\partial v_k}\right)\mathrm{where}& \left(51\right)\\ {\left(\frac{\partial D\left(u\left(T\right)\right)}{\partial v_k}|\right)}_{u^{\left(n\right)}}=-2\gamma {\int }_{\Omega }\left(I_0\left(x-u^{\left(n\right)}\left(x,T\right)\right)-I_1\left(x\right)\right){\nabla }_u{I}_0\left(x-u^{\left(n\right)}\left(x,T\right)\right)·\frac{\partial u^{\left(n\right)}\left(x,T\right)}{\partial v_k}dx& \left(52\right)\end{array}$

[0155] STEP 3: Set n←n+1, and return to Step 1.

[0156] The velocity field can be constructed with various boundary conditions, for example ν(x,t)=0, xεδΩ and tε[0, T], u(x, 0)=ν(x, 0)=0. The differential operator L can be chosen to be any in a class of linear differential operators; we have used operators of the form (−aΔ−b∇∇·cI)^p, p≧1. The operators are 3×3 matrices $\begin{array}{cc}{\nabla }^{T}u=\left(\begin{array}{ccc}\frac{\partial u_1}{\partial x_1}& \frac{\partial u_1}{\partial x_2}& \frac{\partial u_1}{\partial x_3}\\ \frac{\partial u_2}{\partial x_1}& \frac{\partial u_2}{\partial x_2}& \frac{\partial u_2}{\partial x_3}\\ \frac{\partial u_2}{\partial x_1}& \frac{\partial u_2}{\partial x_2}& \frac{\partial u_2}{\partial x_3}\end{array}\right)& \left(53\right)\\ \begin{array}{c}L=-a\Delta -b\nabla {\nabla }^{*}+cI\\ =\left(\begin{array}{ccc}-a\Delta -b\frac{{\partial }^2}{\partial x_1^2}+c& -b\frac{{\partial }^2}{\partial x_1{x}_2}& -b\frac{{\partial }^2}{\partial x_1{x}_3}\\ -b\frac{{\partial }^2}{\partial x_2{x}_1}& -a\Delta -b\frac{{\partial }^2}{\partial x_2^2}+c& -b\frac{{\partial }^2}{\partial x_2{x}_3}\\ -b\frac{{\partial }^2}{\partial x_3{x}_1}& -b\frac{{\partial }^2}{\partial x_3{x}_2}& -a\Delta -b\frac{{\partial }^2}{\partial x_3^2}+c\end{array}\right)\end{array}& \left(54\right)\end{array}$

[0157] 2.3 Small Deformation Solution

[0158] The large deformation computer algorithms can be related to the small deformation approach described in the '212 patent, by choosing ν=u, so that T=δ small, then approximate I−∇u(·,σ)≈I for σε[0, δ), then û(x, δ)=σ(x)δ and defining u(x)≈u(x, δ), so that Eqn. 55 reduces to the small deformation problems described in the '212 patent: $\begin{array}{cc}\hat{u}\left(x\right)=\arg \underset{v_n,k=1,...}{\min }{\int }_{\Omega }{\Vert \mathrm{Lu}\left(x\right)\Vert }^{2}dxdt+D\left(u\right).& \left(55\right)\end{array}$

[0159] 3. Composing Large Deformation Transformation Unifying Landmark and Image Matching

[0160] The approach for generating a hierarchical transformation combining information is to compose the large deformation transformations which are diffeomorphisms and can therefore be composed, h=h_n^o. . . h₂o_h1. Various combinations of transformations may be chosen, including the affine motins, rigid motions generated from subgroups of the generalized linear group, large deformation landmark transformations which are diffeomorphisms, or the high dimensional large deformation image matching transformation (the dimension of the transformations of the vector fields are listed in increasing order). Since these are all diffeomorphisms, they can be composed.

[0161] 4. Fast Method for Landmark Deformations Given Small Numbers of Landmarks

[0162] For small numbers of landmarks, we re-parameterize the problem using the Lagrangian frame, discretizing the optimization over space time Ω×T, into N functions of time T. This reduces the complexity by an order of magnitude given by the imaging lattice [Ω].

[0163] For this, define the Lagrangian positions of the N-landmarks x_i, i=1, . . . , N as they flow through time φ_i(·), i=1, . . . , N, with the associated 3N-vector $\begin{array}{cc}\Phi \left(t\right)=\underset{3N\times 1}{\underset{{\backslash }_{{\mathrm{——}}_{\mathrm{\backslash /}}{\mathrm{——}}^{/}}}{\left(\begin{array}{c}\varphi \left(X_1,t\right)\\ \varphi \left(X_2,t\right)\\ .\\ .\\ .\\ \varphi \left(X_N,t\right)\end{array}\right)}}& \left(56\right)\end{array}$

[0164] The particle flows Φ(t) are defined by the velocities ν(·) according to the fundamental O.D.E. $\begin{array}{cc}\frac{d\varphi \left(t,x\right)}{dt}=v\left(\varphi \left(t,x\right),t\right).& \left(57\right)\end{array}$

[0165] It is helpful to define the 3N×3N covariance matrix K(Φ(t)): $\begin{array}{cc}K\left(\Phi \left(t\right)\right)\underset{\underset{3N\times 3N}{{\backslash }_{{\mathrm{———————————————————}}_{\mathrm{\backslash /}}{\mathrm{————————————————————}}^{/}}}}{\stackrel{.}{=}\left(\begin{array}{cccc}K\left(\varphi \left(t,x_1\right),\varphi \left(t,x_1\right)\right)& K\left(\varphi \left(t,x_1\right),\varphi \left(t,x_2\right)\right)& ...& K\left(\varphi \left(t,x_1\right),\varphi \left(t,x_N\right)\right)\\ K\left(\varphi \left(t,x_2\right),\varphi \left(t,x_1\right)\right)& K\left(\varphi \left(t,x_2\right),\varphi \left(t,x_2\right)\right)& ...& K\left(\varphi \left(t,x_2\right),\varphi \left(t,x_N\right)\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ K\left(\varphi \left(t,x_N\right),\varphi \left(t,x_1\right)\right)& K\left(\varphi \left(t,x_N\right),\varphi \left(t,x_2\right)\right)& ...& K\left(\varphi \left(t,x_N\right),\varphi \left(t,x_N\right)\right)\end{array}\right).}& \left(58\right)\end{array}$