Joint image edge reconstruction and its application in multicontrast MRI
Abstract
We propose a new joint image reconstruction method by recovering edge directly from observed data. More specifically, we reformulate joint image reconstruction with vectorial totalvariation regularization as an minimization problem of the Jacobian of the underlying multimodality or multicontrast images. Derivation of data fidelity for Jacobian and transformation of noise distribution are also detailed. The new minimization problem yields an optimal convergence rate, where is the iteration number, and the periteration cost is low thanks to the closeform matrixvalued shrinkage. We conducted numerical tests on a number multicontrast magnetic resonance image (MRI) datasets, which show that the proposed method significantly improves reconstruction efficiency and accuracy compared to the stateofthearts.
Keywords. Joint image reconstruction, matrix norm, optimal gradient method, multicontrast.
1 Introduction
The advances of medical imaging technology have allowed simultaneous data acquisition for multimodality images such as PETCT, PETMRI and multicontrast images such as T1/T2/PD MRI. Multimodality/contrast imaging integrates two or more imaging modalities (or contrasts) into one system to produce more comprehensive observations of the subject. Such technology combines the strengths of different imaging modalities/contrasts in clinical diagnostic imaging, and hence can be much more precise and effective than conventional imaging. The images from different modalities or contrasts complement each other and generate high spatial resolution and better tissue contrast. However, effectively utilizing sharable information from different modalities and reconstructing multimodality images remain as a challenging task especially when only limited data are available. Therefore, the main approach of joint image reconstruction is to incorporate similarities, such as anatomical structure, between all modalities/contrasts into the process to improve reconstruction accuracy [2, 9, 17].
In this paper, we consider a new approach by jointly reconstructing edges of images, rather than the images themselves, directly from multimodality/contrast imaging data. The final images can be obtained very easily given the reconstructed edges. We show that this new approach results in an type minimization that can be effectively solved by accelerated proxgradient methods with an optimal convergence rate, where is the iteration number. Moreover, the subproblems reduce complex vectorial/joint total variation regularizations to simple matrixvalued shrinkages, which often have cheap closedform solutions. This is in sharp contrast to primaldual based image reconstruction algorithms, such as primaldual hybrid gradient method (PDHG) and alternating direction method of multipliers (ADMM), which only yield convergence rate. Therefore, the proposed method can produce high quality reconstructions with much improved efficiency over the stateofthearts methods in multimodality/contrast image reconstruction problems.
The contributions of this paper can be summarized as follows. (a) We develop a novel twostep joint image reconstruction method that transforms the vectorial TV regularized minimization of image into an minimization of Jacobian. This enables a numerical scheme with optimal convergence rate by employing accelerated gradient descent method (Section 3.2). (b) In the resulting minimization problems, the main subproblem involves matrix norm (instead of vectorial TV) and can be solved easily using generalized shrinkage for matrices. We provide closeform solutions for several cases originated from commonly used vectorial TVs (Section 3.4 and Appendix A). (c) We analyze the noise distribution after transformation, and incorporate it into the data fidelity of Jacobian in the algorithm which significantly improves reconstruction accuracy and efficiency (Section 3.3). (d) We conduct a series to numerical tests on several multicontrast MRI datasets and show the very promising performance of the proposed methods (Section 4). Although we only focus on multicontrast MRI reconstruction where all data is acquired in Fourier space, it is worth noting that the proposed method can be readily extended to other cases, such as those with Radon data, hence is applicable to reconstructions involving other types of imaging modalities.
The remainder of the paper is organized as follows. We first provide an overview of recent literatures in joint image reconstruction in Section 2. Then we present the proposed method and address a number of details in Section 3. In Section 4, we conduct a number of numerical tests on a variety of multicontrast MR image reconstruction datasets. Section 5 concludes our findings.
2 Related Work
There have been a significant amount effort devoted to develop models and algorithms that can effectively take the anatomy structure similarities across modalities/contrasts into account during joint reconstructions. As widely accepted, these structure similarities can be exploited using the locations and directions of the edges [7, 10, 11] characterized by the magnitude and direction of the gradient of an image. The active researches on multimodal/contrast image reconstructions have focused mainly on how to effectively utilize the complimentary information on these structure similarities to improve the accuracy and robustness of the reconstructions.
Inspired by the success of total variation (TV) based image reconstructions for scalarvalued images, many algorithms for joint reconstruction of multimodal images extend the classical TV to vectorial TV for joint multimodal image reconstructions. This extension aims at capturing the sharable edge information and performing smoothing along the common edges across the modalities. There have been several different ways to define the TV regularization for multimodal images (represented by vectorvalued functions). For instance, , where is the (scalarvalued) image of modality/channel/contrast in an modality imaging problem, is proposed in [3]. In [6], a number of variations along this direction, called collaborative TV (CTV), are studied and summarized comprehensively. A specific TV takes a particular form to integrate partial derivatives of the image across the modalities. For example, takes norm of image gradient at every point (pixel) for each modality, then norm over all pixels, and finally norm (direct sum) across all modalities. Another commonly used CTV variation, called joint TV (JTV), is formulated as and has been successfully applied to color image reconstruction [4, 23], multicontrast MR image reconstruction [17, 19], joint PETMRI reconstruction [11], and in geophysics [16]. As one can see, JTV is taking norm in terms of gradients, norm across modalities, and then finally for pixels. In general, the gradient of a vectorvalued image consisting multiple modalities/contrasts is a tensor at each pixel, and the TV takes specific combination of norms regarding partial derivatives (gradients), modalities (channels, or contrasts), and pixels, respectively.
Another joint regularization approach is to extend anisotropic TV regularization from scalar to vectorvalued case. In this approach, anisotropic TV is employed to replace isotropic TV by an anisotropic term that incorporates directional structures exhibited by either the original data or the underlying image. Anisotropic TV been applied in standard singlemodal image reconstruction with successes in [13, 15, 18], where structure tensor is used to provide information of both size and orientation of image gradients. In particular, the model proposed in [15] employs an anisotropic TV regularization of form for image reconstruction. Here, is the gradient of the singlemodal, scalarvalued 2D image , the diffusion tensor is determined by the eigenvalues and corresponding eigenvectors of the structure tensor . The model developed in [13] alternately solves two minimizations: one that estimates structure tensor using the image from previous iteration, and the other one improves the image with an adaptive regularizer defined from this tensor. The idea of anisotropic TV has been extended to vectorvalued images with the anisotropy adopted from the gradients of multimodal/contrast images, which incorporates similarities of directional structures in joint reconstruction. The method developed in [9] projects the gradient in the total variation functional onto a predefined vector field given by the other contrast for joint reconstruction of multicontrast MR images. In [9], a directional TV (DTV) is proposed with the form , where is the residue of projection of to , and and are scalarvalued functions each representing the image of one modality and for all . In other words, this is the anisotropic diffusion by using the structure tensor of given reference image .
With geometric interpretation of gradient tensor above, authors in [5] suggest to consider a vectorvalued image as a parameterized 2dimensional Riemann manifold with metric , where is a vectorvalued image and is the Jacobian of . Then the eigenvector corresponding to the larger eigenvalue gives the direction of the vectorial edge. Based on this framework, several forms of vectorial TV (VTV) have been developed. In [22], a family of VTV formulations are suggested as the integral of over the manifold, where and denote the larger and smaller eigenvalues of the metric , respectively, and is a suitable scalarvalued function. A special choice of , i.e., the Frobenius norm of the Jacobian , reduce to JTV mentioned above. For another special choice of , there is , where is the largest singular value of the Jacobian of . This can be computed by using the dual formulation with where is the standard simplex in .
Besides joint TV or tensor based regularization, the parallelism of level sets across multimodality/contrast images are also proposed as joint image regularization in [8, 11, 12, 16]. The main idea is to exploit the structural similarity of two images and measured by the parallelism of their gradients and at each point using , for some functions and . Then the regularization in joint reconstruction takes form Several different choices of and have been studied in these works. For instance, in [12], and are taken as identities, or and . In [16], is the identity and . In [11] the side information on the level set, namely, the location and direction of the gradients of a reference image, is available to assist the reconstruction.
Another joint reconstruction approach different from aforementioned methods is to recover gradient information for each of the underlying multimodal images from the measurements in Fourier domain, then use the recovered gradients to reconstruct images. This approach is motivated by the idea of Bayesian compressed sensing and applied to joint multicontrast MR image reconstruction in [2]. In [2], gradients of the multicontrasts images are reconstructed jointly from their measurements in Fourier space under a hierarchical Bayesian framework, where the joint sparsity on gradients across multicontrast MRI is exploited by sharing the hyperparameter in the their maximum a posteriori (MAP) estimations. Their experiments show the advantage of using joint sparsity on gradients over conventional sparsity. However, their method requires extensive computational cost. A twostep gradient reconstruction of MR images is also proposed in [20], however, only singlemodality/contrast image is considered. In [20], the authors showed that this twostep gradient reconstruction approach allows to reconstruct image with fewer number of measurements than required by standard TV minimization method.
3 Proposed Method
In this section, we propose a new joint image reconstruction method that first restores image edges (gradients/Jacobian) from observed data, and then assembles the final image using these edges. Without loss of generality, we assume all images are dimensional, i.e., the image domain (all derivations below can be easily extended to higher dimensional images). For simplicity, we further assume is rectangular. For single channel/modality/contrast case, we use function represents the image, such that stands for the intensity of image at . In multimodality case, where is the number of modalities.
It is also convenient to treat an image as a matrix in discretized setting which we mainly work on in practice, and further stack the columns of to form a single column vector in , where is the total number of pixels in the image. For multimodality case, we have to represent the (discretized) image of modality for .
3.1 Vectorial totalvariation regularization
Standard totalvariation (TV) regularized image reconstruction can be formulated as a minimization problem as follows:
(1) 
where represents the data fidelity function, e.g., , for some given data sensing matrix and observed partial/noisy/blurry data . By solving (1), we obtain a solution which minimizes the sum of TV regularization term and data fidelity term with a weighting parameter that balances the two terms. It is shown that TV regularization can effectively recover high quality images with well preserved object boundaries from limited and/or noisy data.
In joint multimodality image reconstruction, the edges of images from different modalities are highly correlated. To take such factor into consideration, the standard TV regularized image reconstruction (1) can be simply replaced by vectorial TV (VTV) regularized counterpart:
(2) 
where is the vectorial TV of . In the case that is continuously differentiable, VTV is a direct extension of standard TV as an “ of gradient” as
(3) 
where is the Jacobian matrix at point , and indicates some specific matrix norm. For example, let be an by matrix, then we may use one of the following matrix norms as :

