Reconstruction method for motion compensated high quality 4D-CBCT image based on bilateral filtering

Abstract

The present disclosure provides a reconstruction method for motion compensated high quality four-dimensional (4D)-Cone Beam Computed Tomography (CBCT) image based on bilateral filtering, including the following steps: step 1, building a motion compensated three-dimensional (3D)-CBCT reconstruction model based on a simultaneous algebraic reconstruction technique (SART) to reconstruct a high quality CBCT image of a reference phase (phase 0%); step 2, creating an iterative optimization procedure and estimating a 4D Deformable Vector Field (DVF) model between phase 0% and other 4D phase images, thereby obtaining an exact 4D-DVF inclusive of an inverse sliding motion on a surface of a locomotive organ; and step 3, successively deforming the high quality phase 0% image in accordance with an optimized 4D-DVF to obtain a sequence of final high quality 4D-CBCT images. This method permits exact reconstruction of high quality 4D-CBCT images without changing the hardware structure of an existing linear accelerator for conventional radiotherapy.

Claims

1. A reconstruction method for motion compensated high quality four-dimensional (4D)-Cone Beam Computed Tomography (CBCT) image based on bilateral filtering, the method allowing for exact modeling reconstruction for an inverse sliding motion on a surface of a locomotive organ and comprising the following steps: building a motion compensated three-dimensional (3D)-CBCT reconstruction model based on a simultaneous algebraic reconstruction technique (SART) to reconstruct a high quality CBCT image of a reference phase, the reference phase being phase 0%; creating an iterative optimization procedure and estimating a 4D Deformable Vector Field (DVF) model between phase 0% and each other 4D phase image, thereby obtaining an exact 4D-DVF inclusive of the inverse sliding motion on the surface of the locomotive organ; and successively deforming the high quality phase 0% image in accordance with an optimized 4D-DVF to obtain a sequence of final high quality 4D-CBCT images; wherein DVF at each 4D phase is optimized by an algorithm that generates the optimized DVF by registration between 1) measured projections of the phase and 2) simulated forward projection of a reference phase 0% image deformed to the phase; an energy function used in DVF optimization by projection registration is as follows:
f.sub.1(v.sup.0.fwdarw.t)=∥p.sup.t−Aμ.sup.0(x+v.sup.0.fwdarw.t)∥l.sup.2.sub.2+βφ(v.sup.0.fwdarw.t)
f.sub.2(v.sup.t.fwdarw.0)=∥p.sup.0−Aμ.sup.t(x+v.sup.t.fwdarw.0)∥l.sup.2.sub.2+βφ(v.sup.t.fwdarw.0)
s.t.v.sup.0.fwdarw.t(x+v.sup.t.fwdarw.0)+v.sup.t.fwdarw.0=v.sup.t.fwdarw.0(x+v.sup.0.fwdarw.t)+v.sup.0.fwdarw.t=0 wherein f.sub.1 and f.sub.2 represent symmetric inverse continuous energy functions, with the superscript 0 representing phase 0% and the superscript t representing other phase t; A is a matrix for providing a forward projection of the deformed phase 0% image; v.sup.0.fwdarw.t represents a forward DVF; v.sup.t.fwdarw.0 represents an inverse DVF; ∥p.sup.t−Aμ.sup.0(x+v.sup.0.fwdarw.t)∥l.sup.2.sub.2 and ∥p.sup.0−Aμ.sup.t(x+v.sup.t.fwdarw.0)∥l.sup.2.sub.2 are data fidelity terms; φ(v.sup.0.fwdarw.t) and φ(v.sup.t.fwdarw.0) are corresponding normalization terms; the inverse continuity of DVF is guaranteed by the third row in the above equations; the smoothness of the smoothing term is controlled by β; and taking modeling of inverse sliding motions on surfaces of locomotive organs into consideration, φ(v) is defined as: $\phi \left(v\right)=\underset{v\in R^3}{\overset{}{.Math.}}\underset{i=1}{\overset{3}{.Math.}}\underset{j=1,y_j\in N\left(x_j\right)}{\overset{3}{.Math.}}\left(G_x\left(x_j,y_j,{\sigma }_x^2\right).Math.G_{\mu }\left({\mu }^t\left(x_j\right),{\mu }^t\left(y_j\right),{\sigma }_{\mu }^2\right).Math.G_{v_i}\left(v_i\left(x_j\right),v_i\left(y_j\right),{\sigma }_v^2\right){\left(\frac{\partial v_i}{\partial x_j}\right)}^2\right)$ where $\frac{\partial v_i}{\partial x_j}$ represents the difference of a neighborhood of each voxel in DVF, j represents a latitude index of the voxel in the image, and i represents an index of the DVF along the three directions of x, y, and z; specific forms of three sub-filter kernels are as follows: $G_x\left(x_j,y_j,{\sigma }_x^2\right)=\exp \left(-\frac{{\left(x_j-y_j\right)}^T\left(x_j-y_j\right)}{2{\sigma }_x^2}\right)$ $G_{\mu }\left({\mu }^t\left(x_j\right),{\mu }^t\left(y_j\right),{\sigma }_{\mu }^2\right)=\exp \left(-\frac{{.Math.{\mu }^t\left(x_j\right)-{\mu }^t\left(y_j\right).Math.}^2}{2{\sigma }_{\mu }^2}\right)$ $G_{v_i}\left(v_i\left(x_j\right),v_i\left(y_j\right),{\sigma }_v^2\right)=\exp \left(-\frac{{\left(v_i\left(x_j\right)-v_i\left(y_j\right)\right)}^T\left(v_i\left(x_j\right)-v_i\left(y_j\right)\right)}{2{\sigma }_v^2}\right)$ G.sub.x is the Gaussian kernel on the spatial domain with the variance of σ.sub.x.sup.2—; G.sub.μ is another image domain based Gaussian kernel with the variance of σ.sub.μ.sup.2; G.sub.v is the DVF domain Gaussian kernel with the variance of σ.sub.v.sup.2; ‘x.sub.j’ is the kernel's center voxel location, and ‘y.sub.j’ is the neighborhood voxel locations within $N\left(x_j\right);\frac{\partial v_i}{\partial x_j}$ represents a difference of a neighborhood of each voxel in DVF; the gradient ∇φ(v)|.sub.v is calculated within the 3×3×3 neighborhood surrounding each voxel of interest; the gradient of φ(v) is given by: $\nabla \phi \left(v\right)=\underset{v\in R^3}{\overset{}{.Math.}}\underset{i=1}{\overset{3}{.Math.}}\underset{j=1,y_j\in N\left(x_j\right)}{\overset{3}{.Math.}}\left(G_x\left(x_j,y_j,{\sigma }_x^2\right).Math.G_{\mu }\left({\mu }^t\left(x_j\right),{\mu }^t\left(y_j\right),{\sigma }_{\mu }^2\right).Math.G_v\left(v_i\left(x_j\right),v_i\left(y_j\right),{\sigma }_v^2\right)\left(v_i\left(x_j\right)-v_i\left(y_j\right)\right)\right).$

