% 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