% deblur3_sb_template.m
% Split-Bregman for sparsity-based edge-preserving image deblurring

if 1 % read data
	load deblur3_data % loads image variable 'y' and PSF variable 'psf'
	[nx ny] = size(y);
	y = double(y);
	nd = nx * ny; % # of data points
	A = Gblur(true(nx,ny), 'psf', psf); % system matrix
	C = Cfd2(nx,ny); % 2D finite differences
end

im plc 2 3
clim = [-10 130]; % display all images on same gray scale
im(1, y, clim);
title('y');
drawnow
y = y(:);

%% general parameters
% beta = 1;
niters = 100; % # of outer (Split-Bregman) iterations
niters_PCG = 5; % # of inner (preconditioned conjugate gradients) iterations

% original function evaluation
% Psi(x) = ||Ax-y||^2/2+beta*||Cx||_1
Psi = @(x) 0.5*norm(col(A * x - y),2)^2 + beta*norm(col(C * x),1);

%% perform Split-Bregman for various mu's

mulist = 1; % for debugging purposes, start with just one "mu"
% mulist = 2.^(-10:5:10); % TODO: adjust mulist to find favorable mu

nmu = length(mulist);
fobjs = zeros(niters+1,1); % vector of objective function values
rs = zeros(niters,1); % vector of residual norms ||w - Cx||_2
fobjs_final = zeros(nmu,1); % final value of objective function for each mu
rs_final = zeros(nmu,1); % final value of residual norm for each mu

for imu = 1:nmu
    mu = mulist(imu);

    % TODO: compute preconditioner for the PCG algorithm
    P = 1;
    
    % initialization
    x = y; % image initialized to the data (note: x is vectorized)
    w = C * x; % initialize the split variable based on initial x
    b = 0; % initially, set scaled dual variable to zero
    
    % Split-Bregman iterations
    fobjs(1) = Psi(x); % initial objective function value

%     pause(1); % in case you want to see more slowly
    
    for ii = 1:niters % outer iterations of Split-Bregman
        ticker(mfilename, ii, niters);

        % TODO: using several iterations of the (preconditioned) conjugate
        % gradient algorithm, approximately solve the subproblem for x
        
%         x = qpwls_pcg3(x,...
%             ,... % TODO: determine remaining arguments for this subproblem 
%             'niter',niters_PCG,'precon',P,'isave','last');
        
        % TODO: solve subproblem for w; this update should have a closed form
%         w = ...;
        
        % TODO: perform scaled dual variable update
%         b = ...;
        
        rs(ii) = norm(C * x - w,2); % residual for this iteration
        fobjs(ii+1) = Psi(x); % objective function value for this iteration
    end
    x = reshape(x,[nx,ny]);
    rs_final(imu) = rs(end);
    fobjs_final(imu) = fobjs(end);

    % plot reconstructed image
    im(4, x, clim);
    title('x*','Interpreter','tex');
%    title('x^*'); % octave
    xlabelf('log_2(\\mu) = %4.1f', log2(mu)); % matlab?
%    xlabelf('log_2(\mu) = %4.1f', log2(mu)); % octave

    im subplot 2
    semilogy(0:niters, fobjs-f_inf), hold all;
    title('\Psi(x) - \Psi(x^{\infty})'); 
    xlabel('n');
        
    im subplot 3
    plot(1:niters, 20*log10(rs)), hold all;
    title('||Cx-w||_2 (dB)'); 
    xlabel('n');

    im subplot 5
    plot(log2(mulist(1:imu)), fobjs_final(1:imu), '-o')
    title('\Psi(x*)'); 
    xlabel('log_2(\mu)');
        
    im subplot 6
    plot(log2(mulist(1:imu)), 20*log10(rs_final(1:imu)), '-o')
    title('||Cx*-w*||_2 (dB)'); 
    xlabel('log_2(\mu)');

    drawnow;
end
