CT IMAGE RECONSTRUCTION ON A LOW DIMENSIONAL MANIFOLD

. The patch manifold of a natural image has a low dimensional structure and accommodates rich structural information. Inspired by the recent work of the low-dimensional manifold model (LDMM), we apply the LDMM for regularizing X-ray computed tomography (CT) image reconstruction. This proposed method recovers detailed structural information of images, signiﬁ-cantly enhancing spatial and contrast resolution of CT images. Both numerically simulated data and clinically experimental data are used to evaluate the proposed method. The comparative studies are also performed over the simultaneous algebraic reconstruction technique (SART) incorporated the total variation (TV) regularization to demonstrate the merits of the proposed method. Results indicate that the LDMM-based method enables a more accurate image reconstruction with high ﬁdelity and contrast resolution.


1.
Introduction. X-ray computed tomography (CT) is a major imaging modality in medical, security, and industrial applications. The filtered back-projection (FBP) is an efficient and robust method for x-ray CT image reconstruction [10], but it generates strong noise and artifacts in the cases of low-dose or incomplete datasets. Extensive efforts have been made to improve image quality for practical purposes [9,19,6]. Iterative methods incorporate prior information of images, and offer distinct advantages over the analytic methods in cases of noisy and few-view data. The statistical iterative methods model the statistics of photons to improve the reconstructed image quality from the low-dose acquisitions [6,20]. Recently, the compressive sensing (CS) approach [3,4] is applied for the image reconstruction from less measurements than that required by the Nyquist-Shannon sampling theorem. Based on the CS theory, image reconstruction algorithms were developed for various problems of CT image reconstruction to improve image quality and reduce radiation dose, such as total variation (TV) regularization [19,20], prior image constrained compressed sensing (PICCS) [5], nonlocal mean (NLM) [1,9], and dictionary learning (DL) [21]. TV is a typical sparse transform for an image, and is a popular regularization form for image reconstruction due to its ability to preserve image edges. However, it is effective only for reconstruction of piecewise constant images and would over-smoothen textured regions, which may sacrifice important details. PICSS could be seen as a generalization of the TV regularization method. It incorporates a preliminary reconstructed image into CS reconstruction technique to achieve more accurate image reconstruction [5]. NLM exploits a high degree of redundancy of an image for de-noising [1]. The similarity is derived from intensity differences between neighboring patches of pixels or voxels. A non-linear filter can be used to reduce image noise by updating each pixel value with a weighted average of its neighbors according to the similarity of involved patches. DL builds adaptive sparse representation elements from a training set of images, and utilizes domain knowledge at a deeper level [21]. The dictionary tends to capture local image features effectively and helps image denoising and feature inference. However, the structural differences between a true image and training images could affect the image reconstruction quality.
The idea of the proposed X-ray CT image reconstruction method is inspired by a recent method called the low-dimensional manifold model (LDMM) [14,16]. Using the image patches discussed in nonlocal methods [2], the LDMM interprets image patches as a point cloud sampled in a low-dimensional manifold embedded in a high dimensional ambient space, which provides a new way of regularization by minimizing the dimension of the corresponding image patch manifold. This can be explained as a natural extension of the idea of low-rank regularization for linear objects to data with more complicated structures [7]. Moreover, the authors in [14] elegantly find that the point-wisely defined manifold dimension can be computed as a Dirichlet energy of the coordinate functions on the manifold, whose corresponding boundary value problem can be further solved by a point integral method proposed in [14]. The LDMM performs very well in image imprinting and super-resolution. The patch manifold of images is generally a low dimensional structure, and accommodates rich structural information [16]. With the LDMM prior knowledge on images, this proposed method significantly enhances contrast resolution of image reconstruction.
In this paper, LDMM-based regularization method is proposed for CT image reconstruction. The LDMM-based reconstruction method is to maximize the data fidelity and minimize the manifold dimensionality, which is performed using the Bregman iteration [15] by updating the patch manifold structure iteratively. By a standard variational approach, the regularization model of image reconstruction can be reduced to the Laplace-Beltrami equations over a point cloud, which is solved using the point integral method [11].
The rest of the paper is organized as follows. In section 2, we provide a detailed description for the CT image reconstruction based on LDMM. A numerical reconstruction algorithm is presented based on Bregman iteration. In section 3, we perform the image reconstruction using numerical simulation data and the clinical raw projection data to evaluate the proposed LDMM-based image reconstruction method. In addition, we also conduct comparative studies with popular simultaneous algebraic reconstruction technique (SART) with TV method to present the merits of the proposed method. After that, we conclude the paper in the last section.
2. Image reconstruction method. In this section, we first review the statistical model of x-ray CT imaging. After that, we discuss the proposed method of CT image reconstruction based on LDMM and its numerical algorithm.
2.1. Statistical model for X-ray CT imaging. In x-ray CT imaging, the number ξ of x-ray photons recorded by a detector element is a random variable, which obeys a Poisson distribution [10]: The expectation value of x-ray photons along a path from x-ray source to i-th detector element obeys Beer-Lambert law: where b i is the number of x-ray photons detected by i-th detector element in the blank scanning (without any object in the beam path), and µ( r) is the linear attenuation coefficient of the object. To implement the numerical computation, Eq.
(2) can be discretized as, where µ is a vector composed of pixel values on image of linear attenuation coefficients, and A i is the weighting coefficients of the pixel values on i-th beam path. Since data are independent between detectors, the likelihood function for x-ray photons probability distribution on detectors is defined by, where Y = (y 1 , y 2 , · · · , y N ) T . According to the Bayesian rule: p(µ|Y )p(Y ) = (Y |µ)p(µ), the image reconstruction task can be implemented by maximizing a posteriori (MAP) distribution p(µ|Y ), which is equivalent to the following minimization problem in term of the monotonic property of the natural logarithm [20,12]: where r(µ) = − ln(p(µ)) is a regularization term expressing the prior knowledge about the attenuation image µ , and N is the total number of x-ray beam paths.
In the context, we propose to use the low-dimension of an image as prior knowledge to conduct the image reconstruction, which is discussed in the next section. After inserting Eq. (3) in Eq. (5), a second-order approximation is applied to simplify the complicated optimization to a quadratic optimization: Image reconstruction based on LDMM. The classical image restoration models mainly focus on local properties of the objective image, such as smoothness and jumps. Image features can be further enhanced due to its possible repetitive patterns non-locally. The nonlocal based image restoration methods [2] extract and match non-local repetitive structures of images using image patches. An essential observation of nonlocal methods is that images can be restored by enhancing similar patterns which may not lie in nearby regions of the original image domain. Therefore, comparing with the direct regularization methods on the image domain, the quality of image restoration can be usually improved using nonlocal methods. For instance, nonlocal based variational methods [2,8,22] and nonlocal based wavelet frame based methods demonstrate outstanding image restoration results [17]. Let µ denote an image containing m × n pixels: µ = {µ(i, j) | 1 ≤ i ≤ m, 1 ≤ j ≤ n} , and P (i0,j0) (I) denotes an image patch centering at (i 0 , j 0 ) with size of 2s 1 s 2 , namely, An image is decomposed into a set of patches. These patches can be overlapping or nonoverlapping. Let P(µ) denote the patch transformation which maps any given image I to be the set of image patches. P(µ) can be also viewed as a point set in R d with a dimension of d = 2s 1 × s 2 . It samples a low dimensional manifold M(µ) embedded in R d , which is called the patch manifold of µ [14,16]. The patch manifold is low dimensional for many natural images, such as X-ray CT images. A patch manifold can be constructed directly from an image using patching size 2s 1 × 2s 2 . As an example illustrated in the right image of Fig. 1, we construct a patch manifold of a CT image using patching size 16 × 16 and color-code its piecewise dimension on the image to show the variation of the manifold dimension. More recently, [14] proposes to regularize the dimension of the patch manifold M(µ) for image restoration. Inspired by [14], we use the dimension of the patch manifold as a regularization term to seek the dimension of its patch manifold as small as possible such that detail structure information can be magnified for the CT image reconstruction. Therefore, the optimization model Eq. (6) is reformulated for the measurement data fidelity and the manifold dimensional quantification: where dim(M(µ))(x) denotes the dimension of the patch manifold M(µ) of an image µ at x, a point in M(µ) ⊂ R d . The patch manifold of a nature image may be a set of several manifolds with different dimensions, corresponding to different patterns of the image. In [14], the authors elegantly demonstrate that the dimension of the patch manifold for an image µ at given patch x can be calculated by the coordinate function, where α i is the embedding coordinate function defined by α i (x) = x i for any x = (x 1 , x 2 , · · · , x d ) ∈ M ⊂ R d . Combining Eqs. (7) and (8), we obtain (9) arg min where M is a manifold embedded in R d , and P(µ) is the patch set of image µ. Note that the dimension formula (8) holds point-wisely for the embedding manifolds. The above regularization is actually the L 1 norm of the local dimension. The optimization (9) can be solved by alternating direction method similar as the one proposed in [14]. In other words, given an estimation of the manifold M n and an estimation of image µ n satisfying P(µ n ) ⊂ M n , the coordinate functions α = {α i } and an updated image µ are computed through the following optimization.
where M = M n , u represents any α i , n is the out normal of M with the boundary ∂M and Ω = P(µ n ). Recently, the point integral method has been proposed to solve Laplace-Beltrami equation over a point cloud [3]. The main idea of the point integral method is to apply following integral approximation for the differential term in Laplace-Beltrami equation: where R t (x, y) are kernel functions given as follows, where C t is a normalizing factor. Using the integral approximation (13), following integral equation can be obtained to approximate the Laplace-Beltrami equation, The integral equation (15) can be further discretized into a matrix equation over the point set using quadrature rule [14]: ,··· ,N is the weight matrix, D = diag({ j w ij } 1,··· ,N ) and V = P(µ k ) − Q k . Thus, the optimization (11a) can be solved based on the matrix equation (16). The detailed formulation and alternating minimization steps for solving Eq. (10) are described in the flowchart for Algorithm 1.

