Image Restoration by Combined Order Regularization with Optimal Spatial Adaptation

TL;DR

This paper introduces a novel non-convex regularization method for image restoration that adaptively combines Hessian-Schatten norm and TV1 with spatially varying weights, improving performance in applications like MRI and microscopy.

Contribution

It proposes a new adaptive non-convex regularization functional with a joint minimization approach and convergence proof, advancing multi-order image restoration techniques.

Findings

Improved image recovery in MRI and microscopy tasks.

Effective adaptive weighting enhances restoration quality.

Convergence of the proposed optimization algorithm is proven.

Abstract

Total Variation (TV) and related extensions have been popular in image restoration due to their robust performance and wide applicability. While the original formulation is still relevant after two decades of extensive research, its extensions that combine derivatives of first- and second-order are now being explored for better performance, with examples being Combined Order TV (COTV) and Total Generalized Variation (TGV). As an improvement over such multi-order convex formulations, we propose a novel non-convex regularization functional which adaptively combines Hessian-Schatten (HS) norm and first order TV (TV1) functionals with spatially varying weight. This adaptive weight itself is controlled by another regularization term; the total cost becomes the sum of this adaptively weighted HS-TV1 term, the regularization term for the adaptive weight, and the data-fitting term. The reconstruction is obtained by jointly minimizing w.r.t. the required image and the adaptive weight. We construct a block coordinate descent method for this minimization with proof of convergence, which alternates between minimization w.r.t. the required image and the adaptive weights. We derive exact computational formula for minimization w.r.t. the adaptive weight, and construct an ADMM algorithm for minimization w.r.t. to the required image. We compare the proposed method using image recovery examples including MRI reconstruction and microscopy deconvolution.

Tables

Table 1. Table 1: Comparison of deconvolution results for TIRF images

IMG	${\gamma }_p$	SSIM						SNR
		TV2	HS	COTV	COHS	COR-I	COR.	TV2	HS	COTV	COHS	COR-I	COR.
Act.1	5	.725	.731	.502	.502	.728	.762	10.53	10.58	9.37	9.37	10.58	10.97
	10	.794	.793	.679	.679	.792	.807	13.44	13.45	12.88	12.89	13.40	13.82
	20	.853	.849	.806	.807	.848	.853	15.62	15.64	15.38	15.43	15.59	15.90
	30	.890	.888	.867	.867	.887	.890	16.88	16.91	16.81	16.81	16.81	17.08
Act.2	5	.735	.744	.516	.516	.741	.780	12.03	12.06	10.74	10.74	12.05	12.34
	10	.828	.829	.712	.712	.829	.844	14.59	14.58	13.96	13.97	14.58	14.74
	20	.883	.881	.839	.839	.880	.882	16.58	16.56	16.38	16.39	16.54	16.55
	30	.904	.900	.880	.880	.897	.899	17.69	17.66	17.58	17.59	17.60	17.57
Tub.1	5	.769	.784	.491	.491	.780	.827	13.42	13.60	11.30	11.30	13.57	13.92
	10	.830	.832	.685	.685	.830	.844	15.70	15.71	14.76	14.76	15.69	15.86
	20	.868	.864	.806	.806	.863	.862	17.46	17.52	17.17	17.18	17.50	17.55
	30	.881	.877	.850	.850	.875	.875	18.27	18.29	18.07	18.11	18.23	18.25

Table 2. Table 2: Comparison of MRI reconstruction results

IMG	Traj.	SR	I/P PSNR	SSIM							SNR
				TV2	HS	TGV	COTV	COHS	COR-I	COR.	TV2	HS	TGV	COTV	COHS	COR-I	COR.
MRI_1	Spiral	.2	10	.972	.976	.979	.983	.982	.976	.983	24.71	25.78	26.68	28.26	28.16	25.72	28.07
			20	.972	.976	.979	.983	.982	.976	.983	24.71	25.78	26.68	28.27	28.19	25.73	28.12
		.1	20	.955	.955	.959	.964	.964	.956	.966	22.07	22.14	22.86	24.20	24.12	22.35	24.09
	Random	.2	10	.902	.927	.927	.926	.934	.946	.962	19.39	20.74	19.72	22.80	22.98	21.84	24.35
			20	.902	.927	.927	.926	.934	.946	.962	19.39	20.73	19.75	22.81	22.99	21.85	24.35
		.1	20	.807	.811	.793	.810	.810	.839	.905	14.42	14.73	14.38	16.38	16.38	15.71	18.52
MRI_2	Spiral	.2	10	.990	.993	.995	.998	.998	.994	.998	22.68	24.53	26.06	30.66	30.66	24.70	30.28
			20	.990	.993	.995	.998	.998	.994	.998	22.68	24.54	26.06	30.77	30.77	24.71	30.45
		.1	20	.978	.978	.982	.993	.993	.982	.992	18.95	19.09	20.20	24.40	24.40	19.93	23.58
	Random	.2	10	.936	.963	.965	.992	.992	.983	.995	16.69	18.87	17.69	27.49	26.60	20.76	27.54
			20	.936	.963	.957	.992	.992	.983	.995	16.69	18.88	16.65	26.65	26.65	20.77	27.60
		.1	20	.853	.856	.846	.887	.904	.895	.959	11.96	12.14	11.77	15.26	15.48	13.44	17.38
MRI_3	Spiral	.2	10	.940	.948	.952	.961	.961	.948	.962	21.50	22.61	23.12	24.55	24.55	22.67	24.56
			20	.940	.948	.952	.961	.961	.949	.962	21.50	22.62	23.12	24.57	24.57	22.67	24.58
		.1	20	.879	.882	.879	.891	.895	.880	.898	18.00	18.18	18.33	19.58	19.58	18.33	19.44
	Random	.2	10	.750	.784	.760	.790	.790	.804	.864	14.88	16.04	15.45	17.32	17.22	16.72	18.33
			20	.743	.780	.753	.783	.783	.796	.854	14.81	15.92	15.31	16.96	16.93	16.49	17.98
		.1	20	.498	.503	.486	.547	.562	.547	.651	8.64	8.79	8.73	9.71	11.02	9.43	11.51
MRI_4	Spiral	.2	10	.992	.994	.990	.996	.996	.994	.996	28.60	29.92	28.13	33.22	32.66	30.12	32.65
			20	.992	.994	.990	.996	.996	.994	.996	28.60	29.94	28.13	33.28	32.72	30.14	32.73
		.1	20	.983	.982	.984	.989	.988	.984	.989	24.69	24.60	25.60	28.01	27.56	25.42	27.21
	Random	.2	10	.971	.977	.978	.978	.979	.980	.985	24.58	26.19	26.22	30.66	30.15	27.38	30.06
			20	.970	.976	.978	.978	.980	.980	.985	24.49	26.10	26.13	30.84	30.20	27.30	30.13
		.1	20	.928	.927	.920	.936	.936	.942	.961	19.48	19.48	19.37	21.72	21.72	20.54	22.54
MRI_5	Spiral	.2	10	.958	.964	.968	.973	.973	.963	.974	21.00	21.98	22.75	24.28	24.28	21.95	24.35
			20	.958	.964	.969	.974	.974	.963	.975	21.00	21.98	22.82	24.31	24.31	21.96	24.38
		.1	20	.922	.923	.926	.936	.936	.927	.938	18.16	18.19	18.50	19.63	19.63	18.59	19.48
	Random	.2	10	.856	.886	.888	.896	.902	.914	.951	17.08	18.11	17.47	19.89	19.89	18.86	21.60
			20	.851	.880	.880	.890	.894	.912	.954	16.97	18.02	17.43	19.88	19.88	18.89	21.91
		.1	20	.737	.736	.732	.763	.759	.779	.871	13.84	13.88	13.43	15.09	14.83	14.59	16.16
MRI_6	Spiral	.2	10	.983	.987	.984	.994	.994	.987	.995	19.36	20.98	19.96	25.89	25.89	21.14	27.25
			20	.984	.987	.984	.994	.994	.987	.996	19.36	20.98	19.96	25.96	25.96	21.14	27.43
		.1	20	.972	.972	.977	.984	.984	.975	.986	17.15	17.18	18.62	20.96	20.96	17.76	21.21
	Random	.2	10	.956	.969	.967	.978	.978	.976	.989	16.61	18.26	17.00	22.09	22.09	18.95	23.84
			20	.952	.965	.968	.976	.976	.976	.990	16.78	18.54	17.22	22.42	22.42	19.18	24.35
		.1	20	.896	.896	.892	.911	.915	.925	.954	12.82	12.89	13.09	15.30	15.30	13.83	16.67

Table 3. Table 3: Comparison of MRI reconstruction results with DAGAN and DL-MRI (Noiseless Measurements)

IMG	SR	SNR					SSIM
		DL- MRI	DAGAN			COR.	DL- MRI	DAGAN			COR.
			10%	50%	100%			10%	50%	100%
I1	.1	18.29	27.71	28.85	28.52	28.30	.722	.953	.964	.961	.991
	.2	22.73	31.23	32.35	31.70	37.78	.949	.971	.978	.975	.999
I2	.1	18.13	27.25	28.18	28.22	27.44	.735	.949	.958	.960	.989
	.2	22.20	30.66	32.04	31.28	36.68	.944	.967	.976	.972	.998
I3	.1	18.40	26.88	27.70	27.81	24.98	.775	.940	.950	.950	.978
	.2	22.09	29.87	30.67	30.59	34.62	.934	.963	.969	.969	.997
I4	.1	17.21	25.96	26.95	27.15	23.72	.770	.933	.944	.946	.973
	.2	21.93	28.56	29.29	29.27	32.92	.935	.956	.964	.962	.996
I5	.1	16.89	24.67	25.88	25.67	22.77	.750	.918	.931	.932	.963
	.2	21.90	28.56	29.45	29.14	31.98	.928	.954	.961	.960	.994
I6	.1	16.92	24.32	25.50	25.40	22.37	.751	.912	.925	.927	.958
	.2	21.58	27.91	28.65	28.46	31.09	.924	.948	.956	.955	.993

