# Block Compressed Sensing of Images Using Adaptive Granular Reconstruction

Authors: Ran Li; Hongbing Liu; Yu Zeng; Yanling Li
Journal: Advances in Multimedia
 (2016)

---

## Abstract

In the framework of block Compressed Sensing (CS), the reconstruction algorithm based on the Smoothed Projected Landweber (SPL) iteration can achieve the better rate-distortion performance with a low computational complexity, especially for using the Principle Components Analysis (PCA) to perform the adaptive hard-thresholding shrinkage. However, during learning the PCA matrix, it affects the reconstruction performance of Landweber iteration to neglect the stationary local structural characteristic of image. To solve the above problem, this paper firstly uses the Granular Computing (GrC) to decompose an image into several granules depending on the structural features of patches. Then, we perform the PCA to learn the sparse representation basis corresponding to each granule. Finally, the hard-thresholding shrinkage is employed to remove the noises in patches. The patches in granule have the stationary local structural characteristic, so that our method can effectively improve the performance of hard-thresholding shrinkage. Experimental results indicate that the reconstructed image by the proposed algorithm has better objective quality when compared with several traditional ones. The edge and texture details in the reconstructed image are better preserved, which guarantees the better visual quality. Besides, our method has still a low computational complexity of reconstruction.


---

## Body

## 1. Introduction

Nyquist frequency sampling theorem is the theoretical basis of traditional image coding (such as JPEG and JPEG 2000), which requests that the number of image transformations is the total number of pixels at least. The full transformation of image results in a large amount of calculations, but most of transformation coefficients are discarded during encoding, which causes the waste of energy consumption. Due to a high computational complexity and a low utilization rate of energy consumption, the traditional image coding is not suitable for the application requiring a light load, for example, wireless sensor network node with the limited energy consumption [1] and remote sensing imaging [2]. In addition, since a few transformation coefficients contain most of useful information, the quality of image reconstruction can be degraded severely once these important coefficients lose; thus it is also challenging for the fault-tolerant mechanism in the wireless communication [3]. Compressed Sensing (CS) [4, 5] transforms the signal with the sub-Nyquist rate, and it can still accurately recover the signal, which motivates the image CS as a new image coding scheme [6]. The CS random sampling can be regarded as the partial image transformation, and it compresses image by dimensionality reduction while transforming image. Due to the advantage of saving mostly the costs of image coding, the image CS attracts lots of academic interests.Many researchers are dedicated to improving the rate-distortion (RD) performance of image CS; the popular method explores the adaptive sparse representation model to improve the convergence performance of minimuml1-norm reconstruction [7, 8]; for example, Chen et al. [9] use the spatial correlation to establish Multiple-Hypotheses (MH) prediction model and recover the more sparse residual to improve the reconstruction quality; Becker et al. [10] proposed the NESTA algorithm based on the first-order analysis to solve the minimum Total Variation (TV) model, which ensures the speediness and robustness of the sparse decomposition; Zhang et al. [11] exploit the nonlocal self-similarity existing in the group of blocks to express the concise sparse representation model; Wu et al. [12] introduce the local autoregressive (AR) model to trace the nonstationary statistical property of image. For above-mentioned methods, the improvement of RD performance is at the expense of the high computational complexity, which results in the fact that the reconstruction time rises rapidly as the dimensionality of image increases; for example, the AR algorithm proposed by Wu et al. [12] requires 1 hour to reconstruct an image with size of 512 × 512 pixels, which makes this method lose the value of application. Compared with the above methods which ignore the computational complexity but pursue the high reconstruction quality, Gan [13] and Mun and Fowler [14] proposed the Smoothed Projected Landweber (SPL) algorithm with a low computational complexity to ensure the better RD performance. However, the SPL algorithm adopts a fixed sparse representation basis (e.g., DCT and wavelet basis) which cannot change adaptively according to the image content in the process of iteration; therefore the potential to improve the reconstruction quality is not developed completely. The Principal Component Analysis (PCA) [15] is the optimal orthogonal transformation matrix to remove spatial redundancy of image; our previous work [16] uses PCA to update continually the sparse representation basis at each SPL iteration. Adapted by the image statistical property, the PCA-based SPL algorithm guarantees the performance improvement. The defect of our previous work is that we ignore the stationary local structural characteristic of image when learning PCA matrix. This defect suppresses the performance improvement of PCA decorrelation, which results in the degradation of reconstructed image by the SPL algorithm.Aimed at the problem that PCA cannot exploit the stationary local structural characteristic of image to update the sparse representation basis at each SPL iteration, this paper proposes the adaptive reconstruction algorithm based on Granular Computing (GrC) theory [17]. The proposed method divides the image into several granules, in which any granule is a set of image patches with the similar structural characteristic. We use PCA to learn the corresponding optimal representation basis of each granule so as to obtain the efficient hard-thresholding shrinkage. Experimental results show that the proposed algorithm can improve the RD performance of image CS and achieve the better subjective visual quality.This paper will explore the adaptive granular reconstruction in block CS of images. First of all, we briefly introduce the background of block CS and GrC-based clustering, and then we describe the adaptive reconstruction algorithm in detail. Afterward, the experimental results are shown and discussed, and finally we make a conclusion on the paper.


## 2. Background

### 2.1. Image Block CS

