PCM-TV-TFV: A Novel Two Stage Framework for Image Reconstruction from Fourier Data

TL;DR

This paper introduces a two-stage image reconstruction framework from Fourier data that combines projection and edge-guided regularization, outperforming traditional methods in accuracy and efficiency.

Contribution

The novel PCM framework integrates a multi-scale projection with edge-guided regularity, enhancing image reconstruction from Fourier measurements.

Findings

Superior reconstruction quality in numerical tests

Effective handling of non-uniform Fourier data

Outperforms classical one-stage methods

Abstract

We propose in this paper a novel two-stage Projection Correction Modeling (PCM) framework for image reconstruction from (non-uniform) Fourier measurements. PCM consists of a projection stage (P-stage) motivated by the multi-scale Galerkin method and a correction stage (C-stage) with an edge guided regularity fusing together the advantages of total variation (TV) and total fractional variation (TFV). The P-stage allows for continuous modeling of the underlying image of interest. The given measurements are projected onto a space in which the image is well represented. We then enhance the reconstruction result at the C-stage that minimizes an energy functional consisting of a fidelity in the transformed domain and a novel edge guided regularity. We further develop efficient proximal algorithms to solve the corresponding optimization problem. Various numerical results in both 1D signals and 2D images have also been presented to demonstrate the superior performance of the proposed two-stage method to other classical one-stage methods.

Tables

Table 1. Table 1: Numerical comparison of the different models on different noise level.

$\sigma$	Model	ssim	psnr	snr	rela_err
\hlineB3 0.1	Proposed	0.9283	39.14	33.58	0.0161
	${\ell }^1$ regularization	0.8478	36.45	30.88	0.0220
	SS-TV	0.8466	35.72	30.16	0.0239
	Tikhonov	0.7640	32.17	26.60	0.0360
\hlineB3 0.4	Proposed	0.8505	29.41	23.84	0.0495
	${\ell }^1$ regularization	0.7092	26.21	20.64	0.0715
	SS-TV	0.6830	23.81	18.24	0.0942
	Tikhonov	0.3041	21.09	15.52	0.1289
\hlineB3 0.7	Proposed	0.7096	24.87	19.30	0.0834
	${\ell }^1$ regularization	0.6221	21.93	16.36	0.1170
	SS-TV	0.6816	22.60	17.03	0.1083
	Tikhonov	0.1414	15.51	9.951	0.2445

Table 2. Table 2: Numerical details for the 2D reconstruction in Fig.   7 with bias error, σ = 0.6 𝜎 0.6 \sigma=0.6 .

\hlineB3 Model	psnr	snr	rela_err
\hlineB3 SS-TV	33.8196	20.87	0.086
PCM-TV	35.6589	22.7092	0.069
\hlineB3

Table 3. Table 3: Numerical results for denoising with different regularities.

\hlineB3 Model	psnr	snr	rela_err
\hlineB3 Noisy	20.9199	6.0853	0.3914
TV	31.0007	16.1661	0.1226
TFV	31.6910	16.8565	0.1133
TV-TFV	32.1948	17.3602	0.1069
\hlineB3

Equations

$$\min \mathcal{L}(x)+\mathcal{R}(x),$$

$$\hat{f}=\mathcal{S}\mathcal{F}f+ϵ,$$

$$\mathcal{F}f(w)={\int }_{{\mathbb{R}}^2}f(x)\exp (-2\pi iw\cdot x)dx.$$

$$\tilde{f}=\sum _{j=1}^n({\mathbf{c}}_f{)}_j{\varphi }_j,$$

$${\mathbf{c}}_f={\mathrm{argmin}}_{\mathbf{c}}\Vert \mathcal{S}\mathcal{F}\Phi \mathbf{c}-\hat{f}{\Vert }_2^2,$$

$$g=\sum _{j=1}^n({\mathbf{c}}_g{)}_j{\varphi }_j,$$

$$\underset{{\mathbf{c}}_g}{\min }\frac{1}{2}\Vert {\mathbf{c}}_g-{\mathbf{c}}_f{\Vert }^2+\mathcal{R}({\mathbf{c}}_g),$$

$${}_a{D}_x^{\alpha }f(x)=\frac{1}{\Gamma (n-\alpha )}\frac{d^n}{d{x}^n}{\int }_a^x(x-\tau {)}^{n-\alpha -1}f(\tau )d\tau ,$$

$${}_x{D}_b^{\alpha }f(x)=\frac{(-1{)}^n}{\Gamma (n-\alpha )}\frac{d^n}{d{x}^n}{\int }_x^b(\tau -x{)}^{n-\alpha -1}f(\tau )d\tau ,$$

