Motion field generation for non-rigid image registration

• US 20060002632A1
• Filed: 06/30/2004
• Published: 01/05/2006
• Est. Priority Date: 06/30/2004
• Status: Active Grant

• Current Assignee: Accuray Incorporated
• Original Assignee: Accuray Incorporated
• Inventors: Kuduvalli, Gopinath, Fu, Dongshan
• Granted Patent: US 7,426,318 B2
• Time in Patent Office: Days
• US Class Current: 382/294
• CPC Class Codes: G06T 2207/10081 Computed x-ray tomography [CT]
• CPC Class Codes: G06T 2207/10121 Fluoroscopy
• CPC Class Codes: G06T 2207/10124 Digitally reconstructed radiograph [DRR]
• CPC Class Codes: G06T 2207/30012 Spine; Backbone
• CPC Class Codes: G06T 7/344 involving models
• CPC Class Codes: G06V 10/754 involving a deformation of the sample pattern or of the reference pattern; Elastic matching

Abstract

A method and system is presented for generating a full motion field of an object, when performing non-rigid image registration between a first and a second image of the object. A mesh grid generator generates in the first image a mesh grid having a plurality of mesh nodes, for each of a hierarchy of mesh resolution levels. A nodal motion estimator estimates, for each mesh node, at least one nodal motion vector that describes a matching of the mesh node with a corresponding mesh node in the second image. At each mesh node, multi-level block matching is performed, using a similarity measure based on pattern intensity. A smoothness constraint is imposed, in order to reconstruct the nodal motion vector for those mesh nodes at which mismatching occurs. A motion field interpolator determines a local motion vector for any desired point of interest within the first image, by interpolating from the nodal motion vectors of the surrounding mesh nodes.

Claims

1. An apparatus for generating a motion field of an object when registering a first image of the object with a second image of the object, the apparatus comprising:

a mesh grid generator configured to generate in the first image a mesh grid having a plurality of mesh nodes;

a nodal motion estimator configured to estimate, for each mesh node in the first image, at least one nodal motion vector that describes a matching of said mesh node with a corresponding mesh node in the second image; and

a motion field interpolator configured to determine a local motion vector for each of a plurality of points of interest within the first image, by interpolating from the nodal motion vectors of the mesh nodes that surround each point of interest.

2. An apparatus in accordance with claim 1,

wherein the mesh grid generator is configured to repeat, at each of a plurality of mesh resolution levels, the act of generating of said mesh grid;

wherein nodal motion estimator is configured to repeat, at each of said plurality of mesh resolution levels, the act of estimating the at least one nodal motion vector for each mesh node; and

wherein the nodal motion vectors from which the motion field interpolator interpolates to determine said local motion vector comprise nodal motion vectors that have been estimated at a final one of said plurality of mesh resolution levels.

3. An apparatus in accordance with claim 2,

further comprising a motion field reconstructor configured to reconstruct the nodal motion vector for a mesh node if any nodal motion vector estimated for any mesh node is unreliable.

4. An apparatus in accordance with claim 2, wherein at each mesh resolution level, the nodal motion estimator is configured to pass onto a subsequent mesh resolution level one or more nodal motion vectors that have been estimated at the current mesh resolution level.

5. An apparatus in accordance with claim 2, wherein said plurality of mesh levels comprise successively increasing mesh resolution levels, the mesh grid at each successive mesh resolution level having a number of mesh nodes that is greater compared to the number of mesh nodes at each previous mesh resolution level.

6. An apparatus in accordance with claim 5,

wherein for the first one of the plurality of mesh resolution levels, the nodal motion estimator is configured to estimate a global translation for the first image and to use the global translation as an estimate for the nodal motion vector for each mesh node in the mesh grid generated for the first mesh resolution level;

and

wherein the global translation represents a translation of an image center of the first image.

7. An apparatus in accordance with claim 2,

wherein at each mesh resolution level subsequent to the first mesh resolution level, the mesh nodes in the first image comprise a first subset of mesh nodes and a second subset of mesh nodes;

wherein said first subset of mesh nodes comprise mesh nodes generated at a previous mesh resolution level; and

