A Deep Learning Reconstruction Framework for Differential Phase-Contrast Computed Tomography with Incomplete Data

TL;DR

This paper introduces a deep learning-based reconstruction framework for differential phase-contrast computed tomography with incomplete data, effectively estimating complete sinograms and improving image quality in challenging scenarios.

Contribution

It proposes a novel deep learning approach tightly integrated with DPC-CT reconstruction, handling incomplete data without complex parameter tuning.

Findings

Validated with synthetic and experimental data

Outperforms traditional algorithms in handling incomplete data

Encapsulates physical imaging model for better generalization

Abstract

Differential phase-contrast computed tomography (DPC-CT) is a powerful analysis tool for soft-tissue and low-atomic-number samples. Limited by the implementation conditions, DPC-CT with incomplete projections happens quite often. Conventional reconstruction algorithms are not easy to deal with incomplete data. They are usually involved with complicated parameter selection operations, also sensitive to noise and time-consuming. In this paper, we reported a new deep learning reconstruction framework for incomplete data DPC-CT. It is the tight coupling of the deep learning neural network and DPC-CT reconstruction algorithm in the phase-contrast projection sinogram domain. The estimated result is the complete phase-contrast projection sinogram not the artifacts caused by the incomplete data. After training, this framework is determined and can reconstruct the final DPC-CT images for a given incomplete phase-contrast projection sinogram. Taking the sparse-view DPC-CT as an example, this framework has been validated and demonstrated with synthetic and experimental data sets. Embedded with DPC-CT reconstruction, this framework naturally encapsulates the physical imaging model of DPC-CT systems and is easy to be extended to deal with other challengs. This work is helpful to push the application of the state-of-the-art deep learning theory in the field of DPC-CT.

Tables

Table 1. TABLE I: Parameters in the DenseBlock

Layers	Input Size	Output Size
BN-Relu-$5\times 5$ Conv.	$H\times W\times C_I$	$H\times W\times k$
Concatenate	$H\times W\times C_I$	$H\times W\times \left(C_I+k\right)$
	$H\times W\times k$
BN-Relu-$5\times 5$ Conv.	$H\times W\times \left(C_I+k\right)$	$H\times W\times k$
Concatenate	$H\times W\times \left(C_I+k\right)$	$H\times W\times \left(C_I+2k\right)$
	$H\times W\times k$
BN-Relu-$5\times 5$ Conv.	$H\times W\times \left(C_I+2k\right)$	$H\times W\times k$
Concatenate	$H\times W\times \left(C_I+2k\right)$	$H\times W\times \left(C_I+3k\right)$
	$H\times W\times k$
BN-Relu-$5\times 5$ Conv.	$H\times W\times \left(C_I+3k\right)$	$H\times W\times k$
Concatenate	$H\times W\times \left(C_I+3k\right)$	$H\times W\times \left(C_I+4k\right)$
	$H\times W\times k$

Table 2. TABLE II: Parameters in the Neural Network

Stages	Layers	Output Size
		($k=16$)
Encoder-1	$3\times 3$ Conv. (strides=1)	$720\times 731\times 2k$
	-Relu-BN
	$3\times 3$ Conv. (strides=1)	$720\times 731\times 2k$
	-Relu-BN
Encoder-2	$3\times 3$ Conv. (strides=2)	$360\times 366\times 2k$
	DenseBlock-1	$360\times 366\times 6k$
Encoder-3	$3\times 3$ Conv. (strides=2)	$180\times 183\times 2k$
	DenseBlock-2	$180\times 183\times 6k$
Encoder-4	$3\times 3$ Conv. (strides=2)	$90\times 92\times 2k$
	DenseBlock-3	$90\times 92\times 6k$
Encoder-5	$3\times 3$ Conv. (strides=2)	$45\times 46\times 2k$
	DenseBlock-4	$45\times 46\times 6k$
Decoder-1	$3\times 3$ ConvTranspose (strides=2)	$90\times 92\times 6k$
	-Relu-BN
	Concatenate	$90\times 92\times 12k$
	$3\times 3$ Conv. (strides=1)	$90\times 92\times 6k$
	-Relu-BN
