Возвращение проекции с синограммы в MATLAB
Я работаю над реконструкцией рентгеновских томографических изображений и пытаюсь восстановить изображение в MATLAB, используя его синограмму, но оно не работает должным образом. Здесь вы можете найти код, который я использовал. Первая часть предназначена для расчета синограммы, которая работает хорошо. Вторая часть (начиная со строки, в которой написано% backprojection) предназначена для расчета обратной проекции (без какого-либо фильтра), и проблема в этой части. Восстановленное изображение не является правильным. Кто-нибудь может мне помочь с кодом или ввести другой код для расчета обратной проекции? Кстати, я не хочу использовать встроенные коды MATLAB, такие как радон или ирадон.
clc
clear all
I=phantom('Modified Shepp-Logan',600);
[N,M]=size(I);
% center of the image
m=round(M/2);
n=round(N/2);
% total number of rhos = number of pixels on the diagonal
rhomax=ceil(sqrt(M^2+N^2));
rc=round(rhomax/2);
%matrix to store the results
res=cast(zeros(rhomax+1,180),'double');
tic
for t=1:45 % below 45 degree use y as variable
costheta=cos(t*pi/180);
sintheta=sin(t*pi/180);
a=-costheta/sintheta; % y=ax+b
for r=1:rhomax
rho=r-rc;
b=rho/sintheta; % y=ax+b
ymax=min(round(-a*m+b),n-1);
ymin=max(round(a*m+b),-n);
for y=ymin:ymax
x=(y-b)/a;
xfloor=floor(x); % the integer part of x
xup=x-xfloor; % the decimal of x
xlow=1-xup;
x=xfloor;
x=max(x,-m);
x=min(x,m-2);
res(rhomax-r+1,180-t)=res(rhomax-r+1,180-t)+xlow*I(y+n+1,x+m+1);
res(rhomax-r+1,180-t)=res(rhomax-r+1,180-t)+xup*I(y+n+1,x+m+2);
end
end
end
for t=46:90
costheta=cos(t*pi/180);
sintheta=sin(t*pi/180);
a=-costheta/sintheta; % y=ax+b
for r=1:rhomax
rho=r-rc;
b=rho/sintheta; % y=ax+b
xmax=min(round((-n-b)/a),m-1);
xmin=max(round((n-b)/a),-m);
for x=xmin:xmax
y=a*x+b;
yfloor=floor(y);
yup=y-yfloor;
ylow=1-yup;
y=yfloor;
y=max(y,-n);
y=min(y,n-2);
res(rhomax-r+1,180-t)=res(rhomax-r+1,180-t)+ylow*I(y+n+1,x+m+1);
res(rhomax-r+1,180-t)=res(rhomax-r+1,180-t)+yup*I(y+n+2,x+m+1);
end
end
end
for t=91:135
costheta=cos(t*pi/180);
sintheta=sin(t*pi/180);
a=-costheta/sintheta;
for r=1:rhomax
rho=r-rc;
b=rho/sintheta;
xmax=min(round((n-b)/a),m-1);
xmin=max(round((-n-b)/a),-m);
for x=xmin:xmax
y=a*x+b;
yfloor=floor(y);
yup=y-yfloor;
ylow=1-yup;
y=yfloor;
y=max(y,-n);
y=min(y,n-2);
res(rhomax-r+1,180-t)=res(rhomax-r+1,180-t)+ylow*I(y+n+1,x+m+1);
res(rhomax-r+1,180-t)=res(rhomax-r+1,180-t)+yup*I(y+n+2,x+m+1);
end
end
end
for t=136:179
costheta=cos(t*pi/180);
sintheta=sin(t*pi/180);
a=-costheta/sintheta; % y=ax+b
for r=1:rhomax
rho=r-rc;
b=rho/sintheta;
ymax=min(round(a*m+b),n-1);
ymin=max(round(-a*m+b),-n);
for y=ymin:ymax
x=(y-b)/a;
xfloor=floor(x); % the integer part of x
xup=x-xfloor; % the decimal of x
xlow=1-xup;
x=xfloor;
x=max(x,-m);
x=min(x,m-2);
res(rhomax-r+1,180-t)=res(rhomax-r+1,180-t)+xlow*I(y+n+1,x+m+1);
res(rhomax-r+1,180-t)=res(rhomax-r+1,180-t)+xup*I(y+n+1,x+m+2);
end
end
end
for t=180
rhooffset=round((rhomax-M)/2);
for x=1:M
r=x+rhooffset;
r=rhomax-r+1;
for y=1:N
res(r,t)=res(r,t)+I(y,x);
end
end
end
toc
rhoaxis=(1:rhomax+1)-rc;
figure
imshow(res,[]);
%backprojection
res1=cast(zeros(600),'double');
[N1,M1]=size(res1);
m1=round(M1/2);
n1=round(N1/2);
rhomax1=ceil(sqrt(m1^2+n1^2));
rc1=round(rhomax1/2);
tic
for t=1:45 % below 45 degree use y as variable
costheta=cos(t*pi/180);
sintheta=sin(t*pi/180);
a=-costheta/sintheta; % y=ax+b
for r=1:rhomax1
rho=r-rc1;
b=rho/sintheta; % y=ax+b
ymax=min(round(-a*m1+b),n1-1);
ymin=max(round(a*m1+b),-n1);
for y=ymin:ymax
x=(y-b)/a;
xfloor=floor(x); % the integer part of x
xup=x-xfloor; % the decimal of x
xlow=1-xup;
x=xfloor;
x=max(x,-m1);
x=min(x,m1-2);
res1(y+n1+1,x+m1+1)=res1(y+n1+1,x+m1+1)+res(rhomax1-r+1,t);
end
end
end
for t=46:90
costheta=cos(t*pi/180);
sintheta=sin(t*pi/180);
a=-costheta/sintheta; % y=ax+b
for r=1:rhomax1
rho=r-rc1;
b=rho/sintheta; % y=ax+b
xmax=min(round((-n1-b)/a),m1-1);
xmin=max(round((n1-b)/a),-m1);
for x=xmin:xmax
y=a*x+b;
yfloor=floor(y);
yup=y-yfloor;
ylow=1-yup;
y=yfloor;
y=max(y,-n1);
y=min(y,n1-2);
res1(y+n1+1,x+m1+1)=res1(y+n1+1,x+m1+1)+res(rhomax1-r+1,t);
end
end
end
for t=91:135
costheta=cos(t*pi/180);
sintheta=sin(t*pi/180);
a=-costheta/sintheta;
for r=1:rhomax1
rho=r-rc1;
b=rho/sintheta;
xmax=min(round((n1-b)/a),m1-1);
xmin=max(round((-n1-b)/a),-m1);
for x=xmin:xmax
y=a*x+b;
yfloor=floor(y);
yup=y-yfloor;
ylow=1-yup;
y=yfloor;
y=max(y,-n1);
y=min(y,n1-2);
res1(y+n1+1,x+m1+1)=res1(y+n1+1,x+m1+1)+res(rhomax1-r+1,t);
end
end
end
for t=136:179
costheta=cos(t*pi/180);
sintheta=sin(t*pi/180);
a=-costheta/sintheta; % y=ax+b
for r=1:rhomax1
rho=r-rc1;
b=rho/sintheta;
ymax=min(round(a*m1+b),n1-1);
ymin=max(round(-a*m1+b),-n1);
for y=ymin:ymax
x=(y-b)/a;
xfloor=floor(x); % the integer part of x
xup=x-xfloor; % the decimal of x
xlow=1-xup;
x=xfloor;
x=max(x,-m1);
x=min(x,m1-2);
res1(y+n1+1,x+m1+1)=res1(y+n1+1,x+m1+1)+res(rhomax1-r+1,t);
end
end
end
toc
rhoaxis1=(1:rhomax1+1)-rc1;
res1=res1./90;
figure
imshow(res1,[]);