积分线性变换可以将信号或图像在更适合的域内表达，并且使得解决相关问题更容易，在图像分析中最常用的积分显示变换是傅里叶变换、离散余弦变换与小波变换。
 
1d傅里叶变换由傅里叶(Fourier)提出，1d傅里叶变换将函数f(x)变换到频率域F(t)表达。F称作频谱，可以显示不同频率的相对成分。用i(根号下-1)表示虚数单位，1d连续傅里叶变换表达为：
 
 
逆变换为
 
 
相应地，对于下标从0开始的离散变量f(x)，1d离散傅里叶变换和逆变换分别为：
 
 
离散傅里叶变换(DFT)的时间复杂度为O(n2)，使用快速傅里叶变换(FFT)会降低为 O(nlogn)
 
从1d傅里叶变换可容易推广到2d，对于连续空间f，2d连续傅里叶变换和逆变换如下：
 
 
(x,y)表示图像坐标，(i,j)表示空间频率。
 
离散傅里叶变换和逆变换如下：
 
 
在matlab中，可使用fft2进行2d快速变换(1d使用fft)，如果用基本公式实现，复杂度为O（n4）。为了验证效果一下裁剪了局部：
 
I=imread('lena.jpg');
I=I(1:50,1:50);
[m,n]=size(I);
F=zeros(m,n);
%% FT
for u=1:m
    for v=1:n
        for x=1:m
            for y=1:n
                F(u,v)=F(u,v)+double(I(x,y))*exp(-2*pi*i*(u*x/m+v*y/n));
            end
        end
    end
end
I2=zeros(m,n);
for x=1:m
    for y=1:n
        for u=1:m
            for v=1:n
                I2(x,y)=I2(x,y)+double(F(u,v))*exp(2*pi*i*(u*x/m+v*y/n))/(m*n);
            end
        end
    end
end
I2=uint8(I2);
subplot(121),imshow(I);
subplot(122),imshow(I2);
F2=fft2(I);%直接调用
figure,
F=abs(F);
F2=abs(F2);
subplot(121),imshow(F,[]);
subplot(122),imshow(F2,[]);
F=fftshift(F);
F2=fftshift(F2);
F=log(F+1);
F2=log(F2+1);
figure,
subplot(221),imshow(F,[]);
subplot(222),imshow(F2,[]);
subplot(223),mesh(F);
subplot(224),mesh(F2);
 
其中，I为原图像，F为I的频谱，F2为直接fft2的频率，I2为F的逆变换结果(也可通过ifft2得到)，i为系统默认设置的常数跟下-1。直接运算或fft得到的F（如figure2）的频率高点都分布在四个角，直接使用fftshift可以把四个角的高频信息移动到最中间(imshow(I,[])在显示图像时会自动调整归一化，让图像正确显示，或者手动归一化实现后inshow(I))。之后abs将负实数和虚数部分调整为正实数，由于高频低频差距过大，所以用log(F+1)缓和高低。运行结果分别如下：
 
 
 
 
 figure1为裁剪后的图像与进行变换和逆变换后的图像，figure2为通过代码实现和直接调用库得到的未调整前的频谱(左图和右图角度不同，但内容一样)，figure3为调整后的频谱以及他们的mesh面。
 
将图像由空间域转为频率域在图像压缩以及去噪等方面有很好的应用。
 
之后文章：DCT WT
 
 

 

 

  
 
 
在上一部分当中，得到了利用三角函数表示周期函数的方法，但是对于非周期函数就...凉了。所以有什么办法吗？没办法（划掉）。这时候我们就需要拿出来我们的黑科技——傅里叶变换。
 
一、傅里叶级数的推广
 
当然这东西肯定不是凭空脑补出来的，而是将傅里叶级数进一步推广到非周期函数上。现在已经得到了周期函数的情况，一种很自然的想法就是将非周期函数化归到周期函数上，那么就可以继续套用傅里叶级数了。
 
如果要强行描述非周期函数的周期性，那它的周期就应该是无穷大，整个定义域都在它的一个周期内，以至于它不可能再重复这一周期。
 
把这个想法用形式化的语言表示出来，就是 
 

   的周期 
 
 

   。因为 
 
 

   ，那么 
 
 

   。接下来观察一下此时的傅里叶级数 
 
 

   。不大容易观(xuan)察(xue)，三角形式有点复杂，不如采用指数形式 
 
 

   。
 
 
当 
 

   时， 
 
 

   从原本的离散变化变成了连续变化， 
 
 

   也就可以表示为关于 
 
 

   的函数 
 
 

   。
 
 
