CN112989081A - Method and device for constructing digital reconstruction image library - Google Patents

Abstract

The invention discloses a method and a device for constructing a digital reconstruction image library, which comprise the following steps: acquiring a pair of intraoperative X-ray images and preoperative CT images, wherein the X-ray images are acquired by an intraoperative X-ray machine, and the CT images are reconstructed 3D CT images; acquiring projection positioning part information when an X-ray machine projects a human body to generate an X-ray image, and calculating to obtain an X-ray projection vector field when the X-ray image is generated; sampling unknown projection parameters of the X-ray projection vector field, and combining the X-ray projection vector field with the sampled unknown projection parameters to obtain a series of projection vector fields; simulating the process of generating an X-ray image, converting the CT value in the CT image into an attenuation coefficient, projecting the CT image one by using the series of projection vector fields to obtain a series of DRR images, and realizing the construction of a DRR image library. The invention obviously reduces the DRR image library, accelerates the library building, improves the registration performance of the X-ray and CT images, and is effectively compatible with the existing single rapid DRR technology.

Description

Method and device for constructing digital reconstruction image library

Technical Field

The invention relates to the field of 2D/3D medical image registration, in particular to a method and a device for constructing a digital reconstruction image library in the registration process of an intraoperative X-ray image and a preoperative CT image.

Background

In the operation process of intervention, radiotherapy, endoscope and the like, a doctor performs the operation by means of an X-ray image generated by irradiating a patient with X-rays, however, the X-ray image in the operation is fuzzy and rich in noise, and tissues such as blood vessels and the like cannot be displayed in the X-ray image, so that the accurate operation of the operation is influenced; therefore, it is often necessary to overlay the intraoperative 2D X line image and the preoperatively reconstructed 3D CT image, i.e. align the 2D X line image and the 3D CT image, so as to help the doctor observe the lesion, navigate the surgical instrument, and achieve precise treatment. The registration process is generally as follows: firstly, projecting a Reconstructed 3D CT image to a two-dimensional space by using a projection transformation and Digital Reconstructed Radio (DRR) method, and obtaining a 2D DRR image by one-time projection; traversing all projection transformations of the projection space to obtain a complete digital reconstruction image library (DRR image library for short); then, searching in a DRR image library, and searching for a DRR image which is most matched with the X-ray image; and finally, overlapping the X-ray image and the best matching DRR image to realize the registration of the X-ray image and the CT image. The projective transformations correspond one-to-one to the DRR images, and the best matching DRR images are searched in the DRR image library, namely the projective transformations of the registered X-ray and CT images are searched. It can be seen that the construction of the DRR image library is a crucial link for registering X-ray and CT images.

However, building a DRR image library is very time consuming. Several seconds are required to generate a DRR image using the most common ray tracing DRR method; constructing a complete DRR image library can take hundreds or even thousands of hours, which is difficult to withstand. The projection space is sampled finely, a huge DRR image library is generated, and a DRR image which is best matched with the X-ray image can be searched from the DRR image library; if the projection space sampling frequency is reduced and the DRR image library is reduced, only DRR images which are matched with the X-ray images in a second best mode can be searched from the DRR images, and then the registration accuracy of the X-ray images and the CT images is reduced. To balance registration accuracy and library construction time, the prior art accelerates library construction mainly from two directions: 1. hardware acceleration is used, or a DRR algorithm is improved, so that the single DRR projection process is accelerated, and the generation time of a single DRR image is shortened; such techniques require traversing all of the projective transformations in the projection space without reducing the resulting image library. 2. Dividing projection parameters into in-plane parameters and out-of-plane parameters, sampling and projecting the out-of-plane parameters to generate a partial DRR image library, and performing in-plane deformation processing on each DRR to replace a DRR image generated by sampling and projecting the in-plane parameters so as to obtain a complete image library; the database building time is greatly reduced because only DRR is sampled and projected on the out-of-plane parameters, but the 3D CT area penetrated by X rays in the process of sampling and post-projecting the out-of-plane parameters is different from the 3D CT area penetrated by the X rays in the process of sampling and post-projecting the in-plane parameters, the content of the generated DRR image is also different, the former is processed by in-plane deformation and cannot effectively replace the latter, and the DRR image information is lost and mixed, so that the registration precision is reduced; such techniques require indirect traversal of all projective transformations in the projection space, and the resulting image library is not scaled down.

Disclosure of Invention

We find that the intraoperative X-ray machine and intraoperative X-ray image dicom data acquired by the intraoperative X-ray machine provide projection positioning part information when generating an intraoperative X-ray image, and the information is fully utilized, so that a DRR image library can be obviously reduced, the construction of the DRR image library is greatly accelerated, and the registration accuracy and speed are improved to a certain extent.

We consider the registration process of X-ray and CT images as a registration process of the intraoperative human body and preoperative 3D CT images. Assuming that the intraoperative human body and preoperative 3D CT are completely overlapped in spatial position, if the projection transformation of an X-ray image generated by intraoperative projection human body is known, the projection 3D CT image is transformed, the generated DRR image is the DRR image which is best matched with the X-ray image, and the projection transformation is the projection transformation of registering the X-ray and the CT image; if only the projection positioning part information when generating the X-ray image is known, then combining the known information, traversing the unknown projection transformation in the projection space, searching in the DRR image library constructed in this way, and also searching the DRR image which is best matched with the X-ray image, and generating the projection transformation thereof, namely the projection transformation of the registration X-ray and CT images. In the second case, when constructing the DRR image library, by combining with the projection positioning part information when generating the X-ray image, only the unknown projection transformation in the projection space needs to be traversed, and the constructed DRR image library is significantly reduced compared with traversing all the projection transformations in the projection space, and the library construction time can be greatly shortened.

In view of this, the present invention provides a method for constructing a digital reconstructed image library, including:

acquiring a pair of intraoperative X-ray images and preoperative CT images, wherein the intraoperative X-ray images are acquired by an intraoperative X-ray machine, and the CT images are 3D CT images formed by reconstructing CT slices;

acquiring projection positioning part information when the X-ray machine projects a human body to generate the intraoperative X-ray image; according to the known information, calculating to obtain an X-ray projection vector field when the X-ray machine generates the intraoperative X-ray image in a projection mode, and determining a known part in a projection space;

sampling unknown projection parameters in the unknown parameter space of the X-ray projection vector field, combining the X-ray projection vector field and the sampled unknown projection parameters to obtain a series of projection vector fields, and traversing the unknown projection transformation in the projection space;

simulating the process of generating an X-ray image by projecting the human body by X-rays, converting the CT value in the CT image into an attenuation coefficient, projecting the CT image one by using the series of projection vector fields to obtain a series of DRR images, and realizing the construction of a DRR image library.

Preferably, the X-ray machine is a C-shaped arm X-ray machine.

Based on the above construction method, the present invention further provides a construction apparatus for a digital reconstructed image library, comprising:

a data acquisition unit for: acquiring a pair of intraoperative X-ray images and preoperative CT images, wherein the intraoperative X-ray images are acquired by an intraoperative X-ray machine, and the CT images are 3D CT images formed by reconstructing CT slices;

an X-ray projection vector field calculation unit for: acquiring projection positioning part information when the X-ray machine projects a human body to generate the intraoperative X-ray image, and calculating to obtain an X-ray projection vector field when the X-ray machine projects to generate the intraoperative X-ray image according to the known information, so as to determine a known part in a projection space;

