由投影重建图像:滤波反投影、FDK、TFDK三维重建算法理论基础

本文主要是介绍由投影重建图像:滤波反投影、FDK、TFDK三维重建算法理论基础,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

1. 基础理论从:
[1] RafaelC.Gonzalez, RichardE.Woods, Gonzalez,等. 阮秋琦等译.数字图像处理(第三版)[M]. 电子工业出版社, 2011.P232

[2] RafaelC.Gonzalez, RichardE.Woods, StevenL.Eddins. 阮秋琦译.数字图像处理:MATLAB版:本科教学版[M]. 电子工业出版社, 2014. (第二版)P102


参考文献[1]从P228-P245

了解:

1.    由投影重建图像

1.     计算机断层(CT)的原理

2.     投影和Radon(雷登)变换(radontransform)

3.     正弦图(sinogram)和Shepp-Logan幻影

4.     傅里叶切片定理

5.     使用平行射线束滤波反投影的重建

6.     Ram-Lak滤波器与Hamming、Hann等窗函数

7.     使用扇形射线束滤波反投影的重建


1. 由投影重建图像

利用Matlab来复现参考文献[1]中的处理流程:

首先生成一幅原始图像:

%代码如下:
I=zeros(512,512);
for i=1:512for j=1:512if ((i-256)*(i-256)+(j-256)*(j-256))<1024I(i,j)=1;endend
end
imshow(I,[]);


得到1°和90°两个方向上的投影(一维),并回抹得到反投影图像(二维)

         

这两个方向的反投影相叠加:


%代码如下:(这里参考了MATLAB自带的函数说明输入:doc iradon打开该参考说明)
%就是怎么得到特定角度下的一个投影,并形成反投影图像:
R=radon(I,0:179);
r1=R(:,2);
II1=iradon([r1 r1],[1 1])/2; %得到1°角度下的投影图像
imshow(II1);r90=R(:,91);
II90=iradon([r90 r90],[90 90])/2; %得到90°角度下的投影图像
imshow(II90);II1_90=II1+II90;%二者叠加
imshow(II1_90);

同理,得到45°和135°方向的反投影图像,并将这四幅反投影图像相叠加:

%代码如下:
r45=R(:,46);
II45=iradon([r45 r45],[45 45])/2;
imshow(II45);r135=R(:,136);
II135=iradon([r135 r135],[135 135])/2;
imshow(II135);II_1_90_45_135=II1+II45+II90+II135;%四个反投影图像相叠加!
imshow(II_1_90_45_135);

