% FUNCTION PROGRAM A = mktomomat(N,M,K)
%----------------------------------------------------------------------------
% The program computes the sparse matrix relating the density distribution
% matrix to the Radon transform.
% 
% INPUT  :  N  -  number of pixels in vertical and horizontal directions
%           M  -  number of illumination directions, 0<= phi < pi
%           K  -  number of rays in each illumination direction
% OUTPUT :  A  -  sparse matrix of size (M*K)x(N*N). The numbering of the
%                 pixels is (x(1),y(1)), (x(2),y(1)), ... , (x(N),y(1)),
%                 (x(1),y(2)), ... , (x(N),y(N)). The numbering of the
%                 tomography lines is (s(1),phi(1)), (s(2),phi(1)), ...
%                 (s(K),phi(1)), (s(1),phi(2)), ... (s(K),phi(M)). 
%----------------------------------------------------------------------------
% CALLS TO : None
% Erkki Somersalo
%----------------------------------------------------------------------------

function A = mktomomat(N,M,K)

x = -1 + 2/N*(0:N); 
y = x;
xd = kron(ones(1,N),x(1:N));
yd = kron(y(1:N),ones(1,N));
xu = kron(ones(1,N),x(2:N+1));
yu = kron(y(2:N+1),ones(1,N));
phi = pi/M*(0:M-1); 
s = -1 + 2/(K-1)*(0:K-1);

I = []; J = []; val = [];
for m = 1:M
  m
  cc = cos(phi(m)); ss = sin(phi(m));
  for k = 1:K
    if ss < eps  % implying cc = 1
      aux = (xd<s(k)).*(xu>=s(k));
      I = [I,((m-1)*K + k)*ones(1,sum(aux))];
      J = [J,find(aux)]; 
      val = [val,(yu(aux) - yd(aux))];
    else  % ss positive
      if abs(cc) < eps % implying ss = 1
        aux = (yd<s(k)).*(yu>=s(k));
        I = [I,((m-1)*K + k)*ones(1,sum(aux))];
        J = [J,find(aux)];
        val = [val,(xu(aux) - xd(aux))];
      elseif cc > 0
        lb = max((1/ss*(cc*s(k) - xu)),1/cc*(yd - ss*s(k)));
        ub = min((1/ss*(cc*s(k) - xd)),1/cc*(yu - ss*s(k)));
        aux = (lb<ub);
        I = [I,((m-1)*K + k)*ones(1,sum(aux))];
        J = [J,find(aux)]; 
        val = [val,(ub(aux) - lb(aux))];
      elseif cc < 0
        lb = max((1/ss*(cc*s(k) - xu)),1/cc*(yu - ss*s(k)));
        ub = min((1/ss*(cc*s(k) - xd)),1/cc*(yd - ss*s(k)));
        aux = (lb<ub);
        I = [I,((m-1)*K + k)*ones(1,sum(aux))];
        J = [J,find(aux)]; 
        val = [val,(ub(aux) - lb(aux))];
      end
    end
  end
end

A = sparse(I,J,val,M*K,N*N);