When applied in image acquisition, CS has high time and space complexity, which may be computationally expensive for some real-world systems. Hence, most of existing image CS approaches [9, 13, 14, 16] split the images into nonoverlapping blocks, and, at the decoder, each block is recovered independently in the recovery algorithm. To the best of our knowledge, the pioneering work of block CS was proposed by Gan [13]. This work was further extended by Gan [13] and Mun and Fowler [14] to improve the performance of block CS. Block CS framework is described as follows. First of all, an Ic×Ir image x with N=IcIr pixels is divided into n small blocks with size of B×B pixels each. Let xi represent the vectorized signal of the ith block (i=1,2,…,n; n=N/B2) through raster scanning. Construct B2×B2 random Gaussian matrix, and get orthogonal matrix Ω by orthonormalizing the Gaussian matrix. Then, we randomly pick MB(≪B2) rows from Θ to produce the block measurement matrix ΦB. Finally, the corresponding observation vector yi of each block is obtained by measuring xi with ΦB as follows:(1)yi=ΦBxi.There are M CS observed values for the original image (M=n×MB), and the total measurement matrix Φ is a block diagonal matrix, in which each diagonal element is ΦB; that is,(2)Φ=ΦBΦB⋱ΦB.The block CS can compute an initial image to accelerate the reconstruction process by adopting the linear Minimum Mean Square Error (MMSE) criterion. The MMSE solution of image block can be computed as follows:(3)x^i0=RxxΦBTΦBRxxΦBT-1,i=1,…,nRxx=1ρ⋯ρ2B-1ρ1⋯ρB-22+B-12⋮⋮⋱⋮ρ2B-1ρB-22+B-12⋯1,in which ρ ranges from 0.9 to 1. We set ρ and B to be 0.95 and 32, respectively, by experience.When reconstructing image, we can construct the relationship model between the block observation vectoryi and a whole image x: (4)y=y1y2⋮yn=ΦBΦB⋱ΦBx1x2⋮xn=Φ·Ex=Θx,in which E is the elementary matrix to rearrange the column vectors block by block to a raster scanning column vector of image, and Θ=Φ·E. The minimum l1-l2 norm reconstruction model can be constructed by using (4) as follows: (5)x^=arg⁡minx⁡y-Θx22+λΨx1,in which ·1 and ·2 denote l1 and l2 norm, respectively, Ψ is the sparse representation basis, and λ is the fixed regularization factor. Equation (5) can be solved by using the SPL algorithm with a low computational complexity; in particular Ψ can be updated flexibly at each SPL iteration. Therefore, the PCA can be used to learn the adaptive sparse representation basis.


### 2.2. GrC-Based Clustering

GrC is proposed by Zadeh in 1979, and he said that the information granule exists widely in our daily life and it can be viewed as an abstraction of actual object. Granulating, organization, and causality are three important concepts of GrC [18]; that is, the granulating divides a whole of the object into several parts, the organization combines several parts into a whole by some specific operators, and the causality is used to model the relationship between the cause and its result. In the above, there is a natural relation between GrC and clustering. The training set can be regarded as the object, and each sample in training set is defined as a single granule. From the above views, we can see that the granulating and organization correspond to the clustering and its inverse, respectively, and the causality describes the internal relation between samples. Different from the traditional clustering method (e.g., K-means [19]), the GrC-based clustering chooses a variety of shapes to flexibly represent a granule, and the relationship between granules can be modeled by using the more mature algebraic system; therefore the GrC-based clustering has better generalization ability.Although the group of points in training set has irregular shape, their borders can be better distinguished by GrC-based clustering, since GrC can express the granule as various shapes, for example, hyperdiamond, hypersphere, and hyperbox. The vector form of granule is denoted asG=(c,gr), in which c is the center of granule and gr is the granularity. Denote ci as a point contained in the granule G, and the granularity of G is computed by(6)gr=maxic-cip,in which ·p is p-norm to measure the distance. The different distance measure represents the different shape of granule; for example, p=1 expresses the corresponding granule as hyperdiamond granule, p=2 expresses the corresponding granule as hypersphere, and p=∞ expresses the corresponding granule as hyperbox granule. The operators between granules include the merge operator ∨ and the decomposing operator ∧. The merge operator ∨ is to merge two smaller granules into a big one, and the decomposition operator ∧ is to decompose a big granule into two smaller ones, in which the merge operator ∨ is used to cluster the samples in the training set. The merge granules of two granules G1=(c1,r1) and G2=(c2,r2) are merged as a new granule:(7)G1∨G2=12P+Q,P-Qp,in which(8)P=c1-r1c12p,Q=c2+r2c12p,c12=c2-c1.Figure 1 illustrates the merging of two spherical granules G1=(1.0, 2.0, 2.0) and G2=(2.5, 2.5, 2.5) in the two-dimensional space, and their merging result is G1∨G2=(1.75, 2.25, 1.75).Figure 1
Merging of spherical granules in two-dimensional space.According to the vectorized representation of granule, the relation between granules inevitably involves the relation between vectors, but the partial order relation between vectors is not consistent with inclusion relation between granules. Therefore, we regard that the fuzzy inclusion relation exists between granules. When merging granules, we measure the fuzzy inclusion relation as follows:(9)μG1,G2=vG2vG1∨G2,in which v(G):G→R is a mapping from the space of granule to the space of real numbers, and this mapping v(G) uses the following formula:(10)vx=∑isin⁡πxi2.Aimed to the training setS in the clustering problem, we construct an algebra system G,μ, which is made up of the granule set G, the fuzzy inclusion relation μ(·,·), and merge operator ∨. By the fuzzy inclusion relation between granules, we control the granule merging and divide gradually S into a number of granules in order for each granule to be corresponding to a category. Figure 2 shows the GrC-based clustering results of spherical and square granules in the two-dimensional space; we can see that the GrC performs the granule merging by using the fuzzy inclusion relation and the granularity, which guarantees the better clustering performance; in particular the control of inclusion relation is very beneficial to distinguish the ownership of the training samples in the border. In a natural image signal, some relations exist between image blocks in the spatial domain. Depending on this experience, the GrC-based clustering divides automatically a whole image into several regions including blocks according to the image content. In each cluster, the structure feature of any block is very similar to others; thus it is more conducive to find the sparse representation basis adapted to the fixed data pattern.Figure 2
GrC clustering results of spherical and square granules in two-dimensional space.
(a)
Original training set (b)
Clustering results