a DRR image library generation unit for: sampling unknown projection parameters in the unknown parameter space of the X-ray projection vector field, combining the X-ray projection vector field and the sampled unknown projection parameters to obtain a series of projection vector fields, and traversing the unknown projection transformation in the projection space; simulating the process of generating an X-ray image by projecting the human body by X-rays, converting the CT value in the CT image into an attenuation coefficient, projecting the CT image one by using the series of projection vector fields to obtain a series of DRR images, and realizing the construction of a DRR image library.

Preferably, the X-ray machine is a C-shaped arm X-ray machine.

The invention fully utilizes the acquired projection positioning part information of the X-ray image in the X-ray machine generation operation, only needs to traverse unknown projection transformation in the projection space, greatly reduces the generated DRR image library, and obviously accelerates the library building time. As mentioned above, the image library generated by the two existing technologies is not reduced, and relatively, the present invention has the following beneficial effects:

(1) the invention obviously reduces the image library, accelerates the registration speed of the X-ray and CT images while accelerating the DRR library construction, and improves the registration precision.

Firstly, the prior art needs to directly or indirectly traverse all projection transformations in the projection space, the number of images in the image library is not reduced, but the invention only traverses unknown projection transformations in the projection space, and the number of images in the image library is greatly reduced, that is, compared with the prior art, the invention greatly reduces the number of DRR images which need to be searched during registration, thereby obviously accelerating the registration speed.

Secondly, the prior art needs to directly or indirectly traverse all projection transformations in the projection space, so that all projection parameters in the projection space have influence on the registration accuracy; in contrast, the registration error is not generated by the projection parameters obtained by using the known information of the projection positioning, and only the unknown projection parameters influence the registration accuracy, so that the registration accuracy can be improved to a certain extent although the DRR image library is reduced.

Moreover, compared with the existing database building technology based on the in-plane and out-plane parameters, the method has the advantage of improving the registration accuracy. The former replaces the DRR image generated by in-plane parameter sampling projection with the DRR image generated by out-of-plane parameter sampling projection through in-plane distortion processing to accelerate the construction of an image library, but as mentioned above, this replacement loss and confusion of DRR image information lead to reduced registration accuracy; the invention combines the known information of projection positioning, reduces the projection space in the unknown parameter space, and then carries out DRR for each projection transformation one by one, and the generated DRR image has no information loss and confusion.

(2) The method can be effectively compatible with the existing single rapid DRR technology, and further accelerates the construction of the DRR image library on the premise of ensuring that the registration performance of the X-ray and CT images is improved to a certain extent.

The invention aims at reducing the DRR image library, the existing single-time rapid DRR technology aims at shortening the generation time of a single DRR image, the two DRR images are complementary and effectively combined with each other, and the construction of the DRR image library can be further accelerated on the premise of ensuring that the registration performance of X-ray and CT images is improved to a certain extent.

In contrast, the existing database building technology based on in-plane and out-of-plane parameters can be combined with the existing single-pass fast DRR technology, but cannot guarantee that the registration performance of the X-ray and CT images is improved. The number of images in an image library finally obtained by the technology is not reduced, that is, the number of DRR images to be searched during registration is not reduced, and the registration speed is not improved. Moreover, as mentioned above, the DRR image generated by the out-of-plane parameter sampling projection is processed by the in-plane deformation, which cannot effectively replace the DRR image generated by the in-plane parameter sampling projection, but the DRR image information is lost and mixed up, resulting in the reduction of the registration accuracy.

Drawings

In order to clearly understand the technical solution of the present invention, the drawings used in describing the embodiments of the present invention will be briefly described. It should be apparent that these drawings are only a part of the present invention, and that other drawings may be obtained by those skilled in the art without inventive effort.

Fig. 1 is a schematic flow chart of a method for constructing a digital reconstructed image library according to embodiment 1 of the present invention;

FIG. 2 is a schematic diagram of a process for calculating an X-ray projection vector field according to embodiment 1 of the present invention;

FIG. 3 is a schematic diagram of information of a projection positioning portion of an X-ray apparatus according to embodiment 1 of the present invention

FIG. 4 is a schematic diagram of a DRR image library generation flow provided in embodiment 1 of the present invention;

FIG. 5 is a schematic view of a series of X-ray projection vector field calculation processes provided in embodiment 1 of the present invention

Fig. 6 is a schematic structural diagram of a device for constructing a digital reconstructed image library according to embodiment 2 of the present invention;

fig. 7 is a graph illustrating the registration of a set of intraoperative X-ray and preoperative CT image pairs using embodiment 3 of the present invention and a conventional DRR method, respectively.

Detailed Description

In order to make the objects, technical means and advantages of the present invention more apparent and complete, embodiments of the present invention will be described below with reference to the accompanying drawings.

Example 1

The embodiment 1 of the present invention provides a method for constructing a digital reconstructed image library, a flow diagram of which is shown in fig. 1, and the method includes the following steps:

s1, acquiring a pair of intraoperative X-ray images and preoperative CT images, wherein the intraoperative X-ray images are acquired by an intraoperative X-ray machine, and the CT images are 3D CT images formed by reconstructing CT slices.

S2, acquiring information of a projection positioning part when the X-ray machine projects a human body to generate the X-ray image in operation, and calculating to obtain an X-ray projection vector field when the X-ray machine projects to generate the X-ray image in operation according to the information, thereby determining a known part in a projection space when the X-ray image is generated.

Fig. 2 shows an implementation process of step S2, which specifically includes the following steps:

s21, acquiring projection positioning part information:

as shown in fig. 3, from the dicom data of the intraoperative X-ray image and the intraoperative X-ray machine specification, some parameters of the intraoperative X-ray image generated by the X-ray machine are read, and these parameters include:

s211, the height between the isocenter and the ground is obtained from an X-ray machine specification;

s212, acquiring the initial value of the height of the treatment couch from the specification of the X-ray machine, and reading the change value of the height of the treatment couch from the dicom data tag of the X-ray image, wherein the sum of the initial value and the change value is the height of the treatment couch from the ground;

s213, acquiring the distance from the point light source S to the isocenter C along the X-ray projection direction from the X-ray machine specification;

s214, the point light source S projects to the central point D of the flat panel detector along the X-ray projection direction_cRead from the X-ray image dicom data tag;

and rotation angles alpha, beta and gamma for driving the X-ray machine to rotate around the isocenter C around three axial directions by the rack read from the X-ray image dicom data tag, and the size of the flat panel detector and the resolution of the X-ray image acquired from the specification of the X-ray machine.

The above is mostly the acquisition mode of the positioning parameters of the large C-shaped arm X-ray machine; for small and medium-sized C-arm X-ray machines, the height of the isocenter from the ground is often adjusted during surgery, and the isocenter is often read from the X-ray image dicom data tag.

S22, determining the coordinate (x) of the isocenter C_c, y_c, z_c) In (1)x _c ：

As shown in fig. 3, the coordinate system uses the vertical ground as the X-axis direction, the head of the human body points to the foot end as the z-axis direction, the y-axis direction is determined by adopting the right-hand rule, and the section of the uppermost end of the human body in the X-ray irradiation area parallel to the treatment couch is taken as the plane yoz, and S215 represents the thickness of the human body in the X-ray irradiation area; coordinate (x) of isocenter C_c, y_c, z_c) In (1),x _c = H _c - (H _t + H _CT) WhereinH _cRepresenting the height of the isocenter C from the ground,H _tin order to ensure that the treatment bed is at the height from the ground,H _CTthe height of the CT slice is shown, namely the thickness of the human body in the X-ray irradiation area; y is_cAnd z_cIs unknownA parameter;

