A three-dimensional dual-domain deep network for high-pitch and sparse helical CT reconstruction

Computed tomography (CT) has been an important tool in medical diagnosis and the helical scan with multi-row detector is the most used scanning modality in hospitals. When scanning, the X-ray source and detector rotate around the central axis (i.e., the z-direction) while the table with a patient remaining stationary on it is translated along the z-direction at a constant velocity. Therefore, the trajectory of the X-ray source relative to the patient forms a helical curve, where the pitch of the helical curve is related to the moving velocity of the table.

The patient needs to stay stationary during the scanning process, otherwise severe artifacts will be caused in the reconstructed CT images. Thus, the high pitch scanning that can complete the data acquiring in small time has many merits and is very desirable in clinical diagnosis. For example, by utilizing the high pitch helical CT, the whole lung can be scanned within a very small interval and so the patient needs not to hold the breath. This is very useful for patients that can’t perform breath-holding for a long time, such as infant, old patients and unconscious patients. For cardiac CT, scanning the whole heart in one heart beat is a reasonable request.

In [1], Noo et al. used the Tam-Danielsson window [2][3] and Katsevich’s inverse formula [4] to derive the numbers of required detector rows for perfectly reconstructing the helical CT images with a curved detector or a flat detector. When other parameters are fixed, the higher the pitch is, the greater number of the detector rows is needed. However, as the number of the detector rows increasing, the price of the CT scanner will also rise significantly. One way to reduce the price is to use sparse detector [5] with fixed detector row numbers. However, CT images reconstructed from sparse sinograms usually have streak artifacts and lose fine details especially when the scanning pitch is high.

In the literature, many algorithms were proposed to reconstruct the helical CT images. These algorithms can be mainly categorized into four classes: (a) rebinning algorithms, (b) FDK-like algorithms, (c) model-based iterative methods, (d) exact algorithms.

The rebinning algorithms [6][7] first converted the measured helical sinograms to 2D fan-beam or parallel-beam sinograms and then reconstructed the CT images by using any 2D reconstruction methods. To reduce the interpolation error and accurately estimate the 2D sinogram, the sampling rate of the helical scan is usually high.

The FDK-like algorithms [8][9][10] are extended from the circle cone-beam reconstruction algorithm proposed by Feldkamp, Davis and Kress [11]. The FDK-like algorithms usually involve three steps: 1) filtering the sinograms, 2) weighting the filtered data, 3) Backprojecting the weighted data. By utilizing the fast Fourier transform (FFT), the FDK-like algorithms can be implemented efficiently. However, since the FDK-like algorithms were proposed heuristically and lack of theory supports, the artifacts in the reconstruction is usually hard to avoided especially when the cone angles are large.

Model-based methods usually first propose the minimizing functionals by combining the statistical property of the sinogram with different priori hypotheses, and then convert the functionals to iterative reconstruction algorithms by using optimizing methods. In the literature, total variation (TV) [12][13], nonlocal means (NLM) [14], dictionary learning [15][16] are the commonly used priori hypotheses. In [7], Stierstorfer et al. proposed an iterative weighted filtered backprojection (WFBP) method to reconstruct helical CT images. In [17], Sunneg$\mathring{a}$rdh proposed an improved regularized iterative weighted filtered backprojection (RIWFBP) method for the helical cone-beam CT. In [18], Yu and Zeng developed a fast TV-based iterative reconstruction algorithm for the limited-angle helical cone-beam CT reconstruction. Model-based iterative methods are very suitable for ill-posed problems, such as limited angle or sparse angle CT reconstructions. But the computational burdens of the model-based iterative methods are too high to reconstruct the images in real time especially for 3D helical CT.

