医学图像重建1 | Radon变换，滤波反投影算法，中心切片定理

转自微信公众号：机器学习炼丹术
笔记：陈亦新

概述

医学图像重建的目的就是得到上图的f(x,y)的图像。我们只能获取到投影的数据，也就是右边的sensor检测到的强度信息。当然上图来看，是把一个2D的图像投影成了1D的数据，那么这样肯定是无法复原的。

在投影的过程中，并不是上述这一个角度。上述的投影角度为0，是水平从左到右的。假设我们对角度旋转180度，每旋转1度都进行一次投影。那么最终我们会有180个1D的投影数据，然后如何从这些1D投影数据还原2D的原始图像就是我们所说的重建算法。

Radon变换

这个变换讲述的就是将2D物体投影成1D的过程。2D的两个维度记作x和y，1D的数据只有1个维度，我们记作s。但是我们还需要考虑这个radon变成的1D其实是在某一个特定投影角度下的1D数据，所以其实上是还要加上角度的变量 $\theta$ .

用作数学的表示，那就是radon变换就是将f(x,y)变成 $R(\theta,s)$ 的过程。这个公式的推导我是这样理解的，在某一个特定角度下*（水平方向投影）的s和f(x,y)存在这个关系：

$R(s)=\iint_{(x,y)要在水平线上} f(x,y) dxdy$

如果考虑到角度，那么就变成：

$R(\theta,s)=\iint_{(x,y)要在垂直sensor的直线上} f(x,y) dxdy$

现在我们来理解什么是垂直sensor的直线上。

这个图展示了一个投影角度为 $\theta$ 的通用情况。 $R(\theta,s)$ 的s的物理含义是直线与2D坐标原点的垂直距离，也是1D投影距离1D投影坐标原点的距离，就是上图中的p。

现在我们需要吧垂直sensor的直线上约束转换成数学的表示，直线的表示常见就是y-kx-b=0的形式，但是我们需要用p和 $\theta$ 来表示直线，因此直线可以表示为: $xcos\theta+ysin\theta=p$

那么上面公式就变成

$R(\theta,s)=\iint_{xcos\theta+ysin\theta-p=0} f(x,y) dxdy$

但是这还不够好，你在积分下面加个约束，那后续做什么运算都不行。就像是带约束的优化问题有拉格朗日算子法一样，这里我们引入狄拉克 $\delta$ 分布函数。看下百度的结果：

再看下书中的定义：

再看下我的理解：

狄拉克函数是一个概念，理想情况下就是一个类似在0的地方无穷大，非0地方等于0的脉冲函数。当然这种函数可能不可导或者没有什么很好的性质。因此我们可以用一些性质比较好的函数来模拟这种效应。也就是高斯函数。公式为：

$\delta(x)=(\frac{n}{\pi})^{0.5} e^{-nx^2}$

其中n越大，曲线越来越窄，图像如下：

为什么这样定义狄拉克函数呢？首先，这样的定义是符合狄拉克函数的广义概念的，所以他是狄拉克函数的一种定义。然后就是这样定义会有很多很方便的性质，后续用到的时候我再做介绍。

我们积分约束那里，现在我们可以写成这样的形式：

$R(\theta,s)=\iint f(x,y) \cdot \delta(xcos\theta+ysin\theta-p) dxdy$

因为当xcos\theta+ysin\theta-p不等于0的时候，其实 $\delta(xcos\theta+ysin\theta-p)$ 可以看作是0(我的理解哈).

至此，我们理解了什么是radon变换，是一个多角度投影的正向过程。

中心切片定理

中心切片定理是断层扫描成像的理论基础。这个定理还可以叫做：投影切片定理和傅里叶中心切片定理。二维图像的中心切片定义指出：二维图像f(x,y)的 $\theta$ 角度的投影 $p(s)$ 的傅里叶变换 $p(\omega)$ 等于函数f(x,y)的傅里叶变换 $F(\omega cos\theta,\omega sin\theta)$ 同样沿着\theta角度过原点的片段

