TraditionCV_4: 直线拟合，找边，霍夫直线

presentation: 这算是理论知识吧，整理有一段时间了，里面很多东西来自抄袭， 希望各位见到可以留名。其实很明显，这就是去拟合直线的一些方法，可以去定位轮廓边。本文只是最基本的算法原理以及实现，在复杂场景下定位边，这些只是底层代码中的部分原理公式，需要更多的优化和其他增强方法。

最小二乘法：

几何意义：高维空间中的一个向量在低维子空间的投影

小二乘法是由勒让德在19世纪发现的，形式如下式：

观测值就是我们的多组样本，理论值就是我们的假设拟合函数。目标函数也就是在机器学习中常说的损失函数，我们的目标是得到使目标函数最小化时候的拟合函数的模型。

可以加入正则化项，降低过拟合。

回归问题中，损失函数是平方损失，正则化可以是参数向量的L2范数,也可以是L1范数。

L1: regularization*abs§
L2: 0.5 * regularization * np.square§

RANSAC:

Random Sample consensus: 随机采样一致算法

最小二乘法的问题

参考：https://www.cnblogs.com/washa/p/3164212.html

Homography矩阵:单应性矩阵

基本算法思想与流程：

1、选择出可以估计出模型的最小数据集（直线拟合：每次选出两个点；Homography矩阵：每次4个点）

2、用1选出的数据集计算出数据模型；（直线：（x1,y1),(x2, y2)求解a, b;

3、将所有数据带入这个模型（y = ax + b)，计算出内点的数量；（在一点误差范围内）

4、比较当前模型和之前推出的最好的模型的“内点”数量，更新最大“内点”数的模型参数（a, b) 和“内点”数。

5、重复1-4步，直到迭代结束或者当前模型达到要求（内点数目大于一定数量）

迭代次数推导：

iters = math.log(1 - P) / math.log(1 - pow(total_inlier / (SIZE * 2), 2))

推算：

假设内点占比：t = num_in / (num_in + num_out)

每次计算模型使用N个点（直线2个），选取的点至少有一个外点：1 - pow(t, N)

也就是说，在迭代k次的情况下，pow( 1- pow(t, N), k)就是k次迭代计算模型都至少采样到一个“外点”去计算模型的概率。那么能采样到正确的N个点去计算出正确模型的概率：

P = 1 - pow( 1- pow(t, N), k)

求得：

k = log(1 - P) / log(1 - pow(t, N))

内点”的概率 t 通常是一个先验值。然后 P 是我们希望RANSAC得到正确模型的概率。如果事先不知道 t 的值，可以使用自适应迭代次数的方法。也就是一开始设定一个无穷大的迭代次数，然后每次更新模型参数估计的时候，用当前的“内点”比值当成 t 来估算出迭代次数

霍夫变换：

霍夫变换常用于检测直线特征，经扩展后的霍夫变换也可以检测其他简单的图像结构。

已经进行初步的轮廓检测之后，才进行直线检测

参考：https://www.cnblogs.com/php-rearch/p/6760683.html

在霍夫变换中我们常用公式

ρ = xcosθ + ysinθ

表示直线，其中ρ是圆的半径（也可以理解为原点到直线的距离），θ是直线与水平线所成的角度（0~180°），确定了它们，也就确定一条直线了，和下图略有出入的是实际的原点定在图片左上角。

cv::HoughLines(midImage, lines, 1, CV_PI / 180, 150);  // 输入的时二值图像，输出vector向量

原理是对于输入的二值图像中的像素点（有值的），按照步长（参数三参数四对应rho和theta的步长）分别计算出每个点上的所有可能的直线。记录下每条直线经过的点数（即存在多个点计算出的直线有交集），按照阈值（参数五）筛选符合条件的图像。

算法原理：

一条直线可由两个点A=(X1,Y1)和B=(X2,Y2)确定(笛卡尔坐标)；

另一方面，y = kx + q 可以写成关于（k, q)的函数表达式；

对应的变换可以通过图形直观表示：

变换后的空间成为霍夫空间。即：笛卡尔坐标系中一条直线，对应霍夫空间的一个点。

反过来同样成立（霍夫空间的一条直线，对应笛卡尔坐标系的一个点）：

再来看看A、B两个点，对应霍夫空间的情形：

一步步来，再看一下三个点共线的情况：

可以看出如果笛卡尔坐标系的点共线，这些点在霍夫空间对应的直线交于一点：这也是必然，共线只有一种取值可能。

如果不止一条直线呢？再看看多个点的情况（有两条直线）：

