下载文档原格式
实验一投影数据采集与滤波反投影重建一、实验目的及要求：用程序模拟X射线的投影，获得Shepp-Logan模型的投影数据。
对获得的投影数据进行滤波反投影重建，获得Shepp-Logan模型的重建图像。
二、实验基本原理：X射线穿过人体时，人体的各种组织对X射线有不同程度的衰减，即不同的组织有不同的线性衰减系数μ。
假设强度为I0的X射线穿过均匀分布衰减系数为μ的物体，行进了x的距离，强度变为I，按Beer定理有或①若物体时分段均匀的，系数分别是μ1、μ2、μ3、...,相应的长度为x1，x2，x3，...,则下式成立：②更一般的可用下面的积分式表示：③由于只是模拟X射线的投影过程，我们简化了问题。
假设断面的结构如图1.1（Shepp-Logan）所示（各图元均为椭圆），各个椭圆表示了人体的不同的组织（内部是均匀的），分别有不同的线性衰减系数μ。
那么，就可利用公式②来求某条X射线投影值。
各个椭圆（组织）的线性衰减系数μ（Shepp-Logan图的各椭圆的位置、大小和线性衰减系数参见表1.1）是已知的，问题就是球X射线穿过椭圆时的行进距离。
设椭圆的长短轴为a，b；X射线与b的夹角为Φ；椭圆中心到X射线的距离为d。
如图1.2所示。
这样可由④⑤求得X射线穿过椭圆时的行进距离。
再乘上各个椭圆的线性衰减系数μ后累加起来就可得到X射线的投影值。
④⑤图1.1 图1.2表1.1Shepp-Logan头部模型表1.1 r的单位是角度，μ为负值时表示削弱原有椭圆的衰减系数利用滤波反投影重建算法，实现对Shep-Logan头模型的重建。
要用到的原理有：傅立叶切片定理、快速傅立叶变换FF T以及滤波函数的设计。
三、主要仪器设备及实验耗材：具有XP或2000系统，并装有MATLAB系统的PC机。
四、实验内容或步骤：%1．用MATLAB图像处理工具箱的phantom生成Shep-Logan头模型；P=phantom(256);imshow(P)%2．用MATLAB中的radon函数获得Shepp-Logan模型的投影数据；theta1=0:10:170;[R1,xp]=radon(P,theta1); %计算Shep-Logen头模型18个角度theta2=0:5:175;[R2,xp]=radon(P,theta2); % 36个角度theta3=0:2:178;[R3,xp]=radon(P,theta3); % 90个角度的投影%显示投影数据：%18个角度figure,imagesc(theta1,xp,R1);xlabel('\theta');ylabel('x\prime');% 36个角度figure,imagesc(theta2,xp,R2);xlabel('\theta');ylabel('x\prime');% 90个角度figure,imagesc(theta3,xp,R3);xlabel('\theta');ylabel('x\prime');%3．用MATLAB中的iradon函数对获得的投影数据进行滤波反投影重建，获得Shepp-Logan 模型的重建图像。
滤波反投影法：
滤波反投影法根据附件三所给接收信息，采用先修正、后投影重建图像的做法，可得到原始图像的吸收率信息。
其原理为：在得到某一角度下的投影函数（一维函数）后，对此函数做滤波处理，得一修正后的滤波函数，再对修正后的滤波函数做反投影运算，得待检测介质吸收率在正方形托盘中的每一点的分布密度函数。
图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命令实现。
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 命令实现。
第26 卷第11 期2006 年11 月光学学报ACTA OP TICA SIN ICAVo l. 26 ,No. 11November , 2006文章编号: 025322239 (2006) 112165729偏折层析的滤波反投影算法及误差分析宋张斌贺安之(南京理工大学信息物理与工程系, 南京210094)摘要: 对偏折层析投影转换为相位层析投影的转换关系迚行了分析,给出明晰的数学关系,幵针对偏折层析的滤波反投影算法重建的结果迚行误差分析。
分析结果表明投影噪声对重建场的作用体现在与由偏折层析滤波反投影算法的滤波器有关的倾斜函数上。
因此提出了改迚的偏折层析滤波反投影算法,数值模拟表明,改迚算法在有效抑制倾斜现象的同时,对重建结果不会造成明显的失真。
在此基础上改迚的算法被用于真实火箭燃气射流密度场的三维重建中。
关键词: 信息光学; 偏折层析; 重建算法; 误差分析中图分类号: O438 文献标识码: AFil t e red B ac k2P r oject i o n Al gori t h m of Def lect i on To m og r ap h ya n d Er r or A n al ys isSong Y a ng Zhang Bin He Anzhi( Dep a r t me n t of I nf or m a t i on Physics & Engi neeri ng Na n ji ng U n iversit y of Scie nce & Tech nology , Na n ji ng 210094)Abs t r act : The conversion f rom deflection tomography p rojection to phase tomograp hy p rojection is analyzed , and an explicit exp ression correspondin g to the conversion is p r esented. An er ror analysis is made to the reconst r ucted fields by f iltered back2p rojection ( DFB P ) algorithm of deflection tomography. Results show that the effect of p rojection noise on the reconst ructed fields is rep resented by a slope f unction related to the filter used in deflection tomograp hic f iltered back2p rojection algorithm. So the deflection tomographic f ilter ed back2p rojection algorithm is modified. Numerical simulation shows that the modified algorithm dep resses the slope ph enomena efficiently , while no obvious distortion is int roduced to the reconst r uction. Based on the modified algor ithm , the three2dimensional reconst ruction for den sity field of the real rocket exhausted plumes is carried out .Key w or ds : information optics ; deflection tomograp hy ; reconst ruction algorithm ; error analysis1 引言光学层析技术( Optical Comp uterized Tomograp hy ,O CT)是以光波为载体, 由加载了被测场信息的多方向投影数据重建待测场物理量分布的技术。
平行束滤波反投影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’经过滤波后，反投影后的图像。
基于滤波反投影算法的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变换的一个中心切片。
地理与生物信息学院2012 / 2013 学年第二学期实验报告课程名称：医学图像处理和成像技术实验名称：CT反投影滤波重建算法设计班级学号: B10090405学生姓名: 陈洁指导教师: 戴修斌日期：2013 年 5 月一、实验题目：CT反投影滤波重建算法设计二、实验内容：1.显示图像；2.获得仿真投影数据；3.基于获得的仿真投影数据重建图像。
三、实验要求：1.Shepp-Logan头模型：画出Shepp-Logan头模型，简称S-L模型，头模型尺寸设定为128×128；2.仿真投影数据的获得：从头模型中获得投影数据，投影数据格式为180×185，即[0，179°]范围内角度每隔1°取样，每个角度下有185个探测器；3.卷积反投影重建算法的实现：基于获得的仿真投影数据重建图像，使用R-L卷积函数，重建尺寸为128×128。
四、实验过程：实验1. Shepp-Logan头模型①算法实现流程：I. S-L头模型由10个位置、大小、方向、密度各异的椭圆组成，象征一个脑断层图像。
Shepp-Logan头模型中的椭圆参数：II. 使用循环语句给像素赋值：for i=1:10for x….for y…..判断点(x, y)是否在第i个椭圆内;如是,则将第i个椭圆折射指数赋给点(x, y);endendendIII. 显示仿真头模型：使用imshow(f,[])函数显示出图像。
②实验代码：clear all;p=[0 0 0.92 0.69 pi/2 10 -0.0184 0.874 0.6624 pi/2 20.22 0 0.31 0.11 72/180*pi 0-0.22 0 0.41 0.16 108/180*pi 40 0.35 0.25 0.21 pi/2 50 0.1 0.046 0.046 0 60 -0.1 0.046 0.046 0 7-0.08 -0.605 0.046 0.023 0 80 -0.605 0.023 0.023 0 80.06 -0.605 0.046 0.023 pi/2 8];N=256;x=linspace(-1,1,N);y=linspace(-1,1,N);f=zeros(N,N);for i=1:Nfor j=1:Nfor k=1:10A=p(k,3);B=p(k,4);x0=p(k,1);y0=p(k,2);x1=(x(i)-x0)*cos(p(k,5))+(y(j)-y0)*sin(p(k,5));y1=-(x(i)-x0)*sin(p(k,5))+(y(j)-y0)*cos(p(k,5));if((x1*x1)/(A*A)+(y1*y1)/(B*B)<=1) %判断条件f(i,j)=p(k,6);endendendendf=rot90(f);imshow(f,[])③运行结果：实验2. 获得仿真投影数据：①算法实现流程：I. θ∈ [00, 10, ..., 1790]， s ∈[-92, -91, ..., 91，92]；II. 对于第i 个椭圆求出对应θ和s 的仿真投影数据：其中，(x 0, y 0)为中心坐标，A 为长轴，B 为短轴，a 为旋转角度，ρ为折射指数。
《滤波反投影程序设计报告》课程名称：生物医学图像处理2院系：生物医学工程姓名：学号：完成日期： 2017年4月23日一、设计目的用Matlab实现平行束滤波反投影算法，比较不同滤波函数的效果。
二、实验原理（一）图像重建模型——shepp Logan头模型是图像重建标准体模，由10个位置、大小、方向、密度各异的椭圆组成，代表一个脑部断层。
（二）重建理论推导中心切片定理是从投影图像重建图像的理论基础，表述为：某断层图像f(x,y)在视角为θ时得到的平行投影的一维傅里叶变换等于f(x,y)二维傅里叶变换F（w1,w2）过原点的一个垂直切片，且切片与轴w1相交成θ角。
如下图所示。
公式表述为：F（wcos,wsin）=P(w,) ①将在-坐标系中表达的F（w1,w2）引入新的极坐标系中，新坐标系与原坐标系原点重合，有w1=wcos,w2=wsin.面元换算为dw1dw2=wdwd.有 f(x,y)=== +②注意到在其傅里叶变换存在如下关系：P(w,将上式代入②式，有f(x,y)=③令③式内积分等于g(xcos+ysin),则有g(xcos+ysin)=t=xcosθ+ysinθdw实际上，g(xcos+ysin)就是投射角度为时的滤波投影，在t-s坐标系中表达时，转化为g（t,）=h(t)*p(t,),h(t)是传递函数H（w）=|w|的傅里叶逆变换，p(t,)是P（w,）的傅里叶逆变换。
所以③式可写成f(x,y)=θ④在图中注意到Xr=rcos()=xcos是从原点出发的通过点（r,）的射线方程，④式可写为：f(x,y)=④⑤两式表明：f(x,y)在（x,y）处的重建，等于通过该点的所有角度下滤波投影的累加所得到的像素值，而Xr=rcos()=xcos的变化代表了所有平行投影射线。
（三）Radon变换一个无限薄的切片内相对线性衰减系数的分布是由它的所有线积分的集合唯一决定，揭示了函数和投影之间的关系，若函数为f(x,y),则不同角度下的投影可写为P(t,)=⑥（四）滤波函数由于直接反投影法把取自有限空间的投影均匀回抹到了射线所及的无限空间的各个像素上，使得原来像素值为0的点不为0，从而产生星状伪迹，滤波反投影算法用人为设计的一维滤波函数对所得投影数据进行卷积，而后进行反投影和累加时，由于正负抵消，可一定程度上消除星状伪迹。
滤波反投影法的实施步骤引言滤波反投影法（Filtered backprojection）是一种重建计算机断层成像（Computed Tomography，CT）中常用的算法。
它通过滤波和反投影的过程，将射线在被测物体中的吸收信息转化为图像。
本文将介绍滤波反投影法的实施步骤，以帮助读者了解该算法的实现过程。
步骤一：数据采集1.准备一个旋转的X射线源和一个旋转的探测器，旋转范围通常为180度或360度。
2.将被测物体置于射线源和探测器之间，以获取多个不同角度下的投影数据。
3.每个角度下，记录探测器上的射线强度，得到一系列投影数据。
步骤二：滤波1.将获取的投影数据进行滤波处理，以增强图像的边缘信息。
2.使用一维高斯滤波器或锥形滤波器对每个投影数据进行滤波。
3.滤波的目的是减少伪影和噪音，同时保留有用的图像信息。
步骤三：反投影1.对滤波后的投影数据进行反投影计算。
2.对每个投影数据，将其按照相应的角度反投影到一个平面上。
3.反投影的过程是将每个点的投影值分布到相应的位置上，形成二维空间中的像素点。
步骤四：重建1.对反投影得到的所有平面进行叠加，以得到三维空间中的体素值。
2.使用重建算法，如Feldkamp算法或迭代重建算法，将二维平面上的像素点转化为三维体素值。
3.通过重建算法，可以还原出被测物体的内部结构信息。
步骤五：图像处理和显示1.对重建得到的体素值进行图像处理，如去噪、增强对比度等。
2.使用图像显示工具，如图像处理软件或CT设备，将处理后的图像进行可视化显示。
3.检查显示的图像是否符合预期，如有必要可以进行后续的图像处理和调整。
结论滤波反投影法是一种常用的断层成像重建算法，通过数据采集、滤波、反投影、重建和图像处理等步骤，可以将射线在被测物体中的吸收信息转化为可视化的图像。
本文介绍了滤波反投影法的实施步骤，希望能对读者的研究和实践工作提供有益的参考。
基于滤波反投影算法的CT系统参数标定及成像研究摘要本文主要针对CT系统的参数标定以及成像进行研究，在大致介绍CT设备及其成像原理的基础上，对题目所给问题进行分析计算解答，并采用基于滤波反投影重建算法的逆Radon 变换等算法对数据进行处理分析，得出结果，经过误差分析后，与CT系统自身特点相结合，结果比较理想。
问题一:根据X射线与标定模板的特殊位置如图（4）和图（8），两位之间角度90o，分析特殊位置图形，然后用MATLAB软件可画出CT系统绕某固定的旋转中心逆时针旋转180次的图，找到与特殊位置图形描述相近的图形，发现两特殊位置旋转次数差为90次，由此确定X射线的180个方向，探测器单元之间的距离=椭圆的长轴长/长轴长投影在探测器上单元的数量，由于旋转中心的位置和椭圆中心的位置相对不变，所以根据x和y的偏移量来确CT系统旋转中心在正方形托盘的位置。
问题二：建立CT系统滤波反投影重建算法，用iradon函数对附件2接收信息数据逆变换，即可得到出未知介质的吸收率，未知介质的几何形状，位置等信息。
问题三：由题意可知，该题同问题二相似，只是要求求另一种未知介质的相关信息，所以此题所用数学方法和问题二相同。
问题四：题目所给标定模型存在一定的规律性，并且容易受外界环境的影响，因此计算得到的CT系统参数存在误差。
所以我们自行设计了新的模板来改进标定精度和稳定性。
关键词：特殊位置、滤波反投影重建算法、逆Radon变换一、问题重述CT可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。
一种典型的二维CT系统是平行入射的X射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。
X射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。
对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。
图像重建是CT技术的一个研究热点。
重建算法的现实是对算法研究的一个重要环节。
由Math Works公司推出的MATLAB工具软件具有强大的数学计算和图像处理功能，并为算法提供了
一个方便有效的研究和实现的平台。
本文在图像重建分析的基础上，运用MATLAB实现了基于扇
束的滤波反投影重建算法的计算机模拟。
引言
图像重建技术在诸多领域中发挥着重要作用，在重建算法的研究
和实现过程中，存在着是一系列极其复杂的图像处理问题和数学计算
图５ １２８×１２８的Ｓｈｅｅｐ
Ｌｏｇａｎ头模型图像
图６ 扇束射角增量为０．３°投影值
2. 选取滤波函数，并离散化处理，如：R-L滤波函数，则离散化形式为：
（式1-8）其中：
（式1-9）
图７ 重建效果图
总结
本文在分析基于扇束滤波反投影算法的基础上，详细介绍了该算法“模拟产生投影数据——修正投影——加权滤波——反投影重建”整个计算机现实过程，并充分利用Matlab强大的图像处理功能，无需大量的编程，现实了图像重建算法的计算机模拟。
高效的工程计算语言，它从本质上提供了对图像的支持，使用它可以对数字图像形成的离散数据矩阵进行一次性的处理，较其他高级语言。
图像滤波实验报告分析的三个基本流程下载温馨提示:该文档是我店铺精心编制而成，希望大家下载以后，能够帮助大家解决实际的问题。
文档下载后可定制随意修改，请根据实际需要进行相应的调整和使用，谢谢!并且，本店铺为大家提供各种各样类型的实用资料，如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等，如想了解不同资料格式和写法，敬请关注!Download tips: This document is carefully compiled by theeditor.I hope that after you download them,they can help yousolve practical problems. The document can be customized andmodified after downloading,please adjust and use it according toactual needs, thank you!In addition, our shop provides you with various types ofpractical materials,such as educational essays, diaryappreciation,sentence excerpts,ancient poems,classic articles,topic composition,work summary,word parsing,copy excerpts,other materials and so on,want to know different data formats andwriting methods,please pay attention!图像滤波实验报告分析：三大基本流程在图像处理领域，图像滤波是一种常见的技术，用于改善图像质量、消除噪声或突出特定特征。
《滤波反投影程序设计报告》课程名称：生物医学图像处理2院系：生物医学工程姓名：学号：完成日期：2017年4月23日一、设计目的用Matlab实现平行束滤波反投影算法，比较不同滤波函数的效果。
二、实验原理（一）图像重建模型——shepp Logan头模型是图像重建标准体模，由10个位置、大小、方向、密度各异的椭圆组成，代表一个脑部断层。
（二）重建理论推导中心切片定理是从投影图像重建图像的理论基础，表述为：某断层图像f(x,y)在视角为θ时得到的平行投影的一维傅里叶变换等于f(x,y)二维傅里叶变换F（w1,w2）过原点的一个垂直切片，且切片与轴w1相交成θ角。
如下图所示。
公式表述为：F （wcos θ,wsin θ）=P(w,θ) ① 将在ω1-ω2坐标系中表达的F （w 1,w 2）引入新的极坐标系ω−θ中，新坐标系与原坐标系原点重合，有w1=wcos θ,w2=wsin θ. 面元换算为dw1dw2=wdwd θ.有 f(x,y)= ∫∫F （wcosθ，wsinθ）e j2πw(xcosθ+ysinθ)wdwdθ∞02π0 =∫∫P(w,θ)e j2πw(xcosθ+ysinθ)wdwdθ∞02π0 =∫∫P(w,θ)e j2πw(xcosθ+ysinθ)wdwdθ∞0π0 +∫∫P(w,θ+π)e −j2πw(xcosθ+ysinθ)wdwdθ∞0π0 ② 注意到在θ与θ+π两个角度下的平行束射线投影的投影存在如下关系：p (t,θ)=p (−t,θ+π)其傅里叶变换存在如下关系：P(w,θ+π)=P(−w,θ)将上式代入②式，有 f(x,y)= ∫∫P (w,θ)|w|e j2πw(xcosθ+ysinθ)dwdθ∞−∞π0 ③ 令③式内积分等于g(xcos θ+ysin θ),则有g(xcos θ+ysin θ)=∫|w |P (w,θ)e j2πwt |+∞−∞t= xcosθ+ysinθdw实际上，g(xcos θ+ysin θ)就是投射角度为θ时的滤波投影，在t-s 坐标系中表达时，转化为g （t,θ）=h(t)*p(t,θ),h(t)是传递函数H （w ）=|w|的傅里叶逆变换，p(t, θ)是P （w, θ）的傅里叶逆变换。
所以③式可写成f(x,y)= ∫g(x π0cosθ+ysinθ)d θ ④ 在图3.5中注意到Xr=rcos(θ−φ)=xcos θ+ysinθ是从原点出发的通过点（r,θ）的射线方程，④式可写为：f(x,y)= ∫g[rcos (θ−φ),φ]π0dφ ⑤ ④⑤两式表明：f(x,y)在（x,y ）处的重建，等于通过该点的所有角度下滤波投影的累加所得到的像素值，而Xr=rcos(θ−φ)=xcos θ+ysinθ的变化代表了所有平行投影射线。
（三）Radon 变换一个无限薄的切片内相对线性衰减系数的分布是由它的所有线积分的集合唯一决定，揭示了函数和投影之间的关系，若函数为f(x,y),则不同角度下的投影可写为P(t,θ)=∬f (x,y )δ(xcosθ+ysinθ−t )dxdy +∞−∞ ⑥（四）滤波函数由于直接反投影法把取自有限空间的投影均匀回抹到了射线所及的无限空间的各个像素上，使得原来像素值为0的点不为0，从而产生星状伪迹，滤波反投影算法用人为设计的一维滤波函数对所得投影数据进行卷积，而后进行反投影和累加时，由于正负抵消，可一定程度上消除星状伪迹。
现在常用的滤波函数有R-L、S-L、Hanning、Hamming、Parzen、Sigwin窗函数，本次设计比较了R-L、S-L和Parzen 滤波函数的效果，发现Parzen滤波效果最好。
1.R-L滤波函数R-L滤波函数频率响应为：R-L滤波函数滤波计算简单，避免了大量的正弦、余弦计算，所重建图像轮廓清楚，空间分辨率高，但有Gibbs现象，表现为明显的振荡响应，特别是工件的边缘和衰减系数变化剧烈（即密度变化）时，尤为明显。
S-L滤波器对投影中的高频成分具有抑制作用，进而使重建图像的振荡响应减小，对含噪声的投影数据，重建质量较用RL的好，但在低频段不及R-l滤波函数好，这是由于H（w）在低频段也偏离了|w|的缘故3.Parzen滤波函数三、程序实现程序包含FBP_GUI.m和dps.m两个m文件，其中投影、滤波和反投影过程均为自己编写，没有使用Matlab自带函数，附录中对过程有详细注释，程序大致流程为:四、运行结果（1）sheep Logan256*256图像重建（2）自定义灰度图像重建（3）自定义RGB图像重建（滤波函数进行滤波处理。
五、总结本次程序设计完成了设计目的，投影和反投影过程没有使用自带函数，但是在处理像素比较大的图像时程运行时间比较长，是一个缺点，图像也存在一定误差，效果不是特别完美，不过也算是锻炼了自己的能力，对图像重建中对滤波反投影算法有了深刻的了解，增加了对这门课程的兴趣，期待在以后的学习中能进一步完善程序，缩短运行时间，减小图像误差。
六、附录——程序代码①dps.m文件代码如下：function [a]=dps(P)tic;%P=phantom(256);%P=imread('gray1.jpg');%P=rgb2gray(P);[N,N]=size(P);subplot(2,3,2);imshow(P);title([int2str(N),'*',int2str(N),'原始图像']);%先进行自定义radon变换------------------------------------------------------------thm=45; %45度时会出现最大尺寸pre = imrotate(P,thm);[mmax,nmax] = size(pre);s=1;%创建一个180*nmax的空白图片，用以存储投影后的线状图片Final = zeros(180/s,nmax);%这里180代表180角度，每个角度投影成为一条线t = 1;for theta = 1:s:179Protate = imrotate(P,theta); %对原图旋转一个角度，求和（线积分）Pf = sum(Protate,1);[mreal,nreal]=size(Pf); %计算实际尺寸%确定起始点if (nmax - nreal)/2-floor((nmax - nreal)/2) == 0From = floor((nmax - nreal)/2 + 1);%总点数为偶数时elseFrom = floor((nmax - nreal)/2) + 1;%总点数为奇数时end%确定结束点End = floor((nmax-nreal)/2) + nreal;%将一个角度Radon变换后的线状图存入结果图像的某一行Final(180/s-t,From:End) = Pf; %从最底下一行开始存起%上移一行t = t + 1;end%再逆时针旋转R=imrotate(Final,90);subplot(2,3,3);imshow(R,[]);title('自定义投影后图像');z=2*ceil(norm(size(P)-floor((size(P)-1)/2)-1))+3;% radon变换默认平移点数/角度e=floor((z-N)/2)+2;R1=R(e:(N+e-1),:);[mm,nn]=size(R1);subplot(2,3,4);imagesc(R1);title([int2str(mm),'*',int2str(nn),'正弦图']);colormap(gray);colorbar;%滤波函数构造------------------------------------------------------------g=1-N:N-1;% d=1;% R-L滤波函数% for i=1:2*N-1% if g(i)==0% h(i)=1/(4*(d^2));% else if mod(g(i),2)==0% h(i)=0;% else% h(i)=(-1)/(pi^2*d^2*(g(i)^2));% end% end% end%S-L滤波函数% d=1;% for i=1:2*N-1% h(i)=2/(pi^2*d^2*(1-4*g(i)^2));% end%Parzen滤波函数for i=1:2*N-1h(i)=(24*pi*g(i)*cos(pi*g(i))-96*sin(pi*g(i))-48*pi*g(i)*cos(pi*g(i)/2)+384*sin(pi*g(i)/2)-2*(pi^3)*(g(i)^3)-72*pi*g(i))/(4*(pi^5)*(g(i)^5));endh(N)=0.0435;%将投影与滤波函数卷积----------------------------------------------------G=zeros(N,180);for m=1:180for i=1:Nb=0;for n=1:Nb=b+h(N+n-i)*R1(n,m);G(i,m)=b;endendend%投影滤波后线性内插进行图像重建----------------------------------------------a=zeros(N); %重建图像初始化，每个像素点像素值为0ns=(N+1)/2;for m=1:180 %遍历每个投影角度r=pi/180; %将角度换为弧度for i=1:Nfor j=1:N %遍历原图的每一个像素点Xrm=(i-ns)*cos(m*r)+(j-ns)*sin(m*r)+ns-1; %坐标转换if Xrm<1n=1; %内插运算整数值t=abs(Xrm)-floor(abs(Xrm)); %内插运算小数值elsen=floor(Xrm);t=Xrm-floor(Xrm);endif n>N-1n=N-1;endc=(1-t)*G(n,m)+t*G(n+1,m); %内插后的滤波投影值a(N+1-j,i)=a(N+1-j,i)+c; %像素值的叠加endendend%将P、a的像素值变换到0-1之间P=double(P);P=P-min(min(P));kk=max(max(P));for i=1:Nfor j=1:NP(i,j)=P(i,j)/kk;endenda=double(a);a=a-min(min(a));kkk=max(max(a));for i=1:Nfor j=1:Na(i,j)=a(i,j)/kkk;endendsubplot(2,3,5);imshow(a,[]); %重建图形的显示title('滤波反投影重建图像');%重建图像质量评价--------------------------------------------------------%计算归一化均方距离判据d;ave=0;for x=1:Nfor y=1:Nave=ave+P(x,y);endendave=ave/(N*N);x1=0;x2=0;for x=1:Nfor y=1:Nx1=x1+(P(x,y)-a(x,y))^2;x2=x2+(P(x,y)-ave)^2;endendevaluate_d=sqrt(double(x1/x2));%计算归一化平均绝对距离判据r;x3=0;x4=0;for x=1:Nfor y=1:Nx3=x3+abs(P(x,y)-a(x,y));x4=x4+abs(P(x,y));endendevaluate_r=x3/x4;%计算最坏情况距离判据ens=floor(N/2)-1;g=zeros(ns);for x=1:nsfor y=1:nsT=(P(2*x,2*y)+P(2*x+1,2*y)+P(2*x,2*y+1)+P(2*x+1,2*y+1))/4;R=(a(2*x,2*y)+a(2*x+1,2*y)+a(2*x,2*y+1)+a(2*x+1,2*y+1))/4;g(x,y)=abs(T-R);endendevaluate_e=max(max(g));%结果文本显示------------------------------------------------------------o=ones(500,1000);subplot(2,3,6);imshow(o,[]);s_title=['图像重建精度判据如下:'];text(0,0,s_title,'Fontsize',14);s=num2str(toc);s_one=['run time = ' s ' s;'];text(0,100,s_one,'FontSize',10);s=num2str(evaluate_d);s_two=['归一化均方距离判据d=' s ';'];text(0,200,s_two,'Fontsize',10);s=num2str(evaluate_r);s_three=['归一化平均绝对距离判据r=' s ';'];text(0,300,s_three,'Fontsize',10);s=num2str(evaluate_e);s_four=['最坏情况距离判据e=' s ';'];text(0,400,s_four,'Fontsize',10);toc②FBP_GUI.m文件代码如下：function varargout = FBP_GUI(varargin)% FBP_GUI MATLAB code for FBP_GUI.fig% FBP_GUI, by itself, creates a new FBP_GUI or raises the existing% singleton*.%% H = FBP_GUI returns the handle to a new FBP_GUI or the handle to% the existing singleton*.%% FBP_GUI('CALLBACK',hObject,eventData,handles,...) calls the local% function named CALLBACK in FBP_GUI.M with the given input arguments. %% FBP_GUI('Property','Value',...) creates a new FBP_GUI or raises the% existing singleton*. Starting from the left, property value pairs are% applied to the GUI before FBP_GUI_OpeningFcn gets called. An% unrecognized property name or invalid value makes property application % stop. All inputs are passed to FBP_GUI_OpeningFcn via varargin.%% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one% instance to run (singleton)".%% See also: GUIDE, GUIDATA, GUIHANDLES% Edit the above text to modify the response to help FBP_GUI% Last Modified by GUIDE v2.5 15-Apr-2017 14:24:03% Begin initialization code - DO NOT EDITgui_Singleton = 1;gui_State = struct('gui_Name', mfilename, ...'gui_Singleton', gui_Singleton, ...'gui_OpeningFcn', @FBP_GUI_OpeningFcn, ...'gui_OutputFcn', @FBP_GUI_OutputFcn, ...'gui_LayoutFcn', [] , ...'gui_Callback', []);if nargin && ischar(varargin{1})gui_State.gui_Callback = str2func(varargin{1});endif nargout[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); elsegui_mainfcn(gui_State, varargin{:});end% End initialization code - DO NOT EDIT% --- Executes just before FBP_GUI is made visible.function FBP_GUI_OpeningFcn(hObject, eventdata, handles, varargin) % This function has no output args, see OutputFcn.% hObject handle to figure% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % varargin command line arguments to FBP_GUI (see VARARGIN)% Choose default command line output for FBP_GUIhandles.output = hObject;% Update handles structureguidata(hObject, handles);% UIWAIT makes FBP_GUI wait for user response (see UIRESUME)% uiwait(handles.figure1);% --- Outputs from this function are returned to the command line. function varargout = FBP_GUI_OutputFcn(hObject, eventdata, handles) % varargout cell array for returning output args (see VARARGOUT);% hObject handle to figure% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)% Get default command line output from handles structure varargout{1} = handles.output;% --- Executes on button press in pushbutton1.function pushbutton1_Callback(hObject, eventdata, handles)% hObject handle to pushbutton1 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)I=phantom(256);[A]=dps(I);% --- Executes on button press in pushbutton2.function pushbutton2_Callback(hObject, eventdata, handles)% hObject handle to pushbutton2 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) [fn,pn,fi]=uigetfile('*.*','选择图片');lujing=[pn fn];%得到待重建图像存储路径I=imread(lujing);[A]=dps(I);% --- Executes on button press in pushbutton3.function pushbutton3_Callback(hObject, eventdata, handles)% hObject handle to pushbutton3 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) [fn,pn,fi]=uigetfile('*.*','选择图片');lujing=[pn fn];%得到待重建图像存储路径I=imread(lujing);I=rgb2gray(I);[A]=dps(I);% --- Executes on button press in pushbutton4.function pushbutton4_Callback(hObject, eventdata, handles)% hObject handle to pushbutton4 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) close all。
