图像复原之由投影重建图像

简述部分

引言

• 最初接触由投影重建图像这块内容的时候是在车牌识别中。

• 上图是在90°的投影下的结果。
• 下面我们开一个简单图像的特定角度下的投影
  
• 当我们收集到各个角度的投影后，并希望通过这些投影的图像重建原图像。上图的最后一张图片就是直接重建的图像。可见，其有非常明显的“晕环”现象。

雷登变换

• 雷登变换阐述了一幅图像与其在各个角度下投影的具体表示。
  

• 下面我们在matlab上展示一副图像的雷登变换
  

• 图像在各个角度下的雷登变换的集合我们称之为正弦图。

• 以第一行的图像为例。明显的，其像素为白色在90°的投影下是最多的，在0°或180°的投影下是最少的。对应于右侧的正弦图也能够体现出来。

• 试想这样一个问题，如果我们把各个角度下的投影经过一次反投影在求和是否会复原图像呢？答案是肯定的。

*上图显示了由正弦图直接得到的反投影图像。如引言所述，可见其有非常明显的“晕环现象”。有人可能会想到，如果将每次投影的角度间隔选的小一点是否还会存在这样的问题呢？当然了，增大采样次数是一个消耗资源的方法。而我们这里还有更好的解决这个问题的办法，这个方法是建立在傅里叶切片定理上的。

傅里叶切片定理

• 傅里叶切片定理用一句话表示就是：一个投影的一维傅里叶变换就是得到该投影原图的二维傅里叶变换的一个切片，其切片角度就是投影的角度θ。

• 由此定理我们就可以在投影的频域做做文章来消除晕环现象了，这种重建方法称之为滤波反投影法。

• 总结由此方法得到反投影的步骤如下：
  

  1. 计算每个投影的一维傅里叶变换
     
  2. 用一个滤波函数|w|乘以每个傅里叶变换，就是加窗。
     
  3. 得到每个滤波后的一维反傅里叶变换。
     
  4. 将3得到的求和

• 这就是滤波反投影法（家汉明窗）复原的图像。可见，已经很好的消除了“晕环现象”。

扇形射线束滤波反投影的重建

• 试想，当我们用扇形射线束代替上述的平行射线束自然会有更加不错的效果。
  

• 扇形射线束是当前CT系统使用的方法，具有高分辨率，高SNR和更快的扫描时间。

---

• 下面用matlab的fanbeam实现基于扇形射线束的投影图像

• 使用ifanbeam实现图像重建
  

• 左图是直接重建的效果，右图加汉明窗并且将传感器间隔缩小到原来1/10的效果。

示例源码

clc;clear;close all;g1 = zeros(600, 600);g1(100:500, 250:350) = 1;g2 = phantom('Modified Shepp-Logan', 600);subplot(2,2,1);imshow(g1);subplot(2,2,3); imshow(g2);theta = 0:0.5:179.5;%执行雷登变换[R1, xp1] = radon(g1, theta);[R2, xp2] = radon(g2, theta);%显示投影图像(正弦图)r1_show = flipud(R1');r2_show = flipud(R2');subplot(2,2,2);imshow(r1_show, [], 'XData', xp1([1 end]), 'YData', [179.5, 0]);axis xy;axis on;xlabel('\rho');ylabel('\theta');subplot(2,2,4);imshow(r2_show, [], 'XData', xp2([1 end]), 'YData', [179.5, 0]);axis xy;axis on;xlabel('\rho');ylabel('\theta');%从正弦图得到反投影图f1 = iradon(R1, theta, 'Hamming');f2 = iradon(R2, theta, 'Hamming');figure;subplot(1,2,1);imshow(f1, []);subplot(1,2,2);imshow(f2, []);%}%使用扇形射线束D = 1.5*hypot(size(g1,1), size(g2,2))/2;b1_line=fanbeam(g1, D, 'FanSensorGeometry', 'line', 'FanSensorSpacing', 1, 'FanRotationIncrement', 0.5);b1_line_s=flipud(b1_line');b2_line=fanbeam(g2, D, 'FanSensorGeometry', 'line', 'FanSensorSpacing', 1, 'FanRotationIncrement', 0.5);b2_line_s=flipud(b2_line');b1_src=fanbeam(g1, D, 'FanSensorGeometry', 'arc', 'FanSensorSpacing', .08, 'FanRotationIncrement', 0.5);b2_src=fanbeam(g2, D, 'FanSensorGeometry', 'arc', 'FanSensorSpacing', .08, 'FanRotationIncrement', 0.5);figure;subplot(2,2,1);imshow(g1);subplot(2,2,2); imshow(b1_line_s, [], 'XData', [0, 850], 'YData', [0, 360]);axis xy;axis on;ylabel('扇形旋转角度');xlabel('传感器个数');subplot(2,2,3);imshow(g2);subplot(2,2,4); imshow(b2_line_s, [], 'XData', [0, 850], 'YData', [0, 360]);axis xy;axis on;ylabel('扇形旋转角度');xlabel('传感器个数');%扇形反投影滤波重建B1 = fanbeam(g2, D);fB1=ifanbeam(B1, D);B2 = fanbeam(g2, D, 'FanSensorSpacing', .05, 'FanRotationIncrement', .5);fB2=ifanbeam(B2, D, 'FanSensorSpacing', .05, 'FanRotationIncrement', .5, 'filter', 'Hamming');figure;subplot(1,2,1);imshow(fB1, []);subplot(1,2,2);imshow(fB2, []);