|
1.IntroductionThe ability of spectral computed tomography (CT) to distinguish photons in terms of their different energies allows this technique to decompose scanned objects into their basis materials.1–5 Compared with linear attenuation coefficients from traditional single energy CT, the distribution of basis materials can effectively reflect inner information. For example, in clinical practice, liver fat quantification requires the calculation of the components of fat, tissue, blood, and contrast agent.6 With the upgrading of product applications and the popularity of spectral CT, many material decomposition algorithms have been proposed. These algorithms can be divided into three types, namely, projection-domain processing methods, image-domain processing methods, and direct iterative reconstruction methods. Projection-domain methods use the projections from detectors as the input data and bring the decomposition model into nonlinear polychromatic projection process.7–9 Therefore, they are sensitive to the mismatch of the projection data in different energy-bins. Image-domain methods assume that the effective attenuation coefficients of each pixel of the spectral CT images can be decomposed into the linear combination of several basis functions.10–16 They are widely using in current preclinical and clinical applications. However, the performance of image-domain method is usually affected by the ring artifacts and beam-hardening artifacts within the CT images especially for the spectral CT using the photon-counting detectors. The direct iterative reconstruction method incorporates the models of the material decomposition and the physics of spectral CT transmission.17–20 For example, in Ref. 20, the iterative reconstruction method can even perform MMD of traditional integral detectors using angle-dependent filters. Theoretically, it has the potentials of decreasing noise and artifacts by taking the advantages of accurate modeling. However, its computational cost of the forward- and backward-projections between the material images and the projections in different energy-bins is expensive. In this paper, we focus on the MMD problem in image-domain. Traditional material decomposition algorithms in spectral CT use the linear combination of a set of basis functions to represent the linear attenuation coefficient of location and then apply these functions in material decomposition.21 Different basis functions lead to different material decomposition methods. For example, Ref. 22 uses the linear attenuation curves of preselected materials as the basis functions, leading an analytical solution of MMD problem. However, regardless of the choice of basis functions, these traditional material decomposition algorithms are voxelwise based, using no more than dozens of pixels around the target pixel. Therefore, their performances usually depend on the reconstructed image quality of spectral CT. They are often affected by the beam-hardening artifacts and ring artifacts in the spectral CT using the photon-counting detectors. An idea to avoid these problems is using the data-driven methods, for example, machine learning, instead of the analytical solution such as Ref. 22. In this way, we can expand the receptive field of data-driven algorithms to provide more information for their learning. Wang23 published an outlook article on deep learning and presented a roadmap of CT imaging methods in the deep learning framework. Inspired by his work, we imported the deep learning technique into spectral CT and proposed a multimaterial decomposition (MMD) method in image-domain.24 Deep learning has made considerable progress in many fields. In image analyses, deep neural networks can extract the features of images without manual marking and are thus effective in image classification.25–27 Recently, deep learning has shown its effectiveness in medical imaging,28 such as auxiliary diagnosis,29–31 lesion localization,32,33 and image noise reduction.34,35 By using deep learning, we can break through the restriction of sampling theorem and realize reconstruction using few views.36 Using deep learning to acquire priori knowledge of images can also enhance the results of traditional MMD algorithm.37 As a universal approximator,38 neural networks can be trained to solve the MMD problem by given adequate and suitable training data. In this paper, we propose a convolutional material decomposition (CMD) algorithm on the basis of patches instead of voxels of the spectral CT images by building two convolutional neural networks (CNN) based on existing networks. Section 2 presents the CMD algorithm. Then, the results of simulation and physical experiment are presented in Sec. 3 to validate it. Sections 4 and 5 are the discussion and conclusion. 2.Method2.1.MMD ModelAs an image-domain method, the reconstructed image of spectral CT is used as the input data. The reconstructed image is the x-ray linear attenuation coefficient of an object scanned at different spectra. Here, is the spatial position, and is the spectrum or energy-bin.35 Knowing by the spectral CT reconstruction algorithm, we choose the attenuation coefficient of the basis material at different spectra as the basis function. Thus, we have where denotes the volume fraction of basis material at position , is the number of basis materials, and is the number of spectra. As proven in Ref. 4, we apply the constraints on volume and mass conservation, all of which lead to the following constraints:The constraints on each have the same form as the result of the image classification, that is, the possibility of each image type. Therefore, a CNN designed for image classification can be transformed to fit the material decomposition problem. 2.2.Deep Learning Based on CNNEquations (1) and (2) describe the relationship between the reconstructed image and the decomposition result. Traditional material decomposition algorithms in image-domain are based on these voxelwise equations. They often suffer from the spatial-specific errors or artifacts such as beam-hardening artifacts, metal artifacts, ring artifacts, and so on. Although some decomposition algorithms add regularization terms based on a priori knowledge, e.g., total variation, into the cost function to improve image quality, the results are still unsatisfactory. These spatial-specific features are difficult to describe using explicit mathematical formulas; nevertheless, they can be learned using deep neural networks. As a universal approximator, neural networks have been proved effective if enough neurons are given.38 This paper will show that neural networks can determine the patterns of material decomposition if the neural networks are trained with sufficient samples. Enlarging the receptive field of each pixel can improve material decomposition performance as well. As shown in Fig. 1, we train the CNN to fit each patch and the exact decomposition result of this patch’s central voxel. Then, by moving this receptive field voxel-by-voxel among the reconstructed CT image, we obtain the decomposition result. We refer to this method as the CMD algorithm. CNN has become one of the most important methods of image analysis.26,38 A CNN consists of several simple operations or so-called layers. By sending an image into the CNN and going through several layers, such as convolutional, activation, pooling, and fully connected layers, we can transform the image into the result we want, including classification or material decomposition. The image can be considered as a three-rank tensor , whose first two parameters and describe the position of a voxel and is the channel.39 For a chromatic image, ranges from 1 to 3, and for the reconstructed image of spectral CT, ranges from 1 to . All the layers mentioned above are the derivable operations on this tensor, transforming tensor into another tensor. Because these parameters are derivable, most of parameters in the layers can be trained automatically by stochastic gradient descent (SGD) or other training algorithms.38 Details of these layers can be found in the Appendix. 2.3.CNN Used in CMDWe build two CNNs for the decomposition task. Both CNNs are simplified versions of the CNNs that work effectively in image classification. Let and . is the size of the input patch , which is the receptive field with size , and channel number . 2.3.1.Visual geometry group networkThe first CNN is based on a visual geometry group (VGG) network.25 The VGG network is stacked by convolution layers with small kernel. Its structure is shown in Fig. 2. All the activation layers are ReLUs except for the last one, which is Softmax. Fig. 2Structure of the first CNN based on VGG used in this work. The left parts of the rectangles are the layers’ names, and the right parts are the sizes of the input tensors and output tensors . The numbers after the convolutional layers are , which refer to the size of the convolutional kernels. Parameter of convolutional layers is determined by the input tensor. ![]() 2.3.2.Deep residual networkTo compare the performances of the different networks, we build the second CNN based on a deep residual network (DRN),27 which is more complicated than a VGG network. The network comprises three types of blocks that are shown in Fig. 3. Block A does not change the size of the input tensor. Block B halves the size of the input image by downsampling the first two parameters. Block C halves the size while doubling its channels. The full structure of the DRN is shown in Fig. 4. 2.4.Loss Function and SGD with MomentumAs shown in Fig. 1, is the training sample, where is a three-rank tensor and is a one-rank tensor. represents the material decomposition result of the central voxel of . The CNN is initialized randomly and trained using the minibatch SGD method. To establish a regression problem, we use mean square error (MSE) as the loss value, where is the output of the CNN. is a set of , who is the minibatch in the training step, the validation set in validation step, and test set in test step. is the number of tensors in . As defined above, every parameter in the CNN layers is derivable or assigned with a gradient. According to the chain rule of derivation,38 the partial derivatives of to the parameters in the CNN are calculable. For a parameter of the CNN, the updating function is where lr is the learning rate or the step length of the gradient descent algorithm. To improve the stability of the convergence, we add the momentum term to Eq. (4)However, no specific method can be used for training CNNs. In this study, a learning rate decrease policy is adopted to train the CNN and achieve satisfactory convergence. Initially, we set and patience = P and then pick VAL_NUM patches randomly from the validation phantoms as the validation set. We then pick BATCH_NUM samples from the training set randomly and use them to train the CNN through SGD. For every ITER_NUM iteration (the training patches are reselected for every iteration), we send the validation set to calculate the MSE of the CNN output. If the MSE is greater than the lastMSE for more than patience times sequentially, we divide lr into two, add P_STEP to patience, repick the validation set, and continue the iteration until the MAX_ITERATION_TIMES or while lr is less than . The pseudocode is shown in Fig. 5. 2.5.Convergence Property of Neural NetworkThe SGD with momentum method is widely used in deep learning optimization and turns out to be effective among various applications. In the case of a local minimum nearing . At time step , Eq. (5) can be rewritten as Define momentum operator Equation (11) is equivalent to Reference 40 proved that when is close to the local minimum , if we may get where is the spectral radius of . When is a constant matrix, which ensures Eq. (8) converges. Actually, may be changed in the above iterations. In this case, Ref. 40 verified that this momentum method could robustly yield linear convergence with rate empirically around a local minimum .The difference between the method in Ref. 40 and ours is that we evaluate the gradient by sampling a minibatch from the training set instead of calculating over the whole training set, which is a widely used trick to accelerate the convergence without losing too much accuracy. In this way, the SGD with momentum method is converged in the expectation sense. However, because the fitting ability of deep neural network is nearly infinite, the search for the global optimal solution in most cases will lead to overfitting.41 On the other hand, critical point is quite rare in a deep neural network. The SGD method tends to stop in some “flat minima.”42 The loss of flat minima that SGD could reach is related with the structure of neural network.43 As stated in the paper, the CNNs we used are based on VGG and DRN, which are widely verified in image classification. As a task immigration, these CNNs perform well according to our experiments. 3.Numerical Simulation and Physical Experiment ResultsTo train the CNN to solve the MMD problem, we need organize three datasets, namely, the training dataset for SGD to optimize the parameters in the CNNs, the validation dataset for adjusting the hyperparameters (such as the number of mini-batches, lr, and so on), and the test dataset for evaluating the performance of the VGG-CMD and DRN-CMD. The three datasets should satisfy the same distribution.38 The ground truth values for the three datasets are also necessary as the supervised signal for training or the evaluation criterion. The phantom generation, projection, and reconstruction are all coded in MATLAB R2014a. The construction, training, and testing of CNN are coded by deep learning library Keras 2 in Python.44 3.1.Numerical Simulation ResultsIn numerical simulations, we used 115 digital phantoms based on the 2-D Shepp–Logan phantom to train and evaluate the proposed CMD algorithm. All of these phantoms comprise ten ellipses and five materials such as soft tissue, lung, bone, blood, and air. Their linear attenuation coefficients are obtained from Ref. 5, as shown in Table 1. One of the phantoms is shown in Fig. 6. The concentrations of these ellipses are shown in Table 2. “×” in Table 2 means the concentration is 0. “∼” means the concentration is 1 minus the other materials’ concentrations for meeting the constraints of Eq. (2). “” means that the concentration of that ellipse satisfies 0.1 normal distribution among different phantoms. Then, we randomly changed the position, size, and direction of the ellipses, and the concentrations of the materials to generate other similar 114 phantoms, as shown in Fig. 7. Table 1ID of materials used this paper and in Ref. 5.
Fig. 6An example of the digital phantoms. (a) One of the phantoms who comprises 10 ellipses. (b)–(f) The concentrations of five materials: soft tissue, lung, bone, blood, and air, respectively. ![]() Table 2Mean concentrations of materials in each ellipse of the generated phantoms.
Fig. 7Some phantoms used in the simulation. The size, position, and direction of every ellipse and the concentrations of the five materials are different (but close). ![]() Note that the linear attenuation coefficients of materials 1, 4 and 2, 5 in this paper are both relatively close, which worsens the ill-condition of MMD problem. The x-ray beam spectra used in numerical simulations were generated with a Siemens simulator45 at 75, 135, 105, and 95 kilovolts peak (kVp), respectively. A 12-mm thickness aluminum filter was added. To simulate five-energy-bin photon-counting detectors, each spectrum is divided into five energy-bins as shown in Fig. 8. The x-ray fan-beam covered a 20-mm-diameter field of view with 256 detector elements. The number of views over 360 deg is 360. The reconstructed images were discretized into a grid. We use the method in Ref. 46 to obtain the projection, which can simulate the statistical fluctuation and scattering in the imaging process, and can simulate the unsatisfactory situation such as pile up and charge sharing of the detectors based on real detectors. The spectral CT images in different energy-bins were reconstructed with the ASD-POCS algorithm with 25 iterations.47 Fig. 8In numerical simulations, four spectra were generated using a Siemens simulator. Each spectrum was divided into five different energy-bins. Twenty energy-bins of the four spectra used in this paper. The maximum energies of (a–d) are 75, 135, 105, and 95 kVp, respectively. The scales of the energy-bins are marked in each figure. ![]() Two different simulation studies were carried out. The first one only used the reconstructed spectral CT images of all of 115 phantoms with the 75 kVp spectrum of Fig. 8(a). We divided them randomly into three parts with the number of 85, 15, and 15, and used them to generate image patches, which were used as training sets, verification sets, and test sets, respectively. The patches we used has . As stated in Sec. 2.4, we set , , P_STEP = 0.5, VAL_NUM = 16384, BATCH_NUM = 128, ITER_NUM = 64, , and . The MSE of decomposition results of VGG-CMD and DRN-CMD algorithms are shown in Fig. 9 for the 15 test phantoms. As a contrast, we used an extended direct inversion algorithm in image-domain algorithm to calculate the MMD voxel-by-voxel (EDI-MMD). It can be seen as a spectral CT version of the original one proposed by Mendonca et al.4,22 The procedure of EDI-MMD is shown in Fig. 10. The basis material triplets set is in Table 3. Table 3Triplets set for phantom in simulation.
On our GPU (GTX 1070), training VGG-CMD and DRN-CMD took about 3 and 12 h, respectively. After separation, the average MSE of all the 15 testing phantoms are for VGG-CMD, for DRN-CMD, and for EDI-MMD, respectively. The MSE of EDI-MMD algorithm of different testing phantoms are very close because the illcondition of MMD problems of the test phantoms are similar with the same five basis materials as shown in Figs. 6(b)–6(f). There are some differences among the MSE of VGG-CMD and DRN-CMD algorithms for different phantoms due to their different structures. However, compared to the results of EDI-MMD, the MSE are reduced more than two orders in general. The decomposition results of one phantom are shown in Fig. 11. The rows from top to bottom are the decomposition ground truth and the results of VGG-CMD, DRN-CMD, and EDI-MMD, respectively. The columns from left to right are the concentrations of five materials: soft tissue, lung, bone, blood, and air, respectively. The EDI-MMD algorithm decomposes the soft tissue, lung, and blood wrong due to the close linear attenuation coefficients. Fig. 11Ground truth and decomposition results of three algorithms. The rows from top to bottom are ground truth, results of VGG-CMD, DRN-CMD, and EDI-MMD algorithms. The columns from left to right are the concentrations of soft tissue, lung, bone, blood, and air. ![]() We chose the 44th, 69th, and 102nd rows of the decomposition result of lung to draw profiles, which divide the ellipses halves. The positions of these profiles are showed in the second column of Fig. 11, and the section lines are drawn in Fig. 12. The MSE of the three methods on the profiles are presented in Table 4. Fig. 12Profiles of the decomposition results of lung. (a–c) Mean the different positions of the 44th, 69th, and 102nd rows. ![]() Table 4MSE of the decomposition section lines of material 2.
Equation (1) indicates that the decomposition result should only depend on single voxel if we ignore the unideal factors in spectral CT image. Although it uses patches as the input data, the output of CMD algorithm should rely mainly on the central voxel of every patch. We chose a test phantom. A square in the center of each patch of its reconstructed images is set to zero. The decomposition results are shown in Fig. 13. Although the zero region only occupies of the whole patch, the degradation of the output of the CMD is drastic. This result indicates that the deep learning procedure of CMD algorithm can insure that the central voxel of each patch plays a more important role. Fig. 13Decomposition results of (a) VGG-CMD and (b) DRN-CMD. The top row of each figure is the results with normal patches. The bottom row of each figure is the results with the patches setting the central square to zero. The columns from left to right are the concentrations of soft tissue, lung, bone, blood, and air. ![]() In the second simulation, we studied the CMD performance on the spectral CT images from an unknown or never trained spectrum. Traditional material decomposition algorithms need to measure of the basis materials and spectral CT images of the scanned object at the same x-ray spectrum. relates to not only the physical properties of material itself but also the spectrum of incident x-ray beam and the detection efficiency of the detector. Therefore, their decomposition results are sensitive to the variation of the effective spectrum of spectral CT systems. Figure 14 shows the decomposition results of EDI-MMD algorithm. The first, second, and third rows are the results when the x-ray spectra and of the basis materials within the corresponding energy-bins are all exactly provided. If the x-ray spectrum used in CT scanning is changed from the one measuring of the basis materials, the decomposition results of EDI-MMD will deteriorate. As showed in the last row of Fig. 14, of the basis materials used in EDI-MMD algorithm were generated by performing interpolations among the of the basis materials of the spectra of Figs. 8(a)–8(c). The average MSE of the results of spectra Figs. 8(a)–8(c) is , while the MSE of the results of spectrum Fig. 8(d) is . Fig. 14Decomposition results of EDI-MMD algorithm. The rows from top to bottom are the results using spectrum (a), (b), (c), and (d) of Fig. 8, respectively. The columns from left to right are the concentrations of soft tissue, lung, bone, blood, and air. ![]() Trained by patches from multiple different spectra rather than one specific spectra, CMD shows its potential to solve the MMD problem of unknown x-ray spectra. Because training for every CT system and its configuration on kVp/filtering is sometimes uneconomical. In this simulation, we randomly picked 85 phantoms and reconstructed their spectral CT images at the 15 energy-bins of Figs. 8(a)–8(c). Then, from each of these 15 channels of patch, five channels were randomly selected to form the input of CNN. That is, the input tensor still had . The other 15 phantoms were randomly picked to generate the validation dataset in exactly the same way. We use the same learning rate decrease policy as in the first simulation to train the CNNs. As the test dataset, the reconstructed images of the rest 15 phantoms in all of 20 energy-bins of Fig. 8(a)–8(d) were used to test the decomposition results of the CMD algorithm. The results are shown in Fig. 15. For the first three spectra, the average MSE on 15 testing phantoms are for VGG-CMD and for DRN-CMD, respectively. And for the fourth spectra, which was not used in training, the MSE of the above two CMD algorithms are 9.7165e-4 and , respectively. As shown in Fig. 15, the CMDs’ performance did not turn worse when the x-ray spectrum was changed from the one calibrating the basis materials. This result means that we can use the proposed CMD algorithm to solve the MMD problems in different spectral CT systems with different hardware configurations even if the spectrum of a CT system is unknown. 3.2.Physical Experiment ResultsIn the physical experiment, we used the reconstructed images of our spectral CT system to validate CMD algorithm and compare it with the EDI-MMD algorithm. The experiments were completed on our spectral CT system as shown in Fig. 16. It uses a conventional Hamamatsu L12161-07 x-ray tube and a linear photon-counting detector array (eV3500, eV PRODUCTS, Saxonburg, Pennsylvania).46,48 The detector array has 256 elements and 5 energy-bins of each one. The energy-bins were set at [26, 33], [33, 40], [40, 50], [50, 60], and [60, 80] keV. The size of each element is . The x-ray tube’s max energy was set 75 kVp, and the tube current . Over a 360-deg scan, 360 projections were acquired, and every projection radiated for 5 s. The distance between the x-ray source and the detector was 70 cm, and the distance between the x-ray source and the center of rotation was 44 cm. The reconstructed image had whose size was . The image was reconstructed by the filtered-backprojection algorithm for each energy-bin. The phantom used in the experiments was a 50-mm diameter resinic cylinder with eight 11-mm-diameter round holes to insert centrifugal tubes, as shown in Fig. 17. We used three types of solutions, namely, NaI, , and , with different concentrations to test the performance of the CMD algorithm and compare it with EDI-MMD algorithm. Four basis materials were used. NaI and represented the typical components of the contrast agents used in clinics, and represented the typical components of the human body. According to the concentrations in Refs. 11 and 49, all the concentrations (weight/volume) of the three solutes we used are listed in Table 5. In each spectral CT scanning, we used eight solutions with different concentrations in one row of Table 5. After six scans, we obtained reconstructed images of different solutions and energy-bins. Fig. 17(a) All of the 48 centrifugal tubes filled with different concentration solutions. (b) The resinic cylinder phantom. ![]() Table 5Concentrations in centrifugal tubes.
Figures 18(a)–18(f) show the six reconstructed images in the first energy-bin of [26, 33] keV, and Figs. 18(g)–18(k) are the reconstructed images of the first phantom in all five energy-bins, whose display window are . It shows that the images at higher energy-bins are darker than the ones at lower energy-bins, which means high energy x-ray photons were less absorbed than low energy x-ray photons. The diameter of centrifuge tubes contains 36 voxels in reconstructed images. Note that there are some ring artifacts in the images which were caused by the different detection response functions element-by-element, especially the abnormally high detection efficiency of the elements on the edges of the 32-element modules that comprised the whole detector array. Moreover, these elements suffered from severer pile-up distortion compared with other elements, resulting in further deviation from the Beer–Lambert law. Fig. 18(a)–(f) The spectral CT reconstructions with the first energy-bin [26, 33] keV. (a), (b) The reconstructed image of the NaI solutions, (c), (d) for solutions, and (e), (f) for solutions. The display window is [0, 1]. (g)–(k) The reconstructed images of the NaI solutions in all five energy-bins. The display window is . ![]() As shown in Fig. 18, the concentrations within the red circles were used to generate the validation dataset. Those within the green circles were used to generate the testing dataset. The rest were used for the training dataset. Because it was difficult to exactly obtain the effective spectrum of the x-ray beam used in experiment, we could not get the exact linear attenuation coefficients of every solutes. Thus, we used the maximum concentration solutions of three solutes and water as the basis materials to validate the decomposition algorithms. The regions within the centrifugal tubes were extracted to validate our CMD algorithm and compare it with EDI-MMD algorithm. To improve the generalization of the CNNs with limited training data, we randomly overturned the input patches.35 We obtained 36,324 samples for the training dataset, and 6054 samples for the validation and testing datasets, respectively. Because there always had some noise and artifacts, we used a smaller learning rate than that used in the simulation. We set , , P_STEP = 0.5, VAL_NUM = 6054, BATCH_NUM = 128, ITER_NUM = 64, , and . We also calculated the decomposition results by EDI-MMD algorithm as a contrast. The triplets set for EDI-MMD is in Table 6. Table 6Triplets set for phantom in experiment.
Figure 19 shows the decomposition results of six different concentration solutions highlighted in italics in Table 5. Figures 19(a)–19(f) are the results of NaI, NaI, , , , and , respectively. In each figures [Figs. 19(a)–19(f)], the columns from left to right are the results of four basis materials, NaI, Gd , , and , respectively. The rows from top to bottom are the ground truth, decomposition results of VGG–CMD, DRN-CMD, and EDI-MMD algorithms, respectively. We can find that both the VGG-CMD and DRN-CMD algorithms can provide good decomposition results. They are much better than the results of EDI-MMD algorithm which are seriously suffered from the ring artifacts existing in the spectral CT images. As shown in the last rows of Figs. 19(a)–19(d), the decomposition results are far from desire for the low concentration solutions. Moreover, the deeper network used in DRN-CMD algorithm can provide a smoother and more uniform results. In addition, compared with Ref. 37, the decomposition of traditional materials cannot achieve good results even if there is a clear prior knowledge (the solution consists of two basic materials: water and solvent) in the case of serious ring artifacts. Fig. 19Decomposition results of italic cells in Table 5, namely, (a) NaI, (b) NaI, (c) , (d) , (e) , and (f) . In each figure of (a–f), there are images. The rows from top to bottom are the ground truth, the decomposition results of VGG-CMD, DRN-CMD, and EDI-MMD algorithm. The columns from left to right are the results of NaI, , , and . ![]() Figure 20 shows the decomposition MSE of three algorithms. The average MSE of six testing images are for VGG–CMD, for DRN–CMD, and for EDI-MMD. In general, the MSE of CMD algorithms are reduced more than one order comparing to the MSE of EDI-MMD. Figure 21 shows the profiles of the decomposition results of NaI locating at the colorful lines in Fig. 19(a). MSE of the decomposition section line are 1.0645e-2, 2.7356e-3, and for three methods. The conclusion is similar to the numerical simulations. Although the CMD algorithm has never been trained by the testing solutions, the generalization of the CNNs guarantees its accuracy. Compared with the direct inversion decomposition method, CMD can resist the effects of ring artifacts that provide much more accurate and uniform decomposition results. 4.DiscussionCompared with traditional direct inversion algorithms, the proposed algorithm improves the decomposition results significantly especially for the cases that artifacts, e.g., ring and beam-hardening artifacts, existing in the CT images or the spectrum of x-ray source is distorted. For example, in current spectral CT systems using photon-counting detectors, the ring artifacts often appear due to the property inconsistency among the detector units (e.g., detection response function, dead time, pulse pile-up, etc.). As an image-domain method, no matter how the images are reconstructed, the proposed CMD algorithm can provide better decomposition results than direct inversion methods as shown in Fig. 19. It works well even the linear attenuation coefficients of the basis materials are close. Due to the robustness of CMD, it is also expected to be used in other incomplete-data or sparse-data CT cases in the future, e.g., limited-angle CT, few-view CT, and low-dose CT.46–48,49–52 Because the forward-propagation process in the CNN can be configured on GPU, the CMD algorithm can complete the material decompositions in a few seconds once it has been well trained. However, CMD also has some disadvantages. Traditional MMD algorithms in image domain only need the linear attenuation coefficients of the basis materials. When environment changes, we can update the algorithms easily. However, CMD method needs to gather training data and be retrained. The simulation 2 shows that the trained neural network is compatible with different imaging conditions, but it is almost impossible to generate every data for all the situations. For CMD, the biggest problem is that it is difficult to get the ground truth in practical situations for training. Even though we know all the material distributions of the designed phantoms, it is still difficult to ensure the diversity of data. We hope to train the neural network with simulation data as real as possible, and then decompose the materials directly on the real samples in our further experiments. Incorporating the models of the material decomposition and the physics of spectral CT, the direct iterative reconstruction and decomposition method sometimes provide good results. In the future, it is worth combining CMD algorithm with iterative reconstruction algorithm. 5.ConclusionIn this study, we proposed a deep-learning-based CMD algorithm to solve the MMD problem of spectral CT in image-domain. We redesigned two CNNs that worked well in the image classification field. Given the property of the Softmax layer, the output of the CMD can satisfy the mass and volume conservation constraints automatically. Simulation and experimental results proved the CMD algorithm can robustly solve the MMD problem of spectral CT. 6.Appendix6.1.Layers Used in This Work6.1.1.Convolutional layerThe convolutional layer is the most important layer in a CNN. The advanced features in an image are detected by convolutional layers. A convolutional layer receives a three-rank tensor and outputs a three-rank tensor as well. A convolutional layer can be divided into the weight tensor and the bias matrix . Receiving a tensor , the output tensor is Two types of convolutional operations are available for the edge of an image, i.e., using Eq. (11) to directly lead the size reduction or supplementing 0 around the input tensor to keep the size unchanged. 6.1.2.Activation layerAs shown above, no matter how many convolutional layers are piled up, the effect is the same as the case of only one convolutional layer because the operation in the convolutional layer is linear.38 To enhance the capability of the CNN, nonlinearity should be introduced. We explain two types of activation layers, each of which has its own features. Rectified linear unit (ReLU) is a simple activation function that effectively introduces nonlinearity into CNNs and has thus been widely used since its discovery. ReLU sets the negative value of the output tensor of a convolutional layer to 0, that is, ReLU simplifies the process of derivation and lightens the vanishing gradient problem of traditional artificial neural networks, resulting in the rapid convergence of networks. Softmax is the last layer of a CNN. It receives a one-rank tensor and outputs a one-rank tensor:38 Obviously, the output of the Softmax activation layer satisfies the constraints of Eq. (2). In image classification, the output of Softmax (in common, the Softmax activation layer is the last layer of the CNN, so it is also the output of a CNN) is considered as the possibility of each image type predicted by CNN. However, in the material decomposition problem, we train the output as the volume fraction of each material of the central pixel of the input patch. 6.1.3.Batch-normalization layerSimilar to the ordered subset used in CT reconstruction, which can accelerate convergence, we send a certain number of samples to update the CNN at a time during training. This set of samples is called a minibatch. BN is set between the convolutional and activation layers to scale and translate the values in a batch uniformly and to improve the following activation layers.38,53 BN accelerates convergence, relaxes the requirements of the initial values of networks, and makes very deep networks possible. 6.1.4.Max pooling layerThe pooling layer is a screening operation of the input three-rank tensor. Given the strides of the and directions and , the output is As shown in Eq. (14), the max pooling layer shrinks the size of the image and introduces nonlinearity into the CNN. The max pooling layer also preserves only one value of a region of the image. Therefore, small translations in the image will be ignored by the CNN. In another word, the translation invariance is improved. 6.1.5.Fully connected layerAn inchoate neural network consists of fully connected layers and activation layers. This composition is ineffective.38 After extracting features through the convolutional layer, the fully connected layers can work effectively. The fully connected layer receives a 1-rank tensor and outputs a one-rank tensor: In a CNN, the first two parameters of the image’s size, and , will shrink when flowing through the convolutional and pooling layers. The third parameter enlarges when operated by the convolutional layer. Therefore, the image becomes increasingly small while the number of its channels increases. Finally, we ignore the difference in the first two parameters, flatten the three-rank tensor image into a one-rank tensor, and send it into the fully connected layer. 6.2.Codes AvailableOur trained codes, some synthetic and real datasets used in this paper can be downloaded from Github repository available at: https://github.com/ZhengyangChen/Convolutional_material_decomposition. AcknowledgmentsThis work was supported in part by NSFC 11775124, 61571256, 11525521, 81427803, and the National Key Research and Development Program of China, 2017YFC0109103 and 2018YFC0115502. ReferencesY. Long and J. A. Fessler,
“Multi-material decomposition using statistical image reconstruction for spectral CT,”
IEEE Trans. Med. Imaging, 33
(8), 1614
–1626
(2014). https://doi.org/10.1109/TMI.2014.2320284 ITMID4 0278-0062 Google Scholar
M. Patino et al.,
“Material separation using dual-energy CT: current and emerging applications,”
Radiographics, 36
(4), 1087
–1105
(2016). https://doi.org/10.1148/rg.2016150220 Google Scholar
J. Chu et al.,
“Combination of current-integrating/photon-counting detector modules for spectral CT,”
Phys. Med. Biol., 58
(19), 7009
–7024
(2013). https://doi.org/10.1088/0031-9155/58/19/7009 PHMBA7 0031-9155 Google Scholar
P. R. Mendonça, P. Lamb and D. V. Sahani,
“A flexible method for multi-material decomposition of dual-energy CT images,”
IEEE Trans. Med. Imaging, 33
(1), 99
–116
(2014). https://doi.org/10.1109/TMI.2013.2281719 ITMID4 0278-0062 Google Scholar
L. Li et al.,
“Spectral CT modeling and reconstruction with hybrid detectors in dynamic-threshold-based counting and integrating modes,”
IEEE Trans. Med. Imaging, 34
(3), 716
–728
(2015). https://doi.org/10.1109/TMI.2014.2359241 ITMID4 0278-0062 Google Scholar
Y. Xue et al.,
“Statistical image-domain multimaterial decomposition for dual-energy CT,”
Med. Phys., 44
(3), 886
–901
(2017). https://doi.org/10.1002/mp.12096 MPHYA6 0094-2405 Google Scholar
R. E. Alvarez and A. Macovski,
“Energy-selective reconstructions in x-ray computerised tomography,”
Phys. Med. Biol., 21
(5), 733
–744
(1976). https://doi.org/10.1088/0031-9155/21/5/002 PHMBA7 0031-9155 Google Scholar
L. Lehmann et al.,
“Generalized image combinations in dual KVP digital radiography,”
Med. Phys., 8
(5), 659
–667
(1981). https://doi.org/10.1118/1.595025 MPHYA6 0094-2405 Google Scholar
H. Xue et al.,
“A correction method for dual energy liquid CT image reconstruction with metallic containers,”
J. X-Ray Sci. Technol., 20
(3), 301
–316
(2012). https://doi.org/10.3233/XST-2012-0339 JXSTE5 0895-3996 Google Scholar
L. Li et al.,
“A dynamic material discrimination algorithm for dual MV energy x-ray digital radiography,”
Appl. Radiat. Isot., 114 188
–195
(2016). https://doi.org/10.1016/j.apradiso.2016.05.018 ARISEF 0969-8043 Google Scholar
P. Sukovic and N. H. Clinthorne,
“Penalized weighted least-squares image reconstruction for dual energy x-ray transmission tomography,”
IEEE Trans. Med. Imaging, 19
(11), 1075
–1081
(2000). https://doi.org/10.1109/42.896783 ITMID4 0278-0062 Google Scholar
K. Taguchi et al.,
“Image-domain material decomposition using photon-counting CT,”
Proc. SPIE, 6510 651008
(2007). https://doi.org/10.1117/12.713508 PSISDG 0277-786X Google Scholar
Y. Zou and M. D. Silver,
“Analysis of fast kV-switching in dual energy CT using a pre-reconstruction decomposition technique,”
Proc. SPIE, 6913 691313
(2008). https://doi.org/10.1117/12.772826 PSISDG 0277-786X Google Scholar
C. Maaß, M. Baer and M. Kachelrieß,
“Image-based dual energy CT using optimized precorrection functions: a practical new approach of material decomposition in image domain,”
Med. Phys., 36
(8), 3818
–3829
(2009). https://doi.org/10.1118/1.3157235 MPHYA6 0094-2405 Google Scholar
T. Niu et al.,
“Iterative image-domain decomposition for dual-energy CT,”
Med. Phys., 41
(4), 041901
(2014). https://doi.org/10.1118/1.4866386 MPHYA6 0094-2405 Google Scholar
M. Petrongolo, X. Dong and L. Zhu,
“A general framework of noise suppression in material decomposition for dual-energy CT,”
Med. Phys., 42
(8), 4848
–4862
(2015). https://doi.org/10.1118/1.4926780 MPHYA6 0094-2405 Google Scholar
R. F. Barber et al.,
“An algorithm for constrained one-step inversion of spectral CT data,”
Phys. Med. Biol., 61
(10), 3784
–3818
(2016). https://doi.org/10.1088/0031-9155/61/10/3784 PHMBA7 0031-9155 Google Scholar
Y. Zhao, X. Zhao and P. Zhang,
“An extended algebraic reconstruction technique (E-ART) for dual spectral CT,”
IEEE Trans. Med. Imaging, 34
(3), 761
–768
(2015). https://doi.org/10.1109/TMI.2014.2373396 ITMID4 0278-0062 Google Scholar
T. Zhao, L. Li and Z. Chen,
“K-edge eliminated material decomposition method for dual-energy x-ray CT,”
Appl. Radiat. Isot., 127 231
–236
(2017). https://doi.org/10.1016/j.apradiso.2017.06.018 ARISEF 0969-8043 Google Scholar
J. W. Stayman and S. Tilley,
“Model-based multi-material decomposition using spatial-spectral CT filters,”
(2018) http://aiai.jhu.edu/papers/CT2018_stayman.pdf Google Scholar
A. M. Alessio and L. R. MacDonald,
“Quantitative material characterization from multi-energy photon counting CT,”
Med. Phys., 40
(3), 031108
(2013). https://doi.org/10.1118/1.4790692 MPHYA6 0094-2405 Google Scholar
P. R. Mendonça et al.,
“Multi-material decomposition of spectral CT images,”
Proc. SPIE, 7622 76221W
(2010). https://doi.org/10.1117/12.844531 PSISDG 0277-786X Google Scholar
G. Wang,
“A perspective on deep imaging,”
IEEE Access, 4 8914
–8924
(2016). https://doi.org/10.1109/ACCESS.2016.2624938 Google Scholar
Z. Chen and L. Li,
“Preliminary research on multi-material decomposition of spectral CT using deep learning,”
in Fully3D,
523
–526
(2017). Google Scholar
K. Simonyan and A. Zisserman,
“Very deep convolutional networks for large-scale image recognition,”
(2014). Google Scholar
J. Schmidhuber,
“Deep learning in neural networks: an overview,”
Neural Networks, 61 85
–117
(2015). https://doi.org/10.1016/j.neunet.2014.09.003 NNETEB 0893-6080 Google Scholar
K. He et al.,
“Deep residual learning for image recognition,”
in Computer Vision and Pattern Recognition,
770
–778
(2016). Google Scholar
H. Greenspan, B. van Ginneken and R. M. Summers,
“Guest editorial deep learning in medical imaging: overview and future promise of an exciting new technique,”
IEEE Trans. Med. Imaging, 35
(5), 1153
–1159
(2016). https://doi.org/10.1109/TMI.2016.2553401 ITMID4 0278-0062 Google Scholar
P. Hu et al.,
“Automatic 3D liver segmentation based on deep learning and globally optimized surface evolution,”
Phys. Med. Biol., 61
(24), 8676
–8698
(2016). https://doi.org/10.1088/1361-6560/61/24/8676 PHMBA7 0031-9155 Google Scholar
P. Moeskops et al.,
“Automatic segmentation of MR brain images with a convolutional neural network,”
IEEE Trans. Med. Imaging, 35
(5), 1252
–1261
(2016). https://doi.org/10.1109/TMI.2016.2548501 ITMID4 0278-0062 Google Scholar
M. Anthimopoulos et al.,
“Lung pattern classification for interstitial lung diseases using a deep convolutional neural network,”
IEEE Trans. Med. Imaging, 35
(5), 1207
–1216
(2016). https://doi.org/10.1109/TMI.2016.2535865 ITMID4 0278-0062 Google Scholar
H.-C. Shin et al.,
“Deep convolutional neural networks for computer-aided detection: CNN architectures, dataset characteristics and transfer learning,”
IEEE Trans. Med. Imaging, 35
(5), 1285
–1298
(2016). https://doi.org/10.1109/TMI.2016.2528162 ITMID4 0278-0062 Google Scholar
Q. Dou et al.,
“Automatic detection of cerebral microbleeds from MR images via 3D convolutional neural networks,”
IEEE Trans. Med. Imaging, 35
(5), 1182
–1195
(2016). https://doi.org/10.1109/TMI.2016.2528129 ITMID4 0278-0062 Google Scholar
J. M. Wolterink et al.,
“Generative adversarial networks for noise reduction in low-dose CT,”
IEEE Trans. Med. Imaging, 36
(12), 2536
–2545
(2017). https://doi.org/10.1109/TMI.2017.2708987 ITMID4 0278-0062 Google Scholar
H. Chen et al.,
“Low-dose CT via convolutional neural network,”
Biomed. Opt. Express, 8
(2), 679
–694
(2017). https://doi.org/10.1364/BOE.8.000679 BOEICL 2156-7085 Google Scholar
Y. Han, J. Kang and J. C. Ye,
“Deep learning reconstruction for 9-view dual energy CT baggage scanner,”
(2018). Google Scholar
S. Dorn et al.,
“Prior-based multi material decomposition for duel energy CT,”
(2018) https://www.dkfz.de/en/roentgenbildgebung/ct/ct_conference_contributions/Prior-Based-Multi-Material-Decomposition-for-DECT_SabrinaDorn.pdf?m=1527056293 Google Scholar
I. Goodfellow et al., Deep Learning, MIT Press, Cambridge
(2016). Google Scholar
L. Li et al.,
“A tensor PRISM algorithm for multi-energy CT reconstruction and comparative studies,”
J. X-Ray Sci. Technol., 22
(2), 147
–163
(2014). https://doi.org/10.3233/XST-140416 JXSTE5 0895-3996 Google Scholar
J. Zhang, I. Mitliagkas and C. Ré,
“Yellowfin and the art of momentum tuning,”
(2017). Google Scholar
A. Choromanska et al.,
“The loss surfaces of multilayer networks,”
in AISTATS,
192
–204
(2015). Google Scholar
S. Hochreiter and J. Schmidhuber,
“Flat minima,”
Neural Comput., 9
(1), 1
–42
(1997). https://doi.org/10.1162/neco.1997.9.1.1 NEUCEB 0899-7667 Google Scholar
D. Soudry and E. Hoffer,
“Exponentially vanishing sub-optimal local minima in multilayer neural networks,”
(2017). Google Scholar
R. Li, L. Li and Z. Chen,
“Spectrum reconstruction method based on the detector response model calibrated by x-ray fluorescence,”
Phys. Med. Biol., 62
(3), 1032
–1045
(2017). https://doi.org/10.1088/1361-6560/62/3/1032 PHMBA7 0031-9155 Google Scholar
E. Y. Sidky and X. Pan,
“Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,”
Phys. Med. Biol., 53
(17), 4777
–4807
(2008). https://doi.org/10.1088/0031-9155/53/17/021 PHMBA7 0031-9155 Google Scholar
L. Li et al.,
“Full-field fan-beam x-ray fluorescence computed tomography with a conventional x-ray tube and photon-counting detectors for fast nanoparticle bioimaging,”
Opt. Eng., 56
(4), 043106
(2017). https://doi.org/10.1117/1.OE.56.4.043106 Google Scholar
H. Q. Le and S. Molloi,
“Least squares parameter estimation methods for material decomposition with energy discriminating detectors,”
Med. Phys., 38
(1), 245
–255
(2011). https://doi.org/10.1118/1.3525840 MPHYA6 0094-2405 Google Scholar
Z. Chen et al.,
“A limited-angle CT reconstruction method based on anisotropic TV minimization,”
Phys. Med. Biol., 58
(7), 2119
–2141
(2013). https://doi.org/10.1088/0031-9155/58/7/2119 PHMBA7 0031-9155 Google Scholar
M. Chang et al.,
“A few-view reweighted sparsity hunting (FRESH) method for CT image reconstruction,”
J. X-Ray Sci. Technol., 21
(2), 161
–176
(2013). https://doi.org/10.3233/XST-130370 JXSTE5 0895-3996 Google Scholar
J. Wang et al.,
“Penalized weighted least-squares approach to sinogram noise reduction and image reconstruction for low-dose x-ray computed tomography,”
IEEE Trans. Med. Imaging, 25
(10), 1272
–1283
(2006). https://doi.org/10.1109/TMI.2006.882141 ITMID4 0278-0062 Google Scholar
S. Ioffe and C. Szegedy,
“Batch normalization: accelerating deep network training by reducing internal covariate shift,”
(2015). Google Scholar
|