when the X-ray machine projects the human body, the height information of the human body is givenH _c 、H _t 、H _CTThe exact position of the body in the yoz plane is not given, so that only the isocenter coordinate can be determinedx _c And y cannot be determined_cAnd z_c(ii) a Next, in the calculation steps S23 to S26, the parameter y is assumed first_cAnd z_cAre known.

S23, calculating a transformation matrix of the X-ray machine rotating around the isocenter around three axial directions:

firstly, translating the coordinate system, moving the origin to the isocenter C, and calculating a translation matrixT _t (ii) a Then calculating the rotation matrix obtained after the X-ray machine rotates alpha, beta and gamma angles around the X, y and z axes in sequenceR(ii) a Finally, the coordinate system is translated, the original point is moved back to the position of the original coordinate system from the isocenter C, and the translation matrix isT _t Inverse matrix ofT _t ^-1 (ii) a Transformation matrixTThe calculation formula of (a) is as follows:

wherein

，

Namely the x-axis, of the optical system,

to be wound around

Rotating the rotation matrix corresponding to the alpha angle; vector quantity

Is the vector corresponding to the y axis after the coordinate system rotates around the x axis by an angle alpha,

indicating winding

Rotating the rotation matrix corresponding to the angle beta; vector quantity

The vector corresponding to the z-axis is formed by rotating the coordinate system around the x-axis by an angle alpha and then rotating the coordinate system around the y-axis by an angle beta,

indicating winding

And rotating the rotation matrix corresponding to the gamma angle. Rotation matrix

Calculated according to the following formula:

here, the

Representing a vector of windings

A rotation matrix rotated by an angle theta

Of and

the parameters are in one-to-one correspondence, and then each matrix value can be obtained through calculation;

no matter how the sequence of the rotation of the X-ray machine around the three axial directions in actual operation is, the finally obtained rotation matrix is the same as the rotation matrix R calculated according to the method.

S24, calculating coordinates of the point light source of the X-ray machine:

point light source S coordinate (x)_s, y_s, z_s) The calculation formula of (a) is as follows:

wherein the content of the first and second substances,

showing the distance from the point light source S to the isocenter C along the X-ray projection direction, the default point light source S is positioned below the treatment couch at the initial position, and the X-ray

The initial position is perpendicular to the treatment couch.

S25, calculating coordinates of all points of the flat panel detector:

all points D of flat panel detector_mnCoordinates of the object

The calculation method of (2) is as follows:

firstly, calculating the central point D of the flat panel detector according to the formula (6)_cThe coordinates of (a):

wherein

Representing the point light source S to the central point D of the flat panel detector along the X-ray projection direction_cThen all points D of the flat panel detector are calculated according to the formula (7)_mnThe coordinates of (a):

wherein the content of the first and second substances,

indicating point D_mnThe initial position coordinates of the optical sensor (c),

，

，S_L*S_Wis the flat panel detector size, R_L*R_WIs the X-ray image resolution; when in use

，

When the temperature of the water is higher than the set temperature,

from which it is derived

(ii) a Will be provided with

Substituting the formula (7) to obtain all points D of the flat panel detector_mnCoordinates of (2)

。

S26, calculating X-ray projection vector fields of all points projected to the flat panel detector by the point light source of the X-ray machine:

x-ray projection vector field

The calculation formula of (a) is as follows:

wherein

。

Calculation at the above-mentioned S22 to S26Parameter y during the entire process of generating the projection vector field of the X-ray image_cAnd z_cUnknown; parameter(s)

Has been obtained in step S21 as a known parameter; parameter(s)H _CT = R_CT-h*P_h，R_CT-hAnd P_hRespectively, the height resolution and pixel height of the CT slice, are read from the dicom data of the CT slice, and thus the parametersH _CTIs also a known parameter.

Assuming that the parameter y is known_cAnd z_cThrough the steps, the projection vector field when the X-ray image is generated can be accurately calculated, and the projection vector field is used for projecting the 3D CT image, so that the DRR image which is optimally matched with the X-ray image can be obtained; in other words, once the parameter y is determined_cAnd z_cThe registration of the X-ray and CT images can be achieved. However, as previously mentioned, in practice the parameter y_cAnd z_cUnknown, to achieve X-ray to CT image registration, the parameter y needs to be determined_cAnd z_cThe value of (c). Considering that the DRR images correspond to the projection transformation one by one, if the DRR image which is best matched with the X-ray image is found, the corresponding projection transformation is the projection transformation for registering the X-ray and the CT image, namely the projection vector field of the X-ray image is generated by projection; knowing the projection vector field, the parameter y can be derived_cAnd z_cThe value of (c). Therefore, to determine the parameter y_cAnd z_cTo realize the registration of X-ray and CT images, the traversal (y) in the X-ray irradiation area is needed_c, z_c) All possible projection vector fields are obtained through all combinations of the three-dimensional CT image, and the 3D CT image is subjected to projection DRR one by one to obtain a DRR image library; (y) corresponding to DRR image in library that best matches X-ray image_c, z_c) In combination, that is the parameter y_cAnd z_cThe value of (c).

Thus, having determined the known portion in projection space via step S2, the next step is traversal (y)_c, z_c) All possible combinations of (a) to (b), generate a DRR image library.

S3, sampling unknown projection parameters in the unknown parameter space of the X-ray projection vector field, combining the X-ray projection vector field and the sampled unknown projection parameters to obtain a series of projection vector fields, and traversing the unknown projection transformation in the projection space; simulating the process of generating an X-ray image by projecting the human body by X-rays, converting the CT value in the CT image into an attenuation coefficient, projecting the CT image one by using the series of projection vector fields to obtain a series of DRR images, and realizing the construction of a DRR image library.

Fig. 4 shows a specific implementation process of step S3, which includes the following three steps:

s31, sampling unknown projection parameters in the unknown parameter space of the X-ray projection vector field, combining the X-ray projection vector field and the sampled unknown projection parameters to obtain a series of projection vector fields, and traversing the unknown projection transformation in the projection space.

As shown in fig. 5, step S31 further includes the following steps:

s311, specifying the origin of the coordinate system:

a first read voxel point in the CT dicom data is used as a coordinate system origin;

s312, sampling unknown parameters of the projection vector field:

in the projection vector field

Unknown parameter y_cAnd z_cIn the yoz plane, two parameters are sampled at high frequency and equal interval in a full coverage manner to obtain a series (y)_ci, z_cj) A value; wherein

A represents the area covered by the 3D CT in the yoz plane, namely the irradiation area of the X-ray in the yoz plane; i = 0,1, …, I-1, J = 0,1, …, J-1, I and J representing the number of sample points on the y-axis and z-axis, respectively;

s313, obtaining a series of projection vector fields:

sampling value (y)_ci ,z_cj) Are substituted into a formula (8) one by one to obtain a series of projection vector fields

Where i and j correspond to sample values (y)_ci ,z_cj) I = 0,1, …, I-1, J = 0,1, …, J-1; m and n correspond to flat panel detector point D_mn，m = 0,1,…，R_L-1，n = 0,1,…，R_W-1; for example, the sampled value (y)_c0 ,z_c0) Correspondingly obtaining the point light source S of the X-ray machine to all the points D of the flat panel detector_mnProjected vector field of

Value of the sample (y)_c2 ,z_c3) Correspondingly obtaining a point light source S to all points D of the detector_mnProjected vector field of

(ii) a Traverse (y)_ci ,z_cj) All combinations, all projective transformations unknown in the projection space are also traversed.