Algorithm 1
Initialize an initial image µ, Q 0 = 0 and set parameters λ and β. while the current solution is not converged do 1.Compute the weight matrix W = {w ij } from the patch image P (µ k ) and let weight w ij = R t (x i , x j ), and the matrices 2.Solve the linear systems: Update µ by solving the problem: 4.Update: Q k+1 = Q k + U − P (µ k+1 ) end while 3. Numerical experiments. In the section, we test the proposed LDMM-based reconstruction method with numerical simulation and real experimental datasets obtained on a GE clinical CT scanner. We constructed a patch manifold of an image with a patch size of 1616 by shifting 8 pixels. In addition, we also perform the comparative studies with the state of the art method, the simultaneous algebraic reconstruction technique (SART) with a total variation (TV) [18], to demonstrate the merits of the proposed method. All numerical computations in this section are implemented by MATLAB in a PC with 16G RAM and 2.8GHz CPU.
3.1. Simulation data. A phantom was adapted from a CT slice to evaluate the proposed algorithms. Data acquisitions are simulated with polychromatic x-ray source operated at 120 kVp/5mAs dose for the x-ray imaging. The radius of the scanning trajectory was 53.852cm. The source-to-detector distance was 94.6746 cm. 540 projections are uniformly acquired over a 360-degree angular range. 765 detector elements are equiangular distributed on each projection view. The phantom is discretized into a 512 × 512 matrix. We choose the regularization parameters as λ = 0.3 and β = 0.15. The x-ray imaging process was simulated to generate projection data according to the x-ray propagation forward model, the polychromatic Beer-Lambert law. The projection datasets were corrupted by Poisson noise to simulate real x-ray imaging experiments. We performed the image reconstruction using the proposed LDMM-based image reconstruction and the simultaneous algebraic reconstruction technique (SART) with a total variation (TV), respectively. Results show that the LDMM-based method is able to produce more accurate image reconstruction with high fidelity and detailed features than the SART+TV method that would over-smoothen textural pattern in reconstructed image. We calculated the peak-to-noise ratio (PSNR) and structural similarity (SSIM) for the reconstructed image, and obtained PSNR of 52.42 and SSIM of 0.9940 for the LDMM-based reconstruction, and PSNR of 39.53 and SSIM of 0.9911 for the SART+TV method. Fig. 2 shows a comparison between the reconstructed images and the ground truth. Fig. 3 presents the comparison of profiles along the horizotal and vertical midlines in the phantom and reconstructed images.  3.2. Low dose data. A realistic phantom adapted from a human CT slice is used to evaluate the proposed algorithms. We use an computer-assisted tomography simulation (Cat-Sim) software [13], which was developed by GE Global Research Center, to simulate x-ray imaging for the phantom. CatSim incorporates polychromaticity, realistic quantum and electronic noise models, finite focal spot size and shape, finite detector cell size, and detector cross-talk for the simulation of real x-ray imaging. All acquisitions are simulated with polychromatic x-ray source operated at 120 kVp and 0.2mSv dose for the low dose x-ray imaging. The radius of the scanning trajectory is 54.1cm. Source-todetector distance is 94.9cm. 888 detector elements are equiangular distributed on each projection view. 984 projections are uniformly acquired over a 360-degree angular range, generating the sinogram, as shown in Fig. 6. The phantom is discretized into a 512 × 512 matrix. We choose the regularization parameters as λ = 0.3 and β = 0.3. We performed the image reconstruction respectively using the proposed LDMM-based image reconstruction and SART +TV, which achieve PSNR of 18.72 and SSIM of 0.79 for the LDMM-based reconstruction method, and PSNR of 18.69 and SSIM of 0.74 for the SART+TV method. Results show that the LDMM-based image reconstruction method had a better contrast resolution of the reconstructed image than the SART+TV method, as shown in Fig. 5.  . The sinogram measured from a clinical x-ray CT scanner.