wherein said second subset of mesh nodes comprise mesh nodes that are newly added at the current mesh resolution level.

8. An apparatus in accordance with claim 7,

wherein at each mesh resolution level subsequent to the first mesh resolution level, said nodal motion estimator is configured to use the nodal motion vectors passed on from a previous mesh resolution level as estimates for nodal motion vectors for the first subset of mesh nodes;

wherein said nodal motion estimator comprises an interpolator configured to interpolate from the nodal motion vectors estimated for the first subset of mesh nodes, to estimate the nodal motion vectors for the second subset of mesh nodes; and

wherein said nodal motion estimator further comprises a nodal motion refiner configured to refine, for each mesh node in both the first and the second subsets, the nodal motion vector that was estimated for said mesh node.

9. An apparatus in accordance with claim 8,

wherein the nodal vector refiner is configured to:

define a block centered on a mesh node in the first image;

search for a matching mesh node in the second image that maximizes a similarity measure between the block in the first image and another block centered around the matching mesh node in the second image; and

modify the nodal motion vector so that the nodal motion vector describes a mapping of said mesh node in the first image onto said matching mesh node in the second image.

10. An apparatus in accordance with claim 9,

wherein the nodal vector refiner is configured to repeat, for each of a plurality of successively increasing image resolution levels, the acts of defining a block centered around a mesh node in the first image, searching for a matching mesh node in the second image, and refining the nodal motion vector.

11. An apparatus in accordance with claim 10, wherein the size of said block decreases at each increasing image resolution level.

12. An apparatus in accordance with claim 9,

wherein said block matching subsystem is configured to define a search window around said mesh node, when searching for said matching mesh node, and to search for the location of said maximum within said window.

13. An apparatus in accordance with claim 12, wherein said block matching subsystem is configured to review a plurality of search candidates when searching for the location of said maximum within said search window.

14. An apparatus in accordance with claim 9, wherein the similarity measure is based on a pattern intensity function.

15. An apparatus in accordance with claim 14,

wherein the first and second images are discretized images each characterized by an array of pixels, each pixel having an associated pixel value representative of the intensity of the image at a surface unit area corresponding to the pixel;

and

wherein a mathematical formulation of the pattern intensity function comprises;

$\underset{i,j}{\sum <br><br>}<br><br>\underset{k,l⋐<br><br>R}{\sum <br><br>}<br><br>\frac{{\mathrm{\sigma <br><br>}}^2}{{\mathrm{\sigma <br><br>}}^2+{\left(I_{\mathrm{diff}}<br><br>\left(i,j\right)-I_{\mathrm{diff}}<br><br>\left(i+k,j+l\right)\right)}^2},$

wherein I_diff(i,j)=I_second(i,j)−

I_first(i,j),

wherein I_diff(i,j) represents an array of pixel values for a difference image formed by subtracting the corresponding pixel values of the second image from each pixel value of the first image,

wherein I_first(i,j) represents the (i,j)-th pixel value of the first image of the object, I_second(i,j) represents the (i,j)-th pixel value of the second image of the object,

wherein σ

is a weighting constant for weighting the pattern intensity function, and

wherein R is a neighborhood defined around the (i,j)-the pixel as a center point.

16. An apparatus in accordance with claim 3,

wherein the nodal motion estimator is configured to determine whether a mismatch occurred between any mesh node and its corresponding mesh node.

17. An apparatus in accordance with claim 16, wherein the nodal motion estimator is configured to compute a difference between a nodal motion vector for a mesh node, and a median nodal motion vector for said mesh node, said median nodal motion vector comprising the mathematical median of the nodal motion vectors estimated for the mesh nodes that surround said mesh node;

wherein the nodal motion estimator is configured to determine that the nodal motion vector for said mesh node is reliable, if said difference is less than a predefined threshold;

wherein the nodal motion estimator is configured to determine that the nodal motion vector for said mesh node is unreliable, if said difference is greater than or equal to said predefined threshold; and

wherein the nodal motion estimator is configured to replace any unreliable nodal motion vector with said median nodal motion vector.

18. An apparatus in accordance with claim 2,

wherein the motion field reconstructor is configured to impose a smoothness constraint on the nodal motion vectors estimated by the nodal motion estimator.

