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);