Exact helical CT reconstruction can be obtained by using Katsevich’s algorithms [19][4], Zou and Pan’s PI-line algorithms [20][21]. These algorithms have mathematical inverse formulas for helical CT reconstructions and so are theoretically exact. In [22][23], by extending Katsevich’s algorithms and the PI-line algorithms, the exact methods for cone-beam CT reconstructions with general scanning curves were proposed. In [24], Yan et al. accelerated the implementation of the Katsevich algorithm [4] by utilizing the graphics processing unit (GPU). In their implementation, the Tam–Danielsson window rather than the PI line is used to determine whether a sinogram corresponding to a scanning angle needs to be backprojected to reconstruct a slice of CT image. Their method doesn’t utilize the periodic properties of the PI-line and so needs to calculate the backprojection positions on the detector for pixels in all the CT images. Exact helical algorithms can reconstruct good quality CT images even if the helical pitch is high. But when the sinograms are sparse, streak artifacts in the CT images reconstructed by the exact algorithms can’t be avoided.

Recently, convolutional neural network (CNN) based methods have also been used for CT reconstructions. These methods can be roughly divided into three categories. The first category uses CNNs only to pre-process the sinograms [25][26] or post-process the CT images [27][28][29]. This type of methods needs an extra step to convert the sinograms to CT images. The second category utilizes end-to-end CNNs to reconstruct CT images directly from measured sinograms [30][31][32][5][33]. Since the reconstruction algorithm is modelled as network layers embedded in CNNS, no extra step is needed when testing or deploying these trained CNNs. However, embedding reconstruction algorithm in CNNs is GPU-memory expensive and may result in GPU-memory exhausting error when training this type of CNNs especially for 3D CT. The third category uses CNNs to approximate the solution of the unroll iterative algorithms [34][35], where each iteration of the algorithms was approximated by a subnetwork. This type of algorithms exhausts more GPU-memory than those in the second category since in the networks the operators of projection and backprojection need to be performed $n$ times, where $n$ is the iteration numbers.

In this paper, we focus on embedding the Katsevich algorithm [4] into CNNs to propose an end-to-end deep network to reconstruct helical CT images with high pitch and sparsely spaced detectors. Due to the limited GPU-memory, to the best of our knowledge, in the literature no exact algorithm (Katsevich’s algorithms or PI-line algorithms) was embedded in CNNs to reconstruct helical CT images. In [5], an end-to-end deep network was proposed for helical CT, but the sinograms measured from helical scans were converted to 2D fan-beam sinograms and the embedded reconstruction algorithm was also for fan-beam geometry. As the helical pitch increases, the image volume involved to generate the rebinned fan-beam sinogram and the error of the estimated fan-beam sinogram will both increase. Therefore, for helical CT with high pitch, the CNN using rebining methodology may not work well. The Katsevich algorithm [4] can reconstruct good CT images even if the helical pitch is high. Our network uses Katsevich algorithm as the domain transfer layer and so can also reconstruct good CT images directly from high pitch helical sinograms. Moreover, because our network utilizes the features extracted from both sinograms and CT images, it can simultaneously reduce the streak artifacts caused by the sparsity of sinograms and reconstruct fine details in the CT images.

The sinogram measured from the helical scan usually has large size in the z-direction. If directly input the whole sinogram into CNN to reconstruct all the CT images, huge GPU-memory is needed. Besides, most of the parameters of Katsevich algorithm [4] depend on the position of the CT images and computing them is also time consuming and GPU-memory consuming. By observing that the parameters of Katsevich algorithm [4] are periodic, we divide the sinograms into parts and reconstruct the CT images cycle-by-cycle. By this way, the parameters of Katsevich algorithm for all cycles can be pre-calculated and the GPU-memory loading can be reduced.

The rest of the paper is organized as follows. Section II gives our new GPU implementation of the Katsevich algorithm [4]. Section III describes the detailed structure of our proposed network. Section IV presents the experimental results. Discussion and conclusion are given in Section V.

In this section, we give a new GPU implementation of the Katsevich algorithm [4], which reconstructs the helical CT images pitch by pitch and has low computational and GPU-memory burdens. Our implementation of Katsevich algorithm [4] is very suitable for deep learning since its paramerters are the same for all pitches and so need only to calculate once.

Katsevich [4] proposed an exact reconstruction formula for helical cone-beam CT and Noo et al. [1] researched how to implement it efficiently and accurately. Let

be the helical trajectory formed by the X-ray source position with radius $R$ and pitch $p$, where

is the start point of the helix,

is a cylinder inside the helix, $0<r<R$ and $\underline{x}=(x,y,z)$.

Then, for any $\underline{x}\in U$, the reconstruction formula can be written as

