Reconstructing a 3D Medical Image from a Few 2D Projections Using a B-Spline-Based Deformable Transformation

Abstract

(1) Background: There was a need for 3D image reconstruction from a series of 2D projections in medical applications. However, additional exposure to X-ray projections may harm human health. To alleviate it, minimizing the projection numbers is a solution to reduce X-ray exposure, but this would cause significant image noise and artifacts. (2) Purpose: In this study, a method was proposed for the reconstruction of a 3D image from a minimal set of 2D X-ray projections using a B-spline-based deformable transformation. (3) Methods: The inputs of this method were a 3D image which was acquired in previous treatment and used as a prior image and a minimal set of 2D projections which were acquired during the current treatment. The goal was to reconstruct a new 3D image in current treatment from the two inputs. The new 3D image was deformed from the prior image via the displacement matrixes that were interpolated by the B-spline coefficients. The B-spline coefficients were solved with the objective function, which was defined as the mean square error between the reconstructed and the ground-truth projections. In the optimization process the gradient of the objective function was calculated, and the B-spline coefficients were then updated. For the acceleration purpose, the computation of the 2D and 3D image reconstructions and B-spline interpolation were implemented on a graphics processing unit (GPU). (4) Results: When the scan angles were more than 60°, the image quality was significantly improved, and the reconstructed image was comparable to that of the ground-truth image. As the scan angles were less than 30°, the image quality was significantly degraded. The influence of the scan orientation on the image quality was minor. With the application of GPU acceleration, the reconstruction efficiency was improved by hundred times compared to that of the conventional CPU. (5) Conclusions: The proposed method was able to generate a high-quality 3D image using a few 2D projections, which amount to ~ 20% of the total projections required for a standard image. The introduction of the B-spline-interpolated displacement matrix was effective in the suppressing noise in the reconstructed image. This method could significantly reduce the imaging time and the radiation exposure of patients under treatment.

1. Introduction

A 3D image, such as computed tomography (CT) or cone-beam CT (CBCT), can be generated from a series of 2D projections acquired in a single X-ray scan. They have been widely used in radiation therapy for target localization [1,2,3,4]. CBCT imaging delivers a radiation dose to a large volume of the patient’s body and collects the residual dose on the detectors. As the accumulated imaging doses to patients could be significant and increase the risk of secondary cancer induction [5,6,7], there is a demand for low-dose CBCT with an adequate image quality in clinic. For reducing the imaging dose, two parameters (exposure level and projection number) were needed to be adjusted. However, reducing the exposure level will lead to increased noise in the reconstructed images, while reducing the projection number will violate the Nyquist–Shannon sampling theorem, leading to serious streak artifacts.

To solve these issues, compressed sensing methods based on the total variation (TV) have been developed to improve the image quality of low-dose CT/CBCT [8,9,10,11,12]. Although these methods successfully reduce the noise and streak artifacts, they also tend to over-smooth the edge information. Alternatively, other compressed-sensing-based methods used prior images for low-dose CBCT [13]. In the prior-image-constrained compressed sensing (PICCS) method, prior images were used as an additional constraint to minimize the image’s total variation [14]. Adaptive-prior-image-constrained compressed sensing (APICCS) was then proposed to enhance the image quality of PICCS [15,16]. The rigid or deformable registration between the prior image and the reconstructed volume was then incorporated to improve the image quality [17,18]. Moreover, a motion model was introduced to further reduce the scan angle needed for CBCT reconstruction [19,20,21].

Based on the compressed sensing theory, images can be reconstructed from limited-view and limited-angle projections [22]. However, these methods generally require the sparseness prior of the image, which may be not valid. Based on the deformable transform, images can be reconstructed from prior images [23,24,25,26,27]. Lei et al. proposed a deformable field map method to reconstruct 3D images from limited-angle 2D projections. The missing information was made up of highly correlated prior images. However, due to the high fluctuation of the gradient derived from the objective function, the resulted deformation field map comes with a higher deviation, which could result in the failure of the image transformation [28,29]. To solve this issue, the gradient or deformation field map should be smoothed in a rational way.

Several smoothing methods for deformable transformations can be found in image registration algorithms, including demons registration, viscous fluid registration, B-spline registration and thin-plate splines. Among them, B-spline is the most popular one because of its flexibility and robustness which provide the ability to represent the complex spatial distribution of the data and image in 3D [30,31]. B-spline is used to define the displacement matrixes which map voxels in a moving image to those in a reference image. The individual voxel movement between the two images is parameterized in terms of the uniformly spaced control points that are aligned with the voxel grid, and the displacement vectors are obtained via the interpolation of the control point coefficients using piecewise continuous B-spline basis functions. B-spline interpolation is computationally intensive and requires hours for higher image resolutions [32].

