% Reconstruct phantom using TV regularization
%
% Samuli Siltanen Feb 2017

% Regularization parameter
alpha = 100;

% Choose the size of the unknown. The image has size (2^j)x(2^j);
% let's denote n = (2^j)*(2^j).
j = 6;
M = 2^j;
n = M*M;

% Choose the angles 
Nang = 33;

% Load precomputed results of the routine tomo01_RadonMatrix_comp.m
loadcommand = ['load data/RadonMatrix_', num2str(j), '_', num2str(Nang), ' A angles'];
eval(loadcommand)
% Load precomputed results of the routine tomo03_NoCrimeData_comp.m
loadcommand = ['load data/tomo03_NoCrime_', num2str(M),'_',num2str(Nang), ' mncn f'];
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 (M^2)x(M^2). This implements difference
% between consecutive values assuming periodic boundary conditions.
LH = -eye(M^2,M^2)+diag(ones(M*(M-1),1),M)+diag(ones(M,1),-M^2+M);
LV = zeros(M^2,M^2);
L = eye(M);
L = L-[L(:,end),L(:,1:end-1)];
for jjj = 1:M
    LV((jjj-1)*M+[1:M],(jjj-1)*M+[1:M]) = -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/tomo07 M recn Nang mncn angles

tomo07_TV_plot

