Rigid motion artifact reduction in CT using frequency domain analysis

Abstract

BACKGROUND:

It is often unrealistic to assume that the subject remains stationary during a computed tomography (CT) imaging scan. A patient rigid motion can be decomposed into a translation and a rotation around an origin. How to minimize the motion impact on image quality is important.

OBJECTIVE:

To eliminate artifacts caused by patient rigid motion during a CT scan, this study investigated a new method based on frequency domain analysis to estimate and compensate motion impact.

METHODS:

Motion parameters was first determined by the magnitude correlation of projections in frequency domain. Then, the estimated parameters were applied to compensate for the motion effects in the reconstruction process. Finally, this method was extended to helical CT.

RESULTS:

In fan-beam CT experiments, the simulation results showed that the proposed method was more accurate and faster on the performance of motion estimation than using Helgason-Ludwig consistency condition method (HLCC). Furthermore, the reconstructed images on both simulated and human head experiments indicated that the proposed method yielded superior results in artifact reduction.

CONCLUSIONS:

The proposed method is a new tool for patient motion compensation, with a potential for practical application. It is not only applicable to motion correction in fan-beam CT imaging, but also to helical CT.

Keywords

CT fan-beam frequency domain motion artifacts

1 Introduction

The artifact reduction is one of the most important problems in CT reconstruction. Many different factors may lead to increased artifacts, such as the limitations of CT equipment and reconstruction algorithm and the motion of the patients. And one of the major sources of image artifact in CT is patient motion, which is the most difficult to correct [1 –4]. Wagner et al. reported all the patients were subject to head motion during CT even though they were instructed not to do so before scanning. The averages of the mean translation excursion (MTE) and the mean rotation excursion (MRE) were found to be 2.5 mm and 1.06°, respectively [5]. Unfortunately, these are only average values and the motion range is much greater in some cases. If a moving object is reconstructed without motion correction, eminent streaks or even “double images” could be produced [6 –8]. Due to the relatively high radiation dose associated with CT scanning that can provide injury to humans, it is undesirable to repeat the scan if motion occurs. Hence it is very important to develop algorithms and techniques to compensate motion artifacts in CT [9 –12].

Rigid motion can be divided into translational motion and rotational motion, which are coexistent. However, correction method is different according to different types of motion [13 –15]. Several approaches for motion artifact reduction in CT scan have been investigated. Generally speaking, they can be classified into three categories. The first is to minimize motion effects, Li and his colleagues improved temporal resolution by increasing rotation speed of the gantry. This technique is restricted by the scanner’s mechanical limits, and has the disadvantage of reducing the image signal-to-noise (SNR) ratio [16, 17]. Multi-acquisition or overlapping half-scan methods were also carried out by discarding motion corrupted image datasets and reconstructing an image from the remaining non-motion corrupted projections [18 –21]. This method is not effective for constant motion, and increases patient dose. The second is to correct the motion using external motion monitoring devices. Their main drawback is that they require additional hardware [22 –24]. The third is to estimate the motion parameters and compensate for the motion effect in the image reconstruction [25 –27]. The method was proposed in [28] to estimate the motion parameters only from sinogram based on the Helgason-Ludwig consistency condition method (HLCC). However, this method needs a lot of computation due to dependence on a searching algorithm, and has a disadvantage of poor estimation accuracy when the motion curve is complex. The authors of [29] proposed a method to estimate and correct the motion in CBCT, but only for CBCT. Sun et al. [30] proposed a method to iteratively estimate and compensate the motion during the reconstruction. But the iterative processing requires long computation. The third category has still great research potential.

In this paper, the proposed method is in the third category to estimate the motion parameters based on frequency domain analysis. First, the relationships are derived between the rigid motion parameters and projections in fan-beam. The rigid motion can be estimated by the magnitude correlation of projections in frequency domain. It realizes the estimation of rotation and translation simultaneously. Then, the estimated parameters can be introduced into an improved fan-beam filtered back-projection reconstruction algorithm to correct motion artifacts. Finally, this method can be extended to helical CT. The experimental results demonstrate the efficacy of the proposed method.

2 Materials and methods

2.1 Motion estimation in fan-beam

2.1.1 Fan-Beam projection with general in-plane motion

Let us denote a bounded object function in the x-y plane by f (x, y), which is defined on a compact support R², and the corresponding fan-beam projections can be denoted by g_F (β, γ). As shown in Fig. 1, the object is assumed to be rigid with a displacement d = (d_x (t) , d_y (t)) and a rotation $\bar{θ} = θ + φ_{z} (t)$ at time t.

Fig.1

Fan-beam geometry and variables.

According to the geometric relationship defined in Fig. 1, we have

