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

专栏地址：『youcans 的 OpenCV 例程 300篇 - 总目录』

【第 7 章：图像复原与重建】
110. 投影和雷登变换
111. 雷登变换反投影重建图像
112. 滤波反投影重建图像

【youcans 的 OpenCV 例程 300 篇】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-28