2. The reconstruction method for motion compensated high quality 4D-CBCT image based on bilateral filtering according to claim 1, wherein during modeling in step 1, the optimized 4D-DVF between phase 0% and other 4D phases is substituted into the SART reconstruction procedure to obtain an improved motion compensated SART algorithm, thereby obtaining a high quality phase 0% initial image based on all 4D projections; and a mathematical model of motion compensated SART is described below: letting p.sup.t=(p.sub.1.sup.t, p.sub.2.sup.t, . . . , p.sub.I.sup.t) represent a logarithmically compressed line integral of 4D-CBCT projections of phase t and μ.sup.t=(μ.sub.1.sup.t, μ.sub.2.sup.t, . . . , μ.sub.J.sup.t) represent an attenuation coefficient for a phase t image; and the mathematical model of motion compensated SART is expressed as: ${\mu }_j^{0,\left(k+1\right)}={\mu }_j^{0,\left(k\right)}+\frac{\underset{\mathrm{jn}}{.Math.}{d}_{\mathrm{jn}}^{t.fwdarw.0}\underset{i}{.Math.}\left[a_{\mathrm{in}}\frac{p_i^t-\underset{n}{.Math.}{a}_{\mathrm{in}}{\mu }_n^{t,\left(k\right)}}{\underset{n=1}{\overset{J}{.Math.}}{a}_{\mathrm{in}}}\right]}{\underset{i}{.Math.}{a}_{in}}$ wherein k is the number of iterations; j is a voxel index for the phase 0% image; n is a voxel index for the phase t image; a.sub.in is an intersection length of forward projection simulation i with image voxel n, which is obtained by Digitally Reconstructed Radiograph (DRR); d.sub.jn.sup.t.fwdarw.0 represents an inverse DVF for deforming the phase t image to phase 0; introducing DVF reflects the idea of motion compensation; and the phase 0% initial image μ.sub.j.sup.0,(0) is obtained via Total Variation (TV) minimization reconstruction.