As there are many computation intensity components in CBCT reconstruction, it is more efficient to implement them on a graphics processing unit (GPU) which can handle multiple similar tasks through parallel computation [33]. Its effectiveness was already demonstrated in accelerating image reconstruction tasks, including DRR [34], CBCT [35], 4DCT [36], etc. The integrated GPU-CPU platform can potentially be applied to different medical image applications, such as medical image segmentation and registration [37,38]. Also, GPUs paved the way for various promising real-time and on-line clinical applications, such as tumor motion tracking [39], respiration motion tracking and cardiac motion tracking [40].

In this study, a B-spline-based deformable transformation method was proposed for 3D CBCT reconstructions. Instead of the real-value form of the displacement matrixes in traditional methods, the B-spline coefficients of the displacement matrixes were employed. These coefficients were optimized by minimizing the objective function and finally interpolated the real-value displacement matrixes. For acceleration purposes, the computation of 2D and 3D reconstructions were implemented on a GPU. The proposed method was evaluated through a clinical head-and-neck case and compared with the traditional method. Finally, the influence of the scan angle and scan orientation on the image quality was analyzed, and the advantage of this method was discussed.

2. Methods

2.1. Deformation-Based Reconstruction Method

The traditional method used the real-value displacement matrixes to transform the prior image into the target image. In the following, V0 and V are the prior image and target image, and D is the displacement matrix. The transformation from V0 to V can be formulated as

$$V=V\left(D\right)=V0\left(i+D_{ijk}^x,j+D_{ijk}^y,k+D_{ijk}^z\right),i,j,k=1,\dots ,n,$$

(1)

where $D_{ijk}^x,D_{ijk}^y,D_{ijk}^z$ are the displacements matrixes along x, y, and z axes, respectively. n is the image size, and ijk are the 3D indexes of voxel. V was transformed from V0 through trilinear interpolation as described below.

$$\begin{array}{c}\overline{i}=\lfloor i+D_{ijk}^x\rfloor ,\overline{j}=\lfloor j+D_{ijk}^y\rfloor ,\overline{k}=\lfloor k+D_{ijk}^z\rfloor \\ \widehat{i}=i+D_{ijk}^x-i_0,\widehat{j}=j+D_{ijk}^y-j_0,\widehat{k}=k+D_{ijk}^z-k_0\hfill \\ u_1=(1-\widehat{k})\ast V{0}_{\overline{i},\overline{j},\overline{k}}+\widehat{k}\ast V{0}_{\overline{i},\overline{j},\overline{k}+1},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em} u_2=(1-\widehat{k})\ast V{0}_{\overline{i},\overline{j}+1,\overline{k}}+\widehat{k}\ast V{0}_{\overline{i},\overline{j}+1,\overline{k}+1}\\ v_1=(1-\widehat{k})\ast V{0}_{\overline{i}+1,\overline{j},\overline{k}}+\widehat{k}\ast V{0}_{\overline{i}+1,\overline{j},\overline{k}+1},\hspace{1em}\hspace{1em}\hspace{1em}{v}_2=(1-\widehat{k})\ast V{0}_{\overline{i}+1,\overline{j}+1,\overline{k}}+\widehat{k}\ast V{0}_{\overline{i}+1,\overline{j}+1,\overline{k}+1}\\ w_1=(1-\widehat{j})\ast u_1+\widehat{j}\ast u_2,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} w_2=(1-\widehat{j})\ast v_1+\widehat{j}\ast v_2\\ V_{ijk}=(1-\widehat{i})\ast w_1+\widehat{i}\ast w_2\end{array}$$

(2)

Here, $\overline{i},\overline{j},\overline{k}$ are the three-dimensional indices of grid adjacent to the voxel coordinates, and $\widehat{i},\widehat{j},\widehat{k}$ are the corresponding fractional parts. $(u_1,u_2)$,$(v_1,v_2)$ and $(w_1,w_2)$ are pairs of weights along three axes.

The fidelity of target image could be guaranteed by approximating its projections to the ground-truth projections. The 2D projections Y can be calculated from V as

$$Y=P(V)=P\left(V\left(D\right)\right),i,j,k=1,\dots ,n,$$

(3)

where P is the projection operator implemented by ray-tracing algorithm. As there are not enough equations in (3), to solve this ill-posed problem, the regularization term was introduced. Here, the bending energy defined as below was employed to enforce the flatness of D.