Frobenius norm: . This essentially treats as an vector.

Induced 2norm: where is the largest singular value of . This norm is advocated in [14] with a geometric interpretation when used in VTV.

Nuclear norm: where are singular values of . This is a convex relaxation of matrix rank.
Obviously, there are a number of other choices for VTV due to the many variations of matrix norms. However, in this paper, we only focus on these three norms as they are the mostly used ones in VTV regularized image reconstructions.
It is also worth noting that, in general, the VTV norm with matrix norm is defined for any function (not necessarily differentiable) as
(4) 
where is the dual norm of . Although we would not always make use of the original VTV definition (4) in discrete setting, we show that they can help to derive closedform softshrinkage with respect to the corresponding matrix norm as in Appendix A.
3.2 Joint edge reconstruction
The main idea of this paper is to reconstruct edges (gradients/Jacobian) of multimodality images jointly, and then assemble the final image from these edges. To that end, we let denote the “gradients” (or Jacobian) of . If is differentiable, then is the Jacobian matrix of at . For example, assuming there are three modalities and at every point , then is a matrix
(5) 
where at every for and ( here). As a result, the VTV of in (3) simplifies to . If is not differentiable, then may become the weak gradients of (when ) or even a RadonNikodym measure (when ). However, in common practice using finite differences for numerical implementation to approximate partial derivatives of functions, such subtlety does not make much differences. Therefore, we treat as the Jacobian of throughout the rest of the paper and in numerical experiments.
Assume that we can derive the relation of the Jacobian and the original observed data and form a data fidelity of , as an analogue of data fidelity of image (justification of this assumption will be presented in Section 3.3). Then we can reformulate the VTV regularized image reconstruction problem (2) about to a matrix norm regularized inverse problem about as follows:
(6) 
Now we seek for reconstructing instead of . This reformulation (6) has two significant advantages compared to the original formulation (2):