3. The reconstruction method for motion compensated high quality 4D-CBCT image based on bilateral filtering according to claim 1, wherein on the basis of the high quality phase 0% CBCT image and the optimized 4D-DVF thus obtained, a sequence of final high quality 4D-CBCT images is obtainable by deforming the phase 0% image; and the mathematical description of this procedure is presented below:
μ.sub.n.sup.t,(k)=Σ.sub.j d.sub.jn.sup.0.fwdarw.tμ.sub.j.sup.0,(k) wherein μ.sub.n.sup.t,(k) represents that a phase t CBCT image is obtainable by deforming μ.sub.j.sup.0,(k) using 4D-DVFs, i.e., d.sub.jn.sup.0.fwdarw.t.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a comparison diagram of section curve dose rates of an FFF beam and a flattened beam in the prior art.

(2) FIG. 2 is a block diagram of a technical route of motion compensation for high quality 4D-CBCT based on bilateral filtering.

(3) FIG. 3 is a comparison diagram of reconstruction results of 4D Non-Uniform Rational B-Splines (NURBS)-based Cardiac-Torso (NCAT) phantoms.

(4) FIG. 4 is a comparison diagram of reconstruction results of 4D motion trajectory (in z-axis direction) in ROI of NCAT phantom.

(5) FIG. 5 is a diagram of verification results of algorithm with clinical patient data.

(6) FIG. 6 is a comparison diagram of reconstruction results from clinical data of a second case.

(7) FIG. 7 is a comparison diagram of reconstruction results from clinical data of a third case.

DETAILED DESCRIPTION

(8) A reconstruction method for motion compensated high quality 4D-CBCT image based on bilateral filtering includes the following steps.

(9) Step 1, Collection of 4D-CBCT Projections

(10) Firstly, 4D-CBCT projections are collected. A waveform of amplitudes and phases of the respiratory motion of a patient over time during treatment can be obtained by using an existing mature camera-based respiratory motion detection system. With the collection of the 4D-CBCT projections, each frame of collected projections may be labeled with a corresponding time tag by a monitored respiratory signal. After the completion of the collection of all projections, the CBCT projections with time and amplitude tags of respiratory signals can be subjected to 4D phase segmentation by respiratory amplitude or phase signal, thereby obtaining 4D-CBCT projections.

(11) This step is a data pre-processing procedure. In specific implementation, the 4D-CBCT projections may be obtained using Varian Respiration Positioning Management (RPM) system. The system includes a globular infrared light reflector holder placed on a patient's abdomen and an infrared camera fixed to the ceiling of a treatment room to monitor the motion amplitudes of a globular infrared light reflector. Then, the patient's respiratory amplitude and phase signals over time can be recorded by the infrared camera. The recorded respiratory signals are distributed to the projections inclusive of the patient's respiratory motion collected by the CBCT system, equivalent to labeling each collected projection with the time tag. Subsequently, 4D-CBCT projections can be obtained by 4D grouping of the projections according to the time tags of a respiratory motion curve.

(12) Step 2, Initial 4D-CBCT Image Reconstruction and Initial 4D-DVF Generation

(13) Initial 4D-CBCT reconstruction and initial 4D-DVF generation are carried out using Total Variation (TV) minimization and Demons registration algorithms, respectively.