$$E\left(D\right)={\displaystyle \sum _{k=1}^n{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{m=1}^3\left({\left(\frac{\partial D_{ijk}^x}{\partial x}\right)}^2+{\left(\frac{\partial D_{ijk}^y}{\partial y}\right)}^2+{\left(\frac{\partial D_{ijk}^z}{\partial z}\right)}^2\right)}}}}$$

(4)

The objective function was then established as

$$\min E\left(D\right)\begin{array}{cc}& \end{array}s.t.\begin{array}{cc}& \end{array}\Vert P(V)-Y\Vert =0$$

(5)

This problem could be transformed to an unconstraint problem with a weight µ

$$f\left(D\right)=\Vert P(V)-Y\Vert +\mu \ast E\left(D\right)$$

(6)

The gradient of f(D) could be derived as

$$\nabla f\left(D\right)=2P^{-1}\left(P\left(V\left(D\right)-Y\right)\right)\nabla V\left(D\right)+\mu \ast \nabla E\left(D\right)$$

(7)

where P⁻¹ is the inverse projection operator which backprojects 2D projections onto 3D volume. $\nabla V\left(D\right)$ could be calculated with respect to D along three axes.

$$\begin{array}{c}\frac{\partial V_{ijk}}{\partial D_{ijk}^z}=\frac{\partial w_1}{\partial D_{ijk}^z}*(1-\widehat{i})+\frac{\partial w_2}{\partial D_{ijk}^z}*\widehat{i}\\ =\left(\frac{\partial u_1}{\partial D_{ijk}^z}*\left(1-\widehat{j}\right)+\frac{\partial u_2}{\partial D_{ijk}^z}*\widehat{j}\right)*\left(1-\widehat{i}\right)+\left(\frac{\partial v_1}{\partial D_{ijk}^z}*\left(1-\widehat{j}\right)+\frac{\partial v_2}{\partial D_{ijk}^z}*\widehat{j}\right)*\widehat{i}\\ =\left(\begin{array}{l}\left(-V{0}_{\overline{i},\overline{j},\overline{k}}+V{0}_{\overline{i},\overline{j},\overline{k}+1}\right)*\left(1-\widehat{j}\right)\\ +\left(-V{0}_{\overline{i},\overline{j}+1,\overline{k}}+V{0}_{\overline{i},\overline{j}+1,\overline{k}+1}\right)*\widehat{j}\end{array}\right)*\left(1-\widehat{i}\right)\hspace{0.17em}\hspace{0.17em}+\left(\begin{array}{l}\left(-V{0}_{\overline{i}+1,\overline{j},\overline{k}}+V{0}_{\overline{i}+1,\overline{j},\overline{k}+1}\right)*\left(1-\widehat{j}\right)\\ +\left(-V{0}_{\overline{i}+1,\overline{j}+1,\overline{k}}+V{0}_{\overline{i}+1,\overline{j}+1,\overline{k}+1}\right)*\widehat{j}\end{array}\right)*\widehat{i}\end{array}$$

(8a)

$$\frac{\partial V_{ijk}}{\partial D_{ijk}^y}=\frac{\partial w_1}{\partial D_{ijk}^y}\ast \left(1-\widehat{i}\right)+\frac{\partial w_2}{\partial D_{ijk}^y}\ast \widehat{i}=\left(-u_1+u_2\right)*\left(1-\widehat{i}\right)+\left(-v_1+v_2\right)*\widehat{i}$$

(8b)

$$\frac{\partial V_{ijk}}{\partial D_{ijk}^x}=-w_1+w_2$$

(8c)

$\nabla E$ could be calculated with respect to D as below.

$$\begin{array}{l}\frac{\partial E}{\partial D}=-2(D_{i+1,j,k}^x-D_{i,j,k}^x)+2(D_{i,j,k}^x-D_{i-1,j,k}^x)-2(D_{i,j+1,k}^x-D_{i,j,k}^x)\\ \hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+2(D_{i,j,k}^x-D_{i,j-1,k}^x)-2(D_{i,j,k+1}^x-D_{i,j,k}^x)+2(D_{i,j,k}^x-D_{i,j,k-1}^x)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=12D_{i,j,k}^x-2D_{i+1,j,k}^x-2D_{i-1,j,k}^x-2D_{i,j+1,k}^x-2D_{i,j-1,k}^x-2D_{i,j,k+1}^x-2D_{i,j,k-1}^x\end{array}$$

(9)

