冈萨雷斯《数字图像处理》学习笔记（4）--图像复原与重建(含傅里叶切片定理推导)

一、图像复原模型

若H是线性的，空间不变的过程，则退化图像在空间域通过下式给出:
g(x,y)=h(x,y)∗f(x,y)+δ(x,y)g(x,y)=h(x,y)*f(x,y)+\delta(x,y)g(x,y)=h(x,y)∗f(x,y)+δ(x,y)其中，h(x, y)是退化函数的空间表示，符号∗*∗表示卷积。空间域的卷积和频域的乘法组成了一个傅里叶变换对，所以可以用等价的频域表示来写出前面的模型:
G(u,v)=H(u,v)F(u,v)+N(u,v)G(u,v)=H(u,v)F(u,v)+N(u,v)G(u,v)=H(u,v)F(u,v)+N(u,v)

二、图像重建模型

CT成像原理的实质是衰减系数成像。主要涉及朗伯比尔定律和Rodon变换

1. 朗伯比尔定律
具体的可以通过百度了解。简单来说，就是每种物质的核外电子数目都不一样，对X光的吸收也不一样。通过检测前后能量的差异可以求出物体对X光的吸收能力的大小，也就是衰减系数。把不同的衰减系数对应到不同的像素值，就得到X光照片。

2. Rodon变换
笛卡尔坐标系中的一条直线可以由 y=ax+by=ax+by=ax+b 来描述，可以写出它的法线式（即垂直于 y=ax+by=ax+by=ax+b 的垂线段，该垂线段所在直线的倾斜角为 θ\thetaθ，ρ\rhoρ 是该线段的长度）：xcosθ+ysinθ=ρxcos\theta+ysin\theta=\rhoxcosθ+ysinθ=ρ

平行射线束的投影可由这样的一组直线建模，看下图：

