Radon 变换原理和应用

前言： 承接 Hough 变换，从简单的直线检测说起，推广到 Radon 变换。介绍了 Radon 变换的基本原理和应用。主要是Hough直线检测的拓展，深度较浅。

上启自：详解 Hough 变换

1. 从直线检测说起

1.1 基本问题

在 详解 Hough 变换 中，对直线检测折腾这么久，已经是个成熟的”直线检测者”了，所以基本问题就不从点分析，直接从图出发。即

如何检测下图（a）中的直线？

当然可以使用 Hough 变换，如直接套用那篇博文中自编的代码就可以得到结果，展示如下：

再简单概述下原理：即，将原空间中的点变换到参数空间，每个原空间的点对应参数空间一条直线。原空间中有N个点就对应参数空间中N个直线。参数空间中直线相交就是表明对应原空间的点共线。参数空间中交点重复次数(亮度)越多，就表明原空间中共线的点越多。这些点连起来就是要求的直线。

注：如果使用 Hough变换 博文中的代码，将 \(\theta\) 范围修改下，其实不必 \([0,2\pi]\)，修改为\([0,\pi]\) 即可。

1.2 另辟蹊径

检测之前一般都是进行边缘提取，我们的目标是检测 图(b) 边缘二值图中是否存在直线，以及，如果存在，表示出这条直线。

先抛弃 Hough 变换的思路，考虑另外想法。

（1） 第一种想法：考虑从不同方向把 (b) 图拍扁（术语叫投影）。比如从两个角度（下图 \(A1,A2\) ）拍扁（投影）它，示意如下：

这里使用两个方向 (\(\alpha=0\) 和 \(\alpha=45°\) )的光线 \(A1,A2\) 对边缘进行投影。假设投影到投影面的位置点的数量越多，投影结果越亮。那么图中，灰色箭头标志的地方会非常亮。因为整条直线的像素点都会投影到该处。

不仅如此，事实上，考察所有投影方向，即 \(\alpha\in [0,\pi]\) （注意：因为不分正负，所以 \([0,\pi]\) 和 [\(\pi,2\pi]\) 结果是一样的），最亮的点依旧会是图示中灰色标志的部分。这是直线的特点。

所以，是否我们可以考察所有方向下，对边缘图像进行投影，找到最亮的点（或者达到规定亮度阈值的点，比如存在多条直线）和此时对应的光线方向。 那这条直线不就检测出来了，直线方程不就确定了吗？

先不着急编程实践，再思考另一种想法。

（2）第二种想法：Hough 变换中，有提到可以用 \(\theta\) 和 \(d\) 表示一条直线，如下图示意（不懂为何，可以参考 Hough 变换博文）：

其中， \(d\) 为圆心到直线距离，\(\theta\) 是垂线与 \(x\) 轴夹角。那么这一想法就是使用 穷举法 的思想，即

因为 \(\theta\) 是有范围的，同时因为在图像中检测直线，所以 \(d\) 也是有范围，最大为图像的斜对角距离。于是，我们就考察“所有”的 \(\theta\) 和 \(d\) ，这里的”所有“取决于你的精度，即 \(\theta\) 以及 \(d\) 各自的离散值取值间隔。来表示图像空间内”所有“的直线。

表示出”所有“直线后，下面我们要做的就是，将直线一条一条的与原图相匹配。记录该直线的点与边缘图中点重叠的个数，把该个数记为该线与图的匹配度。比如取所有线的其中四条，作个示意：

其中，\(L1,L2\) 和 \(L3,L4\) 故意取了\(\theta\) 相同，但 \(d\) 不同所表示的直线（图中L1,L2没有标角度和距离，怕线看起来太乱了）。上面所说的匹配度即该条直线与边缘图中像素重叠数，如\(L1\) 与 边缘线没有像素重叠，匹配度就为0，\(L4\)同理。其中这举例的四条直线中，匹配度最高的是\(L2\)。因为它与边缘图中包含直线的所有像素都重合。

不仅如此，事实上，考察所有 \(\theta\)，距离 \(d\) 所表示的直线，匹配度数值最高的依旧会是\(L2\)。

