Fiducial-less tracking with non-rigid image registration

• US 7,327,865 B2
• Filed: 06/30/2004
• Issued: 02/05/2008
• Est. Priority Date: 06/30/2004
• Status: Active Grant

• Current Assignee: Accuray Incorporated
• Original Assignee: Accuray Incorporated
• Inventors: Kuduvalli, Gopinath, Fu, Dongshan
• Primary Examiner(s): AZARIAN, SEYED H
• Application Number:
  US10/880,486
• Publication Number:
  US 20060002630A1
• Time in Patent Office: 1,315 Days
• Field of Search: 382/128, 382/129, 382/130, 382/131, 382132-133, 382/154, 382/168, 382/181, 382/209, 382/219, 382/232, 382/255, 382/260, 382/274, 382/285, 382/291, 382/294, 382/305, 382/318, 600/411, 600/130, 600/407, 375/240.16, 378/8, 378/65, 378/28
• US Class Current: 382/128
• 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/24 Aligning, centring, orientation detection or correction of the image
• CPC Class Codes: G06V 2201/03 Recognition of patterns in medical or anatomical images

Abstract

A method and system is presented for tracking patient motion in image guided surgery, using skeletal reference structures instead of fiducial markers. A non-rigid 2D/3D image registration is performed between pre-operative 3D scan data and intra-operative 2D images. The CT numbers are modified to compensate for the difference between the ratio of bone-to-tissue attenuation at the CT scan energy level, and the x-ray projection energy level. DRRs are generated from the modified CT numbers. An ROI, within which the image registration is restricted, is defined within the DRRs by calculating a modified Shannon entropy for each image, and seeking the region in which the image entropy is maximized. The DRRs and 2D images are enhanced by applying a top hat filter to increase the visibility of the skeletal structures. The 2D/3D transformation parameters are determined by generating a full motion field that accounts for non-rigid motions and deformations. The full motion field is derived by estimating a plurality of local motion fields, using multi-level block matching and a similarity measure based on pattern intensity.

Claims

1. A method of registering a first image of an object with data representative of a second image of the object, the object including one or more structures, the method comprising:

reconstructing from the data at least one DRR (digitally reconstructed radiograph);

enhancing the DRR and the first image so as to increase the visibility in the DRR and the first image of at least one structure;

selecting an ROI (region of interest) within the DRR, wherein the ROI includes at least one structure in the object;

estimating one or more 2D local motion fields within the ROI, and constructing a full 3D motion field from the local motion fields;

and determining from the full 3D motion field a set of non-rigid transformation parameters that represent the difference in the position and orientation of at least one structure as shown in the first image, as compared to the position and orientation of the at least one structure as shown in the DRR;

wherein the first image comprises a 2D x-ray projection image, the data comprise 3D scan data, the object comprises an anatomical region, and the one or more structures comprises at least one reference structure and at least one target.

2. A method in accordance with claim 1, wherein the act of estimating the one or more 2D local motion fields within the ROI comprises maximizing a similarity measure between the DRR and the x-ray image, for each of a plurality of image resolution levels.

3. A method in accordance with claim 2, wherein the measure of similarity is based on pattern intensity.

4. A method in accordance with claim 3,

wherein the pixel intensity values at the i-th row and j-th column of each array of pixels for the near real time 2D x-ray projection image and the DRR can be represented as I_Live(i,j) and I_DRR(i,j), respectively;

and

wherein determining the measure of similarity between the DRR and the near real time 2D x-ray projection image comprises comparing the pixel intensity values I_Live(i,j) and I_DRR(i,j).

5. A method in accordance with claim 4, wherein the step of determining a measure of similarity between the DRR and the x-ray image comprises:

a) forming a difference image by subtracting a corresponding pixel value of the DRR from each pixel value of the x-ray image, the pixel value at the i-th row and j-th column of the array of pixel values for the difference image being given by;

I_diff(i,j)=I_Live(i,j)−

I_DRR(i,j);

and

b) applying upon each pixel of the difference image a pattern intensity function, wherein the pattern intensity function is defined by summing asymptotic functions of the gradients of the difference image over the pixels within a neighborhood R, and

wherein the neighborhood R is defined so that the gradients of the difference image can be considered in at least four directions.

6. A method in accordance with claim 5, wherein the pattern intensity function is characterized by a mathematical formulation given by:

$\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 σ

