John’s Equation-based Consistency Condition and Corrupted Projection Restoration in Circular Trajectory Cone Beam CT

In transmitted X-ray tomography imaging, the acquired projections may be corrupted for various reasons, such as defective detector cells and beam-stop array scatter correction problems. In this study, we derive a consistency condition for cone-beam projections and propose a method to restore lost data in corrupted projections. In particular, the relationship of the geometry parameters in circular trajectory cone-beam computed tomography (CBCT) is utilized to convert an ultra-hyperbolic partial differential equation (PDE) into a second-order PDE. The second-order PDE is then transformed into a first-order ordinary differential equation in the frequency domain. The left side of the equation for the newly derived consistency condition is the projection derivative of the current and adjacent views, whereas the right side is the projection derivative of the geometry parameters. A projection restoration method is established based on the newly derived equation to restore corrupted data in projections in circular trajectory CBCT. The proposed method is tested in beam-stop array scatter correction, metal artifact reduction, and abnormal pixel correction cases to evaluate the performance of the consistency condition and corrupted projection restoration method. Qualitative and quantitative results demonstrate that the present method has considerable potential in restoring lost data in corrupted projections.


Method
John's equation. In 1938, Fritz John derived the ultra-hyperbolic PDE, which can solve the problem of line integrals on a characteristic surface. Let ε and η denote the X-ray source and detector cell, respectively. John's equation 24 can be expressed as follows: i j j i 2 2 where u denotes a set of line integrals of the object function f, which represents a sufficiently differentiable function of an object along the line through ε and η, and can be expressed as R Geometry configuration of circular trajectory CBCT. In this section, we introduce circular trajectory CBCT geometry, which is utilized to convert the ultra-hyperbolic PDE into an ODE. As shown in Fig. 1, o-xyz is a global coordinate system, where o is the rotation isocenter; (α 1 , α 2 ) is the local coordinate system of the flat panel detector, where α 2 is parallel to the z-axis and α 1 is orthogonal to α 2 . The rotation angle of the X-ray source and detector around the x-axis is θ; ρ and d denote source-to-isocenter distance and isocenter-to-detector distance, respectively.

Consistency condition in circular trajectory CBCT.
In this section, we would like to explain further why Patch's consistency condition method cannot be directly utilized in circular trajectory CBCT. To calculate the unmeasured data in the projections acquired via helical CBCT geometry, Patch 27 initially used geometry parameters of a third-generation helical CT system to rewrite the variables in John's equation as follows: where (k 1 , k 2 ) is the corresponding frequency domain of (α 1 , α 2 ), and the superscript * denotes a Fourier transform. The subscript denotes a derivative; for example, u θ represents the derivative of u with respect to θ. Note that Eq. (4) is a first-order PDE for θ and z, and for fixed frequency component (k 1 , k 2 ), it is referred to as a family of first-order ODEs. Although Patch converts the complex John's equation into an ODE, the helical coordinates in Eq. (4) consider z as a numerical variable. The change in variable z along the vertical axis destroys the homogeneity in a helical system 28 ; thus, Patch's equation is difficult to apply in practice. In theory, helical CT can degenerate into circular trajectory CBCT when z = 0. However, Eq. (4) contains a derivative of z; thus, this equation will not work when z is a constant. So, Eq. (4), which was deduced by Patch, could not be utilized in circular trajectory CBCT. Therefore, we have derived a John's equation-based consistency condition (JECC) and the derivation is described in detail in this section. To make John's equation subject to the circular trajectory CBCT geometry, we firstly rewrite ε and η in terms of five parameters, namely, α 1 , α 2 , ρ, d, and θ. Then ε and η are defined by Eqs (5) and (6) can be transformed into the equivalent equations as follows: Combined Eqs (5), (6) with (7), the derivatives of the CBCT parameters with respect to ε and η are listed in Table 1.
According to the derivative of the compound function and the derivatives shown in Table 1, we can get  Figure 1. Geometry configuration of circular trajectory CBCT.  3 2 In order to facilitate writing, define the differential operator L Then John's equation can be transformed as follows: Due to ρ and d are constants, we would like to express ∂/∂r and ∂/∂d using other variables. Therefore two gradient differential equations are derived firstly, where l denotes the line through two points ε and η According to integration by parts, Eq. (18) can be expressed as That is, Similarly Scientific RepoRts | 7: 4920 | DOI:10.1038/s41598-017-05249-5 Eqs (20) and (21) can be rewritten with α 1 , α 2 , ρ, d and θ as follows:  After transposition and combining like terms, the expressions of ∂/∂ρ and ∂/∂d are as follows: With Eqs (24) and (25), Eq. (17) can be expressed as Finally the equation we have derived is The Fourier transform of Eq. (28) is derived with respect to α 1 and α 2 . Let u* represent the Fourier transform of u and (k 1 , k 2 ) be the corresponding frequency form of (α 1 , α 2 ). Consequently, Eq. (28) is transformed into Let U denote the right side of (29), in terms of derivative definition, Eq. (29) can be simply expressed as Eq. (30) is the JECC equation in the frequency domain. Note that the partial derivatives contained in U are consistent with u*(θ), e.g., when u*(θ) denotes the Fourier transform of the projection at angle θ, ∂u*/∂k 2 means ∂u*(θ)/∂k 2 . To avoid any confusion, we add an angle label to symbol U, therefore Eq. (30) can also be written as follows: 2 As shown in Fig. 2, JECC works well when k 2 ≠ 0 in Eq. (30) and it can restore the area in the blue zone, which represents high-frequency components. However, the zero-frequency component represented by the white horizontal line is beyond the capability of JECC when k 2 = 0 in Eq. (30). To address this issue, a corrupted projection u is first interpolated into uʹ via spatial interpolation (SI) in the projection domain, and then uʹ is converted into the frequency domain via fast Fourier transform (FFT) to restore the low-frequency components, which are represented by the purple rectangle. Thus, the proposed method can combine the advantages of JECC and SI to restore corrupted data in the frequency domain. In the following section, we refer to the proposed restoration method as the JECC method. Figure 3 shows the flowchart of the JECC method, or the corrupted projection restoration method. The black dots in the projection represent the pixels that should be restored. Let u(θ) denote the projection at view angle θ, and u(θ + dθ) and u(θ − dθ) be the adjacent projections at angle θ + dθ and θ − dθ, respectively.  In Step 1, horizontal 1D cubic spline interpolation induces few artifacts because ofthe horizontal shift-invariant weighting in the Feldkamp-Davis-Kress (FDK) algorithm 29 . Thus, cubic spline interpolation is utilized to interpolate the projection in each view to obtain the initial projection denoted by u I (θ). The purpose of this step is to restore the low-frequency components of the corrupted projection.

