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. http://abel.ee.ucla.edu/cvxopt/userguide/coneprog.html) 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: http://users.ece.gatech.edu/justin/Publications_files/gurbuz08co.pdf
%pylab inline
from numpy import *
#from scipy import ndimage
from scipy import misc
from numpy.linalg import matrix_rank #for debug only
Getting a hand drawn image
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]
plt.set_cmap('gray')
plt.figure(figsize=(8,3.3));
plt.subplot(121)
imshow(lines);
plt.subplot(122)
#threshold
lines = where(lines>128,255,0)
imshow(lines);
Generating the Hough basis
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]
#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=50
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:
PI.append([ri,th])
#pic +=new
#plt.gca().set_aspect(1.)
#imshow(pic)
#plt.set_cmap('gray')
PI = array(PI)
M = shape(PI)[0] #length of decision variable x
print M
Constructing A and b
#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?
#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
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\)
#now, let's try solving the noiseless case as a linear program.
#to recast as LP, we use derivation of
#http://users.ece.gatech.edu/~justin/l1magic/downloads/l1magic.pdf
#guide to cvxopt.solvers.lp is available: http://cvxopt.org/userguide/coneprog.html#linear-programming
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):
plt.figure(figsize=(12,4))
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]]
plt.subplot(121)
plot(pi_sol)
plot([0,M],[thresh,thresh],linestyle='--',color='g')
gca().set_title("Solution")
plt.subplot(122)
imshow(image)
gca().set_title("Input image and Reconstruction")
for PII in PI[XX]:
endpts = h_endpts(PII[0],PII[1],N,N)
plot(endpts[0],endpts[1],c='y',lw=2)
gca().set_xlim([0,N])
gca().set_ylim([0,N])
#imshow( sum(dstack([h_k(PII) for PII in PI[XX]]),axis=-1) )
#gca().set_title("Reconstruction")
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)
plot_sol(sol_perfect,perfect_lines)
print "%d pixels, %d measurements (%.3f), %d hough shape functions"%(N*N,K,float(K)/N/N,M)
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)
plot_sol(sol_perfect,perfect_lines)
print "%d pixels, %d measurements (%.3f), %d hough shape functions"%(N*N,K,float(K)/N/N,M)
Dantzig Selector
Didn't have time to make sure this is right. Looks like redundant constraints...
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)
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.)
#plot_sol(sol_perfect,perfect_lines)
print "%d pixels, %d measurements (%.3f), %d hough shape functions"%(N*N,K,float(K)/N/N,M)
Using hand-drawn lines
Need to get the noise model working before expecting this to work.
b = sample_image(lines)
sol = noiseless_hough_cs(b)
plot_sol(sol,lines)
Think about noise model!