Table 4. Table 4: Comparison of MRI reconstruction results with DAGAN and DL-MRI (Noisy Measurements)

IMG	SR	SNR					SSIM
		DL- MRI	DAGAN			COR.	DL- MRI	DAGAN			COR.
			10%	50%	100%			10%	50%	100%
I1	.1	16.16	19.61	19.42	19.50	19.64	.349	.738	.837	.820	.922
	.2	19.91	24.47	24.67	25.21	22.72	.475	.845	.890	.888	.956
I2	.1	16.43	18.67	18.80	19.07	19.40	.363	.718	.819	.806	.911
	.2	19.80	24.22	24.06	24.89	22.60	.481	.847	.883	.890	.952
I3	.1	16.77	18.52	17.91	18.37	18.60	.393	.717	.801	.793	.905
	.2	19.48	22.13	20.92	21.61	22.16	.452	.795	.784	.769	.938
I4	.1	16.08	17.48	16.75	17.52	18.19	.398	.702	.783	.777	.902
	.2	19.38	20.69	19.39	19.68	21.56	.448	.746	.715	.698	.932
I5	.1	15.41	16.81	16.25	16.75	17.70	.394	.696	.769	.766	.889
	.2	19.36	20.90	19.32	19.42	21.20	.457	.757	.706	.696	.925
I6	.1	14.87	16.84	16.09	16.53	17.34	.391	.686	.757	.749	.880
	.2	19.52	20.54	19.28	19.26	21.17	.473	.762	.729	.708	.924

Equations

$$\hat{s}(\mathbf{r})=\arg \underset{s}{\min }F(s,h,m)+\lambda R(s),$$

$$F(s,h,m)=\sum _{\mathbf{r}}|(h\ast s)(\mathbf{r})-m(\mathbf{r}){|}^2$$

$$m(\mathbf{r})=\mathcal{P}((h\ast s)(\mathbf{r}))+\eta (\mathbf{r}),$$

$$\hat{\mathbf{m}}={\mathcal{T}}_{s}M\mathcal{F}s(\mathbf{r})+\hat{\mathbf{\eta }},$$

$$\begin{array}{c}{s}_{opt}=\underset{s}{\mathrm{argmin}}F(s,h,m)+\lambda \underset{R_1(s)}{\underbrace{\sum _{\mathbf{r}}\Vert ({\mathbf{d}}_1\ast s)(\mathbf{r}){\Vert }_2}}\\ \text{mbox}with{d}_1(\mathbf{r})=[d_x(\mathbf{r}),d_y(\mathbf{r}){]}^T,\end{array}$$

$$\begin{array}{c}{s}_{opt}=\underset{s}{\mathrm{argmin}}F(s,h,m)+\lambda \underset{R_2(s)}{\underbrace{\sum _{\mathbf{r}}\Vert ({\mathbf{d}}_2\ast s)(\mathbf{r}){\Vert }_2}},\\ \text{mbox}with{\mathbf{d}}_2(\mathbf{r})=[d_{xx}(\mathbf{r})d_{yy}(\mathbf{r})\sqrt{2}{d}_{xy}(\mathbf{r}){]}^T,\end{array}$$

$$s_{opt}=\underset{s}{\mathrm{argmin}}F(s,h,m)+\lambda \sum _{\mathbf{r}}{\Vert \mathbf{\zeta }(({\mathbf{H}}_2\ast s)(\mathbf{r}))\Vert }_p.$$

$$s_{opt}=\underset{s}{\mathrm{argmin}}F(s,h,m)+\lambda {\alpha }_1{R}_1(s)+\lambda {\alpha }_2{R}_2(s),$$

$$\begin{array}{rl}(s_{opt},{\mathbf{p}}_{opt})& =\underset{s,\mathbf{p}}{\mathrm{argmin}}F(s,h,m)\\ & +\lambda {\alpha }_1\sum _{\mathbf{r}}\Vert ({\mathbf{d}}_1\ast s)(\mathbf{r})-\mathbf{p}(\mathbf{r}){\Vert }_2\\ & +\lambda {\alpha }_2\frac{1}{2}\sum _{\mathbf{r}}\Vert {\mathbf{d}}_1(\mathbf{r})\ast {\mathbf{p}}^T(\mathbf{r})+\mathbf{p}(\mathbf{r})\ast {\mathbf{d}}_1^T(\mathbf{r}){\Vert }_F,\end{array}$$

$$\begin{array}{rl}(s_{opt},{\beta }_{opt})& =\underset{s,\beta }{\mathrm{argmin}}F(s,h,m)+\lambda R_{sa}(s,\beta ,p)+L(\beta ,\tau )+\mathcal{B}(s),\end{array}$$

$$\begin{array}{rl}{R}_{sa}(s,\beta ,p)& =\sum _{\mathbf{r}}\beta (\mathbf{r})\Vert ({\mathbf{d}}_1\ast s)(\mathbf{r}){\Vert }_2+\sum _{\mathbf{r}}(1-\beta (\mathbf{r})){\Vert \mathbf{\zeta }(({\mathbf{H}}_2\ast s)(\mathbf{r}))\Vert }_p,\\ & \text{mbox}subjectto0\le \beta (\mathbf{r})\le 1,\end{array}$$

$$L(\beta ,\tau )=-\sum _{\mathbf{r}}\tau \log (\beta (\mathbf{r})(1-\beta (\mathbf{r}))).$$

$$J_{sa}(s,\beta ,\tau ,h,m)=F(s,h,m)+\lambda R_{sa}(s,\beta ,p)+L(\beta ,\tau )+\mathcal{B}(s).$$

$$F^{(j)}(s,h,m)=\sum _{\mathbf{r}}|(h\ast (E^{(j)}s))(\mathbf{r})-m(\mathbf{r}){|}^2.$$

$$\begin{array}{rl}{R}_{sa}^{(j)}(s,\beta ,p)& =\sum _{\mathbf{r}}\beta (\mathbf{r})\Vert ({\mathbf{d}}_1\ast (E^{(j)}s))(\mathbf{r}){\Vert }_2\sum _{\mathbf{r}}(1-\beta (\mathbf{r})){\Vert \mathbf{\zeta }(({\mathbf{H}}_2\ast (E^{(j)}s))(\mathbf{r}))\Vert }_p.\end{array}$$

$$\begin{array}{rl}{J}_{sa}^{(j)}(s,\beta ,\tau ,h,m)& =F^{(j)}(s,h,m)+\lambda R_{sa}^{(j)}(s,\beta ,p)+L(\beta ,\tau )+\mathcal{B}(E^{(j)}s).\end{array}$$

$$\begin{array}{rl}{\hat{s}}^{(K)}(\mathbf{r})& =\underset{s}{\mathrm{argmin}}{J}_{sa}^{(K)}(s,\beta ,\tau ,h,m)\\ & \equiv \underset{s}{\mathrm{argmin}}{F}^{(K)}(s,h,m)+\lambda R_{sa}^{(K)}(s,\beta ,p)+\mathcal{B}(E^{(K)}s).\end{array}$$

$$f$$

$$\bar{\beta }(\mathbf{r})$$

$$\equiv \underset{\beta }{\mathrm{argmin}}{R}_{sa}(f,\beta ,p)+L(\beta ,\tau ),$$

$$\text{mbox}subjectto0\le \beta (\mathbf{r})\le 1.$$

$${\hat{s}}^{(j)}(\mathbf{r})$$

$$\equiv \underset{s}{\mathrm{argmin}}{F}^{(j)}(s,h,m)+\lambda R_{sa}^{(j)}(s,\bar{\beta },p)+\mathcal{B}(E^{(j)}s).$$

$$d(\mathbf{r})=\Vert ({\mathbf{d}}_1\ast f)(\mathbf{r}){\Vert }_2-{\Vert \mathbf{\zeta }(({\mathbf{H}}_2\ast f)(\mathbf{r}))\Vert }_p.$$

$$\bar{\beta }(\mathbf{r})=\frac{1}{2}\left(1-sign(d(\mathbf{r}))\left(\sqrt{\frac{4{\tau }^2}{d^2(\mathbf{r})}+1}-\frac{2\tau }{|d(\mathbf{r})|}\right)\right).$$

$${\beta }_{(k+1)}$$

$$=\underset{\beta }{\mathrm{argmin}}{R}_{sa}(s_{(k)},\beta ,p)+L(\beta ,\tau ).$$

$$s_{(k+1)}$$

$$=\underset{s}{\mathrm{argmin}}F(s,m,h)+\lambda R_{sa}(s,{\beta }_{(k+1)},p)+\mathcal{B}(s).$$

$$F^{(j)}(\mathbf{s},\mathbf{H},\mathbf{m})={\Vert \mathbf{H}{\mathbf{E}}^{(\mathbf{j})}\mathbf{s}-\mathbf{m}\Vert }_2^2.$$

$$\Vert ({\mathbf{d}}_1\ast (E^{(j)}s))(\mathbf{r}){\Vert }_2={\Vert {\mathbf{P}}_{i(\mathbf{r})}{\mathbf{D}}_f{\mathbf{E}}^{(j)}\mathbf{s}\Vert }_2$$

$${\Vert \mathbf{\zeta }(({\mathbf{H}}_2\ast (E^{(j)}s))(\mathbf{r}))\Vert }_p={\Vert \mathbf{\zeta }(\mathcal{S}({\mathbf{Q}}_{i(\mathbf{r})}{\mathbf{D}}_s{\mathbf{E}}^{(j)}\mathbf{s}))\Vert }_p$$

$$R_{sa}^{(j)}(\mathbf{s},\mathbf{\beta },p)=\sum _{i=1}^{N^2}{\bar{\beta }}_i{\Vert {\mathbf{P}}_i{\mathbf{D}}_f{\mathbf{E}}^{(j)}\mathbf{s}\Vert }_2+\sum _{i=1}^{N^2}(1-{\bar{\beta }}_i){\Vert \mathbf{\zeta }(\mathcal{S}({\mathbf{Q}}_i{\mathbf{D}}_s{\mathbf{E}}^{(j)}\mathbf{s}))\Vert }_p,$$

