Enhanced Compressed Sensing Recovery with Level Set Normals

by

 

Code for compressed sensing recovery of images.

Dgraph(image, graphfile)
function Dgraph(image, graphfile)

verb = 0;
windowing = 1;

NIds_read = 0;
[Ny,Nx,L] = size(image);
szW = round(sqrt(Ny*Nx)/15); % the window will ne of size 2*szW+1 
file = ['graph_data/' graphfile '.mat'];
if nargin<2
    graphfile = 'temp';
end;
if( exist('graph_data')<1 )
    mkdir('graph_data');
end;
if windowing
%     %%%%%%%% read neighbours from previous iteration and update value on w, but not the neighbours 
%     if( exist(['graph_data/' graphfile '.mat'])==2 )
%         disp('read previous NIds');
%         load(file,'NIds');
%     end;
%     if( exist(['graph_data/NIds_' int2str(Ny) '_' int2str(Nx) '_' int2str(szW)])==2 )
%         load(['graph_data/NIds_' int2str(Ny) '_' int2str(Nx) '_' int2str(szW)]);
%         NIds_read = 1;
%     end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% collect patches
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verb
    disp('collecting patches ...');
    tic
end;

% patch size = 2m+1 x 2m+1
m = 2; 
%m = 2; 
s = 2*m+1;

% extend image to deal with boundaries
u = image;
uext = zeros(Ny+2*m,Nx+2*m,L);
uext(m+1:m+Ny,m+1:m+Nx,:) = u;
uext(1:m,:,:) = uext(m+m:-1:m+1,:,:);
uext(m+Ny+1:m+Ny+m,:,:) = uext(m+Ny:-1:m+Ny+1-m,:,:);
uext(:,1:m,:) = uext(:,m+m:-1:m+1,:);
uext(:,m+Nx+1:m+Nx+m,:) = uext(:,m+Nx:-1:m+Nx+1-m,:);

if not(exist(['graph_data/Tp_ids_' int2str(Ny) '_' int2str(Nx) '_' int2str(m)])==2)    
    Tp_ids = zeros(Ny,Nx,s^2);
    % first m rows, the window goes from row 1:m
    [X,Y] = meshgrid(1:s,1:s);
    window = sub2ind([Ny+2*m Nx+2*m],X(:),Y(:));
    window = reshape(window,1,s^2);
    for x=1:Nx
        for y=1:Ny
            Tp_ids(y,x,:) = window + y - 1;
        end     
        window = window + Ny + 2*m;
    end
    Tp_ids = reshape(Tp_ids,Ny*Nx,s^2); 
    save(['graph_data/Tp_ids_' int2str(Ny) '_' int2str(Nx) '_' int2str(m)],'Tp_ids');
    clear window
else
    load(['graph_data/Tp_ids_' int2str(Ny) '_' int2str(Nx) '_' int2str(m)],'Tp_ids');
end;
Tp = zeros(Ny*Nx,L*s^2);
for l=1:L
    ul = reshape(uext(:,:,l),Ny+2*m,Nx+2*m);
    Tp(:,(l-1)*(s^2) + (1:s^2)) = reshape(ul(Tp_ids(:)),Nx*Ny,s^2);
end;
save(file,'Tp','m','s','image');