S32, converting the CT value in the CT image into an X-ray attenuation coefficient:

in order to simulate the process of X-ray projection human body generation operation X-ray image, firstly, the CT value of each voxel in the 3D CT image is converted into the attenuation coefficient of corresponding tissue to X-ray, and the conversion formula is as follows:

wherein k represents that the X-ray passes through the kth tissue of the human body,

represents the attenuation coefficient of the k-th tissue to the X-ray,

the attenuation coefficient of water to X-ray is shown, and HU represents CT value.

S33, projecting the CT images one by one to a flat panel detector by the series of projection vector fields to obtain a series of DRR images, and realizing the construction of a DRR image library:

analog X-rayProcess for projecting an intraoperative X-ray image of the human body, starting from a point source S, with the series of projection vector fields

(I = 0,1, …, I-1, J = 0,1, …, J-1) projecting the 3D CT image one by one onto all points D of the flat panel detector_mn（m = 0,1,…，R_L-1，n = 0,1,…，R_W-1) generating a series of DRR images using ray tracing; point D_mnThe gradation value DRR (m, n) of (d) is calculated as follows:

wherein the content of the first and second substances,I ₀represents the intensity of the X-ray light source, and may be designated as 1; k denotes the X-ray crossing the kth tissue of the body,

represents the attenuation coefficient of the k-th tissue;

represents X-ray

The length of the kth tissue is obtained by calculating the distance between two points of the X-ray penetrating into and out of the kth tissue;

projection vector field

Traversing m and n to obtain corresponding sampling parameters (y)_ci ,z_cj) Projecting the generated DRR image; for example, a projected vector field

Traversing m and n to obtain a sampling parameter (y)_c0 ,z_c0) Projecting the generated DRR image; projection vector field

Traversing m and n to obtain a sampling parameter (y)_c2 ,z_c3) Projecting the generated DRR image; by projecting vector fields in series

And (I = 0,1, …, I-1, J = 0,1, … and J-1) projecting the 3D CT images one by one to obtain a series of DRR images, thereby realizing the construction of a DRR image library.

The above is a specific implementation manner of the method for constructing a digital reconstructed image library provided in embodiment 1 of the present invention. Compared with the prior art, the method has the following advantages: (1) the image library is obviously reduced, the DRR image quantity of registration search of the X-ray and CT images is greatly reduced while the DRR library building is accelerated, and the registration speed can be obviously accelerated; moreover, the registration is not influenced by the known parameters in the projection space and is only influenced by the unknown parameters in the projection space, so that the registration accuracy can be improved to a certain extent; compared with the in-plane and out-of-plane parameter database building technology with information loss and confusion, the method has the advantage of improving the registration accuracy. (2) The invention aims at reducing the DRR image library, the existing single-time rapid DRR technology aims at shortening the generation time of a single DRR image, the two DRR images are complementary and effectively combined with each other, and the construction of the DRR image library can be further accelerated on the premise of ensuring that the registration performance of X-ray and CT images is improved to a certain extent.

Example 2

Based on the method for constructing a digital reconstructed image library provided in embodiment 1, embodiment 2 of the present invention further provides a device for constructing a digital reconstructed image library, a schematic structural diagram of which is shown in fig. 6, and the device includes the following units:

a data acquisition unit 41 for: acquiring a pair of intraoperative X-ray images and preoperative CT images, wherein the intraoperative X-ray images are acquired by an intraoperative X-ray machine, and the CT images are 3D CT images formed by reconstructing CT slices;

an X-ray projection vector field calculation unit 42 for: acquiring projection positioning part information when the X-ray machine projects a human body to generate the intraoperative X-ray image, and calculating to obtain an X-ray projection vector field when the X-ray machine projects to generate the intraoperative X-ray image according to the information, so as to determine a known part in a projection space when the X-ray image is generated;

a DRR image library generating unit 43 for: sampling unknown projection parameters in the unknown parameter space of the X-ray projection vector field, combining the X-ray projection vector field and the sampled unknown projection parameters to obtain a series of projection vector fields, and traversing the unknown projection transformation in the projection space; simulating the process of generating an X-ray image by projecting the human body by X-rays, converting the CT value in the CT image into an attenuation coefficient, projecting the CT image one by using the series of projection vector fields to obtain a series of DRR images, and realizing the construction of a DRR image library.

Compared with the prior art, the digital reconstructed image library construction device provided by the embodiment 2 has the following advantages: (1) the image library is obviously reduced, the DRR image quantity of registration search of the X-ray and CT images is greatly reduced while the DRR library building is accelerated, and the registration speed can be obviously accelerated; moreover, the registration is not influenced by the known parameters in the projection space and is only influenced by the unknown parameters in the projection space, so that the registration accuracy can be improved to a certain extent; compared with the in-plane and out-of-plane parameter database building technology with information loss and confusion, the method has the advantage of improving the registration accuracy. (2) The invention aims at reducing the DRR image library, the existing single-time rapid DRR technology aims at shortening the generation time of a single DRR image, the two DRR images are complementary and effectively combined with each other, and the construction of the DRR image library can be further accelerated on the premise of ensuring that the registration performance of X-ray and CT images is improved to a certain extent.

The X-ray projection vector field calculation unit 42 includes:

an acquire projection information subunit 421 configured to: acquiring the height from an isocenter point to the ground, the distance from a point light source S to an isocenter point C along the X-ray projection direction, the size of a flat panel detector and the resolution of an X-ray image from an X-ray machine specification; reading point light source S from X-ray image dicom data label to flat panel detector central point D along X-ray projection direction_cThe frame drives the X-ray machine to rotate around the isocenter C around three axial rotation angles alpha, beta and gamma; and obtaining the initial value of the height of the treatment bed from the X-ray machine specification, and reading the height change of the treatment bed from the dicom data tag of the X-ray imageThe sum of the two values is the height of the treatment bed from the ground.

The above is mostly the acquisition mode of the positioning parameters of the large C-shaped arm X-ray machine; for small and medium-sized C-arm X-ray machines, the height of the isocenter from the ground is often adjusted during surgery, and the isocenter is often read from the X-ray image dicom data tag.

A calculate X-ray projection vector field subunit 422 for: calculating an X-ray projection vector field when generating the intraoperative X-ray image, specifically including: the system comprises an isocenter point determining subunit 4221, a transformation matrix calculating subunit 4222, a point light source coordinate calculating subunit 4223, a detector coordinate calculating subunit 4224 and a projection vector field calculating subunit 4225.

The determine isocenter subunit 4221 is configured to: determining the isocenter C coordinate (x)_c, y_c, z_c) In (1)x _c ；

The coordinate system takes the vertical ground surface upward as the direction of an X axis, the head pointing to the foot end of the human body as the direction of a z axis, the direction of the y axis is determined by adopting a right hand rule, and the section of the uppermost end of the human body in the X-ray irradiation area parallel to the treatment bed is taken as a plane yoz; coordinate (x) of isocenter C_c, y_c, z_c) In (1),x _c = H _c - (H _t + H _CT) WhereinH _cRepresenting the height of the isocenter C from the ground,H _tin order to ensure that the treatment bed is at the height from the ground,H _CTthe height of the CT slice is shown, namely the thickness of the human body in the X-ray irradiation area; y is_cAnd z_cAre unknown parameters.

The calculation transformation matrix subunit 4222 is configured to: calculating transformation matrix of X-ray machine rotating around three axial directions by using isocenter as centerT；

