CT图像伪影MATLAB仿真

在去金属伪影过程中，为了研究金属伪影的去除方法，需要获取相应的图像。对此，本文采用了MATLAB对该过程进行仿真。

处理流程

1. 获取原始图像
2. 在原始图像中添加金属
3. 正投影，得到有金属轨迹的投影图
4. 获取金属投影图
5. 对金属轨迹的投影值进行变换
6. 反投影得到有金属伪影的图像

前三步的结果图如下所示：

投影图变换

获取投影图像，并对金属区域采用如下变换函数，即可得到目标图像。

最终的图像如下所示，可以看到金属周围的明暗条纹。

不进行变换的情况下，直接进行反投影，

完整实现代码如下：

% 金属伪影模拟
clc; 
clear; 
close all
%% 建立原始图像并加入金属球
nSize = 402; 
Img0=phantom('Modified Shepp-Logan',nSize);
figure;imshow(Img0),title('原始无金属图像');

metalMask = zeros(size(Img0));
vec = -200:201;
[xx,yy] = meshgrid(vec,-vec);
%加入金属球
x1 = 30;y1 = 50;r1 = 80;
Img0((xx-x1).^2+(yy-y1).^2<=r1) = 2;
metalMask((xx-x1).^2+(yy-y1).^2<=r1) = 2;
x1 = 60;y1 = 100;r1 = 40;
Img0((xx-x1).^2+(yy-y1).^2<=r1) = 2;
metalMask((xx-x1).^2+(yy-y1).^2<=r1) = 2;

R0 = radon(Img0);
srcImg = iradon(R0,0:179, nSize);

figure;imshow(Img0),title('原始有金属图像-无伪影');
figure;imshow(R0, []),title('原始有金属投影图-理想投影');
figure;imshow(srcImg, []),title('理想投影重建');

metalImg = metalMask;
notMetalImg = srcImg;
notMetalImg(metalImg>0) = 0;
%% 投影与反投影
%Radon变换的角度范围
radon_theta=0:179;
% 金属的Radon变换
R_metal=radon(metalImg,radon_theta);

figure,imshow(R_metal,[]),title('金属的Radon变换');
xlabel('\theta');

%非金属部分的Radon变换
R_notMetal=radon(notMetalImg,radon_theta);

figure,imshow(R_notMetal,[]),title('非金属部分的Radon变换');
xlabel('\theta');

%% 金属投影区域进行模拟变换
% prj_metalSim = R_notMetal;
prj_metalSim = R0;
[m, n] = size(prj_metalSim);
a = 10;
c = 0.5;
for i=1:m
    for j=1:n
        if R_metal(i,j) > 0
            x = prj_metalSim(i,j);
            h =x- (2*c*(x-a) +1 - sqrt(4*c*(x-a) + 1))/(2*c)*sign(x-a);
            prj_metalSim(i,j) = h*4;
        end
    end
end
metal_image=iradon(prj_metalSim, radon_theta, nSize);

figure;imshow(prj_metalSim,[]);title('模拟金属的投影图')
figure;imshow(metal_image,[]);title('模拟金属伪影的CT图')

x = 10:200;
h =x- (2*c.*(x-a) +1 - sqrt(4*c.*(x-a) + 1))/(2*c).*sign(x-a);
figure; plot(x, h), title('变换函数')

save('data.mat', 'metal_image', 'metalMask', 'prj_metalSim')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

如有错误，敬请指正！
如果觉得本文对你有帮助，欢迎一键三连。