where $\lambda_{i}(\underline{x})$ and $\lambda_{o}(\underline{x})$ are the extremities of the $\pi$-line passing through $\underline{x}$ with $\lambda_{i}(\underline{x})<\lambda_{o}(\underline{x})$.

is the Hilbert filter, $g(\lambda,\underline{\theta})$ is the measured projection data at X-ray source position $\underline{a}(\lambda)$ in the direction $\underline{\theta}$, vector $\underline{m}(\lambda,\underline{\theta})$ is normal to the $\kappa$-plane $\mathcal{K}(\lambda,\psi)$ of the smallest $|\psi|$ value that contains the line of direction $\underline{\theta}$ through $\underline{a}(\lambda)$ and the $\kappa$-plane $\mathcal{K}(\lambda,\psi)$ is any 2D plane that has three intersections $\underline{a}(\lambda)$, $\underline{a}(\lambda+\psi)$ and $\underline{a}(\lambda+2\psi)$ with the helix, where $\psi\in(-\pi,\pi)$. In [1], to numerically calculate $g^{F}(\lambda,\underline{\theta})$, the integral over $\gamma$ on the $\kappa$-plane with the minimum value $\left|\psi\right|$ is converted to the line integral along the curve that the $\kappa$-plane intersects with the detector plane.

Let $g(\lambda,\alpha,w)$ be the sinogram measured from the helical scan by curved detector, then the implementation of (4) in [1] is described as follows:

1. 1.

   Taking the derivative at constant direction via chain rule:

   $$g_{1}(\lambda,\alpha,w)=\left.\left(\frac{\partial g(q,\alpha,w)}{\partial q}+\frac{\partial g(q,\alpha,w)}{\partial\alpha}\right)\right|_{q=\lambda}$$

   (8)

2. 2.

   Length-correction weighting:

   $$g_{2}(\lambda,\alpha,w)=\frac{D}{\sqrt{D^{2}+w^{2}}}g_{1}(\lambda,\alpha,w)$$

   (9)

3. 3.

   Forward height rebinning: computing

   $$g_{3}(\lambda,\alpha,\psi)=g_{2}\left(\lambda,\alpha,w_{\kappa}(\alpha,\psi)\right)$$

   (10)

   by interpolation for all $\psi\in\left[-\pi/2-\alpha_{m},\pi/2+\alpha_{m}\right]$, where $\alpha_{m}=\arcsin(r/R)$ is the half fan angle,

   $$w_{\kappa}(\alpha,\psi)=\frac{DP}{2\pi R_{o}}\left(\psi\cos\alpha+\frac{\psi}{\tan\psi}\sin\alpha\right)$$

   (11)

4. 4.

   1-D Hilbert transform in $\alpha$: at constant $\psi$, compute

   $$g_{4}(\lambda,\alpha,\psi)=\int_{-\pi/2}^{+\pi/2}\mathrm{~{}d}\alpha^{\prime}h_{H}\left(\sin\left(\alpha-\alpha^{\prime}\right)\right)g_{3}\left(\lambda,\alpha^{\prime},\psi\right)$$

   (12)

   where $h_{H}$ is the Hilbert filter defined in (7).

5. 5.

   Backward height rebinning: compute

   $$g_{5}(\lambda,\alpha,w)=g_{4}(\lambda,\alpha,\hat{\psi}(\alpha,w))$$

   (13)

   where $\hat{\psi}(\alpha,w)$ is angle $\psi$ of the smallest absolute value that satisfies the equation

   $$w=\frac{DP}{2\pi R_{o}}\left(\psi\cos\alpha+\frac{\psi}{\tan\psi}\sin\alpha\right)$$

   (14)

6. 6.

   Post-cosine weighting:

   $$g^{F}(\lambda,\alpha,w)=\cos\alpha g_{5}(\lambda,\alpha,w)$$

   (15)