$${\mathcal{N}}_f(\mathbf{u},\mathbf{v})$$

Full text

K****eywords Total Variation, Image Restoration, Multi-Order Regularization, Hessian-Schatten norm, Spatially Adaptive Regularization, Magnetic Resonance Imaging, Total Internal Reflection Fluorescence Microscopy.

00footnotetext: © 2020 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.

1 Introduction

Regularization plays an important role in image reconstruction/restoration by stabilizing the inversion of the imaging forward model against noise and other distortions in modalities such as photography [1], microscopy [2], astronomical imaging [3] and magnetic resonance imaging [4]. Here we consider the measurement model in which the measured image is expressed as a convolution of the underlying image with a blurring function. In regularized reconstruction, the required image is the solution of the following minimization problem:

[TABLE]

where $R(s)$ is the regularization functional, and $F(s,h,m)$ is the data fidelity term with $m$ being the measured image and $h$ being the blurring kernel of the imaging system. We restrict the data fidelity term to be of the form,

[TABLE]

Strictly speaking, the data-fidelity term should be the negative log-likelihood of the noise model. However, our focus is on developing improved regularization, and hence we use the above form of data-fidelity as an approximation, irrespective of the noise model. In our evaluation of the regularized image reconstruction methods, we consider two types of imaging forward models with respect to our restoration experiments. The first measurement model is related to total internal reflection fluorescence (TIRF) microscopy and is given by

[TABLE]

where $h(\mathbf{r})$ is the PSF of TIRF microscope, ${\mathcal{P}}(\cdot)$ represents the Poisson process, and $\eta(\mathbf{r})$ is additive white Gaussian noise. The second measurement model involved in our experiments corresponds to magnetic resonance imaging (MRI) and is given by

[TABLE]

where ${\cal F}$ represents Fourier transformation operation, $M$ represents the under-sampling mask of ones and zeros, ${\cal T}_{s}$ represents the operations of picking samples from non-zero locations of $M$, and $\hat{\boldsymbol{\eta}}$ is the complex noise vector where each component comes from a Gaussian distribution. For the first case, equation (2) is not the exact negative log-likelihood, and for the second case, it is the exact negative log-likelihood. For the MRI model in the second case, $h$ is the inverse Fourier transform of $M$ and $m(\mathbf{r})={\cal F}^{-1}({\cal T}^{\dagger}_{s}\hat{\bf m})$, with ${\cal T}^{\dagger}_{s}$ denoting the adjoint of ${\cal T}_{s}$, which is essentially the operation of embedding the Fourier samples into image form. Note that, in this case, $m(\mathbf{r})$ can have complex values and it is taken care by the operator $|\cdot|$ in equation (2).

As mentioned before, our focus is on developing an improved regularization method. The quality of restoration is mainly determined by the ability of $R(s)$ to discriminate between characteristics of the underlying image and noise. While priors can be defined based on general image characteristics such as sparsity of the roughness [5, 6, 7, 8], there are also priors which are tailored to specific classes of images [9, 10, 11]. The latter type of regularization approach utilizes learning paradigms including deep learning based techniques [12, 13, 14, 15], dictionary learning [9, 10, 16, 17] and model fitting [11]. Here the functionals are built from training images and then applied for restoration involving images from the same class. Such methods are believed to out-perform general priors, when suitable training sets are available. At the same time, the necessity of training set and computational complexity limits their applicability.

On the other hand, general priors such as Tikhonov regularization [18] and total variation (TV) [19] do not need training samples and have been applied with robust performance in multiple domains. Among these general priors, total variation (TV) [19] has been widely applied [5, 6, 7, 8] because of its ability to recover sharp image features in the presence of noise and in the cases of undersampling. First order TV (TV1) restoration is given by

[TABLE]

where $d_{x}(\mathbf{r})$ and $d_{y}(\mathbf{r})$ are filters implementing first order derivatives, $\frac{\partial}{\partial x}$, and $\frac{\partial}{\partial y}$ respectively. While TV1 is able to retain edges [20] in the reconstruction as compared to standard $\ell_{2}$ norm based Tikhonov regularization [18], it presents drawbacks such as staircase artifacts [21, 22]. Higher order extensions of TV [23, 24, 25, 26] have been proposed to avoid staircase artifacts and they deliver better restoration, albeit at the cost of increased computations. Second order TV (TV2) [24] restoration was proposed as

[TABLE]

where $d_{xx}(\mathbf{r}),d_{yy}(\mathbf{r})$, and $d_{xy}(\mathbf{r})$ are discrete filters implementing second order derivatives $\frac{\partial^{2}}{\partial x^{2}}$, $\frac{\partial^{2}}{\partial y^{2}}$ and $\frac{\partial^{2}}{\partial x\partial y}$ respectively. Another second-order derivative based formulation is Hessian-Schatten (HS) norm regularization [27] [26], which has been proposed as a generalization of the standard TV2 regularization. It is constructed as an $\ell_{p}$ norm of Eigen values of the Hessian matrix, which becomes the standard TV2 for $p=2$. HS norm with $p=1$ has been proven to yield best resolution in the reconstruction, since this better preserves Eigen values of the Hessian [26]. Let ${\bf H}_{2}({\bf r})$ be the matrix filter composed of $d_{xx}({\bf r})$, $d_{yy}({\bf r})$, and $d_{xy}({\bf r})$ and let ${\boldsymbol{\zeta}}(\cdot)$ be the operator that returns the vector containing the Eigen values of its matrix argument. Then HS norm regularization of order $p$ can be expressed as

[TABLE]

Since the Eigen values are actually directional second derivatives taken along principle directions, setting $p=1$ better preserves the local image structure. It has to be noted that the costs given in the equations (5), and (6) are often minimized using gradient based approaches with smooth approximations of the form $R_{k}(s)={\sum\limits_{\mathbf{r}}}\sqrt{\epsilon+||d_{k}(\mathbf{r})\ast s(\mathbf{r})||_{2}^{2}},\ k=1,2$ where $\epsilon$ is a small positive constant [19], [28]. This approach has been proven to converge to the minimum of the exact form as $\epsilon\rightarrow 0$ [28]. Approaches to minimize the cost without smooth approximation include primal-dual method [29], and alternating direction of multiplier method (ADMM) [30]. A detailed comparison of such approaches has been provided in [31].

It has been demonstrated that combining first- and second-order derivatives is advantageous in accurately restoring image features [32, 33, 34, 35]. In this regard, the combined order TV [33] uses scalar relative weights for combining first- and second-order variations, with the relative weights left as user parameters and the solution is estimated by means of optimization problem of the form

[TABLE]

where $\alpha_{1}$ and $\alpha_{2}$ determine the relative weights. Although HS norm can be in principle combined with TV1 by the same way as standard TV2 is combined as given above, this possibility has not been explored. A generalization for total variation to higher order terms, named as total generalized variation (TGV) has also been proposed [34] [36]. It is generalized in the following ways: it is formulated for any general derivative order, and for any given order, it is generalized in the way how the derivatives are penalized. In the literature, only the second order TGV form has been well explored for image reconstruction, which takes the following form:

[TABLE]

where ${\bf p}({\bf r})$ is an auxiliary $2\times 1$ vector image. The TGV functional is able to spatially adapt to the underlying image structure because of the minimization w.r.t. auxiliary variable ${\bf p}$. Near edges, ${\bf p}({\bf r})$ approaches zero leading to TV1-like behavior which allows sharp jumps in the edges. On the other hand, in smooth regions, ${\bf p}({\bf r})$ approaches ${\bf d}_{1}*s({\bf r})$ leading to TV2-like behavior which will avoid staircase artifacts. However, the drawback with TGV functional is that the weights $\alpha_{1}$ and $\alpha_{2}$ have to be chosen by the user. In summary, existing combined order methods can bring the advantages of multi-order derivatives albeit at the cost of additional tuning.

We propose a novel spatially adaptive regularization method, in which, the weights involved in combining first- and second-order derivatives are determined from the measured image without user-intervention. Here the relative weight between first- and second-order terms becomes an image, and this weight is determined without user-intervention through minimization of a composite cost function. Our contributions can be summarized as follows:

• •

We construct a composite regularization functional containing two parts: (i) the first part is constructed as the sum of first- and second-order derivative magnitudes with spatially varying relative weights; (ii) the second part is an additional regularization term for preventing rapid spurious variations in the relative weights. For the first order term, we use the norm of the gradient, and for the second order term, we use the Schatten norm of the Hessian [27] [26]. The composite cost functional is convex with respect to either the required image or the relative weight, but it is non-convex jointly.

• •

We construct a block coordinate descent method involving minimizations w.r.t. the required image and the relative weight alternatively with the following structure: the minimization w.r.t. the required image is carried out using ADMM approach [37, 38] and the minimization w.r.t. the relative weight is carried out as a single step exact minimization using a formula that we derive in this paper.

• •

Since the total cost is non-convex, the reconstruction results are highly dependent on the initialization for block-coordinate descent method. We handle this problem using a multi-resolution approach, where, a series of coarse-to-fine reconstructions are performed by minimization of cost functionals defined through upsampling operators. Here, minimization w.r.t. the relative weight and the required image is carried out alternatively, as we progress from coarse to final resolution levels. At the final resolution level, the above-mentioned block coordinate descent method is applied.

• •

Note that the sub-problem of minimization w.r.t. to the required image involves spatially varying relative weights. Further, this sub-minimization problem in the above-mentioned multi-resolution loop involves upsampling operators. Hence, the ADMM method proposed by Papafitsoros et al. [33] turns out to be unsuitable. We propose improved variable splitting method and computational formulas to handle this issue.

• •

We prove that the overall block coordinate descent method converges to a local minimum of the total cost function by using Zangwill’s convergence theorem.