(14) Firstly, a sequence of initial 4D-CBCT images is reconstructed from the collected 4D-CBCT projections by TV minimization. The TV minimization reconstruction method is adopted because effective noise removal from reconstructed images can be achieved by this approach with image boundaries well retained. The number of projections at each phase is about 50-60 after the 4D phase segmentation of the projections. When initial 4D-CBCT reconstruction is carried out with such sparse projections by the traditional analytical method, strong streak artifacts and extremely poor image quality are caused. The TV minimization method can prevent images from being polluted by the streak artifacts, allowing for further improvement on the quality of subsequent 4D-CBCT images.

(15) After the initial 4D-CBCT images are obtained, initial 4D-DVF at reference phase % versus other 4D phases is generated using Demons registration algorithm. Demons registration algorithm is a classical deformable image registration algorithm based on similarity between image pixel values, which can generate a 3D-DVF between two matching images. The initial 4D-DVF obtained using Demons registration algorithm can be used to trigger further 4D-DVF optimization, thereby being conducive to obtaining the optimized 4D-DVF rapidly.

(16) Step 3, Iterative 4D-DVF Optimization

(17) The DVF between a reference phase 0% image and one of other 4D phases is generated by registration between measured projections of the target phase and the forward projection of the deformed phase 0%. The procedure is as shown by formula (2): p.sup.t represents the measured projections of the target phase. μ.sup.0(x+v.sup.0.fwdarw.t) represents DVF pixel v.sup.0.fwdarw.t, causing the phase 0% image μ.sup.0 to be deformed from the reference phase 0% to the target phase t. A forward projection matrix A is added in front of the deformed phase 0% image to obtain the forward projection of the deformed image. ∥p.sup.t−Aμ.sup.0(x+v.sup.0.fwdarw.t)∥l.sup.2.sub.2 represents l.sub.2 norm of differences between the measured projections of the target phase t and the forward projection of the phase 0% image deformed to the target phase t, namely difference square. This forms the fidelity term of a DVF optimization object function. The design of the penalty term βφ(v.sup.0.fwdarw.t) of the DVF optimization object function focuses on natural DVF smoothing and maintaining of extraction of locally non-uniform sliding motion of a locomotive organ. The specific smoothing term is designed as equation (3). The smoothing term mainly includes three Gaussian filter kernels based on bilateral filtering, namely spatial domain-based Gaussian kernel G.sub.x, image pixel domain-based Gaussian kernel G.sub.μ, and DVF domain-based Gaussian kernel G.sub.v. During optimization, the effect of the whole bilateral filtering can be adjusted by setting variances corresponding to different Gaussian kernels, such as variance σ.sub.x.sup.2 corresponding to G.sub.x, variance σ.sub.μ.sup.2 corresponding to G.sub.μ, and variance σ.sub.v.sup.2 corresponding to G.sub.v. The estimation of the globally uniform motion vector field of the locomotive organ and the capture of the locally nonuniform inverse sliding motion occurring on the surface of the locomotive organ can be eventually enabled based on DVF.

(18) As described above, the DVF between two phases is adjusted by registration between the projections at the target phase and the forward projection of the deformed reference phase. In this process, value assignment for different variances of the bilateral filtering kernels is required. In specific implementation, algorithm verification is firstly performed on a 4D Non-Uniform Rational B-Splines (NURBS)-based Cardiac-Torso (NCAT) phantom in ideal conditions. It was found that the optimal result can be obtained at σ.sub.x=3 mm. When σ.sub.x is too large or too small, DVF would be over-smoothed or no local spatial feature would be captured. For the variance σ.sub.μ on the image pixel domain, let σ.sub.μ be equal to the difference between a mean of pixels at the lung tissue and a mean of pixels in the vicinity of the thoracic wall so as to enable natural transition of image pixels from the thoracic wall to the lung tissue. For the NCAT phantom, σ.sub.μ=0.03 mm.sup.−1 satisfies this condition and the reconstruction result is ideal. Finally, the key is how to determine the suitable variance σ.sub.v on the DVF domain. In theory, σ.sub.v should be greater than the DVF difference between any two points within the spatial range of σ.sub.x unless σ.sub.x includes a sudden pixel change region of the thoracic wall. In the sudden pixel change region of the thoracic wall, σ.sub.v should be smaller than the DVF pixel difference of this region. Thus, small pixel motions can be filtered out while large pixel motions at the thoracic wall can be retained. To avoid over-segmentation, only large sharp and discontinuous motions (e.g., inverse sliding motion on the surface of a locomotive organ) are set to be retained. By observing clinical patient data, the amount of movement of 10 mm can be relatively regarded as a large sliding motion threshold. Under this assumption, it was found that a good reconstruction result based on the NCAT phantom data can be obtained at σ.sub.v=2.5 mm.

