111. 雷登变换反投影重建图像

7. 图像复原与重建索引

一、投影重建图像

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

二、雷登变换反投影重建图像（Radon transform back projection）

• 空间点 X 通过摄像机 P 被作用到图像平面的图像点 m = PX ， 获得采集到的图像数据，这种投影关系称为摄像机的正向投影 （forward projection） ，简称投影。

• 反向投影是针对图像平面的基本几何元素而言的，图像平面点 m 的反投影是指在摄像机 P 的作用下具有像点 m 的所有空间点的集合。

• 反投影一个点形成部分图像的过程，是将直线 复制到图像上，直线上每点的灰度值是 。遍历投影信号中的每个点，得到：
  

  
  对所有反投影图像积分，得到最终的图像
  
  其离散形式为：
  
  反投影的图像有时称为层图，可以理解为投影图像的一个近似。

三、例程

• 9.25：雷登变换反投影重建图像

    # 9.24: 雷登变换反投影重建图像
    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]
        res = np.zeros((steps, channels, channels))
        for s in range(steps):
            expandDims = np.expand_dims(image[:, s], axis=0)
            repeat = expandDims.repeat(channels, axis=0)
            res[s] = ndimage.rotate(repeat, s * 180/steps, reshape=False).astype(np.float32)
        invRadon = np.sum(res, axis=0)
        return invRadon


    # 读取原始图像
    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])

    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("Radon transform"), plt.imshow(imgRadon1, 'gray')  # 绘制 sinogram 图
    plt.subplot(233), plt.axis('off'), plt.title("Inv-Radon transform"), plt.imshow(imgInvRadon1, 'gray')
    plt.subplot(234), plt.axis('off'), plt.title("origin image"), plt.imshow(img2, 'gray')
    plt.subplot(235), plt.axis('off'), plt.title("Radon transform"), plt.imshow(imgRadon2, 'gray')
    plt.subplot(236), plt.axis('off'), plt.title("Inv-Radon transform"), plt.imshow(imgInvRadon2, 'gray')
    plt.tight_layout()
    plt.show()

四、资料

youcans_的博客：
https://blog.csdn.net/youcans/article/details/123088322