## 2.1. Image Block CS

When applied in image acquisition, CS has high time and space complexity, which may be computationally expensive for some real-world systems. Hence, most of existing image CS approaches [9, 13, 14, 16] split the images into nonoverlapping blocks, and, at the decoder, each block is recovered independently in the recovery algorithm. To the best of our knowledge, the pioneering work of block CS was proposed by Gan [13]. This work was further extended by Gan [13] and Mun and Fowler [14] to improve the performance of block CS. Block CS framework is described as follows. First of all, an Ic×Ir image x with N=IcIr pixels is divided into n small blocks with size of B×B pixels each. Let xi represent the vectorized signal of the ith block (i=1,2,…,n; n=N/B2) through raster scanning. Construct B2×B2 random Gaussian matrix, and get orthogonal matrix Ω by orthonormalizing the Gaussian matrix. Then, we randomly pick MB(≪B2) rows from Θ to produce the block measurement matrix ΦB. Finally, the corresponding observation vector yi of each block is obtained by measuring xi with ΦB as follows:(1)yi=ΦBxi.There are M CS observed values for the original image (M=n×MB), and the total measurement matrix Φ is a block diagonal matrix, in which each diagonal element is ΦB; that is,(2)Φ=ΦBΦB⋱ΦB.The block CS can compute an initial image to accelerate the reconstruction process by adopting the linear Minimum Mean Square Error (MMSE) criterion. The MMSE solution of image block can be computed as follows:(3)x^i0=RxxΦBTΦBRxxΦBT-1,i=1,…,nRxx=1ρ⋯ρ2B-1ρ1⋯ρB-22+B-12⋮⋮⋱⋮ρ2B-1ρB-22+B-12⋯1,in which ρ ranges from 0.9 to 1. We set ρ and B to be 0.95 and 32, respectively, by experience.When reconstructing image, we can construct the relationship model between the block observation vectoryi and a whole image x: (4)y=y1y2⋮yn=ΦBΦB⋱ΦBx1x2⋮xn=Φ·Ex=Θx,in which E is the elementary matrix to rearrange the column vectors block by block to a raster scanning column vector of image, and Θ=Φ·E. The minimum l1-l2 norm reconstruction model can be constructed by using (4) as follows: (5)x^=arg⁡minx⁡y-Θx22+λΨx1,in which ·1 and ·2 denote l1 and l2 norm, respectively, Ψ is the sparse representation basis, and λ is the fixed regularization factor. Equation (5) can be solved by using the SPL algorithm with a low computational complexity; in particular Ψ can be updated flexibly at each SPL iteration. Therefore, the PCA can be used to learn the adaptive sparse representation basis.


## 2.2. GrC-Based Clustering

GrC is proposed by Zadeh in 1979, and he said that the information granule exists widely in our daily life and it can be viewed as an abstraction of actual object. Granulating, organization, and causality are three important concepts of GrC [18]; that is, the granulating divides a whole of the object into several parts, the organization combines several parts into a whole by some specific operators, and the causality is used to model the relationship between the cause and its result. In the above, there is a natural relation between GrC and clustering. The training set can be regarded as the object, and each sample in training set is defined as a single granule. From the above views, we can see that the granulating and organization correspond to the clustering and its inverse, respectively, and the causality describes the internal relation between samples. Different from the traditional clustering method (e.g., K-means [19]), the GrC-based clustering chooses a variety of shapes to flexibly represent a granule, and the relationship between granules can be modeled by using the more mature algebraic system; therefore the GrC-based clustering has better generalization ability.Although the group of points in training set has irregular shape, their borders can be better distinguished by GrC-based clustering, since GrC can express the granule as various shapes, for example, hyperdiamond, hypersphere, and hyperbox. The vector form of granule is denoted asG=(c,gr), in which c is the center of granule and gr is the granularity. Denote ci as a point contained in the granule G, and the granularity of G is computed by(6)gr=maxic-cip,in which ·p is p-norm to measure the distance. The different distance measure represents the different shape of granule; for example, p=1 expresses the corresponding granule as hyperdiamond granule, p=2 expresses the corresponding granule as hypersphere, and p=∞ expresses the corresponding granule as hyperbox granule. The operators between granules include the merge operator ∨ and the decomposing operator ∧. The merge operator ∨ is to merge two smaller granules into a big one, and the decomposition operator ∧ is to decompose a big granule into two smaller ones, in which the merge operator ∨ is used to cluster the samples in the training set. The merge granules of two granules G1=(c1,r1) and G2=(c2,r2) are merged as a new granule:(7)G1∨G2=12P+Q,P-Qp,in which(8)P=c1-r1c12p,Q=c2+r2c12p,c12=c2-c1.Figure 1 illustrates the merging of two spherical granules G1=(1.0, 2.0, 2.0) and G2=(2.5, 2.5, 2.5) in the two-dimensional space, and their merging result is G1∨G2=(1.75, 2.25, 1.75).Figure 1
Merging of spherical granules in two-dimensional space.According to the vectorized representation of granule, the relation between granules inevitably involves the relation between vectors, but the partial order relation between vectors is not consistent with inclusion relation between granules. Therefore, we regard that the fuzzy inclusion relation exists between granules. When merging granules, we measure the fuzzy inclusion relation as follows:(9)μG1,G2=vG2vG1∨G2,in which v(G):G→R is a mapping from the space of granule to the space of real numbers, and this mapping v(G) uses the following formula:(10)vx=∑isin⁡πxi2.Aimed to the training setS in the clustering problem, we construct an algebra system G,μ, which is made up of the granule set G, the fuzzy inclusion relation μ(·,·), and merge operator ∨. By the fuzzy inclusion relation between granules, we control the granule merging and divide gradually S into a number of granules in order for each granule to be corresponding to a category. Figure 2 shows the GrC-based clustering results of spherical and square granules in the two-dimensional space; we can see that the GrC performs the granule merging by using the fuzzy inclusion relation and the granularity, which guarantees the better clustering performance; in particular the control of inclusion relation is very beneficial to distinguish the ownership of the training samples in the border. In a natural image signal, some relations exist between image blocks in the spatial domain. Depending on this experience, the GrC-based clustering divides automatically a whole image into several regions including blocks according to the image content. In each cluster, the structure feature of any block is very similar to others; thus it is more conducive to find the sparse representation basis adapted to the fixed data pattern.Figure 2
GrC clustering results of spherical and square granules in two-dimensional space.
(a)
Original training set (b)
Clustering results