is a weighting constant for weighting the pattern intensity function, and where R is a neighborhood defined around a pixel (i,j) as a center point.

7. A method in accordance with claim 5, wherein the at least four directions comprise:

a) a substantially horizontal direction;

b) a substantially vertical direction;

c) a substantially diagonal direction at about 45 degrees; and

d) a substantially diagonal direction at about −

45 degrees.

8. A method in accordance with claim 1,

wherein the 2D x-ray image and the DRR are each characterized by an array of pixels, each pixel having an associated pixel value that represents the image intensity of a corresponding unit of the anatomical region.

9. A method in accordance with claim 8, wherein the step of defining the ROI comprises:

a) calculating an entropy measure H of the DRR; and

b) selecting an ROI so as to maximize the entropy measure H.

10. A method in accordance with claim 9, wherein the entropy measure H is a modified Shannon entropy measure 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.

11. A method in accordance with claim 8,

wherein the 3D scan data comprise a plurality of CT numbers representing the x-ray attenuation intensity of corresponding 3D CT voxels, each CT voxel represents a corresponding volume element of the anatomical region;

wherein each CT voxel is disposed within one of a plurality of axial voxel slices, each axial voxel slice representing a corresponding axial slice of the anatomical region;

wherein each CT number represents the attenuated intensity of an x-ray beam that has been generated at a CT scan energy level and that has traversed a corresponding volume element of the anatomical region;

wherein the 2D x-ray projection image is generated in near real time by transmitting through the target an x-ray imaging beam having a known intensity and a known origination position and angle relative to the target; and

wherein the DRR is a synthetic image representing the radiographic image of the target that would be obtained with the x-ray imaging beam at the known intensity, origination position and angle, if the target were positioned in accordance with the 3D scan data.

12. A method in accordance with claim 11, wherein reconstructing the DRR comprises modifying each of the plurality of CT numbers to compensate for the difference between the ratio of bone-to-tissue attenuation at the scan energy level, and the ratio of bone-to-tissue attenuation at the known intensity level of the x-ray beam.

13. A method in accordance with claim 12, wherein reconstructing the DRR further comprises:

casting a plurality of hypothetical rays through the plurality of 3D CT scan voxels from the known intensity and the known origination position and angle, so that each hypothetical imaging beam traverses a subset of the 3D CT voxels;

for each hypothetical ray, integrating along the ray the CT numbers for the subset of CT voxels; and

projecting the integrated value of the CT numbers onto an imaging plane.

14. A method in accordance with claim 13, wherein the step of integrating the CT numbers along the hypothetical ray comprises:

performing a bi-linear interpolation for the voxels encountered by the ray; and

for each voxel of interest, performing a one-dimensional polynomial interpolation over the voxel of interest and for voxels on each voxel slice adjacent to the voxel slice on which the voxel of interest is situated.

15. A method in accordance with claim 13, wherein the step of modifying the CT numbers comprises performing a mathematical operation derived from a non-linear x-ray attenuation model and given by:

C(x,y,z)=a C₀(x,y,z)e^bC0(x,y,z)

where C(x,y,z) represents the modified CT number of a 3D CT voxel having a location (x,y,z);

a and b represent weighting coefficients; and

C₀(x,y,z) represents the un-modified CT number, based on a linear attenuation model, of a 3D CT voxel having a location (x,y,z).

16. A method in accordance with claim 1, wherein the act of enhancing the DRR and the 2D image comprises applying a filter operator to each image, and

wherein the filter operator comprises a top-hat filter configured to select the pixel having the optimal pixel value from each of at least two different neighborhoods within the image, and to eliminate the remaining pixels in the neighborhoods, thereby transforming each image into an enhanced image.

17. A method in accordance with claim 16, wherein the filter operator is defined mathematically by:

f_e=f+w×

[f−

γ

_B(f)]−

b×

[φ

_B(f)−

f],

where f_e represents an enhanced image,

f represents an original image,

w and b represent weighting coefficients,

γ

_B(f) represents a structural element for the opening of f, and

φ

_B(f) represents a structural element for the closing of f.

18. A method in accordance with claim 1, wherein the reference structure comprises a skeletal structure, and the target comprises at least one of a tumor and a lesion.

19. A method in accordance with claim 1, wherein the non-rigid transformation parameters comprise three translational parameters (x, y, and z), and three rotational parameters (r, p, w);