19. An apparatus in accordance with claim 18,

wherein the motion field reconstructor is further configured to formulate a cost function, and to minimize said cost function in order to impose said smoothness constraint.

20. An apparatus in accordance with claim 19,

wherein a mathematical expression characterizing the cost function comprises:

E(d)=∫

∫

β

(d−

u)² dxdy+λ

∫

∫

(d_,x²+d_,y²)dxdy,

wherein

d represents a desired local estimate for a nodal motion vector at coordinates (x,y),

u represents a locally estimated nodal motion vector at coordinates (x,y),

β

represents a reliability constant that ranges from 0 to 1, 0 indicating an unreliable estimation, and 1 indicating a reliable estimation, and

E(d) represents said cost function.

21. An apparatus in accordance with claim 20, wherein a discretized form of said mathematical expression for said cost function comprises:

E(d_i,j)=Σ

Σ

β

_i,j(d_i,j−

u_i,j)²+λ

Σ

Σ

└

(d_i,j−

d_i-1,j)²+(d_i,j−

d_i,j)²┘

wherein

u_i,j represents the locally estimated nodal motion vector for the ij-th mesh node,

d_i,j represent the desired local estimate for the nodal motion vector for the ij-th mesh node,

the value of β

_i,j is 1 if the estimation is reliable,

the value of β

_i,j is 0 if the estimation is unreliable, and

wherein the summations over the indices i and j are over the total number of mesh nodes.

22. An apparatus in accordance with claim 21, wherein the motion field reconstructor is configure to minimize said cost function by solving a system of simultaneous linear equations.

23. An apparatus in accordance with claim 22, wherein a mathematical expression characterizing said system of simultaneous linear equations comprises:

$\begin{array}{c}\frac{\mathrm{\delta <br><br>}<br><br>\text{\hspace{1em}<br><br>}<br><br>E<br><br>\left(d_{i,j}\right)}{\partial <br><br>d_{i,j}}=\left({\mathrm{\beta <br><br>}}_{i,j}+4<br><br>\mathrm{\lambda <br><br>}\right)<br><br>d_{i,j}-\mathrm{\lambda <br><br>}<br><br>\left(d_{i-1,j}+d_{i+1,j}+d_{i,j-1}+d_{i,j+1}\right)-{\mathrm{\beta <br><br>}}_{i,j}<br><br>u_{i,j}\\ \text{\hspace{1em}<br><br>}<br><br>=0,\end{array}$

where the indices i and j range from 0 to the total number of mesh nodes.

24. An apparatus in accordance with claim 1,

wherein a mathematical expression of the local motion vector for a point of interest comprises (dx,dy), dx and dy representing translations along mutually orthogonal x- and y-directions, respectively;

wherein a mathematical expression of the nodal motion vector for the i-th node that surrounds said point of interest comprises(dx_n(i),dy_n(i)); and

wherein said motion field interpolator is configured to interpolate the local motion vector from the nodal motion vectors estimated for the mesh nodes that surround said point of interest by performing a summation described by;

$\begin{array}{c}\mathrm{dx}=\underset{i=1}{\overset{m}{\sum <br><br>}}<br><br>N<br><br>\left(i\right)<br><br>\text{\hspace{1em}<br><br>}<br><br>\mathrm{dx}<br><br>\left(i\right),\\ \mathrm{dy}=\underset{i=1}{\overset{m}{\sum <br><br>}}<br><br>N<br><br>\left(i\right)<br><br>\text{\hspace{1em}<br><br>}<br><br>\mathrm{dy}<br><br>\left(i\right),\end{array}$

where N(i) is a shape function for the i-th mesh node surrounding said point of interest, and m is the total number of mesh nodes that surround said point of interest.

25. An apparatus in accordance with claim 24,

wherein said point of interest is characterized by coordinates (ξ

,η

) along mutually orthogonal ξ

- and η

-axes, respectively,

wherein the total number m of mesh nodes surrounding the point of interest is 9; and

wherein the shape functions N(i), for i=1, . . . , 9, are given by;

