Image reconstruction through compressive sampling matching pursuit and curvelet transform

ABSTRACT


INTRODUCTION
Image reconstruction improves the quality of images for better visual interpretation, understanding, and analysis.The applications of image reconstruction are robotics, entertainment, augmented reality, human-computer interaction, satellite image processing, astrological studies, and surveillance.Image reconstruction involves the entire image formation process.The aim of image reconstruction is to recover image information that has been lost during image formation.It is a method to reconstruct a proper image from degraded image by using various models.Image reconstruction is broadly categorized into three types based on input data.The three types of reconstruction involve image denoising, image restoration and image super-resolution.Image denoising can be defined as the problem of mapping from a noisy image to a noise-free image.Identifying the noise parameters is the crucial step of denoising.Next challenging task is to remove noise while preserving edge information.Image restoration is described as the reconstruction of the uncorrupted image from blurred and noisy image.Essentially, it attempts to accomplish the inverse of the imperfections in the image formation system.Recovering high-resolution image from one or more low-resolution images is called super-resolution.
In iterative reconstruction, an image is reconstructed by iteratively optimizing an objective function, which normally includes a data fidelity term and an edge-preserving regularization term [1].Greedy algorithms adapt the iterative reconstruction.Matching pursuit (MP) and orthogonal matching pursuit (OMP)  ISSN: 2088-8708 Int J Elec & Comp Eng, Vol. 13, No. 6, December 2023: 6277-6284 6278 are the initial algorithms that evolved from greedy techniques.A dictionary is used to represent the input signal using as few components as possible.Each column of the dictionary is known as an atom.For each iteration, the best optimal single atom is selected from the dictionary and the value of the corresponding coefficient is updated.Here, only one atom is added to the support set in each iteration [2].Regularized orthogonal matching pursuit (ROMP) has less run time but the quality of reconstruction is low [3].
The present work addresses these problems through the CoSaMP algorithm.The proposed method identifies 2s (s is the sparsity level of the signal) atoms and adds it to the support set.The restricted isometry property (RIP) is developed for analyzing l1 minimization and it plays a key role in CoSaMP.If the matrix Φ satisfies the RIP of order k, this means that every subset of k columns of the matrix is approximately orthonormal.This property is the main criteria for strong convergence results of CoSaMP technique.The reconstruction guarantee of CoSaMp is uniform and stable.Also, an important feature of CoSaMP is that new indices in a signal estimate can be added or deleted from the current set of chosen indices.In contrast to this, MP, OMP and ROMP suffer from the fact that a chosen index remains in the signal representation until the end.
The organization of the remainder of this paper is as follows.The related work regarding image reconstruction is discussed in section 2. Section 3 presents the mathematical foundation of CoSaMP.Section 4 discusses the reconstruction by the CoSaMP algorithm.Section 5 describes experimental evaluation and summarizes the results and section 6 concludes the paper.