This work is an extension of the work presented in the conference paper [35], where we only considered a differentiable approximation of the total variation functionals, and did not incorporate the joint optimization of the relative weight and the required image. In other words, only the multi-resolution loop is applied without the block-coordinate descent method in the final resolution. This method was named Spatially Adaptive Multi-order TV (SAM-TV). Here we name the improved approach more appropriately as Combined Order Regularization with Optimal Spatial Adaptation (COROSA). The rest of the paper is organized as follows: section 2 deals with the formulation of COROSA functional and the multi-resolution framework, and section 3 presents the ADMM formulation and iterations for the optimization problems associated with COROSA. Section 4 presents the simulation results and comparisons.

2 COROSA Image Restoration

2.1 COROSA formulation

In the proposed COROSA approach, the restoration problem is formulated as given below:

[TABLE]

where

[TABLE]

$L(\beta,\tau)$ is a regularization term for $\beta$, which will be specified soon, and ${\cal B}(\cdot)$ is the indicator function for constraining the restored image to a particular range of positive values. From this formulation, it is clear that the relative weight image $\beta$ is also considered as a minimization variable and the optimal image of weights is determined jointly with the required image. Since $R_{sa}(s,\beta,p)$ is linear in $\beta$, minimizing with respect to $\beta$ means that it will essentially act as a switching between first- and second-order terms. In this context, the role of $L(\beta,\tau)$ is to prevent spurious switching; in other words, its role is to prevent rapid switching between the first- and second-order terms caused by insignificant differences in their magnitudes. We set $L(\beta,\tau)$ as

[TABLE]

Here, a lower value of $\tau$ will cause a more rapid switching between first- and second-order terms and vice versa. We denote the overall cost by $J_{sa}(s,\beta,\tau,h,m)$, and we write

[TABLE]

An assumption that is implicitly made by most of the image restoration algorithms is that there is no $s$ for which $F(s,h,m)$, $R_{1}(s)=\sum_{\bf r}\|({\bf d}_{1}*s)({\bf r})\|_{2}$, and $R_{2}(s)=\sum_{\bf r}\|({\bf H}_{2}*s)({\bf r})\|_{2}$ will have zero value simultaneously. We will also use this assumption for proving the convergence of the iterative method that we propose in the following sections.

2.2 Multiresolution method

The regularization functional, $R_{sa}(s,\beta,p)$, is non-convex jointly with respect to $\beta$ and $s$, although it is convex with any one of them alone. Hence, the reconstruction result becomes sensitive to initialization and finding an efficient initialization becomes crucial. To this end, we adopt multiresolution approach for the initialization. To describe the multi-resolution approach, we define the following:

[TABLE]

[TABLE]

In the above, $E^{(j)}s$ denotes the image obtained by interpolating $s$ by a factor $2^{j}$ along both axes. We will defer the description of the implementation of $E^{(j)}$ to the end. Here, $(h*(E^{(j)}s))({\bf r})$ denotes convolving the interpolated image $E^{(j)}s$ with $h$ followed by accessing the pixel at position $\mathbf{r}$. Further, $({\bf d}_{1}*(E^{(j)}s))({\bf r})$ and $({\bf H}_{2}*(E^{(j)}s))({\bf r}))$ have similar interpretation except that the first expression will be a vector, and the second one will be a matrix. Note that the variable $s$ is an $\frac{N}{2^{j}}\times\frac{N}{2^{j}}$ image in both the functionals $F^{(j)}(s,h,m)$, and $R^{(j)}_{sa}(s,\beta,\tau,p)$. On the other hand, $\beta$ always has size ${N}\times{N}$, which is the size of the measurement $m(\mathbf{r})$. We denote the resulting scale-$j$ cost by

[TABLE]

It should be noted that, size of the variable in scale-$j$ cost is $\frac{N}{2^{j}}\times\frac{N}{2^{j}}$; however, the cost is always evaluated on a ${N}\times{N}$ grid through interpolation by $E^{(j)}$. This will help to ensure a better convergence in the multi-resolution method to be described below.

Let $K$ denote the user-defined number of multi-resolution levels. To initialize the multi-resolution loop, we set $\beta({\bf r})=0$ and perform the following minimization:

[TABLE]

With $\hat{s}^{(K)}$ as the initialization, we iterate for $j=K-1,\ldots,0$ with the following minimizations:

[TABLE]

The resulting restored image at the end of the multi-resolution loop, $\hat{s}^{(0)}$, can be an initialization for the joint minimization problem given in equation (10). Note that, the minimization problem of equation (19) has to be done iteratively and requires an initialization. We use $E^{(1)}\hat{s}^{(j+1)}({\bf r})$ as the initialization. In this regards, the way the multi-scale costs $\{J_{sa}^{(j)}(s,\beta,\tau,h,m),j=0,\ldots,K\}$ are constructed significantly helps to make the initialization $E^{(1)}\hat{s}^{(j+1)}({\bf r})$ to be very close the the required minimum $\hat{s}^{(j)}({\bf r})$. In other words, since each $J_{sa}^{(j)}$ is constructed through $2^{j}$-fold interpolation on the original reconstruction grid, the initialization $E^{(1)}\hat{s}^{(j+1)}({\bf r})$ is typically very close the the required minimum $\hat{s}^{(j)}({\bf r})$.

In the multi-resolution method described above, if we set $p=2$, the result $\hat{s}^{(0)}$ will be equivalent to final reconstruction of SAM-TV approach [35] except the fact that SAM-TV uses smooth approximations for the first- and second-order TV terms, whereas, here, the exact non-differentiable form of TV functionals are used. Now we consider solving the sub-problem of determining the adaptive weight $\bar{\beta}$ in equation (18). The exact solution for $\bar{\beta}$ is given in the following Proposition.

Proposition 1**.**

Let $d(\mathbf{r})$ be defined as

[TABLE]

If $d(\mathbf{r})=0$, the solution for the problem in equation (18) is $\bar{\beta}(\mathbf{r})=0.5$. When $d(\mathbf{r})$ is non-zero, the solution $\bar{\beta}$ is unique and given by

[TABLE]

The proof of Proposition 1 is given in Appendix. Next, unlike equation (18), the subproblem of equation (19) cannot be solved exactly and has to be solved iteratively. We will develop an ADMM based method to solve this problem in Section 3.

Now it remains to specify the implementation of $E^{(j)}$. It can be implemented by $j$ stages of $2$-fold interpolation. The $2$-fold interpolation can be implemented by inserting a zero next to each pixel along both axes (which is called as the two-fold expansion) and then filtering by an appropriate interpolation filter. In our implementation we use $u(\mathbf{r})=\frac{1}{64}[1\;4\;6\;4\;1]^{T}[1\;4\;6\;4\;1]$ as the interpolation filter, which is the two-scale filter of cubic B-spline [39].

2.3 Obtaining the final restoration by block coordinate

descent method

By using the result of the above multiresolution method as the initialization, the final reconstruction has to be obtained by minimizing the cost of equation (10) jointly with respect to $\beta$ and $s$. We propose to use a simple block coordinate descent method. Let $s_{(0)}=\hat{s}^{(0)}$, where $\hat{s}^{(0)}$ is the result of the multi-resolution loop described before. With $k=0,\ldots,N_{b}$, the block coordinate descent method involves the following series of minimizations with respect to $\beta$ and $s$. Let $s_{(k)}$ and $\beta_{(k)}$ be the current estimate of the minimum at cycle $k$. Then the next refined estimate is computed as the following set of minimizations:

[TABLE]

[TABLE]

As evident, the iterations given above are similar to the iterations given in the multi-resolution method of section 2.2. The difference is that the minimization with respect to $\beta$ and $s$ for each of the cost functions in the series $\{J_{sa}^{(j)}(s,\beta,\tau,h,m),j=0,\ldots,K\}$ is done only once in the multiresolution method. On the other hand, the minimizations in the block coordinate decent method (BCD) represented by the equation (22) and (23) are done alternatively on the same cost function $J_{sa}(s,\beta,\tau,h,m)$ until convergence. The functional $J_{sa}(s,\beta,\tau,h,m)$ is convex with respect to either of $s$ and $\beta$, and the BCD method represented by equations (22) and (23) converges to the solution of the problem given in (10), provided that each of the minimizations is exact as per the convergence theorem of Bertsekas [40].

Now we consider solving the sub-problems. The subproblems of determining the adaptive weight $\bar{\beta}$ in the block coordinate descent method of equation (22) is identical to the sub-problem of the multi-resolution method (equation (18)), and hence can be solved exactly. On the other hand, the sub-problem of equation (23) is similar to the problem of the equation (19), and cannot be solved exactly. Hence the convergence result of Bertsekas [40] will not be applicable. However, it is easy to show that BCD iteration converges to the minimum if $J_{sa}(s_{(k+1)},\beta_{(k+1)},\tau,h,m)<J_{sa}(s_{(k)},\beta_{(k+1)},\tau,h,m)$ using Zangwill’s global convergence theorem. We will provide the convergence statement along with the proof after describing the ADMM method for solving the problems of equations (19) and (23) in the next section.

3 Image recovery with fixed relative weight for first- and

second-order derivatives

The main computational task in the block coordinate descent iteration represented by the equations (22) and (23), is the computation of $s_{(k+1)}$. Similarly, in the multi-resolution method represented by equations (18) and (19), the main task is the computation of $\hat{s}^{(j)}({\bf r})$. Note that the cost in (19) becomes algebraically identical to the cost of (23) for $j=0$. Hence the cost in equation (23) can be considered as a special case of the cost in the equation (19). So we consider only the description of the minimization of the cost in (19). We will use the ADMM approach for solving the minimization problem given in the equation (19). Here, the result of previous level $j+1$, denoted by $\hat{s}^{(j+1)}$, can be used for initializing after interpolating by factor of two. Although ADMM is well-known and its application for total variation based image restoration is not new [30, 41, 42, 43], implementation of standard ADMM causes some numerical problems because of the spatially varying relative weight $\bar{\beta}$. In the following, we will first describe the formulation that will lead to standard ADMM and then describe the modification necessary to handle the associated numerical issues. The first step in constructing an ADMM algorithm for minimizing composite functionals is to define an equivalent constrained optimization such that the sub-functionals act on different set of variables that are related by means of linear equality constraints. Then writing the augmented Lagrangian [40] for the constrained problems leads to the required ADMM algorithm.

