An X-ray imaging system, as shown in
Two-dimensional projection images may be acquired from different gantry angles, and multiple two-dimensional projection images acquired from different gantry angles may be used collectively to build a three-dimensional reconstructed image of the patient anatomy. Accurate reconstruction of the three-dimensional image requires an accurate expression of a geometric relationship between the three-dimensional volume in which the imaging system resides (e.g., a treatment room) and the two-dimensional images acquired at a given gantry angle and FPD position.
Referring to
The actual trajectory of the X-ray source may differ from the above ideal description. For example, the X-ray source may move in and out of the vertical plane by a marginal amount. Moreover, as shown in
The relationship between the X-ray source and the FPD may also vary from the above-described ideal. For example, the relation between the detector assembly and the gantry might not be very rigid. Therefore, at each gantry angle, and due to gravity, the detector assembly may sag a different amount with respect to the central axis (CAX) of the beam (i.e., the line joining the X-ray source and the imaging isocenter), as shown in
However, it has been noted that the out of plane rotations η and σ are quite difficult to determine with reasonable accuracy and these two angles have only a small influence on the image quality compared to other parameters. In practical implementations, these two angles can be kept small (≤10) through good mechanical design and high-accuracy machining. It is therefore reasonable to assume that η=σ=0.
In summary, according to the above model, the plane of the detector is perpendicular to the CAX, but the trajectory of the X-ray source is not a perfect circle and the values of SID f and SAD D in
A projection matrix may be used to describe the geometric relationship between any three-dimensional point (xf,yf,zf) in the imaging room and its projection pixel coordinate (u,v) on the two-dimensional FPD. The projection matrix corresponds to a given position of the source/detector (i.e., X-ray source and FPD) pair. For an ideally circular source trajectory, the projection matrix may be written as:
where the symbol “*” denotes multiplication throughout this document.
However, as noted, the source trajectory is typically not a perfect circle of radius D in the vertical xfzf plane. [T]θ and [R]θ are the unknown translation and rotation, respectively, from the frame of reference of the imaging/treatment room to the frame of reference of the FPD. The subscript θ indicates that these transformations change with the gantry angle θ. Therefore, the projection matrix Pθ for a non-ideal source trajectory may be written as:
Calibration of an imaging system may include determination of the above projection matrix Pθ for the non-ideal source trajectory. Conventionally, this calibration involves imaging of a geometry calibration phantom. Geometry calibration phantoms typically consist of radio-opaque beads at known three-dimensional locations with respect to some known frame of reference. This reference frame may be within the phantom itself. The elements of the projection matrix are found by solving equations which relate the known three-dimensional locations within the phantom with the detected two-dimensional pixel locations in a projection image.
Others calibration methods that do not involve a phantom. For example, calibration may include an iterative reconstruction method of perturbing a model projection matrix and generating forward projections to best-match observed projections. Another method uses redundant projections over 180° and optimization to determine misalignment parameters with respect to an ideal projection matrix.
The existing methods for determining a non-ideal projection matrix are inefficient, time-consuming, based on inaccurate assumptions, and/or incompatible with treatment workflow.
The present inventors have recognized that, for many imaging systems, the assembly which mounts the imaging source 1 to the gantry 3 is often quite rigid, and particularly in a case where the X-ray imaging source 1 is a radiation therapy delivery system (i.e., a linac). As a result, even though the trajectory of the X-ray source 1 does not circumscribe a perfect circle, the actual trajectory may be quite stable and reproducible over a substantial time frame (e.g., several months or years). Moreover, the inventors have recognized that the geometric relationship between the gantry and the detector panel is typically not as rigid (e.g., in order to reduce bulk and/or costs) and therefore varies to a greater extent over time.
The inventors have further discovered a system to decompose the non-ideal projection matrix Pθ into a “flexible” projection sub-matrix and a “rigid” projection sub-matrix. The terms “flexible” and “rigid” are used for convenience and are not intended to denote any particular degree of rigidity or flexibility. Decomposition of the projection matrix Pθ into a first projection sub-matrix PθFlexible and a second projection sub-matrix PθRigid according to some embodiments may be illustrated mathematically as shown below:
According to the embodiments described below, the flexible projection sub-matrix PθFlexible and the rigid projection sub-matrix PθRigid are determined, and then used to convert coordinates of a two-dimensional coordinate system of a detector to a three-dimensional coordinate system, using calibrated projection matrix Pθ=PθFlexible*PθRigid. Thereafter, the flexible projection sub-matrix PθFlexible is periodically re-determined according to a first calibration schedule, and the rigid projection sub-matrix PθRigid is periodically re-determined according to a second calibration schedule, which is less frequent than first calibration schedule.
Calibration of the flexible projection sub-matrix PθFlexible may be faster and less-intrusive than conventional calibration of the non-ideal projection matrix Pθ. According to some embodiments, the flexible projection sub-matrix PθFlexible is re-determined (i.e., calibrated) daily or weekly, and the rigid projection sub-matrix Pa is calibrated every 6-12 months during scheduled maintenance service. The foregoing may result in increased uptime and cost savings.
According to some embodiments, calibration of the flexible projection sub-matrix PθFlexible does not require a phantom, thereby providing further cost savings. As such, calibration of the flexible projection sub-matrix PθFlexible may occur on-the-fly and, if desired, during all scans.
According to some embodiments, most or all of the steps of process 600 are embodied within executable code of a system control program stored in and executed by one or more processing units (e.g., microprocessors, microprocessor cores, execution threads) of control system 10 of
Initially, at S610, a first sub-matrix of a projection matrix is determined. The projection matrix describes a geometrical relationship between a point of a three-dimensional coordinate system and a point of a two-dimensional coordinate system. For example, the projection matrix may describes the geometrical relationship between a point of a three-dimensional coordinate system of an imaging room and a point of a two-dimensional coordinate system of a flat panel detector located in the imaging room.
The first sub-matrix may be associated with in-plane rotation ϕ of the flat panel detector and variation in the principal point (u0,v0) as described above. According to some embodiments, the first sub-matrix represents the first projection sub-matrix PθFlexible described above.
In some embodiments of S610, the first projection sub-matrix PθFlexible is determined as:
Accordingly, determination of the first projection sub-matrix PθFlexible according to some embodiments includes determination of pixel dimensions pw, ph, co-plane rotation ϕ, and principal point (u0,v0). Pixel dimensions pw and ph may be directly obtained from the design specification of the flat panel detector, and several embodiments for determining co-plane rotation ϕ, and principal point (u0,v0) will be described below.
Some embodiments for determining co-plane rotation ϕ, and principal point (u0,v0) and, as a result, the first projection sub-matrix PθFlexible) may be automated, performed on-the-fly and/or performed without use of a phantom. Advantageously, some on-the-fly calibration procedures may be conducted during acquisition of patient images.
A second sub-matrix of the projection matrix is determined at S620. The first sub-matrix determined at S610 and the second sub-matrix determined at S620 comprise a decomposition of the above-described projection matrix. The second sub-matrix is associated with system translations and rotations other than in-plane rotation ϕ of the flat panel detector and variation in the principal point (u0,v0)
According to some embodiments, the second sub-matrix comprises the rigid projection sub-matrix PθRigid described above, and is computed at S620 according to:
where [T]θ represents translation quantity, [R]θ represents rotation quantity, and f represents a source to imager distance (SID).
S620 may comprise determination of the complete non-ideal projection matrix Pθ a using a fiducial phantom or other techniques known in the art, and determination of the second projection sub-matrix PθRigid based on the projection matrix Pθ and the first projection sub-matrix PθFlexible is as follows: PθRigid=[PθFlexible]−1*Pθ. An example of a conventional technique for determining a non-ideal projection matrix Pθ using a fiducial phantom will be provided below.
Next, at S630, a coordinate of the two-dimensional coordinate system is converted to a coordinate of the three-dimensional coordinate system, based on the first sub-matrix and the second sub-matrix. Since the first sub-matrix and the second sub-matrix comprise a decomposition of the projection matrix between the two coordinate systems, a coordinate of the two-dimensional system (e.g., (u, v)) may be converted to a coordinate of the other system (e.g., (xf,yf,zf)) using the transformation:
Operation at S630 may comprise any one or more uses of the projection matrix consisting of the first sub-matrix and the second sub-matrix. For example, S630 may comprise reconstruction of a three-dimensional image based on two-dimensional projection images acquired by the imaging system at various gantry angles, image-assisted patient positioning procedures, radiation field verification, etc. S630 may span several treatment fractions of one or more patients, and any period of time.
S640 and S660 comprise monitoring agents of process 600. Specifically, execution continues at S630 to convert coordinates as desired until it is determined at S640 that the first sub-matrix is to be re-determined. The determination at S640 is based on a first calibration schedule associated with the first sub-matrix. For example, S640 may comprise determining whether a particular amount of time (e.g., one day) has passed since a last determination of the first sub-matrix. The schedule may be event-based. In one example, the first sub-matrix is to be re-determined prior to imaging each patient. The first calibration schedule may comprise a combination of time- and event-based rules.
Flow proceeds to S650 if it is determined at S640 that the first sub-matrix is to be re-determined. The first sub-matrix is re-determined at S650 as described with respect to S610 and flow returns to S630. Accordingly, subsequent coordinate conversion at S630 uses the newly-determined first sub-matrix and the existing second sub-matrix (i.e., Pθ=PθnewFlexible*PθcurrentRigid).
Similarly, execution continues at S630 until it is determined at S660 that the second sub-matrix is to be re-determined. The determination at S650 is based on a second calibration schedule associated with the second sub-matrix. The second calibration schedule may be less frequent than the first calibration schedule according to some embodiments. The second calibration schedule may indicate a re-determination of the second sub-matrix every six months, at each service appointment, or any combination thereof.
If it is determined at S660 to re-determine the second sub-matrix, flow returns to S620 to re-determine the second sub-matrix and then proceeds as described above. In this regard, subsequent coordinate conversions at S630 use the currently-existing first sub-matrix and the newly-determined second sub-matrix (i.e., Pθ=PθcurrentFlexible*PθnewRigid).
Some embodiments for determining the co-plane rotation ϕ, and the principal point (u0,v0) for determination of the first projection sub-matrix at S610 will now be described in more detail. Some of the described embodiments provide automated phantom-less on-the-fly determination of the first projection sub-matrix.
Initially, a MLC of the imaging device is set to a predefined leaf configuration. The predefined configuration may be symmetric about the Y axis. Next, at a single gantry angle, a first image is acquired by the detector while the collimator is at a 0 degree rotational position and a second image is acquired by the detector while the collimator is at a 90 degree rotational position. The two images are superimposed to obtain a rectangular-shaped MLC pattern.
According to some embodiments, only one projection image is acquired, while the collimator is at 0 degrees, to determine the co-plane rotation ϕ, the translation vector {tilde over (S)}, and the principal point (u0,v0).
According to some embodiments, and also illustrated in
As described above, the planned and measured positions of the end (longitudinal) edges of the leaves may be used to identify the Yr axis. However, since the leaves are typically prevented from significantly restricting the open-field, the end positions of these protruded leaves lie near the boundary. While the lateral edges can be detected with better accuracy, the effect of beam penumbra near the boundary may affect the accuracy of detecting the end (i.e., longitudinal) edge of the leaves if those edges are near the boundary. Therefore, in another variation of on-the-fly determination of the first sub-matrix, the MLC may use Y-jaws with a center “notch” such as the notches of
A notch as described above may be simulated by positioning an opposing pair of leaves near the upper Y boundary close together, leaving a narrow opening in the center between the leaves, and similarly positioning an opposing pair of leaves near the lower Y boundary. The Yr axis may then be identified by identifying the openings in the projection image and joining them with a straight line.
The following is a description of a conventional technique for determining the non-ideal projection matrix Pθ. As described above, the non-ideal projection matrix Pθ may be used in some embodiments, in conjunction with the first sub-matrix determined at S610, to determine the second sub-matrix at S620.
Initially, a geometry phantom is placed approximately at the imaging isocenter of the imaging system, with the axis of the phantom being approximately aligned with the expected IEC axis. Images of the phantom are acquired at each of many gantry angles (e.g., at 0.5 or 1 degree intervals). Fiducials are detected in the images and linear equations are formulated which relate three-dimensional locations of the fiducials and the corresponding two-dimensional projections, with the parameters of the projection matrix Pθ as unknowns.
The overdetermined system of linear equations is solved to obtain Pθphantom in the phantom coordinate system, and the location of the imaging source is determined in the phantom coordinate system. Next, a circle is fit through the source points, the center of rotation and plane of rotation are determined, and the isocenter is located at the center of rotation. The IEC Y axis is determined to be normal to the plane of rotation and the source location is projected at gantry 0 degree position on the vertical plane. The IEC Z axis is the unit vector from the center to the projected source location and the IEC X axis is the cross product of the Y and Z axis. Using the above information, Pθphantom is transformed to the projection matrix Pθ.
Avoidance of unnecessary determination of the second sub-matrix may be desirable because this determination (as described above) requires the use of phantoms with embedded fiducials at known three-dimensional locations with respect to a frame of reference. A quality assurance procedure according to some embodiments only requires assurance phantoms at static locations in the imaging room. Moreover, knowledge of the three-dimensional locations of the assurance phantoms is not required. The quality assurance procedure may therefore be performed on-the-fly and with significantly less disruption than re-determination of the second sub-matrix.
A quality assurance procedure according to some embodiments is intended to determine whether the current source trajectory suitably matches the trajectory modeled by the currently-determined second sub-matrix. According to some embodiments, it is first ensured that the first sub-matrix (e.g., PθFlexible) is current for all gantry angles. Next, spherical radio-opaque beads (e.g., ball bearings) are attached to the couch, projections images of the beads are acquired (with or without a patient positioned on the couch), and a three-dimensional image of the beads is reconstructed from the projection images, using the projection matrix composed of the current first sub-matrix and the current second sub-matrix.
Next, the dimensions and/or blur of the ball bearings are determined in three-dimensional space from the three-dimensional image. If the determined dimensions are within a suitable threshold of the actual dimensions, and/or if the detected blur is not significant blur, it is determined that the current second sub-matrix does not require re-determination. Notably, the foregoing procedure may be performed quickly and/or during patient imaging, and may eliminate unnecessary determination of the second sub-matrix (i.e., PθRigid).
An alternative quality assurance procedure attempts to match computed locations of fiducials from different pairs of gantry angles. After ensuring that the first sub-matrix (e.g., PθFlexible) is current for all gantry angles, two-dimensional projections) (u1,v1), . . . , (uN,vN) of spherical radio-opaque beads are acquired from gantry angles θ1, . . . , θN, respectively. The beads are attached to the couch in locations (possibly near the edge of open field) such that their projection is visible even while acquiring portal images or CBCT projection images of the patient.
Assuming the unknown three-dimensional location of such a bead to be (xf,yf,zf), and if the calibration of the rigid component is still holding properly,
is obtained for 1≤i≤N;
wherein in this equation, (ui,vi), Pθi are known but (xf,yf,zf) and λi are unknown.
Next, a known three-dimensional point (xfknown,yfknown,zfknown) is chosen, and its expected two-dimensional projection is computed for the gantry angles under the assumption that the calibration of the rigid component is still valid; and as a result
is obtained, for 1≤i≤N, where in this equation (ui
Corresponding equations are subtracted to obtain:
where (kx,ky,kz) is the unknown constant three-dimensional vector from the unknown bead location (xf,yf,zf) to the known three-dimensional point (xfknown,yfknown,zfknown).
For a given gantry angle θi, there are three equations and four unknowns, λi,kx,ky,kz. However, the three-dimensional unknown vector (kx,ky,kz) is common to all gantry angles. Therefore, under the assumption that the calibration of the rigid component is still holding properly, six equations from two gantry angles (θp, θq) can be used to find five unknowns values (λp, λq,kx,ky,kz).
Thus, the assumption of whether the second sub matrix (i.e., the calibration of rigid components) is sufficiently valid can now be checked. From the projections corresponding to gantry angles θ1, . . . , θN, many pairs of projections, preferably orthogonal to one another, are obtained. (kx,ky,kz) is computed from each such pair. In other words, the unknown location (xf,xf,zf) of each bead can be computed from each such projection pair assuming that the second sub-matrix is still valid. If sufficient agreements in the computed values of the unknown location of the beads are found, it can be concluded that the second sub-matrix is valid.
It should be understood that the embodiments presented above are examples rather than limitations. Those in the art can modify and vary the embodiments without departing from their spirit and scope.
Number | Date | Country | Kind |
---|---|---|---|
201410024798.3 | Jan 2014 | CN | national |
The present application is a continuation of U.S. patent application Ser. No. 15/598,331, filed on May 18, 2017, which a continuation of U.S. patent application Ser. No. 14/598,877, filed on Jan. 16, 2015 and claims priority to Chinese Patent Application No. 201410024798.3, filed on Jan. 20, 2014, and entitled “Calibration Method for X-Ray Imaging System.” The disclosures of each of the above-referenced applications are expressly incorporated herein by reference to their entireties.
Number | Name | Date | Kind |
---|---|---|---|
5442674 | Picard | Aug 1995 | A |
5963612 | Navab | Oct 1999 | A |
5963613 | Navab | Oct 1999 | A |
6044132 | Navab | Mar 2000 | A |
6049582 | Navab | Apr 2000 | A |
6731283 | Navab | May 2004 | B1 |
6753828 | Tuceryan | Jun 2004 | B2 |
20020105484 | Navab | Aug 2002 | A1 |
20020113756 | Tuceryan | Aug 2002 | A1 |
20070172033 | Gorges | Jul 2007 | A1 |
20130121457 | Maltz | May 2013 | A1 |
20130161520 | Jansen | Jun 2013 | A1 |
20150085979 | Zheng | Mar 2015 | A1 |
Number | Date | Country |
---|---|---|
101074869 | Nov 2007 | CN |
101515370 | Aug 2009 | CN |
101582161 | Nov 2009 | CN |
102743184 | Oct 2012 | CN |
103006251 | Apr 2013 | CN |
Entry |
---|
Wein et al. (“Self-calibration of geometric and radiometric parameters for cone-beam computed tomography”, Proc. 11th Int. Meeting on Fully 3D Image Reconstruction in Radio1og-y and Nuclear Medicine, Dec. 31, 2011). |
Rika Baba et al., Comparison of flat-panel detector and image-intensifier detector for cone-beam CT, Computerized Medical Imaging and Graphics, 26:153-158, 2002. |
S. M. Johnston et al., GPU-based iterative reconstruction with total variation minimization for micro-CT, Medical Imaging 2010: Physics of Medical Imaging, 7622(762238): 1-11, 2010. |
D Panetta et al., An optimization-based method for geometrical calibration in cone-beam CT without dedicated phantoms. Physics in Medicine and Biology, 53:3 841-3361, 2008. |
B. P. Selby et al., Geometry calibration for X-ray equipment in radiation treatment devices, retrieved from https://www.researchgate.net/publication/5271860 on Feb. 22, 2016. |
Wolfgang Wein et al., Self-Calibration of Geometric and Radiometric Parameters for Cone-Beam Computed Tomography, Proc. 11th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, 2011. |
Yang Kai et al., A geometric calibration method for cone beam CT systems, Medical Physics, 33(6): 1695-1706, 2006. |
Richard Hartley et al., Camera Models, Multiple View Geometry in Computer Vision, Cambridge University Press, Chapter 6: 153-165, 2004. |
M. J. Daly et al., Geometric calibration of a mobile C-arm for intraoperative cone-beam CT, Medical Physics, 35(5): 2124-2136, 2008. |
The third Office Action in Chinese Application No. 201410024798.3 dated Jul. 10, 2019, 20 pages. |
Number | Date | Country | |
---|---|---|---|
20200022669 A1 | Jan 2020 | US |
Number | Date | Country | |
---|---|---|---|
Parent | 15598331 | May 2017 | US |
Child | 16587630 | US | |
Parent | 14598877 | Jan 2015 | US |
Child | 15598331 | US |