由投影重建图像:滤波反投影、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

相关文章

深入理解Mysql OnlineDDL的算法

《深入理解MysqlOnlineDDL的算法》本文主要介绍了讲解MysqlOnlineDDL的算法,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小... 目录一、Online DDL 是什么?二、Online DDL 的三种主要算法2.1COPY(复制法)

Spring的基础事务注解@Transactional作用解读

《Spring的基础事务注解@Transactional作用解读》文章介绍了Spring框架中的事务管理,核心注解@Transactional用于声明事务,支持传播机制、隔离级别等配置,结合@Tran... 目录一、事务管理基础1.1 Spring事务的核心注解1.2 注解属性详解1.3 实现原理二、事务事

Java中最全最基础的IO流概述和简介案例分析

《Java中最全最基础的IO流概述和简介案例分析》JavaIO流用于程序与外部设备的数据交互,分为字节流(InputStream/OutputStream)和字符流(Reader/Writer),处理... 目录IO流简介IO是什么应用场景IO流的分类流的超类类型字节文件流应用简介核心API文件输出流应用文

从基础到高级详解Python数值格式化输出的完全指南

《从基础到高级详解Python数值格式化输出的完全指南》在数据分析、金融计算和科学报告领域,数值格式化是提升可读性和专业性的关键技术,本文将深入解析Python中数值格式化输出的相关方法,感兴趣的小伙... 目录引言:数值格式化的核心价值一、基础格式化方法1.1 三种核心格式化方式对比1.2 基础格式化示例

redis-sentinel基础概念及部署流程

《redis-sentinel基础概念及部署流程》RedisSentinel是Redis的高可用解决方案,通过监控主从节点、自动故障转移、通知机制及配置提供,实现集群故障恢复与服务持续可用,核心组件包... 目录一. 引言二. 核心功能三. 核心组件四. 故障转移流程五. 服务部署六. sentinel部署

从基础到进阶详解Python条件判断的实用指南

《从基础到进阶详解Python条件判断的实用指南》本文将通过15个实战案例,带你大家掌握条件判断的核心技巧,并从基础语法到高级应用一网打尽,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一... 目录​引言:条件判断为何如此重要一、基础语法:三行代码构建决策系统二、多条件分支:elif的魔法三、

Python WebSockets 库从基础到实战使用举例

《PythonWebSockets库从基础到实战使用举例》WebSocket是一种全双工、持久化的网络通信协议,适用于需要低延迟的应用,如实时聊天、股票行情推送、在线协作、多人游戏等,本文给大家介... 目录1. 引言2. 为什么使用 WebSocket?3. 安装 WebSockets 库4. 使用 We

从基础到高阶详解Python多态实战应用指南

《从基础到高阶详解Python多态实战应用指南》这篇文章主要从基础到高阶为大家详细介绍Python中多态的相关应用与技巧,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录一、多态的本质:python的“鸭子类型”哲学二、多态的三大实战场景场景1:数据处理管道——统一处理不同数据格式

MySQL数据类型与表操作全指南( 从基础到高级实践)

《MySQL数据类型与表操作全指南(从基础到高级实践)》本文详解MySQL数据类型分类(数值、日期/时间、字符串)及表操作(创建、修改、维护),涵盖优化技巧如数据类型选择、备份、分区,强调规范设计与... 目录mysql数据类型详解数值类型日期时间类型字符串类型表操作全解析创建表修改表结构添加列修改列删除列

Python 函数详解:从基础语法到高级使用技巧

《Python函数详解:从基础语法到高级使用技巧》本文基于实例代码,全面讲解Python函数的定义、参数传递、变量作用域及类型标注等知识点,帮助初学者快速掌握函数的使用技巧,感兴趣的朋友跟随小编一起... 目录一、函数的基本概念与作用二、函数的定义与调用1. 无参函数2. 带参函数3. 带返回值的函数4.