3.1 Constrained formulation and variable splitting

For notational convenience, we switch to vector based notations. Let the $N\times N$ image $s({\bf r})$ be represented by scanned vector ${\bf s}$ in $\mathbb{R}^{N^{2}}$, such that its $i$th element $s_{i}$ is given by $s_{i({\bf r}^{\prime})}=s({\bf r}^{\prime})$ with ${\bf r}^{\prime}=[r_{1}\;r_{2}]$ satisfying $i({\bf r}^{\prime})=r_{2}N+r_{1}$. Let ${\bf m}$ and $\bar{\boldsymbol{\beta}}$ also be defined from $m({\bf r})$ and $\bar{\beta}(\mathbf{r})$ in a similar way with the components denoted by $m_{i}$ and $\bar{\beta}_{i}$. Let ${\bf H}$ be the matrix equivalent of convolving an image with $h({\bf r})$, such that the scanned vector of $(h*s)({\bf r})$ is given by ${\bf H}{\bf s}$. Next, let ${\bf E}^{(j)}$ be the matrix equivalent of interpolation by a factor $2^{i}$. With this, the scanned vector of $E^{(j)}s({\bf r})$ is given by ${\bf E}^{(j)}{\bf s}$. In terms of the new notational scheme, the data fidelity term can be written as

[TABLE]

Similarly, let ${\bf D}_{x}$, ${\bf D}_{y}$, ${\bf D}_{xx}$, ${\bf D}_{yy}$, and ${\bf D}_{xy}$ be the matrices defined from $d_{x}({\bf r})$, $d_{y}({\bf r})$, $d_{xx}({\bf r})$, $d_{yy}({\bf r})$, and $d_{xy}({\bf r})$ for representing convolution operations. Let ${\bf D}_{f}=[{\bf D}_{x}^{T}\;{\bf D}_{y}^{T}]^{T}$ and let ${\bf D}_{s}=[{\bf D}_{xx}^{T}\;{\bf D}_{yy}^{T}\;{\bf D}_{xy}^{T}]^{T}$. Let ${\cal S}({\bf v}):\mathbb{R}^{3}\rightarrow\mathbb{R}^{4}$ be the mapping that returns \left[\begin{array}[]{cc}v_{1}&v_{3}\\ v_{3}&v_{2}\end{array}\right], where ${\bf v}=\{v_{1},v_{2},v_{3}\}\in\mathbb{R}^{3}$ represents the three second order derivatives. Let ${\bf P}_{i}$ be a $2\times 2N^{2}$ matrix having ones at locations $(1,i)$, and $(2,N^{2}+i)$ and zeros at all other locations. Let ${\bf Q}_{i}$ be the $3\times 3N^{2}$ matrix having ones at locations $(1,i)$, $(2,N^{2}+i)$, and $(3,2N^{2}+i)$ and zero at other locations. Then we can use the following substitutions in Eq.(15):

[TABLE]

[TABLE]

With these, the regularization functional can be expressed as

[TABLE]

where we have replaced the images $s$ and $\bar{\beta}$ by the components of their vectorial form. Next, to simplify the task of expressing and comparing the two forms of ADMM, we introduce two more definitions as given below:

[TABLE]

Note that $\mathbf{u}\in\mathbb{R}^{N^{2}}$, $\mathbf{v}\in\mathbb{R}^{2N^{2}}$, and $\mathbf{w}\in\mathbb{R}^{3N^{2}}$. With these, the regularization can be expressed as

[TABLE]

where $\mathbf{1}$ is the vector of ones. Now the minimization of $F^{(j)}(\mathbf{s,H,m})+\lambda R^{(j)}_{sa}({\bf s},{\boldsymbol{\beta}},p)+{\cal B}({\bf E}^{(j)}{\bf s})$ can be equivalently expressed as

[TABLE]

subject to the conditions that

[TABLE]

This constrained formulation of the reconstruction problem leads to the ADMM algorithm, which is essentially a series of cyclic minimization of individual sub-functionals of the above cost. However, through some reconstruction trials, we found that ADMM method obtained from this formulation leaves some artifacts in the reconstruction. These artifact disappear only with very large number of iterations.

Here we present an alternative formulation that leads to a better converging ADMM algorithm. To this end, we first recognize that the cost given in the equation (30) can also be expressed as

[TABLE]

where

[TABLE]

[TABLE]

with ${\bf B}=diag({\boldsymbol{\beta}})$. The corresponding constrained problem becomes

[TABLE]

subject to condition that

[TABLE]

From our experiments, we found that the ADMM steps constructed based on the above constrained formulation leads to better converging algorithm.

3.2 Augmented Lagrangian and the ADMM steps

Writing the ADMM steps for the above problem is straightforward and well-known in the literature. However, for proving the convergence of block coordinate descent method represented by equations (22) and (23), we need to specify the steps here. Further, the constraint of the modified formulation given by equation (37) involves non-circulant matrices and hence, it requires some special consideration. To proceed further, we use the symbol $\bar{\mathbf{s}}$ in the place of ${\mathbf{s}}$ to avoid notational conflict with the iterations of equations (18) and (19). To construct the ADMM algorithm from the above constraint form of the problem, we define

[TABLE]

where $\mathbf{w}$ represents Lagrangian multipliers [40]. Let

[TABLE]

and

[TABLE]

Let $\bar{\mathbf{s}}^{(j)}$ be the vector of length $N^{2}/2^{2j}$ obtained by scanning $E^{(1)}\hat{s}^{(j+1)}(\mathbf{r})$, where $\hat{s}^{(j+1)}(\mathbf{r})$ is the result of previous iteration in the multi-resolution loop of the equations (18) and (19). Then ADMM iterations proceed as follows for $k=0,1,2,\ldots$:

[TABLE]

We will use the partitioned notation \mathbf{w}^{(k)}=\left[\begin{array}[]{c}{\bf w}_{f}^{(k)}\\ {\bf w}_{s}^{(k)}\\ {\bf w}_{0}^{(k)}\end{array}\right] for defining the solution for the above minimization problems.

First note that $L_{a}(\bar{\mathbf{s}},\mathbf{d}_{f},\mathbf{d}_{s},\mathbf{d}_{0},\mathbf{w},\lambda)$ is non-differentiable w.r.t. $\mathbf{d}_{f}$, $\mathbf{d}_{s}$, and $\mathbf{d}_{0}$, and we use the notion of proximal map to carry out the minimization defined in the equations (41), (42), and (43). The minimization w.r.t. $\mathbf{d}_{f}$ given in the equation (41) is essentially the minimization of $\frac{\gamma}{2}\left\|{\bf D}_{f}^{\prime}\mathbf{E}^{(j)}\bar{\mathbf{s}}-{\bf d}_{f}\right\|_{2}^{2}+(\mathbf{w}^{(k)}_{f})^{T}({\bf D}_{f}^{\prime}{\bf E}^{(j)}\bar{\mathbf{s}}-{\bf d}_{f})+\mathcal{N}_{f}(\lambda{\boldsymbol{1}},\mathbf{d}_{f})$, w.r.t. ${\bf d}_{f}$. This can be rewritten as $\frac{\gamma}{2}\left\|\bar{\mathbf{d}}_{f}^{(k)}-{\bf d}_{f}\right\|_{2}^{2}+\mathcal{N}_{f}(\lambda{\boldsymbol{1}},\mathbf{d}_{f})+const.$, where $\bar{\mathbf{d}}_{f}^{(k)}=\mathbf{D}_{f}^{\prime}{\bf E}^{(j)}\bar{\mathbf{s}}^{(k)}+(1/\gamma){\bf w}_{f}^{(k)}$, and $const$ is a term that is independent of ${\bf d}_{f}$. The minimum of this cost is defined as the proximal map of $\mathcal{N}_{f}(\lambda{\boldsymbol{1}},\cdot)$ applied on $\bar{\mathbf{d}}_{f}^{(k)}$ [44]. This map is expressed as [44]

[TABLE]

where $\mathcal{T}(\mathbf{x},t)$ denotes the soft-threshold operator given by

[TABLE]

Similarly, solution to the minimization problem given in equation (42) is the proximal map of $\mathcal{N}_{s}(\lambda{\bf 1},\cdot,p)$ applied on $\bar{\mathbf{d}}_{s}^{(k)}=\mathbf{D}_{s}^{\prime}{\bf E}^{(j)}\bar{\mathbf{s}}^{(k)}+(1/\gamma){\bf w}_{s}^{(k)}$. Let $|||\cdot|||_{t}$ denote the operator that applies soft-thresholding on the Eigen values of its matrix arguments and returns the resulting matrix. Then the proximal map of $\mathcal{N}_{s}(\lambda{\bf 1},\cdot,p)$ can be expressed as [26]

[TABLE]

where

[TABLE]

Finally, the minimum defined in equation (43) is the proximal map of ${\cal B}({\bf d}_{0})$ applied on $\bar{\mathbf{d}}_{0}^{(k)}={\bf E}^{(j)}\bar{\mathbf{s}}^{(k)}+(1/\gamma){\bf w}_{0}^{(k)}$. This is denoted as $\mathbf{d}_{0}^{(k+1)}=\mathcal{P}_{u}(\bar{\mathbf{d}}_{0}^{(k)})$, where $\mathcal{P}_{u}(\cdot)$ denotes the clipping of components of the its vectors onto the range $[0,b]$, with $b$ as the user-defined upper bound.

The cost $L_{a}(\bar{\mathbf{s}},\mathbf{d}_{f},\mathbf{d}_{s},\mathbf{d}_{0},\mathbf{w},\lambda)$ is differentiable w.r.t. $\bar{\bf s}$, and the minimization defined in the equation can be obtained by equating gradient to zero. Let \mathbf{M}=\left[\begin{array}[]{c}{\bf D}_{f}^{\prime}\\ {\bf D}_{s}^{\prime}\\ \bf I\end{array}\right] and let \mathbf{d}^{(k+1)}=\left[\begin{array}[]{c}{\bf d}_{f}^{(k+1)}\\ {\bf d}_{s}^{(k+1)}\\ {\bf d}_{0}^{(k+1)}\end{array}\right]. Then, the solution of the last minimization is given by the following equation,

