Abstract

The analysis model has been previously exploited as an alternative to the classical sparse synthesis model for designing image reconstruction methods. Applying a suitable analysis operator on the image of interest yields a cosparse outcome which enables us to reconstruct the image from undersampled data. In this work, we introduce additional prior in the analysis context and theoretically study the uniqueness issues in terms of analysis operators in general position and the specific 2D finite difference operator. We establish bounds on the minimum measurement numbers which are lower than those in cases without using analysis model prior. Based on the idea of iterative cosupport detection (ICD), we develop a novel image reconstruction model and an effective algorithm, achieving significantly better reconstruction performance. Simulation results on synthetic and practical magnetic resonance (MR) images are also shown to illustrate our theoretical claims.

1. Introduction

Sparse sampling theory [13] plays a key role in a broad spectrum of techniques involved in signal and image processing over the past decade. It states that an unknown signal can be recovered from a small number of random linear measurements given that the signal is sparse. Consider a sparse signal or image, vectorized as , which has very few nonzero elements in the sense that , where is the count of the nonzeros in . We expect to reconstruct by solving the following -norm minimization problem,where denotes the measurement matrix that produces the measurement vector . However, this is an NP-hard problem which prompts us to look for alternatives to solve it in an approximate fashion. A popular and effective way is to rewrite (1) as the basis pursuit problem:As a matter of fact, this is not the general case since signals or images do not exhibit sparsity directly but have sparse representations in specific transform domains. Broadly speaking, there are two data models for describing signals. The first one is the sparse synthesis model [4, 5], wherein is assumed to admit sparse representation in a fixed dictionary . Put differently, can be viewed as a linear combination of very few dominant atoms from . Thus, the reconstruction process (2) is reformulated asAs we know, a tremendous surge of effort has been devoted to studying the sparse synthesis model, and much progress has been made ranging from theoretical foundations [6, 7] to appealing applications, including denoising [8], inpainting [9], and more. Additionally, a long series of algorithms [1013] together with provable guarantees [14, 15] have been put forward.

While the synthesis model has gained widespread attention, a similar alternative was proposed modeling signals from an analysis perspective [1620]. Mathematically, given an analysis operator , a signal belonging to the analysis model is supposed to admit a sufficiently sparse analysis representation . In particular, we are interested in cases where , so that various information within can be captured [17]. The commonly used analysis operators include the finite difference [21], overcomplete wavelet transforms [22], and the curvelet transform [23]. Just in contrast to the “sparsity” in the synthesis setting, the analysis model concentrates on the zero-valued coefficients. A fundamental notion measuring the quantity of the zeros is defined as cosparsityThe -cosparse analyzed vector is assumed to have unknown locations of the zeros, referred to as the cosupportRecent studies have theoretically shown that the analysis model has its own advantage over the synthesis one [17]. Moreover, adopting the analysis model leads to a collection of successful applications, such as denoising [23, 24], inpainting [25], and medical imaging [26, 27]. Models mentioned above simply exploit the (co)sparsity prior that is implicit in signals or images. However, (co)sparsity alone is insufficient for making reasonable inferences from these models, and the minimum measurement requirement for reconstruction is improved limitedly. Obviously, the trade-off between (co)sparsity and measurements can be further improved by imposing a priori knowledge. A common way to incorporate prior knowledge is through the use of the signal support, which has been extensively studied in the synthesis context [2830]. However, relatively little research has been devoted to imposing prior knowledge in the analysis sense. In this work, we wish to investigate the consequences of exploiting analysis model prior. We concentrate on image reconstruction given the a priori cosupport knowledge. First, we formulate the reconstruction problem in terms of the analysis operator in general position, which means that every set of rows from are linearly independent [17]. In such a case, we derive the minimum number of measurements for uniquely determining a cosparse solution to the inverse problem . The resulting number is smaller than that of the standard analysis model without using prior cosupport. Second, we dive into details of the model associated with the 2D finite difference operator whose rows show significant linear dependencies [17]. We also provide an improved minimum measurement requirement for guaranteeing the uniqueness given the prior cosupport. Having a theoretical foundation for the analysis uniqueness properties, we develop a novel image reconstruction method based on the idea of iterative cosupport detection (ICD). This two-stage scheme proceeds by alternatingly calling its two key components: image reconstruction and cosupport detection. At each iteration, the image is reconstructed using the cosupport knowledge extracted from the previous iteration. After the acquiring of the image estimate, one can identify even better cosupport to be used in the next iteration. Moreover, when performing the cosupport detection, a multidirectional finite difference is considered. Therefore, the detection and the reconstruction parts work together enabling us to gradually obtain reliable cosupport and a reasonable image estimate. Simulation results on synthetic and practical magnetic resonance (MR) images demonstrate the effectiveness and show considerable improvement of the proposed method compared with other regularization based methods. This consequently indicates that, through the use of the analysis model prior, we can achieve a given reconstruction quality with fewer measurements or alternatively obtain a better reconstruction under the same measurement requirement.

