function [k_NCP,Nk] = NCPforCT(X,A,b,angles)
%NCPforCT The NCP criterion for 2D CT reconstruction.
%
% Implements the NCP criterion for 2D CT reconstruction problems where
% the solution vector represents a 2D image and the right-hand side
% represents a sinogram. The output k_NCP is the index for which the
% residual sinogram is closest to white noise.
%
% Input: X - Array in which each solumn is an iteration vector.
% A - The system matrix.
% b - The right-hand side.
% angles - A vector with the projection angles.
%
% Output: k_NCP - The index for which the residual is closest to white noise.
% Nk - The NCP function that measures the distance to white noise.
% Per Christian Hansen, DTU Compute, October 1, 2021.
% P. C. Hansen and J. S. Jorgensen, "AIR Tools II: algebraic iterative
% reconstruction methods, improved implementation," Numerical Algorithms,
% 79 (2018), pp. 107--137.
% Initialization.
nt = length(angles);
np = size(A,1)/nt;
q = floor(np/2);
c_white = (1:q)'./q;
C = zeros(q,nt);
Nk = zeros(size(X,2),1);
% Loop over iteration vectors in X.
for i=1:size(X,2)
rk = b - A*X(:,i); % Residual vector.
R = reshape(rk,np,nt); % Reshape as residual sinogram.
for j=1:nt % Loop over projection angles.
RKH = fft(R(:,j));
pk = abs(RKH(1:q+1)).^2;
c = cumsum(pk(2:end))/sum(pk(2:end));
C(:,j) = c;
end
Nk(i,1) = norm(mean(C,2)-c_white);
end
[~,k_NCP] = min(Nk);