(19) On the basis of the verification of the algorithm performance based on the NCAT phantom, the algorithm is further verified using the clinical patient data. Since the clinical data is collected in real environment, there may exist a series of interfering factors including noise, scattering, beam hardening, etc. The corresponding filter kernel variance that fits the characteristics of preliminary clinical data slightly differs from the aforementioned variance that fits the NCAT phantom. It was found that σ.sub.x=3 mm, σ.sub.μ=0.02 mm.sup.−1 and σ.sub.v=2 mm can well fit the clinical data, and a good clinical reconstruction result can be obtained.

(20) Iterative optimization estimation of DVF is carried out with the above variances. A non-linear conjugate gradient operator is used in the optimization of the object function, and therefore, the gradient of the object function needs to be calculated. The gradient of the penalty term in the object function is derived as shown by equation (4).

(21) Step 4, Continuously Iterative Updating Reconstruction at Reference Phase %

(22) After each 4D-DVF iteration, updating reconstruction of the reference phase % image is carried out, so that the quality of the reference phase image can be continuously improved. The reconstruction procedure is as described by equation (1), where the DVF pixel d.sub.jn.sup.0.fwdarw.t is updated with each 4D-DVF update.

(23) In specific implementation, the optimal estimation of 4D-DVF and the updating of the reference phase % image are carried out alternately. In other words, 4D-DVF is firstly optimized by n iterations (n can be set to 10-20), and a median 4D-DVF at this time is saved. The median 4D-DVF at this time is then substituted into a motion compensated SART algorithm to reconstruct the phase 0% image after the 4D-DVF optimization by n iterations. After next 4D-DVF optimization by n iterations, a second motion compensated SART reconstruction is carried out to obtain the phase 0% image after the second update. After each iterative 4D-DVF optimization, the value of the object function at each phase is recorded. After the mth motion compensated SART reconstruction iteration, if the 4D-DVF optimization object function at each phase trends to converge, the overall 4D-DVF optimization and motion compensated SART reconstruction would be stopped. Moreover, the current optimal 4D-DVF and a high quality phase 0% CBCT image obtained by m motion compensated SART reconstructions are saved.

(24) Step 5, Generation of a Sequence of Final High Quality 4D-CBCT Images

(25) As the object function decreases continuously and eventually trends to be stable during 4D-DVF optimization, the reference phase 0% image reconstructed by continuous iterations is deformed with the resulting optimized 4D-DVF. A sequence of high quality 4D-CBCT images is finally obtained. The whole procedure is as shown by equation (5).

(26) The block diagram of the technical route of the above steps is as shown in FIG. 2.

EXAMPLE 1

(27) The validity of the algorithm was verified on the 4D NCAT phantom first. The respiration period of the 4D phantom was 4 seconds. The superior internal (SI) motion at the edge of the thoracic diaphragm had a maximum of 20 mm, and the anteroposterior (AP) motion in the chest had a maximum of 12 mm. The phantom had ten phases, and 20 projections were simulated at each phase for DVF estimation and 4D-CBCT reconstruction. The phantom image size was set to 256*256*150, and the voxel size was set to 1*1*1 mm.sup.3. The projection size was set to 300*240*20, and the projection voxel size was set to 1*1*1 mm.sup.3.