$\begin{array}{c}N<br><br>\left(1\right)=\left(1-\mathrm{\xi <br><br>}\right)<br><br>\left(1-\mathrm{\eta <br><br>}\right)/4-\left(N_8+N_5\right)/2,\\ N<br><br>\left(2\right)=\left(1-\mathrm{\xi <br><br>}\right)<br><br>\left(1-\mathrm{\eta <br><br>}\right)/4-\left(N_5+N_6\right)/2,\\ N<br><br>\left(3\right)=\left(1+\mathrm{\xi <br><br>}\right)<br><br>\left(1+\mathrm{\eta <br><br>}\right)/4-\left(N_6+N_7\right)/2,\\ N<br><br>\left(4\right)=\left(1-\mathrm{\xi <br><br>}\right)<br><br>\left(1+\mathrm{\eta <br><br>}\right)/4-\left(N_7+N_8\right)/2,\\ N<br><br>\left(5\right)=\left(1-{\mathrm{\xi <br><br>}}^2\right)<br><br>\left(1-\mathrm{\eta <br><br>}\right)/2,\\ N<br><br>\left(6\right)=\left(1-\mathrm{\xi <br><br>}\right)<br><br>\left(1-{\mathrm{\eta <br><br>}}^2\right)/2,\\ N<br><br>\left(7\right)=\left(1-{\mathrm{\xi <br><br>}}^2\right)<br><br>\left(1+\mathrm{\eta <br><br>}\right)/2,\\ N<br><br>\left(8\right)=\left(1-\mathrm{\xi <br><br>}\right)<br><br>\left(1-{\mathrm{\eta <br><br>}}^2\right)/2,\\ N<br><br>\left(9\right)=\left(1-{\mathrm{\xi <br><br>}}^2\right)<br><br>\left(1-{\mathrm{\eta <br><br>}}^2\right).\end{array}$

26. An apparatus in accordance with claim 1, wherein said one or more points of interest are contained within an ROI (region-of-interest) within said first image.

27. An apparatus in accordance with claim 26, wherein said ROI comprises an ROI that has been selected to maximize an image entropy H within said first image.

28. An apparatus in accordance with claim 27, wherein said image entropy H is a modified Shannon entropy given by:

H=−

Σ

_I I P(I) log P(I),

where I is the value of the intensity of the image, at each pixel of the image, and P(I) is the probability of an image intensity value I occurring within the ROI.

29. An apparatus in accordance with claim 1, wherein the motion field is representative of at least one of a translation, a rotation, a rigid deformation, and a non-rigid deformation of the object, between the acquisition of the first image and the acquisition of the second image.

30. An apparatus in accordance with claim 1, wherein said local motion vector comprises at least one of a local translation vector and a local rotation vector.

31. An apparatus in accordance with claim 1, wherein said first image comprises a pre-operative DRR (digitally reconstructed radiograph) of an anatomical region, and said second image comprises an intra-operative x-ray image of said anatomical region.

32. An apparatus for generating a motion field for an object when registering a first image of the object with a second image of the object, the apparatus comprising:

means for generating in the first image a mesh grid having a plurality of mesh nodes;

means for estimating, for each mesh node in the first image, at least one nodal motion vector that describes a matching of said mesh node with a corresponding mesh node in the second image;

means for reconstructing the nodal motion vector for a mesh node if a mismatch occurs between said mesh node and its corresponding mesh node; and

means for determining a local motion vector for each of a plurality of points of interest within the first image, by interpolating from the nodal motion vectors of the mesh nodes that surround each point of interest.

33. An apparatus for generating a motion field representing a motion of an object between an acquisition of a first image of the object and an acquisition of a second image of the object, the apparatus comprising:

a nodal motion estimator configured to estimate, for each of a plurality of mesh nodes and at each of a plurality of mesh resolution levels, at least one nodal motion vector that describes a matching of each node in the first image with its corresponding mesh node in the second image; and

a motion field interpolator configured to compute a local motion vector for each of a plurality of points of interest within the first image, by interpolating from the nodal motion vectors that are estimated at the final mesh resolution level for those mesh nodes that surround said point of interest.

34. An apparatus in accordance with claim 33, further comprising a motion field reconstructor configured to reconstruct the nodal motion vector for any mesh node in which a mismatch occurs between said mesh node and its corresponding mesh node.

