Cone-beam CT functional imaging method by using volume integral model

In order to improve the precision of functional imaging of cone beam computed tomography (CBCT), this paper firstly uses the dynamic contrast enhancement tomography (DCE-CT) of the white rabbit as the measured object and establishes volume integral model to obtain the projection data. Then the optimization method is used to solve the optimal parameter pairs of the voxel time density curve (TDC). Finally, the results of the perfusion are obtained by the deconvolution method. The results show that the TDC correlation coefficient is 83.99% after, and the maximum of Spearman correlation coefficient of the perfusion parameter is 0.5125, and the projection time consumption is 7.633 seconds through the volume integral model. It can be seen that the volume integral model is closer to the real projection and it can obtain more accurate perfusion data.


INTRODUCTION
The essence of tumor radiotherapy is to make stronger radioactivity parity radiation, such as X and Y ray, focus on one point, and utilize the energy contained in the focus to destroy the tumor cell structure [1].Radiotherapy techniques include anatomical imaging guided therapy and functional imaging guided therapy.Anatomical imaging guidance mainly uses the form of tumor to identify malignant tumor and benign tumor.The image generated by anatomical image guidance technology has good resolution, but it cannot reflect the physiological information of the tissues accurately, which is not conducive to early cancer detection.The functional imaging guidance technology that uses the information of tumor physiology can make up for this deficiency [2].The principle of CT perfusion imaging is that firstly injecting a proper dose of contrast agent in the vein, then continuously dynamic scanning on the selected level and obtaining the curve of each voxel with time in the region of interest, finally using different mathematical models to find out the functional parameters of the organization [3].
CBCT system is an important equipment for radiation therapy.The ray source and the area detector are synchronized 360 rotation scanning around the measured object.It has the characteristics of small radiation and high imaging precision [4].
In the study of CBCT perfusion imaging, paper [5] proved the possibility of CBCT functional imaging by projection data.Paper [6] used the deconvolution algorithm to establish the physiological model of blood supply, and solved the perfusion parameters in the perfusion imaging.Paper [7][8][9] designed different projection models to obtain projection data and gained the better perfusion results.
But these projection models were based on the two-dimensional plane, which is not consistent with the actual scanning process of CBCT.Therefore, on the basis of previous research, this paper will use the volume integral model which is closer to the CBCT scanning process to obtain the projection data and complete the CBCT functional imaging.

CBCT SAMPLING PROBLEM AND PERFUSION SIMULATION PROCESS
The scanning time for CBCT is usually tens of seconds, and it is impossible to obtain the density curve of each voxel at every moment.But according to the kinetic theory, there is no significant change in hemodynamic parameters in 2-3min.Therefore, we can use approximate time density curve to complete perfusion imaging.In the CT imaging system, the reconstruction algorithm plays a vital role [10].According to the discreteness of voxel and the coverage of ray beam and voxel, the ray coverage model can be divided into point integral model, line integral model, area integral model and volume integral model.

Obtain the simulated projection data
In the volume integral model, the detector is arranged by a square probe.A ray beam is a four pyramid consisting of four vertices of a ray source and a probe.The volume weight of each voxel can be expressed: wi,j=Vi,j/a 3   (1) Vi,j is the intersecting volume value of i ray beam and j voxel, a 3 is the volume value of voxel, Vi,j is an irregular convex polyhedron, which can be solved by the coordinates of the vertex of the space.

Simplified projection model
Taking the center plane of the measured object as the XOY plane, the radiation source at the rotation angle θ=0 is the origin, and the connection between the ray source and the center of the object is the Y axis, and the space rectangular coordinate system for the central slice is established(Fig.1).
Taking the center of the object as the center of rotation and rotating around one circle, we can see that the three-dimensional figure has two characteristics: first, the figure is symmetrical about the plane; second, every 90 degrees per rotation, the intersecting situation is the same.So the model can be simplified, only to consider the part y≥0 the rotation angle θ∈[0,90].

