From: Subject: Date: Tue, 11 Dec 2007 15:48:30 -0800 MIME-Version: 1.0 Content-Type: text/html; charset="Windows-1252" Content-Transfer-Encoding: quoted-printable Content-Location: http://www.stanford.edu/class/ee263/matlab/tmeasure.m X-MimeOLE: Produced By Microsoft MimeOLE V6.00.2900.3198

% 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=3D30;   % image size is n_pixels by n_pixels
Nd=3D35;        % number of parallel lines for each angle
Ntheta=3D35;    % number of angles (equally spaced from 0 to pi)
N=3DNd*Ntheta;  % total number of lines (i.e. number of measurements

sigma=3D0.7;   % noise level (standard deviation for normal dist.)

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

x=3DX(:);  % 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=3Dzeros(N,1);            % will store the N measurements
lines_d=3Dzeros(N,1);      % will store the position of each line
lines_theta=3Dzeros(N,1);  % will store the angle of each line
i=3D1;
for itheta=3D1:Ntheta
 for id=3D1:Nd
  lines_d(i)=3D0.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=3Dpi/4)
  lines_theta(i)=3Dpi*itheta/Ntheta;
           % equally spaced angles from 0 to pi
  L=3Dline_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=3DL(:);
           % make matrix L into a vector, as for X
  y(i)=3Dl'*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=3Di+1;
           % for next measurement line
 end
end