Designing sparse sensing matrix for compressive sensing to reconstruct high resolution medical images

: Compressive sensing theory enables faithful reconstruction of signals, sparse in domain Ψ , at sampling rate lesser than Nyquist criterion, while using sampling or sensing matrix Φ which satisfies restricted isometric property. The role played by sensing matrix Φ and sparsity matrix Ψ is vital in faithful reconstruction. If the sensing matrix is dense then it takes large storage space and leads to high computational cost. In this paper, effort is made to design sparse sensing matrix with least incurred computational cost while maintaining quality of reconstructed image. The design approach followed is based on sparse block circulant matrix (SBCM) with few modifications. The other used sparse sensing matrix consists of 15 ones in each column. The medical images used are acquired from US, MRI and CT modalities. The image quality measurement parameters are used to compare the performance of reconstructed medical images using various sensing matrices. It is observed that, since Gram matrix of dictionary matrix ( ΦΨ ) is closed to identity matrix in case of proposed modified SBCM, therefore, it helps to reconstruct the medical images of very good quality.

ABOUT THE AUTHOR Vibha Tiwari received BE and ME degree in 1999 and 2002, respectively. She is currently pursuing PhD. Her current research interests include medical image compression, super-resolution image processing, colourization, image restoration and enhancement. The work presented in the paper aims to design sensing matrix so as to transfer low size medical image requiring less storage space and low transmission bandwidth and at remote end reconstruct high fidelity image

PUBLIC INTEREST STATEMENT
The article titled 'DESIGNING SPARSE SENSING MATRIX FOR COMPRESSIVE SENSING TO RECONSTRUCT HIGH RESOLUTION MEDICAL IMAGES'. This article aims to give contribution to global Health Care system.
In case of intricate disease, patient located at another place from where his ultrasound, computed tomography or any other related modality images may be taken and transmitted remotely where more specialized doctors are located for better consultation at affordable cost.
Generally, medical images are of huge size, occupies huge space and bandwidth to transmit remotely. So, if low resolution images are captured then radiation dosage given to patient will reduce and it will take lesser bandwidth for transmission. This low resolution image is then used to reconstruct high resolution image using sparse recovery algorithms in which sensing matrix design plays a vital role. The sensing matrix is designed in this paper for reconstructing high resolution images.

Introduction
Technological growth has reduced the gap between domestic and global Health Care system. Telemedicine is one such field which promises to provide affordable and quality health care facilities to people at large. The successful deployment of Telemedicine depends heavily on storage and transfer of huge size clinical image/video. Various compression techniques have been developed which transforms the signal in some suitable domain and then achieve compression by discarding less significant coefficients. For reconstruction, signal must be sampled at the Nyquist rate. However, compressive sensing (CS) technique can recover signals from lesser number of samples or measurements of sparse signal (Donoho, 2006). Generally, real time signals are sparse in transform domain Ψsuch as wavelet or discrete cosine transform. Work done by Candes, Tao and Romberg (2006) proves that if the sparse sensing matrix Φ consists of small number (k) of non-zero elements drawn from iid Gaussian or Bernoulli random variables then it is possible to recover the k-sparse signal of length N, from lesser number of M measurements, using linear program, on condition that k ≲ M log( N M ) (Candes & Romberg, 2007;Candes, Romberg, & Tao, 2006a). However, computational cost and storage space of using such dense sensing matrix is very high. Therefore, lot of work is going on to design efficient sensing matrix.
In Zelnik-Manor and Eldar (2011), sensing matrix is optimized by minimizing a weighted sum of the interblock coherence and subblock coherence of the equivalent dictionary. The nearly orthonormal blocks of equivalent dictionary is obtained by minimizing mostly the subblock coherence resulting in slight increase in total interblock coherence. Equiangular tight frame (ETF) approach to Gram matrix of equivalent dictionary matrix is used in Shuang, Zhihui, Li, Chang, and Li (2013) to minimize weighted sum of both interblock and subblock coherence. In Shuang et al. (2013), best recovery results were obtained when total subblock coherence is minimized in contrast to Zelnik-Manor and Eldar (2011).
Statistical restricted isometry property (STRIP) is used in Gan, Cong, Thong, and Trac (2009) to prove that deterministic matrices reconstructs signal better comparative to random matrices. The deterministic sensing matrices used were partial FFT matrix, partial Walsh-Hadamard matrix and matrix constructed from kasami code. Deterministic sensing matrix was also used in Lorne, Stephen, Stephen, and Calderbank (2009), which was based on chirp signal. The validity of sensing matrix was proved empirically and by bounding the eigenvalues of submatrices. In Iwen (2014), compressed sensing and recovery algorithm were proposed which performs recovery in O((k log k) logN)-time and provides error bounds as well. The sensing matrix is constructed by randomly subsampling rows from incoherent matrix.
The goal of this paper is to design sensing matrix which will help in reconstructing the image faithfully based on SBCM (Sun, Shu, & Dong, 2013). In this paper, low resolution medical images: Ultrasound (US), Computed Tomography (CT) and Magnetic Resonance Imaging (MRI) of size 128 × 256 have been considered for simulation and analysis. Block Sparse Bayesian Learning Bounded Optimization (BSBL_BO) (Zhang & Rao, 2013) and Smooth approximation to L 0 norm (SL0) (Mohimani, Babaie-Zadeh, & Jutten, 2007) CS algorithms are used for high resolution image reconstruction.
The rest of the paper is organized as follows. In Section 2, design of sparse sensing matrix is discussed. Simulation and results are discussed in Section 3.