Workflow of corrupted projection restoration.
In Step 2, u I (θ), u(θ + dθ), and u(θ − dθ) are converted into frequency domain via Fourier transform. After this step, the corresponding frequency domains u I* (θ), u*(θ + dθ), and u*(θ − dθ) can be obtained. In Step 3, u*(θ) can be restored by substituting u*(θ + dθ) and u*(θ − dθ) separately into Eqs (31) and (32). Then, u R* (θ), which is the updated frequency data at angle θ, is a weighted result that combines Eqs (31) and (32) via a factor ω. ω is a weighting factor which weights the two corresponding results of (31) and (32). The whole derivation is based on the assumption that u is differentiable, which is not true in reality. So the method need to combine more restoration information at different direction. In the paper, we set ω to 0.5 in order to make two adjacent projections have equal contribution. Subsequently, the initial result u I* (θ) in Step 2 is updated to u R* (θ). The JECC equation only holds when k 2 ≠ 0; thus, u I* (θ) is utilized to extract the low-frequency components when k 2 = 0. In conclusion, this step can be expressed as In Step 4, the inverse Fourier transform of u R* (θ) is utilized to obtain the restoration result by FDK algorithm in the first iteration.
Steps (2)-(4) are repeated until the termination criterion is satisfied. In this study, we set the iteration number as the termination criterion because the restored value remains stable after a certain number of iterations.