Discuss the coverage of slice and ray beam on XOY plane
The coverage of the ray and voxels on the plane depends on the coverage of two edge rays and voxels.The coverage of single ray and voxel can be divided into full coverage (dark part in Fig. 2), partial coverage (light part in Fig. 2) and no cover (white part in Fig. 2).Then the voxel index of the full covered is obtained by using the value of the part covered voxel index set and the value of the ray slope k.Set i=1,enter the cycle condition, when i≤m, all the voxels of m row are traversed, and the first voxel (i,j) is taken if the line has the voxel passing through the ray.(i,1), (i,2)… (i,j-1) will be stored in the index set {index_full}.Then set i=2, continue the loop, when i>mend the loop.If there is no voxel which is crossed by ray beam in this row, first solved the minimum number of row imin,then set i-1, enter the loop condition.When imin<i≤m and k>0 or i≤ imin ≤m and k≤0, all the voxel index will be stored in the index set {index_full},then set i=2, continue loop, when i>m, end loop.

Calculate voxel coordinate set
The coordinates of fully covered voxel are the coordinates of four vertices.And the coordinates of partial covered voxel need to be discussed.Fig. 3 shows a diagram of the relation between a single voxel and a ray line in a partial covered case.p1, p2, p3, p4 are the four vertices of the voxel in the XOY plane.Voxel (m,n) has two intersection points (x1,y1), (x2,y2).Set δx=|x1-x2|, δy=|y1-y2|.When δx=a and δy=a , the ray line coincides the diagonal line (as in Fig. 3 (a)).
If k>a, record p3,else p1; When δx=a and δy=a, the ray line cross over the voxel from X direction (as in Fig. 3 (b)).If k>0, record p3 and p4,else p1 and p2; When δx≠a and δy=a, the ray line cross over the voxel from Y direction (as in Fig. 3 (c)).If k>0, record p1and p3 ,else p2 and p4; When δx≠a and δy≠a, the ray line cross over the adjacent edges of the voxel (as in Fig. 3 (d)), it needs to discuss voxel (m-1,n).When (m-1,n) is crossed over by ray line, if k>0, record p3,else p1;When (m-1,n) does not exist or is not crossed over by ray line, if k>0,record p1, p3, p4 else p2, p3, p4.

Calculate coverage of ray beam and voxel
Calculating the coverage of the ray beam and the voxel by the set operation.If a vertex set of l1 is A , and a vertex set of l2 is B, the result is C=(A+B)-(A-B).

Discuss the coverage of slice and ray beam in three-dimensional space
In three-dimensional space, the coverage can be divided into three types: no cutting, partial cutting and complete cutting.When no cutting (as shown in Fig. 4 (a)) , the height of the upper slope of the beam is higher than the upper surface of the voxel, so the coordinates on the upper surface of the beam are all on the z=a/2 plane; When complete cutting (as shown in Fig. 4 (b)), the height of the upper slope of the beam is lower than the upper surface of the voxel, so the coordinates on the upper surface of the voxel are all on the plane where the slope is located; When partial cutting (as shown in Fig. 4 (c)), the intersection of the upper slope of beam and upper surface of the voxel is a cutting line.So partial cutting can be split into no cutting and complete cutting.
In the case of complete cutting, the z value in the surface coordinates of the voxel can be solved by the principle of similar triangle.Set point M'(xM',yM',0) which is on the XOY plane is the projection point M(xM,yM,zM) which is on the plane of the slope.The length of probe is d, the central point of detector is D(xD,yD,0), source is S(xS,yS,0), distance between source and detector is l.Linking N and M': In the case of partial cutting, the expression of the cutting line and the voxel index of the cutting line are first required.Set projection line of cutting line is g=f(x), it is perpendicular to angular bisector of rays beam on the XOY plane and crossing point((xd+xs) •a/d,(yd+ys) •a/d), so the g=f(x) is: Then using the method of calculating the partial coverage index and coordinate (as 2.1.2) to obtain the z value.If z<a/2, the real z value is (x,y,z),else (x,y,a/2).

Calculate the intersecting volume of voxel and ray beams
The part of the ray beam and the voxel is an irregular convex polyhedron.The volume of a convex polyhedron in a three-dimensional space is equal to the sum of the sum of each polygon on which each of its surfaces takes a positive determinant.[14] ,that is, for convex polyhedron A0 A1A2… AN, the volume is: (4)