[TABLE]

where $\mathbf{v}_{k+1}=\mathbf{E}^{(j)T}\mathbf{M}^{T}\left(\mathbf{d}^{(k+1)}-(1/\gamma)\mathbf{w}^{(k+1)}\right)$. Note that this equation has to be solved iteratively since $\mathbf{M}$ is composed of non-circulant matrices. We use conjugate gradient method for solving this problem. To speed-up, we use the inverse of the following approximation of the matrix $\mathbf{A}_{j}$ in the equation (56), $\hat{\mathbf{A}}_{j}$, as the preconditioner:

[TABLE]

All the matrices in the above product are circulant except $\mathbf{E}^{(j)}$. However, the product, $\hat{\mathbf{A}}_{j}$, is circulant because of the special structure of $\mathbf{E}^{(j)}$. Hence, the preconditioning, i.e., multiplying by the inverse of $\hat{\mathbf{A}}_{j}$, is equivalent to applying the inverse of a discrete filter. The following proposition gives the expression for this filter.

Proposition 2**.**

Let $u(z_{1},z_{2})$ be the z-transform of $u(\mathbf{r})=\frac{1}{64}[1\;4\;6\;4\;1]^{T}[1\;4\;6\;4\;1]$, and let $u_{j}(z_{1},z_{2})=\prod_{i=0}^{j-1}u(z^{2^{i}}_{1},z^{2^{i}}_{2})$. Let $u_{j}(\mathbf{r})$ be the inverse z-transform of $u_{j}(z_{1},z_{2})$. Then the filter equivalent of $\hat{\mathbf{A}}_{j}$ is the $2^{j}$-fold decimation of $B(\mathbf{r})=u_{j}(\mathbf{r})*u_{j}(-\mathbf{r})*[1+d_{x}(\mathbf{r})*d_{x}(-\mathbf{r})+d_{y}(\mathbf{r})*d_{y}(-\mathbf{r})+d_{xx}(\mathbf{r})*d_{xx}(-\mathbf{r})+d_{yy}(\mathbf{r})*d_{yy}(-\mathbf{r})+d_{xy}(\mathbf{r})*d_{xy}(-\mathbf{r})]$.

Next, applying these ADMM steps described above for solving the minimization problem of equation (23) is nearly identical except the fact that the up-sampling matrix $\mathbf{E}$ is replaced by identity matrix because the cost is not defined through up-sampling. Here the size of the variable is the same as the size of the measured image. The initialization for ADMM iteration, $\bar{\mathbf{s}}^{(0)}$, now simply comes from $s_{(k)}$, which is the result of previous iteration of the BCD loop represented by equations (22) and (23). Next, the following proposition confirms the convergence of the block coordinate descent algorithm specified by the equations (22) and (23).

Proposition 3**.**

The block coordinate descent method represented by the equations (22) and (23) with the problem of equation (23) solved by ADMM method described above, converges to a local solution of the problem (10) if $J_{sa}(s_{(k+1)},\beta_{(k+1)},\tau,h,m)<J_{sa}(s_{(k)},\beta_{(k+1)},\tau,h,m)$ and if there is no $s$ for which $F(s,h,m)$, $R_{1}(s)=\sum_{\bf r}\|({\bf d}_{1}*s)({\bf r})\|_{2}$, and $R_{2}(s)=\sum_{\bf r}\|({\bf H}_{2}*s)({\bf r})\|_{2}$ will have zero value simultaneously.

4 Experimental results

For evaluating the restoration performance of the proposed COROSA approach, we considered deconvolution of Total Internal Reflection Fluorescence (TIRF) microscopy images and the reconstruction of Magnetic Resonance Images (MRI) from under-sampled Fourier data. These problems involve different measurement and noise models and are hence good candidates for evaluating the performance of the proposed approach alongside the state-of-the-art methods. We compare COROSA with second order TV (TV2) [24], Hessian-Schatten norm regularization (HS) [26], combined order TV (COTV) [33] and TGV2 [36] regularization methods. We also implemented the combined order TV formulation with Hessian-Schatten norm regularization replacing the original second order TV term, for the purpose of comparison. We refer to this method as Combined Order Hessian-Schatten (COHS) regularization. For the HS functional, we found that setting $p=1$ yielded the best performance. We also include the result of the multi-resolution loop represented by equations (18) and (19) in the comparison (without BCD iterations of equations (22) and (23)). We denote this by COROSA-I. For objective comparison, we use Signal to Noise Ratio (SNR) [45] and Structural Similarity Index (SSIM) [46] scores. In tables, we use COR. and COR-I to denote COROSA and COROSA-I respectively.

The smoothing parameter $\lambda$ was tuned for best performance in terms of SSIM and SNR by using original reference images as done by most methods that focus on the design of regularization. In the case of COTV, and COHS, additional tuning is required to fix the parameters determining the first and second order TV terms. In this regard, we set the relative weights between first and second order derivatives so as to yield the lowest regularization functional cost. This ensures that only $\lambda$ is required to be tuned using the reference images.

The spatial weight $\beta(\mathbf{r})$ can be determined through optimization problem defined in (18) and the corresponding result is given in Proposition 1, with the parameter $\tau$ chosen to be a scalar parameter. However, we observed that it is advantageous to make $\tau$ spatially variant for the following reasons: (i) in the regions of low intensity, it is advantageous to make $\beta(\mathbf{r})$ less sensitive to variations in the relative magnitude of first- and second-order derivatives, and hence $\tau$ has to be larger; (ii) in the regions of high intensity, it is advantageous to make $\beta(\mathbf{r})$ more sensitive to variations in the relative magnitude of these terms, and hence $\tau$ has to be smaller. In short, it is advantageous to make $\tau$ spatially variant and inversely proportional to some approximate estimate of the required image, say $\bar{f}(\mathbf{r})$. We use the following strategy to get this approximate estimate: (i) for solving the problem given the equation (18), we use $f$ itself as the approximate estimate, i.e., we set $\bar{f}(\mathbf{r})=f(\mathbf{r})$. Next, in the block-coordinate descent loop specified by the equations (22) and (23), we keep $\bar{f}(\mathbf{r})$ fixed as the image used to initialize the loop, that is we set $\bar{f}(\mathbf{r})=s_{(0)}(\mathbf{r})$. From $\bar{f}(\mathbf{r})$, we compute $\tau(\mathbf{r})$ as follows: we compute the image $\exp(-100*\bar{f}^{2}(\mathbf{r}))$ and then rescale it to the range $[0.01,100]$. This scheme worked well for all our test cases, and hence we kept this scheme for determining $\tau(\mathbf{r})$ in all test cases.

In the first experiment, we evaluate the proposed method for deconvolution of TIRF images. To obtain realistic ground-truth models, we measured TIRF images with negligible noise and deconvolved them by a simple inverse filtering. The resulting images became the ground truth models, which were used to simulate noisy blurred images for evaluating the proposed method and for comparing it with the existing methods. The images for generating ground-truth model were acquired from samples containing labelled Actin and Tubulin. The Actin images were acquired by staining with phalloidin-488 and an EMTB-mCherry transgene was used to get the Tubulin image. Wavelength for excitation was 491nm for Actin samples and 561nm for Tubulin sample. All images were acquired using a 100x objective lens with numerical aperture 1.45 NA. The exposure time was set to 300ms, and this was sufficient to get nearly noise-free images with a quality that is adequate for using them as ground-truth models. These are shown in Figure 1.

For simulating noisy blurred images from ground-truth models, we used blurring kernel with a bandwidth that is slightly lower than the bandwidth of the system that measured the model images. This was done by setting NA = 1.4, while computing the required blurring kernel. We used the following model to generate test images:

[TABLE]

Here $h(\mathbf{r})$ represents the simulated 2D TIRF PSF with NA = 1.4, and with other parameters set to be identical to that of the system that measured the models. Next, ${\mathcal{P}}(\cdot)$ refers to Poisson process with $\gamma_{p}$ representing the scale factor for photon count, and $\eta$ is AWGN of variance $\sigma_{\eta}$. We use this mixed noise model to make the test data realistic, although we do not use the log-likelihood functional of the mixed model in the formulation of our method. For each reference model, we set $\sigma_{\eta}=1$ and generated four noisy blurred images with $\gamma_{p}$ = 5, 10, 20, and 30. The four noisy images corresponding to ActinSample1 are shown in Figure 2 as an example. It has to be pointed out that we did not include TGV in this experiment because of its poor performance, as was observed in [35]. In this regard, it is worthwhile to note that TGV has been used notably only for MRI reconstruction. The deconvolution results in terms of both SSIM and SNR scores are presented in Table 1. The scores show that the proposed COROSA outperforms other methods in most cases, with the performance advantage significant when the noise is high. It is to be noted that TV2 gives slightly better scores than COROSA as the measurements become less noisy. This is due to the fact that the effect of spatial adaptiveness is negligible, and the fact that, TV2 can converge to the minimum more accurately because of its simplicity, which takes over the advantage of spatial adaptiveness. However, in such low-noise cases, the difference in the score is much lower than the advantage that COROSA has in noisy cases. Figure 3 shows a set of restored images corresponding to ActinSample1 image and $\gamma_{p}=20$. In terms of SSIM, the difference between TV2 and COROSA is 0.006 and difference in SNR is 0.2dB. However, as evident from the displayed images, there is a clear visual improvement in the result of COROSA. Another observation is that COROSA significantly outperforms COTV and COHS, because of the spatial adaptivity.

