【youcans 的 OpenCV 例程 300 篇】112. 滤波反投影重建图像

7. 投影重建图像

图像重建的基本思想，就是通过探测物体的投影数据，重建物体的实际内部构造。

7.3 滤波反投影重建图像（Filtered back projection）

直接由正弦图得到反投影图像，会存在严重的模糊，这是早期 CT 系统所存在的问题。

傅立叶中心切片定理表明，投影的一维傅立叶变换是得到投影区域的二维傅立叶变换的切片。滤波反投影重建算法在反投影前将每一个采集投影角度下的投影进行卷积处理，从而改善点扩散函数引起的形状伪影，有效地改善了重建的图像质量。

以 $f (x, y)$ 表示需要重建的图像， $F (u, v)$ 的二维傅里叶反变换是：
$\int^{\infty}_{-\infty} \int^{\infty}_{-\infty} F(u,v) e^{j 2 \pi (ux+vy)} dudv$
利用傅里叶切片定理，可以得到：
$\int^{\pi}_{0} \big[ \int^{\infty}_{-\infty} |\omega|G(\omega,\theta) e^{j 2 \pi \omega \rho} d\omega \big] d\theta$
括号 [] 内部是一个一维傅里叶反变换，可以认为这是一个一维滤波器的传递函数。由于 $∣ω∣|\omega|$ 是一个不可积的斜坡函数（Slope function），可以通过对斜坡加窗进行限制，典型地如汉明窗（Hamming window）、韩窗（Hann window）。

该式也可以使用空间卷积来实现：
$\int^{\pi}_{0} \big[ \int^{\infty}_{-\infty} |\omega|G(\omega,\theta) e^{j 2 \pi \omega \rho} d\omega \big] d\theta \\ = \int^{\pi}_{0} \big[ \int^{\infty}_{-\infty} g(\rho,\theta) s(xcos \theta + ysin \theta - \rho) d\rho \big] d\theta \\$
这表明，将对应的投影 $g(ρ,θ)g(\rho, \theta)$ 与斜坡滤波器传递函数 $s(ρ)s(\rho)$ 的傅里叶反变换进行卷积，可以得到角度 $θ\theta$ 的各个反投影，整个反投影图像可以通过对所有反投影图像积分得到。

例程 9.25：滤波反投影重建图像

SL 滤波器是使用 Sinc 函数对斜坡滤波器进行截断产生，其离散形式为：
$hSL(nδ)=−2π2δ2(4∗n2−1),n=0,±1,±2,...h_{SL}(n \delta) = - \frac{2}{\pi ^2 \delta ^2 (4*n^2 -1)}, \ n=0,\pm 1,\pm 2,...$
利用投影值和 SL 滤波器进行卷积，然后再进行反投影，就可以实现图像重建。

   # 9.25: 滤波反投影重建图像from scipy import ndimagedef discreteRadonTransform(image, steps):  # 离散雷登变换channels = image.shape[0]resRadon = np.zeros((channels, channels), dtype=np.float32)for s in range(steps):rotation = ndimage.rotate(image, -s * 180/steps, reshape=False).astype(np.float32)resRadon[:, s] = sum(rotation)return resRadondef inverseRadonTransform(image, steps):  # 雷登变换反投影重建图像channels = image.shape[0]backproject = np.zeros((steps, channels, channels))for s in range(steps):expandDims = np.expand_dims(image[:, s], axis=0)repeat = expandDims.repeat(channels, axis=0)backproject[s] = ndimage.rotate(repeat, s * 180/steps, reshape=False).astype(np.float32)invRadon = np.sum(backproject, axis=0)return invRadondef SLFilter(N, d):  # SL 滤波器, Sinc 函数对斜坡滤波器进行截断filterSL = np.zeros((N,))for i in range(N):filterSL[i] = - 2 / (np.pi**2 * d**2 * (4 * (i-N/2)**2 - 1))return filterSLdef filterInvRadonTransform(image, steps):  # 滤波反投影重建图像channels = len(image[0])backproject = np.zeros((steps, channels, channels))  # 反投影filter = SLFilter(channels, 1)  # SL 滤波器for s in range(steps):value = image[:, s]  # 投影值valueFiltered = np.convolve(filter, value, "same")  # 投影值和 SL 滤波器进行卷积filterExpandDims = np.expand_dims(valueFiltered, axis=0)filterRepeat = filterExpandDims.repeat(channels, axis=0)backproject[s] = ndimage.rotate(filterRepeat, s * 180/steps, reshape=False).astype(np.float32)filterInvRadon = np.sum(backproject, axis=0)return filterInvRadon# 读取原始图像img1 = cv2.imread("../images/Fig0534a.tif", 0)  # flags=0 读取为灰度图像img2 = cv2.imread("../images/Fig0534c.tif", 0)# 雷登变换imgRadon1 = discreteRadonTransform(img1, img1.shape[0])  # Radon 变换imgRadon2 = discreteRadonTransform(img2, img2.shape[0])# 雷登变换反投影重建图像imgInvRadon1 = inverseRadonTransform(imgRadon1, imgRadon1.shape[0])imgInvRadon2 = inverseRadonTransform(imgRadon2, imgRadon2.shape[0])# 滤波反投影重建图像imgFilterInvRadon1 = filterInvRadonTransform(imgRadon1, imgRadon1.shape[0])imgFilterInvRadon2 = filterInvRadonTransform(imgRadon2, imgRadon2.shape[0])plt.figure(figsize=(9, 7))plt.subplot(231), plt.axis('off'), plt.title("origin image"), plt.imshow(img1, 'gray')  # 绘制原始图像plt.subplot(232), plt.axis('off'), plt.title("inverse Radon trans"), plt.imshow(imgInvRadon1, 'gray')plt.subplot(233), plt.axis('off'), plt.title("filter inv-Radon trans"), plt.imshow(imgFilterInvRadon1, 'gray')plt.subplot(234), plt.axis('off'), plt.title("origin image"), plt.imshow(img2, 'gray')plt.subplot(235), plt.axis('off'), plt.title("inverse Radon trans"), plt.imshow(imgInvRadon2, 'gray')plt.subplot(236), plt.axis('off'), plt.title("filter inv-Radon trans"), plt.imshow(imgFilterInvRadon2, 'gray')plt.tight_layout()plt.show()