其实（3，2）与（4，1）也可以组成直线，只不过它有两个点确定，而图中A、B两点是由三条直线汇成，这也是霍夫变换的后处理的基本方式：选择由尽可能多直线汇成的点。

看看，霍夫空间：选择由三条交汇直线确定的点（中间图），对应的笛卡尔坐标系的直线（右图）。

到这里问题似乎解决了，已经完成了霍夫变换的求解，但是如果像下图这种情况呢？

k=∞是不方便表示的，而且q怎么取值呢，这样不是办法。因此考虑将笛卡尔坐标系换为：极坐标表示。

在极坐标系下，其实是一样的：极坐标的点→霍夫空间的直线，只不过霍夫空间不再是[k,q]的参数，而是的参数，给出对比图：

是不是就一目了然了？

给出霍夫变换的算法步骤：

code:

void hough() {
    Mat souImg = imread("建筑.png");
    imshow("原始图片", souImg);
    Mat contour;
    Canny(souImg, contour, 50, 200);
    imshow("轮廓图片", contour);
    int H_row;
    if (contour.cols > contour.rows)
        H_row = contour.cols;
    else
        H_row = contour.rows;
    Mat H(3*H_row, 180, CV_8S, Scalar(0));
    std::cout << H_row << std::endl;
 
    float theta, rho;
    for (int i = 0; i < contour.rows; i++) {
        for (int j = 0; j < contour.cols; j++) {
            if (contour.at<uchar>(i, j) > 0) {
                for (theta = 0; theta < 180; ++theta) {
                    rho = floor(i*cos(theta*CV_PI / 180) + j*sin(theta*CV_PI / 180));
                    try {
                        H.at<uchar>(rho + H_row, theta) += 1;
                    }
                    catch (...) {
                        std::cout << i << j << rho << theta << std::endl;
                        return;
                    }
                }
            }
        }
    }
    imshow("H", H);
    waitKey(0);
}

其实本质上就是：

交点怎么求解呢？细化成坐标形式，取整后将交点对应的坐标进行累加，最后找到数值最大的点就是求解的
也就求解出了直线。

opencv:

1、标准霍夫变换

会计算图像中的每一个点，计算量比较大，另外它得到的是整一条线，并不知道原图中直线的端点。

# 标准霍夫变换
# 参数1：要检测的二值图（一般是阈值分割或边缘检测后的图）
# 参数2：距离 ρ 的精度，值越大，考虑越多的线
# 参数3：角度 θ 的精度，值越小，考虑越多的线
# 参数4：累加数阈值，值越小，考虑越多的线 计算(r,θ) 累加数，累加数超过一定值后就认为在同一直线上（有一个阈值
lines = cv2.HoughLines(edges, 0.8, np.pi/180, 90) #检测出来的是极坐标（rho, theta)

2、统计概率霍夫直线变换

# HoughLinesP直接给出了直线的断点，在画出线段的时候可以偷懒
# minLineLength：最短长度阈值，比这个长度短的线会被排除
# maxLineGap：同一直线两点之间的最大距离
lines = cv2.HoughLinesP(edges, 0.8, np.pi/180, 150, minLineLength=20, maxLineGap=10)

Code:

import cv2
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
from scipy.optimize import leastsq
import random
import math

__all__ = ['HoughLine', 'LeastLine', 'RansacLine']