傅里叶级数中 
 

   ，事实上，这个积分的上下限不一定是 
 
 

   到 
 
 

   ，只需要积 
 
 

   的一个周期就可以了。
 
 
换句话说，对于任意的 
 

   ，系数可以表示为
 
 

 

   。
 
 
这个积分需要积一整个周期，而此时的周期为无穷大，也就是整个定义域上都需要积，所以要从 
 

   积到 
 
 

   。
 
 
只需要让上式中的 
 

   ， 
 
 

   ，便可以得到 
 
 

   的表达式。不妨令 
 
 

   ，就得到了 
 
 

   。
 
 
因为 
 

   ，所以 
 
 

   ？当然不是，右侧的积分可能为无穷大，无穷小与无穷大的积不一定为无穷小。（如果等于零的话岂不是很有毒）
 
 
但是这对无穷大和无穷小的阶并不好比较，我们得不出 
 

   究竟应该等于什么值。既然 
 
 

   这么烦，那不如把它从这里面丢出去，之后用到 
 
 

   的时候再乘回来就好了，
 
 

 

   。
 
 
现在有了傅里叶级数对应的系数，该搞一搞 
 

   这个式子了。把对应的系数 
 
 

   代进去，再代入 
 
 

   ，变形后有
 
 

 

  
 
 
因为 
 

   ，每次 
 
 

   的增量 
 
 

   都是由于 
 
 

   变为 
 
 

   造成的，所以 
 
 

 

   。
 
 
同时 
 

   连续变化，原本的离散意义下的求和就该变为连续意义下的积分，搞出来 
 
 

 

   。
 
 
至此便推导出了傅里叶变换的两个公式
 

 

  
 
 

 

  
 
 
上式称为傅里叶变换，下式称为傅里叶逆变换。
 
还有另一个版本的傅里叶变换是
 

 

  
 
 

 

  
 
 
这两个版本都差不多，不过就是 
 

   这个系数的处理方法不大一样。mathematica上采用的是第二个版本的傅里叶变换，之前算了半天都跟自己手算的不一样，还以为自己算错了（溜
 
 
二、傅里叶变换的条件
 
由于傅里叶变换是从傅里叶级数推导得来的，所以还是狄利克雷条件，不过此时还要加上第三条， 
 

   在一个周期内绝对可积。
 
 
这一个条件在 
 

   为周期函数时，可以由前两个条件推出来，因为周期和函数值均为有限值，所以在一个周期内一定绝对可积。但是推广到傅里叶变换后，这个推导就不成立了，需要单独判定第三个条件。
 
 
三、性质
 
以下均默认 
 

   表示可以进行傅里叶变换的函数， 
 
 

   ， 
 
 

  
 
 
函数的卷积 
 

  
 
 
1、线性性质： 
 

   , 
 
 

  
 
 
2、尺度变化： 
 

   ， 
 
 

  
 
 
3、对称性： 
 

  
 
 
4、时移性： 
 

   ，
 
 

  
 
 
5、频移性： 
 

   ， 
 
 

  
 
 
6、时域卷积定理： 
 

  
 
 
7、频域卷积定理： 
 

  
 
 
8、微分运算： 
 

  
 
 
这些运算性质都是在采取第一种形式的傅里叶变换下的性质，如果使用第二种形式，会在某些性质上带来常数因子上的差别。
 
前面的7种运算性质的证明用积分的性质，再做点变量代换乱搞搞就可以了。这里主要说下微分运算性质的证明，用分部积分。只用证一阶导的情况就可以了，证出来之后使用数学归纳法可以很容易地推广到任意阶导数的情况。
 

 

  
 
 
微分运算的性质使得傅里叶变换能够将复杂的微分运算转化为简单的乘法运算，所以这个性质的常见应用在于解微分方程。通过傅里叶变换使微分方程变为代数方程，解出代数方程后再利用傅里叶逆变换求出原微分方程的解。
 
举个栗子，解物理上的简谐振动方程，除了常用的特征根法，还能够使用傅里叶变换
 

 

  
 
 
令 
 

   ，方程两边同时傅里叶变换
 
 

 

  
 
 
定义
 

 

   ，
 
 
使得 
 

   。
 
 
解出 
 

   ， 
 
 

   为常数
 
 
进行逆变换
 

 

  
 
 
使用辅助角公式合并
 

 

   ， 
 
 

   , 
 
 

   为常数
 
 
傅里叶变换在微分方程上的应用不局限于此，还能够应用于偏微分方程。但是最常用的并不是傅里叶变换，而是它的一般形式拉普拉斯变换
 

 

  
 
 
四、广义傅里叶变换 
 