wherein (x, y, and z) represent the three translations of the treatment target along the directions of three mutually orthogonal x-, y-, and z- axes, respectively; and

wherein (r, p, w) represent the three rotations (roll, pitch, yaw) about the three mutually orthogonal x-, y-, and z-axes.

20. A method in accordance with claim 19, wherein the x-ray projection images are computed for at least two orthogonal projections A and B of the anatomical region onto respective projection image planes.

21. A method in accordance with claim 20, wherein computing the non-rigid transformation parameters comprises:

a) determining a separate set of transformation parameters for each projection x-ray image; and

b) combining the resulting parameters for each projection to obtain the 3D transformation parameters.

22. A method in accordance with claim 21, wherein the 3D transformation parameters are related to the transformation parameters for projections A and B by the following relationship:

x=(x_A+x_B)/2, y=(y_A−

y_B)/√

{square root over (2)}, z=(y_A+y_B)/√

{square root over (2)},

r=θ

_x, p=(θ

_B−

θ

_A)/√

{square root over (2)}, w=(θ

_B+θ

_A)/ √

{square root over (2)}.

23. A method in accordance with claim 1, wherein the act of estimating one or more local motion fields comprises:

for each of a plurality of mesh resolution levels, generating in the first image a mesh grid having a plurality of mesh nodes, and estimating 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; 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 estimated at a final one of the plurality of mesh resolution levels for the mesh nodes that surround the point of interest.

24. An image registration system for registering at least one 2D image of an anatomical region with previously generated 3D scan data of the anatomical region, the anatomical region including at least one treatment target and at least one reference structure, wherein the 2D image is generated in near real time by detecting one or more radiographic imaging beams after the imaging beams have traversed at least a portion of the anatomical region, the imaging beams having known intensities and known positions and angles relative to the anatomical region, the system comprising:

a) means for providing the 3D scan data of the anatomical region;

b) a DRR generator configured to generate at least one DRR (digitally reconstructed radiograph) of the anatomical region, using the 3D scan data and the known locations, angles, and intensities of the imaging beams;

c) an ROI selector configured to select an ROI (region of interest) within the DRR, the ROI including the image of the treatment target and of the reference structure;

d) an image enhancer configured to enhance the DRR and the 2D image by applying a filter operator to the DRR and to the 2D image;

e) a motion field generator configured to generate a 3D full motion field within the ROI by estimating a plurality of local motion fields within the ROI; and

f) a parameter determiner configured to determine from the 3D full motion field a set of non-rigid transformation parameters that represent the difference in the position and orientation of the treatment target as shown in the 2D image, as compared to the position and orientation of the treatment target as shown in the DRR.

25. A system in accordance with claim 24, wherein the 3D scan data comprise at least one of CT scan data, MRI scan data, and PET scan data.

26. A system in accordance with claim 24,

wherein the 3D scan data comprise a plurality of CT numbers representing the image intensity of corresponding 3D CT voxels, each CT voxel representing a corresponding volume element of the anatomical region;

wherein each CT voxel is disposed within one of a plurality of axial voxel slices, each axial voxel slice representing a corresponding axial slice of the anatomical region;

wherein each CT number represents the attenuated intensity of an x-ray CT beam that has been generated at a CT scan energy level and that has traversed the corresponding volume element of the anatomical region; and

wherein the anatomical region comprises at least some bone and at least some tissue.

27. A system in accordance with claim 26,

further comprising a scan data modifier configured to modify the 3D scan data, before the 3D scan data are used by the DRR generator, so as to compensate for the difference between the ratio of bone-to-tissue attenuation at the CT scan energy level, and the ratio of bone attenuation at the known intensity level of the imaging beam.

28. A system in accordance with claim 27,

wherein the scan data modifier includes a processor for performing on each CT number a mathematical operation derived from a non-linear x-ray attenuation model;

and

wherein the mathematical operation is given by;

C(x,y,z)=a C₀(x,y,z)e^bC0(x,y,z)

where C(x,y,z) represents the modified CT number of a 3D CT voxel having a location (x,y,z);

a and b represent weighting coefficients;

and C₀(x,y,z) represents the un-modified CT number, based on a linear attenuation model, of a 3D CT voxel having a location (x,y,z).

29. A system in accordance with claim 27,

wherein the 3D scan data comprise a plurality of CT numbers representing the image intensity of corresponding 3D CT voxels;

and

