Compressed Sensing Image Reconstruction via Recursive Spatially Adaptive Filtering
nonparametric stochastic approximation approach

PAGE UNDER CONSTRUCTION

Introduction

We introduce a new approach to image reconstruction from highly incomplete data. The available data are assumed to be a small collection of spectral coefficients of an arbitrary linear transform. This reconstruction problem is the subject of intensive study in the recent field of "compressed sensing" (also known as "compressive sampling"). Our approach is based on a quite specific recursive filtering procedure. At every iteration the algorithm is excited by injection of random noise in the unobserved portion of the spectrum and a spatially adaptive image denoising filter, working in the image domain, is exploited to attenuate the noise and reveal new features and details out of the incomplete and degraded observations. This recursive algorithm can be interpreted as a special type of the Robbins-Monro stochastic approximation procedure with regularization enabled by a spatially adaptive filter. Overall, we replace the conventional parametric modeling used in CS by a nonparametric one. We illustrate the effectiveness of the proposed approach for two important inverse problems from computerized tomography: Radon inversion from sparse projections and limited-angle tomography. In particular we show that the algorithm allows to achieve exact reconstruction of synthetic phantom data even from a very small number projections. The accuracy of our reconstruction is in line with the best results in the compressed sensing field.

Results

Simulation of Radon reconstruction from sparse projections (approximating Radon projections as radial lines in FFT domain)

We consider three illustrative inverse problems of compressed sensing for computerized tomography. In particular, we show reconstruction examples of the Shepp-Logan phantom from sparse Radon projections, with 22 and 11 radial lines in FFT-domain (i.e., available Radon projections), and reconstruction from limited-angle projections, with a reduced subset of 61 projections within a 90 degrees aperture. The available portions of the spectrum and the initial back-projection estimates are shown below. As the recursive algorithm progresses, the reconstruction error improves steadily until numerical convergence; for all three cases the reconstruction is exact, with the estimate converging to the original unobserved phantom.

- Sparse projections: 22 radial lines


available portion of the spectrum (22 radial lines)	Back-projection estimate		Estimate after convergence (exact reconstruction)

Exact reconstruction from 22 radial lines (ANIMATION)

- Sparse projections: 11 radial lines


available portion of the spectrum (11 radial lines)	Back-projection estimate		Estimate after convergence (exact reconstruction)

Exact reconstruction from 11 radial lines (ANIMATION)

- Limited-angle


available portion of the spectrum (90 degrees aperture, 61 radial lines)	Back-projection estimate		Estimate after convergence (exact reconstruction)

Exact reconstruction from 90 degrees limited-angle with 61 radial lines (ANIMATION)

Notes:
The animation files are rather big (about 6 Mbytes). They are in .GIF format.
Because of their length, it is recommended to download the files (right click, Save target as...) and open them with an external viewer (e.g., Irfanview).
In the animations, the three numbers on the top of the frames are the PSNR(dB), the iteration count, and the maximum absolute error (l-infinity norm).
The animations show the image estimate and its difference from the original (unknown) image. This difference is normalized with respect to the maximum absolute error.

Reconstruction from low-frequency portion of Fourier spectrum

We present the result of reconstruction of the Cameraman image (256x256 pixels) from the low-frequency portion of its Fourier spectrum (a 128×128 square centered at the DC).


available portion of the spectrum (low frequencies)	Back-projection estimate (zero-padding the high frequencies)		Estimate after 62 iterations

For all the above experiments, the block-matching and 3D filtering algorithm (BM3D) served as the spatially adaptive filter used in the recursive stochastic approximation algorithm. The separable 3D Haar wavelet decomposition is adopted as the transform utilized internally by the BM3D algorithm.