PET-MRI Joint Reconstruction by Joint Sparsity Based Tight Frame Regularization

TL;DR

This paper introduces a novel PET-MRI joint reconstruction model using tight frame regularization and joint sparsity, improving image quality by accounting for different regularities of PET and MRI data.

Contribution

It proposes a non-convex joint sparsity model with a proximal alternating minimization algorithm and proves its global convergence, advancing PET-MRI image reconstruction techniques.

Findings

Achieves better reconstruction performance than existing models

Utilizes a non-convex balanced approach for different image regularities

Provides a convergence proof based on Kurdyka-Lojasiewicz property

Abstract

Recent technical advances lead to the coupling of PET and MRI scanners, enabling to acquire functional and anatomical data simultaneously. In this paper, we propose a tight frame based PET-MRI joint reconstruction model via the joint sparsity of tight frame coefficients. In addition, a non-convex balanced approach is adopted to take the different regularities of PET and MRI images into account. To solve the nonconvex and nonsmooth model, a proximal alternating minimization algorithm is proposed, and the global convergence is present based on Kurdyka-Lojasiewicz property. Finally, the numerical experiments show that the our proposed models achieve better performance over the existing PET-MRI joint reconstruction models.

Tables

Table 1. Table 1 : Comparison of relative errors, PSNR, and correlations.

Radial Sampling
Images		Indices	Initial	Analysis	DDTF	QPLS	JAnal	JSTF	JSDDTF
PET-PD	PET	RelErr	$0.3156$	$0.1112$	$0.0939$	$0.1040$	$0.0945$	$0.0937$	0.0905
		PSNR	$18.8555$	$27.9170$	$29.3820$	$28.5011$	$29.3297$	$29.4054$	29.7059
		Corr	$0.9384$	$0.9918$	$0.9942$	$0.9928$	$0.9942$	$0.9943$	0.9946
	MRI	RelErr	$0.3448$	$0.2148$	$0.2129$	$0.2878$	$0.2732$	$0.2022$	0.1929
		PSNR	$24.0252$	$28.1347$	$28.2127$	$25.5958$	$26.0482$	$28.6616$	29.0692
		Corr	$0.9018$	$0.9633$	$0.9648$	$0.9350$	$0.9410$	$0.9678$	0.9707
PET-T1	PET	RelErr	$0.3156$	$0.1112$	$0.0939$	$0.1034$	$0.0952$	$0.0934$	0.0893
		PSNR	$18.8555$	$27.9170$	$29.3820$	$28.5453$	$29.2649$	$29.4358$	29.8250
		Corr	$0.9384$	$0.9918$	$0.9942$	$0.9929$	$0.9941$	$0.9943$	0.9948
	MRI	RelErr	$0.3571$	$0.1992$	$0.1809$	$0.2266$	$0.2146$	$0.1753$	0.1663
		PSNR	$19.9312$	$25.0007$	$25.8395$	$23.8804$	$24.3556$	$26.1104$	26.5667
		Corr	$0.8967$	$0.9695$	$0.9752$	$0.9601$	$0.9641$	$0.9764$	0.9787
PET-T2	PET	RelErr	$0.3156$	$0.1112$	$0.0939$	$0.1069$	$0.0964$	$0.0923$	0.0881
		PSNR	$18.8555$	$27.9170$	$29.3820$	$28.2628$	$29.1616$	$29.5351$	29.9369
		Corr	$0.9384$	$0.9918$	$0.9942$	$0.9924$	$0.9940$	$0.9944$	0.9949
	MRI	RelErr	$0.3907$	$0.2415$	$0.2285$	$0.2678$	$0.2497$	$0.2108$	0.2023
		PSNR	$19.9442$	$24.1204$	$24.6029$	$23.2245$	$23.8311$	$25.3016$	25.6593
		Corr	$0.9005$	0.9635	$0.9676$	$0.9551$	$0.9608$	$0.9723$	0.9745
Random Sampling
Images		Indices	Initial	Analysis	DDTF	QPLS	JAnal	JSTF	JSDDTF
PET-PD	PET	RelErr	$0.3156$	$0.1112$	$0.0939$	$0.1027$	$0.0943$	$0.0933$	0.0906
		PSNR	$18.8555$	$27.9170$	$29.3820$	$28.6088$	$29.3456$	$29.4376$	29.6987
		Corr	$0.9384$	$0.9918$	$0.9942$	$0.9930$	$0.9942$	$0.9943$	0.9946
	MRI	RelErr	$0.3366$	$0.2081$	$0.2082$	$0.2729$	$0.2618$	$0.2003$	0.1919
		PSNR	$24.2355$	$28.4134$	$28.4054$	$26.0581$	$26.4169$	$28.7440$	29.1145
		Corr	$0.9064$	$0.9657$	$0.9663$	$0.9408$	$0.9452$	$0.9683$	0.9709
PET-T1	PET	RelErr	$0.3156$	$0.1112$	$0.0939$	$0.1036$	$0.0950$	$0.0929$	0.0889
		PSNR	$18.8555$	$27.9170$	$29.3820$	$28.5293$	$29.2815$	$29.4772$	29.8588
		Corr	$0.9384$	$0.9918$	$0.9942$	$0.9929$	$0.9941$	$0.9943$	0.9948
	MRI	RelErr	$0.3712$	$0.1952$	$0.1785$	$0.2180$	$0.2093$	$0.1707$	0.1627
		PSNR	$19.5948$	$25.1765$	$25.9556$	$24.2185$	$24.5722$	$26.3422$	26.7584
		Corr	$0.8878$	$0.9709$	$0.9759$	$0.9629$	$0.9658$	$0.9777$	0.9797
PET-T2	PET	RelErr	$0.3156$	$0.1112$	$0.0939$	$0.1072$	$0.0964$	$0.0928$	0.0882
		PSNR	$18.8555$	$27.9170$	$29.3820$	$28.2321$	$29.1553$	$29.4865$	29.9256
		Corr	$0.9384$	$0.9918$	$0.9942$	$0.9924$	$0.9939$	$0.9943$	0.9949
	MRI	RelErr	$0.4052$	$0.2357$	$0.2270$	$0.2551$	$0.2397$	$0.2089$	0.2016
		PSNR	$19.6260$	$24.3331$	$24.6594$	$23.6451$	$24.1873$	$25.3805$	25.6890
		Corr	$0.8925$	0.9654	$0.9680$	$0.9591$	$0.9638$	$0.9728$	0.9746

Equations

$$\underset{u_1\in {\mathcal{C}}_1,u_2\in {\mathcal{C}}_2}{\min }⟨1,A{u}_1+c⟩-⟨f,\ln (A{u}_1+c)⟩+\frac{\kappa }{2}\Vert {\mathcal{F}}_p{u}_2-g{\Vert }_2^2+\mathcal{R}(u_1,u_2),$$

$$\Vert u{\Vert }^2=\sum _{n\in \mathbb{Z}}|⟨u,{\phi }_n⟩{|}^2\text{for all}u\in \mathcal{H}.$$

$$u\in \mathcal{H}\mapsto Wu=\{⟨u,{\phi }_n⟩:n\in \mathbb{Z}\}\in {\ell }_2(\mathbb{Z}),$$

$$v\in {\ell }_2(\mathbb{Z})\mapsto W^{T}v=\sum _{n\in \mathbb{Z}}v[n]{\phi }_n\in \mathcal{H}.$$

$$u=\sum _{n\in \mathbb{Z}}⟨u,{\phi }_n⟩{\phi }_n,$$

$$({\mathcal{S}}_{q}u)[n]=(q\ast u)[n]=\sum _{k\in \mathbb{Z}}q[n-k]u[k]\text{for}u\in {\ell }_2(\mathbb{Z}).$$

$$W$$

$$W^T$$