7. 7.

   Backprojection:

   $$f(\underline{x})=\frac{1}{2\pi}\int_{\lambda_{i}(\underline{x})}^{\lambda_{o}(\underline{x})}\mathrm{d}\lambda\frac{1}{v^{*}(\lambda,\underline{x})}g^{F}\left(\lambda,\alpha^{*}(\lambda,\underline{x}),w^{*}(\lambda,\underline{x})\right)$$

   (16)

   where

   $$v^{*}(\lambda,\underline{x})=R_{o}-x\cos\left(\lambda+\lambda_{0}\right)-y\sin\left(\lambda+\lambda_{0}\right)$$

   (17)

   	$\displaystyle\alpha^{*}(\lambda,\underline{x})=$			(18)
   	$\displaystyle\arctan$	$\displaystyle\left(\frac{\left(-x\sin\left(\lambda+\lambda_{0}\right)+y\cos\left(\lambda+\lambda_{0}\right)\right)}{v^{*}(\lambda,\underline{x})}\right)$

   $$w^{*}(\lambda,\underline{x})=\frac{D\cos\alpha^{*}(\lambda,\underline{x})}{v^{*}(\lambda,\underline{x})}\left(z-z_{0}-\frac{P}{2\pi}\lambda\right)$$

   (19)

In the above steps, the parameters $\lambda_{i}(\underline{x})$, $\lambda_{o}(\underline{x})$, $v^{*}(\lambda,\underline{x})$, $\alpha^{*}(\lambda,\underline{x})$, $w^{*}(\lambda,\underline{x})$, $w_{\kappa}(\alpha,\psi)$ and $\hat{\psi}(\alpha,w)$ need to calculate. The parameters $w_{\kappa}(\alpha,\psi)$ and $\hat{\psi}(\alpha,w)$ can be pre-calculated only once but the others $\lambda_{i}(\underline{x})$, $\lambda_{o}(\underline{x})$, $v^{*}(\lambda,\underline{x})$, $\alpha^{*}(\lambda,\underline{x})$, $w^{*}(\lambda,\underline{x})$ need to be calculated for every $\underline{x}$, which is time consuming. Luckily, these parameters have the following periodic properties:

where $\underline{x}^{k}=(x,y,z+kp)$, $k$ is an integer and $p$ is the pitch. Thus, if we reconstruct the CT images in the z-direction pitch by pitch, the parameters $\lambda_{o}(\underline{x})$, $\lambda_{p}(\underline{x})$, $v^{*}(\lambda,\underline{x})$, $\alpha^{*}(\lambda,\underline{x})$, $w^{*}(\lambda,\underline{x})$ can also be pre-calculated only once, which can reduce the computation and GPU-memory burdens a lot.

In the following, we describe the detailed GPU implementation scheme of the reconstruction formula (4), which reconstruct the CT images pitch by pitch. In our scheme, to reduce the GPU-memory consuming, the backprojection step is performed slice by slice in the z-direction while the other steps are performed on the whole sinogram needed to reconstruct the CT images of one pitch.

Let $f(\underline{x}_{j}^{k})$ be the $j$-th sliced CT image of the $k$-th pitch in the z-direction, where $j=0,1,2,...,N-1$, $k=1,...M$, $N$ is the total number of CT images to be reconstruct in one pitch, $M$ is the number of divided pitches of the scanning helix, $\underline{x}_{j}^{k}=(x,y,z_{j}+kp)$, $(x,y)$ is a point of the mesh grids formed by the x-coordinate and y-coordinate and $z_{j}$ is the relative position on the z-coordinate in one pitch. Then our scheme can be divided into parameters pre-calculating step and reconstruction step.

The pre-calculating step is described as follows:

1. 1.

   Pre-calculate $w_{\kappa}(\alpha,\psi)$ and $\hat{\psi}(\alpha,w)$ by using equation (11) and the method in [1].

