基于Shepp-Logan添加环形伪影
clc
clear all
L = 256;%图像大小
P = phantom(L);
theta = 0:359;
R = radon(P,theta);%获得正弦图
for ii = 120:120
for jj = 1:360
R(ii,jj) = 0.0;
end
end
%figure,imshow(R,[])
I = iradon(R, theta, 'linear', 'Hamming', L);%%获得重建图像
JG = 2*pi/(4*L);%%角度采样频率
width = 2*pi/(JG) + 1;
[m,n] = size(I);
radius = ceil(sqrt(m^2 + n^2)/2);
POL = zeros(radius*2, width);
for gama = 0:JG:(2*pi)
for rad = 0.25:0.25:radius %%半径采样频率
i = round(gama/JG+1);%%round四舍五入取整
j = (rad-0.25)*4 + 1;
x = rad*cos(gama);%%x坐标变换
y = rad*sin(gama);%%y坐标变换
ceil_x = ceil(x);%ceil向下取整
ceil_y = ceil(y);
floor_x = floor(x);
floor_y = floor(y);
chae3 = abs(x - floor_x);
chae4 = abs(y - floor_y);
x1 = floor_x + L/2;
x2 = ceil_x + L/2;
y1 = floor_y + L/2;
y2 = ceil_y + L/2;
if x1>L
x1 = L;
end
if x1<1
x1 = 1;
end
if x2>L
x2 = L;
end
if x2<1
x2 = 1;
end
if y1>L
y1 = L;
end
if y1<1
y1 = 1;
end
if y2>L
y2 = L;
end
if y2<1
y2 = 1;
end
POL(j,i) = (1-chae3)*(1-chae4)*I(x1,y1) + chae3*(1-chae4)*I(x2,y1) + (1-chae3)*chae4*I(x1,y2) + chae3*chae4*I(x2,y2);%%双线性插值
end
end
figure, hold on
subplot(1,3,1),imshow(I)%直接重建图像
subplot(1,3,2),imshow(POL)%极坐标变化后图像
subplot(1,3,3),imshow(R)%正弦图