【图像处理】平行线投影radon变换

　　从一个角度，用光源照射对象物体，屏幕上会形成对象物体的影子；如果物体是半透明的，那么影子便有灰度而不是纯黑的，这说明屏幕上的像可以反映物体内部对可见光的衰减作用。我们从落于[0~$\pi$]的一系列连续角度照射物体，形成一系列的像，这些像包含物体结构特征信息，基本上可以通过这些像还原物体的形状特征，如果物体是半透明的，那么物体内部的结构也可以还原出来。物体原始形状变换生成这些投影像，称为radon变换；从这些像还原物体形态，称为逆radon变换。人体对可见光是不透明的，但对X光是半透明的，因此CT可以发射X光照射人体，生成人体内部结构的图像信息。
　　radon变换的公式是：

$$x\cos{\theta}+y\sin{\theta}=\rho \\ g(\rho,\theta)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}{f(x,y)\delta(x\cos{\theta}+y\sin{\theta}-\rho)dxdy}$$

本质上就是沿着$x\cos{\theta}+y\sin{\theta}=\rho$确定的多条平行射线，建立法线方向为$\theta$的线段上的投影。一个比较简单的实现是，让图像均匀旋转$\theta$角度，然后计算x轴上的投影。假设角度区间[0~$\pi$]被分割为n份，投影线段长m，则最终我们获取$m\times n$大小的二维矩阵，称为投影矩阵。

　　radon反变换的公式是：

$$f(x,y)=\int_0^{\pi}{g(x\cos{\theta}+y\sin{\theta}, \theta)d\theta}$$ 该反变换操作比较简单，但是计算量大，且输出图像模糊有光晕。

　　为了得到清晰的图像，我们需要进行频域滤波。首先我们知道二维傅里叶变换对为：

$$F(u,v)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}{f(x,y)e^{-j2\pi(ux+vy)}dxdy} \\ f(x,y)= \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}{F(u,v)e^{j2\pi(ux+vy)}dudv}$$
这里引入傅里叶切片定理，其中 $\omega$是频率分量： $$G(\omega,\theta)=\int_{-\infty}^{\infty}{g(\rho,\theta)e^{-j2\pi \omega \rho}d\rho} \\ =\int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}{f(x,y)\delta(x\cos{\theta}+y\sin{\theta}-\rho)e^{-j2\pi \omega \rho}dxdyd\rho}} \\ = \int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}f(x,y)\left[ \int_{-\infty}^{\infty}{\delta(x\cos{\theta}+y\sin{\theta}-\rho)e^{-j2\pi \omega \rho}d\rho} \right]dxdy} \\ = \int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}f(x,y)e^{-j2\pi\omega(x\cos{\theta}+y\sin{\theta})}dxdy}\\ =F(\omega \cos{\theta}, \omega \sin{\theta})$$

这说明一个投影的一维傅里叶变换，是二维投影矩阵的二维傅里叶变换的一个切片。上述最后一步执行换元操作。
频域逆变换为：

$$f(x,y)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}{F(\omega \cos{\theta, \omega \sin{\theta}})d\omega\cos{\theta}d\omega\sin{\theta}} \\ \because d\omega\cos{\theta}d\omega\sin{\theta}=\omega d\omega d\theta \\ \therefore f(x,y)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}{F(\omega \cos{\theta, \omega \sin{\theta}})\omega d \omega d \theta} \\ =\int_{0}^{2\pi}\int_{0}^{\infty}{G(\omega,\theta)e^{j2\pi \omega(x\cos{\theta}+y\sin{\theta})}\omega d \omega d \theta} \\ \because \cos{(\theta+\pi）=-\cos{\theta}},\sin{(\theta+\pi）=-\sin{\theta}} \\ \therefore G(\omega,\theta+\pi)=\int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}f(x,y)e^{-j2\pi\omega(-x\cos{\theta}-y\sin{\theta})}dxdy}\\ =\int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}f(x,y)e^{-j2\pi(-\omega)(x\cos{\theta}+y\sin{\theta})}dxdy}=G(-\omega,\theta) \\ \rightarrow f(x,y)=\int_{0}^{\pi}\int_{0}^{\infty}{G(\omega,\theta)e^{j2\pi \omega(x\cos{\theta}+y\sin{\theta})}\omega d \omega d \theta} \\ + \int_{0}^{\pi}\int_{0}^{\infty}{G(-\omega,\theta)e^{j2\pi (-\omega)(x\cos{\theta}+y\sin{\theta})}\omega d \omega d \theta}\\let (-\omega\rightarrow t) \\ \rightarrow f(x,y)=\int_{0}^{\pi}\int_{0}^{\infty}{G(\omega,\theta)e^{j2\pi \omega(x\cos{\theta}+y\sin{\theta})}\omega d \omega d \theta} \\ + \int_{0}^{\pi}\int_{0}^{-\infty}{G(t,\theta)e^{j2\pi t(x\cos{\theta}+y\sin{\theta})}t d t d \theta} \\ =\int_{0}^{\pi}\int_{-\infty}^{\infty}{G(\omega,\theta)e^{j2\pi t(x\cos{\theta}+y\sin{\theta})} |\omega| d \omega d \theta}$$
我们可以看出上述式子推导过程中利用了[0~ $2\pi$]的数据，虽然最后只需要[0~ $\pi$]的投影数据，但是最终反滤波结果多了一项 $|\omega|$。这里面有什么更深层次的数学原理，我确实是不知道的。相比直接空间积分的结果，傅里叶切片重建效果更好的原因，可能有：

• 我们之前进行radon投影时，看到[0~$\pi$]和[$\pi$~ $2\pi$]的$\theta$角度下得到的投影是互为镜像，粗看认为信息是重复的，于是贸然舍弃一半，只计算了[0~$\pi$]角度的投影数据，这可能隐性地造成了数据不全。
• 直接进行空间积分变换，可能造成不同角度投影之间发生混淆，导致产生模糊和光晕

在最终式子中，积分计算$|\omega|$项是不可能的，该项不可积。我们可以通过引入窗函数，截断$|\omega|$项或者使用其他近似的窗函数，计算积分计算式并滤波，从而得到一个相对较好的结果。

滤波反投影重建逆radon变化步骤如下：

• 计算每个投影的一维傅里叶变换
• 使用截断的$|\omega|$项或者类似的窗口函数进行滤波，得到新的一维傅里叶变换数据
• 计算傅里叶反变换，获得原图

简单滤波的结果

Matlab重建结果，注意有条纹

Matlab的radon逆变换效果好，Matlab的牛逼只有研究过的人才懂，自己实现不了的痛也只有研究过的人才懂。。