Firstly, translating the coordinate system, moving the origin to the isocenter C, and calculating a translation matrixT _t (ii) a Then calculating the rotation matrix obtained after the X-ray machine rotates alpha, beta and gamma angles around the X, y and z axes in sequenceR(ii) a Finally, the coordinate system is translated, the original point is moved back to the position of the original coordinate system from the isocenter C, and the translation matrix isT _t Inverse matrix ofT _t ^-1 (ii) a Transformation matrixTThe calculation formula of (a) is as follows:

wherein

，

Namely the x-axis, of the optical system,

to be wound around

Rotating the rotation matrix corresponding to the alpha angle; vector quantity

Is the vector corresponding to the y axis after the coordinate system rotates around the x axis by an angle alpha,

indicating winding

Rotating the rotation matrix corresponding to the angle beta; vector quantity

The vector corresponding to the z-axis is formed by rotating the coordinate system around the x-axis by an angle alpha and then rotating the coordinate system around the y-axis by an angle beta,

indicating winding

And rotating the rotation matrix corresponding to the gamma angle. Rotation matrix

Calculated according to the following formula:

here, the

Representing a vector of windings

A rotation matrix rotated by an angle theta

Of and

the parameters are in one-to-one correspondence, and then each matrix value can be obtained through calculation;

no matter how the sequence of the rotation of the X-ray machine around the three axial directions in actual operation is, the finally obtained rotation matrix is the same as the rotation matrix R calculated according to the method.

The point light source coordinate calculating subunit 4223 is configured to: calculating the coordinate (X) of the point light source S of the X-ray machine_s, y_s, z_s) The calculation formula is as follows:

wherein the content of the first and second substances,

showing the distance from the point light source S to the isocenter C along the X-ray projection direction, the default point light source S is positioned below the treatment couch at the initial position, and the X-ray

The initial position is perpendicular to the treatment couch.

The calculate detector coordinates subunit 4224 for: calculating all points D of flat panel detector_mnCoordinates of (2)

The calculation method is as follows:

firstly, calculating the central point D of the flat panel detector according to the formula (6)_cThe coordinates of (a):

wherein

Representing the point light source S to the central point D of the flat panel detector along the X-ray projection direction_cThen all points D of the flat panel detector are calculated according to the formula (7)_mnThe coordinates of (a):

wherein the content of the first and second substances,

indicating point D_mnThe initial position coordinates of the optical sensor (c),

，

，S_L*S_Wis the flat panel detector size, R_L*R_WIs the X-ray image resolution; when in use

，

When the temperature of the water is higher than the set temperature,

from which it is derived

(ii) a Will be provided with

Substituting the formula (7) to obtain all points D of the flat panel detector_mnCoordinates of (2)

。

A calculate projection vector field subunit 4225, configured to: calculating projection of point light source S of X-ray machine to all points D of flat panel detector_mnX-ray projection vector field of

The calculation formula is as follows:

wherein

。

In the above-mentioned whole process of calculating X-ray projection vector field subunit 422, parameter y_cAnd z_cUnknown; parameter(s)

Already obtained in subunit 421, as a known parameter; parameter(s)H _CT = R_CT-h*P_h，R_CT-hAnd P_hRespectively, the height resolution and pixel height of the CT slice, are read from the dicom data of the CT slice, and thus the parametersH _CTIs also a known parameter.

To determine the parameter y_cAnd z_cTo realize the registration of X-ray and CT images, the traversal (y) in the X-ray irradiation area is needed_c, z_c) All possible projection vector fields are obtained by all combinations of (1), and the 3D CT images are subjected to one-by-one treatmentPerforming projection DRR to obtain a DRR image library; (y) corresponding to DRR image in library that best matches X-ray image_c, z_c) In combination, that is the parameter y_cAnd z_cThe value of (c). Thus, after determining the known portion in projection space, the next element, is a traversal (y)_c, z_c) All possible combinations of (a) to (b), generate a DRR image library.

The DRR image library generating unit 43 includes: a series projection vector field calculation subunit 431, a CT image preprocessing subunit 432, and a generate DRR image library subunit 433.

The series of projection vector field calculation subunit 431, configured to: sampling unknown projection parameters in the unknown parameter space of the X-ray projection vector field, combining the X-ray projection vector field and the sampled unknown projection parameters to obtain a series of projection vector fields, and traversing the unknown projection transformation in the projection space, specifically comprising:

a determine coordinate system origin subunit 4311, configured to: a first read voxel point in the CT dicom data is used as a coordinate system origin;

a sample unknown parameters subunit 4312 to: in the projection vector field

Unknown parameter y_cAnd z_cIn the yoz plane, two parameters are sampled at high frequency and equal interval in a full coverage manner to obtain a series (y)_ci, z_cj) A value; wherein

A represents the area covered by the 3D CT in the yoz plane, namely the irradiation area of the X-ray in the yoz plane; i = 0,1, …, I-1, J = 0,1, …, J-1, I and J representing the number of sample points on the y-axis and z-axis, respectively;

a calculate series of projection vector field subunits 4313 for: sampling value (y)_ci ,z_cj) Are substituted into a formula (8) one by one to obtain a series of projection vector fields

Where i and j correspond to sample values (y)_ci ,z_cj) I = 0,1, …, I-1, J = 0,1, …, J-1; m and n correspond to flat panel detector point D_mn，m = 0,1,…，R_L-1，n = 0,1,…，R_W-1; for example, the sampled value (y)_c0 ,z_c0) Correspondingly obtaining the point light source S of the X-ray machine to all the points D of the flat panel detector_mnProjected vector field of

Value of the sample (y)_c2 ,z_c3) Correspondingly obtaining a point light source S to all points D of the detector_mnProjected vector field of

(ii) a Traverse (y)_ci ,z_cj) All combinations, all projective transformations unknown in the projection space are also traversed.

The CT image preprocessing subunit 432 is configured to: in order to simulate the process of X-ray projection human body generation operation X-ray image, firstly, the CT value of each voxel in the 3D CT image is converted into the attenuation coefficient of corresponding tissue to X-ray, and the conversion formula is as follows:

wherein k represents that the X-ray passes through the kth tissue of the human body,

represents the attenuation coefficient of the k-th tissue to the X-ray,

the attenuation coefficient of water to X-ray is shown, and HU represents CT value.

The generate DRR image library subunit 433 is configured to: simulating the process of projecting X-ray images in human body generation by X-rays, starting from a point light source S, in the series of projection vector fields

(I = 0,1, …, I-1, J = 0,1, …, J-1) projecting the 3D CT image one by one onto all points D of the flat panel detector_mn（m = 0,1,…，R_L-1，n = 0,1,…，R_W-1) generating a series of DRR images using ray tracing; point D_mnThe gradation value DRR (m, n) of (d) is calculated as follows:

wherein the content of the first and second substances,I ₀represents the intensity of the X-ray light source, and may be designated as 1; k denotes the X-ray crossing the kth tissue of the body,

represents the attenuation coefficient of the k-th tissue;

represents X-ray

The length of the kth tissue is obtained by calculating the distance between two points of the X-ray penetrating into and out of the kth tissue;

projection vector field

Traversing m and n to obtain corresponding sampling parameters (y)_ci ,z_cj) Projecting the generated DRR image; for example, a projected vector field

Traversing m and n to obtain a sampling parameter (y)_c0 ,z_c0) Projecting the generated DRR image; projection vector field

Traversing m and n to obtain a sampling parameter (y)_c2 ,z_c3) Projecting the generated DRR image; by projecting vector fields in series

And (I = 0,1, …, I-1, J = 0,1, … and J-1) projecting the 3D CT images one by one to obtain a series of DRR images, thereby realizing the construction of a DRR image library.