When both $\nabla V\left(D\right)$ and $\nabla E\left(D\right)$ were obtained, $\nabla f\left(D\right)$ in Equation (7) could be solved.

2.2. B-Spline Interpolation

The uniform cubic B-spline basis functions were employed in this study. The real-value displacement matrix was interpolated from B-spline coefficients. In 3D case, the deformation field at any given voxel was determined with the 64 control points in the immediate vicinity of the voxel’s housing tile. B-spline interpolation was performed for each voxel with respect to the 64 control point coefficients as follows

$$D_{ijk}^x=\mathsf{\Phi }\left(I_{}^x\right)={\displaystyle \sum _{l=0}^3{\displaystyle \sum _{m=0}^3{\displaystyle \sum _{n=0}^3{\beta }_l\left(u\right){\beta }_m\left(v\right){\beta }_n\left(w\right)I_{lmn}^x}}}$$

(10)

where I^x_lmn is the spline coefficient defining the x component of the displacement vector for one of the 64 control points. N_x, N_y and N_z are voxels per region (VPR) defining the distance between control points in terms of voxels in the x, y and z directions, respectively. The volume was, therefore, segmented by the B-spline control point grid into many equal-sized tiles of dimensions N_x × N_y × N_z. The three-dimensional indices x_t, y_t and z_t of the tile for the voxel at x, y and z was given by $\lfloor i/N_x\rfloor ,\lfloor j/N_y\rfloor ,\lfloor k/N_z\rfloor$, respectively. The coordinates of the voxel within its tile were $u=i/N_x-x_t,v=j/N_y-y_t,w=k/N_z-z_t$.

The uniform cubic B-spline basis functions were given by β_l(u), β_m(v) and β_n(w) along the x, y and z axes, respectively.