在实际问题当中，经常会遇到一些函数并不满足绝对可积的条件，因而它们对应的傅里叶变换积分发散，并不存在傅里叶变换。但是我们又需要它们的傅里叶变换，所以就有了广义上的傅里叶变换。
 
比如刚刚求的简谐振动方程，对应的代数方程解出来后，发现 
 

   是发散的，此时我们通过定义了一个新函数 
 
 

   解决了发散的问题。暂时无视掉函数发散的问题，带着无穷大继续运算，最后逆变换时再作处理，这便是广义傅里叶变换的核心思想。
 
 
考虑正余弦函数，它们严格意义上的傅里叶变换都是不存在的，但是可以表示为
 

 

  
 
 

 

  
 
 
五、几何意义
 
傅里叶变换的几何意义类似傅里叶级数，当 
 

   时，所有的三角函数 
 
 

   ， 
 
 

   ， 
 
 

   ， 
 
 

   ( 
 
 

   )两两正交。换句话说，所有的三角函数都作为基向量，将 
 
 

   向它们投影。
 
 
实际上，无论是傅里叶级数还是傅里叶变换，都是在无穷维的希尔伯特空间中，将函数定义为空间中的向量，通过三角函数这样一组基向量表示空间中的任意函数。
 
六、物理意义
 
emm这一部分跟数学和oi的关系都不是特别大，就大概简略的写一下了，详细的介绍在网上也有很多资料，详细写的话怕是能再写这么长一篇文章（我懒）。
 
傅里叶级数将函数分解到离散的频率之上，而傅里叶变换将函数分解到连续的频域中，这样使原本频域上离散的点变成一条连续的曲线，对应的就是 
 

   的图像。 
 
 

   描述的是 
 
 

   在 
 
 

   这个频率分量上的大小。
 
 
基于这样的物理意义，傅里叶变换在实际问题当中得到大量应用。比如说最常见的是音乐软件上那个疯狂抖动的条，我也不知道这东西叫啥，反正就是下面这个图里进度条上面的那一坨。这个东西实际上是把现在正在播放的音频进行傅里叶变换，画出的频域图。
 

 

  
 
 
还有一种应用是视频以及图片的防伪和防盗版鉴别当中。将画面进行二维傅里叶变换，叠加高频分量，再进行逆变换即可。高频分量带来的差异很小，肉眼难以分辨，而且难以通过简单的截图和p图操作消除高频分量，因而是一种十分有效的“水印”。
 
除此以外，音视频的压缩也可以采用傅里叶变换，只保留强度较高的频率，去除较弱的频率，减少存储的数据量。

 
 
 
   
   
 
     
   
 
   
 
    

      三课时精通matlab傅里叶变换及逆变换和快速傅里叶变换 
    
 
    

      图像和算法等领域有多年研究和项目经验；指导发表科技核心期刊经验丰富；多次指导数学建模爱好者参赛。 
    
 
    
 
      
     宋星星 
    
 
   
 
   
 
    
 
     ¥12.00 
    
 
    
 
     立即订阅 
    
 
   
  
 
 
 
扫码下载「CSDN程序员学院APP」，1000+技术好课免费看
 
 
APP订阅课程，领取优惠，最少立减5元 ↓↓↓
 
  
 
   订阅后：请点击此处观看视频课程   
 
 
 
 
 
 视频教程-三课时精通matlab傅里叶变换及逆变换和快速傅里叶变换-Matlab 
 
 
 学习有效期：永久观看 
 
 
 学习时长：22分钟 
 
 
 学习计划：1天 
 
 
 难度：低 
 
 
   
 
  
 
  「口碑讲师带队学习，让你的问题不过夜」  
 
  
 
  讲师姓名：宋星星  
 
 
 工程师 
 
 
 讲师介绍：图像和算法等领域有多年研究和项目经验；指导发表科技核心期刊经验丰富；多次指导数学建模爱好者参赛。 
 
 
   ☛点击立即跟老师学习☚   
 
 
   
 
  
 
  「你将学到什么？」  
 
 
 精通matlab傅里叶变换及逆变换和快速傅里叶变换 
 
 
   
 
  
 
  「课程学习目录」  
 
 
 
  
 1.精通matlab傅里叶变换 
 2.精通matlab傅里叶逆变换 
 3.精通matlab快速傅里叶变换 
 
 
 
 
   
 
  
 
  「7项超值权益，保障学习质量」  
 
 