Decoder-2	$3\times 3$ ConvTranspose (strides=2)	$180\times 183\times 6k$
	-Relu-BN
	Concatenate	$180\times 183\times 12k$
	$3\times 3$ Conv. (strides=1)	$180\times 183\times 6k$
	-Relu-BN
Decoder-3	$3\times 3$ ConvTranspose (strides=2)	$360\times 366\times 6k$
	-Relu-BN
	Concatenate	$360\times 366\times 12k$
	$3\times 3$ Conv. (strides=1)	$360\times 366\times 6k$
	-Relu-BN
Decoder-4	$3\times 3$ ConvTranspose (strides=2)	$720\times 731\times 6k$
	-Relu-BN
	Concatenate	$720\times 731\times 8k$
	$3\times 3$ Conv. (strides=1)	$720\times 731\times 4k$
	-Relu-BN
Merge	$1\times 1$ Conv. (strides=1)	$720\times 731\times 1$
—	Add	$720\times 731\times 1$

Table 3. TABLE III: Quantitative measurement using different methods for this synthetic phantom.

	FBP	DLFBP	DNN-SS	SN-I	DD-Net
FSIM	0.841	0.995	0.994	0.872	0.986
IW-SSIM	0.866	0.996	0.996	0.892	0.986

Table 4. TABLE IV: Quantitative measurement using different methods for synthetic data sets.

	Methods	60	90	120	180
FSIM	FBP	0.718	0.790	0.831	0.897
	DLFBP	0.977	0.988	0.992	0.998
	DNN-SS	0.974	0.984	0.990	0.997
	SN-I	0.816	0.888	0.885	0.971
	DD-Net	0.953	0.976	0.982	0.993
IW-SSIM	FBP	0.678	0.788	0.848	0.938
	DLFBP	0.975	0.990	0.995	0.998
	DNN-SS	0.980	0.993	0.996	0.998
	SN-I	0.758	0.840	0.895	0.974
	DD-Net	0.949	0.981	0.988	0.995

Table 5. TABLE V: Quantitative measurement using different methods for the experimental data.

	FBP	DLFBP	DNN-SS	SN-I	DD-Net
FSIM	0.801	0.990	0.990	0.930	0.987
IW-SSIM	0.835	0.990	0.989	0.931	0.987

Table 6. TABLE VI: Quantitative measurement using different methods for experimental data sets.

	Methods	60	90	120	180
FSIM	FBP	0.712	0.782	0.802	0.904
	DLFBP	0.969	0.983	0.988	0.996
	DNN-SS	0.958	0.976	0.984	0.992
	SN-I	0.862	0.912	0.927	0.978
	DD-Net	0.965	0.978	0.983	0.989
IW-SSIM	FBP	0.705	0.793	0.832	0.926
	DLFBP	0.961	0.982	0.988	0.996
	DNN-SS	0.941	0.970	0.981	0.990
	SN-I	0.833	0.901	0.928	0.978
	DD-Net	0.957	0.978	0.985	0.990

Table 7. TABLE VII: Parameters and computational cost in DNN-SS and DLFBP.

	DNN-SS	DLFBP
Number of parameters (million)	37.50	1.36
Inference Time (seconds per sinogram)	3.56	0.21

Equations

$${\alpha }_{\theta }\left(s\right)\approx \frac{p\left(j+1,m\right)-p\left(j-1,m\right)}{2}$$

$$\delta \left(x,y\right)=\frac{1}{2}{\int }_0^{2\pi }U\times {\alpha }_{\theta }\left(s\right)\ast h\left(v\right)d\theta$$

$$h\left(v\right)=\frac{1}{2\pi }i\mathrm{sgn}\left(v\right)$$

$$P\left(\delta \left(x,y\right)\right)={\int }_0^{2\pi }\delta \left(x,y\right)dl$$

$$\hat{Y}=F(X)$$

$$\hat{Y}=H\left(X\right)+X$$

$$loss=\frac{1}{2}{\Vert Y-\hat{Y}\Vert }_2^2+1-MS\_SSIM(Y,\hat{Y})$$

$$I_n=\frac{I-mean(I)}{std(I)}$$