It reduces to an type minimization (6) and can be solved effectively by accelerated gradient method. For example, using Algorithm 1 based on FISTA [1], we can attain an optimal convergence rate of to solve (6), where is the iteration number. This is in sharp contrast to the best known rate of primaldual based methods (including the recent, successful PDHG and ADMM) for solving (2).
3.3 Formulating fidelity of Jacobian
In an extensive variety of imaging technologies, especially medical imaging, the image data are acquired in transform domains. Among those common transforms, Fourier transform and Radon transform are the two most widely used ones in medical imaging. In what follows, we show that the relation between an image and its (undersampled) data can be easily converted to that between the Jacobian matrix and in the Fourier domains. This allows us to derive the data fidelity of Jacobian and reconstruct it properly for MRI reconstruction problems. Similar idea can be applied to the case with Radon transforms, subject to modification of data fidelity and noise distribution in a straightforward manner. In this paper, we focus on the Fourier case only as our data are multicontrast MRIs.
Data transformation. Imaging technologies, such as magnetic resonance imaging (MRI) and radar imaging, are based on Fourier transform of images. The image data are the Fourier coefficients of image acquired in the Fourier domain (also known as the space in MRI community). The inverse problem of image reconstruction in such technologies often refers to recovering image from partial (i.e., undersampled) Fourier data.
In discrete setting, let denote the image to be reconstructed, the (discrete) Fourier transform matrix (hence a unitary matrix), and the diagonal matrix with binary values ( or ) as diagonal entries to represent the undersampling pattern (also called mask in space), and be the observed partial data with at unsampled locations. Therefore, the relation between the underlying image and observed partial data is given by .
The gradient (partial derivatives) of an image can be regarded as a convolution. The Fourier transform, on the other hand, is well known to convert a convolution to simple pointwise multiplication in the transform domain. In discrete settings, this simply means that is a diagonal matrix where is the discrete partial derivative (e.g., forward finite difference) operator along the direction (). This amounts to a straightforward formulation for the data fidelity of : Let be the partial derivative of in direction, i.e., then we have
(10) 
where we used the facts that both and are diagonal matrices so they can commute. Therefore, the data fidelity term in (6) can be, for example, formulated as
(11) 
As long as the Fourier data is concerned, the data fidelity of the corresponding gradient can be formulated in a similar way as above.
Noise transformation. We now consider the noise distribution in image reconstruction from Fourier data. Suppose the noise is due to acquisition in Fourier transform domain such that
(12) 
for each frequency value in Fourier domain. Here is the Fourier transform of image . By multiplying on both sides of (12), we see that this fidelity of can be obtained by
Suppose that where and are real and imaginary parts of respectively, and where and are independent for all . Then we know that the transformed noise are independent for different and distributed as bivariate Gaussian
(13) 
where and . Therefore, we can readily obtain the maximum likelihood of and hence the fidelity of . In particular, if , then and are two i.i.d. . i.e., . In this case, we denote , then data fidelity (i.e., negative loglikelihood) of simply becomes , as in contrast to the fidelity term of . In the remainder of this paper, we assume that the standard deviation of real and imaginary parts are both . Furthermore, in numerical experiments, we can precompute and only perform an additional pointwise division of after computing in each iteration.
3.4 Closed form solution of matrixvalued shrinkage
As we can see, the algorithm step (7) calls for solution of type
(14) 
for specific matrix norm and given matrix . In what follows, we provide the closeform solutions of (14) when the matrix norm is Frobenius, induced 2norm or nuclear norm, as mentioned in Section 3.1. The derivations are provided in Appendix A.