大咖讲解
 
 
技术专家系统讲解传授编程思路与实战。
 
 
答疑服务
 
 
专属社群随时沟通与讲师答疑，扫清学习障碍，自学编程不再难。
 
 
课程资料+课件
 
 
超实用资料，覆盖核心知识，关键编程技能，方便练习巩固。（部分讲师考虑到版权问题，暂未上传附件，敬请谅解）
 
 
常用开发实战
 
 
企业常见开发实战案例，带你掌握Python在工作中的不同运用场景。
 
 
大牛技术大会视频
 
 
2019Python开发者大会视频免费观看，送你一个近距离感受互联网大佬的机会。
 
 
APP+PC随时随地学习
 
 
满足不同场景，开发编程语言系统学习需求，不受空间、地域限制。
 
 
   
 
  
 
  「什么样的技术人适合学习？」  
 
 
想进入互联网技术行业，但是面对多门编程语言不知如何选择，0基础的你
掌握开发、编程技术单一、冷门，迫切希望能够转型的你
想进入大厂，但是编程经验不够丰富，没有竞争力，程序员找工作难。
 
 
   
 
  
 
  「悉心打造精品好课，1天学到大牛3年项目经验」  
 
 
 【完善的技术体系】 
 
 
技术成长循序渐进，帮助用户轻松掌握
 
 
掌握Matlab知识，扎实编码能力
 
 
 【清晰的课程脉络】 
 
 
浓缩大牛多年经验，全方位构建出系统化的技术知识脉络，同时注重实战操作。
 
 
 【仿佛在大厂实习般的课程设计】 
 
 
课程内容全面提升技术能力，系统学习大厂技术方法论，可复用在日后工作中。
 
 
   
 
  
 
  「你可以收获什么？」  
 
 
 熟练matlab傅里叶变换及逆变换和快速傅里叶变换 
 
 
 掌握matlab傅里叶变换及逆变换和快速傅里叶变换 
 
 
 
 
 
 
 
 
   
 
 
 
    
   

    订阅课程 开始学习 
   
  
 
 

 
 
绪
 
      在做CT图像处理的时候遇到很多问题，对于滤波反变换有许多细节存在疑问，经过多天查找资料和利用MATLAB程序一步步实现后终于豁然开朗，于是想要总结成文，作为笔记方便今后查看。文中若有错误欢迎指出！
 
 
 
目 录
 
（1）傅里叶逆变换法
 
（2）直接反投影法
 
（3）滤波反投影法
 
（4）MATLAB程序
 
 
 

      CT图像的形成和重构，在数学上的描述分别为拉东变换和拉东反变换。简单点说，拉动变换是用于描述通过X射线扫描物体形成CT图像这个过程的一种数学表示，而拉东反变换描述的则是对CT投影数据进行重构，还原成物体图像的一种数学方法。本文章主要讲的是CT图像重构的方法，也就是拉东反变换的具体实现方法，当然可能也有小伙伴想了解一下什么是拉东变换，这里贴上传送门：拉东变换。
 
其中拉东反变换在数学上的表达式如下：
 
 
由公式可看出拉东逆变换在数学上实际上是一个二维傅里叶逆变换，但由于二维傅里叶逆变换计算量大，耗时长，因此计算机对于拉东逆变换的具体实现方法主要为反投影法，反投影法又分为直接反投影法和滤波反投影法。由此可见CT重构方法主要为以下三种：
 
（1）傅里叶逆变换法
 
（2）直接反投影法
 
（3）滤波反投影法
 
下面将具体对这三种图像重构方法进行解释。
 
 
 

（1）傅里叶逆变换法
 
该方法基于一个重要的定理：中心切片定理。
 
该定理简单地理解就是：通过角度为θ扫描得到的投影，该投影的一维傅里叶变换，与对整个图像二维傅里叶变换后，二维频域中对应θ角度的一个切片信号是相同的，下面两个图理解起来更直观。
 
 
 

  
  
 
 
 
 

  
  
 
 
根据该理论，傅里叶逆变换法可以简单分成以下步骤：
 
① 假设每旋转1°就扫描一次，当对物体扫描了180°之后，我们就能得到180个投影信号（就是180根投影线）→在临床上，若使用平行扫描CT，我们拿到手的数据就是这个（在数学上，就是对图像进行拉东变换）
 
② 对180个投影信号进行一维傅里叶变换
 
③ 对②得到的180个一维频域信号，根据相对应的扫描角度，在空间中旋转排列，拼成一个二维频域空间（如图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。
 
      
 
 
 

  
  
 
 
 
 
这种先反投影，后滤波的方法需要用到两次二维傅里叶变换，计算量太大，因此也不是一个特别理想的矫正方法。
 
而滤波反投影法，顾名思义就是先滤波，后反投影的重构方法。具体步骤如下：
 
①对180个投影信号先分别进行一维傅里叶变换
 