$$relI=\frac{Met-MetonFBP}{MetonFBP}$$

Full text

I Introduction

The invention of X-ray computed tomography (CT) has led to a revolution in many fields such as medical imaging, nondestructive testing and materials science. However, since traditional X-ray contrast is generated by the difference in attenuation, weak-absorbing materials are not imaged satisfactorily, limiting the range of application of X-ray CT. To overcome this problem, X-ray phase-contrast CT (PC-CT) uses the phase shift that X-rays undergo when passing through matter as the imaging signal to nondestructively provide the internal physical and biomedical properties of the specimens. Over the last years, several phase-contrast imaging techniques have been developed [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]. One of the recent developments is differential PC-CT (DPC-CT), based on a grating interferometer [12, 13, 14, 15, 16, 17, 18, 19, 20]. It has become more and more popular as a powerful analysis tool for soft-tissue and low-atomic-number samples.

Image reconstruction plays an important role for the development of DPC-CT. Filtered back-projection (FBP) with an imaginary Hilbert filter is generally preferred since it keeps a good balance between reconstrucion speed and image quality when applied to complete data. However, constrained by implementation conditions, DPC-CT with incomplete projections occurs quite frequently. The corresponding FBP reconstruction will have quite visible artifacts and noise.

Reconstruction with incomplete data has attracted more and more interests. Guna Kim et al investigated analytic CT reconstruction in sparse-angular sampling using a new interpolation method [25] to reduce patient radiation dosage. According to the results, the quality of the images reconstructed by their method was considerably improved over the cases using cubic interpolation method. Deriving a maximum likelihood (ML) reconstruction algorithm with regularization [26] for differential phase-contrast imaging, Thomas et al used spherically symmetric basis functions and differential footprints in forward and back-projection to avoid the need for numerical differentiation. Their results showed that sparsely sampled data could be handled efficiently. Based on a novel spling-based discretization of the forward model and an iterative reconstruction algorithm using the alternating direction method of multipliers, Masih et al gave out an iterative reconstruction method [27] for DPC-CT with fewer angular views, and the results suggested that their method allows to reduce the number of required views by a factor of four. Jian Fu et al developed an algebraic iteration reconstruction technique [28] for incomplete data DPC-CT. By minimizing the image total variation, their work could permits accurate tomographic imaging with less data. These reconstruction techniques could be better than FBP, but they still have some limits such as expensive time consumption for the succesive iterative steps and the complicated parameter selection.

A more recent trend is the application of deep learning (DL). It has led to a series of breakthroughs for image classification[29, 30] and segmentation[31] and also demonstrated impressive results on signal denoising[32] and artifacts reduction[33, 34].

There is currently a scarcity of researches on applying DL to DPC-CT reconstruction, and there remains an important need to develop the relative techniques. Nevertheless, several works have applied DL to absorption-based incomplete CT reconstruction. They can be grouped into two categories. The first category can be classified as post-processing. Its key idea is to reduce artifacts on the CT image domain. A method presented by Cierniak [35] uses hopfield-type neural network on the back-projection to solve the image deburring problem. This approach bypasses the problem of too many parameters by fixing them in the back-projection step and degenerates the reconstruction to an image-based filtering approach. Based on a persistent homology analysis, Han et al developed a deep learning residual architecture [36] for sparse-view CT reconstruction. The input of this architecture is the initial corrupted reconstruction image from FBP or other algorithm. It firstly estimates topologically simpler streaking artifacts from the input image and then subtracts the estimated result from the input image to get artifact-free image. Obviously this method is independent on X-CT reconstruction and works in an indirect way. Using multi-scale wavelet, they extended their work to limited angle CT reconstruction [37]. Jin et al also proposed a deep convolutional neural network [38] for inverse problem in imaging. It is similarly independent on X-CT reconstruction, but the estimated result is the final CT image not the artifacts. With dialted convolutions, Pelt et al introduced an architecture [39] to capture features at different image scales and densely connect all feature maps with each other. Their method is also independent on CT reconstruction, but is able to achieve accurate results with fewer parameters, which reduced the risk of overfitting the training data. Zhicheng Zhang et al proposed a method [40], which takes full advantages of DenseNet [41] and deconvolution [42] to remove streaking artifacts from the sparse-view CT images. The second category can be classified as sinogram completion. Its key idea is to complement the incomplete projections before reconstructing them with analytical algorithms. Hoyeon Lee et al developed an interpolation method [43] using convolutional neural network (CNN) to interpolate the missing data in sinogram from sparse-view CT. Their work shows better result than other interpolation methods, like linear interpolation method. Changed the network to U-Net [31] and combined with residual learning [44], their extended work [45] outperformed the existing interpolation methods and iterative image reconstruction approach. Donghoong Lee et al employed DL on hybrid domain of sparsely sampled CT to restore high quality images [46]. Firstly, they apply DL on the linear interpolated sinogram to get full sampled sinogram, then, DL is utilized again on the CT image reconstructed from the full sampled sinogram to obtain the final CT image. Their experiments shows it is able to restore images to a quality similar to fully sampled images.

