滤波反投影
滤波反投影法的实施步骤
滤波反投影法的实施步骤1. 简介滤波反投影法（Filtered Back Projection, FBP）是一种重建计算机体层成像（Computed Tomography, CT）图像的方法。
它通过对一系列飞行时间或射线投影数据进行滤波和反投影操作，来恢复被测物体的内部结构信息。
本文将详细介绍滤波反投影法的实施步骤。
2. 数据获取首先，我们需要获取一系列的投影数据。
投影数据是通过将被测物体从不同角度进行透射扫描而得到的，通常使用X射线或者其他形式的射线。
3. 数据预处理在进行滤波反投影之前，我们需要对投影数据进行一些预处理操作。
这些操作包括：去除背景噪声、校准投影数据以及执行各种纠正操作，以确保数据的准确性和一致性。
4. 滤波操作滤波操作是滤波反投影法中的关键步骤，它通过对投影数据进行频域滤波来增强图像的对比度和细节信息。
常用的滤波方法包括：Ram-Lak、Shepp-Logan、Hann窗等。
5. 反投影操作在经过滤波操作后，我们需要进行反投影操作来恢复原始图像。
反投影是将滤波后的投影数据逆向投影回图像空间的过程。
反投影操作涉及到几何重建算法和数学运算，可以使用快速反投影算法（FFT-based Back Projection）等高效算法来加快图像重建速度。
6. 重建图像优化得到初步的重建图像后，我们可以对其进行一些优化处理，以提高图像质量和视觉效果。
常见的优化方法包括：去除伪影、降噪、增加对比度、增强图像细节等。
7. 结果评估最后，我们需要对重建图像进行结果评估。
这包括比较重建图像和原始物体的差异，评估图像质量和准确性，并根据需要进行后续的调整和改进。
8. 总结滤波反投影法是一种常用的重建CT图像的方法，它通过对投影数据进行滤波和反投影操作，实现了对被测物体的内部结构信息的恢复。
本文介绍了滤波反投影法的实施步骤，包括数据获取、数据预处理、滤波操作、反投影操作、重建图像优化和结果评估等。
通过遵循这些步骤，可以有效地进行CT图像的重建工作，并获得高质量的图像结果。
滤波反投影法
滤波反投影法：
滤波反投影法根据附件三所给接收信息，采用先修正、后投影重建图像的做法，可得到原始图像的吸收率信息。
其原理为：在得到某一角度下的投影函数（一维函数）后，对此函数做滤波处理，得一修正后的滤波函数，再对修正后的滤波函数做反投影运算，得待检测介质吸收率在正方形托盘中的每一点的分布密度函数。
图1给出了滤波反投影法重建原始图像的流程图。
图1滤波反投影法流程图
反投影法重建原始图像的步骤：
（1）在对应于投影函数的角度下对投影函数做一维Fourier变换；
（2）对（1）得到的变换结果乘以权重因子；
（3）对（2）加权后得到的结果做一维傅立叶；
（4）对（3）所得函数做直接反投影；
（5）改变投影角度，得到180个不同的投影角度，对每一角度，重复上述步骤（1）~（4）。
R-L（Ram-Lak）滤波函数：
此函数的基本条件是二维图像函数的频率是有界的，显然，此题所得附件五的所有数据满足此条件。
故频域中的滤波函数可表示为：
其函数图像如图1.
图1R-L滤波函数图像
连续的R-L卷积函数所得结果为：
离散的R-L卷积函数所得结果为：
根据上述滤波原理，在本题中，对附件五中数据的具体滤波过程可用Matlab内置的Ram-Lak命令实现。
基于2017数学建模的滤波反投影算法应用
Image & Multimedia Technology •图像与多媒体技术Electronic Technology & Software Engineering 电子技术与软件工程• 91【关键词】CT 重构 randon 变换 滤波反投影1 CT图像重建原理的知识背景CT 系统基本过程是：平行入射的X 射线垂直于探测器平面发射，形成一个发射-接收CT 系统，每个探测器单元都看做是一个接收点，且间隔距离相等。
计算机断层成像图像重建的过程是按照一定的算法将已经检测到的投影数据进行数学运算，最终得到断层图像。
Radon 变换及其逆变换：物体断层被射线扫描后需要用重建算法计算才能得到CT 图像，图像重建的基础是Radon 变换及其逆变换。
假设每条射线相互平行，对于一个二维平面进行射线检测可得到一条投影数据，该投影数据称为二维平面的一个Radon 变换；如果检测中该平面旋转180度，同时将对应的投影数据进行组合，则得到类似正弦分布形式的图像，从正弦图获取二维平面图像的变换称为Radon 反演。
用公式可分别描述为：，由于matlab 中封装有radon 函数，使用时直接调用函数：R=radon （I ，theta ）。
2 滤波反投影算法radon 函数使用的算法是滤波反投影法，反投影算法因为引入“星”状伪影而导致重建的图像失真，为了消除这个伪影，在进行反投影重建之前将数据修正，最后对修正后的投影数据进行反投影，这样就获得没有伪影的重建图像。
该方法是在空间域中把投影的数据直接反向投射到需要重建的图像中，然后将逐个的反投影图像累加起来。
滤波反投影法基本实现步骤：对数据作一维傅里叶变换→滤波函数：R-L 函数→对滤波后的数据作傅里叶逆变换→反投影求图像函数。
本文简要介绍推导傅里叶变换的过程：令为f 的二维傅里叶变换.单变量函G φ(ω)F(ω cosφ,ω sin φ )为通过φ角的F 切片，并记g φ (p)基于2017数学建模的滤波反投影算法应用文/李春梅为由合成方程 确定的函数，则 (Ff φV )(ω)=F(ω cos φ,ω sin φ)，其中F 是单变量傅里叶变换算子，它建立了Radon 变换和傅里叶变换的联系.然后采用极坐标u=ω cos φ,v=ω sin φ表示傅里叶合成公式得将这个积分分解成两个积分式，通过变换、合并，最后使用投影切片定律重写这个积分形式为：f(x,y)= d ω d φ由此得到合成方程。
滤波反投影法
滤波反投影法：
滤波反投影法根据附件三所给接收信息，采用先修正、后投影重建图像的做法，可得到原始图像的吸收率信息。
其原理为：在得到某一角度下的投影函数（一维函数）后，对此函数做滤波处理，得一修正后的滤波函数，再对修正后的滤波函数做反投影运算，得待检测介质吸收率在正方形托盘中的每一点的分布密度函数f(x,y)。
图1给出了滤波反投影法重建原始图像的流程图。
图1滤波反投影法流程图
反投影法重建原始图像的步骤：
（1） 在对应于投影函数的角度下对投影函数做一维Fourier 变换；
（2） 对（1）得到的变换结果乘以权重因子|ρ|；
（3） 对（2）加权后得到的结果做一维傅立叶；
（4） 对（3）所得函数做直接反投影；
（5） 改变投影角度，得到180个不同的投影角度，对每一角度，重复上述步骤（1）
~（4）。
R-L （Ram-Lak ）滤波函数：
此函数的基本条件是二维图像函数的频率是有界的，显然，此题所得附件五的所有数据满足此条件。
故频域中的滤波函数可表示为：
G (ρ)={|ρ|， |ρ|≤ρ0 0， 其它
其函数图像如图1.
图1R-L 滤波函数图像
连续的R-L 卷积函数所得结果为：
g (R )=ρ02[2sin c (2ρ0R )−sin c 2(ρ0R )]
离散的R-L 卷积函数所得结果为：
g (nT )={ 14T 2 ， n =0 0 ， n 为偶数−1n 2π2T 2，n 为奇数
根据上述滤波原理，在本题中，对附件五中数据的具体滤波过程可用Matlab 内置的Ram-Lak 命令实现。
ct成像fbp算法公式
ct成像fbp算法公式CT（计算机断层扫描）成像中的FBP（滤波反投影）算法是一种常用的重建方法。
在FBP算法中，首先对投影数据进行滤波，然后对滤波后的数据进行反投影，从而得到重建的图像。
以下是FBP算法的简要步骤：1. 采集投影数据：在CT扫描过程中，X射线源沿某一方向旋转，通过物体照射到探测器上，得到一系列投影数据。
2. 滤波：对投影数据进行滤波，以消除噪声和伪影。
常用的滤波方法有理想低通滤波、高斯滤波等。
3. 反投影：将滤波后的投影数据按照一定的角度间隔进行重新采样，然后对采样数据进行反投影运算，得到重建的图像。
4. 图像重建：对反投影结果进行图像重建，得到最终的CT图像。
关于FBP算法的公式，以中心切片法为例，可以分为以下几个部分：1. 滤波：对投影数据P(x,y)进行滤波，得到滤波后的投影数据P'(x,y)。
滤波公式如下：P'(x,y) = ∫P(x-Δx, y-Δy) * h(Δx, Δy) dx dy其中，h(Δx, Δy)是滤波函数，Δx和Δy分别表示x和y方向上的偏移量。
2. 反投影：将滤波后的投影数据P'(x,y)按照一定的角度间隔进行重新采样，得到采样数据P''(u,v)。
反投影公式如下：u = x -Δx * cos(θ)v = y -Δy * cos(θ)其中，θ是投影方向与水平方向的夹角。
3. 图像重建：对采样数据P''(u,v)进行插值，得到重建的图像I(x,y)。
插值方法有线性插值、双线性插值等。
需要注意的是，FBP算法在实际应用中可能会受到噪声、伪影等因素的影响。
为了提高图像质量，可以采用其他重建方法，如迭代重建算法（如ART、MLEM等）或模型驱动的重建方法（如MBIR等）进行优化。
滤波反投影法迭代方程
滤波反投影法是一种用于图像重建的算法，其迭代方程通常由以下步骤组成：
1. 对当前投影图像进行滤波操作，以去除噪声和伪影。
2. 将滤波后的投影图像进行反投影，得到重建图像的更新值。
3. 将更新值与前一次迭代的重建图像进行叠加，得到新的重建图像。
4. 重复步骤1-3，直到达到预设的迭代次数或满足收敛条件。
具体来说，滤波反投影法的迭代方程可以表示为：
\(I^{k+1} = I^k + \lambda \left( \text{滤波后的投影图像} - \text{反投影图像} \right)\)
其中，\(I^{k+1}\)表示第\(k+1\)次迭代的重建图像，\(I^k\)表示第\(k\)次迭代的重建图像，\(\lambda\)是控制迭代的步长，\(\text{滤波后的投影图像}\)是滤波后的投影图像，\(\text{反投影图像}\)是反投影得到的图像。
需要注意的是，具体的迭代方程可能会因不同的滤波器和反投影方法而有所不同。
基于滤波反投影算法的CT系统成像研究
基于滤波反投影算法的CT系统成像研究摘要：CT系统的安装会使得旋转中心发生偏离，从而影响成像质量，因此需要借助于已知结构的样品来标定CT系统的参数，并且利用标定的参数对未知结构的样品进行图像重建。
首先根据直接反投影算法和滤波反投影算法对收集到的数据中的接收信息分别进行图像重建，通过成像图像可知，滤波反投影算法更优；旋转中心可能发生偏移以及CT系统具有初始角度，依次进行旋转、平移、裁剪和残影去除操作，来校正投影图像，从而得到较高质量的图像。
关键词：CT成像原理（影像医学与核医学）；滤波反投影法；图像重建；吸收率引言CT(Computed Tomography)是用X线束从多个方向对人体检查部位具有一定厚度的层面进行扫描，由探测器而不用胶片接收透过该层面的X线，转变为可见光后，由光电转换器转变为电信号，再经模拟/数字转换器转为数字，输人计算机处理。
数字矩阵中的每个数字经数字/模拟转换器转为由黑到白不等灰度的小方块，称之为像素，并按原有矩阵顺序排列，即构成CT图像。
所以，CT图像是由一定数目像素组成的灰阶图像，是数字图像，是重建的断层图像。
首先根据直接反投影算法和滤波反投影算法对收集到的数据中的接收信息分别进行图像重建，将图像重建[4-6]的两种结果进行对比，得出效果较好的模型；然后，旋转中心可能发生偏移以及CT系统具有初始角度，通过旋转、平移、裁剪和残影去除等操作来校正投影图像，最后对图像进行标准化调整，从而提高了成像质量。
1 模型的准备与建立1.1 CT成像的数学基础Rand变换如图1所示，直线g是xOy平面内任意一条直线，t是原点到直线g的距离，φ为原点到直线g的垂线与x轴的夹角。
对于xOy平面内任意一条直线可以由(t,φ)唯一确定。
二维平面中函数f(x,y)沿着直线的积分等于其Rand变换。
中心切片定理中心切片定理是CT图像重建算法的基础，在非衍射源情况下，含义是图像在某个视角下平行投影的一维Fourier变换等同于该图像二维Fourier变换的一个中心切片。
滤波反投影
滤波反投影法重建CT 图像实验指导书一、 实验目的1． 了解傅立叶变换法、直接反投影法重建CT 图像的原理；2． 掌握滤波反投影法重建CT 图像的原理和基本方法。
二、 实验器材装有MATLAB 程序的PC 机，滤波反投影法图像重建演示软件，投影数据。
三、 实验原理CT 图像重建问题实际上就是如何从投影数据中解算出成像平面上各像素点的衰减系数。
图像重建的算法有多种，如反投影法、傅立叶变换法、迭代法、滤波反投影法等。
在介绍算法前，有必要先介绍从投影重建图像的重要依据，即中心切片定理。
1． 中心切片定理密度函数(,)f x y 在某一方向上的投影函数()g R θ的一维傅立叶变换函数()g θρ是原密度函数(,)f x y 的二维傅立叶变换函数(,)F ρθ在(,)ρθ平面上沿同一方向且过原点的直线上的值。
图1 中心切片定理2．傅立叶变换法如果在不同角度下取得足够多的投影函数数据，并作傅立叶变换，根据中心切片定理，变换后的数据将充满整个(,)u v 平面。
一旦频域函数(,)F u v 或(,)F ρβ的全部值都得到后，将其做傅立叶反变换，就能得到原始的密度函数(,)f x y ，即所要重建的图像。
上述图像重建算法称为傅立叶变换法，图2给出了傅立叶变换重建方法的流程图。
图中指出，对于每次测得的投影数据先作一维傅立叶变换。
根据中心切片定理，可将此变换结果看成二维频率域中同样角度下过原点的直线上的值。
在不同投影角下所得的一维变换函数可在频域中构成完整的二维傅立叶变换函数，将此二维变换函数做一次逆变换，就得到了所要求的空间域中的密度函数。
为了在二维逆变换中采用快速傅立叶变换算法，通常在逆变换前要将极坐标形式的频域函数变换成直角坐标形式的数据。
图2 傅立叶变换重建图像的过程采用傅立叶变换法重建图像时，投影函数的一维傅立叶变换在频域中为极坐标形式，把极坐标形式的数据通过插补运算转换为直角坐标形式的数据时，计算工作量较大。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
平行束滤波反投影1100500121 赵伟伦 准备知识:一维Fourier 变换：dt et f f f F t i ⎰+∞∞--⋅==πωω2)()(~)( 一维逆Fourier 变换： ωωπωd e f f F x f t i ⎰+∞∞--⋅==21)(~)~()( 且有：)~(~),(11f F F f f F F f --⋅=⋅=重要的性质：（卷积特性）)(~)(~)*(ωωgf g f F ⋅=； )(~)(~)(ωωgf g f F *=⋅ 二维Fourier 变换： dX e x x f f f F x x i R ),(),(22121221212),(),(~)(⋅-⎰==ωωπωω； 逆二维Fourier 变换： Ω==⋅-⎰d e f f F x x f x x i R ),(),(221122121212),(~)~(),(ωωπωω； 中心切片定理：),)(ˆ()(2ϕωωfF f F r =Φ, 其中),(ˆϕr f 是),(21x x f 的Radon 变换： 解释：一个二元函数的Radon 变换关于r 的一维Fourier 变换与这个二元函数的二维Fourier 变换形式相等。
滤波反投影：思路：)(),(121f F F x x f ⋅=-()()[][]ϕϕωωϕωϕωϕωωϕωϕωϕωωωϕωωϕωϕωωϕωϕωωωϕωωωππωωππωωππωωππωωπd r f F r d fF F d d e fF x x r d d e fF d d e f F d d e f d d e f F X r x x r r r r i r x x i r x x i rx x i x x i R Φ⋅=-Φ⋅=-∞+∞-⋅∞+∞-⋅∞+⋅∞+⋅*⇔=⋅⇔⇔Φ⋅=Φ=⇔⇔⇔⇔⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰)(H ),(ˆfourier fourier ),()(H ),)(ˆ(]),)(ˆ([),),(),(),(),)(ˆ(),)(ˆ()(~)(1),(1202121),(),(20),(),(2200),(),(2200221),(),(222121212121212121212变化变化等于函数点乘后的个函数的卷积的并根据卷积的性质：两设旋转角为为坐标映射到探测器上，设为用极坐标方式表示出来（把，可知），（由于中心切片定理)(),(~),(r H r f r G *=ϕϕ)(r H 是滤波器总结：ϕϕϕωωϕωππωπd r H r fd def F X f X r X r r i r Φ⋅=Φ⋅=+∞∞-⎰⎰⎰=⎥⎦⎤⎢⎣⎡⋅=)(*),(ˆ),)(ˆ()(020 解释为：投影数据),(ˆϕr f 先进行滤波)(*),(ˆr H r f ϕ 在对滤波数据进行投影ϕϕπd r H r f X r Φ⋅=⎰)(*),(ˆ0简单例子：（大圆与小圆）通过已得到的正投影‘round.dat’经过滤波后，反投影后的图像。
正投影数据：滤波图像：反投影后的图像：总的滤波反投影过程：1，得到图像的正投影2，滤波（投影与滤波器卷积）3，反投影（一个像素所有角度的滤波函数值的和再乘以pi/views ）views r H r P d r H r P x x f viewsπϕϕϕϕπ∆=*=∑⎰=)(*),()(),(),(0021例子：GUI界面实现：修改前：程序讲述：主程序：读入一张512*512的黑白图片，展示图片；求正投，展示正投结果；求滤波后的结果并展示；求反投影后的的结果并展示。
view=360;bins=512;U=1.5; %U=view_radiusP=zeros(view,bins);clc;clear;A = imread('angle.jpg');width = 512;height = 512;A = im2double(A);figure,imshow(A)t1 = -1;t2 = 1;n= 2048;dt = (t2-t1)/n;bins = 512;view = 360;du = 3/(bins-1);P = zeros(view,bins);time = clock;for i = 1:viewphi = (i - 1)*(pi/180);for j = 1:bins%将视野范围512等分（定义512个接受器）r = (j - 1)*du - 1.5 + 0.5*du;%设定第一个接受器的与探测器中心的距离int = integral(r,phi,t1,t2,dt,width,height,A);P(i,j) = int;endendouttime = etime(clock,time);s_png2=sprintf ('angle.dat') ;s_png2=strcat ('D:\文件\ct\zwl\',s_png2) ;fp = fopen(s_png2,'wb');fwrite(fp,P','float64');fclose(fp);s_png = sprintf('angle.dat');%angles_png = strcat('D:\文件\ct\zwl\',s_png);fp = fopen(s_png,'r');P = fread(fp,view*bins,'float64');fclose(fp);P = reshape(P,bins,view);P = P';U=1.5;cellsize=2*U/bins;theta0=2*pi/view;figure, imshow(P);integ=RL_filter(P,view,bins,cellsize);figure, imshow(integ,[]);Image=back_projection_circle(integ,cellsize,view,theta0,U);figure,imshow(Image,[])s_png2=sprintf ('R_L_projection1.dat') ;s_png2=strcat ('',s_png2) ;fp = fopen(s_png2,'wb');fwrite(fp,Image','float64');fclose(fp);子函数1：备注：将图像积分区域转化为二维图像中的点的方法：function int = integral(r,phi,t1,t2,dt,width,height,A)int = 0;for t = t1:dt:t2x = r*cos(phi) - t*sin(phi);y = r*sin(phi) + t*cos(phi);%探测器与中心距离为r 时，这条探测器上从-1到1共2048个点， %转变为旋转角phi 下的坐标%将积分区间上的点转化为二维图像中的点%把坐标轴上的点的坐标转化为图像上的位置ϕϕϕϕϕϕϕϕcos sin )sin (cos )cos sin ()sin (cos t r y t r x t r l t r l +=-+=+-++=Φ+Φ=⊥tx = 256*(sqrt(2)*x + 1) + 0.5; ty = 256*(1 - sqrt(2)*y) + 0.5;5122*5.0225122*)1(5122*5.0225122*)1(-+--=+--=dy y dx x fx = floor(tx);fy = floor(ty);%把矩阵上的信息映射到图像上 ϕ%主要用差值法：%非边界点时用差值法，对于一个像素里的点用4的顶点的值去近似此点的值%if (tx<1 || ty<1 || tx>width || ty>height)mu = 0;elseif (tx == width && ty ~= height)mu = (fy+1 - ty)*(fx+1 - tx)*A(fx,fy) + (ty - fy)*(fx+1 - tx)*A(fx,fy+1);elseif (tx ~= width && ty == height)mu = (fy+1 - ty)*((fx+1 - tx)*A(fx,fy) + (tx - fx)*A(fx+1,fy));elseif (tx == width && ty == height)mu = A(fx,fy);elsemu1 = (fx+1 - tx)*A(fx,fy) + (tx - fx)*A(fx+1,fy);mu2 = (fx+1 - tx)*A(fx,fy+1) + (tx - fx)*A(fx+1,fy+1);mu = (fy+1 - ty)*mu1 + (ty - fy)*mu2;endint = int + dt*mu;End%设定足够长的滤波函数%P(i,:)*filter的长度分别为bins和2*bins-1；卷积的结果为的长度为3*bins-2,。
选取其中% [bins,2*bins]的范围内的值子函数2：function integ=RL_filter(P,view,bins,cellsize)filter_length = 2 * bins - 1;for i = 0:filter_length-1tmp = i - (filter_length - 1) / 2;if mod(tmp,2)==0filter(i+1) = 0.0;elsefilter(i+1) = -1.0 / (tmp * tmp * pi * pi) / cellsize;filter((filter_length - 1) / 2 + 1) = 0.25 / cellsize;endendfor i = 1 : viewP_filtered(i,:)=conv(P(i,:), filter);endinteg=P_filtered(:,bins:2*bins);子函数3：function Image=back_projection_circle(integ,cellsize,view,theta0,U)for height=1:512for width=1:512ee=0;for i=1:viewrr=(-U+(width-1)*cellsize+0.5*cellsize)*cos(i*theta0)+(U-(height-1)*cellsize-0.5*cellsize)*sin(i*theta0);%将图像上的点投影到接受器上，在接受器上的坐标（地址：有正有负）。