$\bar{γ} = arcsin (sin γ - \frac{d_{x} (t) cos θ + d_{y} (t) sin θ}{R})$ (1) $\bar{β} = \bar{θ} - \bar{γ} = \bar{θ} - arcsin (sin γ - \frac{d_{x} (t) cos θ + d_{y} (t) sin θ}{R})$ (2) $θ = β + γ$ (3) where the parameter β (β ∈ [0, 2π)) is the fan-beam projection view angle, γ is the angle between any fan-beam and the center of ray source, and R is the source path radius.

Assuming γ is a small value, we can make use of the replacement of same value infinitesimal to simplify procedure and solve the problem easily. Hence, Equations (1) and (2) can be simplified as $\bar{γ} = γ - \frac{d_{x} (t) cos θ + d_{y} (t) sin θ}{R}$ (4) $\bar{β} = \bar{θ} - \bar{γ} = β + \frac{d_{x} (t) cos θ + d_{y} (t) sin θ}{R} + φ_{z} (t)$ (5)

Finally, we obtain the projections for a rigid object in the fan-beam geometry $g_{F} (\bar{β}, \bar{γ}) = g_{F} (β + \frac{d_{x} (t) cos θ + d_{y} (t) sin θ}{R} + φ_{z} (t), γ - \frac{d_{x} (t) cos θ + d_{y} (t) sin θ}{R})$ (6) where the subscript F indicates the fan-beam geometry.

2.1.2 Frequency domain characteristics

According to the Fourier shift property, for any projection view angle $\bar{β}$ , the 1D Fourier transform of the projections $g_{F} (\bar{β}, \bar{γ})$ can be defined as [31]. $F_{1} [g_{F} (\bar{β}, \bar{γ})] = e^{- 2 π i ω α (t)} F_{1} [g_{F} (\bar{β}, γ)]$ (7) where F₁ [•] represents the 1D Fourier transform, and $α (t) = \frac{d_{x} (t) cos θ + d_{y} (t) sin θ}{R}$ .

When the X-ray source rotates a circle, each projection can be detected twice by detector $g_{F} (β, γ) = g_{F} (β + π + 2 γ, - γ)$ (8)

Thus, for each projection view angle β₁ ∈ [0, π - 2γ_max), γ_max is the maximum angle of γ_, there must be a projection view angle β₂ ∈ [π, 2π - 2γ_max) which satisfies the equation $| F_{1} [g_{F} (β_{1}, γ)] | = | F_{1}^{*} [g_{F} (β_{2} + 2 γ, γ)] |$ equality. We define the g_F (β₁, γ) and g_F (β₂ + 2γ, γ) as a pair of conjugate projection series. Considering the object with a rigid, we have $| F_{1} [g_{F} ({\bar{β}}_{1}, \bar{γ} (t_{1}))] | = | F_{1}^{*} [g_{F} ({\bar{β}}_{2} + 2 \bar{γ} (t_{1}), \bar{γ} (t_{2}))] |$ (9) where * denotes the complex conjugate, t₁ and t₂ represent the corresponding scanning times for ${\bar{β}}_{1}$ and ${\bar{β}}_{2}$ .

To ensure the existence of the conjugate projection series for a rigid object, we have the view angle ${\bar{β}}_{1} \in [0, π - 2 γ_{max} - μ)$ and ${\bar{β}}_{2} = π + {\bar{β}}_{1} \in [π, 2 π - 2 γ_{max} - μ)$ (μ is the maximum rotating angle of a rigid object).

2.1.3 Motion estimation

Assuming the X-ray source rotates at a constant velocity ω = 2πrad/ - s along a circle, and there are n fan-beam projection view angles β (β ∈ [0, 2π)) and p fan-beam projections $g_{F} (\bar{β}, \bar{γ})$ per projection view angle. The magnitude correlation of projection data in frequency domain is defined as $r = \frac{\sum_{λ = 1}^{p} (x_{λ} - \bar{x}) \cdot (y_{λ} - \bar{y})}{\sqrt{\sum_{λ = 1}^{p} {(x_{λ} - \bar{x})}^{2} \cdot \sum_{λ = 1}^{p} {(y_{λ} - \bar{y})}^{2}}}$ (10) where $x_{λ} = | F_{1} {g_{F} [{\bar{β}}_{1}, \bar{γ} {(t_{1})}_{λ}]} |$ , $y_{λ} = | F_{1}^{*} {g_{F} [{\bar{β}}_{2} + 2 \bar{γ} {(t_{1})}_{λ}, \bar{γ} {(t_{2})}_{λ}]} |$ , $\bar{x}$ and $\bar{y}$ are the average for x_λ and y_λ.