In the second experiment, we considered reconstruction of MRI images from undersampled k-space measurements. The source MRI images are shown in Figure 5. It has to be mentioned that some of these reference images have been obtained using traditional reconstruction methods and hence contain artifacts. However, this does not affect the validity of our comparisons done here, because, all the methods evaluated in this paper including the proposed method remove these artifacts. Hence, we focus on how well the actual structures are reproduced by various methods in the presence of noise and undersampling, and ignore the removal of artifacts present in the reference images. For sampling, we generated random and spiral trajectories using the MATLAB code given by Chauffert et al. [47]. In addition, we modified the spiral trajectory by filling in low frequency region, since the default spiral had low sample density in the low frequency region. Both sampling trajectories with 10% and 20 % densities are shown in Figure 5. For simulating the thermal noise, we added white Gaussian noise to the k-space data. To control the noise, we use the strategy of [10]: the Gaussian variance is adjusted such that, if added to both real and imaginary parts of full Fourier transforms of the image, it results in PSNR of 10dB and 20dB upon the Fourier inversion. We generated six sample sets from each MRI image as follows: for each sampling trajectory, we generated 20% sample set with two noise levels (10dB and 20dB), and 10% sample set with one noise level (20dB). This make a total of 36 sample sets. The comparisons in terms of SSIM and SNR scores for all cases are given in Table 2.

The reconstruction scores in Table 2 show that, in most cases, the proposed COROSA approach is better than all competitive methods, with COTV and COHS being the nearest in terms of scores. In terms of SSIM, COROSA is always better, and in terms of SNR, COROSA is better in most of the cases. In general, combined order methods perform better in MRI reconstruction, when compared to TIRF restoration. Figure 6 shows the reconstruction results for MRI 6 image with 10% sampling and 20dB PSNR samples. It is clear from the figure that COROSA suppresses background artifacts with the help of spatial adaptation better than TV2 and HS. At the same time, COROSA also avoids staircase artifacts that is normally caused if TV1 is used without spatial adaptation as present in COTV and COHS reconstructions. A selected region of this reconstruction is shown in Figure 7 to further emphasize this fact, along with the relative weight $\beta(\mathbf{r})$. In the figure, white regions correspond to locations where first order TV term is dominant, and black regions correspond to locations where second order TV term is dominant. Grey regions represent locations where both orders are weighted equally. It is clear from the figure that spatial weights follow the local intensity structure, leading to the elimination of staircase artifacts. On the other hand, such artifacts are evident in COTV reconstruction, which uses similar cost function and optimization framework, albeit with constant global relative weights. Recall that SSIM of COROSA for this selected region is 0.7663 and that of COTV is 0.7324. Figure 8 shows the evolution of adaptive weight $\beta({\bf r})$ across different stages of multi-resolution and block coordinate descent iterations. It can be noted that $\beta({\bf r})$ follows fine variations in the relative magnitudes of first- and second-order derivatives that are not obvious from the visual inspection of the displayed intensity images.

We have also given the reconstruction results for MRI 1 image with 10% sampling and 20dB PSNR samples in Figure 9, confirming that the visual improvement with COROSA is consistent with the higher scores. Further, COTV and COHS introduce prominent artifacts in the reconstruction as seen in zoomed-in view of the same set of images displayed in Figure 10 along with $\beta(\mathbf{r})$. It is also clear that these artifacts are not the remnants of any artifacts from the original model, but were produced by COTV and COHS methods. On the other hand, COROSA results do not have any artifacts and the adaptive weights clearly pick out image boundaries from uniform intensity regions. There are two exceptions in the results that do not confer to the pattern exhibited by other test cases. In the first case, COTV was giving better SNR than COROSA, with comparable SSIM scores for reconstruction of MRI 1, MRI 2 and MRI 4 images from spiral trajectory. This is because of the modification in the spiral trajectory, where low frequency regions are filled and hence the advantage of multiresolution framework utilized in COROSA becomes less significant. In the second case, COTV was giving higher SNR score for MRI 4 image reconstruction from random sampling, while having lower SSIM score when compared to COROSA. When we evaluated the corresponding results visually, we found that the reference MRI 4 image itself has block-like piecewise constant regions. Because of this, COTV results have higher SNR compared to COROSA because of the fact that COROSA keeps a minimum amount of second-order smoothing even in the piece-wise constant regions. However, in terms of SSIM score and visual quality, the difference is insignificant. Next, it is worthwhile to note that, TGV has the ability to be spatially adaptive because of the auxiliary variable, $\mathbf{p}$, and it also retains convexity, which makes it quite attractive. However, our method outperforms TGV significantly. Further, surprisingly, even basic non-adaptive methods such as COTV and COHS outperform TGV in many test cases. A possible reason that we inferred based on some reconstructions trials, is that, the inferior performance of TGV is due to the lack of efficient optimization method to handle the auxiliary variable. Specifically, the convergence of all known optimization methods proposed for TGV is highly dependent on the value of the smoothing parameters $\alpha_{1}$, $\alpha_{2}$ and $\lambda$. This leads to the inferior performance of TGV although it is based on rich and elegant mathematical formulation.

Recently, there has been significant interest in deep learning based image restoration [12, 13, 14, 15] and dictionary learning methods [10, 16, 17]. Although these methods require training samples, unlike our method which does not need any training, we make comparison with representative methods from these categories to obtain a perspective. We chose the method of Deep de-aliasing GAN (DAGAN) proposed by Yang et al. [12] and the dictionary learning method of Ravishankar et al. [10] for this purpose, and utilized the code provided by the authors. For DAGAN we chose the best Pixel-Frequency-Perceptual-GAN-Refinement (PFPGR) topology of the network for comparison, and also used the same dataset used by the authors with identical settings. Since our method and the dictionary learning methods do not require training samples, we evaluated DAGAN with varying number of training images: we evaluated this method with 10$\%$, 50$\%$ and 100$\%$ of the 15972 training images used by the authors. We selected 6 images from the set of 50 test images provided by the authors for evaluation. These images are shown in Figure 11.

As the first part of this experiment, we tested reconstruction from under-sampled data without noise using random trajectory. It has to be noted that DAGAN was not designed to take complex measurements. Instead, it takes the absolute of the inverse Fourier transform of the zero-filled k-space values. To eliminate this disadvantage for DAGAN, we made the random sampling symmetric in Fourier domain, yielding real images. The reconstruction results for 10% and 20% sampling ratios are given in Table 3. As the results indicate, COROSA performs better than DAGAN and DL-MRI consistently with $20\%$ samples, in terms of both SNR and SSIM. With $10\%$ sampling, COROSA performs better than the other methods in most cases; however, in some cases, DAGAN has higher SNR scores than COROSA although COROSA is still better in terms of SSIM scores. This is due to the fact that DAGAN reconstructions have artifacts, and SSIM score is more sensitive to artifacts than SNR measure. Further, on visual inspection, we found that DL-MRI reconstructions are over-smoothed, especially the reconstruction with 10% samples. We also noted that the reconstructions using COROSA do not create artifacts in all test cases including the cases where DAGAN has better SNR, as seen in Figure 12. On the other hand, DAGAN reconstruction shows spurious structures introduced during reconstruction. The main point that we want to emphasize in this regard is that COROSA is able to adapt to image structure, due to the robust initialization from multi-resolution framework. In summary, COROSA outperforms both DAGAN and DL-MRI in all cases in terms of SSIM score.

In the next part, we repeated the first part with k-space samples corrupted by AWGN. The noise variance was adjusted such that a specific SNR is achieved with the zero-filled inverse transformed images. These specific SNR’s were chosen to be 0.5dB less than that of noise-free counter-parts for the case of $10\%$ sampling, and to be 1.5dB less for the case of $20\%$ sampling. The results are shown in Table 4, which demonstrate that COROSA outperforms both DAGAN and DL-MRI in terms of both SNR and SSIM scores, except for two cases. In these two cases, SNR of DAGAN is better. It should be re-emphasized that, in terms of SSIM scores, COROSA is better than both DAGAN and DL-MRI in all test cases. Figure 13 shows the reconstruction results corresponding to $20\%$ sampling, which confirms the superiority of COROSA reconstruction.

With regards to the computation time, we ran all MATLAB algorithms on Core i7-3770 CPU with 8 GB RAM. We found that COROSA restoration for $256\times 256$ image is completed in 185.7s. For the same task, TGV2, TV2, HS, COTV, COHS, DAGAN and DL-MRI take 24.2s, 2.3s, 2.9s, 2.8s, 3.6s, 0.2s and 507s respectively. Here, DAGAN restoration was performed using a trained network and the training itself requires approximately 4 hours for 10% training set, 20 hours for 50% training set and 40 hours for the full training set of 15972 images. Regarding overall computational complexity, our algorithm can be considered to be of $O(N^{2}\log(N))$ complexity for an $N\times N$ image, because of FFT operations that dominate over other filtering operations.

For detailed analysis, COROSA can be viewed mainly in two parts: the multiresolution loop involving alternative refinements of the required image and the adaptive weights iterations (equations (18) and (19)) and the block coordinate descent iteration involving alternative refinements of the required image and the adaptive weights in the final resolution (equations (22) and (23)). The first part is essentially similar to the second except for the cost of expansion operation. Ignoring the cost of expansion operation, one round of refinement of the image and the adaptive weight costs 30+83$N_{i}$+44$N_{CG}N_{i}$ additions, 35+78$N_{i}$+43$N_{CG}N_{i}$ multiplications, $N_{i}$+1 square root operations and 2+6$N_{i}$+2$N_{CG}N_{i}$ FFT operations, where $N_{i}$ and $N_{CG}$ are the number of ADMM iterations and the number of CG iterations corresponding to the minimization of equation (19), respectively. This complexity is multiplied by $N_{f}+K$ where $K$ is the number of multiresolution levels, and $N_{f}$ is the number of block-coordinate descent iterations (equations (22) and (23)).

5 Conclusion

We developed a novel form of regularization scheme that combines first- and second-order derivatives in a manner that is adaptive to the image structure. The adaptation is achieved by the fact that the relative weight that combines first- and second-order derivatives is determined by the same cost functional along with an additional regularization for preventing rapid variations in the adaptive weights. We used isotropic TV for the first-order term, and Hessian-Schatten norm for the second order term. We constructed an iterative method for minimizing the resulting non-convex and non-differentiable cost functional, and we proved the convergence of the iterative method. We demonstrated that, the proposed regularization method outperforms notable regularization methods in the literature when the noise is high in the case of deblurring for microscopy. Further, in MRI reconstruction, we demonstrated that the proposed regularization method outperforms existing regularization methods, and a recently proposed Deep GAN method when the noise and/or under-sampling are high.