35. An apparatus for generating a motion field representing a motion of an object between an acquisition of a first image of the object and an acquisition of a second image of the object, the apparatus comprising:

means for estimating, for each of a plurality of mesh nodes within the first image at each of a plurality of mesh resolution levels, at least one nodal motion vector that describes a matching of each node in the first image with its corresponding mesh node in the second image; and

interpolation means for computing a local motion vector for each of a plurality of points of interest within the first image, by interpolating from the nodal motion vectors that are estimated, at a final one of said mesh resolution levels, for those mesh nodes that surround said point of interest.

36. An apparatus for determining a local motion vector for at least one point of interest within a first image of an object when registering the first image onto a second image of the object, the apparatus comprising:

a mesh grid generator configured to generate, for each of a plurality of mesh resolution levels, a mesh grid within the first image, said mesh grid having a plurality of mesh nodes;

a nodal motion estimator configured to estimate, for each of said plurality of mesh resolution levels, at least one nodal motion vector for each mesh node in the first image, said nodal motion vector describing a matching of said mesh node onto a corresponding mesh node in the second image; and

an interpolator configured to interpolate from the nodal motion vectors that are estimated for the mesh nodes that surround the at least one point of interest, so as to determine said local motion vector.

37. A computer-readable medium having stored therein computer-readable instructions for a processor, wherein the instructions, when read and implemented by the processor, cause the processor to:

input and store for a first image of an object data representative of a mesh grid having a plurality of mesh nodes, for each of a plurality of mesh resolution levels;

estimate, for each mesh node in each mesh resolution level, at least one nodal motion vector that describes a matching of said mesh node onto a corresponding mesh node in a second image; and

compute a local motion vector for one or more points of interest in the first image by interpolating from the nodal motion vectors estimated at a final mesh resolution level for the mesh nodes that surround each point of interest.

38. A computer-readable medium in accordance with claim 37, the computer-readable medium having stored therein further computer-readable instructions for the processor,

wherein said further instructions, when read and implemented by the processor, cause the processor to detect any mesh node for which a mismatch occurs between said mesh node and its corresponding mesh node, and to reconstruct the nodal motion vector for said detected mesh node by imposing a smoothness constraint.

39. An apparatus for generating a motion field representing a motion of an object between an acquisition of a first image of the object and an acquisition of a second image of the object, the apparatus comprising:

a local motion determiner configured to determine a local motion vector for one or more points of interest within the first image by interpolating from nodal motion vectors estimated for a plurality of mesh nodes surrounding each point of interest.

40. A block matching system for performing multi-level block matching during image registration between a first image of an object and a second image of the object, the block matching system comprising:

means for defining a block centered around each of a plurality of mesh nodes that delineate a mesh grid in the first image; and

means for attributing, to each mesh node, a matching mesh node in the second image that maximizes a similarity measure between said block in the first image and another block centered around the matching mesh node in the second image.

41. A method of generating a motion field for an object when registering a first image of the object with a second image of the object, the method comprising:

generating in the first image a mesh grid having a plurality of mesh nodes;

for each mesh node in the first image, estimating at least one nodal motion vector that describes a matching of said mesh node with a corresponding mesh node in the second image; and

determining a local motion vector for each of a plurality of points of interest within the first image, by interpolating from the nodal motion vectors of the mesh nodes that surround the point of interest.

42. A method in accordance with claim 41,

further comprising repeating, at each of a plurality of mesh resolution levels, the acts of generating in the first image a mesh grid having a plurality of mesh nodes and estimating at least one nodal motion vector for each mesh node;

and

wherein the nodal motion vectors of the mesh nodes that surround the point of interest comprise nodal motion vectors that have been estimated at a final one of the plurality of mesh resolution levels.

43. A method in accordance with claim 42, further comprising:

if any nodal motion vector that was estimated for any mesh node is unreliable, reconstructing the nodal motion vector for the mesh node.

44. A method in accordance with claim 42, further comprising passing, at each mesh resolution level, one or more nodal motion vectors that have been estimated at the current mesh resolution level onto a subsequent mesh resolution level.

45. A method in accordance with claim 42,

