mlem.m

关于最大似然重建方法的实现

clear
clc
time=cputime;
load test_phantom_zhou;
% I=imread('4.bmp');
imag_size=127;
% p=shepplogan(imag_size,imag_size,1);
I=double(I);
p=imresize(I,[imag_size,imag_size]);
% p=p./max(max(p));
% p=shepplogan(imag_size,imag_size,1);
% % %p=phantom('Modified Shepp-Logan',imag_size);%产生shepp-logan模板
% % %im(p)
% % %imshow(p)
% % %b=radon(p);%产生模板的投影数据
nb=imag_size;%确定投影数据的维数
na=180;
% % %[nb,na]=size(b);%nb为detector得数目，na为所取投影角的数目
% p1=reshape(p,imag_size*imag_size,1);
% load 'ct-sino-s509-P'
a=getG(imag_size,imag_size,nb,na);%产生概率矩阵a
a=a./sum(a(:,1),1);
% load 'phantom-sino-noised-510e6'
% b1=sinoo;
b=a*p(:);%产生投影数据
[b,c1]=psnoise(b, 5*10e6);
b=reshape(b,nb,na);
% b1=b;
b1=zeros(nb,na);
index1=35;%%所知角度的范围：25～155
index2=180-index1;
b1(:,index1:index2)=b(:,index1:index2);

% [b1,c1]=psnoise(b, 10e5);%在投影数据中加入噪声
% fid=fopen('sinogram_2');
% b=fread(fid,'short');
% b1=max(b,0.0);
% b1=b;
% load 'noised-data-128.mat'
% load 'noisedFree-data-128.mat'
% load 'real-data.mat'

itera_num=40;%确定迭代次数
n=imag_size*imag_size;
m=size(a,1);
w=zeros(m,1);
x=ones(n,1);
c=(sum(a))';
result_ml=zeros(n,itera_num);
for j=1:itera_num
    j
    w=a*x;
    mask=w>0.0;%对w中的零点不进行处理
    w(mask)=b1(mask)./w(mask);
    s=a'*w;
    x=x.*s;
    x=x./c;
    x=max(x,0.0);
    result_ml(:,j)=x;
end
% figure(2)
result=reshape(result_ml(:,itera_num),imag_size,imag_size);
time2=cputime-time;
%result=result';
iptsetpref('ImshowBorder','tight')
imshow(result,[])
% im(result1)
%  si_ml=reshape(p,n,1);
%  error_ml=rsse(result_ml,si_ml);
p=p./max(max(p));
re=result./max(max(result));
error_ml=mse(re,p);
%   figure(2)
% plot(error_ml)