% file: tmeasure.m
%
% image reconstruction from line measurements (tomography via least-squares)
%
% this is the matlab script that was used to produce the data in tomodata.m

n_pixels=30;   % image size is n_pixels by n_pixels
Nd=35;        % number of parallel lines for each angle
Ntheta=35;    % number of angles (equally spaced from 0 to pi)
N=Nd*Ntheta;  % total number of lines (i.e. number of measurements

sigma=0.7;   % noise level (standard deviation for normal dist.)

X=  % ...the mystery image goes here... (we hid this part!)

x=X(:);  % write the matrix X (the original image) as one big column vector
         % (first column of X goes first, then 2nd column, etc.)

figure(1)      % display the original image
colormap gray
imagesc(X)
axis image

y=zeros(N,1);            % will store the N measurements
lines_d=zeros(N,1);      % will store the position of each line
lines_theta=zeros(N,1);  % will store the angle of each line
i=1;
for itheta=1:Ntheta
 for id=1:Nd
  lines_d(i)=0.7*n_pixels*(id-Nd/2-0.5)/(Nd/2);
           % equally spaced parallel lines, distance from first to
           % last is about 1.4*n_pixels (to ensure coverage of whole
           % image when theta=pi/4)
  lines_theta(i)=pi*itheta/Ntheta;
           % equally spaced angles from 0 to pi
  L=line_pixel_length(lines_d(i),lines_theta(i),n_pixels);
           % L is a matrix of the same size as the image
           % with entries giving the length of the line over each pixel
  l=L(:);
           % make matrix L into a vector, as for X
  y(i)=l'*x+normrnd(0,sigma);
           % l'*x gives "line integral" of line over image,
           % that is, the intensity of each pixel is multiplied by the
           % length of line over that pixel, and then add for all pixels;
           % a random, Gaussian noise, with std sigma is added to the
           % measurement
  i=i+1;
           % for next measurement line
 end
end
