公司开发个网站,青州营销型网站建设,湛江市研发网站建设,做网站广告联盟赚钱转自#xff1a;CT图像重构方法详解——傅里叶逆变换法、直接反投影法、滤波反投影法_Absolute Zero-CSDN博客_反投影法
绪 在做CT图像处理的时候遇到很多问题#xff0c;对于滤波反变换有许多细节存在疑问#xff0c;经过多天查找资料和利用MATLAB程序一步步实…转自CT图像重构方法详解——傅里叶逆变换法、直接反投影法、滤波反投影法_Absolute Zero-CSDN博客_反投影法
绪 在做CT图像处理的时候遇到很多问题对于滤波反变换有许多细节存在疑问经过多天查找资料和利用MATLAB程序一步步实现后终于豁然开朗于是想要总结成文作为笔记方便今后查看。文中若有错误欢迎指出 目 录
1傅里叶逆变换法
2直接反投影法
3滤波反投影法
4MATLAB程序 CT图像的形成和重构在数学上的描述分别为拉冬变换和拉冬反变换(RANDON:氡。简单点说拉动变换是用于描述通过X射线扫描物体形成CT图像这个过程的一种数学表示而拉东反变换描述的则是对CT投影数据进行重构还原成物体图像的一种数学方法。本文章主要讲的是CT图像重构的方法也就是拉东反变换的具体实现方法当然可能也有小伙伴想了解一下什么是拉东变换这里贴上传送门拉冬变换。(Radon拉冬变换可以用来进行直线识别的算法。它是一种原始灰度图像到p,o二维矩阵的映射映射关系是灰度图像的像素值对参数p,o的直线的线积分或者叫投影可理解为垂直于这个直线方向的像素值累计求和)
其中拉冬反变换在数学上的表达式如下 由公式可看出拉冬逆变换在数学上实际上是一个二维傅里叶逆变换但由于二维傅里叶逆变换计算量大耗时长因此计算机对于拉冬逆变换的具体实现方法主要为反投影法反投影法又分为直接反投影法和滤波反投影法。由此可见CT重构方法主要为以下三种
1傅里叶逆变换法
2直接反投影法
3滤波反投影法
下面将具体对这三种图像重构方法进行解释。 1傅里叶逆变换法
该方法基于一个重要的定理中心切片定理。
该定理简单地理解就是通过角度为θ扫描得到的投影该投影的一维傅里叶变换与对整个图像二维傅里叶变换后二维频域中对应θ角度的一个切片信号是相同的下面两个图理解起来更直观。 图1 图2
根据该理论傅里叶逆变换法可以简单分成以下步骤
① 假设每旋转1°就扫描一次当对物体扫描了180°之后我们就能得到180个投影信号就是180根投影线→在临床上若使用平行扫描CT我们拿到手的数据就是这个在数学上就是对图像进行拉冬变换
② 对180个投影信号进行一维傅里叶变换
③ 对②得到的180个一维频域信号根据相对应的扫描角度在空间中旋转排列拼成一个二维频域空间如图3
④ 由于数据是离散的直接按照角度进行排列难以铺满整个二维空间因此还需要对空缺的地方进行插值一般三次样条插值效果最好但插值会带来一定误差。另外由于中心的信号密集周围的信号稀疏显然会损失一部分高频数据造成高频信号失真这就是采用傅里叶逆变换法重构图像时会使得图像边缘模糊的原因。
⑤ 对④中拼接而成的二维图像进行二维傅里叶逆变换就可重构原图 图3 傅里叶逆变换法
该方法的缺点有
a/高频信号有所失真
b/在插值时还涉及到极坐标和直角坐标的变换计算量大
c/需要用到二维傅里叶逆变换总体耗费时间长
由于计算机处理二维傅里叶逆变换的计算量太大因此很少直接使用该方法实现拉东逆变换。 2直接反投影法
反投影法的原理是将所测得的投影值按照其原投影路径平均地分配到经过的每一个点上把各个方向的投影值都这样反投影后在把每个角度的反投影图像进行累加从而推断出原图。为方便理解下面分别用两张图通过MATLAB编程来还原一下该过程程序在最后
原图如下图1简单图2复杂 ①首先看一下图像经过拉东变换即CT扫描后的数据图中每一列代表对应每一个角度的投影信号在这里只是拼起来成一张图了 ②对每个角度的信号进行反投影选了几个代表性角度 ② 把以上角度发反投影图像叠加起来 当越来越多的反投影图像加起来之后可以看到 到此直接反投影法结束。值得注意的是由于反投影图是离散叠加的显然在中心处信号集中边缘处信号稀疏因此在最后需要在空缺的地方进行插值才能得到最终的原图像。上面演示的过程中没有进行插值仅仅是把反投影图叠加因此能够看出每一张图都不是很光滑。
从以上过程可以很直观地看出直接反投影法出现伪迹的原因
①在“回抹“过程中就是把投影信号平均到每个二维空间点的过程中会把原图像本来是0的像素点也”抹“上一个平均值最终使得重构的图像中存在误差
②插值过程中会带来误差且周围信号稀疏高频信号有所失真导致图像边缘模糊
③在反投影图不断叠加过程中能够看到这种叠加方式会带来明显的“星状伪迹“这是造成图像边缘模糊的最重要因素。在投影数据少的时候更明显现实中投影数据的多少取决于机器。
关于第③点在数学上有研究表明反投影重构后的图像fb(x,y)和真实原图f(x,y)之间存在以下卷积关系 如果要重构成真实的图像就需要把这个由于反投影法本身带来的1/r效应去掉由于1/r会使得图片变得模糊因此1/r也称之为模糊因子。在实际操作中把1/r影响去除的方法这就是下面这种最常用的计算机、商业普遍使用CT图像重构方法滤波反投影法。 3滤波反投影法
从上面的式1可以看到1/r跟原图像是卷积关系通过傅里叶变换转变到频域上后会变成乘积关系这说明直接在频率上处理会更加简便。 在频域中式3中的r称之为权重因子其实就相当于一个滤波器想要去除模糊因子的影响重构成原始图像只需要三步
①把重构图像fb(x,y)进行二维傅里叶变换得到Fb(x,y)
②在频域中把Fb(x,y)与滤波器r相乘
③把r×Fb(x,y)进行二维傅里叶逆变换即能够得到原图f(x,y)
具体过程见下图图中的ρ就是公式中的r。 图4 直接反投影法 这种先反投影后滤波的方法需要用到两次二维傅里叶变换计算量太大因此也不是一个特别理想的矫正方法。
而滤波反投影法顾名思义就是先滤波后反投影的重构方法。具体步骤如下
①对180个投影信号先分别进行一维傅里叶变换
②在频域中对所有投影信号进行滤波也就是乘上权重因子r
③把所有滤波后的信号再进行一维傅里叶逆变换还原到时域
④对每一个已经滤波的投影信号进行反投影最后叠加这一步跟前面的直接反投影法步骤完全相同只是投影信号已经过了滤波
具体流程见下图 图5 滤波反投影法
这个方法的好处是把两次二维傅里叶变换变成了两次一维傅里叶变换计算速度大大提升。于是滤波反投影法的核心问题就变成了如何选择一个合适的滤波器r。数学中的滤波器r是一个理想化的滤波函数现实中不存在因此只能设计一个近似的滤波函数来代替r的作用。
下面常用的滤波函数有R-L滤波函数和S-L滤波函数。
①R-L滤波函数Ram-Lak滤波器 Ramp Filter 从频域图中可以看到该滤波器的作用是把低频信号减少从而突出高频信号因此经过该函数滤波后重建图像的轮廓会更清晰并且函数简单实现起来更方便。但由于该函数在频域上有一个加窗处理就是把高频部分直接截断了因此在时域中会出现有许多振荡的小尾巴截窗振荡效应这将会让重构出来的图像存在Gibb’s现象重构图像存在振荡不连续。
②S-L滤波函数Shepp-Logan滤波器 S-L滤波器在频域上不是直接加窗截断而是通过一些比较平滑的窗函数对函数进行约束。因此经过S-L滤波后重构的图像振荡更少重构的图像质量比R-L滤波器更好一些但是因为S-L在高频段并没有直接截断偏离了理想滤波器的效果因此在高频段轮廓上的重构效果没有R-L滤波器好。
可以在下图更直观地看到两个滤波器之间的区别 下面通过MATLAB编程来还原一下滤波反投影法的过程程序在最后只用复杂那个图进行展示
①对每一列投影信号分别在频域进行R-L滤波(高频信号明显突出了) ②对每个滤波后的投影信号反投影这里开始与直接反投影法一样 ③ 把所有反投影图像叠加起来 为了验证手动重构的效果下面将使用matlab中的iradon函数直接对CT投影数据进行重构把手动重构的图像和matlab自带函数重构的图像进行对比。
iradon是matlab中的拉东逆变换函数但是其实际采用的重构方法并不是对数据进行二维傅里叶逆变换而是通过滤波反投影法实现的且默认插值方法为“线性插值linear“默认滤波方法为”R-L滤波函数“这在matlab的help中可以直接看到。 仔细观察两个重构结果会发现iradon的重构结果更光滑一些那是因为手动方法只是简单粗暴地把全部离散的反变换图片叠加起来而iradon则会在最后对图像进行线性插值使得图像更连续更光滑重构效果更好。
最后把iradon重构的图像放大来看看图像细节 可以看到重构的图像其实并不连续甚至能够看出一些振荡的波纹这就是R-L滤波函数截窗振荡效应所带来的Gibb’s现象。 4MATLAB程序 %% CT反投影法
%%Part 1
%% 1、直接反投影法
% 绘制一张简单的图进行直接反投影法
clc,clear,close all
%
% 简单图
% Izeros(512,512);
% for i1:512
% for j1:512
% if ((i-256)*(i-256)(j-256)*(j-256))1024
% I(i,j)1;
% end
% end
% end
% figure,imshow(I,[]);
% %
% 复杂图
I phantom(512);
figure,imshow(I,[]);
%% 绘制代表角度的反投影图像
Rradon(I,0:179); %对物体CT扫描180°的数据每1°扫描一次
figure,imshow(R,[]);title(CT图像);
% 绘制1、45、90、135、180的投影信号和反投影图像
theta[1 45 90 135];
for i1:length(theta)rR(:,i);IIiradon([r r],[theta(i) theta(i)],linear,none)/2;%回抹反投影之前不滤波II1{i}II;figure,subplot(1,2,1),imshow(r,[]);title([num2str(theta(i)) °投影信号]);subplot(1,2,2),imshow(II,[]);title([往 num2str(theta(i)) °方向反投影]);
end
II2II1{1}II1{2}II1{3}II1{4};
figure,imshow(II2,[]);% 对反投影图像进行叠加
rR(:,1);
II_riradon([r r],[1 1])/2;
k[30 15 10 1];
for j1:4AII_r;for i2:k(j):180rR(:,i);IIiradon([r r],[i i],linear,none)/2;%回抹反投影之前不滤波AAII; end
% figure,imshow(A,[min(min(I)) max(max(I))]);figure,imshow((A),[]);
end%% Part 2
%% 2、滤波反投影法
clc,clear,close all
% 复杂图
I phantom(512);
figure,imshow(I,[]);% 对投影做傅里叶变换
Rradon(I,0:179);
widthlength(R);
% 设计R-L滤波器
filter2*[0:round(width/2-1), width/2:-1:1]/width;
% 展示滤波前的CT投影信号
figure,imshow(R,[]),title(滤波前的CT投影信号);% 每一列做傅里叶变换
r_fftfft(R,729);
% 每一列做傅里叶变换后滤波
r_fft_filterr_fft.*filter;
% 滤波后反变换(real取实部)
r_fft_filter_vreal(ifft(r_fft_filter));
% 展示滤波后的CT投影信号
% figure,imshow(r_fft_filter_v,[]),title(滤波后的CT投影信号);
figure,imshow(r_fft_filter_v),title(滤波后的CT投影信号);% 绘制1、45、90、135、180的投影信号和反投影图像
theta[1 45 90 135];
for i1:length(theta)rr_fft_filter_v(:,i);IIiradon([r r],[theta(i) theta(i)],linear,none)/2;%由于iradon默认有滤波因此关掉默认的滤波选项none直接用上面已经手动滤波的数据进行反投影II1{i}II;figure,subplot(1,2,1),imshow(r,[]);title([num2str(theta(i)) °投影信号]);subplot(1,2,2),imshow(II,[]);title([往 num2str(theta(i)) °方向反投影]);
end
II2II1{1}II1{2}II1{3}II1{4};
figure,imshow(II2,[]);% 对反投影图像进行叠加
rr_fft_filter_v(:,1);
II_riradon([r r],[1 1])/2;
k[30 15 10 1];
for j1:4AII_r;for i2:k(j):180rr_fft_filter_v(:,i);IIiradon([r r],[i i],linear,none)/2;%由于iradon默认有滤波因此关掉默认的滤波选项none直接用上面已经手动滤波的数据进行反投影AAII; end
% figure,imshow(A,[min(min(I)) max(max(I))]);figure,imshow((A),[]);
end% 上面是反投影叠加的分步计算下面直接输入全部投影数据到iradon函数中进行验证
Biradon(r_fft_filter_v,0:179,linear,none);
figure,
imshow(B,[]);title(验证滤波后图像); 以上就是对CT图像重构的整理结果更多是自己在学习过程中的理解若有表述不当的地方欢迎指出。 ------------------------------------------------------- 2021.05.27 笔记补充 ----------------------------------------------------------- 重构图像514×514与原图像512×512尺寸不一致的原因分析
首先需要感谢楼下 study_art 同学提出了这个问题以及 xjch1900 同学的提醒一开始确实没考虑过这个问题后来查看源代码找到了原因在此更新一下笔记。 ①radon函数的matlab文件解释
% Grandfathered syntax
% R RADON(I,THETA,N) returns a Radon transform with the
% projection computed at N points. R has N rows. If you do not
% specify N, the number of points the projection is computed at
% is:
%
% 2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))3
%
% This number is sufficient to compute the projection at unit
% intervals, even along the diagonal.
我们一开始先使用了radon函数来模拟了CT的扫描过程得到的结果是每个角度的投影数据。假设我们扫描的是一张正方形的图像大小为512×512如果我们要把扫描的数据都存储下来就必须设置一个足够长的数组来记录这些数据显然当扫描方向为45°角时得到的扫描数据是最长的如下图所示需要的数组长度至少为 512×√2≈725 。 回到matlab的radon函数从上面的文件解释中可看到radon的语法中的N便是用来定义这个最长长度的。但当我们没有对这个N有专门的限制时该函数就会自动计算出这个“最长长度”为了涵盖图像形状的不同情况如有的图像为长方形它设置了一个公用的计算公式
N2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))3
【其中ceil函数表示“向上取整”floor表示“向下取整”norm表示“范数计算”用在这里其实就是计算正方形斜边的长度】
如果套入这个公式计算出来的N为
size(I)[512 512];
N2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))3N 729 这比自己的理论值725大了一些但实际上这是为了涵盖所有图像情况造成的这并不会对最后的结果带来太大影响。
在正文最后的matlab程序中我并没有专门定义这个N因此radon函数自动把N通过公式计算出来了大小为729同时还定义了180个扫描角度所以得到的所有投影数据的矩阵大小为729×180. ②iradon函数的matlab文件解释
% I IRADON(R,THETA,INTERPOLATION,FILTER,FREQUENCY_SCALING,OUTPUT_SIZE)
% OUTPUT_SIZE is a scalar that specifies the number of rows and columns in
% the reconstructed image. If OUTPUT_SIZE is not specified, the size is
% determined from the length of the projections:
%
% OUTPUT_SIZE 2*floor(size(R,1)/(2*sqrt(2)))
%
% If you specify OUTPUT_SIZE, IRADON reconstructs a smaller or larger
% portion of the image, but does not change the scaling of the data.
在使用iradon函数进行图像重构时同样有一个output_size的参数可以设置表示输出图像的大小我们可以直接设置为原始图像的大小即512这样得到的图像就跟原图像大小一模一样了。但在正文最后部分我同样没有专门设置图像大小所以该函数就会自动套入一条涵盖了所有图像形状的计算公式
OUTPUT_SIZE 2*floor(size(R,1)/(2*sqrt(2)))
由于我们在第一步用radon函数模拟CT扫描时没有专门定义N得到的投影数据的矩阵R大小为729×180因此直接套入这个output_size公式后就直接算出
size(R,1)729;
OUTPUT_SIZE 2*floor(size(R,1)/(2*sqrt(2)))OUTPUT_SIZE 514所以最后得到的重构图像大小就变成了514×514. ③解决方法 使用radon函数时自己定义N。如在正文的例子中直接在radon函数中输入N的理论计算值725这样在重构iradon时就不需要设置output_size了输出结果直接为512×512使用radon函数时不定义N让函数自己计算图像的N然后在重构iradon时设置output_size的大小为原始图像大小512输出结果同样为512×512备注若radon和iradon函数都不专门设置时会出现重构图像与原图像尺寸大小不一致的情况但这并不会改变重构图像本身的缩放比例因此若没有专门的尺寸大小必须前后一致的严格要求可以不用专门设置。