wherein the DRR generator comprises;

a ray casting subsystem configured to cast a plurality of hypothetical rays through the 3D CT voxels, at the known intensities and from the known positions and angles;

a CT number integrator for integrating along each hypothetical ray the CT numbers corresponding to the CT voxels that are traversed by the hypothetical ray; and

a projector for projecting the integrated values of the CT numbers onto an imaging plane.

30. A system in accordance with claim 29,

wherein the CT number integrator includes:

a bi-linear interpolator configured to perform bi-linear interpolation on the voxels encountered by each ray; and

a polynomial interpolator configured to perform, for each voxel of interest within a voxel slice, a one-dimensional polynomial interpolation over the voxel of interest and over voxels on each adjacent voxel slice.

31. A system in accordance with claim 24, wherein the radiation source comprises an x-ray source, and the 2D radiographic image comprises a 2D x-ray projection image.

32. A system in accordance with claim 24,

wherein the filter operator comprises a top-hat filter configured to select the pixel having the greatest pixel value from each of at least two different neighborhoods within the DRR and the 2D image, and to eliminate the remaining pixels in the neighborhoods.

33. A system in accordance with claim 32,

wherein the filter operator is defined mathematically by:

f_e=f+w×

[f−

γ

_B(f)]−

b×

[φ

_B(f)−

f],

where f_e represents an enhanced image,

f represents an original image,

w and b represent weighting coefficients,

γ

_B(f) represents a structural element for the opening of the original image f, and

φ

_B(f) represents a structural element for the closing of the original image f.

34. A system in accordance with claim 32,

wherein the ROI selector comprises an entropy calculator for calculating an entropy measure H of the DRR.

35. A system in accordance with claim 34,

wherein the entropy measure H comprises a modified Shannon entropy measure given by:

H=−

Σ

I P(I)log P(I),

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

36. A system in accordance with claim 35,

wherein the ROI is selected so that the entropy measure H is maximized within the ROI.

37. A system in accordance with claim 24,

wherein the 2D image and the DRR are each characterized by an array of pixels, each pixel having an associated pixel value that represents the image intensity of a corresponding unit of the anatomical region.

38. A system in accordance with claim 24,

further comprising a similarity measure calculator configured to determine a measure of similarity between the ROI and the 2D image;

and

wherein the motion field generator is configured to estimate the 2D local motion fields by maximizing the similarity measure determined by the similarity measure calculator.

39. A system in accordance with claim 38,

wherein the measure of similarity is based on pattern intensity;

wherein the pixel intensity values at the i-th row and j-th column of each array of pixels for the 2D image and the DRR can be represented as I_Live(i,j) and I_DRR(i,j), respectively; and

wherein the measure of similarity is determined by the similarity measure calculator by comparing pixel intensity values I_Live(i,j) and I_DRR(i,j).

40. A system in accordance with claim 39,

wherein the similarity measure calculator is configured to form a difference image by subtracting a corresponding pixel value of the DRR from each pixel value of the ROI, so that the pixel value at the i-th row and j-th column of the array of pixel values for the difference image is given by:

I_diff(i,j)=I_Live(i,j)−

I_DRR(i,j);

wherein the similarity measure calculator is further configured to apply upon each pixel of the difference image a pattern intensity function, the pattern intensity function being defined by summing asymptotic functions of the gradients of the difference image over the pixels within a neighborhood R; and

wherein the neighborhood R is defined so that the gradients of the difference image can be considered in at least four directions.

41. A system in accordance with claim 40, wherein the pattern intensity function

is characterized by a mathematical formulation given by:

$\underset{i,j}{\sum <br><br>}<br><br>\underset{k·<br><br>k⋐<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 σ

is a weighting constant for weighting the pattern intensity function, and where R is a neighborhood defined around a pixel (I,j) as a center point.

42. A system in accordance with claim 40, wherein the at least four directions

comprise:

a) a substantially horizontal direction;

b) a substantially vertical direction;

c) a substantially diagonal direction at about 45 degrees; and

d) a substantially diagonal direction at about −

45 degrees.

43. A system in accordance with claim 24,

wherein the non-rigid transformation parameters comprise three translational parameters (x, y, and z), and three rotational parameters (r, p, w);

wherein (x, y, and z) represent the three translations of the treatment target along the directions of three mutually orthogonal x-, y-, and z-axes, respectively; and