说到这，你可能会发现，其实这就是所提到的第一种想法的另种思路。只是这里的 \(L1,L2,\cdots\) 就是第一种想法中的 一条 光线，如 第二种想法里的 \(L2\) 直线其实就是 第一种想法里 \(A2\) 方向下，最亮点对应的那一条光线。匹配度也就对应于第一种想法里的亮度，也就是重合度。

唯一的区别和需要注意的是，两个想法中的所用角度表示方法不同（但属于无伤大雅的东西）：前者是光线的倾角，后者是光线垂线的倾角。

1.3 编程实现

基于以上想法，编程，进行直线检测：

注：这里使用想法二中的倾角 \(\theta\)，这一倾角的好处和直线表示方法，参看Hough变换博文。

直接给出结果，通过 \(\theta\) 和 \(d\) 确定的直线方程为 $$ d = xcos\theta+ysin\theta $$

完整代码如下：

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取图片
imgPath = "C:\\Users\\zhangwei156\\Desktop\\figure\\waterm.jpeg"
grayImg = cv2.imread(imgPath, 0)
    
# 提取边缘
edgeImg = cv2.Canny(grayImg, 300, 500)

H, W = edgeImg.shape

# 精度
theta_div = 500
d_div = 500
theta_max = np.pi
dMax = np.floor(np.sqrt(H**2 + W**2))

# 离散取值
theta = np.linspace(0, np.pi, theta_div)[:-1]
d = np.linspace(-dMax, dMax, d_div)

# 分辨率
theta_res = np.pi/(theta_div-1)
d_res = 2*dMax/(d_div-1)

# 记录 亮度/匹配度/叠加次数 的模板
forceImg = np.zeros((theta_div, d_div), np.uint16)

mesh = np.meshgrid(theta, d)
for _t, _d in zip(mesh[0].flatten(), mesh[1].flatten()):
    # 分母为 0 情况
    if _t == 0:
        continue
        y = np.arange(0, H)
        x = y * 0 + 1
    else:
        x = np.arange(W)
        y = ((_d-x*np.cos(_t))/np.sin(_t)).astype(np.int64) # theta 和 d 确定的直线
    pixel = 2  # 上下容错像素范围
    # y轴上下容错像素
    for i in range(1, pixel+1):
        x = np.hstack((x,x,x))
        y = np.hstack((y-i,y,y+i))
    # 剔除 y < 0 和 y > H 也即图像外的点
    ind = np.where( (y>=0) & (y<H) )
    if ind[0].shape[0] == 0:
        continue
    y = y[ind[0]]
    x = x[ind[0]]
    index = np.where(edgeImg[y, x]!=0 )[0]
    
    # 取值
    tt = int(_t/theta_res)
    dd = int((_d+dMax)/d_res)
    # 亮度值
    forceImg[tt, dd] = index.shape[0]

# 结果
plt.imshow(forceImg)
plt.show()

# 检查结果正确与否的代码
# 最亮处所代表的 theta 和 d 值
O = np.where(forceImg == np.max(forceImg))
# theta 和 d 的真实值
theta_ = O[0]*theta_res
d_ = O[1]*d_res - dMax
# 在原图中显示检测到的直线
for _o in (zip(theta_, d_)):
    x_ = np.arange(W)
    y_ = ((_o[1]-x_*np.cos(_o[0]))/np.sin(_o[0]))
    ind_ = np.where((y_>0) & (y_<H))
    x_ = x_[ind_[0]]
    y_ = y_[ind_[0]]

plt.imshow(grayImg)    
plt.plot(x_,y_,color = "r")
plt.show()

结果展示：

其中，\(forceImg\) 中最亮点对应的实际 \((\theta, d) = (1.561, 181.383)\) ，即\(\theta = 89.46° ; d=181.4\) 。看来垂线还不是纯竖直的，有一点点偏。

此时，你会有个很“吃惊”的发现。

对比一下前面展示的使用 Hough 变换方法时结果，两者 \(forceImg\) 图是一致的！

想象一下整个过程，实际上，针对上述问题，它就是 Hough 变换的解决方式。它就是 Hough 变换。

换句话说，想法一，想法二以及之前的 Hough 直线检测想法都是同源的！

