Enhanced Compressed Sensing Recovery with Level Set Normals

by

 

Code for compressed sensing recovery of images.

CS_TV(R, f, uclean, alpha, rd, verb)
function [u, snr_out] = CS_TV(R, f, uclean, alpha, rd, verb)

if nargin<6
    verb = 1;
end;

%%%%% parameters
nAL = 10000;
nInner = 1;

%%%%%% normalize parameters to size of image and CS samples
[Ny,Nx] = size(R);
N = Nx*Ny;
alpha = alpha*sqrt(N)/nnz(R);
gamma = alpha/10; 
fftscale = sqrt(N);

%%%%%%%%%%% auxiliary variable for FFT-based minimization
[X,Y] = meshgrid(1:Nx,1:Ny);
auxFFT = cos(2*pi/Ny*Y) + cos(2*pi/Nx*X) - 2;


%%%%%%%%%%%%% initialize split variables
u = zeros(Ny,Nx);
upast = u;

dx = Dx(u);
dy = Dy(u);
ldx = zeros(Ny,Nx);
ldy = zeros(Ny,Nx);

K = gamma + alpha*conj(R).*R - 2*rd*auxFFT;
aux_rhs0 = alpha*ifft2(R.*(f))*fftscale;

%%%%%%%%%%%%% iterate AL algorithm
for iAL=1:nAL
    aux_rhs = aux_rhs0 + Dxt( ldx + rd*dx ) + Dyt( ldy + rd*dy );
    for inner = 1:nInner;
       
        %update u
        rhs =  aux_rhs + gamma*u;
        u  = real(ifft2( fft2(rhs)./K ));

        %update d
        [dx,dy] = shrink2( Dx(u) - ldx/rd, Dy(u) - ldy/rd, 1/rd);
    end;
    
    %update Lagrange multipliers
    ldx = ldx + rd*( dx - Dx(u) );
    ldy = ldy + rd*( dy - Dy(u) );
    
    
    % compute energy and AL residuals
    ux = Dx(u); uy = Dy(u); Rfftu = R.*fft2(u)/fftscale;
    Energy(iAL) = sum(sum( sqrt( ux.*conj(ux) + uy.*conj(uy) ) ))/N  + 0.5/N*alpha*norm( Rfftu - f , 'fro')^2;
    res(iAL) = 1/N*norm( dx-ux, 'fro' ) + 1/N*norm( dy-uy, 'fro' ); 
    var_noise(iAL) = 1/N*norm( Rfftu - f ,'fro');
    SNR(iAL) = snr(real(u),uclean);          
    
    % stopping condition
    stop_cond = [ ( res(iAL) > 1e-3 ) norm(u-upast,'fro')/norm(u,'fro') norm(u-upast,'fro')/N ];
    upast = u;
    if( max(stop_cond)<1e-6 )
        disp(['stop at iteration ' int2str(iAL)])
        break;
    end;      
    
end;
snr_out = SNR(iAL)

if verb
    cpt_fig = 101;
    figure(cpt_fig); clf;
    subplot(221); imshow(abs(u)); title(['TV SNR = ',num2str(SNR(iAL))]);
    subplot(222); semilogy(Energy+1e-9); title('Energy to minimize');
    subplot(223); semilogy(res+1e-9); title('residuals'); legend('grad u = d');
    subplot(224); plot(SNR); title('SNR');    
    pause(0.01);

    nerror = sqrt(conj(u-uclean).*(u-uclean));
    cpt_fig = 102;
    figure(cpt_fig); clf;
    subplot(131);imshow(real(uclean));title(['original image'])
    subplot(132);imshow(real(u));title(['reconstrection SNR=' num2str(SNR(iAL))]);
    subplot(133);imshow(nerror, [min(nerror(:)) max(nerror(:))]); title(['|| uclean - u ||']);
    pause(0.01);
end;


%%%%%%%%%%%% estimate differential functions with finite differences

function d = Dx(u)
[rows,cols] = size(u); 
d = zeros(rows,cols);
d(:,2:cols) = u(:,2:cols)-u(:,1:cols-1);
d(:,1) = u(:,1)-u(:,cols);
return

function d = Dxt(u)
[rows,cols] = size(u); 
d = zeros(rows,cols);
d(:,1:cols-1) = u(:,1:cols-1)-u(:,2:cols);
d(:,cols) = u(:,cols)-u(:,1);
return

function d = Dy(u)
[rows,cols] = size(u); 
d = zeros(rows,cols);
d(2:rows,:) = u(2:rows,:)-u(1:rows-1,:);
ud(1,:) = u(1,:)-u(rows,:);
return

function d = Dyt(u)
[rows,cols] = size(u); 
d = zeros(rows,cols);
d(1:rows-1,:) = u(1:rows-1,:)-u(2:rows,:);
d(rows,:) = u(rows,:)-u(1,:);
return



%%%%%%%%%%%% auxiliary shrinkage functions

function xs = shrink(x,lambda)

s = sqrt(x.*conj(x));
ss = s-lambda;
ss = ss.*(ss>0);

s = s+(s<lambda);
ss = ss./s;

xs = ss.*x;

return;


function [xs,ys] = shrink2(x,y,lambda)

s = sqrt(x.*conj(x)+y.*conj(y));
ss = s-lambda;
ss = ss.*(ss>0);

s = s+(s<lambda);
ss = ss./s;

xs = ss.*x;
ys = ss.*y;

return;

Contact us