$${\beta }_l\left(u\right)=\{\begin{array}{l}{\left(1-u\right)}^3/6:l=0\\ \left(3u^3-6u^2+4\right)/6:l=1\\ \left(-3u^3-3u^2+3u+1\right)/6:l=2\\ u^3/6:l=3\end{array}{\beta }_m\left(v\right)=\{\begin{array}{l}{\left(1-v\right)}^3/6:m=0\\ \left(3v^3-6v^2+4\right)/6:m=1\\ \left(-3v^3-3v^2+3v+1\right)/6:m=2\\ v^3/6:m=3\end{array}{\beta }_n\left(w\right)=\{\begin{array}{l}{\left(1-w\right)}^3/6:n=0\\ \left(3w^3-6w^2+4\right)/6:n=1\\ \left(-3w^3-3w^2+3w+1\right)/6:n=2\\ w^3/6:n=3\end{array}$$

VPR was an important parameter which controlled the approximation quality of B-spline curve for the target image. The fine VPR reduced the approximation error but caused higher computation costs. The coarse VPR increased the approximation error but was more time efficient. The balance between the accuracy and efficiency was investigated in the latter section.

2.3. B-Spline-Based Reconstruction Method

With B-spline-interpolated displacement matrix D, Equations (1) and (3) were reformulated as

$$V=V\left(\mathsf{\Phi }\left(I\right)\right) \mathrm{and} Y=P(V)=P\left(V\left(\mathsf{\Phi }\left(I\right)\right)\right)$$

(11a)

and the objective function of Equation (5) was updated as

$$\min E\left(\mathsf{\Phi }\left(I\right)\right)\begin{array}{cc}& \end{array}s.t.\begin{array}{cc}& \end{array}\Vert P\left(V\left(\mathsf{\Phi }\left(I\right)\right)\right)-Y\Vert =0$$

(11b)

The gradient of the objective function was then expressed as:

$$\nabla f\left(\mathsf{\Phi }\left(I\right)\right)=2P^{-1}\left(P\left(V\left(\mathsf{\Phi }\left(I\right)\right)\right)-Y\right)\nabla V\left(\mathsf{\Phi }\left(I\right)\right)+\mu \ast \nabla E\left(\mathsf{\Phi }\left(I\right)\right)$$

(12)

The gradient $\nabla V$ could be calculated using chain rule as

$$\nabla V\left(\mathsf{\Phi }\left(I\right)\right)=\frac{\partial V}{\partial I}=\frac{\partial V}{\partial \mathsf{\Phi }}\frac{\partial \mathsf{\Phi }}{\partial I}$$

(13)

Combining Equation (2) and Equation (11), $\frac{\partial V}{\partial \mathsf{\Phi }}$ could be calculated as

$$\frac{\partial V_{ijk}}{\partial {\mathsf{\Phi }}_{ijk}^x}=-w_1+w_2$$

(14a)

$$\begin{array}{ll}\frac{\partial V_{ijk}}{\partial {\mathsf{\Phi }}_{ijk}^y}& =\frac{\partial w_1}{\partial {\mathsf{\Phi }}_{ijk}^y}\ast \left(1-\widehat{i}\right)+\frac{\partial w_2}{\partial {\mathsf{\Phi }}_{ijk}^y}\ast \widehat{i}\\ & =\left(-u_1+u_2\right)*\left(1-\widehat{i}\right)+\left(-v_1+v_2\right)*\widehat{i}\end{array}$$

(14b)

$$\begin{array}{cc}\frac{\partial V_{ijk}}{\partial {\mathsf{\Phi }}_{ijk}^z}& =\frac{\partial w_1}{\partial {\mathsf{\Phi }}_{ijk}^z}*(1-\widehat{i})+\frac{\partial w_2}{\partial {\mathsf{\Phi }}_{ijk}^z}*\widehat{i}\\ & =\left(\frac{\partial u_1}{\partial {\mathsf{\Phi }}_{ijk}^z}*\left(1-\widehat{j}\right)+\frac{\partial u_2}{\partial {\mathsf{\Phi }}_{ijk}^z}*\widehat{j}\right)*\left(1-\widehat{i}\right)+\left(\frac{\partial v_1}{\partial {\mathsf{\Phi }}_{ijk}^z}*\left(1-\widehat{j}\right)+\frac{\partial v_2}{\partial {\mathsf{\Phi }}_{ijk}^z}*\widehat{j}\right)*\widehat{i}\\ & =\left(\begin{array}{l}\left(-V{0}_{\overline{i},\overline{j},\overline{k}}+V{0}_{\overline{i},\overline{j},\overline{k}+1}\right)*\left(1-\widehat{j}\right)\\ +\left(-V{0}_{\overline{i},\overline{j}+1,\overline{k}}+V{0}_{\overline{i},\overline{j}+1,\overline{k}+1}\right)*\widehat{j}\end{array}\right)*\left(1-\widehat{i}\right)\hspace{0.17em}\hspace{0.17em}+\left(\begin{array}{l}\left(-V{0}_{\overline{i}+1,\overline{j},\overline{k}}+V{0}_{\overline{i}+1,\overline{j},\overline{k}+1}\right)*\left(1-\widehat{j}\right)\\ +\left(-V{0}_{\overline{i}+1,\overline{j}+1,\overline{k}}+V{0}_{\overline{i}+1,\overline{j}+1,\overline{k}+1}\right)*\widehat{j}\end{array}\right)*\widehat{i}\end{array}$$

(14c)

According to the definition of $\mathsf{\Phi }\left(I\right)$ in Equation (10),$\frac{\partial \mathsf{\Phi }}{\partial I}$ is given by

$$\frac{\partial \mathsf{\Phi }}{\partial I}={\displaystyle \sum _{l=0}^3{\displaystyle \sum _{m=0}^3{\displaystyle \sum _{n=0}^3{\beta }_l\left(u\right){\beta }_m\left(v\right){\beta }_n\left(w\right)}}}$$

(15)

The bending energy function could also be updated as

$$E\left(I\right)={\displaystyle \sum _{k=1}^n{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{m=1}^3\left({\left(\frac{\partial {\mathsf{\Phi }}_{ijk}^x}{\partial x}\right)}^2+{\left(\frac{\partial {\mathsf{\Phi }}_{ijk}^y}{\partial y}\right)}^2+{\left(\frac{\partial {\mathsf{\Phi }}_{ijk}^z}{\partial z}\right)}^2\right)}}}}$$

(16)

and its discrete form is

$$E\left(I\right)={\displaystyle \sum _{k=1}^n{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{m=1}^3\left(\begin{array}{l}{\left({\mathsf{\Phi }}_{i+1,j,k}^x-{\mathsf{\Phi }}_{i,j,k}^x\right)}^2+{\left({\mathsf{\Phi }}_{i,j,k}^x-{\mathsf{\Phi }}_{i-1,j,k}^x\right)}^2\\ +{\left({\mathsf{\Phi }}_{i,j+1,k}^x-{\mathsf{\Phi }}_{i,j,k}^x\right)}^2+\left({\mathsf{\Phi }}_{i,j,k}^x-{\mathsf{\Phi }}_{i,j-1,k}^x\right)\\ +{\left({\mathsf{\Phi }}_{i,j,k+1}^x-{\mathsf{\Phi }}_{i,j,k}^x\right)}^2+{\left({\mathsf{\Phi }}_{i,j,k}^x-{\mathsf{\Phi }}_{i,j,k-1}^x\right)}^2\end{array}\right)}}}}$$