(28) FIG. 3 shows a reconstructed image of phase 40% with and without sliding motion constraint. FIG. 3(a) shows reconstructed sagittal images of phase 40% without bilateral filtering; FIG. 3(b) shows a reconstructed identical sagittal slice with bilateral filtering; FIG. 3(c) shows the phantom ground truth; and FIG. 3(d), FIG. 3(e) and FIG. 3(f) show the Region Of Interest (ROI), where sliding motions occurred at the edges of the heart and the veins. ROI showed clear boundaries of the edges of the heart (indicated by the arrows in figures (d), (e) and (f)) obtained with bilateral filtering. In addition, the veins (marked with yellow dotted lines) were reconstructed more accurately with bilateral filtering (for example, the length of the veins was reconstructed better). FIG. 3(g), FIG. 3(h) and FIG. 3(i) show the difference between reconstructions of ribs with and without bilaterally filtering. In FIG. 3(h), the top edges 1 and 2 of the ribs were matched with the reference positions of the ribs in the ground truth as shown in FIG. 3(i), and bilateral filtering was omitted in FIG. 3(g).

(29) To detect motion errors caused by sliding motion images at the boundaries of a locomotive organ, a 4D motion trajectory in z-axis direction was extracted from the edge of the heart in the above figure (see the dotted line at the marked position in FIG. 3(f)). In each 4D image, the correspondingly detected dotted line positions were used for motion trajectory rendering. FIG. 4 shows that the trajectory extracted with consideration of sliding based on bilateral filtering was much better in degree of matching with the 4D trajectory ground truth than the trajectory with no consideration of this factor. Root Mean Square Error (RMSE) and Maximum motion Error (MaxE) were calculated with each motion from the trajectory. RMSE was defined as follows:

(30) $\begin{array}{cc}\mathrm{RMSE}=\sqrt{\frac{1}{9}\times \underset{\mathrm{ph}=1}{\overset{9}{.Math.}}{\left({\mathrm{pos}}_{\mathrm{ph}}^R-{\mathrm{pos}}_{\mathrm{ph}}^T\right)}^2}& \left(7\right)\end{array}$

(31) where pos.sub.ph.sup.R represents a feature structure position in the ph.sup.th phase image in the reconstructed image; pos.sub.ph.sup.T represents the position of the same feature structure in the corresponding phase in the ground truth image; ph represents a phase index; and there were 10 4D phases here. MaxE was the maximum error among the 10 phases.

(32) The results showed that based on the result of bilateral filtering, the RMSE and MaxE of the 4D trajectory were 0.796 mm and 1.02 mm, respectively. Compared with the result of original reconstruction without bilateral filtering, the RMSE and MaxE of the trajectory were 2.704 mm and 4.08 mm, respectively.

(33) A preliminary clinical patient data test was carried out on the basis of the data of the NCAT phantom. CBCT projections of two patients with lung cancer were collected, and a preliminary test was performed on the sliding compensated 4D-CBCT reconstruction algorithm based on bilateral filtering. To find out the ground truth for the clinical data, quantitative comparison was carried out and the clinical ethics permission by Institutional Review Board (IRB) was submitted. The data of each patient was projections collected in a full fan mode. The scanning time was 4-6 minutes. About 2000 projections were obtained. The obtained projections were divided into 10 phases using corresponding Respiration Positioning Management (RPM) signals during collection. Division in this way resulted in averagely about 200 projections at each phase, and a high quality 4D-CBCT image was reconstructed using the TV minimization reconstruction technique, which could be used as a reference image for clinical result quantification. This would be not allowed under the conventional clinical scan protocol. The purpose was merely to collect sufficient projections, ensuring that the 4D phase segmentation could result in sufficient projections at each phase to enable reconstruction of a high quality 4D-CBCT clinical data image as reference for clinical effect evaluation and verification on the algorithm.

(34) In consideration of only about 670 projections collected at most by rotation of 360 degrees under one-minute CBCT scan protocol in clinical practice and actually only about 60 projections at each 4D phase after 4D phase segmentation, All the collected 4D projections were undersampled so as to simulate the real one-minute CBCT scan scenario. In other words, the average projections at each phase were under-sampled to about 40. Subsequently, a series of parallel control 4D-CBCT reconstructions were obtained using different reconstruction methods and then compared with the proposed motion compensated reconstruction based on bilateral filtering. The used reconstruction methods were traditional Feldkamp-David-Kress (FDK) reconstruction, TV minimization reconstruction, and reconstruction with no consideration of sliding motion modeling in algorithm.