Clinical data.
A clinical CT raw projection dataset obtained from GE Healthcare was used to evaluate the LDMM-based image reconstruction method. After appropriate preprocessing for the projection data, we obtained a set of fanbeam sinogram, as shown in Fig. 6. In the x-ray imaging, the field of view (FOV) is of a 25 cm radius, and the radius of the scanning trajectory is 54.1cm. Source-todetector distance is 94.9cm. 984 projections are uniformly acquired over a 360degree angular range. 888 detector elements with 1.024mm pitch were equiangular distributed on a projection view. The image matrix was of 512 × 512 pixels. We choose the regularization parameters as λ = 0.2 and β = 0.5. We conducted the image reconstruction from the sinogram using the proposed LDMM-based reconstruction method. For comparison, the FBP method and the SART with TV regularization were applied as well to perform the image reconstruction from same projection dataset. Results show that the LDMM-based image reconstruction outperforms the other two reconstruction methods, as shown in Fig. 7. The LDMMbased method well preserves structural information especially texture features of the reconstructed image. The SART with TV iteration method was suitable to reconstruct high contrast images, whereas it would over-smoothen textured regions in medical images, resulting in the loss of details. FBP keeps the structural information, but it produced noisy image. 4. Discussions and conclusion. The major contribution in this paper is to propose an image reconstruction method aided by the regularization of a low dimensional manifold model (LDMM). The patch manifold of a natural image is generally with a low dimensional structure, and accommodates rich structural information. LDMM regularization method well recovers structural information of images, and promises substantially to increase spatial and contrast resolution of the image reconstruction. The comparison between the proposed method and the representative the SART with TV methods has been performed to illustrate the merits of the LDMM-based reconstruction approach. The raw datasets from a clinical CT scanner have been used to evaluate the performance of the image reconstruction methods. Results show that the regularization method of low dimensional manifold is an efficient and robust image reconstruction technique, and well preserves image edges and structural details of the reconstructed image comparing to the state of the art SART with TV image reconstruction. The iterative algorithm also incorporates prior knowledge, and account for photon statistics at a low dose level. Major computational cost is matrix-vector multiplication operations in the LDMM-based image reconstruction. Because matrix-vector multiplication is well computed parallel, the computational speed of the proposed iterative method can be significantly improved by parallel programming. This LDMM-based image reconstruction approach is very promising for medical imaging and other applications.