Appendix

Proof of Proposition 1

Proof: Solution $\bar{\beta}(\mathbf{r})\in[0,1]^{N\times N}$ is given by

[TABLE]

For the case where $d(\mathbf{r})=\|({\bf d}_{1}*(E^{(j)}s))({\bf r})\|_{2}-\left\|{\boldsymbol{\zeta}}(({\bf H}_{2}*(E^{(j)}s))({\bf r}))\right\|_{p}=0$, equating derivative w.r.t. $\beta$ to zero gives $2\tau\bar{\beta(\mathbf{r})}=\tau$ which gives the solution $\bar{\beta}(\mathbf{r})=0.5$. Now we examine the case where $d(\mathbf{r})$ is non-zero. For this, we first define the following:

[TABLE]

Next, we note that $\bar{\beta}(\mathbf{r})$ is less than $0.5$ if $v_{1}(\mathbf{r})>v_{2}(\mathbf{r})$ and greater than $0.5$ if $v_{1}(\mathbf{r})<v_{2}(\mathbf{r})$. Hence, the following transformation will reduce the range of minimization variable to $[0,0.5]$:

[TABLE]

Translating the optimization problem in terms of ${\beta}_{r}(\mathbf{r})$ gives

[TABLE]

where

[TABLE]

Equating derivative w.r.t. $\beta_{r}(\mathbf{r})$ to zero yields the quadratic equation

[TABLE]

The roots are given by

[TABLE]

where $\zeta=\frac{2\tau}{|d(\mathbf{r})|}$. From the form of the expressions, it is clear that $\beta_{r}^{+}(\mathbf{r})$ is positive. Also, as $|d(\mathbf{r})|$ ranges from [math] to $\infty$, $\beta_{r}^{+}(\mathbf{r})$ ranges from [math] to $0.5$. Hence, $\beta_{r}^{+}(\mathbf{r})\in[0,0.5]^{N\times N}$, which is the required solution. Substituting $\beta_{r}^{+}(\mathbf{r})$ in the equation (60) gives the expression of the equation (21). Since $\beta_{r}^{+}(\mathbf{r})\in[0,0.5]$, it is clear that $\beta(\mathbf{r})$ is guaranteed to be in the range $[0,1]$. In fact, it is guaranteed to be in the interval $(0,1)$ because $L$ will be infinity if $\beta(\mathbf{r})\in\{0,1\}$.

Proof of Proposition 2

The matrix in the equation (57) can be expressed as

[TABLE]

Note that $\mathbf{E}^{(j)}$ is the matrix equivalent of $j$-stage implementation of two-fold upsampling; this upsampling is realized as an expansion by a factor of two, which is inserting a zero after each pair of samples along both axes, and then filtering by $u(\mathbf{r})=\frac{1}{64}[1\;4\;6\;4\;1]^{T}[1\;4\;6\;4\;1]$. This operation is also equivalent to the single stage implementation involving $2^{j}$ fold expansion, and then filtering by $u_{j}(z_{1},z_{2})=\prod_{i=0}^{j-1}u(z^{2^{i}}_{1},z^{2^{i}}_{2})$ where $u(z_{1},z_{2})$ is the $z$-transform of $u(\mathbf{r})$ [48]. Let $u_{j}(\mathbf{r})$ be the inverse $z$-transform of $u_{j}(z_{1},z_{2})$. From the structure of $\mathbf{E}^{(j)}$, we infer that multiplication by $\mathbf{E}^{(j)T}$ is equivalent to convolution by $u_{j}(-\mathbf{r})$ followed by decimation by a factor of $2^{j}$ along both axes, which is the operation of skipping $2^{j}-1$ samples for each block of $2^{j}$ samples.

Next, the matrix ${\bf D}_{x}$ represents convolution by $d_{x}(\mathbf{r})$ and ${\bf D}_{x}^{T}$ represents convolution by $d_{x}(-\mathbf{r})$. The other matrices within the square brackets are similarly interpreted.

Hence the operation equivalent to multiplication by $\hat{\mathbf{A}}_{j}$ is the following three stage operations in sequence: (i) $2^{j}$ fold expansion; (ii) filtering by $B(\mathbf{r})=u_{j}(\mathbf{r})*u_{j}(-\mathbf{r})*[1+d_{x}(\mathbf{r})*d_{x}(-\mathbf{r})+d_{y}(\mathbf{r})*d_{y}(-\mathbf{r})+d_{xx}(\mathbf{r})*d_{xx}(-\mathbf{r})+d_{yy}(\mathbf{r})*d_{yy}(-\mathbf{r})+2d_{xy}(\mathbf{r})*d_{xy}(-\mathbf{r})]$; (iii) $2^{j}$ fold decimation. This three stage operation is equivalent to convolving by $2^{j}$ fold decimated version of $B(\mathbf{r})$ [48].

Proof of Proposition 3

Our proof will be based on Zangwill’s global convergence theorem. It states three conditions to be satisfied by the iterates to ensure convergence. These conditions translated for our problem are the following: (i) the sequence $\{(s_{(k)},\beta_{(k)}\}_{i=1,2,\ldots}$ is a descent sequence, i.e., the sequence should satisfy $J_{sa}(s_{(k+1)},\beta_{(k+1)},\tau,h,m)<J_{sa}(s_{(k)},\beta_{(k)},\tau,h,m)$; (ii) the sequence of iterates should be contained in a compact set; (iii) the mapping that generates the iterates should be closed, i.e., if $\mathcal{M}$ is mapping such that $(s_{(k+1)},\beta_{(k+1)})=\mathcal{M}(s_{(k)},\beta_{(k)})$, then it should be a closed mapping.

To verify the second condition, note that $\{(s_{(k)},\beta_{(k)}\}_{i=1,2,\ldots}$ is within the sub-level set satisfying $J_{sa}(s,\beta,\tau,h,m)\leq J_{sa}(s_{(0)},\beta_{(0)},\tau,h,m)$. This is a bounded set because the function $J_{sa}(s,\beta,\tau,h,m)$ is bounded below and has empty null space. Note that a bounded set in Euclidean space is compact.

Note that $J_{sa}(s_{(k)},\beta_{(k+1)},\tau,h,m)<J_{sa}(s_{(k)},\beta_{(k)},\tau,h,m)$. This is because $\beta_{(k+1)}$ is computed by exact of minimization of $J_{sa}(s_{(k)},\beta,\tau,h,m)$ with respect to $\beta$. Next, $s_{(k+1)}$ is computed by iterative minimization of $J_{sa}(s,\beta_{(k+1)},\tau,h,m)$ with respect to $s$ using ADMM. By assumption, $J_{sa}(s_{(k+1)},\beta_{(k+1)},\tau,h,m)<J_{sa}(s_{(k)},\beta_{(k+1)},\tau,h,m)$. Hence we, have $J_{sa}(s_{(k+1)},\beta_{(k+1)},\tau,h,m)<J_{sa}(s_{(k)},\beta_{(k)},\tau,h,m)$. This verifies the first condition.

To verify the third condition, we will first verify that each cycle of ADMM is a continuous mapping. We first consider the mapping $(\mathbf{d}_{f}^{(k+1)},\mathbf{d}_{s}^{(k+1)},\mathbf{d}_{0}^{(k+1)})=\mathcal{K}(\bar{\mathbf{s}}^{(k)},\bar{\beta})$ represented by the equations (41), (42), and (43). Since these equations represent exact single step minimizations of convex sub-functionals, the mapping $\mathcal{K}(\cdot)$ is continuous. Next, consider the mapping $\bar{\mathbf{s}}^{(k+1)}=\mathcal{L}(\mathbf{d}_{f}^{(k+1)},\mathbf{d}_{s}^{(k+1)},\mathbf{d}_{0}^{(k+1)},\bar{\beta})$ represented by the equation (44). Here, we ignore the presence of $\mathbf{w}^{(k)}$ since it is generated by a simple affine transformation. This minimization is implemented by conjugate gradient iterations. Conjugate gradient iteration with any number of steps is equivalent to a minimization of convex quadratic function within a subspace and it is also continuous. Hence the mapping, $\mathcal{M}_{2}(\cdot,\bar{\beta})=\mathcal{L}(\mathcal{K}(\cdot,\bar{\beta}),\bar{\beta})$, which represents one cycle of ADMM, is continuous. If the minimization given in equation (23) is implemented with $N_{a}$ cycles of ADMM with initialization $s_{(k)}$, we can represent this as $s_{(k+1)}=\mathcal{M}_{2}^{N_{a}}(s_{(k)},\bar{\beta})$. Now, the minimization in the equation (22) is continuous operation because it is implemented by exact minimization, and function with respect to $\beta$ alone is convex. The result, $\bar{\beta}$, is a function of $s_{(k)}$. We denote the minimization operation specified by the equation (22), by $(s_{(k)},\bar{\beta})=\mathcal{M}_{1}(s_{(k)})$. Hence one cycle of block coordinate descent can be represented as $s_{(k+1)}=\mathcal{M}_{2}^{N_{a}}(\mathcal{M}_{1}(s_{(k)}))$. The above mapping is also continuous since it is composition of continuous mappings. Since continuity is a special case of closedness, we can say that the mapping that generates the iterates is a closed mapping.

References for MRI Images

MRI 1, https://www.xraygroup.com.au/index.php/our-services/mri

MRI 2, https://www.healthcare.siemens.es/magnetic-resonance-imaging/options-and-upgrades/upgrades/magnetom-trio-upgrade/use

MRI 3, https://radiopaedia.org/images/208569

MRI 4, https://www.usoccdocs.com/contract-service-mri/

MRI 5, https://www.researchgate.net/post/[What_do_you_think_diagnosis_of_this_pediatric_brain_MRI2](What_do_you_think_diagnosis_of_this_pediatric_brain_MRI2)

MRI 6, https://www.healthcare.siemens.co.uk/magnetic-resonance-imaging/0-35-to-1-5t-mri-scanner/magnetom-aera/use