examp.m

这是在网上下的一个东东

%initialize the workspace
clear
rand('state',0);
randn('state',0);
%initialize the design matrix
G = zeros(94,256);

%number of model parameters in the problem
m = 256;

%parameter indices increase going down each column
%m(1) is in the NW corner, m(16) is in the SW corner
%m(241) in the NE corner, and m(256) is in the SE corner

%design matrix entries for the column scan
%ones where each ray hits a block, zeros otherwise
for i=1:16
  for j=(i-1)*16+1:i*16
    G(i,j) = 1.;
  end
end

%design matrix for the row scan
%ones where each ray hits a block, zeros otherwise
for i=1:16
  for j=i:16:240+i
    G(i+16,j) = 1.;
  end
end

%enlarge the g design matrix to accommodate the diagonal scans

%g matrix for the SW to NE diagonal scan, upper part
%because we're now going through blocks on the diagonal, we
%have entries of sqrt(2) where we intersect a block
for i=1:16
for j=0:i-1
G(i+32,i+j*15) = sqrt(2.);
end
end

%g matrix for the SW to NE diagonal scan, lower part
for i=1:15
for j=0:15-i
G(i+48,(i+1)*16+j*15) = sqrt(2.);
end
end

%g matrix for the NW to SE diagonal scan, lower part
for i=1:16
for j=0:i-1
G(i+63,17-i+17*j) = sqrt(2.);
end
end

%g matrix for the NW to SE diagonal scan, upper part
for i=1:15
for j=0:15-i
G(i+79,(i*16)+1+17*j) = sqrt(2.);
end
end

%
% Setup the true model.
%
mtruem=zeros(16,16);
mtruem(9,9)=1;
mtruem(9,10)=1;
mtruem(9,11)=1;
mtruem(10,9)=1;
mtruem(10,11)=1;
mtruem(11,9)=1;
mtruem(11,10)=1;
mtruem(11,11)=1;
mtruem(2,3)=1;
mtruem(2,4)=1;
mtruem(3,3)=1;
mtruem(3,4)=1;
mtruev=reshape(mtruem,256,1);

%
% Compute the data.
%
dtrue=G*mtruev;
%
% Add the noise.
%
rand('state',0);
d=dtrue+0.01*randn(size(dtrue));

%
% First, plot the true model.
%

figure(1);
bookfonts;
imagesc(mtruem,[-0.2 1.2]);
%;
colormap(gray);
H=colorbar;
set(H,'FontSize',18);
%print -deps mtomo.eps

%
% Next, compute the generalized inverse solution, using 87 singular values.
%
t=cputime;
[U,S,V]=svd(G);
UP=U(:,1:87);
VP=V(:,1:87);
SP=S(1:87,1:87);
mg=VP*inv(SP)*UP'*d;
disp('cputime for gen inv solution ');
cputime-t

%
% And plot the solution.
%
figure(2);
bookfonts;
imagesc(reshape(mg,16,16),[-0.2 1.2]);
%;
colormap(gray);
H=colorbar;
set(H,'FontSize',18);
%print -deps mg.eps

%
% Next,compute the Kaczmarz solution, 94 iterations.
%
t=cputime;
mkac=kac(G,d,0.0,200);
disp('Kac time');
cputime-t
%
% And plot the solution.
%
figure(3);
bookfonts;
imagesc(reshape(mkac,16,16),[-0.2 1.2]);
%;
colormap(gray);
H=colorbar;
set(H,'FontSize',18);
%print -deps mkac.eps

%
% Next,compute the ART solution, 200 iterations.
%
t=cputime;
mart=art(G,d,0.0,200);
disp('ART time');
cputime-t
%
% And plot the solution.
%
figure(4);
bookfonts;
imagesc(reshape(mart,16,16),[-0.2 1.2]);
%;
colormap(gray);
H=colorbar;
set(H,'FontSize',18);
%print -deps mart.eps

%
% Next,compute the SIRT solution, 200 iterations.
%
t=cputime;
msirt=sirt(G,d,0.0,200);
disp('SIRT time');
cputime-t
%
% And plot the solution.
%
figure(5);
bookfonts;
imagesc(reshape(msirt,16,16),[-0.2 1.2]);
colormap(gray);
H=colorbar;
set(H,'FontSize',18);
%print -deps msirt.eps

disp('kac relative error')
norm(mkac-mtruev)/norm(mtruev)
disp('ART relative error')
norm(mart-mtruev)/norm(mtruev)
disp('SIRT relative error')
norm(msirt-mtruev)/norm(mtruev)
disp('SVD relative error')
norm(mg-mtruev)/norm(mtruev)