Edge-oriented dual-dictionary guided enrichment (EDGE) for MRI-CT image reconstruction

2.1 Background

CT and MRI have different imaging principles. The CT image reconstruction is to solve the inverse problem of the integrals of the X-ray attenuation coefficient distribution in the human body. X-ray projections or integrals of linear attenuation coefficients is usually described as ∫ μ ( l ) dl = ln ( I 0 I )(1) where I ₀ is the incident intensity of an X-ray beam, I is the transmitted intensity, and μ (l) is the X-ray attenuation coefficient along an X-ray path l. During a CT scan, a number of X-ray projections are acquired. Then, the relationship of Equation (1) can be discretized in the matrix form: Ax = p(2) where the vector x ∈ ℝ N × 1 stands for a CT image, N is the total number of pixels, A ∈ ℝ M × N is a system matrix which is determined by the geometry of the CT system, and the vector p ∈ ℝ M × 1 stands for projection data with M being the product of the number of angle views and the number of detector elements. The CT reconstruction problem is to solve the inverse problem of the linear systemEquation (2).

Different from X-ray CT, each pixel value of an MRI image represents a property relevant to the resonance nuclear density (hydrogen protons) and two relaxation times, with the details depending on a specific pulse sequence. The scan signals, known as k-space data, are recorded by phase and frequency encoding. An image can be reconstructed by the inverse Fourier transform, x = F - 1 y(3) where x is an MRI image, and y is its k-space data.

Although Equations (2-3) represent different models, the two reconstruction problems can be both transfered into similar mixed-norm regularization problems according to compressed sensing (CS) theory [16]. Moreover, though the physical principles and quantities are different for these imaging modalities, the human structural information in these images are correlated [17, 18]. Inspired by the correlative nature of CT and MRI datasets, we are motivated to study whether an MRI image can be reconstructed from less k-space data when the prior CT information of the same patient is available.

Recently, CS theory has served as a powerful tool in producing high quality image restoration and reconstruction from under-sampled data which may be much less than the requirements of the traditional Shannon/Nyquist sampling theory [19 –24]. According to the CS theory, a highly under-sampled image reconstruction problem is to solve an underdetermined system of linear equations F _u x = y by minimizing the l ₀ norm of a sparsified transformΨx. The corresponding optimization problem may be described as min x ∥ Ψ x ∥ 0 s . t . F u x = y(4)

This problem is also known as a sparse coding problem which is an NP-hard problem. Usually, the l ₀ norm can be replaced by the l ₁ norm because the l ₁ norm minimization is a proxy for the l ₀ sparsity constraint when the restricted isometry property (RIP) holds [19]. Then, the problem can be efficiently solved by linear programming in the real domain or second order cone programming in the complex domain.

DL is a promising implementation of CS theory by training a dictionary for efficient sparse representation [25, 26]. It aims to find a proper representation of data in adaptive low-dimensional subspaces. DL related methods are widely used, already giving significant improvements in the field of medical image reconstruction [27 –29].

Given an image of size N × N, it can be decomposed into small patches of size b × b, b << N. Each patch can be expressed as a vector x p ∈ ℝ b × 1. All the patches are extracted from the image according to a designed partition. A dictionary D ∈ ℝ b × K is a matrix consisting of k atoms (columns). An initial dictionary generated from sample images is usually redundant or over-complete. With this dictionary D, each vector x _p of an image can be sparsely represented as ∥ x p - D α ∥ 2 2 < ε(5) where the error bound ε > 0. α ∈ ℝ K × 1 is a sparse representation vector which has few none-zero elements, i.e., ∥α ∥ ₀ << N. Finding a sparse representation of the vector x _p is to solve the following l ₀ norm minimization problem min α ∥ α ∥ 0 s . t . ∥ x p - D α ∥ 2 2 < ε(6)

If an image contains S patches, DL is to find a dictionary that all the patches can be sparsely represented by min D , α ∑ s = 1 S ( ∥ x p _ s - D α s ∥ 2 2 ) s . t . ∥ α s ∥ 0 < T 0(7) where T ₀ is a required sparsity level.

Many algorithms were proposed to solve the dictionary learning problem Equation (7) [25 , 30]. As one of the popular algorithms, the k-means singular value decomposition (K-SVD) algorithm iteratively updates the dictionary atoms to better fit data. It does SVD on the errors, and refines the current dictionary atoms and coefficients simultaneously. As a widely-used method, K-SVD enjoys an excellent performance in convergence and sparsity [25]. The DL method was reported to perform well for MRI reconstruction from under-sampled k-space data.

However, image features may be in large disparities across multimodality images or even within one modality, for example, dual-energy CT. These different images cannot be directly integrated into one dictionary. Thus, two or more dictionaries are required to describe image features in different domains. A number of dual-dictionary learning (DDL) methods were proposed to solve this problem, such as MRI-CT reconstruction, multi-energy CT reconstruction, image reconstruction from sparse data [12 , 31–33]. DDL enables us to establish a connection between different types of images using the prior information of two modalities.

In multimodality image reconstruction, the existing DDL algorithms usually use the Euclidean distances among different patches to choose the dictionary atoms, which work well for the multimodality images in consistent contrasts, or from a high contrast image to a low contrast destination; for example, estimating or reconstructing a CT image of the human brain from an MRI image using DDL [12]. However, these algorithms would perform poorly in reconstructing a high contrast image from a low contrast counterpart; for example, reconstructing an MRI image of the human brain from a CT image. Therefore, we are motivated to investigate an edge-oriented dual-dictionary guided enrichment (EDGE) alternative to solve the MRI reconstruction problem from the training dataset consisting of well-registered MRI and CT images subject to highly under-sampled k-space data and a complete CT dataset. The key is to use the relative positions guided by the edges and contours in the images to establish the connections among the patches instead of the Euclidean distance that can be misled by inconsistent contrast and grey-scalebias.

The EDGE algorithm has two steps. First, an initial MRI image is generated by an edge-oriented DDL algorithm from a training dataset consisting of well-registered MRI and CT images of selected previous patients, and the CT image of a current patient. The EDGE mechanism enables the two dictionaries well mapped and self-adaptive. Then, the estimated MRI image will be used as an initial value to facilitate the final MRI image reconstruction from highly under-sampled k-space data.

Step 1. Initial MRI image estimation using the EDGE strategy.

First, to establish two dictionaries D ^MR and D ^CT, we need to collect a series of MRI and CT images of patients. As a training dataset, these images may be gathered from different CT and MRI scanners. The dual-modality images of each patient require to be well-registered by structure similarity. For ahigh-resolution CT image u CT ∈ ℝ N × N, a patch vector u i CT ∈ ℝ b × 1 (b is the patch size) can be represented by a specific dictionary D CT = [ d 1 CT , d 2 CT , . . . , d K CT ] ∈ ℝ b × K according to DL theory as [27] ∥ u i CT - D CT α i CT ∥ 2 2 < ε(8)

Using the orthogonal matching pursuit (OMP) algorithm, the patches in Equation (8) can be represented as [34], u i CT = D CT α i CT + η i(9)

When MRI and CT data of an object are obtained simultaneously on a hybrid CT-MRI scanner or sequentially in a consistent status, u ^MR has a corresponding CT image u ^CT. We assume that the patches u i CT and u i MR of these two modal images in registration satisfy u i MR = Qu i CT + ε i CT(10) where ε ^CT is the error and Q is the linear transform operator. Combining Equations (9) and (10), we have u i MR = QD CT α i CT + ( Q η i + ε i CT )(11)

Equation (11) may be simplified to ∥ u i MR - QD CT α i CT ∥ 2 2 < δ(12)

Therefore, through sparse representation and dictionary learning, a high-quality patch u i MR can be sparsely coded by the same vector α i CT in the relationship of D ^MR = Q · D ^CT. It other words, the patch of MRI image u i MR can be approximately recovered by multiplying D ^MR by the sparse representation α i CT obtained from the corresponding CT image and the dictionary D ^CT if the atom mapping from the dictionary D ^CT to D ^MR is provided. That is, u i MR = D MR α i CT(13)

Next, the key is to establish an atom-to-atom mapping between these two dictionaries of different imaging modalities. Here we propose an edged-orientated method to establish such a mapping. It was proved that the public edge information between the well-registered MRI and CT images was useful to improve the MRI image quality in MRI-CT hybrid imaging [35]. Different from that method, here the edge information in the low contrast CT images is used to choose the atoms of the CT dictionary to sparsely code a given CT image, instead of the traditional Euclidean distance measure.

For a given CT image u ^CT, some sample CT images are first chosen from the initial dataset whose structure information is close to the target image. Then, a specific dictionary D ^CT may be generated from the patches of these sample images. To represent each patch u i CT of the given CT image, we perform an edge-based searching to select the closest atoms in D ^CT.