In this paper, we report a new deep learning reconstruction framework for DPC-CT with incomplete projections. Different from post-processing and sinogram completion methods, we firstly project the initial reconstructed DPC-CT image into corrupted sinogram, in which the missing information will be completed in a fashion of CT scanning instead of interpolation, and then remove the artifacts on the projection sinogram domain. It is the tight coupling of the DL and FBP algorithm in the phase-contrast projection sinogram domain. The estimated result is the complete projection sinogram not the DPC-CT image or the artifacts. Embedded with DPC-CT reconstruction, it naturally encapsulates the physical imaging model of DPC-CT systems and is easy to be extended to deal with other challengs such as photon starvation and phase wrapping. When training, this framework firstly obtains the forward projections from the initial reconstruction by applying FBP to the original incomplete projection sinogram. Taking the complete projection sinogram as a target, they are then fed into the neural network to get the net parameters by deep learning. After training, this framework is determined and can reconstruct the final DPC-CT image for a given incomplete phase-contrast projection sinogram. Taking sparse-view DPC-CT reconstruction as an application example, this framework has been validated by using synthetic data sets and experimental data sets. This work is helpful to push the application of deep learning in the field of DPC-CT.

II Methods

II-A Framework overview

Fig.1 shows the proposed deep learning reconstruction framework for DPC-CT with incomplete projections. It is based on FBP and a neural network and called DLFBP. This framework consists of five parts. Initial FBP reconstruction of the original incomplete phase-contrast sinogram is the first part. The second part is the forward projection operator which is applied to the initial reconstructed image to obtain corrupted sinogram. Then the differential operation is conducted to get differential phase-contrast sinogram. The fourth part is a neural network and used to execute deep learning. Final FBP reconstruction of the complete phase-contrast sinogram from the fourth part is the last part of the framework.

A two-dimensional object can be described by a complex refractive index distribution $n\left({x,y}\right)=1-\delta\left({x,y}\right)+i\beta\left({x,y}\right)$, where $x$ and $y$ describe the coodinate system of the sample. In differential phase-contrast imaging, one measures the effect of variations of the real part $\delta$ by evaluating the tiny refraction angles of X-rays induced by the specimen with a grating Talbot-Lau interferometer. The differential phase-contrast projection can be expressed by ${\alpha_{\theta}}\left(s\right)=\partial\left({\int_{l}{\delta\left({x,y}\right)dl}}\right)/\partial s$, where $s$ describes the coordinate system of the detector, $\theta$ the rotation view angle of the object, and $l$ the incident ray direction. Using the three point differential method, the differential phase-contrast projection can be discretely re-expressed into Eq.(1).

[TABLE]

Here, $p$ is a two-dimensional matrix and each element in this matrix represents the line integral value detected by one detector channel. It is used in the third part to get differential phase-contrast sinogram.

Eq.(2) is the well-known two dimensional fan beam FBP algorithm which is adopted to reconstruct the DPC-CT image in this article. It has many extension versions for different CT scanning configurations. In this equation, $\delta\left({x,y}\right)$ represents the DPC-CT image, $U$ the geometrical weight factor, ${\alpha_{\theta}}\left(s\right)$ the differential phase-contrast projection sinogram, $h$ the Hilbert filter, Eq.(3), $v$ the frequency variant and $\theta$ the rotation angle.