Optimize volume integral model using GPU parallel computing
The traditional CPU executes the program in a serial way, which is mainly used to control logic.And GPU uses multithreading to calculate large and dense data, and its speed and throughput are higher than CPU.
Since the computation of each ray beam's projection data on detector is independent of each other in volume integral model, we can use GPU parallel computing to optimize volume integral model and reduce projection time.

Obtain the approximate time density curve
Paper [13] confirmed that the time density curve obtained by the contrast agent regiment in functional CT is usually described by a gamma function.In Paper [5], the reverse tangent function is introduced to simulate the trend of the TDC curve reaching its peak.Paper [7] considers the initial density value before the contrast agent is not injected to characterize the actual hemodynamic condition of the organism.The model expression is: (5) X(t) is the density value of voxel changing with time; P indicates the density of objects when no contrast agent is injected; t is the time after injection of contrast agent; T is the time when the contrast agent enters the scanning layer; a, b and c are the model parameters.

Obtain CBCT perfusion parameters
The contrast agent in artery indicates Ca(t).After entering the tissue, the concentration of contrast agent will rise sharply in a short time, then smooth, and finally return to the state of no perfusion.
The concentration changes of the contrast agent in the tissue Q(t) can be calculated by the product of the tissue blood flow BF and the tissue residency function R(t) ,and the TDC of the tissue can be detected by the CT device.
The relation between R(t) and time t is: MTT indicates the average time for contrast agents to flow out of the veins after the artery enters the tissue.Ve is volume of intercellular space.E is associated with surface permeability PS and blood flow BF: Finally, blood volume BV can be obtained by using the blood flow BF and the average of time multiplied by time MTT:

Experimental materials and preparation
The experimental image is a New Zealand white rabbit with VX2 tumor implanted in legs at the princess Margaret Hospital of Toronto.After two weeks of tumor implantation, 4 groups of 120 DICOM format sequence images (Fig. 5) were obtained, and the specific parameters were as Table 1.The distance from the ray source to the axis of rotation(mm) 541 The distance between the ray source and the detector (mm) 949 Image pixel size 512*512 As a result of the calculation of speed, a tumor area (9 x 9 pixel size), such as Fig. 6, is selected as the object of study.Graphics cards NVDIA GeForce GTX 960M

Projection data and cost
The contrast point integral model [7], the line integral model [8] and the area integral model [9] are used to get the CBCT projection sine of the second slices, as shown in Fig. 7.
The figure under the visible volume integral model reduces most of the saw tooth artifact, and the image quality is better than the other three models.Table 3 shows the cost of projection under the four models.The volume integral model brings a lot of computational cost while improving the accuracy of the image.However, after the use of GPU optimization, the calculation of the volume integral model is greatly reduced.

Fig. 2 .
Fig. 2. The coverage of voxel and ray beam on the XOY plane.1.Calculate voxel index The slice is regarded as a m×n grid, and the ray is regarded as an equation y=kx+b, and the two intersection points p1(x1,y1), p2(x2,y2) of the ray and the slicing border.The coordinate range of the ray can be calculated by two intersection points xmin=min(x1,x2), xmax=max(x1,x2), ymin=min(y1,y2), ymax=max(y1,y2).First from the X direction, set x=xmin+a ,enter the loop condition.when x<xmax,(x,kx+b) will be stored in the vertex set {point}.(ceil(x/a),floor(kx+b)/a) will be stored in the index set {index_part}.Then set x=x+a , continue to loop, when x≥xmax end the loop.Then from Y direction, set y=ymin+a, enter the loop condition.When y<ymax, ((y-b)/k,y) will be stored in the vertex set {point}.(ceil(x/a),floor((kx+b)/a)) will be stored in the index set {index_part}.Then set y=y+a , continue to loop, when y≥ymax end the loop.

Fig. 3 .
Fig. 3. Graph of relation between partial coverage voxel and ray beam.

Fig. 4 .
Coverage of ray beams and voxels in space.

Fig. 6 .
Fig. 6.Selected regions of interest in the tumor.The experimental simulation environment parameters are shown in Table2.

7 .
(a) Point integral model (b) Line integral model (c) Area integral model (d) Volume integral model Fig.The sine map of CBCT projection data.

Table 2 .
Simulation environment parameters

Table 3 .
Projection costs of various projection models (units: seconds)