Experiments
Moving BSA scatter correction case. In this section, the JECC method is applied to the simulation of a BSA scatter correction case. As shown in Fig. 4(a) and (b), the BSA dimension is 15 × 7 with each blocker shading 5 × 5 pixels, which are denoted by black dots. In this study, BSA is fixed at the same position in odd-numbered views, as shown in Fig. 4(a), whereas it moves to another location in even-numbered views, as shown in Fig. 4(b). In practice, we divide the projections into two categories: the odd-numbered views use mode I and the even-numbered views use mode II as shown in Fig. 4(c). The movement distance l of BSA between the odd-numbered and even-numbered views is 7 pixels. The phantom we tested is a head phantom called FORBILD, which has rich high-frequency details. The dimension of the flat panel detector is 850 × 200, and its cell size is 1 × 1 mm 2 . Both ρ and d are 500 mm. The number of projections is 1080. Each iteration result is compared with the result of the SI method. These defective cells cause abnormal pixels in raw projections. However, these abnormal pixels can be detected in a projection by repeatedly scanning the projection at various exposure levels. In general, the detector manufacturer provides the abnormal pixel detection and correction protocol. In such protocol, the pixel value that corresponds to the defective cell is linearly interpolated using nearby pixels. However, linear interpolation (LI) typically causes streak artifacts, and artifacts caused by abnormal pixels in one of hundreds of 2D slices are difficult to find among hundreds of 2D slices. The abnormal pixel correction protocol is related to corrupted projection restoration; thus, the JECC method is also utilized in this case.
As shown in Fig. 5, we acquired the raw projections using an in-house bench-top CBCT system, in which the CBCT scanning protocol is presented in Table 2. We then simulated the abnormal pixels induced by the defective detector cells, as illustrated in detail in a previous paper 1 . Furthermore, to further verify the capability of the proposed method, we have tested it in a more challenging case with real data acquired by a detector which has a group of corrupted cells.
Metal artifact reduction case. In terms of CBCT in the circular trajectory geometry, raw 2D projections are firstly acquired by rotating over 360°. Then the transection images are reconstructed by FDK algorithm. For segmenting the metal objects from human tissue, the global threshold method is used in each 2D reconstructed image. In this workflow, all voxels above 3000 HU will be considered as metal voxels. After metal segmentation by threshold, the metal regions are recognized. Subsequently, SI method or the proposed method is selected to restore the corrupted data of the acquired projections with the guidance of metal-only measurement data generated by forward projecting metal images. Then, the corrected image can be reconstructed with the corrected projection data. After the reconstruction, the segmented metal objects are finally reinserted to display the implants. In this case, the projection correlation based view interpolation (PC-VI) method proposed by Yan 16 is also employed to compare the capability of the JECC method.

Performance evaluation. Evaluation by error comparison.
For quantitative measurement, we selected the mean absolute error (MAE) to measure the difference between the conventional SI method and the proposed method. MAE is defined as Figure 5. In-house bench-top CBCT system for raw projection acquisition.  where μ(i) denotes the ith pixel value in the reconstructed image, μ t (i) represents the ith pixel value in the reference image, and N is the number of total pixels. Evaluation by noise reduction. For quantitative comparison, we also used signal-to-noise ratio (SNR) to evaluate noise reduction in the images. SNR is defined as Evaluation by image similarity. In this study, the universal quality index (UQI) 30 was utilized to assess the degree of similarity at the region of interest (ROI) in detail. When the ROI is given, the associative mean, variance, and covariance of the ROI can be defined as: where M represents the number of pixels within the ROI. Then, UQI is defined as UQI shows the similarity level between the two images with a value ranging from 0 to 1. A UQI value closer to 1 indicates that the similarity is high between the reconstructed and reference images.