The remainder of this paper is organized as follows. In Section 2, a detailed description of our proposed method is given from two aspects: analysis operators in general position and the 2D finite difference operator. Uniqueness issues are also explored. Section 3 describes the methods used for comparison and rules for image quality assessment. Section 4 presents simulation results and discussion which validate our theoretical claims. Finally in Section 5, conclusions are summarized.

2. Method

Analysis model prior has been successfully used for many signal processing tasks but has been done with little theoretical justification. In this section, we focus on the analysis-based reconstruction given that the cosupport is known a priori. The number of measurements required for guaranteeing the unique reconstruction is proved to be essentially reduced. We first consider the analysis operators in general position.

2.1. Analysis Operators in General Position
2.1.1. The Cosupport Is Exactly Known

Consider -cosparse image whose cosupport is and a redundant analysis operator in general position; namely, every set of rows from are linearly independent [17]. The cosparse representation analysis model where belongs is related to the recently proposed union-of-subspaces model [3133]:which is a union of all the possible subspaces of dimension . Here, denotes the null-space of the analysis operator indexed by . is the cardinality of . The rows associated with , , define the analysis subspace. Remark that removing rows from for which leaves the subspace unchanged.

Consider the aforementioned noiseless linear inverse problem:where and are assumed to be mutually independent. Provided that the cosupport is exactly known, the cosparse is supposed to satisfy the linear system:The fact that (8) identifies a unique is equivalent to the requirement:This requirement indicates that the minimum number of measurements is:

2.1.2. The Cosupport Is Unknown

Actually, the case we do care about is that we only know the cosparsity level , whereas the cosupport is undetermined. The uniqueness issue in this sense has been explored in [17].

Lemma 1. Let , , be a union of -cosparse analysis subspaces induced by the analysis operator . Then the linear system admits a unique -cosparse solution if and only if for any where

Clearly, for any , we have and . Therefore, we can conclude that the minimum number in terms of the cosparsity level is

2.1.3. The Cosupport Is Inaccurately Known

More generally, we are considerably interested in the case where the cosupport is inaccurately known. In other words, a superset of the true cosupport is available, which comes up in many applications. Even in the absence of available prior knowledge, one can still extract useful information from the current solution and use it subsequently.

We still assume is -cosparse with cosupport . is the prior cosupport containing small errors denoted by . Then, the true cosupport of can be expressed as

We state the condition which enables us to guarantee the uniqueness of the linear inverse problem (7), namely, the minimum measurement requirement.

Proposition 2. Assume that , , is the -cosparse image to be reconstructed. is the prior cosupport with small error . Given and , the minimum number of measurements for identifying a unique solution to the linear system is

Proof. As we have seen, the prior cosupport is not exactly consistent with the true cosupport of ; that is, there exists small error in . It means that we are supposed to consider the invertibility of over the direct sum of any two subspaces . Assume any two -cosparse images , with corresponding cosupports and , respectively. Note that and let . Then, we have . Consider any from as defined in (12). The cosupport of is obtained aswhich yieldsThus, we haveConsequently, we conclude that the minimum number of measurementswhich completes the proof.

Proposition 2 states that the measurement number required for uniquely identifying the solution to (7) tends to be smaller in the sense that the cosupport is known a priori. As approaches , is equivalent to that of the case without prior knowledge; namely, . When the prior cosupport selected is consistent with the ground-truth, the minimum measurement number is achieved as .

2.2. Specific 2D Finite Difference Analysis Operator

In this subsection, we would like to investigate the family of finite difference operators on graphs, (significantly related to total variation norm minimization [21]), which has been proven successful for sparse signal and image recovery [3436]. We explore uniqueness issues in the sense that cosupport knowledge is known a priori. Due to the fact that this class of analysis operators exhibit strong linear dependencies [17], the theoretical results derived above cannot be applied directly. Our analysis is based upon the work [17]. For ease of notation, we will drop the subscript DIF and simply use to represent hereafter when there is no ambiguity. Assume that is defined on a 2D nonoriented graph with vertices and edges . An edge connecting two vertices and can be regarded as a finite difference. is a subset of and the set of vertices connected by at least one edge in is denoted by which is composed of connected components. A connected component is a collection of vertices connected to one another by a walk through vertices within . Thus, the dimension of is given bywhere denotes the number of isolated vertices which have distinct values from all the neighbors. Suppose that the known part of the cosupport is denoted by and the number of the connected components is . Then, in terms of the cosparsity level , we present a concrete bound for , which reveals the minimum measurement number.

Proposition 3. Let be the 2D finite difference analysis operator that computes horizontal and vertical discrete derivatives of an image . The measurement matrix is assumed to be mutually independent from . For a fixed and the known cosupport which corresponds to connected components, the equation admits at most one solution with cosparsity only if

Proposition 3 reveals that the measurements required for uniquely determining the cosparse solution to the inverse problem can be reduced given that some part of the cosupport is known. Additionally, the minimum measurement number decreases monotonically as the number of the connected components increases. When , the minimum number of measurements , which is equivalent to the case that no cosupport information is available. For completeness, the proof can be found in Appendix.

2.3. Proposed Reconstruction Model and Algorithm
2.3.1. Reconstruction Model

Armed with the above-described theoretical analyses, we are now in position to solve the reconstruction problem regularized with cosupport prior. More specifically, we expect to constrain the cosparsity of within the prior cosupport :However, as previously mentioned, the -norm involved in the combinatorial minimization program is an NP-hard penalty and thus is not feasible to be solved for practical applications. An effective and widely used alternative is the -relaxation:in which the -norm enables promoting high cosparsity in the solution. Its desirable convexity facilitates various computationally tractable algorithms, and much recent progress in the theory of analysis -minimization has been made [37, 38]. In this work, we focus on the analysis -recovery which is given in an unconstrained fashion. The objective function is formulated as a linear combination of the data consistency error and a modified cosparsity-inducing penalty:where controls the influence between the fidelity term and the regularization term:Remark that the analysis operator we employed in our method is the 2D finite difference operator involving four directional components; that is, , respectively, compute discrete derivatives in vertical, horizontal, and two diagonal directions. denote the associated four-direction cosupport sets that are detected iteratively.

2.3.2. Algorithm

Having presented the reconstruction model, we now turn to the question how to effectively solve it. We propose a two-stage algorithm based on the idea of iterative cosupport detection (ICD). The proposed ICD allows extracting the reliable information of the underlying solution and enables achieving a reasonable reconstruction. This two-stage scheme proceeds by alternatingly calling its two components: image reconstruction and cosupport detection. In the reconstruction step, we solve a truncated -minimization problem via conjugate gradient method using the cosupport knowledge obtained from the previous iteration. ICD will terminate if the approximate solution is accurate enough. Otherwise, cosupport detection will be performed in light of this inaccurate reconstruction, thereby yielding a better to be used in the next iteration. Consequently, the detection and the reconstruction parts of the proposed ICD work together enabling us to gradually obtain reliable cosupport and a reasonable image estimate.

It is important to note that the proposed ICD requires reliable cosupport detection. However in most cases, it is hard to completely avoid false detections at each iteration. To this end, we address ways to identify the index sets . Actually, the analyzed vectors in each direction are not strictly cosparse but exhibit a strong decay. A conceptually simple detection strategy is to obtain the cosupport in a truncated manner. Without loss of generality, we assume as the four-direction index sets of after sorting in an ascending order. Then, the cosupport set in the th direction is created aswhere the truncated length is fixed and specified in advance. Empirically, a relatively stable reconstruction can be achieved provided that is suitably scaled in a certain range. This strategy is easy to implement but has one drawback; that is, the assumption that the true cosupport is included in the inaccurate prior cosupport cannot be well satisfied when is fixed. To alleviate this drawback, we propose another detection strategy which is more adaptive and effective. The cosupport sets are iteratively learned based on thresholding [28]:where represents the th element of . The threshold is set asHere, is the threshold parameter selected as an exponential function of the iteration number :where is a positive integer. For , we have , indicating that no prior cosupport knowledge is exploited in the first iteration. In the sequel, the threshold decreases with the increase of the iteration so that the cosupport size reduces gradually. However, the cosupport in each direction is not strictly decreasing over the iteration, which allows the current detection to include indices within the true cosupport that are excluded from the detected cosupport in previous iterations. This leads to an attractive self-corrected capacity of the cosupport detection. The proposed algorithm, referred to as ICD, is outlined as follows.