[TABLE]

[TABLE]

Eq.(4) is the forward projection operator. In this equation, $P\left({\delta\left({x,y}\right)}\right)$ the forward projection, ${\delta\left({x,y}\right)}$ represents the initial reconstructed DPC-CT image and $l$ the forward projection path.

[TABLE]

Eq.(5) represents the way how the neural network in this framework produce complete phase-contrast sinogram. Corrupted sinogram $X$ with artifacts is feed into the neural network, then an artifact-free complete phase-contrast sinogram $\hat{Y}$ could be obtained. Finally, FBP is used again on $\hat{Y}$ to get the final DPC-CT image.

[TABLE]

II-B Neural network

The neural network is build upon a more elegant architecture, U-net [31], which has wide applications like semantic segmentation [31] and medical image processing [38], due to its goal is to generate a complete sinogram and the size of input and output are identical. U-net has large number of feature channels to extract high-level features in the extracting path, and large number of feature channels in the expansive path to allow the propagation of context information to higher resolution layers. As a consequence, its symmetric structure is organized in the form of Encoder-Decoder. Recently, research [47] shows many layers contribute very little and can in fact be randomly dropped during training, and DenseNet [41] points out densely-connection can act as a substitute for large number of feature channels because there is no need to relearn redundant feature-maps. Therefore, this neural network takes fully advantages of U-net and DenseNet to generate complete phase-contrast sinogram.

According to the spatial resolution of feature maps, this network can be divided into different stages. As depicted in Fig.1, there are five encoding stages which are surrounded by gray dotted boxes, and four decoding stages surrounded by green dotted boxes. Numbers next to each stage denote the corresponding spatial resolution on that stage.

II-B1 Encoder-Decoder

Encoder is used to extract multi scale features from the input corrupted sinogram. Bigger convolutional (Conv.) filter could provide larger receptive field, but more parameters will be introduced, which may cause overfitting and hardly converging. Consecutive convolution with small filter could provide the same size of receptive field as bigger filter does, with less parameters. Therefore, consecutive $3\times 3$ Conv., Rectified Linear Units (Relu) [48] and BatchNormalization (BN)[49] are used to extract features from the input.

In common use [31, 40], max-pooling is taken to reduce the feature maps’ spatial resolution to increase the size of receptive field and get multi scale features. It computes fastly, but only keeps the maximum value from the pooling window. In order to pay more attention to important values rather than maximum one, strided convolution is used. The parameters in convolution filter will decide which part is important and give it higher weight. Besides, to reduce the feature maps’ spatial resolution as max-pooling does, the stride is set to 2. So $3\times 3$ Conv. (strides=2) is used for reducing spatial resolution.

After high-level features are extracted, which possess high semantic information and are abstract, Decoder should make use of those abstract features to restore a high-quality output with the same size as the input corrupted sinogram. Firstly, deconvolution [42], which can be regarded as the reverse version of convolution, also known as transpose convolution (ConvTranspose), is adopted to recover the spatial resolution, with stride 2. Then high-level features are concatenated with low-level ones, and $3\times 3$ Conv.-Relu-BN are used to refine details. It’s working stage-by-stage. When the feature maps’ spatial resolution are the same as input’s, $1\times 1$ Conv. is used to merge multi channels into one for matching the input’s dimension.

II-B2 Skip Connection

There exists two kinds of skip connections: connection between input and output in the form of residue learning [44] and connections among different stages.

Firstly, input is directly added to the output, so the neural network could be reformulated as (6). Consequently, it only needs to recognize and remove the artifacts from the corrupted sinogram, rather than removing artifacts and building up the whole complete sinogram from those abstract features at the same time. This way can simplify the learning process, which makes the network focus on the artifact reduction.

[TABLE]

Secondly, the output from different stages are connected. Though high-level features contains high semantic information, it is not helpful for recovering the spatial information, while low-level features retain spatial accuracy with less semantic information. So low-level features are concatenated with high-level ones to refine the spatial information, which is similar to [31, 38, 45, 40].

II-B3 DenseBlock

