function [ K, XX, YY ] = ConstructObsMatrixTransTomo(Nx,Ny,theta)
% ConstructObsMatrixTransTomo.m Constructs obsevation matrix for transmission
% tomography using Matlab function radon.m
%
%   Output:
%   K = observation matrix
%   XX = x-cooridinates of the image pixels
%   YY = y-coordinates of the image pixels
%
%   Input:
%   Nx = number of pixels in the image in x-direction
%   Ny = number of pixels in the image in x-direction
%   theta = vector consisting of projection angles
%
% Jyvaskylä Summer School 2014
% Aku Seppanen
% University of Eastern Finland, Kuopio
% Department of Applied Physics

[XX,YY] = meshgrid(1:Nx,1:Ny);
XX = XX(:);
YY = YY(:);

I = zeros(Ny,Nx);
[R,xp] = radon(I,theta);
RR = R(:);
ng_theta = size(R,1);
ng = length(RR);
K = zeros(ng,Ny*Nx); 

for iy = 1:Ny
    for ix = 1:Nx
        ik = (iy-1)*Nx + ix;
        I = zeros(Ny,Nx);
        I(ix,iy) = 1;
        [R,xp] = radon(I,theta);
        RR = R(:);
        K(:,ik) = (1/Nx)*RR;
    end
end

end


