Compressed sensing of parameterized shapes in drawings

Suppose we have an image \(f\), aranged as a column vector of length \(L_1L_2\), where \(L_i\) are the image dimensions. We can write \(f=H p\), where \(H\) is the "inverse" Hough transform, and \(p\) is a parameter space selector. We can form the \(k^{th}\) column of \(H\) by computing the image of the shape given by the parameters \(\pi_k\).

We measure \(y = \Phi f = \Phi H p\). The measurement we take is a projection of \(f\) onto a set of \(\phi_m, m=\{1,2,...,M\}\) basis vectors. In practice, these vectors are random. One choice is using random Gaussian vectors. Another is simply drawing entries from \(\pm 1\) with probability one half. Gurbuz 2008 shows that these two approaches perform significantly better than random pixel measurement (as expected by coherence predictions).

For the noiseless case, we can simply minimize \(|p|_1\) subject to \(y=A p\) (where \(A=\Phi H\)), which can be phrased as a linear program. Any real world measurement will have noise, so we must address this in our solution. Suppose we have \(y=Ap+z\), where \(z \sim \mathcal{N}(0,\sigma^2)\). One way is simply to relax the linear equality constraint by a user-supplied factor \(\epsilon\), so \(|y-Ap|_2 < \epsilon\). This will recover the solution, but is a second order cone problem as opposed to a linear program.

If we wish to continue solving a linear program, we can use the Dantzig Selector, minimizing \(|p|_1\) subject to \(|A^\top (y-Ap)|_\infty < \sigma\sqrt{2\log{N}}\). This constraint says that the residual of a candidate solution must not be very correlated with any column of \(A\). This can be phrased as a linear program.

There is a nice framework for incorporating multiple shape types into a single problem. This essentially amounts to stacking the vectors.