$$\displaystyle\sum_{l=1}^{m}\sum_{k\in{\mathbb{Z}}}q_{l}[n+k]q_{l}[k]=\delta_{n}=\left\{\begin{array}[]{cl}1{}&\text{if}~{}n=0,\vspace{0.4em}\\ 0{}&\text{if}~{}n\neq 0.\end{array}\right.$$

$$\underset{v,W}{\min }\Vert v-Wu{\Vert }_2^2+{\lambda }^2\Vert v{\Vert }_0\text{subject to}{W}^{T}W=I,$$

$$u$$

$$W$$

$$v$$

$$\underset{V,D}{\min }\Vert V-D^{T}G{\Vert }_F^2+{\lambda }^2\Vert V{\Vert }_0\text{subject to}D{D}^T=I$$

$$\begin{array}{rl}& \underset{v,\{W_i{\}}_{i=1}^c}{\min }\sum _{i=1}^c{w}_i\Vert W_i{u}_i-v_i{\Vert }_2^2+\lambda \Vert v{\Vert }_{2,0}\\ & \text{subject to}{W}_i^T{W}_i=I,i=1,\cdots ,c.\end{array}$$

$$P(f|u_1)=\prod _{j=1}^M\frac{(A{u}_1+c)[j{]}^{f[j]}{e}^{-(A{u}_1+c)[j]}}{f[j]!}$$

$$P(g|u_2)\propto \prod _{j=1}^L\exp \left\{-\frac{|({\mathcal{F}}_p{u}_2)[j]-g[j]{|}^2}{2{\sigma }^2}\right\}$$

$$P(f,g|u_1,u_2)=P(f|u_1,u_2)P(g|u_1,u_2)=P(f|u_1)P(g|u_2).$$

$$\begin{array}{rl}& \underset{u_1\in {\mathcal{C}}_1,u_2\in {\mathcal{C}}_2,v}{\min }{\Phi }_1(u_1)+{\Phi }_2(u_2)+\frac{{\mu }_1}{2}\Vert W{u}_1-v_1{\Vert }_2^2+\frac{{\mu }_2}{2}\Vert W{u}_2-v_2{\Vert }_2^2+\lambda \Vert v{\Vert }_{2,0}\end{array}$$

$${\Phi }_1(u_1)$$

$${\Phi }_2(u_2)$$

$$\begin{array}{rl}& \underset{\begin{array}{c}{u}_1\in {\mathcal{C}}_1,u_2\in {\mathcal{C}}_2,\\ W_1,W_2,v\end{array}}{\min }{\Phi }_1(u_1)+{\Phi }_2(u_2)+\frac{{\mu }_1}{2}\Vert W_1{u}_1-v_1{\Vert }_2^2+\frac{{\mu }_2}{2}\Vert W_2{u}_2-v_2{\Vert }_2^2+\lambda \Vert v{\Vert }_{2,0}\\ & \text{subject to}{W}_i^T{W}_i=I,i=1,2,\end{array}$$

$$\underset{u_1\in {\mathcal{C}}_1,u_2\in {\mathcal{C}}_2}{\min }{\Phi }_1(u_1)+{\Phi }_2(u_2)+\lambda \Vert Wu{\Vert }_{2,1}$$

$$\Vert Wu{\Vert }_{2,1}={\Vert {\left(\Vert W{u}_1{\Vert }_2^2+\Vert W{u}_2{\Vert }_2^2\right)}^{1/2}\Vert }_1,u=(u_1,u_2).$$

$$\underset{u_1\in {\mathcal{C}}_1,u_2\in {\mathcal{C}}_2}{\min }{\Phi }_1(u_1)+{\Phi }_2(u_2)+\lambda {\Vert {\left(\Vert \nabla u_1{\Vert }_2^2+\Vert \nabla u_2{\Vert }_2^2\right)}^{1/2}\Vert }_1,$$

$$\begin{array}{rl}\underset{v,W^T{v}_1\in {\mathcal{C}}_1,W^T{v}_2\in {\mathcal{C}}_2}{\min }& {\Phi }_1(W^T{v}_1)+{\Phi }_2(W^T{v}_2)\\ & +\frac{{\mu }_1}{2}\Vert (I-W{W}^T)v_1{\Vert }_2^2+\frac{{\mu }_2}{2}\Vert (I-W{W}^T)v_2{\Vert }_2^2+\lambda \Vert v{\Vert }_{2,0},\end{array}$$

$$\underset{u_1\in {\mathcal{C}}_1,u_2\in {\mathcal{C}}_2}{\min }{\Phi }_1(u_1)+{\Phi }_2(u_2)+\lambda \Vert Wu{\Vert }_{2,0}$$

$$\begin{array}{rl}& \underset{v,W_1,W_2}{\min }{\mu }_1\Vert W_1{u}_1^0-v_1{\Vert }_2^2+{\mu }_2\Vert W_2{u}_2^0-v_2{\Vert }_2^2+\tilde{\lambda }\Vert v{\Vert }_{2,0}\\ & \text{subject to}{W}_i^T{W}_i=I,i=1,2,\end{array}$$

$$\begin{array}{rl}{u}_1^{k+1}& =\underset{u_1\in {\mathcal{C}}_1}{\mathrm{argmin}}{\Phi }_1(u_1)+\frac{{\mu }_1}{2}\Vert W_1^k{u}_1-v_1^k{\Vert }_2^2+\frac{{\alpha }_1^k}{2}\Vert u_1-u_1^k{\Vert }_2^2\\ u_2^{k+1}& =\underset{u_2\in {\mathcal{C}}_2}{\mathrm{argmin}}{\Phi }_2(u_2)+\frac{{\mu }_2}{2}\Vert W_2^k{u}_2-v_2^k{\Vert }_2^2+\frac{{\alpha }_2^k}{2}\Vert u_2-u_2^k{\Vert }_2^2.\end{array}$$

$$\begin{array}{rl}{W}_1^{k+1}& =\underset{W_1^T{W}_1=I}{\mathrm{argmin}}\frac{{\mu }_1}{2}\Vert W_1{u}_1^{k+1}-v_1^k{\Vert }_2^2+\frac{{\beta }_1^k}{2}\Vert W_1-W_1^k{\Vert }_2^2\\ W_2^{k+1}& =\underset{W_2^T{W}_2=I}{\mathrm{argmin}}\frac{{\mu }_2}{2}\Vert W_2{u}_2^{k+1}-v_2^k{\Vert }_2^2+\frac{{\beta }_2^k}{2}\Vert W_2-W_2^k{\Vert }_2^2.\end{array}$$

$$\begin{array}{rl}{v}^{k+1}=\underset{v}{\mathrm{argmin}}\lambda \Vert v{\Vert }_{2,0}& +\frac{{\mu }_1}{2}\Vert v_1-W_1^{k+1}{u}_1^{k+1}{\Vert }_2^2\\ & +\frac{{\mu }_2}{2}\Vert v_2-W_2^{k+1}{u}_2^{k+1}{\Vert }_2^2+\frac{{\gamma }^k}{2}\Vert v-v^k{\Vert }_2^2.\end{array}$$

$$M^j$$

Full text

{AMS}

65K15, 68U10, 90C90, 92C55

1 Introduction

Medical imaging is the technique and process of visualizing the anatomy of a body for clinical analysis and medical intervention, as well as the function of some organs and tissues. For many decades, single medical imaging modality has been available and widely used to image either function or anatomy. For example, positron emission tomography (PET) and single positron emission computed tomography (SPECT) provide information on the distribution of radioisotopes and tracers in living tissue. These distributions allow clinicians to determine certain functions of tissue [33, 35]. On the other hand, magnetic resonance imaging (MRI) and X-ray computed tomography (CT) can provide the information on anatomical structure.

Recent technical advances have allowed for the coupling of PET and MRI, leading to significant advantages over traditional PET-CT and standalone MRI [16]. The resulting PET-MRI modalities are able to acquire functional data and anatomical data simultaneously, while PET-CT modalities acquire data sequentially [16, 18, 26, 52]. Since both PET and MRI data stem from the same underlying object, it would be more desirable to explore the relation between two modality images and develop a joint reconstruction model, rather than restore images independently. Typically, the major drawback of PET is low signal-to-noise ratio and low spatial resolution in the visualized image [35]. On the other hand, MRI provides an anatomical image with superior spatial resolution, whereas the so-called $k$-space data acquisition process is so time consuming that there has been increasing demand for methods reducing the amount of acquired data without degrading the image quality [43]. However, when the $k$-space data is undersampled, the Nyquist sampling criterion is violated, and this inevitably leads to the aliasing in the reconstructed image [43]. Nevertheless, through the joint reconstruction process, the complementary information in MRI and PET images can be borrowed to improve restoration quality, especially in the case of degraded and/or incomplete data. In this paper, we only consider the two dimensional PET-MRI joint reconstruction, as the same modeling can be easily applied to the three dimensional cases, and even to other joint reconstruction problems such as SPECT-CT.

Let $u_{1}$ and $u_{2}$ denote the two dimensional PET and MRI images, and let $A$ and ${\mathscr{F}}_{p}$ denote the acquisition process operator of PET and MRI respectively. The acquired PET data is assumed to be an instance of Poisson distributed vector with expectation $Au_{1}+c$ where $c$ denotes the mean number of background counts [9, 26]. Meanwhile, the noise in MRI can be modeled as a complex white Gaussian noise [26]. Then, by the maximum likelihood principle, the PET-MRI joint reconstruction problem is formulated as the following minimization problem:

[TABLE]

where ${\mathcal{C}}_{1}$, ${\mathcal{C}}_{2}$ are feasible sets, $\kappa$ is a positive constant related to the noise level of MRI data, and ${\mathcal{R}}(u_{1},u_{2})$ is a nonseparable regularization term for $u_{1}$ and $u_{2}$.

In relation to the joint reconstruction, the above model (1) naturally arises the question of what is the proper choice of ${\mathcal{R}}(u_{1},u_{2})$. As the function cannot be independent from structure [26, 45], it is reasonable to assume that the images to be reconstructed have structural correlation [26]. In the literature, the structure of an image is mathematically modeled as a singularity such as edge and ridge, and the idea of structural correlation based joint reconstruction was first explored in the color image processing as the color channels in general share the singularities [8, 25, 29, 36, 47, 50, 53]. The successful examples include the vectorial total variation [8, 29] and the parallel level set (PLS) [25]. Motivated by the works on the color image processing, the authors in [26] recently applied joint total variation (JTV) [26, 31] and PLS priors to the PET-MRI joint reconstruction. Their analysis shows that the Gâteaux derivatives of both JTV and PLS lead to the nonlinear, inhomogeneous diffusion, which promotes the edge-enhancing effect. However, even though both JTV and PLS methods showed significant improvements over the individual reconstruction methods, the regularization term ${\mathcal{R}}(u_{1},u_{2})$ has to be further improved as the artifacts appear from the other modality image when the images have complex structures.

Wavelet frames are known to be effective in capturing image singularities from degraded images. This is due to the multiscale structure of the wavelet frame systems, short supports of the framelet functions with varied vanishing moments, and the presence of both low pass and high pass filters in the wavelet frame filter banks, which are desirable in sparsely approximating images [24]. In addition, due to the redundancy of systems, wavelet frames are more robust to errors than the (bi)orthogonal wavelet bases in [44]. Finally, the special structure of filter banks enables the construction of a so-called data driven tight frame to provide a better sparse approximation of a given input image adaptively [14, 54, 56, 57], with the computational efficiency over the K-SVD method in [1]. The successful applications in the medical image reconstruction can be found in the sparse angle X-ray CT reconstruction (wavelet frames in [21], and data driven tight frames in [56]), the limited view CT reconstruction [19], the sparse MRI reconstruction [40], and the data driven tight frame based CT reconstruction in [57], etc.

In this paper, we propose a new tight frame (wavelet frame and data driven tight frame) based PET-MRI joint reconstruction model. Our proposed model is also based on the observation that two modality images have a strong correlation in the image singularities. As the singularities are sparsely captured by the (data driven) wavelet tight frame, our model enforces the joint sparsity [28] of frame coefficients to exploit such a structural correlation. Related examples can be found in the previous work in [54], where the authors use the vector $\ell_{0}$ norm of frame coefficients to construct a multi channel data driven tight frame enforcing the structural correlation between image channels. Indeed, their applications to color image denoising and joint color-depth image reconstruction show the improvements over the independent reconstruction. Motivated by their success in color image processing, we further explore the application in the PET-MRI joint reconstruction. More precisely, using the vector $\ell_{0}$ norm of frame coefficients, we propose a balanced approach [10, 12, 17] for the joint reconstruction from the data with Poisson noise (PET) and the undersampled $k$-space data (MRI). Finally, as it can be observed that the two modality images have different regularity, we expect the use of balanced approach will help suppress the artifacts which the existing models in [26] may introduce.

The rest of this paper is organized as follows. In Section 2, we review the basic concepts of wavelet frame and data driven tight frame. In Section 3, we introduce our tight frame based PET-MRI joint reconstruction model, followed by the proximal alternating minimization algorithm, and the convergence analysis is given at the end of this section. In Section 4, we present some experimental results, and the concluding remarks are given in Section 5.

2 Preliminaries and Related Works

2.1 Wavelet Tight Frames

In this subsection, we briefly introduce the concept of tight frames and wavelet tight frames. Interested readers may consult [11, 22, 23, 46, 48] for details. Let ${\mathscr{H}}$ be a Hilbert space equipped with the inner product $\langle\cdot,\cdot\rangle$ and the norm $\|\cdot\|$. A sequence $\{\varphi_{n}:n\in{\mathbb{Z}}\}\subseteq{\mathscr{H}}$ is called a tight frame on ${\mathscr{H}}$ if

[TABLE]

Given $\{\varphi_{n}:n\in{\mathbb{Z}}\}\subseteq{\mathscr{H}}$, the analysis operator $W:{\mathscr{H}}\to\ell_{2}({\mathbb{Z}})$ is defined as

[TABLE]

and the synthesis operator $W^{T}:\ell_{2}({\mathbb{Z}})\to{\mathscr{H}}$ is defined as

[TABLE]

Thus, $\{\varphi_{n}:n\in{\mathbb{Z}}\}$ is a tight frame on ${\mathscr{H}}$ if and only if $W^{T}W=I$ with $I:{\mathscr{H}}\to{\mathscr{H}}$ being the identity. This means that for a given tight frame $\{\varphi_{n}:n\in{\mathbb{Z}}\}$, we have the following canonical expression:

[TABLE]

where $Wu=\{\langle u,\varphi_{n}\rangle:n\in{\mathbb{Z}}\}$ is called the canonical tight frame coefficients. Hence, the tight frames are generalizations of orthonormal bases to the redundant systems. In fact, a tight frame is an orthonormal basis if and only if $\|\varphi_{n}\|=1$ for all $n\in{\mathbb{Z}}$.

One of the most widely used class of tight frames is discrete wavelet frame generated by a set of finitely supported filters $\{q_{1},\cdots,q_{m}\}$. In this paper, we only discuss the undecimated wavelet frames, which is also known as the translation invariant wavelet frame transform. Given $q\in\ell_{1}({\mathbb{Z}})$, define a convolution operator ${\mathcal{S}}_{q}:\ell_{2}({\mathbb{Z}})\to\ell_{2}({\mathbb{Z}})$ by

[TABLE]

Given a set of finitely supported filters $\{q_{1},\cdots,q_{m}\}$, define the analysis operator $W$ and the synthesis operator $W^{T}$ respectively by

[TABLE]

Then, the rows of $W$ form a tight frame on $\ell_{2}({\mathbb{Z}})$ if and only if $W^{T}W=I$, i.e. the filters $\{q_{1},\cdots,q_{m}\}$ satisfy one of the so-called unitary extension principle (UEP) condition [32]:

[TABLE]

Once the one dimensional filters generate a wavelet tight frame on $\ell_{2}({\mathbb{Z}})$, the higher dimensional wavelet tight frame could be obtained via the tensor product of the one dimensional filters.

2.2 Data Driven Tight Frames

Even though tensor product wavelet frame based approach is simple to implement and able to achieve sparse representation of a piecewise smooth image, the major disadvantage is that these framelets mostly focus on singularities along the horizontal and vertical directions [14, 39]. For example, when an image has complex geometries such as directional textures, the tensor product wavelet frame coefficients may not be sparse enough [14].

In order to explore the adaptive sparse approximation of images, various data driven methods have been proposed (e.g. [1, 14]). Among these methods, one promising transformation is the so-called data driven tight frame [14]. It aims to construct a tight frame which sparsely approximates a given image $u$ adaptively. The tight frame $W=W(q_{1},\cdots,q_{r^{2}})$ is generated by finitely supported real valued filters $\{q_{1},\cdots,q_{r^{2}}\}$ satisfying (7). More concretely, given an image $u$, the data driven tight frame is obtained via solving the following minimization:

[TABLE]

where the $\ell_{0}$ norm $\|v\|_{0}$ encodes the number of nonzero entries in the coefficient vector $v$. To solve (8), we first reformulate it in the following way. Reshape all $r\times r$ patches of $u$ into $G\in{\mathbb{R}}^{r^{2}\times p}$ where $p$ denotes the number of total patches. Let $D\in{\mathbb{R}}^{r^{2}\times r^{2}}$ be the matrix generated by concatenating filters $\{q_{1},\cdots,q_{r^{2}}\}$ into column vectors $\{\vec{q}_{1},\cdots,\vec{q}_{r^{2}}\}$. Denote $V\in{\mathbb{R}}^{r^{2}\times p}$ as the frame coefficients. Hence, we have

[TABLE]

Under this setting, (8) is equivalent to

[TABLE]

where $\|\cdot\|_{F}$ denotes the Frobenius norm of a matrix. For solving (9), the alternating minimization with closed form solutions is presented in [14], and the proximal alternating minimization (PAM) scheme with global convergence property is proposed in [4].

Recently, such a patch-based adaptive construction of tight frame filters has been extended to the following data driven tight frame for multi channel image in [54]:

[TABLE]

Here, $W_{i}$, $u_{i}$, and $v_{i}$ for $i=1,2,\ldots,c$ respectively denote the tight frame filters, $i$th channel image and frame coefficients, and $\|v\|_{2,0}=\left|\{j:|v_{1}[j]|^{2}+\cdots+|v_{c}[j]|^{2}\neq 0\}\right|$ with $v=(v_{1},\cdots,v_{c})$ is the vector $\ell_{0}$ norm to promote the joint sparsity. This vector $\ell_{0}$ norm $\|\cdot\|_{2,0}$ encodes the structural correlation among the channel images, which demonstrates the advantages in multi channel image restoration [54].

3 Model and Algorithm

3.1 PET-MRI Joint Reconstruction Model

Let $u_{1}$ and $u_{2}$ respectively denote the PET and MRI images to be reconstructed. In the literature, the observed PET data $f\in{\mathbb{R}}_{+}^{M}$ is corrupted by the Poisson noise while the MRI data (or $k$-space data) $g\in{\mathbb{C}}^{L}$ is corrupted by white complex Gaussian noise. Assuming that the noise is independent for each pixel, we have

[TABLE]

where the PET forward operator $A:{\mathbb{R}}^{N}\to{\mathbb{R}}^{M}$ is in general modelled as the discrete attenuated Radon transform [5, 9], and $c$ denotes the mean number of background counts. Meanwhile, for MRI, we have

[TABLE]

where the MRI forward operator ${\mathscr{F}}_{p}:{\mathbb{R}}^{N}\to{\mathbb{C}}^{L}$ is modelled as the unitary discrete Fourier transform followed by a projection onto the measured frequencies, and $\sigma$ is the standard deviation of noise in $k$-space data. Given that $f$ and $g$ are conditionally independent on $u_{1}$ and $u_{2}$, we have

[TABLE]

Then, together with $P(u_{1},u_{2})\propto\exp(-{\mathcal{R}}(u_{1},u_{2}))$ and the Bayes’ rule, the maximum a posteriori is equivalent to the minimization problem (1).

Let $v_{1}$ and $v_{2}$ denote the wavelet frame coefficients of $u_{1}$ and $u_{2}$ under a given wavelet frame transform $W$ respectively. As the singularities of $u_{1}$ and $u_{2}$ are mostly correlated, we assume that the sparsity of $v_{1}$ and $v_{2}$ is correlated. Exploiting the idea of joint sparsity in [54], we propose our joint sparsity tight frame (JSTF) PET-MRI joint reconstruction model as follows:

[TABLE]

where

[TABLE]

and $\kappa=1/\sigma^{2}$. Here, ${\mathcal{C}}_{1}=[0,a_{1}]^{N}$ and ${\mathcal{C}}_{2}=[0,a_{2}]^{N}$ respectively denote feasible sets which reflect the physical properties of PET and MRI images.

As each modality image may contain information on the different features despite the structural correlation, we may have to actively learn the joint sparse approximation of $u_{1}$ and $u_{2}$ [54]. Hence, we also propose the following joint sparsity based data driven tight frame (JSDDTF) joint reconstruction model

[TABLE]

where the tight frames $W_{1}$ and $W_{2}$ are learned from the two modality images $u_{1}$ and $u_{2}$ respectively.

In the literature, the joint sparsity of tight frame coefficients can also be achieved via the following joint analysis (JAnal) based approach

[TABLE]

with

[TABLE]

Indeed, under properly chosen parameters, the above JAnal model (14) can be viewed as a discretization of a variational model including the following JTV model:

[TABLE]

while achieving the restoration results superior to the standard finite difference discretization of JTV (See e.g. [13] for the related details). Meanwhile, our model takes the balanced approach (e.g. [10, 12, 17]) as a special case; by fixing $u_{1}=W^{T}v_{1}$ and $u_{2}=W^{T}v_{2}$, the JSTF model (12) reduces to

[TABLE]

and the similar reasoning can be applied to the JSDDTF model (13) as well, by fixing $u_{i}=W_{i}^{T}v_{i}$ for $i=1,2$.

When $\mu_{1}$, $\mu_{2}=\infty$, the JSTF model (12) can be rewritten as the following nonconvex variant of JAnal model:

[TABLE]

as $v_{i}=Wu_{i}$ for $i=1,2$. Hence, the main difference between (12) and the JAnal model (14) lies in the contribution of $\|Wu_{i}-v_{i}\|_{2}^{2}$, which measures the distance between $v_{i}$ (frame coefficients) and $Wu_{i}$ (canonical coefficients of $u_{i}$). Since the canonical coefficients in general may not be sparse enough due to the complex geometry of human body, we use the term $\|Wu_{i}-v_{i}\|_{2}^{2}$ to provide the flexibility in sparse approximation of two modality images. Most importantly, under some mild conditions on the tight frame $W$, the canonical coefficients not only provide the information on the image singularities, but their magnitude also reflects the regularity of an image (See [7] for details). This means that, when the JAnal model (14) is used, we implicitly assume that the two modality images should have the same regularity, as well as the exact coincidence of singularities. However, even though the singularities of two modality images have a strong correlation, they may not exactly coincide with each other as different modality images reflect different physical properties of a human body. In addition, when two images have different regularity, the JAnal model can cause the singularities in one modality image to affect the smooth region in the other even if their singularities coincide. As a consequence, the JAnal model may introduce the artifacts in the reconstructed images. In contrast, $\mu_{1}$ and $\mu_{2}$ in our proposed model play a role of balancing the structural correlation of two modality images and their different regularity. As we can observe that MRI images are less noisy compared to PET image, in our proposed model, $\mu_{2}$ is chosen to be larger than $\mu_{1}$, and we expect that this different choice of $\mu_{i}$ will help to suppress such artifacts. Finally, the same arguments are applied to the JSDDTF model (13) as it has the same modeling philosophy as the JSTF model (12) except that we learn two tight frames adaptively from $u_{1}$ and $u_{2}$.

3.2 Alternating Minimization Algorithm

We propose the proximal alternating minimization (PAM) algorithms [2] to solve both (12) and (13). As the two schemes are similar except that $W_{1}$ and $W_{2}$ are additionally updated in (13), we only consider the algorithm for solving (13). Given initializations $u_{1}^{0}$ and $u_{2}^{0}$, the initializations of $v$, $W_{1}$, and $W_{2}$ are obtained via

[TABLE]

using the algorithm in [14]. After the initializations, we optimize $\{u_{1},u_{2},v,W_{1},W_{2}\}$ by solving the model (13) alternatively. The full details are described in Algorithm 1. Note that we consider the subproblems separately whenever they are separable.

It is easy to see that each problem in (18) is strongly convex and smooth, and there are numerous algorithms to solve it. In our algorithm, $u_{1}^{k+1}$ is updated using the projected scaled gradient method [34] which is a special case of preconditioned alternating projection algorithm [37]; Let $u_{1}^{k,0}=u_{1}^{k}$. For $j=0$, $1$, $2$, $\cdots$

[TABLE]

Similarly, $u_{2}^{k+1}$ is updated using the projected gradient method [42]; Let $u_{2}^{k,0}=u_{2}^{k}$. For $j=0$, $1$, $2$, $\cdots$

[TABLE]

where ${\mathscr{F}}_{p}^{*}$ is the conjugate transpose of ${\mathscr{F}}_{p}$. In any case, we omit the outer iteration $k$ on $M$, $u_{1}$, and $u_{2}$. Since the constraints ${\mathcal{C}}_{1}$ and ${\mathcal{C}}_{2}$ are convex, the above iterations converge to the global minimizers for appropriately chosen $\rho_{1}^{j}$ and $\rho_{2}^{j}$ [34, 42].

To solve (19) and (20), we introduce

[TABLE]

using $r\times r$ patches of $u_{1}$ and $u_{2}$. Under the above reformulation, (19) and (20) become

[TABLE]

and

[TABLE]

where $V=(V_{1},V_{2})$. Hence, to solve (21), we use the following closed form formulae:

[TABLE]

The closed form solution to (22) is given by

[TABLE]

where ${\mathcal{T}}_{\lambda,\mu}$ with $\mu=(\mu_{1},\mu_{2})$ is the generalized hard thresholding formula [54] defined as

[TABLE]

Here, $({\mathcal{T}}_{\lambda,w}[(U_{1},U_{2})])_{j}$ and $U_{i,j}$ denote the $j$th row vector of ${\mathcal{T}}_{\lambda,w}[(U_{1},U_{2})]$ and $U_{i}$ respectively.

3.3 Convergence Analysis

This section is devoted to the convergence analysis of our alternating minimization algorithm. Due to the same reason to the PAM algorithm, we focus on the convergence analysis of sequence $\{(u_{1}^{k},u_{2}^{k},W_{1}^{k},W_{2}^{k},v^{k}):k\in{\mathbb{N}}\}$ in (13) generated by Algorithm 1. To do this, we assume the followings:

1. A1.

$f>0$ and $A$ satisfies

[TABLE]

where $A[i,j]$ denotes the $(i,j)$th entry of $A$. 2. A2.

$c>0$ is a constant vector. 3. A3.

There exist $0<L<U$ such that $L\leq\alpha_{1}^{k},\alpha_{2}^{k},\beta_{1}^{k},\beta_{2}^{k},\gamma^{k}\leq U$ for all $k\in{\mathbb{N}}$.

We first introduce some basic notation and definitions.

Definition 3.1**.**

Let $f:{\mathbb{R}}^{n}\to{\mathbb{R}}\cup\{\infty\}$ be a proper and lower semicontinuous (lsc) function.

The domain of $f$, denoted as $\mathrm{dom}(f)$ is defined as

[TABLE] 2. 2.

For each $x\in\mathrm{dom}(f)$, the Fréchet subdifferential of $f$ at $x$ is defined as

[TABLE]

If $x\notin\mathrm{dom}(f)$, then we set $\partial_{F}f(x)=\emptyset$. 3. 3.

The (limiting-) subdifferential of $f$ at $x$ is defined as

[TABLE]

The domain of $\partial f$ is defined as \mathrm{dom}(\partial f)=\{x\in{\mathbb{R}}^{n}:\partial f(x)\neq\emptyset\big{\}}. 4. 4.

$x\in\mathrm{dom}(f)$* is a critical point of $f$ if $0\in\partial f(x)$.*

In the literature, the global convergence of the alternating schemes are established on the framework of the Kurdyka-Lojasiewicz (KL) property in [38, 41]. Even though it is challenging to verify whether a given function $f$ satisfies the KL property, there are some functions of special cases which have proven to satisfy the KL property, such as analytic functions and semi-algebraic functions (e.g. [2, 3, 6, 55] ).

To apply the framework of KL property to our model, we let

[TABLE]

for notational simplicity. Recall ${\mathcal{C}}_{1}=[0,a_{1}]^{N}$ and ${\mathcal{C}}_{2}=[0,a_{2}]^{N}$, and let ${\mathcal{D}}=\{W\in{\mathbb{R}}^{r^{2}N\times N}:W^{T}W=I_{r^{2}\times r^{2}}\}$. Then we define

[TABLE]

where $\imath_{{\mathcal{A}}}$ is the indicator function of a set ${\mathcal{A}}$: $\imath_{{\mathcal{A}}}(x)=0$ if $x\in{\mathcal{A}}$, $\imath_{{\mathcal{A}}}(x)=\infty$ otherwise.

Using the above notation, we reformulate the model (13) as

[TABLE]

Note that it is not hard to verify that $\Phi_{1}(u_{1})$ and the indicator functions are analytic, and $r_{5}(x_{5})$ is a semi-algebraic function. Together with that the remaining terms are all polynomials, we can see that the functional $H$ defined as in (28) satisfies the KL property.

Under this reformulation, we can present the following result on the global convergence.

Theorem 3.2**.**

Assume that A1-A3 hold. Let $H(x)$ be the objective function defined as in (28). Then the sequence $\{x^{k}:k\in{\mathbb{N}}\}$ generated by Algorithm 1 converges to a critical point of $H$. Moreover, $\{x^{k}\}_{k\in{\mathbb{N}}}$ has the following finite length property:

[TABLE]

The proof of Theorem 3.2 is similar to the framework given in [3, Theorem 3.7.]. In fact, $H(x)$ defined as in (28) satisfies KL property, and noting that $\Phi_{1}(x_{1})$ is analytic, proper, lsc, strictly convex and coercive [27], it is not hard to see that $H(x)$ satisfies the first condition in [3, Theorem 3.7.]. (See [3] for the details). Hence, the proof is completed provided that we show that the sequence $\{x^{k}:k\in{\mathbb{N}}\}$ generated by Algorithm 1 is bounded. In fact, since the constraints ${\mathcal{C}}_{1}$, ${\mathcal{C}}_{2}$ and ${\mathcal{D}}$ are compact sets, it suffices to show that $v^{k}$ is bounded.

Lemma 3.3**.**

*Assume that A1-A3 hold. Let $H(x)$ be the objective function defined as in (28). For $\{x^{k}:k\in{\mathbb{N}}\}$ generated by Algorithm 1, there exist $R_{1}$, $R_{2}>0$ such that $\|v_{1}^{k}\|_{2}\leq R_{1}$ and $\|v_{2}^{k}\|_{2}\leq R_{2}$ for all $k\geq 0$. *

Proof 3.4**.**

Note that the constraint on $u_{1}$ and $u_{2}$ leads to

[TABLE]

for all $k\geq 0$. If we choose $R_{1}=a_{1}\sqrt{N}$ and $R_{2}=a_{2}\sqrt{N}$, then the proof is completed by mathematical induction. For $k=0$, (27) implies

[TABLE]

where the equality comes from the fact that $W_{1}^{0}$ and $W_{2}^{0}$ are tight frames.

For the mathematical induction, we assume that $\|v_{1}^{k}\|_{2}\leq R$ and $\|v_{2}^{k}\|_{2}\leq R$ for $k\geq 0$. Again, by (27), we have

[TABLE]

*from the fact that $W_{1}^{k}$, $W_{2}^{k}\in{\mathcal{D}}$ for all $k\in{\mathbb{N}}$. This completes the proof. *

4 Numerical Results

In this section, we present some experimental results to compare our JSTF model (12) and JSDDFT model (13) with the several existing methods. We choose to compare with the the following analysis based individual reconstruction models:

[TABLE]

both of which are solved by the split Bregman method [15], the data driven tight frame (DDTF) individual reconstruction models:

[TABLE]

solved by the PAM algorithm, and the quadratic parallel level set (QPLS) method in [26]. We also compare with the JAnal model (14) solved by the split Bregman method, as a replacement of the JTV model (15).

The experiments are conducted with $256\times 256$ PET and MRI images taking values in $[0,1]$, which are available at Harvard Whole Brain Atlas webpage111http://www.med.harvard.edu/AANLIB/home.html. Throughout this paper, we set both ${\mathcal{C}}_{1}$ and ${\mathcal{C}}_{2}$ are $[0,1]^{N}$. Note that MRI images are classified into the positron density (PD) weighted image, the T1 weighted image, and the T2 weighted image according to the pulse sequence design [30] (See Fig. 1). In this sequel, we perform experiments using different MRI images for different $u_{2}$, which we shall refer to PET-PD, PET-T1, and PET-T2 respectively. To measure the quality of restored images, we compute the relative error, the PSNR, and the correlation between the reconstructed image $\widetilde{u}_{i}$ and the true image $u_{i}$ respectively defined as

[TABLE]

where $\overline{u}_{i}$ and $\overline{\widetilde{u}}_{i}$ denote the mean of $u_{i}$ and $\widetilde{u}_{i}$ respectively.

To synthesize the data, we generate the forward PET operator $A$ as described in [5], and we set $c=1$ in all cases. To generate $256\times 256$ PET data $f$ with Poisson noise, we use the MATLAB built-in function “imnoise($\cdot$,‘poisson’)”by scaling $Au_{1}+c$ with a suitable factor before applying “imnoise”, and then scaling it back with the same factor. More precisely,

[TABLE]

and we set factor to be $2*10^{8}$ in all cases. Note that the larger “factor” is, the larger noise level is. To generate $g$, we first generate ${\mathscr{F}}_{p}=R_{\Lambda}{\mathscr{F}}$, where $R_{\Lambda}$ is the projection on the known frequency region $\Lambda$, and ${\mathscr{F}}$ is the unitary discrete Fourier transform. Here, we choose $R_{\Lambda}$ to be the sampling along the $30$ radial lines and the $10\%$ random sampling described in [20]. The gaussian noise with standard deviation $0.05$ is also added to both real and imaginary part of ${\mathscr{F}}_{p}u_{2}$.

In all experiments, we choose the same initializations for the fair comparison; $u_{1}^{0}$ is obtained by the Expectation-Maximization algorithm [49], and $u_{2}^{0}$ is obtained by the inverse discrete Fourier transform with zero filling (See Fig. 2). For the individual reconstruction models (30) and (31), the JAnal Model (14) as well as our JSTF model (12), we use one level piecewise cubic B-spline wavelet frame transformation. For JSDDTF model (13) as well as DDTF individual reconstruction models (32) and (33), we use $8\times 8$ undecimated discrete cosine transform filters [51] for the initial guess of (17). In all cases, we set $\kappa=1$. In (12), (13), (32), and (33), we set $\mu_{1}=0.05$, $\mu_{2}=1$, $\alpha_{1}^{k}=\alpha_{2}^{k}=0.001$, $\beta_{1}^{k}=\beta_{2}^{k}=\gamma^{k}=0.00005$, and $\rho_{1}^{j}=\rho_{2}^{j}=0.5$ in all cases. All parameters related to the (joint) sparsity are manually chosen so that we obtain the optimal restoration results, especially the optimal results in PET image for the joint reconstruction models. In addition, as the wavelet frame systems consist of low pass filter and high pass filters and the low pass filter coefficients are not sparse in general [13], we follow the convention that we do not penalize the frame coefficients corresponding to the low pass filter.

Table 1 summarizes relative errors, PSNR values, and correlations of the aforementioned five restoration models, and Fig. 3-Fig. 8 present visual comparisons of the results. We can see that both JSTF (12) and JSDDTF model (13) consistently outperforms both the individual reconstruction models and the existing joint reconstruction models in [26]. This verifies that there exists a correlation on image singularities between the two modality images, and exploiting this correlation results in the better reconstruction results. In addition, compared to the existing joint reconstruction models, we can see that our proposed models introduce less artifacts in both modality images, leading to the visual improvements over the existing joint reconstruction methods which are consistent with the improvements in indices. It is worth noting that our proposed model can improve the restoration qualities even using a static wavelet tight frame. This demonstrates that the improvements mainly come from simultaneously considering the different regularity of the different modality images and the joint sparsity. Meanwhile, we can see that the QPLS model and the JAnal model (14), which only take the structural correlation into account, show the degradations in MRI restoration results compared to the independent reconstruction methods due to the artifacts. Finally, we list some zoom-in views in Fig. 9 to illustrate that our models (12) and (13) restore structures better than the existing methods.

5 Conclusions

In this paper, we proposed a joint sparsity based tight frame regularization PET-MRI joint reconstruction model, together with a proximal alternating minimization algorithm. The numerical experiments show that our models (12) and (13) both outperform the existing models in [26]. This performance gain of our proposed models mainly comes from taking the different regularity of the different modality images into the consideration, as well as the structural correlation. Finally, our convergence analysis demonstrates that the sequence generated by our algorithm globally converges to a critical point of the proposed model.

Acknowledgments

We would like to thank the anonymous reviewers for their constructive suggestions and comments that helped tremendously with improving the presentations of this paper.