当反投影的角度采样增多时,这里从0°到179°每隔1°共采180个反投影,相叠加后形成反投影图像:(不使用滤波器时的反投影重建图像如下:


%这里用了iradon函数,线性插值,不使用滤波器。
ii2=iradon(R,0:179,'linear','none');
imshow(ii1,[]);

ii1=iradon(R,0:179,'linear','Ram-Lak');%默认使用线性插值,Ram-Lak滤波器!

该结果如下:


将重建后的图像与原图像进行对比:(左边是原图、右边是重建得到的图像,注意模糊与振铃现象)

  


我们再使用MATLAB自带的一头部幻影图像Shepp-Logan:

原图:


未使用滤波器的重建图像:


使用Ram-Lak滤波器  :


使用Shepp-Logan滤波器:


%生成以上图像的代码如下:
P=phantom(512);
theta=0:179;
[R,xp]=radon(P,theta);
I1P=iradon(R,0:179,'linear','none');
I2P=iradon(R,0:179,'linear','Ram-Lak');
I3P=iradon(R,0:179,'linear','Hamming');


以下内容参考文献[1]将滤波反投影的基础知识过一遍:


笛卡尔坐标系的一条直线由它的斜截式描述: ;或由其法线方程来描述:




平行射线束的投影可以由一组这样的直线建模。如图5.37所示,投影信号中的任意一点由沿着直线 的射线和给出。连续变量的情况下,线求个变为线积分,由下式给出:



该式是沿xy平面内任意一条直线的f(x,y)的投影(线积分)的公式,就是雷登(Radon)变换

符号R{f(x,y)}或R{f}有时用于代替(3)式中的 来表示f的雷登变换。雷登变换是由投影重建图像的基石,计算机断层(CT)是其在图像处理领域的主要应用。在离散情况下,(3)式变为


正弦图与雷登变换:


    







正弦图包含了重建图像f(x,y)所需的信息。

正弦图的视觉分析仅限于实际应用、但有时对于算法开发是有帮助的。

 

CT的关键目的是从投影得到物体的三维表示。其方法是反投影每一个投影,然后对反投影求和以产生一幅图像(切片),再堆积所有的结果以产生三维物体的再现。(这里是指用二维切片堆积成三维体,与FDK/TFDK直接就是三维重建是不同的!)


5.11.4傅里叶切片定理

 

傅里叶切片定理即投影的一维傅里叶变换和被投影区域图像的二维傅里叶变换间的关系。

投影 的一维傅里叶变换为:




式(11)就是著名的傅里叶切片定理(或投影切片定理)。它说明了一个投影(一维)的傅里叶是得到这个投影的二维区域f(x,y)的二维傅里叶变换对应角度下的一个切片。正如图5.41所示,任意一个投影的一维傅里叶变换可以沿着一个角度提取一条直线的F(u,v)的值来得到,而该角度就是投影时所用的角度。


下面推导滤波反投影公式,将用到傅里叶切片定理。

F(u,v)的反傅里叶变换为:






当c=0.54时,该函数称为汉明窗(RichardHamming);

当c=0.5时,称为韩窗(Juliusvon Hann)

加了窗函数的滤波器在空域的振铃现象减弱。

我们可以预期,由于使用汉明窗的反投影有较小的振铃,但稍微模糊一点。见下图:


原图  ;使用Ram-Lak滤波器  ;  使用Hamming窗加窗后的滤波器

        

CT的多数应用中(特别是医学上),像振铃这样的人为缺陷有严重的厉害关系,使其最小化是有意义的工作。调整滤波算法(可以做文章的地方,一些硕士论文就这样自创新的滤波器,发现效果有所改进,好,成文!比如文献

[1] 张銮. 基于平板探测器的锥束CT重建技术研究[D]. 中北大学, 2010.)、和硬件制造方面的改进(如,提供探测器的检测细腻度,即采样粒度)




因为斜坡滤波器(甚至在被加窗时)在频率域的直流项为零,故每一幅反投影图像的均值将为零。这将意味着,每一幅反投影图像都将有正像素和负像素值,当所有的反投影图像相加形成最终的重建图像时,一些负像素值位置可能变成正像素,而平均值可能不为零,但是,典型地,最终的图像将还是有负像素值。

当有关平均值的知识未知时,就使用标定的方法将图像的像素值都归一化到一个区间[0,255]。当典型的平均值的知识是可用的时候,可将该值加到频率域的滤波器上,从而抵消斜波并防止直流项为零。当在空域中使用卷积时,截断空间滤波器的长度(斜坡的反傅里叶变换)的真正效果都将防止其有零均值,这样就完全避免的迫零问题。

滤波反投影算法:

1.     计算每一个投影的一维傅里叶变换;

2.     用滤波函数 乘以每一个傅里叶变换;

3.     得到每一个滤波后的变换的一维反傅里叶变换;

4.     对步骤3得到的所以一维反变换积分(求和)








这篇关于由投影重建图像:滤波反投影、FDK、TFDK三维重建算法理论基础的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/286566

相关文章

不懂推荐算法也能设计推荐系统

本文以商业化应用推荐为例,告诉我们不懂推荐算法的产品,也能从产品侧出发, 设计出一款不错的推荐系统。 相信很多新手产品,看到算法二字,多是懵圈的。 什么排序算法、最短路径等都是相对传统的算法(注:传统是指科班出身的产品都会接触过)。但对于推荐算法,多数产品对着网上搜到的资源,都会无从下手。特别当某些推荐算法 和 “AI”扯上关系后,更是加大了理解的难度。 但,不了解推荐算法,就无法做推荐系

基于人工智能的图像分类系统

目录 引言项目背景环境准备 硬件要求软件安装与配置系统设计 系统架构关键技术代码示例 数据预处理模型训练模型预测应用场景结论 1. 引言 图像分类是计算机视觉中的一个重要任务,目标是自动识别图像中的对象类别。通过卷积神经网络(CNN)等深度学习技术,我们可以构建高效的图像分类系统,广泛应用于自动驾驶、医疗影像诊断、监控分析等领域。本文将介绍如何构建一个基于人工智能的图像分类系统,包括环境

康拓展开(hash算法中会用到)

康拓展开是一个全排列到一个自然数的双射(也就是某个全排列与某个自然数一一对应) 公式: X=a[n]*(n-1)!+a[n-1]*(n-2)!+...+a[i]*(i-1)!+...+a[1]*0! 其中,a[i]为整数,并且0<=a[i]<i,1<=i<=n。(a[i]在不同应用中的含义不同); 典型应用: 计算当前排列在所有由小到大全排列中的顺序,也就是说求当前排列是第

csu 1446 Problem J Modified LCS (扩展欧几里得算法的简单应用)

这是一道扩展欧几里得算法的简单应用题,这题是在湖南多校训练赛中队友ac的一道题,在比赛之后请教了队友,然后自己把它a掉 这也是自己独自做扩展欧几里得算法的题目 题意:把题意转变下就变成了:求d1*x - d2*y = f2 - f1的解,很明显用exgcd来解 下面介绍一下exgcd的一些知识点:求ax + by = c的解 一、首先求ax + by = gcd(a,b)的解 这个

Open3D 基于法线的双边滤波

目录 一、概述 1.1原理 1.2实现步骤 1.3应用场景 二、代码实现 2.1关键函数 输入参数: 输出参数: 参数影响: 2.2完整代码 三、实现效果 3.1原始点云 3.2滤波后点云 Open3D点云算法汇总及实战案例汇总的目录地址: Open3D点云算法与点云深度学习案例汇总(长期更新)-CSDN博客 一、概述         基于法线的双边

综合安防管理平台LntonAIServer视频监控汇聚抖动检测算法优势

LntonAIServer视频质量诊断功能中的抖动检测是一个专门针对视频稳定性进行分析的功能。抖动通常是指视频帧之间的不必要运动,这种运动可能是由于摄像机的移动、传输中的错误或编解码问题导致的。抖动检测对于确保视频内容的平滑性和观看体验至关重要。 优势 1. 提高图像质量 - 清晰度提升:减少抖动,提高图像的清晰度和细节表现力,使得监控画面更加真实可信。 - 细节增强:在低光条件下,抖

2024年流动式起重机司机证模拟考试题库及流动式起重机司机理论考试试题

题库来源:安全生产模拟考试一点通公众号小程序 2024年流动式起重机司机证模拟考试题库及流动式起重机司机理论考试试题是由安全生产模拟考试一点通提供,流动式起重机司机证模拟考试题库是根据流动式起重机司机最新版教材,流动式起重机司机大纲整理而成(含2024年流动式起重机司机证模拟考试题库及流动式起重机司机理论考试试题参考答案和部分工种参考解析),掌握本资料和学校方法,考试容易。流动式起重机司机考试技

【数据结构】——原来排序算法搞懂这些就行,轻松拿捏

前言:快速排序的实现最重要的是找基准值,下面让我们来了解如何实现找基准值 基准值的注释:在快排的过程中,每一次我们要取一个元素作为枢纽值,以这个数字来将序列划分为两部分。 在此我们采用三数取中法,也就是取左端、中间、右端三个数,然后进行排序,将中间数作为枢纽值。 快速排序实现主框架: //快速排序 void QuickSort(int* arr, int left, int rig

零基础学习Redis(10) -- zset类型命令使用

zset是有序集合,内部除了存储元素外,还会存储一个score,存储在zset中的元素会按照score的大小升序排列,不同元素的score可以重复,score相同的元素会按照元素的字典序排列。 1. zset常用命令 1.1 zadd  zadd key [NX | XX] [GT | LT]   [CH] [INCR] score member [score member ...]

poj 3974 and hdu 3068 最长回文串的O(n)解法(Manacher算法)

求一段字符串中的最长回文串。 因为数据量比较大,用原来的O(n^2)会爆。 小白上的O(n^2)解法代码:TLE啦~ #include<stdio.h>#include<string.h>const int Maxn = 1000000;char s[Maxn];int main(){char e[] = {"END"};while(scanf("%s", s) != EO