Algorithm 4 (ICD). Consider the following:(1) Input. Consider undersampled measurements , cosupport detection parameter (or ), regularization parameter , and maximum iteration number .(2) Iteration. While the stopping criterion is not reached, do the following:(i) Image Reconstruction. Solve the minimization problem (24) for the image estimate using the cosupport information of the previous iteration .(ii) Cosupport Detection. Update the cosupport information using the detection criterion equation (27) (or (26)) based on .(iii)Consider .(3) Output. The reconstructed image .

3. Evaluation

To evaluate the performance of the proposed method, we performed simulations on both synthetic and practical MR images. Similar to prior work on CS-MRI [34, 39], we simulated the data acquisition by randomly sampling the 2D discrete Fourier transform coefficients of test images according to the patterns. Thus, the measurement matrix in (24) was defined as the undersampled Fourier transform . The number of compressive measurements was measured in terms of the percentage of total number of Fourier coefficients, namely, the sampling ratio (SR). All simulations were carried out under MatLab R2011b running on a PC with a 3.2 GHz processor and 4 GB memory. To assess the effectiveness of the proposed method, we compared it with other potential reconstruction techniques, including (i) SparseMRI [39]: the leading MR image reconstruction approach combining the wavelets and the standard total variation (TV) regularization, regardless of the effect of the wavelets, this approach can be approximately viewed as the proposed method without using cosupport information; thus the contribution of the integration of analysis a priori knowledge can be demonstrated; (ii) NLTV [40]: the well-known method based on wavelet sparsity and nonlocal TV penalty; (iii) ISD-TV [2830]: the method constraining the sparsity of wavelet coefficients over the complement of the known support, which can be regarded as incorporating prior knowledge to the sparse synthesis model; in order to make a fair comparison, we added a TV term; (iv) SDBS-TV [41]: another synthesis-based method combining support knowledge and block-sparse property in the wavelet domain. The proposed method was tested using both truncated and threshold-based strategies, respectively, named as ICD-TR and ICD-TH. We also expected further enhancement by combining our method with the wavelet penalty, named as ICD-WT, which in general solveswhere denotes the wavelet transform operator and . Remark that threshold-based strategy is used for cosupport detection in ICD-WT.

For quantitative evaluation, the reconstruction quality was measured by the relative -norm error (RLNE) which is a standard image quality metric indicating the difference between the reconstruction and the ground-truth As for practical MR images containing more fine features, the quality of the reconstruction is quantified by two other metrics. The first one is the high-frequency error norm (HFEN) [42], defined aswhere LoG is a rotationally symmetric Laplacian of Gaussian filter capturing detailed textures. The filter kernel is of size pixels and with a standard deviation of pixels, the same as that in [42]. The second one is the structural similarity (SSIM) index [43], comparing local patterns of pixel intensities between and that have been normalized for luminance and contrast:Here, and are mean intensities at the th local window of and , while and are the corresponding standard deviations. denotes the covariance and the constants , are included to avoid instability.

4. Results and Discussion

4.1. Shepp Logan Phantom

