% Reconstruct Shepp-Logan phantom using TV regularization
%
% Samuli Siltanen Feb 2015

% Regularization parameter
alpha = 1;

% Choose the size of the image (attenuation function)
N = 64;
n = N^2;

% Choose the angles between 0 and 180 degrees
T  = 15;

% Load precomputed data
loadcommand = ['load data/XRmatrix_',num2str(N),'_',num2str(T),' A th'];
eval(loadcommand)
loadcommand = ['load data/XRMC_NoCrime_', num2str(N),'_',num2str(T), ' mncn target'];
eval(loadcommand)

% Construct the quadratic optimization problem matrix
H = zeros(5*n);
H(1:n,1:n) = 2*A.'*A;

% Construct the vector h of the linear term
h = alpha*ones(5*n,1);
h(1:n) = -2*A.'*mncn(:);

% Construct prior matrix of size (N^2)x(N^2). This implements difference
% between consecutive values assuming periodic boundary conditions.
LH = -eye(N^2,N^2)+diag(ones(N*(N-1),1),N)+diag(ones(N,1),-N^2+N);
LV = zeros(N^2,N^2);
L = eye(N);
L = L-[L(:,end),L(:,1:end-1)];
for jjj = 1:N
    LV((jjj-1)*N+[1:N],(jjj-1)*N+[1:N]) = -L;
end

% Compute reconstruction using quadprog
% Construct input arguments for quadprog.m
Aeq         = [...
    [LH,-eye(n),eye(n),zeros(n),zeros(n)];...
    [LV,zeros(n),zeros(n),-eye(n),eye(n)]];
beq         = zeros(2*n,1);
lb          = zeros(5*n,1);
ub          = repmat(Inf,5*n,1);
AA          = -eye(5*n);
AA(1:n,1:n) = zeros(n,n);
iniguess    = zeros(5*n,1);
b           = [repmat(10,n,1);zeros(4*n,1)];
QPopt   = optimset('quadprog');
QPopt   = optimset(QPopt,'Algorithm','interior-point-convex','Display','iter');
[uvv,val,ef,output] = quadprog(H,h,AA,b,Aeq,beq,lb,ub,iniguess,QPopt);
recn = uvv(1:n);

% Save results to file
save data/XR09 N recn T mncn

XR09_TV_plot