## 3. Adaptive Granular Reconstruction

### 3.1. SPL Iteration

In the framework of block CS, the convex optimization algorithm (e.g., NESTA [10]) can be only used to independently reconstruct each block. However, the varying block sparsity results in the uneven qualities of reconstructed blocks due to the different local structure characteristics of image. The uneven reconstructed qualities bring about lots of blocking artifacts. In addition, the small block size restricts the performance of block reconstruction, which makes the reconstructed image contain a lot of noises. Aimed to the above problem, Gan [13] proposed the image reconstruction algorithm based on SPL iteration by combining Projection Onto Convex Set (POCS), hard-thresholding shrinkage, and wiener filtering. The SPL algorithm guarantees the better reconstruction performance with a low computational complexity. The flow of SPL algorithm is presented as shown in the followingThe Flow of SPL IterationTask. Find the optimal solution x^ of model (5).Input: block CS observation vector yi,  i=1,2,…,n; block measurement matrix ΦB; sparse representation basis Ψ; regularization factor λ for hard thresholding.Initialization: set the initial MMSE linear estimator x^i(0)=RxxΦBTΦBRxxΦBT-1,  i=1,2,…,n.Main iteration: increment k by 1, and apply these steps:
(i)
Weiner filtering: the Winner filter with 3 × 3 sliding window is performed to processx^(k)as follows:(11)x^^k=Wienerx^k,33.(ii)
POCS: the filtered image is projected into the convex setC; that is,(12)x-ik=x^^ik+ΦBTyi-ΦBx^^ik,i=1,…,n.(iii)
Hard-thresholding shrinkage: transformx-(k) by using the sparse representation basis Ψ, and perform the hard-thresholding shrinkage as follows:(13)x_k=Ψx-k.σk=medianx_k0.6745x~k=Thresholdx_k,λ=x_k,x_k>λσk2log⁡N0,else,x⌢k=Ψ-1x~k(iv)
POCS: perform the POCS operator again:(14)x^ik+1=x⌢ik+ΦBTyi-ΦBx⌢ik,i=1,…,n.(v)
Stopping rule: stop whenD(k+1)-D(k)<ε, where D(k)=x^(k)-x-(k-1)2/N; ε is the stopping tolerance.Output: the reconstructed image x^.Due to the presence of blocking artifacts and noise in the MMSE linear estimator, we firstly use the 3 × 3 Wiener filter to reduce the blocking artifacts and smoothen the image. Afterward, the filtered image is projected to the convex setC={xi:yi=ΦBxi,i=1,…,n} again, and the transformation coefficients on the sparse representation basis Ψ are shrunk by the hard thresholding to eliminate further noises. Finally, the POCS is performed again to enforce the image into the convex set C. Repeat the above flow until meeting the stopping criterion. In the process of hard-thresholding shrinkage, the solution is not sparse due to the existence of noises. We can get the more sparse transformation coefficients by enforcing the noise components to 0.


### 3.2. GrC-Based PCA Hard-Thresholding Shrinkage