DenseBlock is the basic module in DenseNet [41]. It distills skip connection into a simple pattern: all layers are directly connected with each other in the same block. As shown in Fig.2, each layer obtains additional inputs from all preceding layers and passes on its own feature maps to all subsequent layers. It exploits the potential of the network through feature reuse instead of extreamly deep or wide architectures, and yield condensed models and highly parameter-efficiency.

Densely connection makes model easy to train, because each layer in the same block has direct access to the gradients from the loss functions and the original input signal, resulting in implicit deep supervision [50]. Further, it proves that small number of feature channels are sufficient to obtain state-of-the-art results, because global information can be reached from everywhere within the network, and unlike in traditional network architectures, there is no need to replicate it from layer to layer. Besides, each layer just adds $k\left({e.g.\,\,k=16}\right)$ feature maps to global information, which regulates how much new information each layer contributes to the global one and reduces redundancy. This pattern also has regularizing effect resulting from highly-parameter-efficiency, which reduces overfitting on tasks with smaller training set sizes.

In this paper, each DenseBlock contains four BN-Relu-$5\times 5$ Conv. (strides=1) layers, and each layer produces $k$ new feature maps. The parametric architecture of one DenseBlock is shown in Table I, where $H$ and $W$ denotes spatial resolution of feature maps’, ${C_{I}}$ denotes the number of input feature channels.

II-B4 Parameters in Neural Network

The whole architecture of this neural network is shown in Table II. All the parameters in this neural network are initialized using a Gaussian distribution with zero mean and standard deviation $\sqrt{\frac{2}{{{n_{in}}}}}$ in which ${n_{in}}$ indicates the number of input units in each layer.

Mean square error (MSE) can force the network to learn difference on pixel value, and Multi-scale structure similarity (MS-SSIM) [51] can make it to learn the difference on structure. In this paper, MSE and MS-SSIM are combined to be the loss function, Eq.(7). Here, $\hat{Y}$ denotes output, and $Y$ represents corresponding grouth truth. The initial learning rate is set to $1\times{10^{-4}}$, and gradually reduced to $1\times{10^{-5}}$.

[TABLE]

II-C Artifact reduction in the sinogram domain

The output of a neural network can not be entirely correct, and even the latest optimization algorithm [52] can not improve the accuracy to 100%. When DL is applied in the CT image domain, the imprecise value of the output will be directly reflected in the CT image, resulting in the structural deformation. However, Eq. (2) tells that when leveraging FBP on sinogram to obtain CT image, the object information on CT image is obtained by weighted combination of information from all angles on sinogram. Thus, the final CT image reconstructed from the sinogram in which the artifact is reduced by DL could be more tolerant to inaccurate output value, and could get more accurate structure. The experimental section also confirms that the results of artifact reduction conducted in the sinogram domain are better than those whose artifact reduction is conducted on CT image domain.

II-D Running modes

This framework has two running modes. One is training mode and another is working mode. Training mode has following steps:

A set of incomplete phase-contrast sinograms is firstly matched with the corresponding complete sinogram into many pairs of training data which include an incomplete sinogram and a complete sinograms. 2. 2.

These data is fed into the framework depicted in Fig.1 one pair by one pair and the network parameters are updated iteratively. 3. 3.

When all pairs are used once, an outer learning iteration is completed. 4. 4.

Repeat steps 2 and 3 utill the learning converges.

The procedure for working mode is much simpler. When an incomplete phase-contrast sinogram is fed into the framework determined by training mode, the output of the framework will be a high quality DPC-CT image.

III Experiments

Sparse-view DPC-CT is a typical case with incomplete data. Taking sparse-view DPC-CT as an application example, this section validates the proposed reconstruction framework by synthetic data sets and experimental data sets.

III-A Data preparation

III-A1 Synthetic

For synthetic data sets, 500 phantoms are used to obtain the fan beam phase-contrast sinograms with different sampling factors. Each phantom consists of tens of ellipses with random refraction coefficients, size, and location. And the sampling factors are set to be 1, 4, 6, 8 and 12. They correspond to 720, 180, 120, 90 and 60 views, respectively. The size of each phantom is $512\times 512$ pixel. Fan beam sinograms are generated by using the embedded MATLAB function $fanbeam\left({}\right)$. The width of all the sinograms are 731 pixels. The sinogram with sampling factor 1 has a size $720\times 731$ pixel and is treated as complete one. Other sinograms are incomplete. Sinograms of 400 phantoms are used to train the framework and those from another 100 phantoms are used to test the framework.

