2.2. Reconstruction method

To incorporate the trained DAE as prior information, we formulated the following reconstruction model:

x represents the MR image to be reconstructed, and g denotes the acquired signal in k space. F indicates the 2D Fourier transform operator, while S is the undersampling operator. || . || _F denotes Frobenius norm. P _i refers to the operator that extracts a 64-by-64 pixel image patch centered at the i -th pixel from x Ω is the set of all pixels in the image. The first term in Eq. (3) enforced the fidelity between the reconstructed image and the acquired k-space data. The second term incorporated the trained DAE Φ(·|θ^) to regularize the quality of the reconstructed image. β is a regularization parameter balancing the contributions from the two terms. By solving this optimization problem, we expect to reconstruct an MR image of similar image quality as the prior image, while its content is defined by the measurement g , reflecting the anatomy/structure at the time of scanning instead of being biased towards the prior information.

To solve Eq. (3), we employed the forward-backward splitting algorithm. More specifically, the original problem was solved by tackling the following two subproblems alternatively in each iteration until convergence:

The first subproblem is of a quadratic form, which can be solved efficiently using the conjugate gradient least-square algorithm (46) with the initial guess set as x ^(t) , i.e. the solution obtained in the previous iteration step. In this study, we assumed the solution is a real image and hence enforced this at this step of the iterative process by taking the absolute value.

The second subproblem in Eq (4), on the other hand, involves the DAE, which is a highly non-linear and non-convex function. To solve such a complex optimization problem, we proposed to employ a fix-point scheme, i.e. at each iteration, we fix Φ(Pix|θ^) as Φ(Pix(t+12)|θ^) and hence the second subproblem becomes

The solution to this modified optimization problem can be obtained explicitly as

where Pi* is the adjoint operator of P _i , which simply places the extracted patch i back to its location in the MR image. In this case, Pi*Pix=n2x , where n is the patch size. Note that we used one patch for each pixel, so the patches were overlapping and Each pixel was covered by n ² patches. In the above equation, the summation indicated that the adjoint operator added each patch to corresponding locations, and overlapping patches were summed. The normalization was reflected in the denominator term of this equation. In general, if patches were overlapped, the updated pixel value would be the average over all patches covering this pixel, together with the x(t+12) term.

The complete update scheme to solve the reconstruction problem in 3 is summarized in Algorithm 1 .

The interpretability of the reconstruction Algorithm 1 can be illustrated in Figure 1A , which offers a geometric visualization about the iterative reconstruction process. At the beginning of each iteration, with the solution x ^(t) as the initial guess, solving the least square problem in step 1 generates a new solution x(t+12) . For each patch xi=Pix(t+12) , the DAE generates a corresponding one on the manifold by first computing the coordinate in the latent space E(xi|θ^E) and then restoring the patch for this coordinate x^i=D(E(xi|θ^E)|θ^D) . These patches are assembled to form a new image. Finally, the new image and x(t+12) are averaged with weights proportional to β and 1 , respectively, yielding a new solution x ^(t+1) . This process continued until convergence, where a solution on, or sufficiently close to, the manifold prior and meeting the data fidelity was generated.