Moving BSA scatter correction results. Visualization-based evaluation. The experimental results
were compared and analyzed through visual inspection. Figure 6(a-d) show the reference image and the images on the 100th slice reconstructed using the SI method, JECC 1st iteration, and JECC 4th iteration, respectively. Figure 6(e-h) show the reference image and the images on the 80th slice reconstructed using the SI method, JECC 1st iteration, and JECC 4th iteration, respectively. The images reconstructed using the JECC method are visually better than those reconstructed using the SI method. Compared with the SI method, the JECC method achieves good performance in cases of LI-induced streak artifacts. We can obtain a visually satisfactory image quality after the JECC 4th iteration.
Profile-based evaluation. Figure 7(a and b) plot the vertical profiles on the 80th slice and the 100th slice, respectively. The profiles of the results of the JECC method are considerably closer to those of the reference image than those of the results of the conventional SI method. The SI method causes more fluctuations around the reference value, and this result is consistent with the appearance of streak artifacts in the image. The results also indicate that the proposed method can achieve better profiles compared with the SI method.
Quantitative evaluation. The difference between the reference image and the images reconstructed using the conventional SI method and the JECC method are quantitatively evaluated by measuring the MAEs of the reconstructed results in different views. The calculated MAEs with the projection views of 135, 270, 360, 540, and 1080 are shown in Table 3. Compared with the SI method, the first iteration of the JECC method results in an appreciable improvement. Moreover, as the iteration number increases, the MAEs decrease. At the fourth iteration, the JECC method achieves the lowest MAE. After four iterations, the proposed method outperforms the traditional correction method (i.e., SI), demonstrating 72.21% reduction in terms of MAE. The MAEs have not been significantly changed by the fourth iteration onward.
Furthermore, SNR is used to quantitatively evaluate the noise reduction of the present method. The results are presented in Table 4. From 135 views to the 1080 views, the SNR gains of the JECC method compared with that of the SI method are between 7 dB and 3.6 dB. Compared with the SI method, the present method yields higher SNR. This result demonstrates that the present method can achieve noticeable gains in terms of noise and artifact suppression.
To evaluate image similarity between the reconstructed results and the reference image, we selected the ROI marked with a red dashed square in Fig. 6(a) to calculate the UQI scores. The corresponding UQI scores are    shown in Fig. 8. In all the views, the present method yields higher UQI scores (over 0.9), thereby implying that it outperforms the traditional method in terms of UQI by over 0.24 on average.
Abnormal pixel correction results. Visualization-based evaluation. Figure 9 shows the reconstructed images without correction, corrected with SI method and the JECC method on the 190th slice and 220th slice. As shown in the first column of Fig. 9, the abnormal pixels induce ring artifacts. As shown in the second column, the conventional SI method can efficiently remove the ring artifacts. A satisfactory image quality can be acquired using the JECC method after four iterations. The zoomed ROIs in Fig. 9 show that several streak artifacts, which are indicated by the red arrows, are still present. These streak artifacts degrade image quality severely. However, such artifacts are mitigated in the images in the third and fourth columns, and these artifacts even disappear visually in the images in the fifth column, which shows the results obtained after four iterations of the JECC method. Figure 10 displays one more challenging case with ring artifacts. As shown in the first column, uncorrected images suffer from severe bending artifacts caused by a group of blocked pixels on detector. The ring artifacts are so wide that they degrade image diagnosis severely. The second column represents the results of SI method. SI  . Image without correction and images reconstructed using the SI method and the JECC method. The first row presents the images on the 190th slice, and the third row presents the images on the 220th. The second and fourth rows present the magnified ROIs marked with yellow squares in the images in the first and third rows. method can reduce dark ring artifacts at the cost of newly introduced streak artifacts. SI method estimates the missing data by utilizing nearby uncorrupted data, therefore the interpolated data deviates from original true values in the dark shadow. When the corrupted regions are large, the severe deviations will unavoidably lead to new streak artifacts. However, due to the complete utilization of the information of adjacent projections, the proposed method outperforms SI method in both artifact suppression and bone structure preservation, which can be observed in the zoomed ROIs in the Fig. 10.
Profile-based evaluation. Figure 11(a and b) plot the horizontal profiles on the 190th slice and 220th slice. The profiles of the results of the JECC method are considerably closer to those of the reference image than those of the results of the SI method. The profile of the image without any correction does not match the reference well because of the ring artifacts; this result partially proves that ring artifacts seriously degrade the quality of images. The use of the SI method allows for the improvement of image quality; however, the fluctuation in the profile demands for further improvement. By contrast, the JECC method achieves high image quality, as depicted by the blue dashed line in Fig. 11(a and b).
Quantitative evaluation. Table 5  In addition, the SNRs of the results are presented in Table 6. Compared with the SI method and the uncorrected results, the present method yields relatively higher SNRs. From 135 views to 1080 views, the SNRs of the JECC method increased by 5.44 dB to 1.99 dB compared with the conventional SI method. Thus, the present method can achieve noticeable SNR gains over the conventional SI method.
In this case, we selected the ROI marked with a red dashed square in Fig. 9 to calculate the UQI scores. The corresponding UQI scores are shown in Fig. 12. In all the cases, the present method yields UQI scores that are 0.51 higher than those of the uncorrected image and 0.19 higher than those of the SI method. These results partially demonstrate that the JECC method is a better solution for restoring local detailed information.
Metal artifact reduction results. Visualization-based evaluation. To evaluate the capability of JECC method, one anthropomorphic male pelvic phantom (CIRS Inc., Norfolk, VA, USA) was used to simulate bilateral hip prostheses. The reconstructed images without correction, corrected with SI method, PC-VI, and JECC method are shown in Fig. 13, respectively. Transverse and coronal views indicate that the uncorrected image suffers from severe beam hardening artifacts between bilateral metal prostheses. In sagittal view, there are severe streak artifacts reducing image quality. All SI method, PC-VI, and JECC method can reduce those metal artifacts to varying degrees. However, SI method only removes some streak artifacts at the cost of the loss of bone structures surrounding the metal prostheses. As for PC-VI, though the bone structures lost in the result of SI method are preserved, the existence of residual streak artifacts observed in the transverse and coronal views contaminates the tissues between bone structures as seen from the sagittal view. Compared with above two methods, JECC method achieves better performance in both the preservation of bone edge structures and the suppression of metal artifacts in the transverse, sagittal, or coronal views. Figure 11. Image profiles of the results in Fig. 9. (a) Profiles across the 315th to 375th rows in the 145th column of the results indicated by a blue line in the top row in Fig. 9; (b) profiles across the 345th to 385th rows in the 170th column of the results indicated by a blue line in the third row in Fig. 9.   Table 6. SNR measurements from 135 views to 1080 views cases. Figure 14(a and b) show the profiles indicated by the blue dashed line in transverse and coronal views, respectively. Figure 14 indicates that the profile line of JECC is closer to that of the reference than those of SI and PC-VI. More fluctuations around the reference value imply the loss of fine bone edge structures and the appearance of metal artifacts. These profiles also demonstrate that JECC can gain better results compared with other two methods.