class HoughLine:
    def __init__(self):
        pass

    def fitline_houghlines(self, edges):
        """

        :param edges:
        :return: lines: class 'numpy.ndarray': [[[rho, theta]] [[..]] .. [[..]]]
        """
        # 标准霍夫变换
        # 参数1：要检测的二值图（一般是阈值分割或边缘检测后的图）
        # 参数2：距离 ρ 的精度，值越大，考虑越多的线
        # 参数3：角度 θ 的精度，值越小，考虑越多的线
        # 参数4：累加数阈值，值越小，考虑越多的线 计算(r,θ) 累加数，累加数超过一定值后就认为在同一直线上（有一个阈值
        # 检测出来的是极坐标（rho, theta)
        lines = cv2.HoughLines(edges, 0.8, np.pi/180, 90)
        return lines

    def fitline_houghlinesP(self, edges):
        """

        :param edges:
        :return: lines: class 'numpy.ndarray': [[[x1 y1 x2 y2]] [[....]]..[[....]]]
        """
        # 概率霍夫直线变换
        # HoughLinesP直接给出了直线的断点，在画出线段的时候可以偷懒
        # minLineLength：最短长度阈值，比这个长度短的线会被排除
        # maxLineGap：同一直线两点之间的最大距离
        lines = cv2.HoughLinesP(edges, 0.8, np.pi/180, 150, minLineLength=20, maxLineGap=10)
        return lines

    def drawline_fitline_houghlines(self, lines, drawing):
        """

        :param lines: lines: class 'numpy.ndarray': [[[rho, theta]] [[..]] .. [[..]]]
        :param drawing:
        :return:
        """
        for line in lines:
            rho, theta = line[0]
            a = np.cos(theta)
            b = np.sin(theta)
            x0 = a * rho
            y0 = b * rho
            x1 = int(x0 + 1000 * (-b))
            y1 = int(y0 + 1000 * a)
            x2 = int(x0 - 1000 * (-b))
            y2 = int(y0 - 1000 * a)
            drawing = cv2.line(drawing, (x1, y1), (x2, y2), (255, 0, 0), 1)
        return drawing

    def drawline_fitline_houghlinesP(self, lines, drawing):
        """

        :param lines: lines: class 'numpy.ndarray': [[[x1 y1 x2 y2]] [[....]]..[[....]]]
        :param drawing:
        :return:
        """
        for line in lines:
            line = line[0]
            x1 = line[0]
            y1 = line[1]
            x2 = line[2]
            y2 = line[3]
            drawing = cv2.line(drawing, (x1, y1), (x2, y2), (255, 0, 0), 1)
        return drawing

# 目标函数y=sin(2πx) , 加上一个正太分布的噪音干扰，用多项式去拟合
def real_func(x):
    return np.sin(2 * np.pi * x)

# 多项式函数(拟合函数，也就是h(x))
# ps: numpy.poly1d([1,2,3])  生成  $1x^2+2x^1+3x^0$*
def fit_func(p, x):
    f = np.poly1d(p)
    return f(x)

# 残差函数
def residuals_func(p, x, y):
    ret = fit_func(p, x) - y
    return ret

# 添加正则项的残差函数
def residual_func_regularization(p, x, y):
    ret = fit_func(p, x) - y
    ret = np.append(ret, np.sqrt(0.5 * regularization * np.square(p)))
    #print('ret', ret)
    return ret


class LeastLine:
    def __init__(self):
        pass

    def fitline_least_square(self, M):
        """
        使用最小二乘法拟合y=sin(2πx)
        :param M: 多项式的次数
        :return: list 拟合参数
        """
        # 随机初始化多项式参数
        # 生成p+1个随机数的列表，这样poly1d函数返回的多项式次数就是p
        p_init = np.random.rand(M + 1)
        #print('p_init:', p_init)
        # 最小二乘法
        # leastsq（）函数可以很快速地使用最小二乘法对数据进行拟合
        # 三个参数：误差函数、函数参数列表、数据点
        # p_lsp = leastsq(residuals_func, p_init, args=(x, y))
        p_lsp = leastsq(residual_func_regularization, p_init, args=(x, y))
        #print('Fitting Parameters:', p_lsp[0])

        # 可视化
        plt.plot(x_points, real_func(x_points), label='real')
        plt.plot(x_points, fit_func(p_lsp[0], x_points), label='fitted curve')
        plt.plot(x, y, 'bo', label='noise')
        plt.legend()
        plt.show()
        return p_lsp

class RansacLine:
    def __init__(self):
        pass

    def fitline_ransac(self, iters, sigma, SIZE, RANDOM_X, RANDOM_Y, P):
        # 最好模型的参数估计和内点数目
        best_a = 0
        best_b = 0
        pretotal = 0
        while iters > 0:
            # 随机在数据中红选出两个点去求解模型
            sample_index = random.sample(range(SIZE * 2), 2)
            x_1 = RANDOM_X[sample_index[0]]
            x_2 = RANDOM_X[sample_index[1]]
            y_1 = RANDOM_Y[sample_index[0]]
            y_2 = RANDOM_Y[sample_index[1]]

            # y = ax + b 求解出a，b
            a = (y_2 - y_1) / (x_2 - x_1)
            b = y_1 - a * x_1

            # 算出内点数目
            total_inlier = 0
            for index in range(SIZE * 2):
                y_estimate = a * RANDOM_X[index] + b
                if abs(y_estimate - RANDOM_Y[index]) < sigma:
                    total_inlier = total_inlier + 1

            # 判断当前的模型是否比之前估算的模型好
            if total_inlier > pretotal:
                iters = int(math.log(1 - P) / math.log(1 - pow(total_inlier / (SIZE * 2), 2)))
                #print(iters)
                pretotal = total_inlier
                best_a = a
                best_b = b

            # 判断是否当前模型已经符合超过一半的点
            #print("t:", total_inlier)
            if total_inlier > SIZE / 2:
                break
        return best_a, best_b