We first tested our method on an ideal example: the Shepp Logan numerical phantom of size . The undersampling scheme we employed is the 2D radial trajectory. Figure 1 shows the reconstruction results of different methods using radial lines (RL). Reconstruction parameters in (24) were set as follows: and for ICD-TR and for ICD-TH. We also tested the proposed method using and radial lines. The RLNE comparison results are presented in Table 1. The superiority of the proposed method is clearly seen from reconstructed images and error maps which were obtained by subtracting the reconstructions from the original image (shown at the same scale). The results indicate that methods without using any prior knowledge, such as SparseMRI and NLTV, cause serious artifacts. Methods imposing prior support knowledge in the synthesis sense, including ISD-TV and SDBS-TV, improve the reconstruction quality to a certain extent, but the edges are not well preserved. The proposed ICD significantly improves the reconstruction and is capable of suppressing more artifacts and preserving more details. As for different cosupport detection strategies, ICD-TH performs better than ICD-TR since the threshold-based strategy enables correcting the cosupport adaptively. Table 2 shows the number of true and false detections of four directional cosupports at each iteration using ICD-TH. It indicates that the number of false detections in each direction decreases gradually until the true cosupport is detected. However, the number of true detections is not necessarily always increasing.

4.2. Practical MR Images

While encouraged by the results on the phantom image, we conducted further simulations on more realistic images to evaluate the practical effectiveness of the proposed method. The CS measurements were generated by undersampling the Fourier coefficients of the fully sampled (ground-truth) MR images whose intensities were normalized to a maximum magnitude of . The original T1-weighted image (Brain-1) of size (courtesy of Professor N. Schuff at the UCSF School of Medicine) and the sampling mask of 30 sampling ratio are shown in Figure 2. The undersampling scheme we used is the variable density random pattern which is widely used in - plane for 3D imaging enabling removing the aliasing interference without degrading the image quality [39]. Reconstruction parameters were set as , for ICD-TR, and for ICD-TH. Figure 2 gives the reconstructed images and error maps through different methods, and Table 3 presents the HFEN and SSIM results. The capabilities of our method were also demonstrated under different sampling ratios as presented in Table 1. The reconstruction results show considerable improvement of the proposed method with respect to the subjective visual quality and quantitative indices. It has been noted from the error maps and the HFEN and SSIM results that the image reconstructed by our method has more fine features and is the closest to the ground-truth. The superior performance of ICD can be attributed to the use of the analysis model prior, namely, the cosupport knowledge. Additionally, we observe that the proposed ICD combining the wavelet penalty, ICD-WT, performs slightly better than ICD-TR and ICD-TH. However, by adding an extra penalty term, the algorithm is slowed down, which weakens the superiority of ICD-WT.

We also tested the proposed ICD on a T2-weighted image (Brain-2) which was acquired from the Siemens scanner T, SE sequence with imaging parameters: TR = 4000 ms, TE = 91 ms, slice thickness = 5.0 mm, flip angle = 120, and the field of view (FOV) = 176 × 220 mm × mm. Reconstruction parameters were set as , , and . Figure 3 presents the reconstructed images and the error maps of different methods under 30 sampling ratio, while Table 3 gives the HFEN and the SSIM values. RLNE results under different sampling percentages are also given in Table 1. It is straightforward to see that the proposed method performs the best under all sampling ratios.

4.3. Parameter Evaluation

In this subsection, we explore effects of the parameters involved in the proposed method. We begin by considering the regularization parameter , since the selection of the optimal -value is necessary. The reconstructions of Brain-1 and Brain-2 by ICD-TH are performed for different -values and under different sampling ratios. Motivated by the similar decay rates of four directional transform coefficients, we employed the same in each direction. Figure 4 displays the curves of RLNE values as a function of . The selected -values are marked with Stars. From the curves, the optimal -values under different sampling ratios are almost identically selected between 10−4~10−3.

Then, we consider effects of the cosupport detection parameters on the performance of the proposed method. As for ICD-TR, we are supposed to evaluate the truncated length . Table 4 presents the RLNE results of both Brain-1 and Brain-2 under sampling ratio with varying . It can be observed that the proper range of is image-dependent. However, the error undulates slightly when ranges from to , corresponding to 80~90% of the image dimension. The iterative process makes ICD-TR less sensitive to a few errors present in the cosupport. However with a fixed , the true cosupport is not easy to be identified. We then consider ICD-TH which is more effective and adaptive. The selection of the threshold parameter is essential in that it affects the convergence behavior of the reconstruction. The RLNE between the true and the reconstructed images at each iteration was computed and plotted for different values of . From the error-iteration curves shown in Figure 5, we see that a large leads to fast convergence but results in a relatively poor reconstruction. This is due to the fact that a large makes the threshold decrease fast so that a number of indices belonging to the true cosupport are excluded from the detected cosupport. On the other hand, a small allows a slowly changing threshold so that the cosupport can be corrected gradually. However, the convergence rate will not be satisfactory provided that is too small. Based on our simulation results, we set = 2~5 to achieve a trade-off between the convergence rate and the reconstruction quality.