Nuclear norm. This norm promotes low rank and yields a closeform solution of (14) as
(17) where is the singular value decomposition (SVD) of .
The computation of (15) is essential shrinkage of vector and hence very cheap. The computations of (16) and (17) involve (reduced) SVD, however, explicit formula also exists as the matrices have tiny size of by (by if images are 3D), where is the number of image channels/modalities.
Remark. It is worth noting that the computation of (14) is carried out at every pixel independently of others in each iterations. This allows straightforward parallel computing which can further reduce realworld computation time.
3.5 Reconstruct image from gradients
Once the Jacobian is reconstructed from data, the final step is to resemble the image from . Since this step is performed for each modality, the problem reduces to reconstruction of a scalarvalued image from its gradient . In [20, 21], the image is reconstructed by solving the Poisson equation since and are the partial derivatives of . The boundary condition of this Poisson equation can be either Dirichlet or Neumann depending on the property of imaging modality. In medical imaging applications, such as MRI, CT, and PET, the boundary condition is simply since it often is just background near image boundary .
Numerically, it is more straightforward to recover by solving the following minimization
(18) 
with some parameter to weight the data fidelity of . The solution is easy to compute since the objective is smooth, and often times closedform solution may exist. For example, in the MRI case where , the solution of (18) is given by
(19) 
where is diagonal and hence trivial to invert. The main computations are just few Fourier transforms.
In fact, it seems often sufficient to retrieve the base intensity of from to reconstruct from , and hence the result is not sensitive to when solving (18). It is also worth pointing out that, by setting , the minimization (18) is just least squares whose normal equation is the Poisson equation mentioned above.
To summarize, we propose the twostep Algorithm 2 which restores image edge and resembles image for multicontrast MRI reconstruction. The two steps are each executed only once (no iteration). Step 1 itself requires iterations which converge very quickly with rate where is iteration number. Step 2 has closed form solution (19) for multicontrast MRI reconstruction.
4 Numerical Results
In this section, we conduct a series of numerical experiments on synthetic and real multicontrast MRI datasets using Algorithm 2 (for short, we call it ER, standing for Edgebased Reconstruction, without adaptive weighting, and ERweighted for the one with adaptive weighting in Section 3.3). For comparison, we also obtained implementation of two stateofthearts method for joint image reconstruction: multicontrast MRI method Bayesian CS (BCS)^{4}^{4}4BCS code: http://martinos.org/~berkin/Bayesian_CS_Prior.zip [2], and Fast Composite Splitting Algorithm for multicontrast MRI (FCSAMT)^{5}^{5}5FCSAMT code: http://ranger.uta.edu/~huang/codes/Code_MRI_MT.zip [17], which are both publicly available online. We tune the parameters of each methods so they can perform nearly optimally. In particular, FCSAMT seems to be sensitive to noise level and we need to tune a very different parameter for different . Specifically, we set the weight in BCS as for both and cases; the TV term and wavelet term weights in FCSAMT are set to and (resp.) for case, and and (resp.) for case; the weight of ER/ERweighted is set to , and the fidelity weight for (19) to obtain image from edges is set to for both and cases.
The image datasets we used are a 2D multicontrast SheppLogan phantom (size ), a simulated multicontrast brain image (size ) obtained from BrainWeb^{6}^{6}6BrainWeb: http://brainweb.bic.mni.mcgill.ca/brainweb/, and a 2D invivo brain image (size ) included in the FCSAMT code. All images contain three contrasts. In particular, the three contrasts in BrainWeb and invivo brain datasets represent the T1, T2, and PD images. Undersampling pattern and ratio are presented below with results. The experiments are performed in Matlab computing environment on a Mac OS with Intel i7 CPU and 16GB of memory. Gaussian white noise with standard deviation and were added to the kspace for the simulated data. We full Fourier data to obtain reference (ground truth) image . The relative (L2) error of reconstruction for modality is defined by and used to measure error of quantitatively for .
4.1 Comparison of different matrix norms
Our first test is on the performance of three matrix norms, namely Frobenius, Induced 2norm, and Nuclear norm, when used in Algorithm 2. These norms correspond to three commonly used vectorial TV regularization in multimodal/channel/contrast joint image reconstruction. All three norms yield closedform solutions for subproblem (7) as we showed in Section 3.4. We apply the proposed Algorithm 2 with these three norms to different image data and noise level combinations, and observe very similar performance. For demonstration, we show a typical result using BrainWeb data with no noise in left panel of Figure 1. As we can see, all norms yield very similar accuracy in terms of reconstruction error.
In terms of realworld computational cost, however, Frobenius norm treats matrix as vector and hence the shrinkage (15) is very cheap to compute, while both induced 2norm and nuclear norm require (reduced) SVD in (16) and (17) that can be much slower despite that the matrices all have small size ( for 2D images with modalities/contrasts). We show the same trajectory of relative errors but versus CPU time in the right panel of Figure 1, and it appears that Frobenius norm is the most costeffective choice for this test. In the remainder of this section, we only use Frobenius norm (corresponding to the specific vectorial TV norm of form ) for the proposed Algorithm 2, but not the induced 2norm (corresponding to ) and nuclear norm (corresponding to ).
4.2 Comparisons on multicontrast MRI datasets
Now we conduct comparison of the proposed algorithm with BSC and FCSAMT on multicontrast MRI datasets. We first test the comparison algorithms on a multicontrast SheppLogan phantom image. To demonstrate robustness of the proposed algorithm, we use a radial mask of undersampling ratio in Fourier space, but add white complexvalued Gaussian noise with standard deviation (for both real and imaginary parts of the noise) to the undersampled Fourier data. Then we apply all comparison algorithm, namely BCS, FSCAMT, and proposed ER and ERweighted to the undersampled noisy Fourier data to reconstruct the multicontrast image. For each method, we record the progress of relative error and show the relative error vs CPU time in Figure 2. From these plots, we can see that the proposed Algorithm 2 (ER/ERweighted) achieves higher efficiency as they reduce reconstruction error faster than the comparison algorithms BCS and FCSAMT. Moreover, the proposed algorithm with the weight (ERweighted) can incorporate the transformed noise distribution and further improve reconstruction accuracy. We observe that ERweighted sometimes performs slightly slower than ER in terms of CPU time (although always faster in terms of iteration number which is not shown here), since the ERweighted requires an additional pointwise division using . We believe that this is an issue that can be resolved by further optimizing code, because complexity wise the additional computation required by ERweighted is rather low compared to other operations (such as Fourier transforms) in ER/ERweighted.
Besides the plot of relative error vs CPU time, we also show the final reconstruction images using these methods in Figure 3 for the noise level case. In Figure 3, the radial mask of sampling ratio (13.5%) is show on upper left corner (white pixels indicate sampled location). The final reconstructed multicontrast images by the four comparison algorithm are shown in the middle column with algorithm name on left. On the right column, we also show their corresponding error image where is the ground truth multicontrast image obtained using full Fourier data. As we can see, the reconstruction error of proposed algorithm is very small compared other others, especially showing less obvious error on edges. This demonstrates the improve reconstruction accuracy using our method.
We also conduct the same test using the BrainWeb image with a Poisson mask of sampling ratio 25%, also with noise level and . The results are shown in Figure 4. For this “near real” image, our method again shows significant improvement of efficiency compared to other methods. In particular, ERweighted outperforms all methods in both cases, suggesting its superior robustness in multicontrast image reconstruction. The final reconstruction images and error images, similar to those for SheppLogan image, are also plotted in Figure 5 for the case.
Finally, we conduct the same test on the invivo brain image with radial mask of sampling ratio 25%. The relative error vs CPU time plots are shown in Figure 6 ( and ). For this image, our method further shows its promising performance with significant improvement when compared to existing methods. The final reconstruction images and error images are also shown in Figure 7 for the case.
5 Conclusion
We proposed to reconstruct Jacobian of mutlimodal/contrast/channel image from which we can resemble the underlying image. We showed the relation between Jacobian and the observed data when the underlying transform is Fourier, and formulate the reconstruction problem of Jacobian as an minimization. Our new method then exhibits an optimal convergence rate which outperforms the rate of primaldual based algorithms. We also derived closedform solutions for the minimization subproblem as matrixvalued shrinkage. The periteration complexity is thus very low. Numerical results demonstrated the promising performance of the proposed method when compared to the stateofthearts joint image reconstruction methods.
Acknowledgement
The authors would like to thank Yao Xiao and Yun Liang for assisting Fang to run the MATLAB code and collect numerical results. This research is supported in part by National Science Foundation grant DMS1719932 (Chen), IIS1564892 (Fang), DMS1620342 (Ye) and CMMI1745382 (Ye), and National Key Research and Development Program of China No: 2016YFC1300302 (Fang) and National Natural Science Foundation of China No: 61525106 (Fang).
Appendix
A Computation of matrixvalued shrinkage
We show that the minimization problem (14) has closedform solutions when the matrix norm is Frobenius, induced 2norm, or nuclear norm.

Frobenius norm. In the case, all matrices can be considered as vectors and hence the vectorvalued shrinkage formula can be directly applied. We omit the derivations here.

Switching min and max, and solving for with fixed and , we obtain and the dual problem of (20) becomes
(21) where we used the facts that and . Therefore, and are the left and right singular vectors of , and hence the optimal solution of (20) is .
To obtain , one needs to compute the largest singular vectors of , or the SVD of . Note that has size , and hence a reduced SVD for is sufficient. In particular, for , there is a closed form expression to compute SVD of . In general, SVD of a matrix is easy to compute and has at most two singular values.

Nuclear norm This corresponds to standard nuclear norm shrinkage: compute SVD of , and set . Computation complexity is comparable to the 2norm case above.
References
 [1] A. Beck and M. Teboulle. A fast iterative shrinkagethresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1):183–202, 2009.
 [2] B. Bilgic, V. K. Goyal, and E. Adalsteinsson. Multicontrast reconstruction with bayesian compressed sensing. Magnetic Resonance in Medicine, 66(6):1601–1615, 2011.
 [3] P. Blomgren and T. F. Chan. Color tv: total variation methods for restoration of vectorvalued images. IEEE transactions on image processing, 7(3):304–309, 1998.
 [4] X. Bresson and T. F. Chan. Fast dual minimization of the vectorial total variation norm and applications to color image processing. Inverse problems and imaging, 2(4):455–484, 2008.
 [5] S. Di Zenzo. A note on the gradient of a multiimage. Computer vision, graphics, and image processing, 33(1):116–125, 1986.
 [6] J. Duran, M. Moeller, C. Sbert, and D. Cremers. Collaborative total variation: a general framework for vectorial tv models. SIAM Journal on Imaging Sciences, 9(1):116–151, 2016.
 [7] M. Ehrhardt. Joint reconstruction for multimodality imaging with common structure. PhD thesis, UCL (University College London), 2015.
 [8] M. J. Ehrhardt and S. R. Arridge. Vectorvalued image processing by parallel level sets. IEEE Transactions on Image Processing, 23(1):9–18, 2014.
 [9] M. J. Ehrhardt and M. M. Betcke. Multicontrast mri reconstruction with structureguided total variation. SIAM Journal on Imaging Sciences, 9(3):1084–1106, 2016.
 [10] M. J. Ehrhardt, P. Markiewicz, M. Liljeroth, A. Barnes, V. Kolehmainen, J. S. Duncan, L. Pizarro, D. Atkinson, B. F. Hutton, S. Ourselin, et al. Pet reconstruction with an anatomical mri prior using parallel level sets. IEEE transactions on medical imaging, 35(9):2189–2199, 2016.
 [11] M. J. Ehrhardt, K. Thielemans, L. Pizarro, D. Atkinson, S. Ourselin, B. F. Hutton, and S. R. Arridge. Joint reconstruction of petmri by exploiting structural similarity. Inverse Problems, 31(1):015001, 2015.
 [12] M. J. Ehrhardt, K. Thielemans, L. Pizarro, P. Markiewicz, D. Atkinson, S. Ourselin, B. F. Hutton, and S. R. Arridge. Joint reconstruction of petmri by parallel level sets. In Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC), 2014 IEEE, pages 1–6. IEEE, 2014.
 [13] V. Estellers, S. Soatto, and X. Bresson. Adaptive regularization with the structure tensor. IEEE Transactions on Image Processing, 24(6):1777–1790, 2015.
 [14] B. Goldluecke and D. Cremers. An approach to vectorial total variation based on geometric measure theory. In Computer Vision and Pattern Recognition (CVPR), 2010 IEEE Conference on, pages 327–333. IEEE, 2010.
 [15] M. Grasmair and F. Lenzen. Anisotropic total variation filtering. Applied Mathematics & Optimization, 62(3):323–339, 2010.
 [16] E. Haber and M. H. Gazit. Model fusion and joint inversion. Surveys in Geophysics, 34(5):675–695, 2013.
 [17] J. Huang, C. Chen, and L. Axel. Fast multicontrast mri reconstruction. Magnetic resonance imaging, 32(10):1344–1352, 2014.
 [18] F. Lenzen and J. Berger. Solutiondriven adaptive total variation regularization. In International Conference on Scale Space and Variational Methods in Computer Vision, pages 203–215. Springer, 2015.
 [19] A. Majumdar and R. K. Ward. Joint reconstruction of multiecho mr images using correlated sparsity. Magnetic resonance imaging, 29(7):899–906, 2011.
 [20] V. M. Patel, R. Maleh, A. C. Gilbert, and R. Chellappa. Gradientbased image recovery methods from incomplete fourier measurements. IEEE Transactions on Image Processing, 21(1):94–105, 2012.
 [21] E. Sakhaee, M. Arreola, and A. Entezari. Gradientbased sparse approximation for computed tomography. In Biomedical Imaging (ISBI), 2015 IEEE 12th International Symposium on, pages 1608–1611. IEEE, 2015.
 [22] G. Sapiro. Vectorvalued active contours. In Computer Vision and Pattern Recognition, 1996. Proceedings CVPR’96, 1996 IEEE Computer Society Conference on, pages 680–685. IEEE, 1996.
 [23] G. Sapiro and D. L. Ringach. Anisotropic diffusion of multivalued images with applications to color filtering. IEEE transactions on image processing, 5(11):1582–1586, 1996.