$${}_a{D}_b^{\alpha }f(x)=\frac{1}{2}{(}_a{D}_x^{\alpha }f(x)+{}_x{D}_b^{\alpha }f(x)\cdot (-1{)}^n),$$

$$\Gamma (\alpha )={\int }_0^{\infty }{t}^{\alpha -1}{e}^{-t}dt.$$

$$x_i=\frac{(i-1)(b-a)}{n}+a,i=1,2,...,n.$$

$$L^{(\alpha)}_{n}=n\left(\begin{array}[]{cccccc}w_{0}^{\alpha}&w_{1}^{\alpha}&\cdots&\cdots&w_{n-1}^{\alpha}&w_{n}^{\alpha}\\ 0&w_{0}^{\alpha}&w_{1}^{\alpha}&\ddots&\ddots&w_{n-1}^{\alpha}\\ 0&0&\ddots&\ddots&\ddots&\vdots\\ \vdots&\ddots&\ddots&w_{0}^{\alpha}&w_{1}^{\alpha}&\vdots\\ 0&\ddots&\ddots&0&w_{0}^{\alpha}&w_{1}^{\alpha}\\ 0&0&\cdots&0&0&w_{0}^{\alpha}\end{array}\right)$$

$$R^{(\alpha)}_{n}=n\left(\begin{array}[]{cccccc}w_{0}^{\alpha}&0&0&\cdots&0&0\\ w_{1}^{\alpha}&w_{0}^{\alpha}&0&\ddots&\ddots&0\\ \vdots&\ddots&\ddots&\ddots&\ddots&\vdots\\ \vdots&\ddots&\ddots&w_{0}^{\alpha}&\ddots&0\\ w_{n-1}^{\alpha}&\ddots&\ddots&w_{1}^{\alpha}&w_{0}^{\alpha}&0\\ w_{n}^{\alpha}&w_{n-1}^{\alpha}&\cdots&\cdots&w_{1}^{\alpha}&w_{0}^{\alpha}\end{array}\right),$$

$$w_j^{\alpha }=(1-\frac{1+\alpha }{j})w_{j-1}^{\alpha },j=1,2,...,n.$$

$$C_n^{(\alpha )}=\frac{1}{2}(L^{(\alpha )}+R^{(\alpha )}).$$

$$u_x^{(\alpha )}=(I_n\otimes C_n^{(\alpha )})U,$$

$$u_y^{(\alpha )}=(C_n^{(\alpha )}\otimes I_n)U.$$

$$w_0=1,w_1=-1.3,w_2=0.195,w_3=0.0455,$$

$$w_4=0.0493,w_5=0.01,w_6=\mathrm{0.006...}$$

$$\underset{{\mathbf{c}}_g}{\min }\frac{1}{2}\Vert {\mathbf{c}}_g-{\mathbf{c}}_f{\Vert }^2+{\mu }_t\Vert \nabla g{|}_{\Gamma }{\Vert }_1+{\mu }_f\Vert {\nabla }^{\alpha }g{|}_{{\Gamma }^c}{\Vert }_1$$

$$s.t.g=\Phi {\mathbf{c}}_g.$$

$$\Gamma_{k}=\texttt{round}\bigl{(}\frac{1}{k-1}\sum_{i=1}^{k-1}\Gamma_{i}\bigr{)},$$

$$\Gamma_{k}=\Bigg{(}{\sum_{i=q}^{k-1}w_{i}\Gamma_{i}>t}\Bigg{)},$$

$$\begin{array}{rl}\text{P}:& {\mathbf{c}}_f={\mathrm{argmin}}_{\mathbf{c}}\Vert \mathcal{S}\mathcal{F}\Phi \mathbf{c}-\hat{f}{\Vert }_2^2,\\ \text{C}:& \underset{{\mathbf{c}}_g}{\min }\frac{1}{2}\Vert {\mathbf{c}}_g-{\mathbf{c}}_f{\Vert }^2+{\mu }_t\Vert \nabla g{|}_{\Gamma }{\Vert }_1+{\mu }_f\Vert {\nabla }^{\alpha }g{|}_{{\Gamma }^c}{\Vert }_1\\ & s.t.g=\Phi {\mathbf{c}}_g.\end{array}$$

$$\text{P}:$$

$$\text{C}:$$

$$\mathbf{prox}_{\lambda\mathcal{R}}({\boldsymbol{v}})=\mathop{\rm argmin}_{{\boldsymbol{c}}}\biggl{\{}\mathcal{R}({\boldsymbol{c}})+\frac{1}{2\lambda}\|{\boldsymbol{c}}-{\boldsymbol{v}}\|_{2}^{2}\biggr{\}}.$$

$${\mathbf{prox}}_{\lambda \mathcal{R}}(\mathbf{v})=(A+1/\lambda I{)}^{-1}\mathbf{v}/\lambda .$$

