医学图像重建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(,s)$的进程。这个公式的推导我是这样了解的，在某一个特定视点下*（水平方向投影）的s和f(x,y)存在这个联系：

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

假如考虑到视点，那么就变成：

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

现在咱们来了解什么是笔直sensor的直线上。

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

现在咱们需求吧笔直sensor的直线上束缚转换成数学的表明，直线的表明常见便是y-kx-b=0的方式，可是咱们需求用p和\theta来表明直线，因而直线能够表明为: $xcos+ysin=p$

那么上面公式就变成

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

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

再看下书中的界说：

再看下我的了解：

• 狄拉克函数是一个概念，抱负情况下便是一个相似在0的当地无穷大，非0当地等于0的脉冲函数。当然这种函数可能不行导或许没有什么很好的性质。因而咱们能够用一些性质比较好的函数来模仿这种效应。也便是高斯函数。公式为：

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

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

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

咱们积分束缚那里，现在咱们能够写成这样的方式：

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

由于当xcos\theta+ysin\theta-p不等于0的时分，其实$(xcos+ysin-p)$能够看作是0(我的了解哈).

至此，咱们了解了什么是radon改换，是一个多视点投影的正向进程。

中心切片定理

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

看下图：

• 右边黑乎乎的是f(x,y)的傅里叶改变图画。

• f(x,y)沿着某一个方向投影得到绿色的1D散布，这个是radon进程。

• 然后把1D投影散布做傅里叶改换得到赤色1D频域散布。

• 这个中心切片定理关键便是说，这个赤色的1D散布，其实是等于右图傍边红线上的数据。

• 这样，咱们就建立起来了，投影数据和f(x,y)的傅里叶改换图画的联系，之后通过2D反傅里叶改换就能够得到f(x,y)的图画了。这便是重建。

关键在于，中心切片定理是怎么证明的。 首要界说$P(,)$为投影的$p(s,)$的傅里叶改换，$F(cos,sin)$为f(x,y)的傅里叶改换。

咱们需求证明：

$P(,)=F(cos,sin)$

这儿咱们需求先知道狄拉克函数的两个性质：

我是这样了解的，狄拉克函数由于积分和为1，所以能够看成对f(x)求希望的进程。假如x变成了ax，倍数改变。那么意味着希望概率的稀释。所以稀释的倍数越大，终究要乘上稀释倍数的倒数。

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

下面推导进程我手写了：

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

• 1维傅里叶改换

• 1维傅里叶反改换

• 2维傅里叶反改换

FBP filtered backprojection 滤波反投影算法

其实通过上面的进程进行重建，也便是反投影，是会得到比较模糊的图画的。由于$P(w,)$在构建的时分，是会呈现密度不均匀的情况。由于每一个视点都是过原点的，所以越靠近原点的密度就越高，约远离原点的区域的密度越低。或许能够说，在傅里叶空间，低频区域的密度高，权重就会过大。因而为了消除这种模糊的效果，要对傅里叶空间进行加权纠正，使其密度均匀。

做法便是，在得到的投影傅里叶改换空间中，乘上wx2+wy2\sqrt{w_x^2+w_y^2}wx2​+wy2​​。

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

咱们需求得到的f(x,y)便是通过反傅里叶改变的方法:

f(x,y)=∫02∫0∞F(w,)e2iw(xcos+ysin)∣w∣dwdf(x,y)=\int_0^{2\pi} \int_{0}^{\infin} F(w,\theta)e^{2\pi i w(xcos\theta+ysin\theta)}|w|dwd\thetaf(x,y)=∫02​∫0∞​F(w,)e2iw(xcos+ysin)∣w∣dwd

依据中心切片定理，原始图画的傅里叶改换$F(w,)$等于投影的傅里叶改换$P(w,)$，投影的$P(w,)$由于密度不均匀需求通过$|w|$进行纠正，这个w的绝对值其实便是wx2+wy2\sqrt{w_x^2+w_y^2}wx2​+wy2​​。纠正后的投影的傅里叶改换咱们写作$Q(w,)$，公式如下：

f(x,y)=∫02∫0∞Q(w,)e2iw(xcos+ysin)dwdf(x,y)=\int_0^{2\pi} \int_{0}^{\infin} Q(w,\theta)e^{2\pi i w(xcos\theta+ysin\theta)}dwd\thetaf(x,y)=∫02​∫0∞​Q(w,)e2iw(xcos+ysin)dwd

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

f(x,y)=∫02q(xcos+ysin,)df(x,y)=\int_0^{2\pi} q(xcos\theta+ysin\theta,\theta)d\thetaf(x,y)=∫02​q(xcos+ysin,)d

这个便是终究的成果了。从这儿其实能够看到有两种方法来做反投影：

1. 向咱们之前说的，咱们对投影p做傅里叶改换得到P，然后对P做加权纠正得到Q，然后由于Q和F依据中心切片定理是相等的，所以对F做2维反傅里叶改换得到f；
2. 依据FBP最后的公式推导，咱们能够对投影p做傅里叶改换得到P，对P做加权纠正得到Q，然后做1维的反傅里叶改变重新得到通过加权纠正的p’也便是公式中的q。然后依据FBP的公式推导，能够发现f和q存在必定的联系，能够得到f。