Whether the hard-thresholding shrinkage can eliminate the noises effectively is important for the reconstruction performance of SPL iteration. Gan [13] performs the hard-thresholding operator in the overlapping block DCT and the undecimated wavelet domain, but they cannot better capture the directional characteristics (e.g., point, line, and edge) so as to degrade the reconstructed quality of image. To overcome the above defects, Mun and Fowler [14] proposed the bivariate shrinkage based on directional transform (e.g., contour and double-tree wavelet), which achieves some performance improvements. The above-mentioned methods focus the useful information as much as possible on a few transformation coefficients by constructing some specific wavelets and remove lots of transform coefficients with little information to suppress noises. According to [20], when performing the hard-thresholding shrinkage, it is more effective to capture the important information by learning the dictionary adapted by image content than by constructing the specific wavelet. However, the dictionary learning algorithms (e.g., K-SVD algorithm [21]) have a high computational complexity; therefore the dictionary learning algorithm is not suited to the SPL iteration requiring a low computational complexity. Our previous work [16] proposed the PCA-based hard-thresholding shrinkage to suppress noises. These methods extract the image patches from the noisy image as samples and perform the PCA to train the orthogonal transformation matrix which removes the spatial correlation. Since the PCA matrix can better separate useful components and useless noises existing in image, the hard-thresholding shrinkage gets some performance improvements. In addition, the PCA-based hard-thresholding shrinkage has a low computational complexity; thus it is so suitable to SPL iteration. Due to the fact that the image has only locally stationary statistics and varying structural characteristics globally, the training set composed by global patches will generate a single PCA matrix with the poor ability to extract the sparse coefficients. Depending on the local stationary statistics of image, it is more practical to divide the global patches into several subsets which contain the patches with similar structural characteristics, and PCA is used to generate the sparse representation basis of each subset. These PCA matrices are specialized to some particular structure characteristics; thus they guarantee that all samples in subset have the more sparse representations.According to the above-mentioned viewpoints, we firstly cluster image patches before using PCA to learn the sparse representation basis.K-means [19] is a popular clustering algorithm, but it requires setting the number of clusters in advance, so that the training set cannot be classified adaptively in terms of the data pattern. Due to the nonstationary statistical characteristics of image, a large difference exists between the structure features of patches; thus K-means is not appropriate for the classification of image patches. The GrC does not need to set a fixed number of clusters, and it can achieve the suitable number of clusters by controlling the threshold of granularity to merge the granules depending on the fuzzy inclusion relation between granules. Therefore, we use the GrC to cluster image patches according to the data pattern of training set. The steps of GrC-based PCA hard-thresholding shrinkage are described as follows.Step 1.
We extract allm patches with the size of d=Wp×Wp from the noisy image x-(k) and denote these patches as pi(i=1,…,m). The upper-left pixel of patch pi is corresponding to the ith pixel; thus each pixel is corresponding to a patch. There are overlapping pixels between patches, and those pixels beyond the border of image cannot be extracted; that is, m=(Ic-Wp)×(Ir-Wp).Step 2.
We usepi as the finest atom to constitute the training set S and set the control parameter ρ of granularity. The merge operator ∨ is used to merge the granule conditionally until all training samples are contained in some granule, and finally each granule corresponds to a category. The flow of clustering is shown in the following.The Flow of Dividing the Training Set S Using GrC ClusteringTask. Merge the samples in S into several granules according to the specified conditions.Input: the training set S={p1,p2,…,pm} and the threshold T of granularity.Initialization: use the m samples in S to induce m granule atoms; that is, the sample pl corresponds to the hyperbox granule atom Gl={pl,pl,0}, in which the first entry is the starting point, the second entry is the terminal point, and the third entry is the granularity. In the above, we obtain the initial hyperbox granule set GS={G1,G2,…,GL},L=m.Main iteration: increment i by 1 to m, and apply these steps:
(i)
Compute the fuzzy inclusions betweenGi and other granules of GS; that is, μil=μGi,Gl,l=1,…,L.(ii)
Select the indexQ of granule with the maximum fuzzy inclusion; that is, Q=arg⁡minl⁡μil.(iii)
Compute the granularityρi=Gi-GQ2, and judge whether to merge granules Gi and GQ: if ρi≤T, Gi and GQ are merged into a new granule Gi∨GQ; otherwise Gi is still a member in GS.(iv)
Update the training setS, that is, to remove pi from S.(v)
Stop the above iteration onceS is empty.Output: hyperbox granule set GS.By this flow, we can get the granule setGS={G1,G2,…,GL} composed of L granules.Step 3.
We compute thed×d covariance matrix Σl of the lth granule Gl as follows:(15)Σl=1Gl∑pi∈GlpipiT-p-p-T,p-=1Gl∑pi∈Glpi,in which Gl is the number of samples in Gl.Step 4.
We compute the eigenvaluesλ1≥⋯≥λd of Σl and their corresponding normalized eigenvectors e1,e2,…,ed.Step 5.
The eigenvectors are used to construct the orthogonal transformation matrixΓl=[e1,e2,…,ed], and we transform each patch into Γl domain as follows:(16)p^i=ΓlT·pi.Step 6.
The hard-thresholding shrinkage is performed as follows:(17)p^^i=p^i,p^i≥λσ0,else,σ=medianp^0.6745,p^=p^1T,p^2T,…,p^GlT.Step 7.
The thresholding coefficientsp^^i are transformed inversely into the pixel domain; that is,(18)p~i=Γl·p^^i.Step 8.
We combine all patches into a whole image after shrinking all patches inGl. Because of the existence of overlapped area between patches, these pixels in overlapped area are computed by averaging all pixel values at the same position.By the above flow, each granule corresponds to a PCA matrix, and the additional computational complexityO(m·Wp2) is from the GrC clustering. Therefore, the GrC-based PCA learning algorithm maintains a low computational complexity when compared with the computations of learning global PCA matrix.


## 3.1. SPL Iteration

In the framework of block CS, the convex optimization algorithm (e.g., NESTA [10]) can be only used to independently reconstruct each block. However, the varying block sparsity results in the uneven qualities of reconstructed blocks due to the different local structure characteristics of image. The uneven reconstructed qualities bring about lots of blocking artifacts. In addition, the small block size restricts the performance of block reconstruction, which makes the reconstructed image contain a lot of noises. Aimed to the above problem, Gan [13] proposed the image reconstruction algorithm based on SPL iteration by combining Projection Onto Convex Set (POCS), hard-thresholding shrinkage, and wiener filtering. The SPL algorithm guarantees the better reconstruction performance with a low computational complexity. The flow of SPL algorithm is presented as shown in the followingThe Flow of SPL IterationTask. Find the optimal solution x^ of model (5).Input: block CS observation vector yi,  i=1,2,…,n; block measurement matrix ΦB; sparse representation basis Ψ; regularization factor λ for hard thresholding.Initialization: set the initial MMSE linear estimator x^i(0)=RxxΦBTΦBRxxΦBT-1,  i=1,2,…,n.Main iteration: increment k by 1, and apply these steps:
(i)
Weiner filtering: the Winner filter with 3 × 3 sliding window is performed to processx^(k)as follows:(11)x^^k=Wienerx^k,33.(ii)
POCS: the filtered image is projected into the convex setC; that is,(12)x-ik=x^^ik+ΦBTyi-ΦBx^^ik,i=1,…,n.(iii)
Hard-thresholding shrinkage: transformx-(k) by using the sparse representation basis Ψ, and perform the hard-thresholding shrinkage as follows:(13)x_k=Ψx-k.σk=medianx_k0.6745x~k=Thresholdx_k,λ=x_k,x_k>λσk2log⁡N0,else,x⌢k=Ψ-1x~k(iv)
POCS: perform the POCS operator again:(14)x^ik+1=x⌢ik+ΦBTyi-ΦBx⌢ik,i=1,…,n.(v)
Stopping rule: stop whenD(k+1)-D(k)<ε, where D(k)=x^(k)-x-(k-1)2/N; ε is the stopping tolerance.Output: the reconstructed image x^.Due to the presence of blocking artifacts and noise in the MMSE linear estimator, we firstly use the 3 × 3 Wiener filter to reduce the blocking artifacts and smoothen the image. Afterward, the filtered image is projected to the convex setC={xi:yi=ΦBxi,i=1,…,n} again, and the transformation coefficients on the sparse representation basis Ψ are shrunk by the hard thresholding to eliminate further noises. Finally, the POCS is performed again to enforce the image into the convex set C. Repeat the above flow until meeting the stopping criterion. In the process of hard-thresholding shrinkage, the solution is not sparse due to the existence of noises. We can get the more sparse transformation coefficients by enforcing the noise components to 0.


