Date: 2025.4.20
代码获取方式
需要相关代码可以在关注博主和订阅本专栏《实用毕业设计》后加文章最后的QQ名片咨询博主。【强烈推荐】
前言
本文主要讲述了基于层析成像技术的CT三维重建方法,包括空域法(滤波反投影,FBP)和频域法(傅里叶切片定理)。
1、方法介绍
-
频域法(傅里叶切片定理、傅里叶衍射投影定理):频域法算法简单,成像速度快,但由于在半圆弧内非等间隔的采样,因此要在频域内插值,一定程度影响了重构图像的质量。但是如果改进迭代算法可以提高图像重建的质量,例如应用非均匀快速傅里叶变换。频域法重建的一般方法就是网格算法,将非均匀数据插值到直角坐标网格,如多项式插值法,然后,用傅立叶逆变换的快速算法来重建图像。
-
空域法(FBP,FBPP):空域法相比频域法重建的图像精度较高,由于采用了一个依赖于几何深度的滤波器,所以能够校正正反方向传播时的衍射效应。但它的缺点也显而易见,它需要有大量的投影,计算量大,经常会由于投影数据不足而产生混叠误差。投影重建的过程是,先把投影由线阵探测器上获得的投影数据进行一次一维傅立叶变换,再与滤波器函数进行卷积运算,得到各个方向卷积滤波后的投影数据;然后把它们沿各个方向进行反投影,即按其原路径平均分配到每一矩阵单元上,进行重叠后得到每一矩阵单元的CT值;再经过适当处理后得到被扫描物体的断层图像。
-
层析成像技术是在物体外部发射物理信号并使之穿过待重建物体,在接收端接收到带有物体内部信息的物理信号。利用计算机图像重建技术,重现物体内部的几何结构的一维或者三维图像。这种方法最大的优点是不会破坏物体内部的几何形态和物理参数,并且得到的图形清晰直观。层析成像技术应用广泛,例如:生物医学工程、无损检测、地球物理、模式识别等领域。(1)CT技术的应用,CT技术就是利用X射线扫描人体,并利用计算机重建出人体内部的三维模型,从而判断人体内部是否有病变或者肿瘤等的产生;(2)地震层析成像技术的应用,利用地震波在不同方向投影的波场信息,对地下介质精细结构进行成像,以其分辨率高、解析成果直观等特点广泛应用于建筑、公路、铁路、环境等方面工程的地质勘测中;(3)无损检测和质量检测应用,例如军事工业中,层析成像用于对炮弹、火炮的质量检查;在石油开发中被用于岩心分析和油管损伤检测。
2、Matlab部分实现
本文只提供了傅里叶切片定理的证明实现过程和Radon变换的使用方法,滤波反投影重建的实现可以通过左侧公告栏中的QQ图标或者直接QQ(2963033731)联系我。
I=phantom(256);
figure, imshow(I, []),title('256*256 Shepp-logan体模');
% Radon变换的使用
[R, xt] = radon(I, [0 45 90]); % 在0、45、90度方向做radon变换
subplot(2, 2, 2);
plot(xt, R(:, 1));
title('水平方向的radon变换曲线');
subplot(2, 2, 3);
plot(xt, R(:, 2));
title('45度方向的radon变换曲线');
subplot(2, 2, 4);
plot(xt, R(:, 3));
title('垂直方向的radon变换曲线');
P=I;
theta1=0:10:170;[R1,xp]=radon(P,theta1); %存在18个角度投影
theta2=0:5:175;[R2,xp]=radon(P,theta2); %存在36个角度投影
theta3=0:2:178;[R3,xp]=radon(P,theta3); %存在90个角度投影
figure,imagesc(theta3,xp,R3);colormap(hot);colorbar;
xlabel('\theta');ylabel('x\prime'); % 定义坐标轴
I1=iradon(R1,10);
I2=iradon(R2,5);
I3=iradon(R3,2);
figure,imshow(I1),title('I1');
figure,imshow(I2),title('I2');
figure,imshow(I3),title('I3');
%1.傅立叶切片定理
%获取 theta=0时候的投影数据
%projection=sum(I,1);
[projection,xt] = radon(I,0);
%figure,plot(1:256,projection),title('投影数据');
figure,plot(xt,projection),title('投影数据');
%对投影进行一维傅立叶变换
f=fft(projection);%对投影进行一维傅立叶变换
fab=f.* conj(f);%由于得到的是复数,我们通过共轭得到模
f3=fftshift(fab);%对所得一维傅立叶变换进行频谱转移
%figure,plot(1:256,f3);title('垂直线束投影的一维傅立叶图像');%显示图像
figure,plot(xt,f3);title('垂直线束投影的一维傅立叶图像');%显示图像
% 直接对原始图像进行二维傅立叶变换
f2D=fft2(I);%对原始图像进行二维傅立叶变换。
fab2D=f2D.* conj(f2D);%取结果的模
f3=fftshift(fab2D); %对所得二维傅立叶变换进行频谱转移
figure,plot(f3);title('当theta=0对原始图像进行二维傅立叶变换') %显示图像
P = phantom(256);
C = radon(P,0:179);
I = iradon(C, 0:179);
figure, imshow(C, []),title('重建图像');
figure, imshow(I, []),title('重建投影图');
%% 基于反投影及滤波反投影的CT图像重建
%2.滤波反投影重建
I=phantom(256);
figure,imshow(I);title('256*256 Shepp-logan体模');
IMG=double(I);
THETA=0:179;
3、实验效果
THE END!