Within the framework, for each incomplete sinogram, the initial FBP reconstruction is firstly executed with Eq.(2) and Eq.(3) to obtain the initial DPC-CT image. Then the forward projection operator in Eq.(4) and the differential method in Eq.(1) are applied to the initial DPC-CT image to generate the corresponding corrupted phase-contrast sinogram with a size $720\times 731$ pixel. Next the iterative deep learning runs to update the network parameters by making comparison between the corrupted and the complete phase-contrast sinogram.

This framework involves complicated calculations such as image reconstruction, forward projection and convolution. They may lead to unignored computation errors and degrade the training efficiency and accuracy. So we apply normalization operation to input $X$ and ground truth $Y$ to avoid this problem. This normalization operation is expressed in Eq. (8) in which $I_{n}$ represents the normalized image, $I$ the raw image, $mean()$ a operator to obtain mean value and $std()$ a operator to obtain standard deviation value.

[TABLE]

III-A2 Experimental

The experimental CT Lymph Nodes data sets [53, 54, 55] is collected from The Cancer Imaging Achive (TCIA) [56]. This collection consists of CT images of 90 patients’ mediastinum and 86 patients’ abdomen. For the consistency of the data sets’ size on both synthetic and experimental, 50 out of 86 patients who underwent abdominal imaging are randomly selected, and then split into 40 patients used for training and 10 patients for testing. 10 images are randomly chosen from each of these 50 patients.

After that, all the operations and procedures are the same as the ones for synthetic data sets. Fig. 3 is an example of how to prepare the training and testing data.

III-B Implementation

This framework is implemented with Python 3.5.2 and Tensorflow 1.8. It runs in a workstation Advantech AIMB-785 with a CPU i7 6700 and a Graphics Processing Unit (GPU) nVidia GTX 1080Ti 11 GBytes.

For synthetic data sets, datas from 60, 90, 120 and 180 views are used to train a model seperately. So experimental data sets do.

Adam [57] algorithm is used to train the neural network. The mini-batch size is 2. All the models are trained for 100 epochs.

III-C Comparison with other methods

In this study, two methods (deep-neural-network-enabled sinogram synthesis, DNN-SS [45] and sinogram-normalization interpolation, SN-I [25]) based on sinogram completion and one method (DenseNet-Deconvolution network, DD-Net [40]) based on post-processing are conducted to compare the performance with the proposed framework.

III-D Image evaluation

When training, MSE and MS-SSIM have been used as the loss function to guide the network to optimize parameters. Consequently, the network will try to generate output that has low MSE and high MS-SSIM index. So MSE and MS-SSIM are not objective and not suitable for quantitative measurement of the network performance. Besides, Lin Zhang et al [58] compared the performance of several full reference image evaluation methods, such as feature similarity (FSIM) [59], MS-SSIM, SSIM [60] and peak signal to noise ratio (PSNR), and results show that FSIM and information content weighted SSIM index (IW-SSIM) [61] are more accurate than others. Therefore, FSIM and IW-SSIM are used as quantitative metrics.

Qualitative evaluation is also carried out using visual inspections and image intensity profiles.

Relative improvement ratio based on FSIM and IW-SSIM index using different methods compared to FBP is calculated with Eq.(9).

[TABLE]

Where $relI$ represent relative improvement ratio, $Met$ denotes metric value using FSIM and IW-SSIM on different methods. $MetonFBP$ denotes metric value using FSIM and IW-SSIM on FBP images.

III-E Results

III-E1 Synthetic data sets

Figs.4-6 present results of one of the 100 synthetic phantoms with 120 views for testing using different methods. Same regions indicated with yellow box are enlarged for better visualization in Fig.5. The image intensity profile in same position is shown in Fig.6.