②在频域中对所有投影信号进行滤波，也就是乘上权重因子r
 
③把所有滤波后的信号再进行一维傅里叶逆变换，还原到时域
 
④对每一个已经滤波的投影信号进行反投影，最后叠加（这一步跟前面的直接反投影法步骤完全相同，只是投影信号已经过了滤波）
 
具体流程见下图：
 
 
 

  
  
 
 
这个方法的好处是，把两次二维傅里叶变换变成了两次一维傅里叶变换，计算速度大大提升。于是滤波反投影法的核心问题就变成了如何选择一个合适的滤波器r。数学中的滤波器r是一个理想化的滤波函数，现实中不存在，因此只能设计一个近似的滤波函数来代替r的作用。
 
下面常用的滤波函数有：R-L滤波函数和S-L滤波函数。
 
①R-L滤波函数（Ramp-Lak滤波器）
 
 
 
从频域图中可以看到该滤波器的作用是把低频信号减少，从而突出高频信号，因此经过该函数滤波后，重建图像的轮廓会更清晰，并且函数简单，实现起来更方便。但由于该函数在频域上有一个加窗处理（就是把高频部分直接截断了），因此在时域中会出现有许多振荡的小尾巴（截窗振荡效应），这将会让重构出来的图像存在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现象。
 
 
 
（4）MATLAB程序
 
 
 
%% CT反投影法
%%==========================Part 1==========================
%% 1、直接反投影法
% 绘制一张简单的图，进行直接反投影法
clc,clear,close all
% ================================================
% 简单图
% I=zeros(512,512);
% for i=1:512
%     for j=1: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,[]);
%============================================

% 绘制代表角度的反投影图像
R=radon(I,0:179);  %对物体CT扫描180°的数据，每1°扫描一次
figure,imshow(R,[]);title('CT图像');
% 绘制1、45、90、135、180的投影信号和反投影图像
theta=[1 45 90 135];
for i=1:length(theta)
    r=R(:,i);
    II=iradon([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
II2=II1{1}+II1{2}+II1{3}+II1{4};
figure,imshow(II2,[]);

% 对反投影图像进行叠加
r=R(:,1);
II_r=iradon([r r],[1 1])/2;
k=[30 15 10 1];
for j=1:4
    A=II_r;
    for i=2:k(j):180
        r=R(:,i);
        II=iradon([r r],[i i],'linear','none')/2;%回抹（反投影）之前不滤波
        A=A+II;   
    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,[]);

% 对投影做傅里叶变换
R=radon(I,0:179);
width=length(R);
% 设计R-L滤波器
filter=2*[0:round(width/2-1), width/2:-1:1]'/width;
% 展示滤波前的CT投影信号
figure,imshow(R,[]),title('滤波前的CT投影信号');

% 每一列做傅里叶变换
r_fft=fft(R,729); 
% 每一列做傅里叶变换后滤波
r_fft_filter=r_fft.*filter; 
% 滤波后反变换(real取实部)
r_fft_filter_v=real(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 i=1:length(theta)
    r=r_fft_filter_v(:,i);
    II=iradon([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
II2=II1{1}+II1{2}+II1{3}+II1{4};
figure,imshow(II2,[]);

% 对反投影图像进行叠加
r=r_fft_filter_v(:,1);
II_r=iradon([r r],[1 1])/2;
k=[30 15 10 1];
for j=1:4
    A=II_r;
    for i=2:k(j):180
        r=r_fft_filter_v(:,i);
        II=iradon([r r],[i i],'linear','none')/2;
        %由于iradon默认有滤波，因此关掉默认的滤波选项（none），直接用上面已经手动滤波的数据进行反投影
        A=A+II;   
    end
%     figure,imshow(A,[min(min(I)) max(max(I))]);
    figure,imshow((A),[]);
end

% 上面是反投影叠加的分步计算，下面直接输入全部投影数据到iradon函数中进行验证
B=iradon(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有专门的限制时，该函数就会自动计算出这个“最长长度”，为了涵盖图像形状的不同情况（如有的图像为长方形），它设置了一个公用的计算公式：
 
N=2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3
 
【其中ceil函数表示“向上取整”；floor表示“向下取整”；norm表示“范数计算”，用在这里其实就是计算正方形斜边的长度】
 
如果套入这个公式，计算出来的N为：
 
size(I)=[512 512];
N=2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3

N =
   
   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函数都不专门设置时，会出现重构图像与原图像尺寸大小不一致的情况，但这并不会改变重构图像本身的缩放比例，因此若没有专门的尺寸大小必须前后一致的严格要求，可以不用专门设置。