wherein (r, p, w) represent the three rotations (roll, pitch, yaw) of the treatment target about the three mutually orthogonal x-, y-, and z-axes.

44. A system in accordance with claim 24,

wherein the one or more 2D images comprises x-ray projection images that represent at least two orthogonal projections A and B of the anatomical region onto respective projection image planes.

45. A system in accordance with claim 24,

wherein the parameter determiner is configured to determine a separate set of transformation parameters for each projection x-ray image;

and

wherein the parameter determiner is further configured to combine the resulting parameters for each projection to obtain the 3D transformation parameters.

46. A system in accordance with claim 45,

wherein the 3D transformation parameters are related to the transformation parameters for projections A and B by the following relationship:

x=(x_A+x_B)/2, y=(y_A−

y_B)/√

{square root over (2)}, z=(y_A+y_B)/√

{square root over (2)},

r=θ

_x, p=(θ

_B−

θ

_A)/√

{square root over (2)}, w=(θ

_B+θ

_A)/ √

{square root over (2)}.

47. A system in accordance with claim 24, wherein the motion field generator comprises:

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.

48. A system in accordance with claim 47,

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.

49. A system in accordance with claim 48,

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.

50. A system in accordance with claim 48, 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.

51. A system in accordance with claim 48, 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.

52. A system in accordance with claim 51,

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.

53. A system in accordance with claim 48,

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.

54. A system in accordance with claim 53,

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.

55. A system in accordance with claim 54,

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.

56. A system in accordance with claim 55,

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.

57. A system for registering at least one image of an object with scan data of the object, the object including a target region and a reference structure, the system comprising:

a) means for providing the scan data of the object;

b) a radiation source configured to generate one or more radiographic imaging beams, the imaging beams having known intensities, and having known positions and angles relative to the object;

c) an imaging system configured to generate a radiographic image of the object by detecting the radiographic imaging beam after the beam has traversed at least a portion of the object; and

d) a controller, including;

i) a DRR generator configured to generate at least one DRR (digitally reconstructed radiograph) of the object, using the scan data and the known location, angle, and intensity of the imaging beam;

ii) an ROI selector configured to select an ROI (region of interest) within the DRR, the ROI including the image of the target region and of the reference structure;

iii) an image enhancer configured to enhance the DRR and the ROI by applying a filter operator to the DRR and to the ROI;

iv) a similarity measure determiner configured to determine a similarity measure between the DRR and the ROI;

v) a motion field generator configured to estimate, for each of a plurality of resolution levels, one or more 20 local motion fields within the ROI, using the similarity measure, and to derive a 3D full motion field from the estimated 2D local motion fields; and

vi) a parameter calculator configured to compute from the full motion field a set of non-rigid transformation parameters representing the difference in the position and orientation of the target region as shown in the ROI, as compared to the position and orientation of the target region as shown in the DRR.

58. A system in accordance with claim 57, wherein the object comprises an anatomical region, the target region comprises a treatment target, and the reference structure comprises a skeletal structure.

59. An image registration system for registering at least one 2D image of an anatomical region with previously generated 3D scan data of the anatomical region, the anatomical region including at least one treatment target and at least one reference structure, wherein the 2D image is generated in near real time by detecting one or more radiographic imaging beams after the imaging beams have traversed at least a portion of the anatomical region, the imaging beams having known intensities and known positions and angles relative to the anatomical region, the system comprising:

a) means for providing the 3D scan data of the anatomical region;

b) a DRR generator configured to generate at least one DRR (digitally reconstructed radiograph) of the anatomical region, using the 3D scan data and the known locations, angles, and intensities of the imaging beams;

c) an ROI selector configured to select an ROI (region of interest) within the 2D image, the ROI including the image of the treatment target and of the reference structure;

d) an image enhancer configured to enhance the DRR and the ROI by applying a filter operator to the DRR and to the ROI;

e) a similarity measure calculator configured to determine a measure of similarity between the DRR and the ROI;

f) a motion field generator configured to estimate, for each of a plurality of resolution levels, one or more 2D local motion fields within the ROI, using the similarity measure determined by the similarity measure calculator, and to derive a 3D full motion field from the estimated 2D local motion fields; and

g) a parameter determiner configured to compute from the full motion field a set of non-rigid transformation parameters representing the difference in the position and orientation of the treatment target as shown in the ROI, as compared to the position and orientation of the treatment target as shown in the DRR.