As expected, severe artifacts exist in the result using FBP, and DLFBP, DNN-SS, SN-I and DD-Net could reduce the artifacts. Though most artifacts are suppressed, there still remains some in the background using SN-I, furthermore, its central part of the image has the clearest structure and the nearer the boundary, the more blurred it is. DD-Net could also remove the artifacts satisfactorily, but the result is overly smoothed, and there is slight distortion in the edge structure. As shown in Fig.5 and Fig.6, DD-Net loses small image details and the pixel value changes gently and far from the reference value. Through visual inspection, DLFBP and DD-Net could remove the artifacts and retain small structures.

Table III shows the quantitative measurement using FSIM and IW-SSIM index on this synthetic phantom. DLFBP gets higher value on FSIM and IW-SSIM index compared to other methods.

Fig.7 shows the average relative improvement ratio based on FSIM and IW-SSIM index using different methods compared to the results from FBP on the whole 100 synthetic testing datas. Tabel IV lists the average FSIM and IW-SSIM values using different methods on 60, 90, 120, and 180 views’ synthetic data sets. DLFBP and DNN-SS get higher metric values than other methods.

III-E2 Experimental data sets

Figs.8-Fig.10 present results of one of the 100 experimental abdominal slice with 120 views for testing using different methods. Same regions indicated with yellow box are enlarged for better visualization in Fig.9. The image intensity profile in same position is shown in Fig.10.

Still, SN-I could not completely remove the artifacts, and its edge is blurred. As shown in Fig.9 and Fig.10, DLFBP, DNN-SS and DD-Net could suppress the artifacts clearly, while the result from DD-Net is overly smoothed and loses details.

Table V shows the quantitative measurement using FSIM and IW-SSIM index on this experimental slice. DLFBP and DNN-SS get higher value on FSIM and IW-SSIM index compared to other methods.

Fig.11 shows the average relative improvement ratio based on FSIM and IW-SSIM index using different methods compared to the results from FBP on the entire experimental testing data sets. Tabel VI lists the average FSIM and IW-SSIM values using different methods on 60, 90, 120, and 180 views’ experimental testing data sets. DLFBP and DNN-SS get higher metric values than other methods.

III-E3 Parameters and computational cost

The synthetic and experimental results show DLFBP and DNN-SS outperform other methods in artifacts reduction and structure preservation. Using FSIM and IW-SSIM index as quantitative measurement, these two get high value than others. Besides, DLFBP and DNN-SS both corrected the missing information in the phase-contrast sinogram domain. But the number of parameters used in DNN-SS is 27 times that of DLFBP. And when inference with a single sinogram, DNN-SS will spend more time than DLFBP. Table VII shows the number of parameters in DNN-SS and DLFBP, and the inference time using a single sinogram.

IV Conclusion

In this paper, we reported a new deep learning reconstruction framework for DPC-CT with incomplete data. Different from post-processing and sinogram completion methods, the artifacts are firstly projected into corrupted sinogram, in which the missing information will be completed in a fasion of CT scanning instead of numerical interpolation, and then it is reduced by a neural network. Validated with synthetic and experimental data sets, it could remove artifacts clearly and preserve structures at the same time.

Compared with sinogram completion and post-processing methods and taking FSIM and IW-SSIM index as quantitative metrics, the results from which missing information are corrected in the phase-contrast sinogram domain are better than those from post-processing method. Post-processing will lose tiny details and deform structures, which may be caused by the imprecise output from neural network. When artifacts reduction is conducted in the sinogram domain, it is more tolerant to the inaccurate output from the neural network, because the objects are weighted combination of all projections in the sinogram.

In addition, this framework has high computational efficiency. The neural network is constructed based on elegant models, U-net and DenseNet. U-net could extract multi-scale features from the input and then make use of those abstract features to generate complete sinogram. But there are large number of feature channels in original U-net, and DenseNet has demonstrated small number of feature channels is sufficient to achieve the same effect through densely connections. By combination of U-net and DenseNet, the reported framework could get slightly better results while the number of parameters is 27 times less than that of the method which also corrected information in the sinogram domain.

Acknowledgment

We acknowledge support from National Natural Science Foundation of China (12014003), the large scale science facilities joint funed by National Natural Science Foundation of China and Chinese Academy of Science (U1432101, 11179009).