(17)

$\nabla E$ could be calculated using chain rule as

$$\nabla E\left(\mathsf{\Phi }\left(I\right)\right)=\frac{\partial E}{\partial I}=\frac{\partial E}{\partial \mathsf{\Phi }}\frac{\partial \mathsf{\Phi }}{\partial I}.\frac{\partial E}{\partial \mathsf{\Phi }}$$

(18)

where

$$\begin{array}{l}\frac{\partial E}{\partial \mathsf{\Phi }}=-2({\mathsf{\Phi }}_{i+1,j,k}^x-{\mathsf{\Phi }}_{i,j,k}^x)+2({\mathsf{\Phi }}_{i,j,k}^x-{\mathsf{\Phi }}_{i-1,j,k}^x)-2({\mathsf{\Phi }}_{i,j+1,k}^x-{\mathsf{\Phi }}_{i,j,k}^x)\\ \hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+2({\mathsf{\Phi }}_{i,j,k}^x-{\mathsf{\Phi }}_{i,j-1,k}^x)-2({\mathsf{\Phi }}_{i,j,k+1}^x-{\mathsf{\Phi }}_{i,j,k}^x)+2({\mathsf{\Phi }}_{i,j,k}^x-{\mathsf{\Phi }}_{i,j,k-1}^x)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=12{\mathsf{\Phi }}_{i,j,k}^x-2{\mathsf{\Phi }}_{i+1,j,k}^x-2{\mathsf{\Phi }}_{i-1,j,k}^x-2{\mathsf{\Phi }}_{i,j+1,k}^x-2{\mathsf{\Phi }}_{i,j-1,k}^x-2{\mathsf{\Phi }}_{i,j,k+1}^x-2{\mathsf{\Phi }}_{i,j,k-1}^x\end{array}$$

(19)

Combining Equations (13) and (15) into Equation (12), $\nabla f\left(I\right)$ could be finally obtained as

$$\nabla f\left(\mathsf{\Phi }\left(I\right)\right)=\left[2P^{-1}\left(P\left(V\right)-Y\right)\ast \frac{\partial V}{\partial \mathsf{\Phi }}+\mu \ast \frac{\partial E}{\partial \mathsf{\Phi }}\right]\ast {\displaystyle \sum _{l=0}^3{\displaystyle \sum _{m=0}^3{\displaystyle \sum _{n=0}^3{\beta }_l\left(u\right){\beta }_m\left(v\right){\beta }_n\left(w\right)}}}$$

(20)

2.4. Evaluation

The proposed method was implemented in an iterative manner as shown in Figure 1. Initially, the B-spline coefficients and the displacement matrixes were set to zero. The gradient of the objective function was then calculated, and the B-spline coefficients were updated based on the resulting gradient. After certain iterations, the optimization process converged when the value of $\nabla f$ was less than the threshold or the maximum number of iterations reached. The real-value displacement matrixes were interpolated from the optimized B-spline coefficients. Then, the target image was transformed from the prior image via the resulted real-value displacement matrixes. Finally, the approximation errors of B-spline-interpolated displacement matrix and the reconstructed target image were calculated for analysis.

To evaluate the similarity between the reconstructed target image (V^R) and the ground-truth target image (V^T), three intensity-based statistical metrics, normalized cross correlation, mutual information and mean absolute percentage error, were employed. Normalized cross correlation (NCC) was defined as

$$NCC\left(X^R,X^T\right)=\frac{\left(X^R-{\mu }^R\right)\left(X^T-{\mu }^T\right)}{{\sigma }^R{\sigma }^T}$$

(21)

where µ^R and µ^T are means of V^R and V^T and where σ^R and σ^T are standard deviations of V^R and V^T. Mutual information (MI) was defined as

$$MI\left(V^R,V^T\right)={\displaystyle \sum _{x^R\in V^R}^{}{\displaystyle \sum _{x^T\in V^T}^{}p\left(x^R,x^T\right){\log }_2\frac{p\left(x^R,x^T\right)}{p\left(x^R\right)p\left(x^T\right)}}}$$

(22)

where p(x^R) and p(x^T) are the normalized histograms of V^R and V^T and where p(x^R, x^T) is the joint histogram of V^R and V^T. In order to calculate the intensity difference between two images, the mean absolute percentage error (MAPE) was defined as