The above-mentioned parameter evaluation is based on simulations on a large number of test images. The optimal ranges of the parameters we provided yield good empirical results. However, we realize that the parameters were chosen manually to minimize the reconstruction error. In practical applications, the ground-truth is not available. In these cases, one can conduct parameter selection using the effective L-curve strategy [44] or more sophisticated approaches [45, 46]. The discussion of these approaches is beyond the scope of this paper.

We also tested our method (ICD-TH) on Brain-1 and Brain-2 using radial and Cartesian patterns of and sampling ratios. From the results in Table 5, we see that different patterns may affect the reconstruction quality. However, once the pattern is fixed, our method yields the best performance.

4.4. Comparison with Two-Direction Case

In this subsection, we conducted simulations on the Shepp Logan phantom with radial lines, Brain-1 and Brain-2 under sampling ratio, using two-direction ICD. Table 6 presents the RLNE comparison results, which reveals that the two-direction ICD also performs well. However, compared with four-direction ICD, we see that the local features of the images can be better preserved by adding diagonal components.

4.5. Computational Complexity

The proposed method and the methods used for comparison except NLTV were implemented using a nonlinear conjugate gradient method with backtracking line-search. In a MatLab implementation, it took  s in average for the proposed ICD-TR and ICD-TH containing only one regularization term, while it took more than  s for other methods with an additional wavelet penalty term. However, we expect substantial reduction in the reconstruction time with code optimization and graphics processing unit implementation.

5. Conclusion

In this work, we have presented a cosparse analysis model based approach to reconstruct images from highly undersampled data using cosupport constraints. We have demonstrated that the analysis model prior can significantly improve the sparse sampling based image reconstruction. An effective iterative algorithm proceeded by alternating between image reconstruction and cosupport detection. The performance of the proposed method was evaluated through simulations on synthetic and practical MR images. The results indicate that the proposed method yields considerably better performance than methods without using prior knowledge and synthesis model based methods imposing support constraints in terms of both reconstruction accuracy and subjective visual quality.

Appendix

Proof of Proposition 3

In this part, we first prove the following two lemmas that are the key elements for the proof of Proposition 3. Note that our goal is to quantify the upper bound:To simplify the problem, we look at an equivalent quantity:

Lemma A.1. Given the known part of the cosupport with connected components, then for a fixed cosparsity level , the value is achieved when .

Proof. If , there exist extra connected components out of . Without loss of generality, we only consider the case where . Let correspond to the extra component. Thus, we have . Notice that a subgraph can be shifted horizontally or vertically unless it has vertices on all four boundaries of . Since and are disconnected, we can shift towards until they first met. Suppose and are the numbers of the vertices and the edges coincided. Clearly, . The resulting subgraph has vertices, edges, and connected components. Then, we add additional edges that are connected to , thereby yielding . From this we can infer that , , andWe further have thatConsequently, we obtain from that satisfies and . However, a smaller is achieved, which completes the proof. One can consider cases in which in the same vein.

Lemma A.1 implies that quantifying the upper bound in (A.1) is equivalent to quantifying the value . In the next lemma, we provide a lower bound for and prove that the value is achieved when all the connected components share the same cosparsity.

Lemma A.2. For a fixed , the value

Proof. We first denote the edge sets of the connected components by with (). Clearly, we have . Then, we enlarge each edge set aswhere (resp., ) denotes the edge extending downwards (resp., rightwards) from the vertex . For any , we have . According to [17, Lemma 15], the value occurs when is square, and the optimal width of is given asThen, we haveSumming the above over all yieldsFor any ,Substituting (A.10) into (A.9),Then for any , the minimum of the right hand side of (A.11) is achieved at . This reveals that all the connected components share the same cosparsity. Therefore, by denoting and substituting into (A.9), we conclude thatas required.

Proof of Proposition 3. The proof of Proposition 3 is trivial by applying Lemmas A.1 and A.2 to (A.1):which leads to the minimum measurement number: The proof is then completed.

Competing Interests

The authors declare that there are no competing interests regarding the publication of this paper.