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

Posted 2022-03-02 小白YouCans

欢迎关注 『youcans 的 OpenCV 例程 200 篇』 系列，持续更新中
欢迎关注 『youcans 的 OpenCV学习课』 系列，持续更新中

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

7. 投影重建图像

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

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

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

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

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

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

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

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

   # 9.25: 滤波反投影重建图像
    from scipy import ndimage
    def 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 resRadon

    def 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 invRadon

    def 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 filterSL

    def 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()

（本节完）

---

版权声明：

youcans@xupt 原创作品，转载必须标注原文链接：(https://blog.csdn.net/youcans/article/details/123088438)

Copyright 2022 youcans, XUPT
Crated：2022-2-22

欢迎关注 『youcans 的 OpenCV 例程 200 篇』 系列，持续更新中
欢迎关注 『youcans 的 OpenCV学习课』 系列，持续更新中

【youcans 的 OpenCV 例程200篇】01. 图像的读取（cv2.imread）
【youcans 的 OpenCV 例程200篇】02. 图像的保存（cv2.imwrite）
【youcans 的 OpenCV 例程200篇】03. 图像的显示（cv2.imshow）
【youcans 的 OpenCV 例程200篇】04. 用 matplotlib 显示图像（plt.imshow）
【youcans 的 OpenCV 例程200篇】05. 图像的属性（np.shape）
【youcans 的 OpenCV 例程200篇】06. 像素的编辑（img.itemset）
【youcans 的 OpenCV 例程200篇】07. 图像的创建（np.zeros）
【youcans 的 OpenCV 例程200篇】08. 图像的复制（np.copy）
【youcans 的 OpenCV 例程200篇】09. 图像的裁剪（cv2.selectROI）
【youcans 的 OpenCV 例程200篇】10. 图像的拼接（np.hstack）
【youcans 的 OpenCV 例程200篇】11. 图像通道的拆分（cv2.split）
【youcans 的 OpenCV 例程200篇】12. 图像通道的合并（cv2.merge）
【youcans 的 OpenCV 例程200篇】13. 图像的加法运算（cv2.add）
【youcans 的 OpenCV 例程200篇】14. 图像与标量相加（cv2.add）
【youcans 的 OpenCV 例程200篇】15. 图像的加权加法（cv2.addWeight）
【youcans 的 OpenCV 例程200篇】16. 不同尺寸的图像加法
【youcans 的 OpenCV 例程200篇】17. 两张图像的渐变切换
【youcans 的 OpenCV 例程200篇】18. 图像的掩模加法
【youcans 的 OpenCV 例程200篇】19. 图像的圆形遮罩
【youcans 的 OpenCV 例程200篇】20. 图像的按位运算
【youcans 的 OpenCV 例程200篇】21. 图像的叠加
【youcans 的 OpenCV 例程200篇】22. 图像添加非中文文字
【youcans 的 OpenCV 例程200篇】23. 图像添加中文文字
【youcans 的 OpenCV 例程200篇】24. 图像的仿射变换
【youcans 的 OpenCV 例程200篇】25. 图像的平移
【youcans 的 OpenCV 例程200篇】26. 图像的旋转（以原点为中心）
【youcans 的 OpenCV 例程200篇】27. 图像的旋转（以任意点为中心）
【youcans 的 OpenCV 例程200篇】28. 图像的旋转（直角旋转）
【youcans 的 OpenCV 例程200篇】29. 图像的翻转（cv2.flip）
【youcans 的 OpenCV 例程200篇】30. 图像的缩放（cv2.resize）
【youcans 的 OpenCV 例程200篇】31. 图像金字塔（cv2.pyrDown）
【youcans 的 OpenCV 例程200篇】32. 图像的扭变（错切）
【youcans 的 OpenCV 例程200篇】33. 图像的复合变换
【youcans 的 OpenCV 例程200篇】34. 图像的投影变换
【youcans 的 OpenCV 例程200篇】35. 图像的投影变换（边界填充）
【youcans 的 OpenCV 例程200篇】36. 直角坐标与极坐标的转换
【youcans 的 OpenCV 例程200篇】37. 图像的灰度化处理和二值化处理
【youcans 的 OpenCV 例程200篇】38. 图像的反色变换（图像反转）
【youcans 的 OpenCV 例程200篇】39. 图像灰度的线性变换
【youcans 的 OpenCV 例程200篇】40. 图像分段线性灰度变换
【youcans 的 OpenCV 例程200篇】41. 图像的灰度变换（灰度级分层）
【youcans 的 OpenCV 例程200篇】42. 图像的灰度变换（比特平面分层）
【youcans 的 OpenCV 例程200篇】43. 图像的灰度变换（对数变换）
【youcans 的 OpenCV 例程200篇】44. 图像的灰度变换（伽马变换）
【youcans 的 OpenCV 例程200篇】45. 图像的灰度直方图
【youcans 的 OpenCV 例程200篇】46. 直方图均衡化
youcans 的 OpenCV 例程200篇182.基于形态学梯度的分水岭算法

youcans 的 OpenCV 例程200篇结束语

youcans的OpenCV例程200篇总目录

youcans 的 OpenCV 例程200篇179.图像分割之 GrabCut 图割法（掩模图像）

youcans 的 OpenCV 例程200篇201. 图像的颜色空间转换

OpenCV 例程200篇 目录-202205更新