## 3.2. GrC-Based PCA Hard-Thresholding Shrinkage

Whether the hard-thresholding shrinkage can eliminate the noises effectively is important for the reconstruction performance of SPL iteration. Gan [13] performs the hard-thresholding operator in the overlapping block DCT and the undecimated wavelet domain, but they cannot better capture the directional characteristics (e.g., point, line, and edge) so as to degrade the reconstructed quality of image. To overcome the above defects, Mun and Fowler [14] proposed the bivariate shrinkage based on directional transform (e.g., contour and double-tree wavelet), which achieves some performance improvements. The above-mentioned methods focus the useful information as much as possible on a few transformation coefficients by constructing some specific wavelets and remove lots of transform coefficients with little information to suppress noises. According to [20], when performing the hard-thresholding shrinkage, it is more effective to capture the important information by learning the dictionary adapted by image content than by constructing the specific wavelet. However, the dictionary learning algorithms (e.g., K-SVD algorithm [21]) have a high computational complexity; therefore the dictionary learning algorithm is not suited to the SPL iteration requiring a low computational complexity. Our previous work [16] proposed the PCA-based hard-thresholding shrinkage to suppress noises. These methods extract the image patches from the noisy image as samples and perform the PCA to train the orthogonal transformation matrix which removes the spatial correlation. Since the PCA matrix can better separate useful components and useless noises existing in image, the hard-thresholding shrinkage gets some performance improvements. In addition, the PCA-based hard-thresholding shrinkage has a low computational complexity; thus it is so suitable to SPL iteration. Due to the fact that the image has only locally stationary statistics and varying structural characteristics globally, the training set composed by global patches will generate a single PCA matrix with the poor ability to extract the sparse coefficients. Depending on the local stationary statistics of image, it is more practical to divide the global patches into several subsets which contain the patches with similar structural characteristics, and PCA is used to generate the sparse representation basis of each subset. These PCA matrices are specialized to some particular structure characteristics; thus they guarantee that all samples in subset have the more sparse representations.According to the above-mentioned viewpoints, we firstly cluster image patches before using PCA to learn the sparse representation basis.K-means [19] is a popular clustering algorithm, but it requires setting the number of clusters in advance, so that the training set cannot be classified adaptively in terms of the data pattern. Due to the nonstationary statistical characteristics of image, a large difference exists between the structure features of patches; thus K-means is not appropriate for the classification of image patches. The GrC does not need to set a fixed number of clusters, and it can achieve the suitable number of clusters by controlling the threshold of granularity to merge the granules depending on the fuzzy inclusion relation between granules. Therefore, we use the GrC to cluster image patches according to the data pattern of training set. The steps of GrC-based PCA hard-thresholding shrinkage are described as follows.Step 1.
We extract allm patches with the size of d=Wp×Wp from the noisy image x-(k) and denote these patches as pi(i=1,…,m). The upper-left pixel of patch pi is corresponding to the ith pixel; thus each pixel is corresponding to a patch. There are overlapping pixels between patches, and those pixels beyond the border of image cannot be extracted; that is, m=(Ic-Wp)×(Ir-Wp).Step 2.
We usepi as the finest atom to constitute the training set S and set the control parameter ρ of granularity. The merge operator ∨ is used to merge the granule conditionally until all training samples are contained in some granule, and finally each granule corresponds to a category. The flow of clustering is shown in the following.The Flow of Dividing the Training Set S Using GrC ClusteringTask. Merge the samples in S into several granules according to the specified conditions.Input: the training set S={p1,p2,…,pm} and the threshold T of granularity.Initialization: use the m samples in S to induce m granule atoms; that is, the sample pl corresponds to the hyperbox granule atom Gl={pl,pl,0}, in which the first entry is the starting point, the second entry is the terminal point, and the third entry is the granularity. In the above, we obtain the initial hyperbox granule set GS={G1,G2,…,GL},L=m.Main iteration: increment i by 1 to m, and apply these steps:
(i)
Compute the fuzzy inclusions betweenGi and other granules of GS; that is, μil=μGi,Gl,l=1,…,L.(ii)
Select the indexQ of granule with the maximum fuzzy inclusion; that is, Q=arg⁡minl⁡μil.(iii)
Compute the granularityρi=Gi-GQ2, and judge whether to merge granules Gi and GQ: if ρi≤T, Gi and GQ are merged into a new granule Gi∨GQ; otherwise Gi is still a member in GS.(iv)
Update the training setS, that is, to remove pi from S.(v)
Stop the above iteration onceS is empty.Output: hyperbox granule set GS.By this flow, we can get the granule setGS={G1,G2,…,GL} composed of L granules.Step 3.
We compute thed×d covariance matrix Σl of the lth granule Gl as follows:(15)Σl=1Gl∑pi∈GlpipiT-p-p-T,p-=1Gl∑pi∈Glpi,in which Gl is the number of samples in Gl.Step 4.
We compute the eigenvaluesλ1≥⋯≥λd of Σl and their corresponding normalized eigenvectors e1,e2,…,ed.Step 5.
The eigenvectors are used to construct the orthogonal transformation matrixΓl=[e1,e2,…,ed], and we transform each patch into Γl domain as follows:(16)p^i=ΓlT·pi.Step 6.
The hard-thresholding shrinkage is performed as follows:(17)p^^i=p^i,p^i≥λσ0,else,σ=medianp^0.6745,p^=p^1T,p^2T,…,p^GlT.Step 7.
The thresholding coefficientsp^^i are transformed inversely into the pixel domain; that is,(18)p~i=Γl·p^^i.Step 8.
We combine all patches into a whole image after shrinking all patches inGl. Because of the existence of overlapped area between patches, these pixels in overlapped area are computed by averaging all pixel values at the same position.By the above flow, each granule corresponds to a PCA matrix, and the additional computational complexityO(m·Wp2) is from the GrC clustering. Therefore, the GrC-based PCA learning algorithm maintains a low computational complexity when compared with the computations of learning global PCA matrix.