Sparse sensing matrix for CS
Real time signals/images are sparse in transform domain. Considering k-sparse signal s ∈ R N sparse in Ψ basis matrix, then a sensing matrix Φ ∈ R M×N can recover the signal with high probability if it satisfies restricted isometric property (RIP) given as Xu and Xu (2013), where k is restricted isometric constant (RIC) which should be very much smaller than 1.

Defining,
Equivalently (1) can be stated as Candes, Romberg, and Tao (2006b), for every subset T ⊂ {1, … , p}, having no more than k indices, A T a submatrix of size M × |T|obtained from matrix A by considering number of columns corresponding to indices in T. However, evaluating RIP is impractical. The alternative way to find sensing matrix which satisfy RIP and hence guarantee to reconstruct signal is to obtain coherence value between sensing and sparsifying matrix. The two matrices should be incoherent for perfect reconstruction, so the value of coherence should be very small (Ben-Haim, Eldar, & Elad, 2010). The coherence can be obtained by taking inner product between two distinct normalized columns of matrix A as Duarte et al. (2006), Since A is obtained after normalizing columns, so its Gram matrix, given as G = A T A , should have diagonal elements close to one. Further analysis of sensing matrix can be done by obtaining average value of the off-diagonal element of G, which is given as Abolghasemi, Ferdowsi, Makkiabadi, and Sanei (2010), Ideally, G should be close to identity matrix of size N × N and coherence and average absolute off-diagonal value must be zero (Abolghasemi et al., 2010). According to Gershgorin's disc theorem, Gram matrix of sub-matrix A T given as, G = A T T A T , should have eigenvalues bounded in magnitude 1 − k , 1 + k for A to satisfy RIP(k, k ). However, RIP is a satisfactory condition but not necessary condition for perfect reconstruction. Mamaghanian, Khaled, Atienza, and Vandergheynst (2011) proposed three alternatives to design k−sparse sensing matrix satisfying RIP1 property. The non-zero elements consisted of entries equal to ± √ k , √ k or k number of ones in each row. The location of non-zero entries was randomly chosen.
Sensing matrix having k number of ones had lower coherence value with Daubechies db10 wavelet sparsity matrix; therefore, it is used in this paper.
In order to achieve lower computational cost, SBCM for CS, was proposed by Sun et al. (2013). The sensing matrix Φ ∈ R M×N considered was a sparse block circulant matrix. The matrix Φ, was divided into sub-matrices. There were L number of columns (L = N∕d) and D number of . The circulant nature of matrix is given as, Φ ij is a sparse matrix having elements denoted by b q (q = 1, … , d). Each row and column of Φ ij had one non-zero element which belong to L set of iid Gaussian random variables. Each row and column of Φ had different non-zero elements. Thus in a row there were L number of coefficients and in a column D number of coefficients. It had been proved that the diagonal elements of Gram matrix of Φ(G = Φ � Φ) is close to 1 and off-diagonal elements are bounded in magnitude with high probability and hence satisfies RIP.
The SBCM approach for designing sensing matrix is further modified in the paper to reduce computational load without degrading reconstructed image quality. Instead of using iid Gaussian random variables, non-zero elements are generated by using inverse Daubechies db2 wavelet coefficients. Since random matrices are completely unstructured, so they increase computation complexity as well as memory requirement for storage, hence costly to implement in practical applications (Do, Gan, & Tran, 2012). Therefore, deterministic approach is followed in choosing the elements of sensing matrix, using this approach further reduces the off-diagonal elements of matrix G.
The structure of sensing matrix is maintained as earlier and the location of non-zero elements in sensing matrix is kept random. The flowchart of proposed design of sensing matrix is shown in Figure   1. Image matrix of size N × N is transformed using Daubechies db10 wavelet basis matrix, from which the sparse column vector s ∈ R N is extracted. The measurement vector y ∈ R M is obtained by multiplying Φ with sparse vector s as, The reduced M number of samples of s are scaled, to bound the magnitude, before giving it as an input to the CS algorithm for reconstruction. The CS algorithm used is SL0. Here, N is set to 256 and M to 128. The sensing matrix Φ of size 128 × 256, is divided into sub-matrices each of size d × d, To analyse the sensing matrix further, Gram matrix G is obtained for three different sensing matrix containing non-zero elements as: (i) fifteen ones in each column present at random locations; (ii) sparse block circulant matrix with iid Gaussian random values, with block-size 32 × 32; and (iii) modified sparse block circulant matrix with sub-matrix size of 128 × 128. The three cases selected are the ones which gave the best results in the respective sensing matrix design method discussed above.
Gram matrix G, of size 256 × 256, is given as, G =Â TÂ , where Â is a normalized column version of matrix A and Â T is the transpose. The mutual coherence, average value of absolute off-diagonal Binary sensing matrix has least coherence value and modified SBCM with 128 × 128 block size has least av and its eigenvalue is bounded to 2.08. The histogram of off-diagonal elements of G is shown in Figure 2(b), (d) and (f). It shows that in case of (a) sensing matrix having 15 ones, 60% of elements are below 0.1, (b) SBCM matrix has 83% elements below 0.1 and (c) modified SBCM has 97% of elements lying below 0.1. The distribution of matrix G of size 256 × 256, is shown in Figure 2(a), (c) and (e) which indicates that for modified SBCM, G has non-zero elements close to unity present in diagonal and just off-diagonal position and the rest elements are close to zero. So G, in case of modified SBCM, closely resembles identity matrix and this feature is useful for faithful reconstruction of images. Off−diagonal Elements value