看下图：

右边黑乎乎的是f(x,y)的傅里叶变化图像。
f(x,y)沿着某一个方向投影得到绿色的1D分布，这个是radon过程。
然后把1D投影分布做傅里叶变换得到红色1D频域分布。
这个中心切片定理关键就是说，这个红色的1D分布，其实是等于右图当中红线上的数据。
这样，我们就建立起来了，投影数据和f(x,y)的傅里叶变换图像的关系，之后通过2D反傅里叶变换就可以得到f(x,y)的图像了。这就是重建。

关键在于，中心切片定理是如何证明的。首先定义 $P(\omega,\theta)$ 为投影的 $p(s,\theta)$ 的傅里叶变换， $F(\omega cos\theta,\omega sin\theta)$ 为f(x,y)的傅里叶变换。

我们需要证明：

$P(\omega,\theta) = F(\omega cos\theta,\omega sin\theta)$

这里我们需要先知道狄拉克函数的两个性质：

我是这样理解的，狄拉克函数因为积分和为1，所以可以看成对f(x)求期望的过程。如果x变成了ax，倍数变化。那么意味着期望概率的稀释。所以稀释的倍数越大，最终要乘上稀释倍数的倒数。

这个也简单，我理解为是期望的平移，从f(0)移动到了f(a).

下面推导过程我手写了：

微信图片_20221202133915(1).jpg

傅里叶变换与反傅里叶变换公式

1维傅里叶变换

1维傅里叶反变换

2维傅里叶反变换

FBP filtered backprojection 滤波反投影算法

其实通过上面的过程进行重建，也就是反投影，是会得到比较模糊的图像的。因为 $P(w,\theta)$ 在构建的时候，是会出现密度不均匀的情况。因为每一个角度都是过原点的，所以越靠近原点的密度就越高，约远离原点的区域的密度越低。或者可以说，在傅里叶空间，低频区域的密度高，权重就会过大。因此为了消除这种模糊的效果，要对傅里叶空间进行加权矫正，使其密度均匀。

做法就是，在得到的投影傅里叶变换空间中，乘上 $\sqrt{w_x^2+w_y^2}$ 。

现在我们来推导FBP算法，也就是反投影算法。

我们需要得到的f(x,y)就是通过反傅里叶变化的方法:

$f(x,y)=\int_0^{2\pi} \int_{0}^{\infin} F(w,\theta)e^{2\pi i w(xcos\theta+ysin\theta)}|w|dwd\theta$

根据中心切片定理，原始图像的傅里叶变换 $F(w,\theta)$ 等于投影的傅里叶变换 $P(w,\theta)$ ，投影的 $P(w,\theta)$ 因为密度不均匀需要通过 $|w|$ 进行矫正，这个w的绝对值其实就是 $\sqrt{w_x^2+w_y^2}$ 。矫正后的投影的傅里叶变换我们写作 $Q(w,\theta)$ ，公式如下：

$f(x,y)=\int_0^{2\pi} \int_{0}^{\infin} Q(w,\theta)e^{2\pi i w(xcos\theta+ysin\theta)}dwd\theta$

我们对w变量做1维的反傅里叶变换，

$f(x,y)=\int_0^{2\pi} q(xcos\theta+ysin\theta,\theta)d\theta$

这个就是最终的结果了。从这里其实可以看到有两种方法来做反投影：

向我们之前说的，我们对投影p做傅里叶变换得到P，然后对P做加权矫正得到Q，然后因为Q和F根据中心切片定理是相等的，所以对F做2维反傅里叶变换得到f；
根据FBP最后的公式推导，我们可以对投影p做傅里叶变换得到P，对P做加权矫正得到Q，然后做1维的反傅里叶变化重新得到经过加权矫正的p'也就是公式中的q。然后根据FBP的公式推导，可以发现f和q存在一定的关系，可以得到f。