## 4. Experimental Results and Analysis

The performances of proposed algorithm are evaluated on the five 512 × 512 test imagesLenna,Barbara,Parrot,Leaves, andHouse including the varying smooth, edge, and texture components. Its parameter settings are described as follows: the regularization factor λ is 2.5, the stopping tolerance ε of SPL iteration is 0.01, the patch size Wp is 7, and the threshold T of granularity is 0.03. On one hand, we evaluate the performance of hard-thresholding shrinkage based on the different sparse representation bases. The proposed algorithm uses the GrC-based PCA basis, which is abbreviated to GrC-PCA, and we select the comparing bases including DCT basis, Double-Tree Discrete Wavelet Transform (DDWT) basis, and Global PCA basis (GPCA). On the other hand, we compare the reconstruction performance of proposed algorithm with the MH based SPL (MH_SPL) algorithm [9], the NESTA algorithm with the minimum TV (NESTA_TV) [10], and the Group-based Sparse Recovery (GSR) algorithm [11]. In all experiments, the block size B is 32; the preset total measurement rate SR is set from 0.2 to 0.6. The Peak Signal-to-Noise Ratio (PSNR) is used to evaluate the objective performance, and all PSNR values are averaged over 5 trials, since the quality of reconstruction can vary due to the randomness of random transformation matrix ΦB. All experiments run under the following computer configuration: Intel(R) Core(TM) i7 @ 3.40 GHz CPU, 8 GB RAM, Microsoft Windows 7 32 bits, and MATLAB Version 7.6.0.324 (R2008a).


### 4.1. Performance Comparison of Hard-Thresholding Shrinkage

Figure3 presents the average RD performance on the five test images of the SPL algorithm with the hard-thresholding shrinkages based on DCT, DDWT, GPCA, and GrC-PCA. We can see that our GrC-PCA basis obtains the highest PSNR value among all bases at any measurement rate, and particularly it has about 1 dB gains compared to the DCT basis. Our GrC-PCA basis exploits the local stationary structural characteristics of image; thus the PSNR of GrC-PCA basis is 0.2 dB larger than GPCA basis. Figure 4 shows the performance of hard-threshold shrinkage using different clustering algorithms. We can see that the RD performance degrades once using the traditional K-means algorithm. However, the PSNR value obtains the obvious improvement at any measurement rate when using the GrC-based clustering method; particularly when the measurement rate is 0.2, the PSNR gain of GrC-PCA is around 0.6 dB.Figure 3
The comparison of hard-thresholding shrinkage performances under the different sparse representation bases.Figure 4
The comparison of hard-thresholding shrinkage performances under the different clustering methods.


### 4.2. Comparison of Reconstruction Performance

To evaluate the RD performance of our GrC-PCA-based SPL algorithm, we select the popular reconstruction algorithms as benchmarks, and they include MH_SPL [9], NESTA_TV [10], and GSR [11]. Table 1 lists the average PSNR value on the five test images when the measurement rate ranges from 0.2 to 0.4. We can see that our method obviously outperforms MH_SPL and the NESTA_TV, and it obtains the PSNR gain of 0.48 dB and 3.27 dB, respectively. However, when compared with the GSR algorithm, the proposed algorithm in this paper deteriorates 0.64 dB on average, which results from the fact that the GSR uses both the local stationary statistics of image and nonlocal similarities of patches. Therefore, the GSR obtains more sparse representation compared to GrC-PCA. However, it can be observed from Table 2 that the GSR algorithm has a high computational complexity, and it requires 2542.65 s to reconstructLenna on average, but our algorithm requires only 30.20 s. In the above, when considering both the reconstruction quality and computational complexity, the proposed algorithm is more practical than GSR. Compared with the DDWT-based MH_SPL algorithm, our algorithm increases 18.54 s on average, which results from the fact that PCA needs to relearn the sparse representation basis at each iteration, so that the proposed algorithm obtains the PSNR improvement at the cost of computational complexity. Figure 5 shows the visual qualities of reconstructedLenna using the different algorithms when the measurement rate is 0.3; we can see that the reconstructed image by NESTA_TV contains lots of blocking artifacts, and MH_SPL algorithm eliminates the blocking artifacts but introduces the noises. However, our algorithm gets the subjective visual quality similar to the GSR algorithm. Figure 6 shows the visual qualities of reconstructed remote sensing imageKirtland using the different algorithms when the measurement rate is 0.3. It can be observed that our method provides better visual quality for the remote sensing application when compared the other reconstruction algorithms.Table 1
Comparison of average PSNR values of the 5 test images using the different methods.
Reconstruction algorithm Measurement rateSR 0.2 0.3 0.4 Average MH_SPL 31.62 33.52 34.82 33.32 NESTA_TV 28.56 30.89 32.15 30.53 GSR 32.56 34.82 35.95 34.44 Proposed 32.12 34.08 35.21 33.80Table 2
Comparison of consumed times when using the different methods to reconstructLenna image.
Reconstruction algorithm Measurement rateSR 0.2 0.3 0.4 Average MH_SPL 6.59 12.52 15.88 11.66 NESTA_TV 44.98 44.63 44.81 44.81 GSR 2272.15 2519.23 2836.56 2542.65 Proposed 32.05 30.61 27.95 30.20Figure 5
Comparison of visual qualities using the different methods to reconstructLenna image when the measurement rate is 0.3.
(a)
NESTA_TV (b)
MH_SPL (c)
GSR (d)
ProposedFigure 6
Comparison of visual qualities using the different methods to reconstructKirtland image when the measurement rate is 0.3.
(a)
NESTA_TV (b)
MH_SPL (c)
GSR (d)
Proposed