RELATED WORK
Tezcan et al. [4] employed variational autoencoders (VAE) as an explicit prior term in reconstruction.This resulted in a powerful image prior which compensated for missing k-space data which also does not require paired datasets for training.Lal et al. [5] proposed structured illumination microscopy (SIM) algorithm for super-resolution microscopy.This enabled researchers to interactively vary various parameters and analyze their effect on the reconstructed image.
Schlemper et al. [6] developed a technique for dynamic magnetic resonance image reconstruction.A deep cascade of convolutional neural networks (CNNs) is used to speed up the data acquisition process.The proposed method outperformed the dictionary learning-based MR image reconstruction.One of the limitations of deep neural network (DNN) is requirement of huge amount of prior training pairs.On the lines of deep image prior framework, Gong et al. [7] presented a personalized network training technique in which prior training pairs are not needed and only the patient's own prior information was required.The optimization problem of the maximum-likelihood estimation was solved through alternating direction method of multipliers algorithm.Tan et al. [8] opined that due to ill-posed problem of image reconstruction, conventional methods cannot achieve high accuracy.Hence, CNN was proposed for mapping complicated nonlinear function for image reconstruction of electrical resistance tomography (ERT).Zhu et al. [9] discussed a novel technique, automated transform by manifold approximation (AUTOMAP) for image reconstruction which makes use of supervised learning that maps between the sensor and the image domain which is driven by the training data.A DNN with AUTOMAP is implemented for learning reconstruction transforms for different image acquisition strategies.It can be observed that reconstruction artefacts are reduced to a great extent when compared to conventional reconstruction methods and also immunity to noise was higher.
In plug-and-play algorithms, full set of measurements are used at every iteration which is a hindrance for applications which have large datasets.Sun et al. [10] addressed this limitation by introducing online extension called plug-and-play stochastic proximal gradient method (PnP-SPGM).Traditional batch PnP framework was extended with new online algorithm PnP-SPGM which can be used for large-scale image reconstruction [10].Wang et al. [11] discussed a compatibly sampling reconstruction network (CSRNet) in which an initial image with acceptable quality is obtained by deep reconstruction network module.Further, CNN is employed to improve the image by residual reconstruction network module [11].Yao et al. [12] proposed deep residual reconstruction network (DR2) which has two stages.In the first stage, high-quality preliminary image is reconstructed by linear mapping and in the second stage, residual learning is used to further improve the reconstruction quality [12].Liu et al. [13] presented an end-to-end multiscale residual CNN, called as MSRNet.Different convolution kernel sizes are used with three parallel channels to exploit different-scale feature information in the reconstruction stage.Also, residual learning is used to speedup training process [13].
Meenakshi and Budhiraja [14] presented OMP for image reconstruction in which stopping criteria for the algorithm is modified.Stopping criteria is determined by the recovery condition and mutual incoherence property on the low frequency coefficients of the image.Jin and Li [15] discuss an improved OMP called as Fourier generalized orthogonal matching pursuit (FgOMP) algorithm for electron probe Extensive experiments are carried out on various measurement matrices and Fourier measurement matrix has been selected for better results.
Wei et al. [16] discuss an adaptive variational partial differential equation (PDE) model.A novel method is introduced combining L2 energy of the image gradient and the total variation (TV) to preserve the most important image features.Even though conventional neural network models do well in feature representation and classification, their performance is reduced for high level features.Towards this, Sungheetha and Sharma [17] proposed capsule network which is suitable to handle high level features.Capsule consists of a group of neurons in which the instantiation parameters are considered as activity vectors to represent the entity of an object.These networks provide more information in a vector format.This better information leads to proper decision making.One of the drawbacks of DNN-based image reconstruction is it needs a large amount of memory and more time for the training.Towards this, Wu et al. [18] proposed the proximal gradient descent algorithm for finite iterations.The penalty function is replaced with trainable convolutional neural network (CNN).
Conventional image reconstruction for clinical positive emission photography (PET) lacks automation for the optimization of advanced image reconstruction algorithms and the computational expense also is high.For addressing these points, a novel deep architecture named as DeepPET is developed which directly and immediately reconstructs PET sinogram data into high quality images [19].Połap and Srivastava [20] presented a technique in which some important area of image is found by heuristic algorithm.Different types of neural networks are trained until a certain level of entropy of these areas is reached.Average entropy of found areas is used for stopping criteria.In photoacoustic tomography (PAI), hybrid imaging is used where both optical and ultrasonic imaging advantages are combined.Even though conventional algorithms used in PAI provide a speedy solution, many artifacts still exist.To overcome this drawback, Lan et al. [21] proposed a new CNN framework known as Y-Net.In Y-Net, both raw data and beamformed images are optimized to reconstruct the initial PA pressure distribution.In low-dose computed tomography (CT) imaging, both analytic and iterative reconstruction algorithms are prone to different types of image artifacts.These methods exploit the DNN in the image domain of CT. Lee et al. [22] presented CT in the sinogram domain for restoring the missing data by implementing CNN using the Caffe library.Missing data is synthesized in sparse-view sinogram.Generally, signal and system models are used for solving an inverse problem for standard reconstruction methods of photoacoustic imaging.But, in handheld scanners, there will be low image contrast and structure loss due to limited detection view and bandwidth.Kim et al. [23] proposed deep CNN to overcome these problems which is designed for real time clinical applications.Training is carried out by large scale synthetic data imitating typical microvessel networks.
When the contrast or resolution of the acquired images differ considerably from training images used to learn the prior, reconstruction quality decrease [4].Even though OMP guarantees quality of reconstruction, its reconstruction time is more.Also, the reconstruction guarantees are not uniform [24].Uniform guarantees can only be achieved with more number of measurements.ROMP has less run time but at the same time compromising with the quality of reconstruction.Also, it will not work on noisy signal.Generally, the images are not sparse.The choice of transform basis will affect the quality of the reconstructed image.Curvelets result in a representation that is sparser than other wavelet transforms.It also offers exact reconstruction, low computational complexity and stability against perturbation [25].In this paper, we are presenting compressive sampling matching pursuit technique which can deliver higher PSNR and better image quality.

MATHEMATICAL FOUNDATION OF CoSaMP
Let x is a signal which is s-sparse where ‖‖ 0 ≤ .One of the difficult parts of reconstruction of signal is to find the locations of the largest components in the target signal. is used as sampling matrix which is also sometimes called as measurement matrix.It has to satisfy the RIP.Random matrices like Gaussian matrix, Bernoulli matrix, and almost other matrices with independent and identically distributed (i.i.d) entries satisfy RIP.
RIP states that the th restricted isometry constant of a matrix  is the least number   for which.
The signal proxy  for an s-sparse signal x, is a vector which can be formed as (2).
Largest  components of y correspond to largest  components of x.The samples are constructed as (3).
Hence, the proxy can be obtained by multiplying matrix Φ * to the samples as in (4).
The target signal is approximated by successive iteration.At each iteration, the approximation results in a residual.As the iteration is continued, samples are updated based on the current residual.A proxy for the residual is constructed using these samples.By this, the large components in the residual are identified and a tentative support is formed for the next approximation.Using the samples, the approximation on this support is estimated by least squares technique.This progression is continued until the halting criteria is reached [26].

PROPOSED METHOD
The main idea of this paper is to propose CoSaMP algorithm for image reconstruction to yield better image quality with reduced number of samples.The following four inputs are required for the algorithm: i) sampling operator; ii) samples of the unknown signal; iii) sparsity of the output signal; iv) halting criteria.The algorithm is initialized such that the initial residual is made equal to an unknown target signal.At each iteration, five steps are performed.− Proxy of the residual is formed from the current samples and largest components of proxy are located.− The largest components of the proxy signal are identified and merged with the previous support.− The target signal is approximated by solving the least squares − From the approximation, only the largest components are retained, and a new approximation is obtained.− Finally, the samples are updated depending on the residual.The algorithm for CoSaMP is shown in algorithm 1.

Fast digital curvelet transform (FDCT) via wrapping
The idea of using curvelets is explained in this section.Curvelet transform is applied on the image for extracting sparse coefficients [27].Figure 1 shows curvelet tiling of space and frequency.Figure 1(a) represents the induced tiling of the frequency plane.In Fourier space representation, curvelets are supported near a "parabolic" wedge.The shaded area represents such a generic wedge.Figure 1(b) represents the spatial Cartesian grid associated with a given scale and orientation.
FDCTs are faster and less redundant than existing methods.There are two versions of the FDCTs.They are i) unequally-spaced fast Fourier transforms (USFFT) and ii) via wrapping.The difference between FDCT and USFFT is the way in which curvelets at a given scale and angle are translated with respect to each other.In USFFT, the translation grid is tilted to be aligned with the orientation of the curvelet.In wrapping, the grid is the same for every angle within each quadrant but each curvelet is given the proper orientation.The FDCT via wrapping accomplishes better accuracy because of the exact numerical tightness of the digital transform.The wrapping approach considers spatial grid to translate curvelets at each scale and angle.FDCT wrapping is employed for extracting sparse coefficients.Figure 2 depicts the architecture of the FDCT via wrapping.
The parallelogram is the tile  , that contains the frequency support of the curvelet.The gray parallelograms are the replicas resulting from periodization.The rectangle is centered at the origin.The wrapped ellipse appears like "broken into pieces".The flow diagram of the complete processing from input image to reconstructed image is shown in Figure 3.The given input image is resized to 64×64.FDCT is applied to the resized image.Output coefficients which are less than threshold is set to zero to obtain sparse matrix.Sampling matrix is initialized.The product of sparse matrix and sampling matrix is obtained and CoSaMp algorithm is applied on the resultant matrix.Inverse curvelet transform is applied to obtain reconstructed image.