区别在于：Hough 变换中使用 \(\theta,d\) 穷举的是过每个边缘点的所有直线，而以上两个想法通过 \(\theta,d\) 穷举的是图像空间内所有直线。但 Hough 变换有基于先验知识，换句话说，Hough变换更精准。

所以之后用 Hough 检测直线时，也可以使用 randon 变换函数（上面这个就是最基础的randon变换，就是后面要引申的本文主题，但这里不得不先提下）。

2. Radon 变换

2.1 回顾

回顾上述过程，从某个方向对图像投影，一个特定的 \(\theta,d\) 确定某条唯一的光线，投影结果为一个固定值，即亮度（累加度）。也即 \((\theta,d) \longrightarrow 亮度\) ，这是个函数对应关系。

亮度的计算是光线对应图像上的像素值累加。上例的特殊之处在于，边缘图为二值图，即边缘部分像素值为 1，黑色部分像素值为 0。

想象若图像像素值为任意值，事实上也可以累加。推广到更一般情况，即不再如图像像素那般是一格格的、离散的，而是连续的。那求和其实就是求线积分。

把问题推广到这种程度，事实上就是要展示的 Radon 变换内容。

2.2 Radon 变换

如下图所示，不再是简单的离散像素，而是连续的函数值 \(f(x,y)\)，我们暂且称之为”密度“。一条由\(\theta,d\) 确定的唯一光线，经过 \(f(x,y)\) 物体，得到亮度值为 \(R\) 的投影。

根据以上的思路，此时亮度值应该是 \(L\) 线积分，即

且光线 \(L\) 可以由 \(\theta,d\) 唯一确定，即

也就有

该式子也就是一个全新的函数对应关系，表示的是给定自变量 \(\theta,d\) ，会有与之对应的函数值 \(R(\theta,d)\) ，该值的物理意义为亮度（或者说是叠加度、衰减后的值、匹配度等等）。而 \(\theta,d\) 为连续的，因此将该函数用图像画出来就是对应于前面 \(forceImg\) 的样子，即横、纵坐标为 \(\theta,d\)，亮度大小为 \(R\)。更一般情况是对应一个三维空间，\(x,y\) 轴为 \(\theta,d\) ，\(z\) 轴对应为 \(R\)。

下面就是对积分的计算，该部分为数学中线积分求解问题，不详述。

法一：直接求解

当给定 \(\theta,d\) 时，此时 R 为 $$ R = \int_{y = \frac{d_1-x \cos\theta_1}{\sin\theta_1}}f(x,y)ds $$ 将 \(y\) 替换一下，高等数学中的线积分求解问题。这是一种思路。

法二：\(\delta\) 函数 \(\delta\) 函数是一个广义函数。简单形式如下：

结合上述直线方程 \(L(\theta,d):\ d-xcos\theta-ysin\theta = 0\)， 则有

最后，radon 变换方程写为： $$ R(\theta,d) = \iint f(x,y)\cdot \delta(d-xcos\theta+y\sin\theta)dxdy $$ 当给定 \(\theta,d\) 时，同样带入计算。

可以总结一下，Radon 变换就是对应\((\theta,d)\) 确定的光线与空间中密度函数的积分。

2.3 Radon 逆变换

想象一下，如果我从不同角度对物体进行照射，得到了每个方向下的投影，即所有 \(\theta,d\) 对应的 \(R\) 值。那么反过来，已知 \(R(\theta,d)\) ，是否也就唯一确定 \(f(x,y)\) 的情况。

至于如何已知 \(R(\theta,d)\) 反求 \(f(x,y)\) ，现在用不到，挖坑，之后用到再探究其方法。

3. 应用

前面已经提到了一个应用：直线检测（同 Hough 检测一致）

还有个经典应用，即 CT 断层成像的重建。

CT 原理其实就是借助 x 射线照射组织器官，不同组织器官对 x 的射线衰减程度（或系数；对应于上面的 \(f(x,y)\) ）不同，一条 x 射线穿过一系列组织器官后，经衰减后，到达采集设备，得到一个衰减后的强度值（对应于上面的 \(R\)）。

我们采集得到一系列的 \(R\) 情况，就可以通过 \(radon\) 逆变换得到衰减信息，即组织器官信息。

换句话说，\(radon\) 变换可用于三维重建，且是三维重建研究方向的helloworld。