## 4.1. Performance Comparison of Hard-Thresholding Shrinkage

Figure3 presents the average RD performance on the five test images of the SPL algorithm with the hard-thresholding shrinkages based on DCT, DDWT, GPCA, and GrC-PCA. We can see that our GrC-PCA basis obtains the highest PSNR value among all bases at any measurement rate, and particularly it has about 1 dB gains compared to the DCT basis. Our GrC-PCA basis exploits the local stationary structural characteristics of image; thus the PSNR of GrC-PCA basis is 0.2 dB larger than GPCA basis. Figure 4 shows the performance of hard-threshold shrinkage using different clustering algorithms. We can see that the RD performance degrades once using the traditional K-means algorithm. However, the PSNR value obtains the obvious improvement at any measurement rate when using the GrC-based clustering method; particularly when the measurement rate is 0.2, the PSNR gain of GrC-PCA is around 0.6 dB.Figure 3
The comparison of hard-thresholding shrinkage performances under the different sparse representation bases.Figure 4
The comparison of hard-thresholding shrinkage performances under the different clustering methods.


## 4.2. Comparison of Reconstruction Performance

To evaluate the RD performance of our GrC-PCA-based SPL algorithm, we select the popular reconstruction algorithms as benchmarks, and they include MH_SPL [9], NESTA_TV [10], and GSR [11]. Table 1 lists the average PSNR value on the five test images when the measurement rate ranges from 0.2 to 0.4. We can see that our method obviously outperforms MH_SPL and the NESTA_TV, and it obtains the PSNR gain of 0.48 dB and 3.27 dB, respectively. However, when compared with the GSR algorithm, the proposed algorithm in this paper deteriorates 0.64 dB on average, which results from the fact that the GSR uses both the local stationary statistics of image and nonlocal similarities of patches. Therefore, the GSR obtains more sparse representation compared to GrC-PCA. However, it can be observed from Table 2 that the GSR algorithm has a high computational complexity, and it requires 2542.65 s to reconstructLenna on average, but our algorithm requires only 30.20 s. In the above, when considering both the reconstruction quality and computational complexity, the proposed algorithm is more practical than GSR. Compared with the DDWT-based MH_SPL algorithm, our algorithm increases 18.54 s on average, which results from the fact that PCA needs to relearn the sparse representation basis at each iteration, so that the proposed algorithm obtains the PSNR improvement at the cost of computational complexity. Figure 5 shows the visual qualities of reconstructedLenna using the different algorithms when the measurement rate is 0.3; we can see that the reconstructed image by NESTA_TV contains lots of blocking artifacts, and MH_SPL algorithm eliminates the blocking artifacts but introduces the noises. However, our algorithm gets the subjective visual quality similar to the GSR algorithm. Figure 6 shows the visual qualities of reconstructed remote sensing imageKirtland using the different algorithms when the measurement rate is 0.3. It can be observed that our method provides better visual quality for the remote sensing application when compared the other reconstruction algorithms.Table 1
Comparison of average PSNR values of the 5 test images using the different methods.
Reconstruction algorithm Measurement rateSR 0.2 0.3 0.4 Average MH_SPL 31.62 33.52 34.82 33.32 NESTA_TV 28.56 30.89 32.15 30.53 GSR 32.56 34.82 35.95 34.44 Proposed 32.12 34.08 35.21 33.80Table 2
Comparison of consumed times when using the different methods to reconstructLenna image.
Reconstruction algorithm Measurement rateSR 0.2 0.3 0.4 Average MH_SPL 6.59 12.52 15.88 11.66 NESTA_TV 44.98 44.63 44.81 44.81 GSR 2272.15 2519.23 2836.56 2542.65 Proposed 32.05 30.61 27.95 30.20Figure 5
Comparison of visual qualities using the different methods to reconstructLenna image when the measurement rate is 0.3.
(a)
NESTA_TV (b)
MH_SPL (c)
GSR (d)
ProposedFigure 6
Comparison of visual qualities using the different methods to reconstructKirtland image when the measurement rate is 0.3.
(a)
NESTA_TV (b)
MH_SPL (c)
GSR (d)
Proposed


## 5. Conclusion

In this paper, we use the GrC-based clustering to improve the performance of PCA-based hard-thresholding shrinkage in the SPL iteration. Considering the local stationary statistical property, the image is firstly decomposed into some granules. Since GrC can divide adaptively the training set according to the data pattern, the samples in granule have the similar structure. Then, the PCA is performed to learn the sparse representation basis dedicated to the granule, and the patches in granule are shrunk by the hard thresholding using PCA matrix. Finally, all patches merge into a whole image. The experimental results show that the proposed algorithm has better RD performance and guarantees the better subjective visual quality of reconstructed image with a low computational complexity.



---

*Source: 1280690-2016-11-28.xml*