pit --- Pictures: Iterative Thresholding algorithms

A compressed sensing based image recovery demonstration

A simple and fast implementation of the iterative hard and soft threosholding algorithms (IHT and ISTA) for image recovery. It provides:

Reconstructions of images from few of their pixels (masked images) (see the Jupyter notebook demoIHT.ipynb for a working example)
It relies on the images being sparse under
- The fast discrete cosine transform (DCT) (based on scipy.fftpack.dct)
- The fast wavelet transform (WT) (based on pywt.wavedec2)
Possible thresholding operations:
- Hard thresholding (IHT)
- Soft thresholding (ISTA)
A demonstration of the functioning of image compression is included (Jupyter notebook demo_image_compression.ipynb)

How to run it

A simple example can be found in the Jupyter notebook demoIHT.ipynb.

Description of what it does

Given a picture (left image) we remove a large number of its pixels (middle image). Then we can reconstructs the image (right image) from the middle one; see the Jupyter notebook demo_image_compression.ipynb on how it is done.

Given picture	95% of its pixels removed	Reconstruction

How it works (theory)

Images are compressible (see demo_image_compression.ipynb for an illustration). Effectively, the reconstruction algorithms of pit search for an image in the space of compressed images that is compatible with the given pixels.

Mathematical description

Let T be an invertible transformation taking images to vectors so that the vectors are sparse (i.e., many coefficients are zero) and let iT be the inverse of T. Examples for such a T are the DCT and many versions of the WT. Moreover, let TO be a thresholding operator (hard or soft thresholding).

Suppose we are given a subset of pixels Xsub of an image Xorig and the indices mask of these pixels. Then the reconstruction algorithm pit.estimate essentially does the following iteration:

X = Xsub
repeat
    x = TO( T(X) )
    X = iT(x)
until a stopping criterion is met

For certain transformations T this algorithm essentially solves the problem:

minimize     norm( T(X), L1 )
subject to   norm( X(mask) - Xsub, Frobenius) <= eta

In order to rigorously guarantee this procedure to work T needs to satisfy certain properties. If X is sparse and there are enough pixels in Xsub (depending on the sparsity of X) then the minimum is attained for Xrec that is eta-close to X in Frobenius norm.

Acknowledgment

I thank Stephan Wäldchen for discussions. This project has been funded by the National Science Centre, Poland (Polonez 2015/19/P/ST2/03001), i.e., This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 665778.

References

S. Foucart and H. Rauhut, A mathematical introduction to compressive sensing (Birkhäuser, 2013)
J. Bobin, J.-L. Starck, and R. Ottensamer, Compressed Sensing in Astronomy, arXiv:0802.0131
A. Beck and M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM J. Imaging Sci., 2(1), 183–202 (2009)

Better visualization

The demo image in the README suggests that the image can be reconstructed from almost nothing.

This perception is largely caused by the way that the results are presented:

The missing pixels are shown as white, mostly drowning out the signal for a human observer.
The image has a high resolution, but has been downscaled by a large factor, so the sparse pixels are barely visible.
The reconstructed image has been downscaled by a large factor, hiding artifacts.

For example, compare this simple nearest neighbor reconstruction with the wavelet-based reconstruction. Without zooming in, they look about the same.

KNN	Wavelet

KNN sparse reconstruction baseline

from PIL import Image
import numpy as np
from scipy.spatial import KDTree

image = np.array(Image.open("marie_curie.jpg").convert("L"))
h, w = image.shape
x = np.arange(w)
y = np.arange(h)
x, y = np.meshgrid(x, y)
is_known = np.random.rand(h, w) < 0.05
is_unknown = np.logical_not(is_known)
image[is_unknown] = 0
xy = np.stack([x, y], axis=2)
data = xy[is_known]
query = xy.reshape(-1, 2)
distances, indices = KDTree(data).query(query, k=5)
values = image[is_known][indices]
weights = np.exp(-0.5 * distances)
values = np.sum(weights * values, axis=1) / np.sum(weights, axis=1)
reconstructed = values.reshape(h, w)
reconstructed = np.clip(reconstructed, 0, 255).astype(np.uint8)
reconstructed = Image.fromarray(reconstructed)
reconstructed.show()
reconstructed.save(f"marie_curie_reconstructed.jpg")

I would recommend to:

Demonstrate the wavelet reconstruction on an image with a lower resolution. This way, the images will not be downscaled when displayed in a figure.
Compare to a simple baseline to show that wavelet reconstruction is better.

martkl / cs_image_recovery_demo Goto Github PK

cs_image_recovery_demo's Introduction

pit --- Pictures: Iterative Thresholding algorithms

A compressed sensing based image recovery demonstration

How to run it

Description of what it does

How it works (theory)

Mathematical description

Acknowledgment

References

cs_image_recovery_demo's People

Contributors

Stargazers

Watchers

Forkers

cs_image_recovery_demo's Issues

KNN sparse reconstruction baseline

Recommend Projects

Recommend Topics

Recommend Org