if not(NIds_read)    
    wside = 2*szW+1;
    NIds = zeros(Ny,Nx,wside^2);
    % first szW rows, the window goes from row 1:wside
    [X,Y] = meshgrid(1:wside,1:wside);
    window = sub2ind([Ny Nx],X(:),Y(:));
    window = reshape(window,wside^2,1);
    for x=1:(szW+1)
        for y=1:(szW+1)      
            NIds(y,x,:) = window;
        end
        for y=(szW+2):(Ny-szW-1)
            NIds(y,x,:) = window + y - szW;
        end
        for y=(Ny-szW):Ny
            NIds(y,x,:) = window + Ny - wside;
        end        
    end
    % middle rows, window centered at pixel with szW rows up nad szW rows down
    for x=(szW+2):(Nx-szW)
        window = window + Ny;
        for y=1:(szW+1)     
            NIds(y,x,:) = window;
        end
        for y=(szW+2):(Ny-szW-1)
            NIds(y,x,:) = window + y - szW;
        end
        for y=(Ny-szW):Ny
            NIds(y,x,:) = window + Ny - wside;
        end        
    end  
    % last szW rows, the window goes from row (Nx-wside):Nx
    for x=(Nx-szW+1):Nx
        for y=1:(szW+1)  
            NIds(y,x,:) = window;
        end
        for y=(szW+2):(Ny-szW-1)
            NIds(y,x,:) = window + y - szW;
        end
        for y=(Ny-szW):Ny
            NIds(y,x,:) = window + Ny - wside;
        end        
    end    
    NIds = reshape(NIds,Ny*Nx,wside^2);
    save(['graph_data/NIds_' int2str(Ny) '_' int2str(Nx) '_' int2str(szW)],'NIds');
end;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute weight graph
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verb
    toc
    disp('computing the weight on the graph ...');
end;

opts.kNN=7; % opts.kNN is the number of neighbors
opts.alpha=5; % opts.alpha is a universal scaling factor
opts.kNNdelta=round(opts.kNN/2); 
opts.quiet = 1;
% opts.kNNdelta specifies which nearest 
% neighbor determines the local scale
if windowing
    [w NIds NNDist] = fgf(Tp,opts,NIds);
    save(file,'NIds','-append');
else
    w = fgf(Tp,opts);
end;
% distances on weights are given by similarity of patches
% w is a sparse symmetric matrix
save(file,'w','-append');

% wv=dn(w,'gph'); % normalize the matrix
% [V E]=eigs(wv);
% 
% save(file,'wv','V','E','-append');








%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute gradient operator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verb
    toc
    disp('computing gradients ...');
    tic
end;
wt = triu(w,1); % upper triangular part of w
w = wt;
[r c v] = find(w);
n = length(r);
Dr = [1:n 1:n];
Dc(1:n) = r;
Dc(n+1:2*n) = c;
Dv(1:n) = 2*sqrt(v);
Dv(n+1:2*n) = -2*sqrt(v);
D = sparse(Dr,Dc,Dv);
Dt = D';
% D NL gradient operator
% Dt NL divergence operator associated to D

% 
rW = r; cW = c; vW = v;
rx = ceil(rW/Ny); ry = rW - ( rx-1 )*Ny;
cx = ceil(cW/Ny); cy = cW - ( cx-1 )*Ny;
theta = atan2(cx-rx,cy-ry);
save(file,'rW','cW','theta','-append');


Dv2(1:n) = ones(1,n);
Dv2(n+1:2*n) = ones(1,n);
D2 = sparse(Dr,Dc,Dv2);
J = D2';
% sumJ = sum(J,2); minsJ = min(sumJ); maxsJ = max(sumJ); medsJ = median(sumJ);
% disp(['number of neighbours between ' num2str(minsJ) ' and ' num2str(maxsJ) ' with median ' num2str(medsJ) ]);
% J adds the components of the gradient associated to the same node i.e if
% d = D*u , then J*d = \sum_{j neighbour i} d_{i,j}
% the result is a vector of the size of the nodes in particular J*(d.*d) is
% the square of the norm of the NL gradient in each node

Dv3(1:n) = ones(1,n);
Dv3(n+1:2*n) = zeros(1,n);
D3 = sparse(Dr,Dc,Dv3);
T = D3;
% T replicates the value of a function in each node to the arcs associated
% to that node, the resuls is a vector of the same size that the arcs in the graph
%i.e if d are the values associated to the arcs and w is a weighting
%function for each node then d.*(T*w) weights each arc with the weight w
%associated to its node.

save(file,'D','Dt','J','T','-append');

if verb
    toc
    spy(D); title('gradient on the graph');pause(0.1);
end;

Contact us