This question deals with an example of image reconstruction related to this other question on signal reconstruction. I have different issues in both examples but there could be an underlying factor.

I tried to reconstruct an image using compressed sensing just as in Coursera's course "Computational Methods for Data Analysis" and described in some detail in this pdf (page 414.) By the way, this is the low resolution image:

enter image description here

However, I have a question regarding the reconstruction of the image using $f \approx \Psi x$ where $f$ is the original image, $\Psi$ is a basis (in this case, a discrete cosine transform) to express the image and $x$ its coefficients. This is what I have done, but I'm not entirely sure it's correct.

from PIL import Image
from scipy.fftpack import dct, idct
from sklearn.linear_model import Lasso

# Loading image in grayscale and obtaining its dimensions
im = Image.open('coffee-cup.jpg').convert('L')
nx, ny = im.shape
# Number of sample points used to reconstruct image
k = 1000

# Create a permutation from 1 to nx*ny and choose the first k numbers
Atest = zeros((nx, ny)).reshape(1, nx*ny)
r1 = permutation(arange(1, nx*ny))
r1k = r1[1:k+1]

# Suppose the first number in the permutation is 42. Then choose the 42th
# elements in a matrix full of zeros and change it to 1. Apply a discrete
# cosine transformation to this matrix, and add the result to other matrix
# Adelta as a new row.
Adelta = zeros((k, nx*ny))
for i, j in enumerate(r1k):
    Atest[0,j] = 1
    Adelta[i, :] = dct(Atest)
    Atest[0, j] = 0

# Use the same permutation to sample the image to be reconstructed
image1 = im.reshape(nx*ny,1)
b = image1[r1k]

# Solve the optimization problem Adelta * x = b
lasso = Lasso(alpha=0.001)

After this, lasso.coef_ contains the coefficients $x$. This is the part I'm not sure about. I transformed the coefficients using the discrete cosine transform. However, the construction of Adelta seems to suggest other more elaborated $\Psi$. Yet, when I plot the reconstruction using the inverse discrete cosine transform (IDCT) and the discrete cosine transform (DCT), this is what I get:

recons2 = dct(lasso.coef_).reshape((nx,ny))
recons = idct(lasso.coef_).reshape((nx,ny))

fig = figure()
ax = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

enter image description here

It worked reasonably well, but notice that the first image (in which IDCT was used) performed better. I would have assumed that, if something, DCT was more appropriate choice since at least we used the result of a DCT in the construction of Adelta. What is the correct procedure here?


2 Answers 2


This article helped me to understand this problem by using the Kronecker product. Here is an excerpt of the article:

Creating the $A$ matrix for 2D image data takes a little more ingenuity than it did in the 1D case. In the derivation that follows, we’ll use the Kronecker product $\otimes$ and the fact that the 2D discrete cosine transform is separable to produce our operator $A$.

Let $X$ be an image in the spectral domain and $D_i=\textit{idct}(I_i)$, where $I_i$ is the identity matrix of size $i$. Then:

$$\begin{align*} \textit{idct2}(X)&=\textit{idct}\left(\textit{idct}\left(X^T\right)^T\right)\\ &=D_m\left(D_nX^T\right)^T\\ &=D_mXD_n^T \end{align*}$$

If $\textit{vec}(X)$ is the vector operator that stacks columns of $X$ on top of each other, then:

$$\begin{align*} \textit{vec}\left(D_mXD_n^T\right)&=\left(D_n\otimes D_m\right)\textit{vec}(X)\\ &=\left(D_n\otimes D_m\right)x\quad\text{where}\ \text{and}\ x\equiv\textit{vec}(X) \end{align*}$$

In the article the library CVXPY is used, but I solved it with the LASSO regularizer:

import numpy as np
import scipy.fftpack as spfft
import scipy.ndimage as spimg
from sklearn.linear_model import Lasso

# read original image and downsize for speed
X = spimg.imread('image.jpg', flatten=True, mode='L') # read in grayscale
ny, nx = X.shape

# extract small sample of signal
k = round(nx * ny * 0.5) # 50% sample
ri = np.random.choice(nx * ny, k, replace=False) # random sample of indices
b = X.T.flat[ri]
b = np.expand_dims(b, axis=1)

# create dct matrix operator using kron (memory errors for large ny*nx)
A = np.kron(
    spfft.idct(np.identity(nx), norm='ortho', axis=0),
    spfft.idct(np.identity(ny), norm='ortho', axis=0)
A = A[ri,:] # same as phi times kron

Now you can apply the LASSO regression

def idct2(x):
    return spfft.idct(spfft.idct(x.T, norm='ortho', axis=0).T, norm='ortho', axis=0)
lasso = Lasso(alpha=0.001)
lasso.fit(A, b)

Xat = np.array(lasso.coef_).reshape(nx, ny).T # stack columns
# Get the reconstructed image
Xa = idct2(Xat)

I know that this is an old question but I am answering it in case googles this topic and would like to know the answer.

The problem is the compression matrix and the basis you are doing the sparse reconstruction in are coherent since they are both constructed from the DCT matrix. Do the sparse reconstruction in the DCT basis since the image will sparse/compressible in that basis. For the compression matrix, you need to choose something that will not correlate with the DCT matrix. Noiselets, Random Gaussian will work well. You might try JL-Lemma type matrices as well though I never found them to be very reliable.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.