2.1

Motion Compensated Weighted Filtered Backprojection

In this paper we combine the motion compensation method proposed by Schäfer et al. [4] with the WFBP algorithm of Stierstorfer et al. [6], a 3D filtered backprojection algorithm for multislice spiral CT. However, our proposed motion compensation method works with any CT reconstruction algorithm using rebinning and backprojection steps.

In the rebinning step of the WFBP algorithm the projection images acquired in the cone beam geometry are virtually rearranged to form semi-parallel projection images better suited for the backprojection step. The cone beam geometry is described using α as the rotation angle, β as the horizontal opening angle of a ray and q as the cone parameter describing the vertical cone angle of a ray. For the semi-parallel projection images the rotation angle is described as θ and (p, q) denote the detector columns and rows, respectively. The geometry is only semi-parallel, as the rays are only parallel along the detector rows, but not along the columns. In the rebinning step the rays of the initially obtained cone beam projection images P(α, β, q) are rearranged to form semi-parallel projection images P(θ, p, q) using the rebinning formulas

with R describing the distance from source to isocenter.

The backprojection step in the WFBP algorithm iterates through all rebinned and filtered projection images P(θ, p, q) and for each projection image it iterates through all voxels x of the voxel volume to be reconstructed. For each voxel x, the corresponding position (p, q) in the projection image P(θ, p, q) the voxel has been projected to during acquisition is calculated. The value at position (p, q) is then added to voxel x in the voxel volume. This is performed for all voxels and all projection images until the entire voxel volume has been reconstructed.

For our motion compensation method we assume that the motion present during acquisition is known. For every projection image P(α, β, q) a motion state M_α is required. Each motion state M_α holds the information about the position of the scanned object relative to the scanner at the moment of the acquisition of P(α, β, q).

Our proposed algorithm is based on the method by Schäfer et al. [4], extending it to work with reconstruction algorithms containing a rebinning step like the WFBP algorithm. When combining the WFBP algorithm with Schäfer’s method the backprojection step of the reconstruction algorithm is extended by moving every voxel x in the voxel volume to be reconstructed according to the known motion pattern. As described above, the backprojection is performed by iterating for each projection image P(θ, p, q) through all voxels x. Knowing the motion state M_α present during the acquisition of the projection image P(α, β, q) obtained at rotation angle α, we can use this information to move each voxel x accordingly, assuming α = θ. For each projection image P(θ, p, q) each voxel position x is recalculated to x_M,α considering the motion state M_α=_θ. Using the new position x_M,α the corresponding position (p, q) in the projection image P(θ, p, q) is calculated and the value of this position is added to the value at voxel position x. This is then repeated for all voxels in the voxel volume and for all obtained projection images. This leads to a more accurate reconstruction, as the motion present during acquisition is considered. However, when combining the WFBP algorithm with the motion compensation method proposed by Schäfer et al. in this manner, some motion artifacts remain in the final result, as the rebinning process is not considered.

Due to the rebinning the assumption that for every projection image P(θ, p, q) one motion state M_α needs to be considered during backprojection is not accurate. This assumption only holds for the projection images P(α, β, q). The rebinning process however rearranges the rays such that the rebinned projection images each consist of information obtained at different time points and hence depicting different motion states. The available motion state M_α is no longer correct for the entire projection image P(θ, p, q) but only for the column where p = 0, as only there θ = α holds. Using our geometry, the columns right of the column p = 0 contain previous motion states and columns towards the left contain subsequent motion states. Hence using a motion state M_α for an entire projection image P(θ, p, q) in the backprojection will lead to motion artifacts remaining in the final result.

We propose a motion compensation method as described in algorithm 1 extending Schäfer’s method to consider the rebinning process, enabling a combination with the WFBP reconstruction algorithm. Steps 1) and 2) are the steps as described in the method of Schäfer et al., while steps 3)-6) are our proposed extension. Our method uses the motion state M_α=θ to calculate a new position x_M,θ of the voxel x, which is then used to calculate an initial position (p, q) in the projection image P(θ, p, q) just as described by Schäfer et al. This position (p, q) is then used to choose a more accurate motion state M’. Using the column position p and the rebinning formulas described in equation (1) the cone beam projection image P(α, β, q) from which the information in this column p was taken in the rebinning and hence which motion state M’ was present can be calculated. With this new more accurate motion state M’ the voxel position x can be recalculated to x’_M,θ, followed by the calculation of the new position (p’, q’) in the projection image P(θ, p, q). This more accurate position (p’, q’) can then be used to write its values to the voxel x. The steps 3) to 5) can be repeated to achieve even more accuracy, however experiments have shown that even one iteration is able to yield good results.

For a more efficient implementation, we propose a modification of the algorithm stated above. Instead of using the exact motion state M’ that was calculated for every position p in the projection image P(θ, p, q), we suggest choosing 3 motion states per projection image, whilst interpolating the motion states in between. Specifically, we chose to use the motion state at position p = 0 which corresponds to the initial motion state M_α=θ as well as one motion state for one column left and one column right of p = 0. However, we are not choosing the left and rightmost columns in the projection image P(θ, p, q), but the columns furthest to the left and right that are actually reached in a projection. These columns need to be calculated individually, as they depend on several parameters. For the three chosen columns the corresponding motion states are calculated using the rebinning equations (1). For a position (p, q) in the projection image P(θ, p, q) the corresponding motion state is then obtained from the motion states of the two closest chosen columns using linear interpolation of the motion values.

2.2

Motion Compensated Iterative Weighted Filtered Backprojection

The method described above can also be combined with an iterative weighted filtered backprojection algorithm, as the one described by Sunnegårdh et al. [7]. The iterative reconstruction algorithm can be described with the vectors p and f denoting the projection data and the resulting image data, and the matrices Q and P denoting backprojection and forward projection operators, respectively. In our case, Q is the weighted filtered backprojection by Stierstorfer et al. [6] and P is the forward projection by Joseph [8]. The iterative reconstruction process can then be described as

An initial result f₀ = Qp is calculated by backprojecting the rebinned and filtered projection data. The initial result f₀ is then forward projected and the difference to the initial projection data p is calculated. This difference image is backprojected and subtracted from the previous result f₀ using a factor α. This is repeated multiple times, improving the result with every iteration.

Including the proposed motion compensation method into this iterative reconstruction process consists of two parts. First, the motion compensation is included into the backprojection operator Q as described above. Second, the motion also needs to be considered in the forward projection operator P. However, in the forward projection the motion is included again rather than removed, as the forward projected result is compared to the original not motion compensated projection data p. Including the motion in the forward projection operator P is performed similar to the motion compensation in the backprojection operator Q. In the forward projection, we are assuming an implementation that iterates through all positions (p, q) in the projection image P’a(θ, p, q) that should be calculated. For each position (p, q), the corresponding ray through the voxel volume is calculated and the values along its path are accumulated. This is repeated for all projection images P’(θ, p, q). Knowing the position (p, q) the corresponding motion state at this position can be interpolated as described above. The motion state can then be included in the calculations of each ray changing its path through the voxel volume accordingly. The value accumulated along the rays path is written to the initial position (p, q) of the projection image P’(θ, p, q). However, the motion compensation calculations need to be included in every backprojection and every forward projection step of every iteration.

Hence, we propose to use a modification of the iterative reconstruction algorithm used above, as the one described by Sunnegårdh et al. [9], that allows for a more efficient implementation of the proposed motion compensation algorithm. The algorithm stated in equation (2) is adapted by considering the combination of forward projection followed by backprojection to be the sum of a low-pass filter L and a linear operator A that causes artifacts:

Hence the iterative reconstruction equation (2) can be rewritten to

as described in more detail by Sunnegårdh et al. [9]. This modification allows to include the motion compensation algorithm only in the first backprojection step, without the need to repeat the motion compensation calculations in every iteration. Furthermore, there is no need to adapt the forward projection operator to consider motion.