$$({\mathbf{prox}}_{\lambda \mathcal{R}}(\mathbf{v}){)}_i=\{\begin{array}{ll}{v}_i-\lambda & \text{mbox}if{v}_i>\lambda ,\\ v_i+\lambda & \text{mbox}if{v}_i<-\lambda ,\\ 0& |v_i|\le \lambda .\end{array}$$

$$\underset{f}{\min }\frac{1}{2}\Vert f-\tilde{f}{\Vert }^2+{\mu }_t\Vert \nabla f{|}_{\Gamma }{\Vert }_1+{\mu }_f\Vert {\nabla }^{\alpha }f{|}_{{\Gamma }^c}{\Vert }_1,$$

$$\underset{f}{\min }\frac{1}{2}\Vert f-\tilde{f}{\Vert }^2+{\mu }_t\Vert \nabla f{\Vert }_1.$$

Full text

{AMS}

35R11, 65K10, 65F22, 90C25

1 Introduction

Image reconstruction from Fourier measurements has been a fundamental problem in various applications, such as magnetic resonance imaging (MRI)[41, 42, 51, 36, 33, 4], ultrasound imaging [59, 19] and synthetic radar imaging [3, 18]. The reconstruction methods in the literature can be roughly classified into two categories: the discrete models and the continuous ones. The discrete models view the underlying image as a discrete vector with certain fixed resolution and usually obtain its approximation through solving a discrete optimization problem consisting of a fidelity term and a regularity term. There are various regularity terms used in the literature, such as total variation [53, 47], total generalized variation [11, 38], and total fractional variation [16, 17, 63] etc.. Moreover, many efficient algorithms such as alternating direction method of multipliers (ADMM), primal-dual methods have been proposed [9, 31, 49, 66, 60, 7, 48, 8] to solve corresponding optimization problems. On the other hand, the continuous models consider the underlying image as a piece-wise smooth function and recover the image from a function approximation point of view. One of its advantages is the flexibility in setting resolution and it has been successfully employed in the reconstruction of super-resolution images [13, 24]. It has also been shown to have superior performance in generalized/infinite-dimensional compressive sensing reconstruction [1, 2].

Image reconstruction has usually been formulated as an optimization problem that minimizes an energy functional in the following form:

[TABLE]

where $\mathcal{L}(x)$ is a fidelity term depending on an empirical estimation of the distribution of noise, and $\mathcal{R}(x)$ is a regularity term with a prior estimation of the structure of the underlying image. For instance, $\mathcal{L}(x)$ is often a least squares term while the noise is assumed to be Gaussian, and other formulations could also be found in [39]. On the other hand, a widely used regularity term $\mathcal{R}(x)$ is the $l_{1}$ type constraint incorporating certain sparsity prior knowledge. Such sparsity may come from edge estimation [40, 35, 12], wavelet transformation [26, 55, 64, 20], different orders of total variation [11, 63, 15, 21] etc..

Most of existing methods solve the above optimization problem with various fidelity terms and regularity terms in either continuous or discrete settings. We will refer to them as one stage methods in this paper. We will leverage both discrete and continuous models to develop a two-stage Projection Correction Modeling (PCM), in which the first stage (P-stage) employs a continuous model and the second stage (C-stage) imposes a discrete regularization/penalty term on the model. In particular, for image reconstruction from Fourier measurements, we will later show in numerical experiments that the proposed two stage PCM has a superior performance comparing with other popular one stage methods.

We demonstrate the idea of PCM with a reconstruction problem. The given data are some finite uniform or non-uniform Fourier measurements. The goal is to reconstruct the underlying image from these measurements. At the P-stage, we consider the underlying image as a function $f$ in a processing domain that is usually a Hilbert space spanned by some basis such as polynomials or wavelets. We will find an “optimal” approximation (projection) in the processing domain by minimizing a certain data fidelity term. In other words, the P-stage projects the Fourier (k-space) measurements into another processing domain that has an accurate representation of the underlying function. We also point out that the P-stage could also be viewed as a dimension reduction step, since we often choose a projection on a much lower dimension subspace. It will also help to reduce the computational time at the second stage. Due to the noise in the measurements and/or imperfect selection of the basis, the approximation in the P-stage will also contain errors. To further improve the reconstruction, we will impose a discrete regularization at the C-stage.

At the C-stage, we will find a “corrected” approximation in the same processing domain by minimizing the sum of a date fidelity term and a regularity term. The data fidelity is the difference between the corrected approximation and the approximated function obtained at the P-stage. On the other hand, we will employ a regularity term on the discrete vector that is the evaluation of the corrected approximation function on discrete grids. In particular, we will consider a hybrid regularity combining total variation and total fractional variation. The total variation (TV) regularization has very good performance in keeping edges in the reconstruction but suffers from the staircase artifact which causes oil-painted blocks. It is mainly due to the fact that TV is a local operator. On the other hand, the total fractional variation (TFV) is a recent proposed regularization term in image processing and has achieved promising results [16, 17, 63, 65, 58, 44]. Imposing TFV could reduce such artifact due to its non-local nature. However, the edges are damped in the reconstructed image based on TFV regularization. Therefore, we will adopt a hybrid regularity with TV on the “edges” and TFV on the “smooth” part.

We remark that edge detection becomes an important task in our method, since we will not know where the true edges are before we proceed. We would like to mention that there is another research pipeline in multi-modality image reconstruction where computed tomography (CT) and MRI scanning would be run at the same time and CT images are used to enhance the performance of MRI as well as reducing the processing time [40, 10, 22, 37]. The set of edges will be given since the CT is much faster than MRI. However, we consider a more challenging task in this paper where the edges are not known and will be reconstructed recursively in our algorithm.

We briefly summarize our contributions below:

1. (1)

We propose a novel two-stage PCM method of image reconstruction that leverages advantages of both discrete and continuous models. 2. (2)

We employ a precise and efficient edge based variational constraint in the regularity term, while the edge is determined through combining techniques of image morphology and thresholding strategy. 3. (3)

We develop an efficient proximal algorithm of solving the proposed model.

The rest of the paper is organized as follows. We introduce the proposed PCM framework for image reconstruction from Fourier measurements in Section 2. We develop in Section 3 an efficient proximal algorithm for the general PCM model and further employ it to derive an algorithm for a specific PCM-TV-TFV model. Numerical experiments and comparisons are presented in Section 4. We finally make some conclusion remarks in Section 5.

2 Projection Correction Modeling

We will present the proposed PCM framework for image reconstruction from Fourier measurements in this section. To this end, we first give a brief introduction of the image reconstruction problem.

Suppose the underlying image is a real-valued function $f:\Omega\to{\mathbb{R}}$, with $\Omega\subseteq{\mathbb{R}}^{2}$. We are given its Fourier data $\hat{f}$ in the following form:

[TABLE]

• $\bullet$

$\mathcal{S}$ is the sampling operator (might be uniform or non-uniform),

• $\bullet$

$\mathcal{F}$ is the continuous Fourier transform as

[TABLE]

• $\bullet$

$\epsilon\in{\mathbb{C}}^{m}$ is random noise.

Our goal is to recover the underlying image $f$ from the given Fourier data.

The challenges come from the non-uniformness of $\mathcal{S}$ as well as the appearance of the noise $\epsilon$. It is well received that most images are piece-wise smooth with potential jumps around the edges. The inverse Fourier transform will not work directly here. The non-uniformness of the samples in the frequency domain (k-space) will make the inverse process unstable and might bring extra approximation errors when the sampled data contain noises. Moreover, the Fourier basis is amenable to Gibbs oscillations artifacts in representing a piece-wise smooth function. To overcome such challenges, we will propose a two-stage PCM framework below.

We point out that the PCM framework can be viewed as a generalization of the classical multi-scale Galerkin method in finite element analysis. In contrast to the classical triangle or polyhedron mesh segmentation in a single scale, the region has elements in different scales. The advantage of doing so is to increase the numerical stability of the solver. More details can be found in chapter 13 of [29].

We present the general idea of such a PCM framework in Figure 1. We shall next introduce the P-stage and the C-stage in PCM framework with more details.

2.1 The P-stage

Suppose $\{\phi_{1},\phi_{2},...,\phi_{n}\}$ is a basis of a subspace of the processing domain and we will find an “optimal” approximation

[TABLE]

where its coefficients ${\boldsymbol{c}}_{f}$ is obtained through solving the following least squares problem

[TABLE]

where $\Phi{\boldsymbol{c}}=\sum_{j=1}^{n}c_{j}\phi_{j}$. The above least squares problem could be solved efficiently by conjugate gradient solvers.

The performance of the P-stage depends on the selection of the basis $\{\phi_{i}\}_{i=1}^{n}$ according to the prior knowledge of the property of the image. We point out that the above framework has also been used in [55, 28] for approximating the inverse frame operator with admissible frames. Analysis of the approximation error has also been discussed there. In particular, Fourier basis has been used to obtain a stable and efficient numerical approximation of a smooth function from its non-uniform Fourier measurements. However, Fourier basis might not be appropriate to represent a piece-wise smooth function. Instead, we will consider wavelets including classical ones such as Haar, Daubechies wavelet or more recent ones such as curvelet [56, 43], shearlet [27, 34] etc., which have a more accurate representation for piece-wise smooth functions. In this regard, it could also viewed as a generalization of the admissible frame method in [55, 28]. In this paper, 1D and 2D Haar wavelets are used in the numerical experiments, but the general PCM framework could also work with other wavelets.

We remark that even with a reasonable selection of basis $\{\phi_{i}\}_{i=1}^{n}$, ${\boldsymbol{c}}_{f}$ from Eq. 2 might not be accurate due to the noise in the data. We proceed to improve it in the correction stage to alleviate those effects.

2.2 The C-stage

We will find a “corrected” approximation

[TABLE]

in the same processing domain through solving the following regularization optimization problem

[TABLE]

where $\mathcal{R}(c_{g})$ imposes prior knowledge based on the understanding of the processing domain. The idea is to find a new set of coefficients that are close to ${\boldsymbol{c}}_{f}$ found in P-stage such that the new approximation satisfies some regularities. In terms of optimization, it provides the proximal guidance during the algorithm implementation.

2.2.1 TV-TFV Regularity

We next discuss the choice of the regularity term $\mathcal{R}({\boldsymbol{c}}_{g})$. One popular choice in the literature uses $\|c_{g}\|_{1}$ by assuming the underlying function has a sparse representation in the processing domain. However, it might not be the best choice for the problem of image reconstruction from Fourier measurements. It does not incorporate the property of images that are piece-wise smooth and contain many textures/features. Moreover, it might have the bias issue in statistics literature that involves the modeling error brought by imperfect selection of basis $\{\phi_{i}\}_{i=1}^{n}$. That is, a regularity term on the coefficients $\|c_{g}\|_{1}$ might not be the best choice of alleviating the bias issue. Instead, we would impose regularization in image/function values on a discrete grid through some prior knowledge about the underlying images such as the piece-wise smoothness, textures/features. In particular, we shall employ a regularity that combines total variation (TV) and total fractional variation (TFV). To this end, we first review some definitions and notations of fractional order derivatives.

We point out that there are several definitions of fractional order derivatives, such as Riemann-Liouville (RL), Grunwald-Letnikov (GL), Caputo etc. We will employ the RL definition [45, 32] in this paper. The left, right and central RL derivatives of order $\alpha\in(n-1,n)$, $n\in\mathbb{N}$ for a function $f(x)$ supported on an interval $[a,b]$, are defined by

[TABLE]

[TABLE]

and

[TABLE]

where $\Gamma(\alpha)$ is the Euler’s Gamma function

[TABLE]

We next introduce the discretization of the RL fractional derivative. In the simplicity of presentation, we will show the 1D discretization over an interval $[a,b]$. The discretization over a 2D regular domain will be a direct generalization of this procedure along horizontal and vertical directions. We consider the following $n$ equidistant nodes on $[a,b]$:

[TABLE]

Let $L^{(\alpha)}$ and $R^{(\alpha)}$ be the matrix approximations of the left and right-sided Riemann-Liouville $\alpha-$order derivative operator ${}_{a}D_{x}^{\alpha}$ and ${}_{x}D_{b}^{\alpha}$ accordingly. With Dirichlet boundary conditions that $f(a)=f(b)=0$, it follows [50] that $L^{(\alpha)}$ and $R^{(\alpha)}$ are two triangular strip matrices with the following structure:

[TABLE]

and

[TABLE]

where $w_{0}^{\alpha}=1$, $w_{j}^{\alpha}=(-1)^{j}\binom{\alpha}{j}$. These coefficients then can be constructed iteratively:

[TABLE]

Furthermore, when $\alpha\in(1,2)$, the matrix approximation of the central RL derivative $C^{(\alpha)}$ will then become

[TABLE]

Let $U\in{\mathbb{R}}^{N}$ be the given discretized image $u$ under an ordinary xy-coordinate with pixel values lexicographically ordered in a column vector. For simplicity, we assume that the original image $u$ is in square size with $n=\sqrt{N}$ rows and columns. Let $\otimes$ denote the Kronecker product. It follows from applying the central RL derivative $C^{(\alpha)}$ to the images along the $x$-direction that

[TABLE]

where $I_{n}\in{\mathbb{R}}^{n\times n}$ is the identity matrix. Similarly, along the $y-$direction we will have

[TABLE]

The procedure will be the same for other choices of $L_{n}^{(\alpha)}$ and $R_{n}^{(\alpha)}$. The $l_{1}$ norm regularization over $u_{x}^{(\alpha)}$ and $u_{y}^{(\alpha)}$ leads to the total fractional variation models [63].

We remark that even though $C_{n}^{(\alpha)}$ is a dense matrix, the matrix $I_{n}\otimes C^{(\alpha)}_{n}$ will be sparse. Furthermore, we observe that $w_{j}^{\alpha}$ decays very fast. For example, take $\alpha=1.3$, the first few $w_{j}$’s

[TABLE]

To further enhance the sparsity of the operation matrix, one can truncate $w_{j}$’s at a certain level to improve the efficiency of the program.

We are now ready to introduce the C-stage in the following form:

[TABLE]

where $\Gamma$ is an open domain centered around the edge part of the image, $\nabla^{\alpha}$ is the discretized $\alpha$-order fractional differential operator.

We point out that an important question in the above model Eq. 3 is to select the $\Gamma$ set. It denotes an open region centered around the edges rather than the edges themselves. The reason is that the set of edges has zero measure in ${\mathbb{R}}^{2}$ and it will not be very meaningful to consider the 2D total (fractional) variation on that.

2.2.2 Construction of the $\Gamma$ set

We shall next discuss the construction of the $\Gamma$ set. We would first detect the edges of the image reconstructed from P-stage. In order to obtain accurate edges, our approach is to find a rough estimation at the initial step, and then update it iteratively with the reconstructed image during the implementation of the reconstruction algorithm. Specifically, we would use the first few iterations of the total variation model [53] as a warm start to obtain an initial estimate of the edges. Meanwhile, we will take an early termination in such an iterative method to avoid staircase artifacts. We would then apply these filter edge detectors such as Sobel, Canny and Prewitt filters [14, 57, 23] on the initial result to improve the accuracy of edges.

Once the edges are obtained, an initial $\Gamma$ set can be easily constructed via morphologic dilation [54], see Figure 2 for a visual illustration.

We would then update the $\Gamma$ set through a few iterations of the previous results. We next discuss the update of the $\Gamma$ set in more details. We denote the $\Gamma$ set obtained at the $i$th iteration by $\Gamma_{i}$. When estimating $\Gamma_{k}$, we would use the previous $\Gamma_{i}$’s:

[TABLE]

where $\texttt{round}(\cdot)$ is the standard round function that returns the closest integer. This is to ensure that $\Gamma_{k}$ will be a binary matrix. A more general thresholding scheme with a given distribution $\mathcal{D}$ at a certain confidence level $t\in[0,1)$ will be

[TABLE]

for some $1\leq q\leq k-1,q\in N_{+}$, $w_{i}\sim\mathcal{D}$, $w_{q}\leq w_{q+1}\leq\cdots\leq w_{k-1}$, $\sum_{i=q}^{k-1}w_{i}=1$.

We could further reduce the computational cost by thresholding on the number of iterations. We point out that the reconstructed image will become more accurate during the iterations of the algorithm and there will be little variance in the detected edges after certain iterations. That is, it is a reasonable to use a small number of iterations in edge detection.

We summarize the procedure for reconstructing the initial $\Gamma$ set as follows:

Warm up and edge detection.

1. (a)

Run TV model for a few iterations (in our experiment $3\sim 5$ will be enough). This depends on a rough estimate of the noise level. A quick noise variance estimation can be found at [52]. 2. (b)

Implement filer based edge detector such as Canny, Sobel filters to detect the edge set. 2. 2.

Obtain $\Gamma$ through image dilation.

The $\Gamma$ set is updated as the underlying image is iteratively updated.

2.3 PCM-TV-TFV in Image Denoising

We shall summarize the P-stage and the C-stage to present the complete process of the two stage PCM-TV-TFV model as follows:

[TABLE]

with $\Gamma$ is an open domain centered around the edges.

We will demonstrate the advantages of the above proposed two stage model through numerical comparisons with other popular models later.

3 Algorithms

In this section we shall develop a proximal algorithm scheme for solving the general PCM optimization problem. Moreover, we will also introduce a split Bregman scheme for solving the TV-TFV regularity problem. We would then combine them to derive a specific algorithm for the PCM-TV-TFV model.

3.1 General PCM Model Solver

We shall first consider the following general PCM model with a general regularity term $\mathcal{R}$:

[TABLE]

We will introduce a general proximal algorithm for solving the above PCM optimization problem. We remark that the proximal algorithms refer to a class of algorithms that are widely used in modern convex optimization literature, for a comprehensive survey, see [49]. To this end, we first review the definition of the proximal operator : for a convex function $\mathcal{R}$, the proximal operator $\mathbf{prox}$ [49] is defined as

[TABLE]

We point out that computing the proximal operator is equivalent to solve a trust region problem [25, 49]. Moreover, we could derive closed forms of such proximal operators for many popular functions. For example, if $\mathcal{R}({\boldsymbol{c}})={\boldsymbol{c}}^{T}A{\boldsymbol{c}}$ for some symmetric positive semidefinite matrix $A$, then

[TABLE]

If $\mathcal{R}({\boldsymbol{c}})=\|{\boldsymbol{c}}\|_{1}$, we have

[TABLE]

which is the component-wise soft thresholding operator.

Consequently, we present an efficient proximal algorithm in Algorithm 1 for solving the general PCM optimization problem, based on Nesterov’s accelerated gradient method as well as the FISTA algorithm [46, 6, 9].

We observe that Algorithm 1 involves a linear program and an iterative proximal operator evaluation problem. Both of them can be computed very efficiently. The update step for $\lambda^{k}$ is taken from the FISTA algorithm [46, 6, 9], which will accelerate the convergence.

3.2 TV-TFV Regularity Solver

In this subsection we shall consider the following TV-TFV regularity problem:

[TABLE]

where $\Gamma$ again is an open domain centered around the “edge” set. If $\mu_{f}$ is set to be 0 and the width of $\Gamma$ is large enough, the above model will reduce to the classical TV denoising model.

[TABLE]

Similarly, if $\Gamma$ is set to be $\emptyset$, it will become the total fractional variational model,

[TABLE]

We remark that this model is used in the image denoising problem. The work-flow of image denoising problem with TV-TFV regularization is summarized in Figure 3.

We shall introduce an algorithm of solving the above TV-TFV regularization problem. We will start with a warm up procedure to get an estimate of the edge regions. We would then employ the split Bregman method [31, 62] to derive the algorithm. We present it in Algorithm 2, where shrink is the soft shrinkage function defined by

[TABLE]

We will demonstrate the proposed TV-TFV regularization at the C-stage has a superior performance by numerical comparisons with the TV denoising model and the TFV denoising model in Section 4.

3.3 PCM-TV-TFV Solver

We now present an algorithm of solving the PCM-TV-TFV model Eq. 5 for image reconstruction in Algorithm 3. It follows from a direct application of ADMM [9] and an accelerated ADMM algorithm [30].

4 Numerical Experiments

In this section we demonstrate the superior performance of our proposed two stage PCM framework for image reconstruction from Fourier measurements. In particular, we will use numerical experiments to show that

1. (1)

The projection step itself can achieve accurate recovery in the case without noise and bias error; 2. (2)

The projection-correction with TV regularity (PCM-TV) has a better performance than many of the state-of-the-art continuous models; 3. (3)

The TV-TFV regularity leads to better results in image denoising; 4. (4)

The projection-correction with TV-TFV regularity (PCM-TV-TFV) model further improve the results of the PCM-TV.

We will focus on the function reconstruction from non-uniform measurements. In particular, we consider the jittered sampling in the frequency domain. That is, we assume the sampling is taken at the following frequencies:

[TABLE]

We display an example of the jittered sampling in Fig. 4.

4.1 Accurate Recovery of the P-Stage

We will use a simple example to demonstrate the accurate recovery of the P-stage when there is no noise and no bias error in the model. For a selected processing domain $\mathcal{H}_{n}={\text{span}}\{\phi_{i}\}_{i=1}^{n}$, we say it has no bias error if the underlying function $f\in\mathcal{H}_{n}$.

We will use the classical Haar wavelets in the processing domain. In the 1D case, the mother wavelet $\psi(x)$ is

[TABLE]

and its descendants are $\psi_{n,k}=2^{n/2}\psi(2^{n}x-k)$, $k=\mathbb{Z},x\in[0,1)$. The 2D Haar wavelet is formulated by simple cross-product.

We consider the following piece-wise constant test function,

[TABLE]

One can verify that the support of this test function is a subset of that of the selected Haar wavelets. Figure 5 shows that our proposed projection stage achieves a very accurate recovery for noiseless case.

4.2 PCM-TV Model

We consider in this subsection the case where the reconstruction contains bias error. In other words, the underlying function does not lie in the finite dimensional subspace spanned by the chosen basis. We will consider the following piece-wise linear test function $f(x)$:

[TABLE]

We will use the same Haar wavelets to construct the processing domain. Figure 6 displays the bias error.

Suppose we are given $m=256$ non-uniform Fourier measurements with some added Gaussian noise $\epsilon$ at various noise levels $\sigma={\|\epsilon\|_{2}}/{\|\hat{f}\|_{\infty}}$.

We will use TV as the regularity at the C-stage of our two stage PCM method and call it PCM-TV. Moreover, we will compare the proposed two stage PCM-TV with the following popular one stage method with different regularities:

• •

$\ell^{1}$ regularization model

[TABLE]

• •

Tikhonov regularization model,

[TABLE]

• •

Single stage TV (SS-TV) model,

[TABLE]

In particular, we will compare them in terms of a few different performance measurements including structural similarity (ssim) [61], peak signal-to-noise ratio (psnr), signal-to-noise ratio (snr) and relative error (rela_err) defined as follows:

[TABLE]

where $u\in{\mathbb{R}}^{d_{1}\times d_{2}}$ is the true image and $\hat{u}\in{\mathbb{R}}^{d_{1}\times d_{2}}$ is the reconstruction. We present the numerical results in Table 1.

We can observe from Table 1 that our proposed PCM-TV has a better performance than all the other three one-stage methods.

We next consider the 2D case. We will consider the 2D function $f$ with a randomly chosen square support of $[0.25,0.5]^{2}\cup[0.61,0.83]^{2}$ under where the entire image region is defined as $[0,1]^{2}$. We display it in Figure 7. The 2D Haar wavelet used here is a direct product of 1D Haar wavelet leading to $\mathcal{H}_{n}^{2}:=\mathcal{H}_{n}\otimes\mathcal{H}_{n}$.

In particular, we will use a piecewise constant test function whose support does not lie in $\mathcal{H}_{n}^{2}$. It consists of four squares and the white ones do not lie in $\mathcal{H}_{n}$, which will cause the bias error in reconstruction.

We will compare the proposed PCM-TV model with the single stage TV (SS-TV) method. We point out that TV regularization usually yield better results than $\ell^{1}$ regularization and $\ell^{2}$ regularization in 2D imaging problem. We present the corresponding numerical results in Table 2.

We could observe from the above numerical results that the two stage PCM-TV has a better performance than the single stage TV regularization method. We would next explain the differences between them. We point out that the two stage PCM-TV model is equivalent to the following optimization problem,

[TABLE]

where $A=\mathcal{S}\mathcal{F}\Phi$ and $\mathcal{R}$ is the TV operator. The corresponding first order optimality condition (from Fermat’s rule [5]) implies that

[TABLE]

where $\partial\mathcal{R}$ denotes the sub-differential of $\mathcal{R}$. On the other hand, for the single stage SS-TV model

[TABLE]

the first order condition implies that

[TABLE]

Compared with Eq. 7, the descent direction in Eq. 8 is distorted by the factor $(A^{*}A)^{-1}$. It might not only cause extra computational cost, but also bring numerical instability and extra errors in computing the inverse of $A^{*}A$.

4.3 Performance of TV-TFV regularity

In this subsection we demonstrate the advantages of the proposed TV-TFV regularity. We consider the following general image denoising problem

[TABLE]

where $\hat{u}$ is the given noisy image. We will compare these three different regularities in the above model: TV, TFV, and TV-TFV.

We display the noisy image and the reconstructed images from these three denoising methods in Fig. 8 111 Image retrieved from http://radiopaedia.org/ by Frank Gaillard.. To better understand the performance, we zoom in the selected part of the image and display them in Fig. 9.

We observe that the anisotropic TV suffers from the staircase artifact due to the fact that the TV is local operator. On the other hand, the reconstruction with TFV regularity has blurry effect on the edges. This is not surprising because the TFV is a non-local method and it is less edge sensitive than TV. Instead, the TV-TFV regularity avoids such artifacts and has a better reconstruction of both the edges and the overall image.

We also present the numerical results of different performance measurements in Table 3. The TV-TFV regularity shows better results in such measurements as well.

4.4 PCM-TV-TFV vs. PCM-TV

Finally, we will combine the projection stage and the correction stage with TV-TFV regularity instead of TV regularity to further improve the performance.

We display the ground truth image and the reconstructed images from the inverse Fourier method, the PCM-TV method and the PCM-TV-TFV method in Fig. 10.

We point out that the noisy image in Fig. 10 is obtained directly by inverse Fourier transform and we can see that the noise level is quite high in this case. Both the PCM-TV and the PCM-TV-TFV are able to produce more reasonable visual results. To see a deep comparison, we zoom in the red square part of Fig. 10 and present the approximation errors in Fig. 11.

We can see that the PCM-TV-TFV has a much less error than the PCM-TV model for the surface reconstruction.

5 Conclusions

In this paper we propose a general two stage Projection Correction Modeling framework for image reconstruction from Fourier measurements. The projection step alleviates the instability from the non-uniformness of the Fourier measurements and the correction step further reduces the noise and the bias effects in the previous step. A precise edge guided TV-TFV regularity shows its own advantages over the models with single TV or TFV regularity. The numerical experiments demonstrate that such combination enhances the reconstruction and reduces the drawbacks of the TV and TFV themselves. Furthermore, we also show that the proposed PCM-TV-TFV has a superior performance even when the measurements have considerably heavy noise.

We remark that we only use wavelets to demonstrate the advantages of the proposed two stage PCM method in the numerical experiments. However, other basis such as shearlets, curvelets or adaptive wavelets could also be incorporated to further improve the results. Moreover, the proposed two stage framework for Fourier measurements could also be extended to other linear measurements such as Radon transform etc..

Acknowledgment

The authors would like to thank the discussion with Yiqiu Dong.