一 Hough变换（霍夫变化）

1.1 什么是Hough变换

霍夫变换是图像处理中的一种特征提取技术，通过投票算法检测识别具有特定几何形状的物体。我们这篇文章，主要对直线进行检测。

1.2 Hough变换原理

我们看见直角坐标系中有两个点A B，他们构成了一条直线y=kx=b。知道了斜率k和截距b就可以确定一条直线。

我们可以将这个坐标系看成一个笛卡尔空间。往深处思考一下，其实x就是我们的输入，k和b就是参数，y就是输出。就是当我们知道一个x，就一定知道直线上对于的点。

接下来，我们就要将笛卡尔空间转化为霍夫空间。其实就是需要将直角坐标系上的点转化为极坐标上的点。

那现在有了点，我们需要处理直线了。

 当我们知道和就可以确定这一条直线了，由于我们有点P的坐标，通过坐标就知道了和的数值。​​​​​​我们将左边的式子变化为右边的式子。此时我们将右式转换为下式。

 此时，我们给定和时，知道x就可以求出y。与之前的y=kx+b相同，都是两个参数确定一条直线。

也就是说在极坐标下，一个点代表一条直线。

1.2 Hough变换原理

在一张图中，只有边缘才可能是直线，因此，在进行Hough变换之前，需要进行边缘检测。

从这幅图我们可以很清晰看出，在经过边缘检测后，图中所有的边缘都检测出来了，但是很明显并不是所有的边缘都是我们需要的，因此需要对需要的边缘进行提取。 其实最好的思路就是，有多少的像素点可以贡献给这条直线。 

ok，知道思路后，我们先对一个像素点分析，图中一个像素点可以有无数条直线过这一点。

这些直线可以可以在极坐标坐标系中进行表示出来。 

那在图中我们想得到的直线上不同的两个点，都可以在极坐标上用不同的曲线表示出。

我们发现直角坐标上的点转化到极坐标中，可以发现转化到极坐标中的曲线相交在某一个点。

结合之前说到的”有多少的像素点可以贡献给这条直线“，所以我们就可以根据极坐标上相交多的一个点来确定这些曲线在直角坐标系中其实是在一条直线上。

 那此时，可能有很多的边缘其实也相交的很多，但是这些并不是我需要的，此时就需要对其阈值进行设置。

1.3 Hough变换整体流程

①直线检测（边缘检测图像，阈值）阈值其实就是直线上像素点的数量

②建立一个极坐标空间

③遍历边缘图像中的每一个边缘点（通过边缘图像坐标转到极坐标）

④便利极坐标中的

⑤根据求出，在矩阵上计数

⑥返回技术超过阈值的和，就得到直线的函数了

1.4 实战代码

import cv2
import numpy as np


def hough(img, threshold=50):  # 有150个点是在一条直线上，就默认是需要的直线
    thetas = np.deg2rad(np.arange(0, 180, 2))  # 0开始到180（不包括180）步长为2的整数序列，随后将角度转为弧度
    row, cols = img.shape
    diag_len = np.ceil(np.sqrt(row ** 2 + cols ** 2))  # 对角线长度
    rhos = np.linspace(-diag_len, diag_len, int(2 * diag_len))
    cos_t = np.cos(thetas)
    sin_t = np.sin(thetas)
    num_theta = len(thetas)
    # 构造计数矩阵
    vote = np.zeros((int(2 * diag_len), num_theta), dtype=np.uint64)
    y_int, x_int = np.nonzero(img)  # 把白色的点（边缘点）拿出来就可以
    # 计数
    for i in range(len(x_int)):
        x = x_int[i]
        y = y_int[i]
        for j in range(num_theta):
            rho = round(x * cos_t[j] + y * sin_t[j]) + diag_len
            vote[int(rho), j] += 1

    # 拿到计数大于阈值的theta和rho
    indeies = np.argwhere(vote > threshold)
    rhos_idx = indeies[:, 0]
    theta_idx = indeies[:, 1]

    return vote, rhos[rhos_idx], thetas[theta_idx]


img = cv2.imread("D:/desk/road.png")
cv2.imshow("src", img)  # 原图

gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
cv2.imshow('gray', gray)  # 灰度图

edges = cv2.Canny(gray, 320, 360)  # 140，70是maxval与minval
cv2.imshow('edge', edges)  # 边缘图像

vote, rhos, thetas = hough(edges)
vote = np.uint8(vote.T)
cv2.imshow("vote", vote)

for rho, theta in zip(rhos, thetas):
    x_center = img.shape[1] / 2
    x1 = int(x_center + 500)
    x2 = int(x_center - 500)
    a = np.sin(theta)
    b = np.cos(theta)
    if a == 0.0:
        a = 1e-5
    if b == 0.0:
        b = 1e-5
    y1 = int((rho - x1 * b) / a)
    y2 = int((rho - x2 * b) / a)
    cv2.line(img, (x1, y1), (x2, y2), (0, 255, 0), 1)

cv2.imshow("hough", img)
cv2.waitKey(0)

 

 

二 Radon变换

2.1什么是radon变换

radon变换是将二维平面内的函数沿不同方向的直线进行积分投影的过程。

2.2 radon变换原理（可以简单了解原理，如果想直入主题可以直接看代码分析）