For any x_λ, we can find a y_λ corresponding to ${\hat{\bar{β}}}_{2}$ within ${\bar{β}}_{2}$ to obtain the maximum of r. To reduce computation the search can be narrowed to ${\bar{β}}_{2} \in [{\bar{β}}_{1} + π - μ, {\bar{β}}_{1} + π + μ)$ . There are m pairs of conjugate projection series $g_{F} ({\bar{β}}_{1}, \bar{γ} (t_{1}))$ and $g_{F} ({\hat{\bar{β}}}_{2} + 2 \bar{γ} (t_{1}), \bar{γ} (t_{2}))$ . The angular differences of ${\bar{β}}_{1}$ and ${\hat{\bar{β}}}_{2}$ can be determined as $Δ β + π = {\hat{\bar{β}}}_{2} - {\bar{β}}_{1} = π + φ_{z} (t_{2}) - φ_{z} (t_{1}) + α (t_{2}) - α (t_{1})$ (11)

Thus, we can achieve the rigid motion view-by-view. If we assume that φ_z (t), d_x (t) and d_y (t) are smooth and continuous, we can approximate φ_z (t), d_x (t) and d_y (t) with finite basis functions φ₀ (t) , φ₁ (t) , …, φ_N (t), namely ${\begin{matrix} φ_{z} (t) = p_{10} φ_{0} (t) + p_{11} φ_{1} (t) + \dots + p_{1 N_{1}} φ_{N_{1}} (t) \\ d_{x} (t) = p_{20} φ_{0} (t) + p_{21} φ_{1} (t) + \dots + p_{2 N_{2}} φ_{N_{2}} (t) \\ d_{y} (t) = p_{30} φ_{0} (t) + p_{31} φ_{1} (t) + \dots + p_{3 N_{2}} φ_{N_{2}} (t) \end{matrix}$ (12)

Without loss of generality, we select φ_{N
₁} (t) = t^{N
₁} and φ_{N
₂} (t) = t^{N
₂}. Considering the initial time t = 0, there will be φ_z (0) = 0, d_x (0) = 0 and d_y (0) = 0. Hence, Equation (12) can be simplified as ${\begin{matrix} φ_{z} (t) = p_{11} t + p_{12} t^{2} + \dots + p_{1 N_{1}} t^{N_{1}} \\ d_{x} (t) = p_{21} t + p_{22} t^{2} + \dots + p_{2 N_{2}} t^{N_{2}} \\ d_{y} (t) = p_{31} t + p_{32} t^{2} + \dots + p_{3 N_{2}} t^{N_{2}} \end{matrix}$ (13) where p_1i, p_2j and p_3j (i = 1, 2, …, N₁, j = 1, 2, …, N₂) are the rigid motion parameters. Hence, the rigid motion estimation are transformed into determining the unknown motion parameters p_1i, p_2j and p_3j. Substituting Equation (13) into (11), we have ${\begin{matrix} p_{11} (t_{21 -} t_{11}) + \dots + p_{1 N_{1}} (t_{21}^{N_{1}} - t_{11}^{N_{1}}) + \\ p_{21} a_{11} + p_{31} b_{11} + \dots + p_{2 N_{2}} a_{N_{2} 1} + p_{3 N_{2}} b_{N_{2} 1} = Δ β_{1} \\ ⋮ \\ p_{11} (t_{2 m -} t_{1 m}) + \dots + p_{1 N_{1}} (t_{2 m}^{N_{1}} - t_{1 m}^{N_{1}}) + \\ p_{21} a_{1 m} + p_{31} b_{1 m} + \dots + p_{2 N_{2}} a_{N_{2} m} + p_{3 N_{2}} b_{N_{2} m} = Δ β_{m} \end{matrix}$ (14) where a_1m = [cos (β_2m) t_2m - cos (β_1m) t_1m]/ - R, $a_{N_{2} m} = [cos (β_{2 m}) t_{2 m}^{N_{2}} - cos (β_{1 m}) t_{1 m}^{N_{2}}] / - R$ , b_1m = [sin (β_2m) t_2m - sin (β_1m) t_1m]/ - R, $b_{N_{2} m} = [sin (β_{2 m}) t_{2 m}^{N_{2}} - sin (β_{1 m}) t_{1 m}^{N_{2}}] / - R$ , and the corresponding scanning times are $t_{1} = \frac{{\bar{β}}_{1}}{2 π} \in [t_{11}, t_{12}, \dots, t_{1 m}]$ and $t_{2} = \frac{{\hat{\bar{β}}}_{2}}{2 π} \in [t_{21}, t_{22}, \dots, t_{2 m}]$ .

There are N₁ + 2N₂ unknowns that are p_1i, p_2j and p_3j, and m independent equations in Equation (14). For actual applications, m is much larger than N₁ + 2N₂. Therefore, we can determine the unknown variables as follows $P = {(H^{T} H)}^{- 1} (H^{T} G)$ (15) $G = {(Δ β_{1}, Δ β_{2}, \dots, Δ β_{m})}^{T}$ (16) $P = {(p_{11}, \dots, p_{1 N_{1}}, p_{21}, \dots, p_{2 N_{2}}, p_{31}, \dots, p_{3 N_{2}})}^{T}$ (17) $H = [\begin{matrix} t_{21} - t_{11} & \dots & t_{21}^{N_{1}} - t_{11}^{N_{1}} & a_{11} & \dots & a_{N_{2} 1} & b_{11} & \dots & b_{N_{2} 1} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ t_{2 m} - t_{1 m} & \dots & t_{2 m}^{N_{1}} - t_{1 m}^{N_{1}} & a_{1 m} & \dots & a_{N_{2} m} & b_{1 m} & \dots & b_{N_{2} m} \end{matrix}]$ (18)

Equations (15–18) are the main results of the proposed motion estimation method.

2.1.4 Motion compensation

Once both the rotational and translational parameters are estimated, image reconstruction can be performed to compensate for the motion effects. We can introduce the estimated parameters into fan-beam filtered back-projection reconstruction algorithm to correct motion artifacts. According to the Fourier shift property we compensate for the motion ${\hat{g}}_{F} (β, γ) = F_{2}^{- 1} {e^{- 2 π i [ξ ({\hat{φ}}_{z} (t) + \hat{α} (t)) - η \hat{α} (t)]} F_{2} [g_{F} (\bar{β}, \bar{γ})]}$ (19) where $\hat{α} (t) = \frac{{\hat{d}}_{x} (t) cos θ + {\hat{d}}_{y} (t) sin θ}{R}$ , ${\hat{φ}}_{z} (t)$ , ${\hat{d}}_{x} (t)$ and ${\hat{d}}_{y} (t)$ are the estimated values of the rotation and translation, F₂ [•] represents the 2D Fourier transform, $F_{2}^{- 1} [•]$ represents the 2D inverse Fourier transform, and $t = \frac{\bar{β}}{2 π}$ .

Finally, we obtain the corrected image by fan-beam filtered back-projection reconstruction algorithm. The flow of motion compensation is shown in Fig. 2.

Fig.2

Motion artifact reduction flow.

2.2 Motion estimation in helical CT

Let us consider a single-slice helical CT projections can be denoted by g (β, γ, z), where β denotes the angular position of the source, γ denotes the angle between the projection line and the center line of the fan-beam, and z denotes the longitudinal position of the slice being projected [32]. We assume the projections for a rigid object (d_x, d_y, d_z, φ_x, φ_y, φ_z), namely $g (\bar{β}, \bar{γ}, \bar{z}) = g (β - Δ β, γ - Δ γ, z - Δ z)$ (20) where Δβ, Δγ and Δz can be determined by d_x, d_y, d_z, φ_x, φ_y, φ_z.

According to the Fourier shift property, we take a 2D FFT with respect to $\bar{γ}$ and $\bar{z}$ at each $\bar{β}$ , namely $F_{2} [g (\bar{β}, \bar{γ}, \bar{z})] = e^{- 2 π i (ξ Δ γ + η Δ z)} F_{2} [g (\bar{β}, γ, z)]$ (21)

In the continuous case, the fan-beam redundancy is described by the relationship $g_{F} (β, γ, z) = g_{F} (β + π + 2 γ, - γ, z)$ (22)

Thus, for each projection view angle β₁ ∈ [0, π - 2γ_max), γ_max is the maximum angle of γ, there must be a projection view angle β₂ ∈ [π, 2π - 2γ_max) which satisfies the equation F₂ [g_F (β₁, γ, z)] = F₂ [g_F (β₂ + 2γ, - γ, z)] equality. Considering an object with a rigid, we have

$F_{2} [g_{F} ({\bar{β}}_{1}, \bar{γ} (t_{1}), {\bar{z}}_{1} (t_{1}))] = e^{- 2 π i {ξ [Δ γ (t_{2}) - Δ γ (t_{1})] + η [Δ z (t_{1}) - Δ z (t_{2})]}} F_{2} [g_{F} ({\bar{β}}_{2} + 2 \bar{γ} (t_{1}), - \bar{γ} (t_{2}), {\bar{z}}_{2} (t_{2}))]$ (23)

where t₁ and t₂ represent the corresponding scanning times for ${\bar{z}}_{1}$ and ${\bar{z}}_{2}$ .

Therefore, the rigid motion (d_x, d_y, d_z, φ_x, φ_y, φ_z) can be estimated slice-by-slice by magnitude correlation and phase correlation of projection data in frequency domain [33 –36]. Finally, according to the Fourier shift property we compensate for the motion slice-by-slice, and the corrected image can be obtained by fan-beam filtered back-projection reconstruction algorithm [25, 32].

3 Results and discussion

To demonstrate the feasibility of the proposed method, we performed simulation on a modified Shepp-Logan head phantom [37] and human head CT images as shown in Fig. 3. The additional tiny circles in the modified Shepp-Logan head phantom are used to simulate the small tissue as shown in Fig. 3(a).

Fig.3

Head phantoms for simulation. (a) is the modified Shepp-Logan head phantom; (b) and (c) Image1 and Image2 are the human head CT images.

In experiments we used fan-beam sinogram. The phantom was placed inside a disk of radius 128 mm, and the initial scanning time was t = 0. The X-ray source was rotated at a constant angular speed ω = 2πrad/ - s along a circle of radius R = 500 mm. There were 720 projection view angles and 1000 detector cells distributed along the collinear detector array. The size of detector aperture was 0.3 mm and γ_max = 0.174rad/ - s. Poisson noise was added to the projection data such that 10⁶ photons were emitted from the X-ray source towards each detector cell [28]. The basis functions were set to φ_N (t) = t^N and N₁ = N₂ = 3 in our preliminary study.

The mean translation excursion (MTE) and mean rotation excursion (MRE) are used to evaluate the motion range of an object. Additionally, we use the compensative mean translation excursion (CMTE), relative mean translation excursion (RMTE), compensative mean rotation excursion (CMRE) and relative mean rotation excursion (RMRE) to evaluate the accuracy of motion estimation [28].

3.1 Shepp-Logan phantom experiments

First, we conducted simulation experiments using the modified Shepp-Logan phantom. The motion patterns in three typical cases are listed in Table 1. Group I simulated the ideal motionless case, Groups II and III simulated the complex motion patterns that are cosine and exponential form. The tests were conducted in 30 groups to evaluate the performance of our proposed motion estimation method. In Table 2, the motion was estimated based on the two motion models: the translational motion model (TMM) that ignores the rotation and the general motion model (GMM) in this paper. In Group I, both the GMM and TMM do not introduce deviation. In Groups II and III, the GMM is superior to the TMM since the rotation has been compensated. Therefore, the rotation cannot be ignored.

Table 1
Motion patterns

I II III

φ (t) 0 2 (1 - e^t) 5.73 (cos (πt/  - 2) - 1)

d_x (t) 0 3 (1 - cos (πt/  - 2)) 4 (1 - e^t)

d_y (t) 0 3 sin (πt) 5 (cos (πt/  - 2) - 1)

MRE/° 0 1.437 2.0833

MTE/mm 0 2.4163 3.4164

	II	III
φ (t)	2 (1 - e^t)	5.73 (cos (πt/ - 2) - 1)
d_x (t)	3 (1 - cos (πt/ - 2))	4 (1 - e^t)
d_y (t)	3 sin (πt)	5 (cos (πt/ - 2) - 1)
MRE/°	1.437	2.0833
MTE/mm	2.4163	3.4164

Table 2

Motion estimation results of the TMM and GMM

		I	II	III
TMM	CMRE	0	1.437	2.0833
	RMRE	∞	100%	100%
	CMTE	0	0.5550	0.9808
	RMTE	∞	22.97%	28.71%
GMM	CMRE	0	0.1421	0.2561
	RMRE	∞	9.89%	12.29%
	CMTE	0	0.3997	0.7912
	RMTE	∞	16.54%	23.16%

Table 3 shows the comparison of results between the proposed method based on the GMM and the HLCC method in [28]. In Group I, the proposed method does not introduce deviation. In Groups II and III, the proposed method is superior to the HLCC method. The rotation estimation accuracy of the proposed method seems slightly higher than the HLCC method, while the translation estimation accuracy of the proposed method is higher than the HLCC method by about 50%. Comparing the running times of the three groups, we find that the proposed method is much lower than the HLCC method.

Table 3

Motion estimation results and running times of Shepp-Logan phantoms

		I	II	III
CMRE	GMM	0	0.1421	0.2561
	HLCC	0.0067	0.1599	0.3692
RMRE	GMM	∞	9.89%	12.29%
	HLCC	∞	11.13%	17.72%
CMTE	GMM	0	0.3997	0.7912
	HLCC	0.1074	0.7217	1.8089
RMTE	GMM	∞	16.54%	23.16%
	HLCC	∞	29.87%	52.95%
Time/s	GMM	0.6674	0.6731	0.6801
	HLCC	164.8225	185.5083	187.1339

The estimated results for Group III are plotted in Fig. 4 to demonstrate the advantage of our proposed motion estimation method. The motion curves for the left, middle and right columns correspond to φ_z (t), d_x (t) and d_y (t). Clearly, the estimated motion curves of the proposed method are closer to the real motion than the HLCC method. The results indicate that the proposed method has a better performance on motion estimation accuracy and real-time.

Fig.4

Real motion and estimated motion results for Group III.

All the reconstructed images are shown in Fig. 5. The motion parameters for the three rows correspond to Groups I to III in Table 1. First column images were reconstructed from sinogram with rigid motion and noise by the traditional FBP algorithm with no motion compensation. Second and third column images were reconstructed from sinogram with rigid motion and noise by the improved FBP algorithm with the GMM-based compensation and the HLCC-based compensation. Clearly, both two methods have a good performance on artifact reduction in Group I. In Groups II and III, the proposed method is superior to the HLCC method in image quality of circle regions where tiny structures existed since the motion estimation accuracy of the proposed method is higher than the HLCC method.

Fig.5

Motion artifact reduction in Shepp-Logan phantom for Groups I to III. First column images are reconstructed by the traditional FBP algorithm with no motion compensation. Second and third columns are reconstructed by the improved FBP algorithm with the GMM-based compensation and the HLCC-based compensation. The motion parameters for the three rows correspond to Groups I to III.

To further illustrate the performance of the proposed method, we use the root mean square error (RMSE), Pearson correlation coefficient (CC) and mean structural similarity index (MSSIM) to judge disparity between the reconstructed images and the ideal phantom [24, 38]. Table 4 shows the comparison of artifact reduction effect for Group III.

Table 4

Comparison of artifact reduction performance on Shepp-Logan phantom for Group III

	RMSE	CC	MSSIM
Uncorrected	0.1402	0.6741	0.3672
GMM	0.0686	0.9577	0.7652
HLCC	0.1162	0.7496	0.5151

In Table 4, we can see that each RMSE of two methods is less than uncorrected, and the RMSE of the proposed method is the smallest and reduce about 41% than the HLCC method. Each CC and MSSIM of two methods is closer to 1 than uncorrected, and the CC and MSSIM of the proposed method are the greatest and increase about 28% and 49% respectively than the HLCC method. The results indicate that the proposed method has a better performance on motion artifact reduction.

3.2 Human head experiments

The following simulation experiments were taken using the human head as shown in Figs. 3(b) and 3(c) to verify the artifact reduction effect of the human body. The head CT images were acquired from a GE Hi-Speed multi-slice CT scanner. 720 views were uniformly selected over a 360° range under 150 mA X-ray tube current and 120 kVp. All of the reconstructed images are 512×512 pixels. The real head motion obtained by Kim et al. [24] and estimated motion results are plotted in Fig. 6. Figure 6(a) shows the estimated results of a head case as in Fig. 3(b), and Fig. 6(b) shows the estimated results of a second head case as in Fig. 3(c).

Fig.6

Real head motion and estimated results. (a) are the results of the proposed method and HLCC method for Image1; (b) are the results of the proposed method and HLCC method for Image2.

To quantitatively evaluate the accuracy of the real head motion estimation results between the proposed method and the HLCC method, we use the RMRE and RMTE to quantify the results, as shown in Table 5. We can observe that the RMRE and RMTE of the proposed method are much lower than the HLCC method. So the proposed method has better motion estimation accuracy for real human experiments.

Table 5

Real head motion estimation results of human head images

		Image1	Image2
RMRE	GMM	25.31%	26.68%
	HLCC	35.58%	36.91%
RMTE	GMM	45.46%	46.89%
	HLCC	106.9%	108.8%

All the reconstructed images are shown in Figs. 7 and 8. To view the texture of reconstructed images more clearly, the zoomed images of a region of interest (ROI) from the corresponding images of Figs. 7 and 8 are shown in Figs. 9 and 10. We can see that the motion-corrected images obtained by our proposed method are very close to the stationary simulation phantom. Some streaks are observed near bony structures in Figs. 9(d) and 10(d).

Fig.7

Motion artifact reduction in Image1. (a) is stationary simulation phantom; (b) is reconstructed by the traditional FBP algorithm with no motion compensation; (c) and (d) are reconstructed by the improved FBP algorithm with motion compensation based on our proposed method and the HLCC method.

Fig.8

Motion artifact reduction in Image2. (a) is stationary simulation phantom; (b) is reconstructed by the traditional FBP algorithm with no motion compensation; (c) and (d) are reconstructed by the improved FBP algorithm with motion compensation based on our proposed method and the HLCC method.

Fig.9

Zoomed images of a ROI from Fig. 7.

Fig.10

Zoomed images of a ROI from Fig. 8.

To quantitatively evaluate the proposed method, we compared respectively the reconstructed images in Figs. 7 and 8 by the RMSE, CC and MSSIM as shown in Tables 6 and 7. The RMSE of the proposed method is the smallest and reduce about 39% than the HLCC method; the CC and MSSIM of the proposed method are the greatest and increase about 23% and 40% than the HLCC method. Therefore, the proposed method has a better performance on motion artifact reduction for real human experiments.

Table 6

Comparison of artifact reduction performance on human head for Image1

	RMSE	CC	MSSIM
Uncorrected	0.1400	0.6823	0.3849
GMM	0.0786	0.8977	0.7126
HLCC	0.1295	0.7292	0.5017

Table 7

Comparison of artifact reduction performance on human head for Image2

	RMSE	CC	MSSIM
Uncorrected	0.1474	0.6445	0.3633
GMM	0.0794	0.8476	0.6716
HLCC	0.1323	0.6913	0.4772

4 Conclusions

In this study, a rigid motion estimation and compensation method based on frequency domain analysis has been proposed in fan-beam geometry. The proposed method can determine the motion parameters by the magnitude correlation of projection data in frequency domain and compensate this motion in the reconstruction. In fan-beam CT, the experiments results on both simulated phantom and human head images have demonstrated that the proposed method is better than the HLCC method in estimation accuracy, real-time and artifact correction effect. Therefore, this method could be applicable to helical CT. The future development of our work includes further testing and validation on helical CT. Another development will concern the affine motion especially the dilatation will be estimated and compensated.

Footnotes

Acknowledgments

This work was supported by Natural Science Foundation of China (No.61340034), China Postdoctoral Science Foundation (No.2013M530873), Tianjin Research Program of Application Foundation and Advanced Technology (No.13JCYBJC15600), Natural Science Foundation of Tianjin of China (No.16JCYBJC28800) and International Excellent Postdoctor Training Program of Tianjin of China.

References

Schöndube

, Allmendinger

, Stierstorfer

, et al., Temporal resolution and motion artifacts in single-source and dual-source cardiac CT, Medical Physics 40 (2013), 031112.

Chen

, Shi

, Feng

, et al., Artifact suppressed dictionary learning for low-dose CT image processing, IEEE Transactions on Medical Imaging 33 (2014), 2271–2292.

Tang

, Tong

, Dong

B.X.

and Dillenseger

J.L.

, A new stationary gridline artifact suppression method based on the 2d discrete wavelet transform, Medical Physics 42 (2015), 1721–1729.

Alsobou

, Ahmad

and Ali

, A motion algorithm to extract physical and motion parameters of mobile targets from cone-beam computed tomographic images, Journal of X-ray Science and Technology 24 (2016), 599–613.

Wagner

, Schicho

, Kaniberger

, et al., Quantification and clinical relevance of head motion during computed tomography, Investigative Radiology 38 (2003), 733–741.

Boas

F.E.

and Fleischmann

, CT artifacts: Causes and reduction techniques, Imaging in Medicine 4 (2012), 229–240.

and Mueller

, Efficient low-dose CT artifact mitigation using an artifact-matched prior scan, Medical Physics 39 (2012), 4748–4760.

Zeng

L.L.

, Wang

, Fox

M.D.

, et al., Neurobiological basis of head motion in brain imaging, Proceedings of the National Academy of Sciences of the United States of America 111 (2014), 6058–6062.

Smith-Bindman

, Lipson

and Marcus

, Radiation dose associated with common computed tomography examinations and the associated lifetime attributable risk of cancer, Archives of Internal Medicine 169 (2009), 2078–2086.

10.

Chen

, Chen

, Yin

, et al., Improving low-dose abdominal CT images by Weighted Intensity Averaging over Large-scale Neighborhoods, European Journal of Radiology 80 (2011), 42–49.

11.

Reading

, Radiation exposure from CT scans in childhood and subsequent risk of leukaemia and brain tumours: A retrospective cohort study, Lancet 380 (2012), 910–910.

12.

Chen

, Yin

, Shi

, et al., Improving abdomen tumor low-doseCT images using a fast dictionary learning based processing, Physics in Medicine & Biology 58(2013), 5803–5820.

13.

Kim

J.H.

, Nuyts

, Kuncic

, et al., A method of motion tracking during CT for motion correction, Nuclear Science Symposium and Medical Imaging Conference, NSS/MIC, Knoxville, TN, USA, IEEE, 2011, pp. 3816–3819.

14.

Reaungamornrat

, Wang

A.S.

, Uneri

, et al., Deformable image registration with local rigidity constraints for cone-beam CT-guided spine surgery, Physics in Medicine & Biology 59 (2014), 3761–3787.

15.

Vaillant

, Prieto

, Kolbitsch

and Penney

, Retrospective rigid motion correction in k-space for segmented radial MRI, IEEE Transactions on Medical Imaging 33 (2014), 1–10.

16.

, Chen

, Jin

, et al, Experimental measurement of human head motion for high-resolution computed tomography system design, Optical Engineering 49 (2010), 402–402.

17.

Fahmi

, Marquering

H.A.

, Streekstra

G.J.

, et al., Automatic detection of CT perfusion datasets unsuitable for analysis due to head movement of acute ischemic stroke patients, Journal of Healthcare Engineering 5 (2014), 67–78.

18.

Jenkinson

, Bannister

, Brady

and Smith

, Improved optimization for the robust and accurate linear registration and motion correction of brain images, Neuroimage 17 (2002), 825–841.

19.

Sruthi

M.V.

, Sree

V.U.

and Soundararajan

, Specific removal of motion artifacts in medical image processing, International Journal of Electronics Communication and Computer Engineering 3 (2012), 442–444.

20.

Kim

B.S.

, Lee

K.R.

and Goh

M.J.

, New imaging strategies using a motion-resistant liver sequence in uncooperative patients, Biomed Research International 2014 (2014), 142658–142658.

21.

You

, Serag

, Evangelou

I.E.

, et al., Robust motion correction and outlier rejection of in vivo functional MR images of the fetal brain and placenta during maternal hyperoxia, Proc Spie 9417 (2015), 363–371.

22.

Chandler

, Wei

and Herron

D.H.

, Semiautomated motion correction of tumors in lung CT-perfusion studies, Academic Radiology 18 (2011), 286–293.

23.

Kim

J.H.

, Nuyts

, Kuncic

and Fulton

, The feasibility of head motion tracking in helical CT: A step toward motion correction, Medical Physics 40 (2013), 041903–041903.

24.

Kim

J.H.

, Nuyts

, Kyme

, et al., A rigid motion correction method for helical computed tomography (CT), Physics in Medicine & Biology 60 (2015), 2047–2073.

25.

Wang

and Vannier

M.W.

, Preliminary study on helical CT algorithms for patient motion estimation and compensation, IEEE Transactions on Medical Imaging 14 (1995), 205–211.

26.

, Wei

, Hsieh

, et al., Data consistency based translational motion artifact reduction in fan-beam CT, IEEE Transactions on Medical Imaging 25 (2006), 792–803.

27.

Berger

, Müller

, Choi

J.H.

, et al., 2D/3D registration for motion compensated reconstruction in cone-beam CT of knees under weight-bearing condition, World Congress on Medical Physics and Biomedical Engineering, Springer International Publishing, Toronto, Canada, 2015, pp. 54–57.

28.

and Wang

, Data Consistency Based Rigid Motion Artifact Reduction in Fan-Beam CT, IEEE Transactions on Medical Imaging 26 (2007), 249–260.

29.

Dang

, Stayman

J.W.

, Sisniega

, et al., Statisticalreconstruction for cone-beam CT with a post-artifact-correctionnoise model: Application to high-quality head imaging, Physics in Medicine & Biology 60(2015), 6153–6175.

30.

Sun

, Kim

J.H.

, Fulton

, et al., An iterative projection-based motion estimation and compensation scheme for head x-ray CT, Medical Physics, 43 (2016), 5705–5716.

31.

Zhao

, Yang

and Yang

, Fan beam image reconstruction with generalized fourier slice theorem, Journal of X-ray Science and Technology 22 (2014), 415–436.

32.

Riviere

P.J.L.

and Pan

, Fourier-based approach to interpolation in helical CT exploiting redundant fan-beam information. Medical Imaging, International Society for Optics and Photonics 3977 (2000).

33.

Hoge

W.S.

, A subspace identification extension to the phase correlation method, IEEE Transactions on Medical Imaging 22 (2003), 277–280.

34.

Druckmüller

, Phase correlation method for the alignment of total solar eclipse images, Astrophysical Journal 706 (2009), 1605–1608.

35.

Tong

, Ye

, Xu

, et al., A Novel Subpixel Phase Correlation Method Using Singular Value Decomposition and Unified Random Sample Consensus, Sensing 53 (2015), 4143–4156.

36.

Yousef

, Li

and Karim

, High-speed image registration algorithm with subpixel accuracy, IEEE Signal Processing Letters 22 (2015), 1796–1800.

37.

Luo

, Huang

, Pan

, et al., A frequency domain based rigid motion artifact reduction algorithm, International Symposium on Multispectral Image Processing and Pattern Recognition (2009), 74972K–74972K-8.

38.

Sheikh

H.R.

, Bovik

A.C.

and De

V.G.

, An information fidelity criterion for image quality assessment using natural scene statistics, IEEE Transactions on Image Processing 14 (2006), 2117–2128.

Rigid motion artifact reduction in CT using frequency domain analysis

Abstract

BACKGROUND:

OBJECTIVE:

METHODS:

RESULTS:

CONCLUSIONS:

Keywords

1 Introduction

2 Materials and methods

2.1 Motion estimation in fan-beam

2.1.1 Fan-Beam projection with general in-plane motion

Table 1 Motion patterns I II III φ (t) 0 2 (1 - e t ) 5.73 (cos (πt/ - 2) - 1) d x (t) 0 3 (1 - cos (πt/ - 2)) 4 (1 - e t ) d y (t) 0 3 sin (πt) 5 (cos (πt/ - 2) - 1) MRE/° 0 1.437 2.0833 MTE/mm 0 2.4163 3.4164

Footnotes

Acknowledgments

References

Table 1
Motion patterns

I II III

φ (t) 0 2 (1 - e^t) 5.73 (cos (πt/ - 2) - 1)

d_x (t) 0 3 (1 - cos (πt/ - 2)) 4 (1 - e^t)

d_y (t) 0 3 sin (πt) 5 (cos (πt/ - 2) - 1)

MRE/° 0 1.437 2.0833

MTE/mm 0 2.4163 3.4164