LPs can be phrased in common language (c.f. and solved efficiently by available packages. SOCPs also. I've done this process of LP before, so maybe I'll stick to the Dantzig selector.

All that remains in my head is figuring out how to efficiently calculate the matrix \(A\). Each column of \(H\) describes a binary image, so we can just pick elements from the random gaussian vectors...

In background reading, I found a paper doing something similar, and used much of their derivations:

In [12]:
%pylab inline
from numpy import *
#from scipy import ndimage
from scipy import misc
from numpy.linalg import matrix_rank #for debug only
Populating the interactive namespace from numpy and matplotlib

Getting a hand drawn image

In [93]:
lines = misc.imread('lines.JPG')
lines = .2126*lines[...,0] + .7152*lines[...,1] + .0722*lines[...,2]
lines = lines.astype(uint8)
lines = 255-lines #invert so lines are white
N = 160
lines = lines[0:0+N,80:80+N]
lines = where(lines>128,255,0)
<matplotlib.figure.Figure at 0x11231b810>

Generating the Hough basis

In [72]:
def find_valid_coords(r,th,N,M):
    if -pi/4<=th<=pi/4:
        Y = arange(N,dtype=int32)
        X = (r*cos(th) - (Y-r*sin(th))*tan(th)).astype(int32)
        X = hstack((X-1,X))
        Y = hstack((Y,Y))
        #X = hstack((X-1,X,X+1))
        #Y = hstack((Y,Y,Y))
        valid = where(logical_and(M>X,X>=0))[0]
    elif -pi/2<=th<=pi/2:
        X = arange(M,dtype=int32)
        Y = (r*sin(th) + (X-r*cos(th))/tan(-th)).astype(int32)
        X = hstack((X,X))
        Y = hstack((Y-1,Y))
        #X = hstack((X,X,X))
        #Y = hstack((Y-1,Y,Y+1))
        valid = where(logical_and(N>Y,Y>=0))[0]
    return X[valid], Y[valid]

def make_h(r,th,N,M):
    Make the NxM pixel image vector corresponding to Hough parameters r,th
    h = zeros((N,M),int32)
    X,Y = find_valid_coords(r,th,N,M)
    h[Y,X] = 255
    return h
def h_k(p):
    return make_h(p[0],p[1],N,N)
def h_endpts(r,th,N,M):
    X,Y = find_valid_coords(r,th,N,M)
    ends = array([0,-1])
    return X[ends],Y[ends]

def pick_top_n(signal,n):
    #select the top n values from signal, return indices
    #this is wasteful, but fast to write
    return argsort(signal)[::-1][:n]
In [98]:
#generate the Hough grid for the problem.  We need to make sure each shape has adequate
#support.  As a heuristic, we only keep shapes with 50 valid pixels.
N_r = 60
N_th = 60
pi_r = linspace(-sqrt(2)*N,sqrt(2)*N,N_r)
pi_th =  linspace(-pi/2,pi/2,N_th)
#PI = dstack(meshgrid(pi_r,pi_th)).reshape(-1,2)
PI = [] #array of r,th values
#pic = zeros((N,N),uint8)
for th in pi_th:
    for ri in pi_r:
        new =  make_h(ri,th,N,N)
        if shape(where(new!=0)[0])[0]>.2*N:
            #pic +=new
PI = array(PI)
M = shape(PI)[0] #length of decision variable x
print M

Constructing A and b

In [99]:
#now we construct A=phi*H.  Let's use the +-1 approach for phi to keep working with integers
K=400 #number of random projections
A = zeros((K,M),int32)
#Phi = 2*randint(2, size=(K,N*N))-1 # +-1 with uniform probability
Phi = randn(K,N*N) # Gaussian, zero mean, unit variance
for j in range(M):
    A[:,j] = dot( Phi, h_k(PI[j]).reshape(-1) )
from numpy.linalg import matrix_rank
print matrix_rank(A), K #are we full rank?
400 400

In [100]:
#now we measure the image in question
def sample_image(image):
    b = zeros((K))
    for i in range(K):
        b[i] = dot( Phi[i], image.reshape(-1) )
    return b
In [103]:
plot(sample_image(make_h(40,.4,N,N) + make_h(6,-.51,N,N) + make_h(20,1.57,N,N)));
plt.title("%d image measurements"%K);

Recast as Linear Programming

The \(L_1\) minimization subject to a linear constraint can be recast as the linear program:

minimize \(\sum u_i\) subject to - \(x-u \leq 0\) - \(-x-u \leq 0\) - \(Ax=b\)

This is a noiseless optimization, but the Dantzig selector can also be cast in a very similar linear programming form:

minimize \(\sum u_i\) subject to - \(x-u \leq 0\) - \(-x-u \leq 0\) - \(A^\top (Ax-b)-\epsilon \leq 0\) - \(-A^\top (Ax-b)-\epsilon \leq 0\)

In [83]:
#now, let's try solving the noiseless case as a linear program.
#to recast as LP, we use derivation of
#guide to cvxopt.solvers.lp is available:
from cvxopt import matrix as cvxmatrix
from cvxopt import solvers
def noiseless_hough_cs(b):
    #decision variable has length 2*M
    c = cvxmatrix(hstack((zeros(M),ones(M))))
    G = cvxmatrix( vstack((
        hstack(( diag( ones(M)), diag(-ones(M)) )),
        hstack(( diag(-ones(M)), diag(-ones(M)) ))
    h = cvxmatrix(zeros(2*M))
    Acvx = cvxmatrix(hstack((A,zeros((K,M)) )) )
    return solvers.lp(c,G,h,Acvx, cvxmatrix(b) )

def plot_sol(sol,image):
    pi_sol = array(sol['x'][:M])[:,0]**2
    #thresh = .25*amax(pi_sol)
    #XX = where(pi_sol>thresh)[0]
    print shape(pi_sol)
    XX = pick_top_n(pi_sol,3)
    print XX
    thresh = pi_sol[XX[2]]
    gca().set_title("Input image and Reconstruction")
    for PII in PI[XX]:
        endpts = h_endpts(PII[0],PII[1],N,N)
    #imshow( sum(dstack([h_k(PII) for PII in PI[XX]]),axis=-1) )
In [85]:
perfect_lines = make_h(40,.4,N,N) + make_h(6,-.51,N,N) + make_h(20,1.57,N,N)
b_perfect = sample_image(perfect_lines)
sol_perfect = noiseless_hough_cs(b_perfect)
print "%d pixels, %d measurements (%.3f), %d hough shape functions"%(N*N,K,float(K)/N/N,M)
     pcost       dcost       gap    pres   dres   k/t
 0:  0.0000e+00 -0.0000e+00  2e+03  7e+01  1e-16  1e+00
 1:  1.7017e+00  1.6998e+00  3e+02  1e+01  2e-16  2e-01
 2:  7.1883e+00  7.1854e+00  1e+02  5e+00  4e-16  7e-02
 3:  1.4151e+01  1.4149e+01  4e+01  2e+00  5e-16  2e-02
 4:  1.8481e+01  1.8480e+01  2e+01  6e-01  7e-16  7e-03
 5:  1.9902e+01  1.9901e+01  7e+00  3e-01  8e-16  3e-03
 6:  2.0463e+01  2.0463e+01  3e+00  1e-01  9e-16  1e-03
 7:  2.0685e+01  2.0685e+01  2e+00  6e-02  1e-15  6e-04
 8:  2.0825e+01  2.0825e+01  6e-01  2e-02  1e-15  2e-04
 9:  2.0872e+01  2.0872e+01  3e-01  1e-02  2e-15  1e-04
10:  2.0900e+01  2.0900e+01  1e-01  4e-03  3e-15  3e-05
11:  2.0907e+01  2.0907e+01  5e-02  2e-03  6e-15  2e-05
12:  2.0912e+01  2.0912e+01  1e-02  5e-04  7e-15  4e-06
13:  2.0914e+01  2.0914e+01  3e-03  1e-04  2e-14  1e-06
14:  2.0914e+01  2.0914e+01  1e-03  5e-05  3e-14  5e-07
15:  2.0914e+01  2.0914e+01  5e-04  2e-05  3e-14  2e-07
16:  2.0914e+01  2.0914e+01  2e-05  9e-07  5e-14  8e-09
Terminated (singular KKT matrix).
[518 933  34]
2500 pixels, 400 measurements (0.160), 1506 hough shape functions

In [86]:
perfect_lines = make_h(40,.4,N,N) + make_h(6,-.51,N,N) + make_h(20,1.57,N,N) #+ make_h(40,1.,N,N)
perfect_lines += .3*255*randn(N,N) #add noise at .3 of signal
b_perfect = sample_image(perfect_lines)
sol_perfect = noiseless_hough_cs(b_perfect)
print "%d pixels, %d measurements (%.3f), %d hough shape functions"%(N*N,K,float(K)/N/N,M)
     pcost       dcost       gap    pres   dres   k/t
 0:  0.0000e+00 -0.0000e+00  2e+03  7e+01  1e-16  1e+00
 1:  4.2912e+00  4.2891e+00  3e+02  1e+01  3e-16  2e-01
 2:  1.9213e+01  1.9210e+01  1e+02  5e+00  4e-16  6e-02
 3:  3.0961e+01  3.0959e+01  5e+01  2e+00  7e-16  2e-02
 4:  3.4770e+01  3.4769e+01  3e+01  9e-01  7e-16  1e-02
 5:  3.6820e+01  3.6819e+01  1e+01  4e-01  8e-16  4e-03
 6:  3.7607e+01  3.7607e+01  5e+00  2e-01  8e-16  2e-03
 7:  3.8060e+01  3.8060e+01  2e+00  7e-02  1e-15  7e-04
 8:  3.8193e+01  3.8193e+01  1e+00  4e-02  1e-15  4e-04
 9:  3.8289e+01  3.8289e+01  3e-01  1e-02  2e-15  1e-04
10:  3.8319e+01  3.8319e+01  1e-01  4e-03  4e-15  4e-05
11:  3.8330e+01  3.8330e+01  3e-02  1e-03  6e-15  9e-06
12:  3.8333e+01  3.8333e+01  7e-03  3e-04  4e-14  2e-06
13:  3.8333e+01  3.8333e+01  3e-03  1e-04  1e-13  1e-06
14:  3.8334e+01  3.8334e+01  2e-03  6e-05  9e-14  6e-07
15:  3.8334e+01  3.8334e+01  4e-04  1e-05  1e-13  1e-07
16:  3.8334e+01  3.8334e+01  1e-05  4e-07  8e-14  4e-09
Terminated (singular KKT matrix).
[ 518 1471  908]
2500 pixels, 400 measurements (0.160), 1506 hough shape functions

Dantzig Selector

Didn't have time to make sure this is right. Looks like redundant constraints...

In [17]:
def dantzig_hough_cs(b,eps):
    #decision variable has length 2*M
    AA = dot(A.T,A); Ab = dot(A.T,b)
    c = cvxmatrix(hstack((zeros(M),ones(M))))
    G = cvxmatrix( vstack((
        hstack(( diag(ones(M)), -diag(ones(M)) )),
        hstack((-diag(ones(M)), -diag(ones(M)) )),
        hstack(( AA,     zeros((M,M)) )),
        hstack((-AA,     zeros((M,M)) ))
    h = cvxmatrix( hstack((zeros(2*M), Ab+eps*ones(M),-Ab+eps*ones(M))) )
    print matrix_rank(asarray(G)), shape(asarray(G))
    return solvers.lp(c,G,h)
In [19]:
perfect_lines = make_h(10,.5,N,N) + make_h(6,-.51,N,N) + make_h(20,1.57,N,N)
b_perfect = sample_image(perfect_lines)
sol_perfect = dantzig_hough_cs(b_perfect,0.)
print "%d pixels, %d measurements (%.3f), %d hough shape functions"%(N*N,K,float(K)/N/N,M)
3092 (6184, 3092)

ValueError                                Traceback (most recent call last)
<ipython-input-19-e11a5dd2f968> in <module>()
      1 perfect_lines = make_h(10,.5,N,N) + make_h(6,-.51,N,N) + make_h(20,1.57,N,N)
      2 b_perfect = sample_image(perfect_lines)
----> 3 sol_perfect = dantzig_hough_cs(b_perfect,0.)
      4 #plot_sol(sol_perfect,perfect_lines)
      5 print "%d pixels, %d measurements (%.3f), %d hough shape functions"%(N*N,K,float(K)/N/N,M)

<ipython-input-17-c2f3bf5860b6> in dantzig_hough_cs(b, eps)
     11     h = cvxmatrix( hstack((zeros(2*M), Ab+eps*ones(M),-Ab+eps*ones(M))) )
     12     print matrix_rank(asarray(G)), shape(asarray(G))
---> 13     return solvers.lp(c,G,h)

/usr/local/lib/python2.7/site-packages/cvxopt/coneprog.pyc in lp(c, G, h, A, b, solver, primalstart, dualstart)
   3018     return conelp(c, G, h, {'l': m, 'q': [], 's': []}, A,  b, primalstart,
-> 3019         dualstart)

/usr/local/lib/python2.7/site-packages/cvxopt/coneprog.pyc in conelp(c, G, h, dims, A, b, primalstart, dualstart, kktsolver, xnewcopy, xdot, xaxpy, xscal, ynewcopy, ydot, yaxpy, yscal)
    685         try: f = kktsolver(W)
    686         except ArithmeticError:
--> 687             raise ValueError("Rank(A) < p or Rank([G; A]) < n")
    689     if primalstart is None:

ValueError: Rank(A) < p or Rank([G; A]) < n

Using hand-drawn lines

Need to get the noise model working before expecting this to work.

In [180]:
b = sample_image(lines)
sol = noiseless_hough_cs(b)
     pcost       dcost       gap    pres   dres   k/t
 0:  0.0000e+00 -0.0000e+00  8e+02  4e+01  1e-16  1e+00
 1:  4.6093e+00  4.6083e+00  1e+02  6e+00  2e-16  1e-01
 2:  2.2520e+01  2.2518e+01  3e+01  2e+00  1e-15  4e-02
 3:  2.9479e+01  2.9478e+01  1e+01  6e-01  1e-15  1e-02
 4:  3.1583e+01  3.1583e+01  5e+00  3e-01  1e-15  5e-03
 5:  3.2633e+01  3.2632e+01  2e+00  8e-02  1e-15  2e-03
 6:  3.2909e+01  3.2909e+01  6e-01  3e-02  2e-15  5e-04
 7:  3.3002e+01  3.3002e+01  3e-01  2e-02  2e-15  3e-04
 8:  3.3054e+01  3.3054e+01  1e-01  6e-03  4e-15  9e-05
 9:  3.3076e+01  3.3076e+01  3e-02  2e-03  4e-15  3e-05
10:  3.3083e+01  3.3083e+01  1e-02  6e-04  1e-14  8e-06
11:  3.3085e+01  3.3085e+01  5e-03  2e-04  2e-14  4e-06
12:  3.3086e+01  3.3086e+01  1e-03  5e-05  3e-14  8e-07
13:  3.3086e+01  3.3086e+01  2e-04  9e-06  4e-14  1e-07
14:  3.3086e+01  3.3086e+01  7e-06  4e-07  9e-14  5e-09
Terminated (singular KKT matrix).

Think about noise model!

In []: