Radon 变换
这周薛老师讲CT的时候说了一下CT重构的Radon变换，傅立叶变换一下然后按照方向填到相应位置再反傅立叶变换一下就能得到原图，感觉很神奇，当时没听懂，现在来看看到底是怎么个玩法
问题定义
问题定义如上图所示（纯个人手工绘制），中间圆环是我们想要求的物体，左边是发射源能够发射一种射线，右边为感应器能够感应到射线的强弱，假设物体上每一点对射线的衰减程度为函数 f(x,y) ，根据射线穿过物体的衰减程度可以得到物体在该方向上每一点的衰弱强度（intensity），如上图所示，相当于每个方向上都可以测量出该物体的“厚度”，我们的目的就是，根据不同方向上的物体的“厚度”求出物体上每一点的衰减程度
如上图所示，转换为数学表达可以表示为一种线积分，直线 L 穿过 f(x,y) ，所对应的强度（intensity）就是函数 f(x,y) 在直线 L 上的线积分 R_L = \int_L f(x,y)\mathrm{d}s \qquad R_L: \text{Intensity at $L$} \\我们的目标是根据不同的线 L 求出 f(x,y) 的表达式，但定义里面包含 f(x,y) 的积分，想要求 f(x,y) 的具体形式就要把该积分去掉，那么前人就想到了用傅立叶变换
线 L 的表示
前面的描述中说的是“在某一个方向上的每一个点的”衰减程度 \rightarrow 某一个方向上的每一个点都可以表示为一条直线 L ，那么我们需要用该方向来定义直线 L ，先来看看 L 的定义
如上图所示，我们设直线方程为 y=kx+b ，直线 L 远离原点的法线方向为\vec{n}=(\cos\theta,\sin\theta)，原点到直线的距离为 p ，则如上图公式所示，直线 L 可表示为 x\cos \theta + y \sin \theta = p \\其中该直线斜率 k 为负数 \alpha < 0 。
可以看到 \theta 即为上图所说的方向， p 为该方向上的每一点放射源，则根据 (\theta, p) 可以定义 Radon 变换
Radon 变换
Radon变换相当于把函数 f(x,y) 通过线积分表示成了另外的一种直线参数的形式，相当于把二维平面 (x,y) 坐标系映射到了直线参数 (\theta, p) 坐标系 \mathcal{R}_L = \int_L f(x, y) \mathrm{d}s \\ \mathcal{R}(\theta, p) = \int_{(\theta, p)}f(x, y)\mathrm{d}s\\ \mathcal{R}(\theta, p) = \underset{x\cos\theta + y\sin \theta = p}{\int} f(x, y) \mathrm{d}s 由于线积分把所有点 (x,y) 约束在了直线 x\cos \theta + y\sin \theta = p 上，故为了能够展开上式使用了 \delta 函数，则展开如下 \mathcal{R}(\theta, p) = \iint f(x, y) \cdot \delta(p - x\cos\theta + y\sin \theta) \mathrm{d}x\mathrm{d}y \\ 接下来是如何根据不同的 (\theta, p) 求出 f(x,y) ，最上面说的傅立叶变换再反变换就得到结果也就是这个
傅立叶变换
在积分形式下得出函数的表达式第一想到的就是傅立叶变换，因为傅立叶变换及其性质里面就有把积分约掉的方法，这里用到的是傅立叶变换的卷积性质，为了能说清楚，自己也清楚一下
卷积： f(t)*g(t) = \int f(\tau)g(t-\tau)\mathrm{d}\tau \\ 傅立叶变换得
\begin{align*} \mathcal{F}[f(t) * g(t)] & = \int \left[ \int f(\tau)g(t-\tau)\mathrm{d}\tau \right]e^{-jwt} \mathrm{d}t \\ & = \int \left[ \int f(\tau)g(t-\tau)\mathrm{d}\tau \right]e^{-jw(\tau + t - \tau)} \mathrm{d}t \\ & = \int f(\tau)e^{-jw\tau} \left[ \int g(t-\tau) e^{-jw(t-\tau)}\mathrm{d}t \right]\mathrm{d}\tau \\ & = \int f(\tau) e^{-jw\tau} \hat{g}(w)\mathrm{d}\tau \\ & = \hat{f}(w) \hat{g}(w) \end{align*} \\
Radon 逆变换
下面就是根据上面的傅立叶变换方法把积分约掉并求出 f 的具体形式，这里我们对 p 求傅立叶变换，很简单的把上面的式子套一遍
\begin{align*} \mathcal{R}(\theta, p) &= \iint f(x, y) \cdot \delta(p - x\cos\theta - y\sin \theta) \mathrm{d}x\mathrm{d}y \\ \mathcal{F}[\mathcal{R}(\theta, p)] &=\int \left[ \iint f(x,y)\delta(p-x\cos \theta-y\sin\theta)\mathrm{d}x\mathrm{d}y \right]e^{-jwp}\mathrm{d}p \\ & = \int \left[ \iint f(x,y)\delta(p-x\cos \theta-y\sin\theta)\mathrm{d}x\mathrm{d}y \right]e^{-jw(p -x\cos \theta - y\sin\theta + x\cos \theta + y\sin\theta)}\mathrm{d}p \\ & = \iint f(x,y) e^{-jw(x\cos \theta + y\sin\theta)} \left[ \int \delta(p-x\cos \theta - y\sin \theta) e^{-jw(p-x\cos \theta-y\sin\theta)}\mathrm{d}p \right]\mathrm{d}x\mathrm{d}y \\ & = \iint f(x,y)e^{-jw(x\cos\theta+y\sin\theta)}\hat{\delta}(w)\mathrm{d}x\mathrm{d}y \\ & = \hat{\delta}(w)\iint f(x,y)e^{-jw(x\cos \theta + y\sin\theta)} \mathrm{d}x\mathrm{d}y \\ & = \hat{f}(w\cos\theta, w\sin\theta)\hat{\delta}(w) \end{align*} \\ 最后我们得到的等式 \mathcal{F}[\mathcal{R}(\theta, p)] = \hat{f}(w\cos\theta, w\sin\theta)\hat{\delta}(w) \\ 还有 \delta 函数的傅立叶变换为常数 1 ，故最终的等式为
\hat{f}(w\cos\theta, w\sin\theta) = \mathcal{F}[\mathcal{R}(\theta, p)] \\
则根据式子可以得到，对 \theta 方向的检测结果进行傅立叶变换得到的数据填充到过原点的 \theta 角度的那条线上，最后进行反傅立叶变换即可得到原函数形式
实现
搞通了还是感觉蛮简单的，最后还差个实现，改天再写