2. 2.

   Define $\underline{x}_{j}:=\underline{x}_{j}^{1}$, For $j=0,1,2,...,N-1$, pre-calculate the extremities $\lambda_{i}(\underline{x}_{j})$ and $\lambda_{o}(\underline{x}_{j})$ of the $\pi$-line passing through $\underline{x}_{j}$ using the algorithm described in [36]. Note that an extra projection step that constrains $\theta$ on the interval $[-\pi,\pi]$ needs to perform at each iteration of the algorithm, which was not explicitly stated in [36]. Also note that here $\lambda_{i}(\underline{x}_{j})$ and $\lambda_{o}(\underline{x}_{j})$ are two matrices corresponding to the mesh grid points $(x,y)$ of $\underline{x}_{j}$. After obtaining $\lambda_{i}(\underline{x}_{j})$ and $\lambda_{o}(\underline{x}_{j})$, we can calculate the minimal and maximal values of $\lambda$ needed to reconstruct $f(\underline{x}_{j})$:

   	$\displaystyle\lambda_{min}(\underline{x}_{j})$	$\displaystyle=\min_{(x,y)~{}\text{of}~{}x_{j}}(\lambda_{i}(\underline{x}_{j})),$		(21)
   	$\displaystyle\lambda_{max}(\underline{x}_{j})$	$\displaystyle=\max_{(x,y)~{}\text{of}~{}x_{j}}(\lambda_{o}(\underline{x}_{j})).$

3. 3.

   For $j=0,1,2,...,N-1$, pre-calculate $v^{*}(\lambda,\underline{x}_{j})$, $\alpha^{*}(\lambda,\underline{x}_{j})$ and $w^{*}(\lambda,\underline{x}_{j})$:

   $$v^{*}(\lambda,\underline{x}_{j})=R_{o}-x\cos\left(\lambda+\lambda_{0}\right)-y\sin\left(\lambda+\lambda_{0}\right)$$

   (22)

   	$\displaystyle\alpha^{*}(\lambda,\underline{x}_{j})$	$\displaystyle=$		(23)
   	$\displaystyle\arctan$	$\displaystyle\left(\frac{\left(-x\sin\left(\lambda+\lambda_{0}\right)+y\cos\left(\lambda+\lambda_{0}\right)\right)}{v^{*}(\lambda,\underline{x}_{j})}\right)$

   $$w^{*}(\lambda,\underline{x}_{j})=\frac{D\cos\alpha^{*}(\lambda,\underline{x}_{j})}{v^{*}(\lambda,\underline{x}_{j})}\left(z_{i}-z_{0}-\frac{P}{2\pi}\lambda\right),$$

   (24)

   where $\underline{x}_{j}=(x,y,z_{j})$ and $\lambda=[\lambda_{min}(\underline{x}_{j}),\lambda_{max}(\underline{x}_{j})]$. Here $u^{*}(\lambda,\underline{x}_{j})$, $\alpha^{*}(\lambda,\underline{x}_{j})$ and $w^{*}(\lambda,\underline{x}_{j})$ are three tensors of dimension three, where the first dimension corresponds to $\lambda=[\lambda_{min}(\underline{x}_{j}),\lambda_{max}(\underline{x}_{j})]$.

Note that in the above step 3), since

for any positive integer $k$, $v^{*}(\lambda,\underline{x}_{j})$, $\alpha^{*}(\lambda,\underline{x}_{j})$ and $w^{*}(\lambda,\underline{x}_{j})$ can be pre-calculated only once.

After the pre-calculating step, we can compute the minimal and maximal values of $\lambda$, $\alpha$ and $w$ needed to reconstruct the CT images of one pitch. These values can be calculated by

The values $\alpha_{min}$, $\alpha_{max}$, $w_{min}$ and $w_{max}$ can be used to determine the numbers of rows and columns of the detector, and $\lambda_{min}$ and $\lambda_{max}$ are used to determine the required sinogram to reconstruct the CT images of one pitch. To reconstruct the $k$-th pitch CT images, the sinogram corresponding to $\lambda=[\lambda_{min}+2k\pi,\lambda_{max}+2k\pi]$ is needed.

Let $g(\lambda^{k},\alpha,w)$ be the $\lambda$-sliced sinogram from $g(\lambda,\alpha,w)$, where $\lambda^{k}=[\lambda_{min}+2k\pi,\lambda_{max}+2k\pi]$ for $k=1,...M$. Then the reconstruction step can be described as follows:

1. 1.

   Perform Derivative (8), Length-correction weighting (9), Forward height rebinning (10), Hilbert transform (12), Backward height rebinning (13) and Post-cosine weighting (15) on $g(\lambda^{k},\alpha,w)$ and obtain $g^{F}(\lambda^{k},\alpha,w)$.

