【youcans 的 OpenCV 例程 300 篇】111. 雷登变换反投影重建图像

7. 投影重建图像

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

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

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

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

反投影一个点形成部分图像的过程，是将直线 $L(ρj,θk)L(\rho_j,\theta_k)$ 复制到图像上，直线上每点的灰度值是 $g(ρj,θk)g(\rho_j,\theta_k)$ 。遍历投影信号中的每个点，得到：
$fθ(x,y)=g(ρ,θ)=g(xcosθ+ysinθ)f_{\theta}(x,y) = g(\rho,\theta) = g(x cos \theta + ysin \theta)$
对所有反投影图像积分，得到最终的图像：
$\int_0^{\pi}f_{\theta}(x,y)d\theta$
其离散形式为：

$\sum_{\theta=0}^{\pi}f_{\theta}(x,y)$

反投影的图像有时称为层图，可以理解为投影图像的一个近似。

例程 9.24：雷登变换反投影重建图像

    # 9.24: 雷登变换反投影重建图像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]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()