投影平面上任意一点的坐标 (ρj,θk)(\rho_j,\theta_k)(ρj,θk) 由沿着线 xcosθk+ysinθk=ρjxcos\theta_k+ysin\theta_k=\rho_jxcosθk+ysinθk=ρj 的射线和给出。该射线和是一个线积分，由下式给出：
g(ρj,θk)=∫−∞∞∫−∞∞f(x,y)δ(xcosθk+ysinθk−ρj)dxdyg(\rho_j,\theta_k )=\int_{-\infty }^{\infty }\int_{-\infty }^{\infty }f(x,y)\delta (xcos\theta_k+ysin\theta_k-\rho_j)dxdyg(ρj,θk)=∫−∞∞∫−∞∞f(x,y)δ(xcosθk+ysinθk−ρj)dxdy
其中，δ(x)={1,if x=00,if x≠0\delta (x)= \begin{cases} 1,& \text{ if } x=0 \\ 0, & \text{ if } x\neq 0 \end{cases}δ(x)={1,0, if x=0 if x=0
这个函数意味着积分是沿着 xcosθk+ysinθk=ρjxcos\theta_k+ysin\theta_k=\rho_jxcosθk+ysinθk=ρj 这条射线计算的。考虑 ρ\rhoρ 和 θ\thetaθ 的所有值，可以推广到：
g(ρ,θ)=∫−∞∞∫−∞∞f(x,y)δ(xcosθ+ysinθ−ρ)dxdyg(\rho,\theta)=\int_{-\infty }^{\infty }\int_{-\infty }^{\infty }f(x,y)\delta (xcos\theta+ysin\theta-\rho)dxdyg(ρ,θ)=∫−∞∞∫−∞∞f(x,y)δ(xcosθ+ysinθ−ρ)dxdy

这就是Rodon变换。

接下来我们需要得到不同角度的反投影再通过求和来重建图像。为了得到反投影的表达式，我们从固定单个点 g(ρj,θk)g(\rho_j,\theta_k )g(ρj,θk) 开始，只需将线 L(ρj,θk)L(\rho_j,\theta_k )L(ρj,θk) 复制到图像上即可，其中沿该条线的所有的点的值是 g(ρj,θk)g(\rho_j,\theta_k )g(ρj,θk)。我们对投影信号的中的所有 ρj\rho_jρj 都重复这个过程，得到如下表达式：
fθk(x,y)=g(xcosθk+ysinθk,θk)f_{\theta_k}(x,y)=g(xcos\theta_k+ysin\theta_k,\theta_k)fθk(x,y)=g(xcosθk+ysinθk,θk)
很明显，该公式适合所有的角度。推广到一般（即所有角度）：
fθ(x,y)=g(xcosθ+ysinθ,θ)f_{\theta}(x,y)=g(xcos\theta+ysin\theta,\theta)fθ(x,y)=g(xcosθ+ysinθ,θ)
最后对所有反投影图像积分，得到最后的图像：
f(x,y)=∫0πfθ(x,y)dθf(x,y)=\int_{0}^{\pi}f_{\theta}(x,y)d\thetaf(x,y)=∫0πfθ(x,y)dθ
因为0度和180度互为镜像，因此求和操作只需执行到180度之前的最后一个角度增量。但使用这种方法得到的结果比较模糊，不能令人满意。但我们可以用傅里叶切片定理重新表示反投影方法来得到明显增强的结果。

傅里叶切片定理
G(ω,θ)=∫−∞∞g(ρ,θ)e−j2πωρdρ=F(ωcosθ,ωsinθ)G(\omega ,\theta )=\int_{-\infty }^{\infty } g(\rho ,\theta ) e^{-j2\pi \omega \rho }d\rho =F(\omega cos\theta ,\omega sin\theta )G(ω,θ)=∫−∞∞g(ρ,θ)e−j2πωρdρ=F(ωcosθ,ωsinθ)

（中间的式子由一维傅里叶变换给出，推导可见我上篇博客）冈萨雷斯《数字图像处理》学习笔记（3）–频率域滤波(含傅里叶变换推导)

该定理表明，一个投影的傅里叶变换，是得到该投影区域的二维傅里叶变换的一个切片。下图说明了这个结果：

傅里叶切片定理证明：
G(ω,θ)=∫−∞∞g(ρ,θ)e−j2πωρdρ=∫−∞∞∫−∞∞∫−∞∞f(x,y)δ(xcoxθk+ysinθk−ρj)e−j2πωρdxdydρ=∫−∞∞∫−∞∞f(x,y)[∫−∞∞δ(xcoxθk+ysinθk−ρj)e−j2πωρdρ]dxdy=∫−∞∞∫−∞∞f(x,y)e−j2πω(xcosθ+ysinθ)dxdy\begin{aligned} G(\omega ,\theta )&=\int_{-\infty }^{\infty } g(\rho ,\theta ) e^{-j2\pi \omega \rho }d\rho \quad \quad \quad\\ &=\int_{-\infty }^{\infty }\int_{-\infty }^{\infty }\int_{-\infty }^{\infty } f(x,y)\delta(xcox\theta_k+ysin\theta_k-\rho_j) e^{-j2\pi \omega \rho }dxdyd\rho \quad \quad \quad\\ =&\int_{-\infty }^{\infty } \int_{-\infty }^{\infty } f(x,y)\left [ \int_{-\infty }^{\infty } \delta(xcox\theta_k+ysin\theta_k-\rho_j) 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+ysin\theta) }dxdy \end{aligned} G(ω,θ)==∫−∞∞g(ρ,θ)e−j2πωρdρ=∫−∞∞∫−∞∞∫−∞∞f(x,y)δ(xcoxθk+ysinθk−ρj)e−j2πωρdxdydρ∫−∞∞∫−∞∞f(x,y)[∫−∞∞δ(xcoxθk+ysinθk−ρj)e−j2πωρdρ]dxdy=∫−∞∞∫−∞∞f(x,y)e−j2πω(xcosθ+ysinθ)dxdy
然后我们令 u=ωcosθu=\omega cos\thetau=ωcosθ 和 v=ωsinθv=\omega sin\thetav=ωsinθ，得到：
G(ω,θ)=[∫−∞∞∫−∞∞f(x,y)e−j2π(ux+vy)dxdy]u=ωcosθ；v=ωsinθ=[F(u,v)]u=ωcosθ；v=ωsinθ=F(ωcosθ,ωsinθ)\begin{aligned} G(\omega ,\theta )&=\left [ \int_{-\infty }^{\infty}\int_{-\infty }^{\infty}f(x,y) e^{-j2\pi (ux+vy)}dxdy\right ]_{u=\omega cos\theta；v=\omega sin\theta}\\ &=[F(u,v)]_{u=\omega cos\theta；v=\omega sin\theta}\\ &=F(\omega cos\theta,\omega sin\theta) \end{aligned}G(ω,θ)=[∫−∞∞∫−∞∞f(x,y)e−j2π(ux+vy)dxdy]u=ωcosθ；v=ωsinθ=[F(u,v)]u=ωcosθ；v=ωsinθ=F(ωcosθ,ωsinθ)
证毕！

接下来利用滤波反投影重建图像.
给定F(u,v)，使用傅里叶反变换得到f(x,y):
f(x,y)=∫−∞∞∫−∞∞F(u,v)ej2π(ux+vy)dudvf(x,y)=\int_{-\infty }^{\infty }\int_{-\infty }^{\infty }F(u,v)e^{j2\pi (ux+vy)}dudvf(x,y)=∫−∞∞∫−∞∞F(u,v)ej2π(ux+vy)dudv
令 u=ωcosθu=\omega cos\thetau=ωcosθ 和 v=ωsinθv=\omega sin\thetav=ωsinθ
由雅可比行列式：
∣J∣=∣∂u∂ω∂u∂θ∂v∂ω∂v∂θ∣=∣cosθ−ωsinθsinθωcosθ∣=ω|J|=\begin{vmatrix} \frac{\partial u}{\partial \omega } & \frac{\partial u}{\partial \theta }\\ \frac{\partial v}{\partial \omega }& \frac{\partial v}{\partial \theta } \end{vmatrix}=\begin{vmatrix} cos \theta & -\omega sin\theta\\ sin\theta & \omega cos\theta \end{vmatrix}= \omega ∣J∣=∣∣∣∣∂ω∂u∂ω∂v∂θ∂u∂θ∂v∣∣∣∣=∣∣∣∣cosθsinθ−ωsinθωcosθ∣∣∣∣=ω
我们得到 dudv=ωdωdθdudv=\omega d\omega d\thetadudv=ωdωdθ
则可把前面的积分表示成极坐标的形式:
f(x,y)=∫02π∫0∞F(ωcosθ,ωsinθ)ej2πω(xcosθ+ysinθ)ωdωdθf(x,y)=\int_{0 }^{2\pi}\int_{0 }^{\infty}F(\omega cos\theta ,\omega sin\theta ) e^{j2\pi\omega (xcos\theta +ysin\theta)}\omega d\omega d\thetaf(x,y)=∫02π∫0∞F(ωcosθ,ωsinθ)ej2πω(xcosθ+ysinθ)ωdωdθ
由上面的傅里叶切片定理有：

f(x,y)=∫02π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθf(x,y)=\int_{0 }^{2\pi}\int_{0 }^{\infty}G(\omega, \theta)e^{j2\pi\omega (xcos\theta +ysin\theta)}\omega d\omega d\thetaf(x,y)=∫02π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ
又根据 cos(θ+π)=−cos(θ)cos(\theta+\pi)=-cos(\theta)cos(θ+π)=−cos(θ) 和 sin(θ+π)=−sin(θ)sin(\theta+\pi)=-sin(\theta)sin(θ+π)=−sin(θ)，所以有：

G(ω,θ+π)=∫−∞∞∫−∞∞f(x,y)e−j2πω(−xcoxθ−ysin⁡θ)dxdy=∫−∞∞∫−∞∞f(x,y)e−j2π(−ω)(xcoxθ+ysin⁡θ)dxdy=G(−ω,θ)\begin{aligned} G(\omega ,\theta+\pi )&=\int_{-\infty }^{\infty}\int_{-\infty }^{\infty}f(x,y)e^{-j2\pi \omega(-xcox\theta-y\sin\theta)}dxdy\\ &=\int_{-\infty }^{\infty}\int_{-\infty }^{\infty}f(x,y)e^{-j2\pi (-\omega)(xcox\theta+y\sin\theta)}dxdy\\ &=G(-\omega ,\theta ) \end{aligned} G(ω,θ+π)=∫−∞∞∫−∞∞f(x,y)e−j2πω(−xcoxθ−ysinθ)dxdy=∫−∞∞∫−∞∞f(x,y)e−j2π(−ω)(xcoxθ+ysinθ)dxdy=G(−ω,θ)

则上式积分可以拆分成这样：
f(x,y)=∫02π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ=∫0π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ+∫π2π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ=∫0π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ+∫0π∫0∞G(ω,θ+π)ej2πω(−xcosθ−ysinθ)ωdωdθ=∫0π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ+∫0π∫0∞G(−ω,θ)ej2π(−ω)(xcosθ+ysinθ)ωdωdθ→换元t=−ω=∫0π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ+∫0π∫0−∞G(t,θ)ej2π(t)(xcosθ+ysinθ)tdtdθ=∫0π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ−∫0π∫−∞0G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ=∫0π∫−∞∞G(ω,θ)ej2πω(xcosθ+ysinθ)∣ω∣dωdθ\begin{aligned} f(x,y)&=\int_{0 }^{2\pi}\int_{0 }^{\infty}G(\omega, \theta)e^{j2\pi\omega (xcos\theta +ysin\theta)}\omega d\omega d\theta\\ &=\int_{0 }^{\pi}\int_{0 }^{\infty}G(\omega, \theta)e^{j2\pi\omega (xcos\theta +ysin\theta)}\omega d\omega d\theta\\ &\quad +\int_{\pi }^{2\pi}\int_{0 }^{\infty}G(\omega, \theta)e^{j2\pi\omega (xcos\theta +ysin\theta)}\omega d\omega d\theta\\ &=\int_{0 }^{\pi}\int_{0 }^{\infty}G(\omega, \theta)e^{j2\pi\omega (xcos\theta +ysin\theta)}\omega d\omega d\theta\\ &\quad +\int_{0 }^{\pi}\int_{0 }^{\infty}G(\omega, \theta+\pi)e^{j2\pi\omega (-xcos\theta -ysin\theta)}\omega d\omega d\theta\\ &=\int_{0 }^{\pi}\int_{0 }^{\infty}G(\omega, \theta)e^{j2\pi\omega (xcos\theta +ysin\theta)}\omega d\omega d\theta\\ &\quad +\int_{0 }^{\pi}\int_{0 }^{\infty}G(-\omega, \theta)e^{j2\pi(-\omega) (xcos\theta +ysin\theta)}\omega d\omega d\theta\\ \\ &\rightarrow 换元 \ t=-\omega\\ \\ &=\int_{0 }^{\pi}\int_{0 }^{\infty}G(\omega, \theta)e^{j2\pi\omega (xcos\theta +ysin\theta)}\omega d\omega d\theta\\ &\quad +\int_{0 }^{\pi}\int_{0 }^{-\infty}G(t, \theta)e^{j2\pi(t) (xcos\theta +ysin\theta)}t dt d\theta\\ &=\int_{0 }^{\pi}\int_{0 }^{\infty}G(\omega, \theta)e^{j2\pi\omega (xcos\theta +ysin\theta)}\omega d\omega d\theta\\ &\quad -\int_{0 }^{\pi}\int_{-\infty }^{0}G(\omega, \theta)e^{j2\pi\omega (xcos\theta +ysin\theta)}\omega d\omega d\theta\\ &=\int_{0 }^{\pi}\int_{-\infty }^{\infty}G(\omega, \theta)e^{j2\pi\omega (xcos\theta +ysin\theta)}|\omega| d\omega d\theta \end{aligned} f(x,y)=∫02π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ=∫0π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ+∫π2π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ=∫0π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ+∫0π∫0∞G(ω,θ+π)ej2πω(−xcosθ−ysinθ)ωdωdθ=∫0π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ+∫0π∫0∞G(−ω,θ)ej2π(−ω)(xcosθ+ysinθ)ωdωdθ→换元 t=−ω=∫0π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ+∫0π∫0−∞G(t,θ)ej2π(t)(xcosθ+ysinθ)tdtdθ=∫0π∫0∞G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ−∫0π∫−∞0G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ=∫0π∫−∞∞G(ω,θ)ej2πω(xcosθ+ysinθ)∣ω∣dωdθ
前面的公式是平行射线束x射线断层的基本结果。它表明完全的反投影图像f(x,y)是由一组平行射线束投影通过如下步骤得到的：