1.算法描述
从投影重建物体的截面图像是图像处理中非常重要的技术此技术在物体的无损伤性检测其内部缺陷的应用中能起很大作用从投影重建图像的技术早在20世纪中期就已经制成常规医疗诊断设备的商品1917年奥地利数学家J.Radon发表的论文证明了二维物体或三维物体可以从许多投影来重建其内部数据英国的EMI公司中央研究所1972年由G.N.Hounsfield设计出X射线断层分析仪此后于1974年5月在加拿大蒙特利尔第一次国际会议上正式命名为Computerized Tomography简称XCT或CT当时这台仪器可以获得人体的各个部位鲜明清晰的断层图像另外美国塔夫茨大学A.M.Cormack在1963年就提出了精确的非迭代的级数展开法用于图像重建的算法中并指出在诊断医学中图像重建的可能性以上两人由于对CT研制工作的开创性贡献获得了1979年诺贝尔医学奖。
计算机层析成像技术（Computed Tomography，简称CT）是利用具有一定能量的射线源（X射线，γ射线）对物体进行断层扫描，并根据物体外部的探测器获得的物理量（指物质对射线的衰减系数）生成的一维投影数据，通过特定的重建算法，得到所扫描断层的二维图像。CT断层图像具有无影像重叠、空间和密度分辨率高、可直接进行数字化处理等优点，通过近十几年的发展已成为非接触无损检测的主流技术，是关键部件检测、机械仿型设计、安全检查等方面强有力的手段，并广泛应用于航空、航天、机械、汽车、船舶、公安等领域。此外，在医学领域，医用CT机已成为疾病诊断的重要辅助工具，对病灶的位置及病变程度的良好再现已使其成为医学影像领域里的主要医疗器械。CT机有三个组成部分。一是数据采集系统，二是图像重建（主要是图像重建算法），三是图像后处理及显示系统。CT图像重建算法主要分为两类，即变换法和级数展开法。变换法包括滤波反投影算法、傅立叶变换法以及ρ滤波法（rho-filtered layergrams），小波变换法；级数展开法主要是代数重建算法。滤波反投影是比较常用的变换算法，它具有速度快，空间和密度分辨率好的优点。
滤波反投影法是目前CT图像重建领域使用得最为广泛的一种算法，其有效性和准确性也在临床上得到了验证。
通常在计算机中，采用的离散反投影算法，因此离散情形下的卷积反投影算法可以如下描述：
其中O(x，y)代表原始图像，R(x，y)代表重建图像。当然孤立地去看相对误差并不是一个好的衡量标准，它必须同人们的视觉评估结合起来进行考虑。
shepp-Logan模型:
用S-L头模型进行计算机仿真研究的主要好处之一是可以获得该模型投影数据的解析表达式。一个中心位置在原点且未经旋转的椭圆，其长轴与X轴重合，短轴与Y轴重合。假设椭圆内的密度为r、椭圆外密度为零，则该椭圆图像可用以下方程表示：
式中A，B分别为椭圆长轴与短轴的长度。
在研究从投影重建图像的算法时，为了比较客观的评价各种重建算法的有效性，人们常选用公认的Sheep-Logan 头模型作为研究对象。该模型由10个位置、大小、方向、密度各异的椭圆组成，象征一个脑断层图像。
该模型的主要参数如下表所示：
多尺度全局重建算法计算流程:
2.仿真效果预览
matlab2022a仿真结果如下：
3.MATLAB核心程序
function filtPR=lvbo(PR);
n = size(PR,1);
sideSize = n;
a = 1;
[Length, Count] = size(PR);
w = [-pi:(2*pi)/Length:pi-(2*pi)/Length];
rn1 = abs(2/a*sin(a.*w./2));
rn2 = sin(a.*w./2);
rd = (a*w)./2;
r = rn1*(rn2/rd)^2;
f = fftshift(r);
for i = 1:Count
IMG = fft(PR(:,i));
fimg = IMG.*f';
g(:,i) = ifft(fimg);%利用FFT,IFFT进行滤波
end
filtPR= real(g);
............................................................
I=phantom(256); %生产头部模型图
figure(1);
imshow(I); %显示图像
IMG=double(I); %双精度显示
[cod_a,cod_h,cod_v,cod_d,map]=wtest(IMG);%一层小波变换
figure(2);
subplot(221)%显示分解系数1
imshow(cod_a,map);
subplot(222)%显示分解系数2
imshow(cod_h,map);
subplot(223)%显示分解系数3
imshow(cod_v,map);
subplot(224)%显示分解系数4
imshow(cod_d,map);
Z1=idwt2(cod_a,cod_h,cod_v,cod_d,'db1');%重建
figure(3);
imshow(Z1,map);
for i=1:256
for j=1:256
error(i,j)=(I(i,j)-Z1(i,j))^2/I(i,j)^2;%误差求解
end
end
for i=1:256
for j=1:256
if I(i,j)==0
error(i,j)=255;
end
end
end
for i=1:256
for j=1:256
if error(i,j)>255;
error(i,j)=255;
end
end
end