Example 3

The method and the device for constructing the digital reconstruction image library provided by the embodiments 1 and 2 are adopted to construct a DRR image library, mutual information is taken as similarity measurement, Powell is taken as optimization search, and a group of preoperative CT angiography and intraoperative X-ray image pairs of a patient with thoraco-aortic endoluminal repair are registered in the embodiment 3 of the invention; constructing a DRR image library by adopting a traditional 6-degree-of-freedom ray tracing DRR method, taking mutual information as similarity measurement, taking Powell as optimized search, and registering the CT angiography and X-ray image pair; we compared the experimental results of the two methods. The only difference between the two registration methods is that the construction methods of the DRR image library are different: embodiment 3 employs the digital reconstructed image library construction methods and apparatus provided in embodiments 1 and 2, while the latter employs the conventional 6-degree-of-freedom ray tracing DRR image library construction method.

We performed experiments on a computer configured as intel i9-9820X CPU and 112G memory in matlab environment. The two DRR image library construction methods have sampling intervals of 3mm for translation and 1 degree for angle. The CT angiograms to be registered are aligned with the X-ray images, the preoperative CT slice resolution is 512X 512, the pixel size is 0.6836mm X0.6836 mm,H _CTis 350 mm; the intraoperative X-ray image was collected by GE Innova4100-IQ C arm, and the projection part information obtained from the specification and dicom data when the intraoperative X-ray image was generated by the X-ray machine is shown in table 1.

TABLE 1 projection part information when X-ray machine generates intraoperative X-ray image

TABLE 2 comparison of the results

Table 2 compares the results of the two registration methods quantitatively, and fig. 7 compares the registration effect of the two methods visually. Fig. 7 shows that the upper, middle and lower three sub-images respectively correspond to an intraoperative X-ray image to be registered, an overlay effect image of the DRR and the X-ray image after the image pair is registered by the conventional DRR method (the overlaid DRR and X-ray images are alternately displayed in a checkerboard), and an overlay effect image of the DRR and the X-ray image after the image pair is registered in embodiment 3.

As can be seen from table 2, when the CT angiography and X-ray image pair are registered, compared with the conventional 6-degree-of-freedom ray tracing DRR method, embodiment 3 employs the digital reconstruction image library construction method and apparatus provided by the present invention, so that (1) the number of images in the DRR image library is greatly reduced, and the library construction time is suddenly reduced from thousands of hours to less than three hours; (2) the registration speed is increased by more than three times; (3) the similarity between the DRR image and the X-ray image after registration is improved to a certain degree. Comparing the positions indicated by arrows in the figure 7 and the following figure with each other by visual inspection, especially observing the matching degree between the rib edge in the DRR image and the rib edge in the X-ray image, it can also be seen that, compared with the conventional DRR method, the embodiment 3 can improve the matching degree of the rib position, that is, improve the registration accuracy, by using the method and the apparatus for constructing the digital reconstructed image library provided by the present invention.

The foregoing is illustrative of the preferred embodiments of the present invention and is not to be construed as limiting thereof in any way. Although the present invention has been described with reference to the preferred embodiments, it is not intended to be limited thereto. Those skilled in the art can make numerous possible variations and modifications to the present teachings, or modify equivalent embodiments to equivalent variations, without departing from the scope of the present teachings, using the methods and techniques disclosed above. Therefore, any simple modification, equivalent change and modification made to the above embodiments according to the technical essence of the present invention are still within the scope of the protection of the technical solution of the present invention, unless the contents of the technical solution of the present invention are departed.

Claims (8)

1. A method for constructing a digital reconstructed image library, comprising:

acquiring a pair of intraoperative X-ray images and preoperative CT images, wherein the intraoperative X-ray images are acquired by an intraoperative X-ray machine, and the CT images are 3D CT images formed by reconstructing CT slices;

acquiring projection positioning part information when the X-ray machine projects a human body to generate the intraoperative X-ray image; according to the known information, calculating to obtain an X-ray projection vector field when the X-ray machine generates the intraoperative X-ray image in a projection mode, and determining a known part in a projection space;

sampling unknown projection parameters in the unknown parameter space of the X-ray projection vector field, combining the X-ray projection vector field and the sampled unknown projection parameters to obtain a series of projection vector fields, and traversing the unknown projection transformation in the projection space;

simulating the process of generating an X-ray image by projecting the human body by X-rays, converting the CT value in the CT image into an attenuation coefficient, projecting the CT image one by using the series of projection vector fields to obtain a series of digital reconstruction images, namely a series of DRR images, and realizing the construction of a DRR image library.

2. The method of claim 1, wherein the X-ray machine is a C-arm X-ray machine.

3. The method according to claim 2, characterized in that, the projection positioning part information when the X-ray machine projects the human body to generate the X-ray image in the operation is obtained; according to the known information, calculating an X-ray projection vector field when the X-ray machine generates the intraoperative X-ray image in a projection manner, thereby determining a known part in a projection space, specifically comprising:

step A, acquiring projection positioning part information;

reading partial parameters of the X-ray machine for generating the intraoperative X-ray image from the intraoperative X-ray machine specification and dicom data of the intraoperative X-ray image, wherein the parameters comprise: the height of the isocenter point from the ground, the height of the treatment bed from the ground, the distance from the point light source to the isocenter along the X-ray projection direction, the distance from the point light source to the center point of the flat panel detector along the X-ray projection direction, the rotation angle of the X-ray machine driven by the rack to rotate around the isocenter around three axial directions, the size of the flat panel detector and the resolution of an X-ray image;

step B, determining the coordinate (x) of the isocenter C_c, y_c, z_c) In (1)x _c ；

The coordinate system takes the vertical ground surface upward as the direction of an X axis, the head pointing to the foot end of the human body as the direction of a z axis, the direction of the y axis is determined by adopting a right hand rule, the section of the uppermost end of the human body in the X-ray irradiation area parallel to the treatment bed is taken as a plane yoz, and S215 represents the thickness of the human body in the X-ray irradiation area; coordinate (x) of isocenter C_c, y_c, z_c) In (1),x _c = H _c- (H _t+ H _CT) WhereinH _cRepresenting the height of the isocenter C from the ground,H _tin order to ensure that the treatment bed is at the height from the ground,H _CTthe height of the CT slice is shown, namely the thickness of the human body in the X-ray irradiation area; y is_cAnd z_cIs an unknown parameter;

step C, calculating a transformation matrix of the X-ray machine rotating around the isocenter around three axial directionsT；

Firstly, translating the coordinate system, moving the origin to the isocenter C, and calculating a translation matrixT _t (ii) a Then calculating the rotation matrix obtained after the X-ray machine rotates alpha, beta and gamma angles around the X, y and z axes in sequenceR(ii) a Finally, the coordinate system is translated, the original point is moved back to the position of the original coordinate system from the isocenter C, and the translation matrix isT _t Inverse matrix ofT _t ^-1 (ii) a Transformation matrixTThe calculation formula of (a) is as follows:

wherein

，

Namely the x-axis, of the optical system,

to be wound around

Rotating the rotation matrix corresponding to the alpha angle; vector quantity

Is the vector corresponding to the y axis after the coordinate system rotates around the x axis by an angle alpha,

indicating winding

Rotating the rotation matrix corresponding to the angle beta; vector quantity

The vector corresponding to the z-axis is formed by rotating the coordinate system around the x-axis by an angle alpha and then rotating the coordinate system around the y-axis by an angle beta,

indicating winding