wherein said plurality of mesh resolution levels comprise successively increasing mesh resolution levels, the mesh grid at each successive mesh resolution level having a number of mesh nodes that is greater compared to the number of mesh nodes at each previous mesh resolution level.

46. A method in accordance with claim 45,

wherein for the first one of the plurality of mesh resolution levels, the act of estimating the at least one nodal motion vector comprises:

estimating a global translation for the first image, the global translation representing a translation of an image center of the first image; and

using the global translation as an estimate for the nodal motion vector for each mesh node in the mesh grid generated for the first mesh resolution level.

47. A method in accordance with claim 42,

wherein at each mesh resolution level subsequent to the first mesh resolution level, the mesh modes in the first image comprise a first subset of mesh nodes and a second subset of mesh nodes;

wherein said first subset of mesh nodes comprise mesh nodes generated at a preceding mesh resolution level; and

wherein said second subset of mesh nodes comprise mesh nodes that are newly added at a current mesh resolution level.

48. A method in accordance with claim 47,

wherein at each mesh resolution level subsequent to the first mesh resolution level, the act of estimating the at least one nodal motion vector comprises:

i) using the nodal motion vectors passed on from the previous mesh resolution level as estimates for nodal motion vectors for the first subset of mesh nodes;

ii) interpolating from the estimated nodal motion vectors for the first subset of mesh nodes to estimate the nodal motion vectors for the second subset of mesh nodes; and

iii) for each mesh node in both the first and second subsets in the first image, refining the nodal vector that was estimated for the mesh node.

49. A method in accordance with claim 48, wherein refining the estimated nodal motion vector for a mesh node in the first image comprises:

defining a block centered around the mesh node in the first image;

searching for a matching mesh node in the second image that maximizes a similarity measure between the block in the first image and another block centered around the matching mesh node in the second image; and

revising the local motion vector so that the nodal motion vector describes a mapping of said mesh node in the first image onto said matching mesh node in the second image.

50. A method in accordance with claim 49, wherein refining the estimated nodal motion vector for the mesh node in the first image further comprises:

repeating, for each of a plurality of successively increasing image resolution levels, the acts of defining a block centered around the mesh node in the first image, searching for a matching mesh node in the second image, and revising the local motion vector so that the nodal motion vector describes a mapping of said mesh node in the first image onto said matching mesh node in the second image.

51. A method in accordance with claim 50, wherein the size of said block decreases at each increasing image resolution level.

52. A method in accordance with claim 49, wherein the act of searching for the matching mesh node comprises:

defining a search window around said mesh node; and

searching within said search window for the location of a maximum of said similarity measure.

53. A method in accordance with claim 52, wherein searching for the location of the maximum within the search window comprises reviewing a plurality of search candidates.

54. A method in accordance with claim 49, wherein the similarity measure is based on a pattern intensity function.

55. A method in accordance with claim 54,

wherein the first and second images are discretized images each characterized by an array of pixel values, each pixel value being a number representative of the intensity of the image at a corresponding surface unit area;

and

wherein a mathematical formulation of the pattern intensity function comprises;

$\underset{i,j}{\sum <br><br>}<br><br>\underset{k,l⋐<br><br>R}{\sum <br><br>}<br><br>\frac{{\mathrm{\sigma <br><br>}}^2}{{\mathrm{\sigma <br><br>}}^2+{\left(I_{\mathrm{diff}}<br><br>\left(i,j\right)-I_{\mathrm{diff}}<br><br>\left(i+k,j+l\right)\right)}^2},$

wherein I_diff(i,j)=I_second(i,j)−

I_first(i,j),

wherein I_diff(i,j) represents an array of pixel values for a difference image formed by subtracting the corresponding pixel values of the second image from each pixel value of the first image,

wherein I_first(i,j) represents the (i,j)-th pixel value of the first image of the object, I_second(i,j) represents the (i,j)-th pixel value of the second image of the object,

wherein σ

is a weighting constant for weighting the pattern intensity function, and

wherein R is a neighborhood defined around the (i,j)-the pixel as a center point.

56. A method in accordance with claim 43, wherein estimating the at least one nodal motion vector comprises determining whether any estimated nodal motion vector for any mesh node is unreliable.