$$MAPE=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|\frac{x_i^R-x_i^T}{x_i^T}\right|}$$

(23)

where x^R_i and x^T_i are the voxels of images V^R and V^T and where n is the total number of voxels.

The proposed method was examined using a clinical head-and-neck patient. The CBCT image size was 256 × 256 × 256, and the cone-beam projection size was 512 × 384. The full-angle reference CBCT and target CBCT were reconstructed with FDK method using two sets of CB projections that were acquired on two separate days. The CB projections were acquired at every 0.55°. The influence of scan angle (30°, 45°, 60° and 90°) and scan orientation (central angle at 0°, 45° and 90°) were investigated. The accuracy of our method was evaluated by comparing the limited-angle CBCT (target image) to the full-angle CBCT (ground-truth) that was reconstructed with Feldkamp–Davis–Kress (FDK) method. The reconstruction and evaluation codes were built on the open-source platform, Plastimatch, which provides B-spline-based image registration and GPU-based image reconstruction toolkits.

3. Results

The images reconstructed using the proposed and FDK methods on the clinical head-and-neck case is shown in Figure 2. The projections were acquired with a 0° scan orientation and a 60° scan angle. The prior image that was reconstructed using a full-angle scan with the FDK method is shown in Figure 2a. The target image that was reconstructed with the FDK method using a limited-angle scan is shown in Figure 2b. The target image that was reconstructed with the B-spline-based reconstruction method using a limited-angle scan is shown in Figure 2c. As the ground-truth image, the target image that was reconstructed with the FDK method using a full-angle scan is shown in Figure 2d. The similarity between the images in Figure 2c,d was high, which indicated that the proposed method could generate comparable images, using a limited-angle scan, to the ground-truth image, using a full-angle scan.

The images that were reconstructed with the deformation-based reconstruction method as described in Section 2.1 are shown in Figure 3. The projections were acquired with a 0° scan orientation and 30°/45°/60°/90° angles. It was noticed that there was strong noise in the target images that were reconstructed with 30°, 45° and 60° scan angles. The image quality of the reconstructed images became worse as the scan angle decreased. Among the four images, the image that was reconstructed with a 90° scan angle (Figure 3d) showed the best quality and was closer to the ground-truth image (Figure 2d).

The images that were reconstructed with the proposed B-spline-based reconstruction method using a limited-angle scan as described in Section 2.3 are shown in Figure 4. The setting of this reconstruction was the same as the one used in Figure 3. It was noticed that the noise in the images as shown in Figure 4 was greatly suppressed comparing to that in the images in Figure 3. The image quality with the 60° and 90° scan angles was better and was closer to that of the ground-truth image (Figure 2d). The image quality with the 90° scan angle (Figure 4d) was the best but required 50% more projections than that with the 60 ° scan angle (Figure 4c).

The influence of the scan orientation on the image quality was also investigated. The images that were reconstructed with the FDK method and the proposed method with a 0° scan orientation and a 60° scan angle are shown in Figure 2b,c. The images that were reconstructed with both methods with a 45° scan orientation and a 60° scan angle are shown in Figure 5a,b, while the images that were reconstructed with both methods with a 90° scan orientation and a 60° scan angle are shown in Figure 5c,d. The images as shown in Figure 3a–c that were reconstructed with the FDK method with three scan orientations were significantly different. The images as shown in Figure 2c and Figure 5b,d that were reconstructed with the proposed method with three scan orientations were similar to each other. This indicated that the influence of the scan orientation on the image quality of the reconstructed image was minor.

The accuracy of the B-spline-interpolated displacement matrix and the reconstructed target image with respect to the different settings of the VPR and scan angles were investigated. In

Table 1. The accuracy of the B-spline-interpolated displacement matrix.

Scan Angle	VPR = 60			VPR = 30			VPR = 5
	NCC	MI	MAPE	NCC	MI	MAPE	NCC	MI	MAPE
30°	0.62	0.72	0.31	0.83	1.21	0.13	0.88	1.31	0.11
60°	0.76	0.99	0.20	0.89	1.25	0.19	0.93	1.33	0.06
90°	0.85	1.08	0.15	0.94	1.29	0.08	0.97	1.37	0.05

, the final B-spline-interpolated displacement matrix was compared to the ground-truth displacement matrix which was obtained from the deformable transformation between the prior and ground-truth images. For a combination of the scan angle and VPR, three similarity metrics (NCC, MI and MAPE) were calculated. The result showed that, with a smaller VPR and a larger scan angle, the similarity between the resulted and ground-truth displacement matrixes was higher. In