EXPERIMENTAL RESULTS
This section demonstrates the results of the proposed method.The algorithm is implemented in MATLAB 2018a on Intel core i5, 7 th generation processor.The image of size 256×256 is taken for experiment purpose.Image is resized to 64×64 by bicubic interpolation.Out of two versions of the FDCT, FDCT via wrapping (command: fdct_wrapping (Image, 1, 2, 3)) is applied on the image for extracting the curvelet coefficients.A cell array of curvelet coefficients C(1,1), C(1,2) and C(1,3) are obtained.C(1,1), C(1,2) and C(1,3) are converted into single dimension matrix and subsequently all these coefficients are concatenated (command: vertcat).Afterwards, smaller coefficients lesser than threshold is set to zero so that few coefficients only are non-zero and most of the coefficients are zero.The threshold is selected based on experimentation.The coefficients resulted are sparse coefficients.The sampling matrix ∅ is initialized to normally distributed random matrix.The number of measurements of the sampling matrix is denoted by 'm'.Total number of coefficients of curvelet transform is denoted as 'N'.Most signals can be recovered when  =  log  [18].The product of sparse matrix and sampling matrix is obtained.Initially, all values of output sparse matrix are initialized to zeroes.CoSaMP algorithm which is described in section 4 is applied.Afterwards, inverse curvelet transform (command: ifdct_wrapping (D, 1, 64, 64)) is applied to obtain final reconstructed image.PSNR for various images are computed.
The experiments are carried out on standard grayscale images namely Barbara, Boat, Lena, Monarch, Parrot and Peppers.The results are compared with the state-of-the-art methods, including alternating direction method of multipliers (ADMM), accelerated proximal gradient method (APGM), stochastic proximal gradient method (SPGM) [12], complex-valued sparse reconstruction network (CSR-Net) [13], deep residual reconstruction (DR2) [14] and multi-scale residual network (MSRNet) [15] as recorded in Table 1.The results of different methods are taken at measurement rate (MR) of 0.25.From Table 1, it can be seen that CoSaMP outperforms previous methods.The visual results of reconstruction by the proposed method are depicted in Figures 4(a The effectiveness of the algorithm on noisy images are done.The image is subjected to salt and pepper noise of various noise density of 0.01, 0.02, 0.05, 0.08, 0.12, and 0.18.The result of applying the CoSaMP algorithm on the noisy images is shown in Figure 5.It can be observed that noise is removed from the noisy images after reconstruction.As the noise density increases, the quality of the reconstructed image also becomes lesser. Int J Elec & Comp Eng ISSN: 2088-8708  Image reconstruction through compressive sampling matching pursuit … (Shashi Kiran Seetharamaswamy) 6279 images.

ISSN: 2088- 8708 Figure 1 .
Figure 1.Curvelet tiling of space and frequency (a) tiling of the frequency plane and (b) spatial cartesian grid

Figure 2 .
Figure 2. Wrapping of data from a parallelogram to a rectangle by periodicity

Figure 3 .
Figure 3. Flow diagram of the process ) to 4(l).Here, ground truth image and reconstructed images of Barbara, Boat, Lena, Monarch, Parrot and Peppers are shown.It can be seen that the reconstructed images are of good quality.

Table 1 .
Comparison of PSNR for different images for various methods