Rotating a rotation matrix corresponding to the gamma angle; rotation matrix

Calculated according to the following formula:

here, the

Representing a vector of windings

A rotation matrix rotated by an angle theta

Of and

the parameters are in one-to-one correspondence, and then each matrix value can be obtained through calculation;

no matter how the sequence of the rotation of the X-ray machine around the three axial directions in actual operation is, the finally obtained rotation matrix is the same as the rotation matrix R calculated according to the method;

step D, calculating the coordinate (X) of the point light source S of the X-ray machine_s, y_s, z_s) The calculation formula is as follows:

wherein the content of the first and second substances,

showing the distance from the point light source S to the isocenter C along the X-ray projection direction, the default point light source S is positioned below the treatment couch at the initial position, and the X-ray

The initial position is vertical to the treatment bed;

e, calculating all points D of the flat panel detector_mnCoordinates of (2)

The calculation method is as follows:

firstly, calculating the central point D of the flat panel detector according to the formula (6)_cThe coordinates of (a):

wherein

Representing the point light source S to the central point D of the flat panel detector along the X-ray projection direction_cThen all points D of the flat panel detector are calculated according to the formula (7)_mnThe coordinates of (a):

wherein the content of the first and second substances,

indicating point D_mnThe initial position coordinates of the optical sensor (c),

，

，S_L*S_Wis the flat panel detector size, R_L*R_WIs the X-ray image resolution; when in use

，

When the temperature of the water is higher than the set temperature,

from which it is derived

(ii) a Will be provided with

Substituting the formula (7) to obtain all points D of the flat panel detector_mnCoordinates of (2)

;

Step F, calculating projection of point light source S of X-ray machine to all points D of flat panel detector_mnX-ray projection vector field of

The calculation formula is as follows:

wherein

；

In the whole process of generating X-ray image projection vector field through calculation of steps B to F, the parameter y_cAnd z_cUnknown, parameter

Are known parameters.

4. A method according to claim 3, characterized by sampling unknown projection parameters in a parameter space in which the X-ray projection vector field is unknown, combining the X-ray projection vector field with the sampled unknown projection parameters resulting in a series of projection vector fields, thereby traversing the unknown projection transformation in projection space; simulating the process of generating an X-ray image by projecting a human body by X-rays, converting the CT value in the CT image into an attenuation coefficient, projecting the CT image one by using the series of projection vector fields to obtain a series of DRR images, and realizing the construction of a DRR image library, wherein the process specifically comprises the following steps:

step G, sampling unknown projection parameters in the unknown parameter space of the X-ray projection vector field, combining the X-ray projection vector field and the sampled unknown projection parameters to obtain a series of projection vector fields, and traversing the unknown projection transformation in the projection space; the method specifically comprises the following steps:

g1, designating the origin of a coordinate system;

a first read voxel point in the CT dicom data is used as a coordinate system origin;

g2, sampling unknown parameters of the projection vector field;

in the projection vector field

Unknown parameter y_cAnd z_cIn the yoz plane, two parameters are sampled at high frequency and equal interval in a full coverage manner to obtain a series (y)_ci, z_cj) A value; wherein

A represents the area covered by the 3D CT in the yoz plane, namely the irradiation area of the X-ray in the yoz plane; i = 0,1, …, I-1, J = 0,1, …, J-1, I and J representing the number of sample points on the y-axis and z-axis, respectively;

g3, obtaining a series of projection vector fields;

sampling value (y)_ci ,z_cj) Are substituted into a formula (8) one by one to obtain a series of projection vector fields

Where i and j correspond to sample values (y)_ci ,z_cj) I = 0,1, …, I-1, J = 0,1, …, J-1; m and n correspond to flat panel detector point D_mn，m = 0,1,…，R_L-1，n = 0,1,…，R_W-1; traverse (y)_ci ,z_cj) All in combination, i.eTraversing all unknown projection transformations in the projection space;

step H, converting the CT value in the CT image into an X-ray attenuation coefficient:

in order to simulate the process of X-ray projection human body generation operation X-ray image, firstly, the CT value of each voxel in the 3D CT image is converted into the attenuation coefficient of corresponding tissue to X-ray, and the conversion formula is as follows:

wherein k represents that the X-ray passes through the kth tissue of the human body,

represents the attenuation coefficient of the k-th tissue to the X-ray,

represents the attenuation coefficient of water to X-ray, HU represents CT value;

step I, projecting the CT images to a flat panel detector one by one according to the series of projection vector fields to obtain a series of DRR images and realize the construction of a DRR image library;

simulating the process of projecting X-ray images in human body generation by X-rays, starting from a point light source S, in the series of projection vector fields

(I = 0,1, …, I-1, J = 0,1, …, J-1) projecting the 3D CT image one by one onto all points D of the flat panel detector_mn（m = 0,1,…，R_L-1，n = 0,1,…，R_W-1) generating a series of DRR images using ray tracing; point D_mnThe gradation value DRR (m, n) of (d) is calculated as follows:

wherein the content of the first and second substances,I ₀represents the intensity of the X-ray light source, and may be designated as 1; k denotes the X-ray crossing the kth tissue of the body,

represents the attenuation coefficient of the k-th tissue;

represents X-ray

The length of the kth tissue is obtained by calculating the distance between two points of the X-ray penetrating into and out of the kth tissue;

projection vector field

Traversing m and n to obtain corresponding sampling parameters (y)_ci ,z_cj) Projecting the generated DRR image; by projecting vector fields in series

And (I = 0,1, …, I-1, J = 0,1, … and J-1) projecting the 3D CT images one by one to obtain a series of DRR images, thereby realizing the construction of a DRR image library.

5. An apparatus for constructing a digital reconstructed image library, comprising:

a data acquisition unit for: acquiring a pair of intraoperative X-ray images and preoperative CT images, wherein the intraoperative X-ray images are acquired by an intraoperative X-ray machine, and the CT images are 3D CT images formed by reconstructing CT slices;

an X-ray projection vector field calculation unit for: acquiring projection positioning part information when the X-ray machine projects a human body to generate the intraoperative X-ray image, and calculating to obtain an X-ray projection vector field when the X-ray machine projects to generate the intraoperative X-ray image according to the information, so as to determine a known part in a projection space;

a DRR image library generation unit for: sampling unknown projection parameters in the unknown parameter space of the X-ray projection vector field, combining the X-ray projection vector field and the sampled unknown projection parameters to obtain a series of projection vector fields, and traversing the unknown projection transformation in the projection space; simulating the process of generating an X-ray image by projecting the human body by X-rays, converting the CT value in the CT image into an attenuation coefficient, projecting the CT image one by using the series of projection vector fields to obtain a series of DRR images, and realizing the construction of a DRR image library.

6. The apparatus of claim 5, wherein the X-ray machine is a C-arm X-ray machine.

7. The apparatus of claim 6, wherein the X-ray projection vector field computing unit comprises: a projection information acquiring subunit and an X-ray projection vector field calculating subunit, wherein:

the sub-unit for obtaining projection information is configured to: reading partial parameters of the X-ray machine for generating the intraoperative X-ray image from the intraoperative X-ray machine specification and dicom data of the intraoperative X-ray image, wherein the parameters comprise: the height of the isocenter point from the ground, the height of the treatment bed from the ground, the distance from the point light source to the isocenter along the X-ray projection direction, the distance from the point light source to the center point of the flat panel detector along the X-ray projection direction, the rotation angle of the X-ray machine driven by the rack to rotate around the isocenter around three axial directions, the size of the flat panel detector and the resolution of an X-ray image;