MAE (10
Furthermore, the performance of JECC method was further demonstrated by clinical measurement data. As shown in Fig. 15, two different kinds of clinical cases were tested in this work. The first row is a head case with a brain stimulator, and the third row is a case with multiple dental fillings which is referred to as the most challenging in the field of metal artifact reduction due to dental dense characteristic and irregular edge shape. In the head case, it is obviously that these radial streak artifacts tangent to metal stimulator severely degrade the uncorrected image quality. The correction result of SI is even worse than the uncorrected image due to failure to reduce original artifacts and some newly introduced artifacts. However, PC-VI and JECC method can get visually satisfactory image quality as seen from zoomed ROIs in the second row in Fig. 15. The correction results   for dental case which is classified as the worst situation are shown in the bottom two rows. The beam hardening artifacts between four dental fillings and severe streak artifacts protruding from the corner of each filling are significantly present in the uncorrected image. SI method removes dark beam hardening artifacts considerably but the tooth edges are blurred due to the inaccurate information restored by interpolation. Some bright artifacts in the uncorrected image disappear after the correction with PC-VI, but some streak artifacts still surround dental fillings. However, there are mildest artifacts remaining in the result of JECC method. Furthermore, the dental edges blurred by SI method are preserved significantly. All in all, JECC method can render an overall improvement of image quality.

Discussion
In this study, we derive a new consistency condition from John's equation for the circular trajectory CBCT. As shown in the scatter correction, abnormal pixel correction, and metal artifact reduction experiments, the proposed restoration method can obtain a visually satisfactory image quality. The quantitative results of the present method are significantly better than those of the conventional method for different slices and view numbers. This outcome can be ascribed to the complete utilization by the present method of the information of adjacent projections. By contrast, the conventional SI method only utilizes the information within the current projection to restore corrupted data in each projection.
In addition, we wish to discuss other relevant issues. As shown in Eq. (30), the adjacent projections, i.e., u*(θ + dθ) and u*(θ − dθ), can be separately substituted into JECC to calculate the restored projection u*(θ). To explore the Fourier properties of JECC, we compare the changes before and after utilizing JECC within the same image in the frequency domain. Figure 16(a and b) show the frequency domain of the image with the SI method and the proposed method, respectively. Figure 16(c) presents the difference by subtracting Fig. 16(a) from 16(b). Figure 16(c) shows that JECC can restore more frequency information of the corrupted projections, particularly their high-frequency components.
Another crucial point we should pay attention to is that the John's equation related deviation is based on the assumption that the image u is reasonably differentiable. In fact, this is inaccurate, e.g., due to the sharp jump at Figure 15. Two clinical cases are displayed here. The first row is the head case with a brain stimulator and the second row presents the corresponding magnified ROIs which are the blue dashed squares containing the metal implant, the bottom two rows present the dental case with three dental fillings on the back teeth and one dental filling on the front tooth. boundary of tissues. However, in view of the complete utilization of the information of adjacent projections, this drawback of John's equation relevant consistency condition can be ignored.

Conclusion
In this study, a new consistency condition from John's equation is derived to restore corrupted projections in circular trajectory CBCT. The proposed method is tested in moving BSA scatter correction, metal artifact reduction, and abnormal pixel correction. The performance results verify that the proposed JECC method is useful and effective in restoring lost data of corrupted projections in circular trajectory CBCT.