As shown in Fig. 1, (a) is the MRI image of a patient’s brain to be reconstructed, (b) is the corresponding CT image at the same slice, (d) and (e) are the MRI and CT sample images selected from the training dataset. Without loss of generality, let us consider a patch u i CT of the CT image (b) whose center locates at N ₀ (x ₀, y ₀). To find the patches in Fig. 1(e) similar to u i CT, we use their location from the boundaries instead of the Euclidean distance of the grayscale values because the contrast of this CT image of the human brain is rather low. Thus, the high contrast skull edges in Fig. 1 (b) and (e) are first segmented and extracted. Two contours can be respectively outlined by single-pixel-wide edges as shown in Fig. 1(c) and (f). The four purple points in Fig. 1(c) are the intersections of the contour and the horizontal and vertical lines through the point N ₀ (x ₀, y ₀). The distances between N ₀ and these four intersections can be calculated by { d T = y T - y 0 d B = y 0 - y B d L = x 0 - x L d R = x R - x 0(14) where x _R , x _L , x _T and y _B are the x- and y-coordinates of these four intersections, respectively.

We can compute two distance ratios of the above distances along the x- and y-axes respectively, { k x _ N 0 = d L d L + d R = x 0 - x L x R - x L k y _ N 0 = d T d T + d B = y T - y 0 y T - y B(15)

Assume that in the sample image Fig. 1(e) N ₁ (x ₁, y ₁) is the closest patch or atom to N ₀ (x ₀, y ₀). As shown in Fig. 1(f), which is the contour image of (e), the distances between N ₁ and the four intersections on the edge can be calculated using the equation similar to Equation (14). Then, two distance ratios of these distances with respect to N ₁ may be given by { k x _ N 1 = d L ′ d L ′ + d R ′ = x 1 - x L ′ x R ′ - x L ′ k y _ N 1 = d T ′ d T ′ + d B ′ = y T ′ - y 1 y T ′ - y B ′(16)

Let k _{x_N 1} = k _{x_N 0} and k _{y_N 1} = k _{y_N 0} , we have { x 1 = k x N 0 · ( x R ′ - x L ′ ) + x L ′ y 1 = y T ′ - k y N 0 · ( y T ′ - y B ′ )(17)

For any given patch u i CT of the CT image Fig. 1(b), a corresponding patch can be determined in the sample image Fig. 1(e) using Equation (17). That is, Equation (17) establishes an edge-oriented atom mapping from the patient’s CT image to the sample CT image, which selects the atoms having the same distance ratios defined by Equation (15).

However, because the values of x R ′ , x L ′ , x T ′ , and x B ′ related to the location of N ₁, (x ₁, y ₁) cannot be directly determined by one forward calculation of Equation (17), generally we need to use an iterative procedure in Fig. 2 to solve this problem, which is a crucial DDL step in our EDGE algorithm.

In this paper, the central pixel of this atom was determined by being rounded to the nearest integer from the output (x ₁, y ₁). Furthermore, as shown in Fig. 1, once the patch N ₁ (x ₁, y ₁) is selected from the sample CT image Fig. 1(e), its corresponding patch M ₁ (x ₁, y ₁) is also determined at the same position in the sample MRI image Fig. 1(d) because these two sample images are presumably well-registered.

To each patch of the known corresponding CT image Fig. 1(b), u i CT, its dictionary D i CT consists of the atoms selected by the procedure in Fig. 2 from the sample CT images. To improve the accuracy, the atom in each dictionary D i CT triples its length by adding the first-order gradient vector of each patch along x- and y- directions respectively [32]. Then, the sparse representation α i CT of each u i CT with respect to D i CT can be calculated using the OMP algorithm. Given on the coherence between the CT and MRI images, a correlative dictionary D i MR can be generated from the corresponding well-registered MRI image. Thus, each patch of the reconstructed MRI image u i MR can be approximately recovered by Equation (13) from the corresponding spare representation α i CT. Once all the patches are obtained, the final MRI image can be reconstructed by combining these patches. Figure 3 is a general workflow of the EDGE algorithm to reconstruct an initial MRI image u ^MR.

Note that the above EDGE algorithm requires that the single-pixel contour should be a convex and closed curve as shown in Fig. 1. Otherwise, there may have more than four intersections of the contour and those horizontal and vertical lines through the point N ₀ (x ₀, y ₀). Figure 4 shows the general process to get this convex contour. First, the high contrast bone edges, Fig. 4(b), can be extracted using the common threshold method from the original CT image. Then, the largest connected domain is selected with a certain area threshold as shown in Fig. 4(c). Figure 4(d) is the minimum convex set containing the connected domain of Fig. 4(c). Finally, we can extract the convex set and its single-pixel contour shown in Fig. 4(d) and (e) respectively.

Step 2. High-resolution MRI image reconstruction using DL from highly under-sampled k-space data and an initial image.

Once an MRI image is obtained by Step 1 from its counterpart CT image and two dictionaries D ^CTand D ^MR, it will be used as an initial input in the following high-resolution MRI image reconstruction from highly under-sampled k-space data. Based on the CSMRI model described in [27], the image reconstruction problem can be approximately described as min x ∑ ij ∥ R ij x - D MR α ij ∥ 2 2 + ν ∥ Fx - y MR ∥ 2 2(18) where the matrix R _ij represents the operator extracting the patch x _ij from x, the weight ν is the Lagrange multiplier which depends on the standard deviation of the noise in the measured k-space data y ^MR of MRI, F is the Fourier transform matrix. Thus, Fx represents the k-space data calculated from x.

Equation (18) is a common l ₂ norm minimization problem which can be solved directly using a classical iterative optimization technique, such as Gauss-Seidel and conjugate gradient (CG). The detailed solution procedure is summarized in Fig. 5.

In Equation (19), y (k _x , k _y ) denotes the updated value at the location (k _x , k _y ) in the k-space, S ₀ represents the original k-space data which are highly under-sampled, Ω is a subset of k-space that has been sampled, and S = FFT (x ⁽ⁿ⁾) represents the k-space data generated from the current iterative image x ⁽ⁿ⁾.

3.1 Noiseless simulation

The sample CT and MRI images in Fig. 6 were selected to generate the dictionaries of these two modalities. The images in each column have the same size and are well-registered. The images in the 3rd row are the closed single-pixel edges of the skulls in the CT images of the 1st row. Figure 7 shows the test MRI and CT images of the same patient. (a) is the test CT image. (b) is the closed single-pixel edge of the skull extracted from (a). (c) is the target MRI image to be reconstructed. (d) is the initial MRI image reconstructed by the step 1 of EDGE (Algorithm 1) from the training data set consisting of the five pairs of well-registered MRI and CT images from other patients shown in Fig. 6, and the CT image of this same imaged object shown in Fig. 7(a). (e) is the reconstructed MRI image by the traditional DDL algorithm without edge-oriented mapping [12]. According to the true value of the target MRI image, calculate the MSE (mean squared error) and SSIM of these two reconstructed results. The MSE and SSIM of Fig. 6(d) are 8720.3 and 0.8003, respectively. However, the MSE and SSIM of Fig. 6(e) are 9909.5 and 0.7043, respectively. Enhanced by the adaptive edge information, EDGE algorithm can provide much better image quality than the traditional DDL algorithm.

Figure 8 shows the final MRI reconstruction results from the 20 folds under-sampling k-space measurement data. (a) is the test MRI image. (b) is the 20 folds under-sampling matrix in k-space. (c) is the reconstruction by the Fourier transform algorithm with zero-filling of measured data. (d) is the initial MRI image reconstructed by the Step 1 of EDGE. (e) is the reconstruction by the traditional DL-MRI algorithm with the initial image (c) [27]. (f) is the reconstruction by the EDGE algorithm with the initial image (d). The MSE of the reconstructed image (f) is 0.0065 which is much lower than the MSE of the reconstructed image (e), 0.0564. It can also be seen that the EDGE reconstruction has much cleaner background and less errors than the results of DL-MRI.

Figure 9 represents the PSNR as a function of the number of iterations for the EDGE and traditional DL-MRI algorithms respectively. After 15 iterations, for the DL-MRI reconstruction the PSNR and SSIM are 12.5128 dB and 0.6292, respectively. In contrast, the PSNR and SSIM of EDGE reconstruction are 27.7277 and 0.9322. Moreover, according to the same iterative convergence condition, our algorithm can stop when the iteration number equaled 6, while DL-MRI method stopped more than 15 iterations. Figure 10 shows the PSNR curves for different under-sampling folds in k-space measurements. We can find that the result of EDGE has a slow decline. However, the result of DL-MRI has a rapid fall after 14 folds.

Benefiting from the edge-oriented patch mapping between D ^CT and D ^MR, similar structural information are preserved within these two dictionaries. The initial MRI image reconstructed by EDGE can be regard as a highly under-sampling MRI image whose under-sampling factor is about 14. Therefore, as shown in Fig. 10, when the under-sampling factor is higher than 13, EDGE keep a much better PSNR than DL-MRI. The mean PSNR of the EDGE is about 18 dB higher than the DL-MRI result when the k-space is 25 under-sampled.

3.2 Noisy data simulation

In the noisy data simulation, the k-space data were added by the Gaussian noise with zero mean and standard deviation σ. And, the k-space data was 20 folds under-sampled. Figure 11 shows the results with the noise level σ = 10 reconstructed by DL-MRI and EDGE respectively. It may be found that the reconstruction of EDGE has a higher PSNR and cleaner background than that of DL-MRI. Table 1 shows the mean PSNR values with different noise levels for EDGE and DL-MRI reconstructions. It can be seen that the PSNR values for EDGE are much better than the results of DL-MRI from either noise-free or noisy data. In the case of high level noise, increasing the number of dictionary atoms, such as 10 or 15 (each image generates two or more atoms), may provide a better reconstruction result.