2. 2.

   For $j=0,1,2,...,N-1$, backproject $g^{F}(\lambda^{k},\alpha,w)$ slice by slice:

   $$f(\underline{x}_{j}^{kp})=\frac{1}{2\pi}\int_{\lambda_{min}(\underline{x}_{j})+2k\pi}^{\lambda_{max}(\underline{x}_{j})+2k\pi}\mathrm{d}\lambda\frac{1}{v^{*}_{j}}g^{F}\left(\lambda^{k},\alpha^{*}_{j},w^{*}_{j}\right)$$

   (27)

   where

   	$\displaystyle v^{*}_{j}$	$\displaystyle=v^{*}(\lambda,\underline{x}_{j})$		(28)
   	$\displaystyle\alpha^{*}_{j}$	$\displaystyle=\alpha^{*}(\lambda,\underline{x}_{j})$
   	$\displaystyle w^{*}_{j}$	$\displaystyle=w^{*}(\lambda,\underline{x}_{j})$

To reconstruct the helical CT images pitch by pitch, we need to perform the pre-calculating step once and store its outputs in memory, and the reconstruction step $M$ times for $M$ pitches. In this paper, the pre-calculating step was implemented by the central processing unit (CPU) while the reconstruction step by GPU.

In this section, we propose a 3D deep learning network to reconstruct CT images directly from sparse helical sinograms with high pitch, where our GPU implementation of the Katsevich algorithm [4] is used as the domain transfer layer in the network.

In this section, we present some simulated experimental results to demonstrate the effectiveness of our method and compare it with the related methods, the nonlocal-TV iterative algorithm [42] and the post-processing method FBPConvNet [43], where the convolutional filters in FBPConvNet are extended from 2D to 3D.

As described in Subsection III-C, the CT images in “the 2016 NIH-AAPM-Mayo Clinic Low Dose CT Grand Challenge” [40] are used to prepare the experimental data, where the data generated from six patients named ’L192’, ’L286’, ’L291’, ’L310’, ’L333’ and ’L506’ are used as the training data and those from the other four patients named ’L067’, ’L096’, ’L109’ and ’L143’ are used as the test data. The training data involve totally 345 pitches (of size $512\times 512\times 10$) for $p=7\pi$ and 260 pitches (of size $512\times 512\times 13$) for $p=14\pi$ and the test data involve 215 pitches for $p=7\pi$ and 163 pitches for $p=14\pi$, where $10\%$ of the training data is randomly chosen as the validation. Note that since the last CT image of one pitch is override with the first CT image of its next pitch, we exclude the last CT image of every pitch in our training and testing data.

For nonlocal-TV iterative algorithm, we tune its parameters empirically in order to obtain the best results. We use the default parameters of FBPConvNet to train it, where the number of training epoches is 100. The CT images reconstructed from the noisy sinograms by the Katsevich algorithm are used as the inputs of FBPConvNet, where the batch size is $batch=5$. We train our network $60$ epoches with batch size $batch=1$. The size of the learnable parameters of FBPConvNet is about 12.4MB while ours is about 2.3MB.

In this paper, we developed a new GPU implementation of the Katsevich algorithm [4] for helical CT, which reconstructs the CT images pitch by pitch. By utilizing the periodic properties of the parameters of the Katsevich algorithm, our scheme only needs to calculate these parameters once for all pitches and so is very suitable for deep learning. By embedding our implementation of the Katsevich algorithm into the CNN, we proposed an end-to-end deep network for the high pitch helical CT reconstruction with sparse detectors. Experiments show that our end-to-end deep network outperformed the Nonlocal-TV based iterative algorithm and the post-processed deep network FBPConvNet both in subjective and objective evaluations.

One limitation of our method is that it requires the pitch $p$ of the helical scan is constant, which inherits from the Katsevich algorithm. In the literature, there exist some exact algorithms for the helical CT reconstruction with variable pitches. However, the periodic properties of the parameters of these algorithms may not hold. Therefore, we need to compute the parameters for every slice when using these algorithms to reconstruct CT images, which will consume a lot of GPU-memory and isn’t fit for deep learning. In our future work, we will research how to solve this issue.