Simulation and results
Medical images of US, MRI and CT available online in DICOM format are used for simulation (Barre, 2014;DICOM, 2014;DICOM Test Images, 2014;Sample DataSet, 2014). Sparsifying matrix is generated from Daubechies db10 wavelet coefficients. Simulation is performed on 50 different images acquired from each type of modality using three different sensing matrices as defined in Section 2. SL0 and BSBL_BO CS algorithms are used for reconstruction of medical images. In BSBL_BO, input parameters are set to the default values as threshold set for pruning is 1e-4, noise variance is given by std(y) 2 ∕1e2, block size is set to 16 and stopping criterion is set to 1e-8. The input parameter of SL0 algorithm is set to default values as min = 0.001 and decrease_factor = 0.5.
The quality of reconstructed images is evaluated using parameters like Mean Square Error (MSE), Peak Signal to Noise Ratio (PSNR), Structural Similarity Index (SSIM), Normalized Cross-Correlation (NKR) and Normalized Absolute Error (NAE). The formula used to calculate performance metrices are as given below (Grgic, Grgic, & Mrak, 2004). MSE is obtained by using relation, where x i, j and y(i, j) are original and reconstructed images, respectively, of size N × N and i and j are the pixel indices. PSNR is obtained as, where n is the number of bits used to represent a pixel. SSIM can be given as, where mean values of the two images are represented by x , y , standard deviation as x , y respectively, C 1 and C 2 are the two constants used to avoid instability. NKR is obtained as, and NAE is obtained by using formula, The vital parameters amongst these are PSNR, MSE and SSIM. PSNR gives the ratio between maximum possible power of original image to the power of error introduced while reconstruction. MSE measures pixel wise distortion and SSIM is used to compute similarity between the original and reconstructed images based on luminance, contrast and structure. Since in medical images, structure of images should not be distorted, therefore SSIM plays a very important role in deciding the quality of reconstructed images. Computational load is calculated based upon CPU processing time. The simulation tool used is Matlab running on a PC with Intel core i3 processor and 4 GB RAM.
Column vectors of sparse medical image in wavelet domain are multiplied by sensing matrix of size 128 × 256. Thus, lesser number of measurements of image of size 128 × 256 are obtained. This low resolution image is given as input to CS algorithm for reconstruction of 256 × 256 size image.
The mean values of different quality indicator parameters are obtained after running simulation over 50 images of each modality. The values obtained for SBCM sensing matrix is obtained over 50 trials. The sensing matrix with 15 ones is denoted by Bin_Phi and the CS algorithm used is mentioned within the brackets in the table. In case of SBCM and modified SBCM sensing matrix SL0 CS algorithm is used.
US images of resolution 128 × 256 are given as input to CS algorithms SL0 and BSBL_BO which reconstructs 256 × 256 size images. The mean values of quality indicator parameters for US image are given in Table 2. PSNR values indicate that the sensing matrix consisting of 15 ones with BSBL_BO CS algorithm performs better. However, mean values of other parameters like SSIM, MSE and NKR values are far better in case of modified SBCM. Lower value of MSE and NAE with very good SSIM value of about 0.913 and NKR of about 0.9737 indicates that good quality US images are reconstructed using modified SBCM, as also evident from Figure 3.
CPU processing time with SL0 algorithm is least as compared to BSBL_BO. Total number of nonzero elements in modified SBCM is 128, in SBCM 1024 and 3840 elements in binary sparse sensing matrix, therefore, this is one of the cause of increased CPU processing time in case of binary sparse sensing matrix and reduced timing in modified SBCM, as also shown in Figure 6(a). In typically 1.21 s US image is reconstructed using modified SBCM sensing matrix.   Table 3. Good quality MRI images are reconstructed as shown in Figure 4, with highest PSNR of about 37.44 , 0.9948 NKR, very good SSIM value of about 0.8447 and lowest MSE and NAE values in typically 1.23 s.
CT images of resolution 128 × 256 are very difficult to reconstruct. As shown in Figure 5, both BSBL_BO and SL0 with sparse binary sensing matrix are unable to reconstruct the low resolution images. The quality indicator parameter values are shown in Overall performance of sensing matrices can be evaluated from graphs of Figure 6. CPU processing time taken is highest for binary sparse sensing matrix and least for modified SBCM as shown in Figure 6   Sampling images in appropriate manner is essential for faithful reconstruction of images. Thus, sampling or sensing matrix plays a significant role in CS. Since Gram matrix G of modified SBCM is closest to identity matrix and then is the G of sparse binary sensing matrix and last is G of SBCM so accordingly sensing matrices reconstructs images having quality in similar order.
Further quality of reconstructed images can be enhanced by using image denoising filters or by developing algorithms based on optimization techniques such as particle swarm optimization method or fuzzy logic method (Shibin, Qingsong, & Yaoqin, 2013). Speckle noise prominent in US images can be removed by performing adaptive homomorphic or non-homomorphic wavelet filtering (Rizi, Ahmadi, & Setarehdan, 2011), Yuvarani and Kumarasamy (2012). MRI images are affected by noise having RICIAN distribution which is difficult to remove, since it is dependent on signal. Adaptive wavelet domain filter may be used for CT and MRI image denoising.
The proposed method can be applied to colour images as well. The luminance (Y) and chrominance (C b and C r ) parts of image may be extracted and then the method may be adopted on the Y, C b and C b -C r transformed parts. The designed sensing matrix may be used to reduce the luminance and chrominance information. At receiver, YC b C r components can be recovered back by applying CS algorithm and then performing inverse transformation. C r part can be obtained by adding the recovered chrominance information. The YC b C r components are then converted back to RGB domain to recover transmitted colour image. In future, the applicability of this technique may be tested on colour images by evaluating reconstructed image quality measurement parameters.    Thus, low resolution images can be obtained in clinical settings with lesser storage and in lesser time. This enables the medical image/signals to be sent to distant diagnostic centre over low data rate transmission channels. At remote end, high resolution image/signals can be recovered for better clinical assessment by specialized doctors. The proposed method is applied on grey images; the method may be tested with colour images.