(35) Respective results were shown in FIG. 5. FIG. 5(a) and FIG. 5(b) show sagittal images of reconstruction results for the first patient with/without sliding motion constraint based on bilateral filtering. FIG. 5(c) shows a patient reference image reconstructed from all sampled projections by television. As indicated by arrows, there was a tumor next to the thoracic wall with a small hole between the tumor and the thoracic wall. Compared with the reference in FIG. 5(c), the small cavity in FIG. 5(a) was apparently obscure. Compared with the control, the boundary of the tumor as shown in FIG. 5(b) was more apparent than that shown in FIG. 5(a). This indicated that more comparable reference results were obtained by reconstruction based on bilateral filtering than by reconstruction without bilateral filtering. For the sake of comprehensive comparison, FIG. 5(d) and FIG. 5(e) show respective reconstructed slices by TV minimization and FDK methods, respectively. The two results indicated that the structure in the image reconstructed by TV minimization was over-smoothed, and the content in the FDK image was seriously worsened by noise pollution.

(36) On the other hand, FIG. 5(f), FIG. 5(g) and FIG. 5(h) show the boundary of the pulmonary rib cage of this patient, thus allowing for comparison of rib motions. The dotted line indicates the upper edge position of the rib. As indicated by arrows, in reconstruction without bilateral filtering (such as FIG. 5(f), the top edge of the rib was slightly shifted upwards by about 2 mm in FIG. 5(f) compared with the bilaterally filtered reconstruction result (FIG. 5(f)) and the reference rib position in the patient reference image FIG. 5(h). Therefore, the real position of the rib in the reconstructed image was well corrected by bilateral filtering.

(37) FIG. 6 shows a coronal image reconstructed from the clinical data of the second patient. FIG. 6(a) shows an original motion compensated reconstructed image without bilateral filtering. FIG. 6(b) shows a bilaterally filtered reconstruction result. FIG. 6(c) shows a clinical patient reference image reconstructed from all projections by TV minimization. FIG. 6(d), FIG. 6(e) and FIG. 6(f) show zoomed-in ROIs (see ROI box in FIG. 6(a)) at the same position in FIG. 6(a), FIG. 6(b) and FIG. 6(c), respectively. It was clearly indicated that the cardiopulmonary boundary (marked by the arrow) was clearer than that in FIG. 6(d), as compared with FIG. 6(f). This further proved the advantage of bilateral filtering, that is, more accurate reconstruction results can be obtained as compared with original motion compensated results.

(38) FIG. 7 shows comparison of reconstruction results for a third patient. Images were reconstructed of phase 40% for the third patient. As shown in the figure, (a) was original reconstruction without bilateral filtering; (b) was a bilaterally filtered reconstruction result; and (c) was a reference for the patient. The results showed that the fiber like structure indicated by yellow arrows could be more apparent by bilaterally filtered reconstruction than by original reconstruction, as compared with the reference for the patient.

(39) Numerical values of quantitative analysis of relative reconstruction errors caused by 1) traditional FDK reconstruction method, 2) TV minimization reconstruction method, 3) original motion compensated reconstruction method with no consideration of sliding motion modeling and 4) bilaterally filtered sliding motion compensated reconstruction method with three pieces of clinical data, respectively, were summarized in table 1. The results showed the reconstruction error caused by the proposed bilateral filtering method was smallest.

(40) TABLE-US-00001 TABLE 1 Comparison of Relative Errors of Reconstructions from Clinical Data RE FDK TV Original Recon Bilateral Filtering Recon Patient 1 29.63% 9.88% 7.05% 6.61% Patient 2 30.55% 9.01% 7.68% 7.06% Patient 3 30.45% 8.85% 7.21% 7.10%

(41) Finally, it should be noted that the above examples are only intended to explain, rather than to limit, the technical solution of the present disclosure. Although the present disclosure is described in detail with reference to the preferred examples, those ordinarily skilled in the art should understand that modifications or equivalent substitutions may be made to the technical solution of the present disclosure without departing from the spirit and scope of the technical solution of the present disclosure, and such modifications or equivalent substitutions should be included within the scope of the claims of the present disclosure.