我们都知道医院的CT扫描，CT扫描的结果将我们的身体进行分层，每一层都会以一张图像显示。将这些切片图像结合在一起就得到了人体的三维结构。

那为什么可以得到每一层的切片呢？

现在我们将这三个球看成一个整体，接下来我们对这个整体进行扫描。

我们用一束x光照射它， 红色的线是平行的x光束，蓝线是x光接受器。

x光穿过球体时，有一部分会被吸收，吸收的多少和穿透长度成正比，以此光束的穿透长度越长，说明被吸收的光强越多，这样接收器上的光强就反应了该方向上物体的厚度了。

基于此，选转光束，就得到不同方向的物体厚度。

当然，每条射线都得到路径上的物体厚度。基于此，我们就可以完成一张图像上的Radon变换了。

如果想了解更多关于立体的Radon变换可以看：Radon变换原理可视化与讲解——给御坂美琴做CT扫描_哔哩哔哩_bilibili

2.3 Radon变化代码分析

其实Radon变化的代码实现比较简单，但是关于其返回值和参数的相关内容比较抽象，我也是通过实验，分析和deepseek才最终了解了其具体内容。

我们先看代码:

from skimage.transform import radon
sinogram = radon(edge,theta=theta_range) #此时theta_range=np.arange(-10,10,1)

我们可以看出代码是很简单的，我们先分析一下参数：

radon(image, theta=None, circle=True, *, preserve_range=False):

image（必填）：输入的二值图像，支持uint8或者float类型

theta：投影角度数组（单位：度）。默认值为0到180之间的180个等间距角度（即theta=np.arange(180)）。

我们的例子中设theta是从-10°到10°，1为步长，每次转1°。

circle：布尔值，默认为True，只处理图像内切圆内的像素（避免角落噪声的干扰）；若设置为False，处理整个图像区域。

preserve_range：是否保持输入的数据范围，若为Flase，图像会进行归一化处理。

目前看起来就是正常的一些参数。

接着我们看其返回值：

其返回值是一个二维数组，形状为(len(num_angles), num_shifts)

len(num_angles):表示投影角度的总数。例如，若theta_range包含180个角度值，则第一维长度为180，每一行对应一个角度的投影数据.

num_shifts:指每个投影角度下的位移步数（即平行光束的平移次数）。在Radon变换中，这通常等于投影方向上的探测器单元数量或图像对角线长度的整数倍。例如256x256图像，num_shifts可能为256或者更高的数值，以保证可以覆盖整个图像。（记住此时提到的位移步数）

每一行对应theta中的一个角度，每一列表示投影位移（投影方向的偏移量）

[[43.56481161 43.46496013 40.9024631  ... 42.50609515 49.82001402
  40.00432686]
 [46.15971292 40.8922127  53.23695708 ... 39.44164372 52.8321457
  45.9958921 ]
 [52.00749004 51.62998355 47.27055231 ... 40.27544849 36.28792505
  31.98761969]
 ...
 [46.6219134  47.59189496 48.64718933 ... 38.64512074 40.03041139
  42.70812504]
 [46.59756545 47.56163714 47.46676294 ... 37.23127786 46.38930346
  52.04752387]
 [39.24030297 36.64408011 43.19367528 ... 38.29984588 37.34519734
  68.05730034]]

这是我对我处理的一张图片，print（sinogram）出的结果。

比如数组的第一行，其实代表的就是在-10°的情况下，每一个位移步数下，垂直投影线覆盖的像素积分值。（一定要区分开，位移步长和像素积分值（我在很多地方都没有看见把这两个同时拿出来进行对比的，导致我在打印出结果时非常困惑））如果关于这个部分还是有点模糊，看下一节的代码分析补充。

默认现在大家都知道上述内容了，那它是如何对直线进行检测，并且矫正的呢？

此时我们假设有一条45°的白色边缘，如果我们在0°投影的话，投影线只能覆盖边缘的局部（像素积分值就很低），但是如果我们在45°进行投影，投影与边缘完全对其，积分值达到最大，那我们就知道直线其实就是在积分值最大的角度，于是我们进行以下操作，找到积分值最大的角度，利用rotate函数进行旋转就行。

max_index = np.argmax(np.sum(sinogram**2,axis=0))
print(max_index)
return theta_range[max_index]

此时我们需要注意max_index并不是我们旋转的角度，它是sinogram的行序列，从序号0开始。

由于我的是（-10，10，1）我print的结果是8（序号），由此其对应的角度就是-2°，返回值结果就是-2。

有了角度使用ndimage.rotate（）进行角度旋转。

2.4 代码分析补充

如果上面的返回值有些困惑，接下来我会举一个例子。

假设输入图像 edge 是一个 ​3×3像素 的边缘图像（白色边缘，黑色背景），theta_range 设置为 [0°, 45°]，则可能的输出如下：

[[ 1.  0.  2. ]    # theta=0° 时的投影值（不同位移步下的积分）
 [ 0.  3.  0. ]]   # theta=45° 时的投影值

第一行就是0°情况下的所有投影值，第二行就是45°情况下的所有投影值.

最开始提到返回值的形状 此时的形状就是（2，3）.2对应theta_range达到两个角度，3代表每个角度的位移步数。其中第一行的1.0,0.0,2.0代表的就是每个角度下不同位移步长的像素积分值。

后续还会继续更新其他的直线检测的方法，LSD,FLD等