Table 2. The accuracy of the reconstructed target image.

Scan Angle	VPR = 60			VPR = 30			VPR = 5
	NCC	MI	MAPE	NCC	MI	MAPE	NCC	MI	MAPE
30°	0.80	0.91	0.17	0.88	1.31	0.09	0.90	1.33	0.07
60°	0.88	1.22	0.12	0.94	1.51	0.07	0.96	1.55	0.04
90°	0.90	1.28	0.10	0.98	1.58	0.05	0.99	1.65	0.02

, the final reconstructed target image was compared to the ground-truth target image. For a combination of the scan angle and VPR, the same three similarity metrics were calculated. The result showed that, with a smaller VPR and a larger scan angle, the similarity between the reconstructed and ground-truth images was higher.

The running time on the CPU and GPU platforms were compared in

Table 3. Running time on CPU and GPU for different reconstruction tasks.

Tasks	Reconstruction with 30° Scan Angle			Reconstruction with 60° Scan Angle
Iterations	CPU	GPU		CPU	GPU
	Intel Xeon 2609	Nvidia 1080	Nvidia 2080	Intel Xeon 2609	Nvidia 1080	Nvidia 2080
100	11,000s	100s	75s	13,400s	120s	100s
200	19,000s	230s	200s	25,500s	250s	220s

. The inputs were the 2D projections acquired with the 30° or 60° scan angle and the 3D prior image. The dimensions of the 2D projections were 512 × 384, and the dimensions of the 3D image was 256 × 256 × 256. The output was the 3D image in the dimensions of 256 × 256 × 256. The computation of the 2D and 3D image reconstructions and B-spline coefficients were all accomplished on a GPU. Two popular graphics cards, Nvidia Geforce GTX 1080 and Nvidia Geforce GTX 2080, were tested. The computer was equipped with a Intel Xeon 2609 CPU and a 64GB DRR4 memory. It showed that the running time on a GPU was hundred times faster than that on a traditional CPU.

4. Discussion

With a limited-angle scan, the proposed method could reconstruct images comparable to the standard images reconstructed with a full-angle scan. This will significantly reduce the time and dose during multiple data acquisition and would be beneficial for patients under treatment. With the introduction of the B-spline interpolation, the image noise and artifacts were suppressed, and the edge of the anatomical structures was improved. The VPR had a certain effect on the image quality, but it could be minimized with smaller values. With the acceleration of a GPU, the iterative CBCT reconstruction could be mostly accomplished within 10 min. This would make it possible for several clinical scenarios.

There are many factors that could affect the image quality and accuracy of the reconstructed images. Among them, the scan orientation and scan angle are the most important ones. Our experiments showed that the scan orientation had a small influence on the images, but the scan angle had a large influence on the resulting images. When the scan angle was less than 30°, there was strong noise and a larger structure discontinuity in the reconstructed images. A larger scan angle would result in a better image but would take more time for reconstruction computation. Therefore, to balance accuracy and efficiency, a 60° scan angle would be favorable for clinical use.

The introduction of the B-spline interpolation was effective for noise suppression. It could also be feasible in the other types of limited-angle reconstruction methods, such as total-variation-based and optical-flow-based reconstruction methods. In addition, this method could also be used for the reconstruction of 4D-CT and 4D-CBCT with limited-angle projections. In these 4D reconstruction tasks, the 3D volume in one phase could be used as a prior image to predict the 3D volume of another phase. As a large amount of projections is usually required in 4D applications, the reduction of the time and dose resulted by the proposed method on data acquisition would be considerable.

It is also possible to extend this method to the other industrial domains, such as magnetic resonance imaging, radar interferometry, electromagnetism, etc. As these fields differ from medical imaging, it would be important to identify the proper setting and parameters for them. For example, the imaging principle of MRIs would differ from that of CT/CBCTs substantially. The former is generated by detecting the direction of the change of protons in the water, while the latter is generated from the X-ray interaction with the material on its path and finally deposing the residual radiation dose on the detectors. The proper preprocessing and adjustments would be necessary.

5. Conclusions

A new B-spline-based limited-angle reconstruction method was developed and demonstrated its feasibility in reconstructing high-quality images comparable to the standard images. With the introduction of the B-spline interpolation, the image noise was significantly suppressed, and the contour of the anatomical structures was improved. With the GPU acceleration, the limited clinical time of image reconstruction could be greatly alleviated. Potentially, this reconstruction method could be used in the other industrial domains, but proper preprocessing and adjustments would be needed.