57. A method in accordance with claim 56, wherein determining whether an estimated nodal motion vector for a mesh node is unreliable comprises:

computing the difference between the nodal motion vector for the mesh node, and a median nodal motion vector, wherein said median nodal motion vector comprises a mathematical median of the nodal motion vectors estimated for the mesh nodes surrounding the mesh node;

if the difference is less than a predefined threshold, determining that the nodal motion vector for the mesh node is reliable;

if the difference is greater than or equal to said predefined threshold, determining that the nodal motion vector for the mesh node is unreliable; and

replacing any unreliable nodal motion vector with said median nodal motion vector.

58. A method in accordance with claim 43, wherein reconstructing the motion field from the local motion vectors further comprises:

imposing a smoothness constraint on the local motion vectors.

59. A method in accordance with claim 58, wherein imposing the smoothness constraint comprises:

formulating a cost function; and

minimizing the cost function.

60. A method in accordance with claim 59, wherein a mathematical expression characterizing the cost function comprises:

E(d)=∫

∫

β

(d−

u)² dxdy+λ

∫

∫

(d_,x²+d_,y²)dxdy,

wherein

d represents a desired local estimate for a nodal motion vector at coordinates (x,y),

u represents a locally estimated nodal motion vector at coordinates (x,y),

β

represents a reliability constant that ranges from 0 to 1, 0 indicating an unreliable estimation, and 1 indicating a reliable estimation, and

E(d) represents said cost function.

61. A method in accordance with claim 60, wherein a discretized form of said mathematical expression for said cost function comprises:

E(d_i,j)=Σ

Σ

β

_i,j(d_i,j−

u_i,j)²+λ

Σ

Σ

└

(d_i,j−

d_i,j)²⁺⁽d_i,j−

d_i,j)²┘

wherein

u_i,j represents the locally estimated nodal motion vector for the ij-th mesh node,

d_i,j represent the desired local estimate for the nodal motion vector for the ij-th mesh node,

the value of β

_i,j is 1 if the estimation is reliable,

the value of β

_i,j is 0 if the estimation is unreliable,

and wherein the summations over the indices i and j are over the total number of mesh nodes.

62. A method in accordance with claim 61, wherein minimizing the cost function comprises solving a system of simultaneous linear equations.

63. A method in accordance with claim 62, wherein a mathematical expression characterizing the system of simultaneous linear equations comprises:

$\begin{array}{c}\frac{\mathrm{\delta <br><br>}<br><br>\text{\hspace{1em}<br><br>}<br><br>E<br><br>\left(d_{i,j}\right)}{\partial <br><br>d_{i,j}}=\left({\mathrm{\beta <br><br>}}_{i,j}+4<br><br>\mathrm{\lambda <br><br>}\right)<br><br>d_{i,j}-\mathrm{\lambda <br><br>}<br><br>\left(d_{i-1,j}+d_{i+1,j}+d_{i,j-1}+d_{i,j+1}\right)-{\mathrm{\beta <br><br>}}_{i,j}<br><br>u_{i,j}\\ \text{\hspace{1em}<br><br>}<br><br>=0,\end{array}$

where the indices i and j range from 0 to the total number m of nodes that surround the point of interest.

64. A method in accordance with claim 41,

wherein a mathematical expression of the local motion vector for a point of interest comprises (dx,dy), dx and dy representing translations along mutually orthogonal x- and y-directions, respectively;

wherein a mathematical expression of the nodal motion vector for the i-th node that surrounds the point of interest comprises(dx_n(i),dy_n(i)); and

wherein interpolating the local motion vector from the nodal motion vectors estimated for the mesh nodes that surround the point of interest comprises performing a summation described by;

$\begin{array}{c}\mathrm{dx}=\underset{i=1}{\overset{m}{\sum <br><br>}}<br><br>N<br><br>\left(i\right)<br><br>\text{\hspace{1em}<br><br>}<br><br>\mathrm{dx}<br><br>\left(i\right),\\ \mathrm{dy}=\underset{i=1}{\overset{m}{\sum <br><br>}}<br><br>N<br><br>\left(i\right)<br><br>\text{\hspace{1em}<br><br>}<br><br>\mathrm{dy}<br><br>\left(i\right),\end{array}$