the calculating X-ray projection vector field subunit is used for: calculating an X-ray projection vector field when generating the intraoperative X-ray image, thereby determining a known portion in a projection space, specifically comprising: determine isocenter point subunit, calculate transform matrix subunit, calculate pointolite coordinate subunit, calculate detector coordinate subunit, and calculate projection vector field subunit, wherein:

the determine isocenter subunit is to: determining the isocenter C coordinate (x)_c, y_c, z_c) In (1)x _c ；

The coordinate system takes the vertical ground surface as the direction of the x axis and the head of the human body as the pointing direction of the foot end asDetermining the direction of the y axis by adopting a right hand rule in the direction of the z axis, and enabling a tangent plane of the uppermost end of the human body in the X-ray irradiation area, which is parallel to the treatment bed, to be a plane yoz; coordinate (x) of isocenter C_c, y_c, z_c) In (1),x _c = H _c- (H _t+ H _CT) WhereinH _cRepresenting the height of the isocenter C from the ground,H _tin order to ensure that the treatment bed is at the height from the ground,H _CTthe height of the CT slice is shown, namely the thickness of the human body in the X-ray irradiation area; y is_cAnd z_cIs an unknown parameter;

the computation transformation matrix subunit is configured to: calculating transformation matrix of X-ray machine rotating around three axial directions by using isocenter as centerT；

Firstly, translating the coordinate system, moving the origin to the isocenter C, and calculating a translation matrixT _t (ii) a Then calculating the rotation matrix obtained after the X-ray machine rotates alpha, beta and gamma angles around the X, y and z axes in sequenceR(ii) a Finally, the coordinate system is translated, the original point is moved back to the position of the original coordinate system from the isocenter C, and the translation matrix isT _t Inverse matrix ofT _t ^-1 (ii) a Transformation matrixTThe calculation formula of (a) is as follows:

wherein

，

Namely the x-axis, of the optical system,

to be wound around

Rotating the rotation matrix corresponding to the alpha angle; vector quantity

Is the vector corresponding to the y axis after the coordinate system rotates around the x axis by an angle alpha,

indicating winding

Rotating the rotation matrix corresponding to the angle beta; vector quantity

The vector corresponding to the z-axis is formed by rotating the coordinate system around the x-axis by an angle alpha and then rotating the coordinate system around the y-axis by an angle beta,

indicating winding

Rotating a rotation matrix corresponding to the gamma angle; rotation matrix

Calculated according to the following formula:

here, the

Representing a vector of windings

A rotation matrix rotated by an angle theta

Of and

the parameters are in one-to-one correspondence, and then each matrix value can be obtained through calculation;

no matter how the sequence of the rotation of the X-ray machine around the three axial directions in actual operation is, the finally obtained rotation matrix is the same as the rotation matrix R calculated according to the method;

the point light source coordinate calculating subunit is configured to: calculating the coordinate (X) of the point light source S of the X-ray machine_s, y_s, z_s) The calculation formula is as follows:

wherein the content of the first and second substances,

showing the distance from the point light source S to the isocenter C along the X-ray projection direction, the default point light source S is positioned below the treatment couch at the initial position, and the X-ray

The initial position is vertical to the treatment bed;

the calculate detector coordinates subunit to: calculating all points D of flat panel detector_mnCoordinates of (2)

The calculation method is as follows:

firstly, calculating the central point D of the flat panel detector according to the formula (6)_cThe coordinates of (a):

wherein

Representing the point light source S to the central point D of the flat panel detector along the X-ray projection direction_cThen all points D of the flat panel detector are calculated according to the formula (7)_mnThe coordinates of (a):

wherein the content of the first and second substances,

indicating point D_mnThe initial position coordinates of the optical sensor (c),

，

，S_L*S_Wis the flat panel detector size, R_L*R_WIs the X-ray image resolution; when in use

，

When the temperature of the water is higher than the set temperature,

from which it is derived

(ii) a Will be provided with

Substituting the formula (7) to obtain all points D of the flat panel detector_mnCoordinates of (2)

;

The projection vector field calculating subunit is configured to: calculating projection of point light source S of X-ray machine to all points D of flat panel detector_mnX-ray projection vector field of

The calculation formula is as follows:

wherein

；

In the whole process of calculating the X-ray projection vector field sub-unit_cAnd z_cFor unknown parameters, parameters

Are known parameters.

8. The apparatus of claim 7, wherein the DRR image library generating unit comprises: a series projection vector field calculation subunit, a CT image preprocessing subunit, and a generate DRR image library subunit, wherein:

the series of projection vector field calculation subunits are configured to: sampling unknown projection parameters in the unknown parameter space of the X-ray projection vector field, combining the X-ray projection vector field and the sampled unknown projection parameters to obtain a series of projection vector fields, and traversing the unknown projection transformation in the projection space, specifically comprising:

determining a coordinate system origin subunit for: a first read voxel point in the CT dicom data is used as a coordinate system origin;

a sample unknown parameter subunit for: in the projection vector field

Unknown parameter y_cAnd z_cIn the yoz plane, two parameters are sampled at high frequency and equal interval in a full coverage manner to obtain a series (y)_ci, z_cj) A value; wherein

A represents the area covered by the 3D CT in the yoz plane, namely the irradiation area of the X-ray in the yoz plane; i = 0,1, …, I-1, J = 0,1, …, J-1, I and J representing the number of sample points on the y-axis and z-axis, respectively;

a compute series projection vector field subunit to: sampling value (y)_ci ,z_cj) Are substituted into a formula (8) one by one to obtain a series of projection vector fields

Where i and j correspond to sample values (y)_ci ,z_cj) I = 0,1, …, I-1, J = 0,1, …, J-1; m and n correspond to flat panel detector point D_mn，m = 0,1,…，R_L-1，n = 0,1,…，R_W-1; traverse (y)_ci ,z_cj) All combinations, namely all unknown projection transformations in the projection space are traversed;

the CT image preprocessing subunit is used for: in order to simulate the process of X-ray projection human body generation operation X-ray image, firstly, the CT value of each voxel in the 3D CT image is converted into the attenuation coefficient of corresponding tissue to X-ray, and the conversion formula is as follows:

wherein k represents that the X-ray passes through the kth tissue of the human body,

represents the attenuation coefficient of the k-th tissue to the X-ray,

represents the attenuation coefficient of water to X-ray, HU represents CT value;

the generate DRR image library subunit is configured to: simulating the process of projecting X-ray images in human body generation by X-rays, starting from a point light source S, in the series of projection vector fields

(I = 0,1, …, I-1, J = 0,1, …, J-1) projecting the 3D CT image one by one onto all points D of the flat panel detector_mn（m = 0,1,…，R_L-1，n = 0,1,…，R_W-1) generating a series of DRR images using ray tracing; point D_mnThe gradation value DRR (m, n) of (d) is calculated as follows:

wherein the content of the first and second substances,I ₀represents the intensity of the X-ray light source, and may be designated as 1; k denotes the X-ray crossing the kth tissue of the body,

represents the attenuation coefficient of the k-th tissue;

represents X-ray

The length of the kth tissue is obtained by calculating the distance between two points of the X-ray penetrating into and out of the kth tissue;

projection vector field

Traversing m and n to obtain corresponding sampling parameters (y)_ci ,z_cj) Projecting the generated DRR image; by projecting vector fields in series

And (I = 0,1, …, I-1, J = 0,1, … and J-1) projecting the 3D CT images one by one to obtain a series of DRR images, thereby realizing the construction of a DRR image library.

Images