Fiducial-less tracking with non-rigid image registration

• US 7,505,617 B2
• Filed: 12/05/2007
• Issued: 03/17/2009
• 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:
  US11/999,596
• Publication Number:
  US 20080101673A1
• Time in Patent Office: 468 Days
• Field of Search: 382/128, 382/129, 382130-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, 382/100, 382/134, 600/411, 600/130, 600/407, 375/240.16, 378/8.65, 378/21, 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.

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 2D image by processor, wherein the data comprise 3D scan data;

estimating one or more 2D local motion fields within a region shown by the at least one 2D image based on a comparison between the at least one 2D image and the first image,

wherein the first image comprises a 2D x-ray projection image, and wherein the act of estimating the one or more 2D local motion fields within the region comprises maximizing a similarity measure between the at least one 2D image and the 2D x-ray projection image, for each of a plurality of image resolution levels;

constructing a full 3D motion field from the estimated 2D local motion fields; and

determining from the full 3D motion field a set of non-rigid transformation parameters that represent the difference in a 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 at least one 2D image;

wherein the object comprises an anatomical region, and the one or more structures comprise at least one reference structure and at least one target.

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

3. A method in accordance with claim 2,

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

and

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

4. A method in accordance with claim 3, wherein determining a measure of similarity between the at least one 2D image and the x-ray projection image comprises:

a) forming a difference image by subtracting a corresponding pixel value of the at least one 2D image from each pixel value of the x-ray projection 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 gradients of a 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.

5. A method in accordance with claim 4, 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.

6. A method in accordance with claim 4, 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.

7. A method in accordance with claim 1,

wherein the 2D x-ray projection image and the at least one 2D image 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.

8. A method in accordance with claim 7, further comprising:

defining the region by calculating an entropy measure H of the at least one 2D image, and selecting the region so as to maximize the entropy measure H.

9. A method in accordance with claim 8, wherein the entropy measure H is a modified Shannon entropy measure given by:

H=−

Σ

_II 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.

10. A method in accordance with claim 7,

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 at least one 2D image is a synthetic image representing the radiographic image of the target that would be obtained with an x-ray imaging beam at the known intensity, origination position and angle, if the target were positioned in accordance with the 3D scan data.

11. A method in accordance with claim 10, wherein reconstructing the at least one 2D image 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 imaging beam.

12. A method in accordance with claim 11, wherein reconstructing the at least one 2D image 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.

13. A method in accordance with claim 12, wherein 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.

14. A method in accordance with claim 12, wherein 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)=aC₀(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).

15. A method in accordance with claim 1, further comprising:

enhancing the at least one 2D image and the first image so as to increase visibility of the at least one 2D image and the first image, wherein the act of enhancing the at least one 2D image and the first 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.

16. A method in accordance with claim 15, 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.

17. 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.

18. 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.

19. A method in accordance with claim 18, 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.

20. A method in accordance with claim 19, 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.

21. A method in accordance with claim 20, 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)}.