if __name__ == "__main__":
    img = cv2.imread('test.jpg')
    # 已经进行初步的轮廓检测之后，才进行直线检测

    # Hough找直线
    drawing = np.zeros(img.shape[:], dtype=np.uint8)
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    edges = cv2.Canny(gray, 50, 150)
    hough = HoughLine()
    lines = hough.fitline_houghlinesP(edges)
    drawing = hough.drawline_fitline_houghlinesP(lines, drawing)


    plt.subplot(231), plt.imshow(img), plt.title('img')
    plt.subplot(232), plt.imshow(drawing), plt.title('drawing')
    plt.subplot(233), plt.imshow(gray), plt.title('gray')
    plt.subplot(234), plt.imshow(edges), plt.title('edges')
    plt.show()

    # 使用最小二乘法拟合y=sin(2πx)
    ll = LeastLine()
    # 加正则化项 正则化系数lambda不能过高，过高则导致欠拟合现象即此时的惩罚项权重较高
    regularization = 0.0001
    # 10个点 随机选取0-1之间的9个数作为x
    x = np.linspace(0, 1, 10)
    # 画图时需要的连续点
    x_points = np.linspace(0, 1, 1000)
    # 目标函数
    y_ = real_func(x)
    # 加上正态分布噪音的目标函数的值
    y = [np.random.normal(0, 0.1) + y1 for y1 in y_]
    p_lsp_0 = ll.fitline_least_square(M=9)

    # Random Sample consensus: 随机采样一致算法
    rc = RansacLine()
    # 数据量。
    SIZE = 50
    # 产生数据。np.linspace 返回一个一维数组，SIZE指定数组长度。
    # 数组最小值是0，最大值是10。所有元素间隔相等。
    X = np.linspace(0, 10, SIZE)
    #print(X)
    Y = 3 * X + 10

    fig = plt.figure()
    # 画图区域分成1行1列。选择第一块区域。
    ax1 = fig.add_subplot(1, 1, 1)
    # 标题
    ax1.set_title("RANSAC")

    # 让散点图的数据更加随机并且添加一些噪声。
    random_x = []
    random_y = []
    # 添加直线随机噪声
    for i in range(SIZE):
        random_x.append(X[i] + random.uniform(-0.5, 0.5))
        random_y.append(Y[i] + random.uniform(-0.5, 0.5))
    # 添加随机噪声
    for i in range(SIZE):
        random_x.append(random.uniform(0, 10))
        random_y.append(random.uniform(10, 40))
    RANDOM_X = np.array(random_x)  # 散点图的横轴。
    RANDOM_Y = np.array(random_y)  # 散点图的纵轴。

    # 画散点图。
    ax1.scatter(RANDOM_X, RANDOM_Y)
    # 横轴名称。
    ax1.set_xlabel("x")
    # 纵轴名称。
    ax1.set_ylabel("y")

    # 使用RANSAC算法估算模型
    # 迭代最大次数，每次得到更好的估计会优化iters的数值
    iters = 100000
    # 数据和模型之间可接受的差值
    sigma = 0.75
    # 希望的得到正确模型的概率
    P = 0.99
    best_a, best_b = rc.fitline_ransac(iters, sigma, SIZE, RANDOM_X, RANDOM_Y, P)
    Y = best_a * RANDOM_X + best_b
    # 直线图
    ax1.plot(RANDOM_X, Y)
    text = "best_a = " + str(best_a) + "\nbest_b = " + str(best_b)
    plt.text(5, 10, text,
             fontdict={'size': 8, 'color': 'r'})
    plt.show()
ax1.set_xlabel("x")
    # 纵轴名称。
    ax1.set_ylabel("y")

    # 使用RANSAC算法估算模型
    # 迭代最大次数，每次得到更好的估计会优化iters的数值
    iters = 100000
    # 数据和模型之间可接受的差值
    sigma = 0.75
    # 希望的得到正确模型的概率
    P = 0.99
    best_a, best_b = rc.fitline_ransac(iters, sigma, SIZE, RANDOM_X, RANDOM_Y, P)
    Y = best_a * RANDOM_X + best_b
    # 直线图
    ax1.plot(RANDOM_X, Y)
    text = "best_a = " + str(best_a) + "\nbest_b = " + str(best_b)
    plt.text(5, 10, text,
             fontdict={'size': 8, 'color': 'r'})
    plt.show()