where N(i) is a shape function for the i-th mesh node surrounding the point of interest, and m is the total number of mesh nodes that surround the point of interest.

65. A method in accordance with claim 64,

wherein the point of interest is characterized by coordinates (ξ

,η

) along mutually orthogonal ξ

- and η

-axes, respectively,

wherein the total number m of mesh nodes surrounding the point of interest is 9; and

wherein the shape functions N(i), for i=1, . . . , 9, are given by;

$\begin{array}{c}N<br><br>\left(1\right)=\left(1-\mathrm{\xi <br><br>}\right)<br><br>\left(1-\mathrm{\eta <br><br>}\right)/4-\left(N_8+N_5\right)/2,\\ N<br><br>\left(2\right)=\left(1-\mathrm{\xi <br><br>}\right)<br><br>\left(1-\mathrm{\eta <br><br>}\right)/4-\left(N_5+N_6\right)/2,\\ N<br><br>\left(3\right)=\left(1+\mathrm{\xi <br><br>}\right)<br><br>\left(1+\mathrm{\eta <br><br>}\right)/4-\left(N_6+N_7\right)/2,\\ N<br><br>\left(4\right)=\left(1-\mathrm{\xi <br><br>}\right)<br><br>\left(1+\mathrm{\eta <br><br>}\right)/4-\left(N_7+N_8\right)/2,\\ N<br><br>\left(5\right)=\left(1-{\mathrm{\xi <br><br>}}^2\right)<br><br>\left(1-\mathrm{\eta <br><br>}\right)/2,\\ N<br><br>\left(6\right)=\left(1-\mathrm{\xi <br><br>}\right)<br><br>\left(1-{\mathrm{\eta <br><br>}}^2\right)/2,\\ N<br><br>\left(7\right)=\left(1-{\mathrm{\xi <br><br>}}^2\right)<br><br>\left(1+\mathrm{\eta <br><br>}\right)/2,\\ N<br><br>\left(8\right)=\left(1-\mathrm{\xi <br><br>}\right)<br><br>\left(1-{\mathrm{\eta <br><br>}}^2\right)/2,\\ N<br><br>\left(9\right)=\left(1-{\mathrm{\xi <br><br>}}^2\right)<br><br>\left(1-{\mathrm{\eta <br><br>}}^2\right).\end{array}$

66. A method in accordance with claim 41, wherein said one or more points of interest are contained within an ROI (region-of-interest) within said first image.

67. A method in accordance with claim 66, wherein said ROI comprises an ROI that has been selected to maximize an image entropy H within said first image.

68. A method in accordance with claim 67, wherein said image entropy H is a modified Shannon entropy given by:

H=Σ

_I I P(I) log P(I),

where I is the value of the intensity of the image, at each pixel of the image, and P(I) is the probability of an image intensity value I occurring within the ROI.

69. A method in accordance with claim 41, wherein the motion field is representative of at least one of a translation, a rotation, a rigid deformation, and a non-rigid deformation of the object, between the acquisition of the first image and the acquisition of the second image.

70. A method in accordance with claim 41, wherein said local motion vector comprises at least one of a local translation vector and a local rotation vector.

71. A method in accordance with claim 41, wherein said first image comprises a pre-operative DRR of an anatomical region, and said second image comprises an intra-operative x-ray image of said anatomical region.

72. A method of determining a local motion vector for at least one point of interest within a first image of an object when registering the first image onto a second image of the object, the method comprising:

A) for each of a plurality of mesh resolution levels,

a) generating within the first image a mesh grid having a plurality of mesh nodes, the point of interest being located in the interior of the mesh grid;

b) for each mesh node in the first image, estimating at least one nodal motion vector, the nodal motion vector describing a matching of said mesh node with a corresponding mesh node in the second image; and

c) passing one or more estimated nodal motion vectors onto a subsequent mesh resolution level; and

B) deriving the local motion vector at the point of interest by interpolating from the nodal motion vectors that have been estimated for the mesh nodes that surround the point of interest at a final one of said plurality of mesh resolution levels.