直线段检测算法LSD：a Line Segment Detector

直线段检测算法LSD：a Line Segment Detector

微信公众号：幼儿园的学霸

目录

简介

影像中的直线特征是视觉感知的重要线索和解释图像的基本依据,相比于点特征，直线段特征对视角和光照变化具有很好的鲁棒性。常用的直线检测有Hough直线检测和LSD等.Hough直线检测将直线检测映射为参数空间中点的检测，而LSD通过查找近似矩形区域来获得直线。LSD是由Rafael Grompone von Gioi等人在2012年提出的一种局部直线快速提取的算法, 它能在线性时间内得出亚像素级精度的检测结果,LSD算法较为稳定且存在错误控制, 可以在无人工调节任何参数的情况下快速得到满意的线特征提取结果.

流程

LSD的核心是像素合并与误差控制。

整体流程大致可以分为如下几个步骤：
1、将输入图像通过高斯核采样，缩放到0.8;
2、计算图像梯度强度，以及每个像素点的梯度角：$atan2(gx,-gy)$，其中$gx、gy$分别为水平和垂直方向梯度;
3、设置图像梯度强度范围到[0, 1023] (大于1023的梯度强度，强制设置为1023)；
4、创建1024个链表，遍历整个梯度图，根据梯度强度，相同梯度值像素坐标放入同一张链表中；
5、将1024个链表，按从大到小顺序，合成一张大链表(首部为1023链表，尾部为0)；
6、取出链表头部存储图像坐标位置，作为种子像素，根据梯度角方向相似，进行区域扩散。(每扩散一个像素，将该像素坐标从链表中删除，并且做标记，之后新的区域扩散，无法再扩散到该像素)；
7、将扩散区域进行矩形拟合；
8、检测拟合矩形内梯度角，根据对应NFA公式(来自论文)，计算出该拟合矩形精度误差，如果不满足，则丢弃；如果满足要求，则将该矩形记录存储，表示为一条检测到的直线；
9、继续回到第6步，从链表中，找到下一个种子点，从剩下图像进行区域扩散，至到遍历完全图，得到所有检测到的直线。

详细步骤

1.尺度缩放:

按照原文，这一步可选择跳过（s =1.0），算法的默认选择是s=0.8。目的是为了消除数字离散图像的锯齿效应。缩放完毕之后用高斯下采样的方式对输入图像进行操作，计算时，先进行高斯滤波，然后进行下采样。
其中高斯核的标准差为$\sigma =\zeta /s$,其中$\zeta =0.6$以此来保持锯齿效应和图像模糊之间的平衡。高斯滤波核的半径为$h=ceil(\sigma \ast \sqrt{2\ast 3\ast ln(10)})$,那么高斯滤波核的窗口大小为$w=1+2\ast h$.
举例说明，如果$s=0.8$,那么$\sigma =0.75,h=ceil(2.787)=3,w=7$

2.梯度和方向计算：

梯度计算中论文用到的是2*2的模板。主要是为了用contrario model的时候保证相邻点方向分布的独立性,减少梯度计算过程中的彼此依赖,同时小模板时计算相对较快.

$$\begin{array}{r}{G}_x=\frac{1}{2}\left[\begin{array}{cc}-1& 1\\ -1& 1\end{array}\right]G_y=\frac{1}{2}\left[\begin{array}{cc}-1& -1\\ 1& 1\end{array}\right]\end{array}\text{(1)}$$

$$\begin{array}{r}{G}_x(x,y)=\frac{1}{2}[f(x+1,y)-f(x,y)+f(x+1,y+1)-f(x,y+1)]\\ G_y(x,y)=\frac{1}{2}[f(x,y+1)-f(x,y)+f(x+1,y+1)-f(x+1,y)]\end{array}\text{(2)}$$

代码中使用的梯度幅值和梯度幅角求解公式为:
$$\begin{array}{rl}& G=\sqrt{G_x^2+G_y^2}\\ & \theta =atan2(Gx,-Gy)\end{array}\text{(3)}$$

这种方法计算的梯度其实是$（x+0.5，y+0.5）$处的，在最后的步骤会对这里的计算偏差进行补偿。同样的，由于计算的梯度是$（x+0.5，y+0.5）$，代码中由于未对图像进行边界扩充，因此图像的最后一行及最后一列未进行梯度方向和幅值的计算。
需要注意的是，从黑到白跟从白到黑是不同的，两者在梯度上具有180°的差异。这意味着，LSD的直线分割结果是有方向的，并非起始点和终止点任意的。举例来说，当灰度图进行亮度翻转时（黑的像素变白，白的像素变黑），那么所有LSD检测到的直线是不变的，但每根直线的起始点和终止点都交换了。

3.梯度伪排序(Pseudo-Ordering)：

LSD是一种贪婪算法,因此像素处理的顺序将对结果有一定的影响. 具有较高梯度幅值的像素点所在的区域通常具有较强的边缘,(而在边缘中,中间的像素一般具有最高的梯度幅值).因此,梯度值越大，越是显著的边缘点，更适合作为种子点,从其入手进行直线分割检测.

但是对梯度值进行完全排序是一个时耗性很高的工作。因此简单的将梯度值划分为1024个等级（bins），这1024个等级涵盖了梯度由0~255的变化范围，这种排序是一个线性的时耗。种子点从梯度值最高的bin开始搜索，依次往下，直至所有点标记为USED。
BINS划分原理:设求解的图像中每个点的梯度幅值为$modgra{d}_{ij}$,其中的梯度幅值最大值为$max\_grad$，需要划分为$bins$个等级，那么求解等级系数为$bin\_coef=\frac{bins-1}{max\_grad}$,因此每个点量化后的梯度幅值为$\widehat{modgra{d}_{ij}}=bin\_coef\ast modgra{d}_{ij}$

小梯度值抑制：
小梯度值点往往出现在平滑区域，或者仅仅是噪声。不在关注的范围内，但是他们的存在往往会严重影响直线角度的计算。考虑到区域增长时有角度承受范围$[-\tau ,\tau ]$。因此小梯度“脏点”的最大影响角度需要小于$\tau$。
一个像素点梯度幅值为g与某方向一条线段夹角为$\alpha$，则$g\ast sin(\alpha )$代表像素点在与该方向垂直的方向梯度变化量$q_0>=1$，才代表该点像素拥有的在该直线方向容忍度范围内level-line可以被后续区域生长接受。
$$\rho =\frac{q}{sin\tau }$$
其中,$\tau =22.5$为论文中区域增长的梯度方向角容忍度值，$q=2$为经验值，因此梯度阈值为$\rho =2/sin(22.5°)=5.226$。

在图像中，直线通常由一系列具有相似梯度方向的像素组成。这意味着直线上的像素在梯度方向上会有较小的差异。
当直线方向与梯度方向接近时，即夹角较小，直线上的像素与梯度方向更为一致。在这种情况下，我们希望保留更多的梯度值，因为这些梯度值可能代表直线上的有效信息。因此，梯度阈值会较小，以保留这些梯度值。
相反，当直线方向与梯度方向差异较大时，即夹角较大，直线上的像素与梯度方向差异较大。在这种情况下，我们认为这些梯度值可能是由噪声或非直线特征引起的。因此，为了过滤掉这些不相关的梯度值，梯度阈值会较大。
通过根据直线的几何特征自适应地调整梯度阈值，LSD直线检测算法可以提高对直线的检测灵敏度和准确性。在直线方向与梯度方向更接近的情况下，保留较小的梯度阈值可以捕捉到细微的直线特征。而在直线方向与梯度方向差异较大的情况下，采用较大的梯度阈值可以过滤掉噪声和非直线特征，从而提高直线检测的鲁棒性。
总之，通过使用梯度方向和角度来计算梯度阈值，LSD直线检测算法能够根据直线的几何特征自适应地调整阈值，并在直线和非直线像素之间进行有效的区分，从而提高直线检测的准确性和鲁棒性。

4.直线（矩形）区域增长:

该步骤的目的是利用区域生长算法合并场里方向近似一致的像素，得到一系列线支撑域.
从被标记为UNUSED的有序像素中取一个像素作为种子点，由种子点搜索角度满足并且状态为UNUSED的点（八邻域）形成的区域称为line-support region，递归地计算区域角度region_angle，如果搜索的点角度$angle$和region_angle的差值在角度承受范围$\tau$内(即:$|angle- region_angle|\leqslant \tau $)，则将点加入到这个区域。每次增加一个点，像素点加入到该区域后，重新按照下面公式计算区域方向角，同时其将被标记为USED。

reigon_angle会更新：

$$region\_angle=arctan(\frac{\sum sin(angl{e}_i)}{\sum cos(angl{e}_i)})$$
其中,i为region内所有的点.该公式可以按照计算区域的梯度方向的角度理解，根据矢量相加：总梯度=梯度1+梯度2+…,那么上式实际就是根据梯度的分量从而求解最终的梯度.

直线支撑区域生长过程如下图所示。

5.矩形拟合

为方便后续计算，将得到的直线支撑区域用拟合矩形覆盖，矩形示意图如下所示。

在LSD算法中，每一个直线支撑区域对应一个拟合矩形区域，每一个矩形区域对应一条直线。矩形区域参数求解如下。
通过将线直线支撑区域看做一个实体，其内每个像素点的梯度值作为像素点的质量，整个实体的质心作为最小外接矩形的中心点。记线支持区域为$r$，其内像素点$j$的坐标为$(x_j,y_j)$，梯度值为$G(x_j,y_j)$，则矩形的中心点$(c_x,c_y)$为
$$c_x=\frac{{\sum }_{j\in r}G\left(x_j,y_j\right)\cdot x_j}{{\sum }_{j\in r}G\left(x_j,y_j\right)},c_y=\frac{{\sum }_{j\in r}G\left(x_j,y_j\right)\cdot y_j}{{\sum }_{j\in r}G\left(x_j,y_j\right)}$$

上面公式可以理解为，利用梯度值作为权重，对位置进行加权处理

确定矩形中心后，确定矩形的角度。矩形的主惯性轴方向为矩阵M的最小特征值对应的特征向量的角度
矩形的参数则通过以下公式求解.
$$\mathbf{M}=\left[\begin{array}{ll}{m}^{xx}& m^{xy}\\ m^{xy}& m^{yy}\end{array}\right]$$
其中,
$$\begin{array}{c}{m}^{xx}=\frac{{\sum }_{j\in r}G\left(x_j,y_j\right)\cdot {\left(x_j-c_x\right)}^2}{{\sum }_{j\in r}G\left(x_j,y_j\right)},\\ m^{yy}=\frac{{\sum }_{j\in r}G\left(x_j,y_j\right)\cdot {\left(y_j-c_y\right)}^2}{{\sum }_{j\in r}G\left(x_j,y_j\right)},\\ m^{xy}=\frac{{\sum }_{j\in r}G\left(x_j,y_j\right)\cdot \left(x_j-c_x\right)\cdot \left(y_j-c_y\right)}{{\sum }_{j\in r}G\left(x_j,y_j\right)}\end{array}$$
记特征值λ和2维列向量$x=(x_1,x_2{)}^T$，使:
$$(\mathbf{M}-\lambda \mathbf{E})\mathbf{x}=0$$
当$\lambda ={\lambda }_{min}|$时，解上式可得：
$$\{\begin{array}{ll}\frac{x_2}{x_1}=\frac{{\lambda }_{\min }-m^{xx}}{m^{xy}}& \left(\left|m^{xx}\right|>\left|m^{yy}\right|\right)\\ \frac{x_2}{x_1}=\frac{m^{xy}}{{\lambda }_{\min }-m^{yy}}& \left(\left|m^{xx}\right|⩽\left|m^{yy}\right|\right)\end{array}$$
进而得到矩形的主惯性轴方向角$\theta =arctan(x_2/x_1)$。
其中${\lambda }_{min}=0.5\ast (m^{xx}+m^{yy}-\sqrt{(m^{xx}-m^{yy}{)}^2+{(2\ast m^{xy})}^2})$

上面部分也可以通过求解协方差矩阵来求解特征向量来代替，得到特征向量后，求特征向量与x轴的夹角，这样理解会比较简单点,该尝试在代码中有添加，可以参考代码部分的求解进行理解。

确定了矩形的中心和矩形的朝向后，包含区域所有点的最小矩形即为逼近得到的矩形。

该部分应为求解最小外接矩形（旋转矩形）,但代码中的公式来源未能找到、明确理解.

6.类内点密度/矩形改进(改善)：

有时候，容忍角度$\tau$的处理方法会产生错误的判断，如下图所示。当两条直线边相交形成了一个夹角，并且夹角小于$\tau$时，就会出现一个包含两条直线的直线支撑域和其对应的矩形，应该将这个矩形分割成这两条直线对应的两个矩形。根据直线支撑域和它对应的矩形是非常匹配的，对应的aligned points 密度非常高，LSD算法采用aligned points密度来检测这个角度问题.如果像素点的角度与矩形主方向角度的夹角在容忍角度$\tau$内，就将该像素点称之为aligned point.
aligned points密度为aligned points数目和矩形面积的比值,定义如下:
$d=k/n$。k为类内点个数，n为R的length*width。设同性点密度阈值为D,若d > D。接受这个矩形。否则，需要将R截断。在这里设置D=0.7。文章认为这个数字既能保证同一个R中的类内点属性相近，也能保证R不会被过分的分割为小的矩形。

R截断的方法有两种：1）缩小角度误差阈值；2）缩小区域半径。在这两种方法中，区域中的一部分像素被保留，而另一些则重新标记为NOT USED，因此它们可以在将来的区域生长中再次使用。

缩小角度容忍阈值：简单的将τ值缩小，再次从当前R的seed开始搜寻，看是否满足要求。
缩小区域半径：逐渐去除远离种子点的像素（点到种子点的距离大于区域半径,需要注意的是，该区域半径会不断进行缩小），直到区域密度满足标准或区域太小而被拒绝为止。当线支持区域对应于一条曲线时，该方法效果最好，并且在满足密度准则之前需要减少区域，通常意味着对曲线有一定程度的近似。区域半径缩小时，每次将半径缩小75%。

具体实现时,
算法会选择第1个点作为种子点，并计算该点的角度(ang_c)，然后会遍历区域R内的所有点，并判断距离种子点的距离是否小于矩形的宽度，如果满足条件，则计算该点的角度与种子点角度(ang_c)之间的差异(ang_d)，并将角度差异累加，同时统计满足条件的点的数量。计算完成后，算法会计算平均角度(mean_angle)和标准差的2倍，这里2倍的标准差将作为重新确定区域的角度容差。
接下来，算法会以种子点重新生长成新的区域（生长过程中的角度容差为上面的2倍标准差），如果生成的新区域点数量小于2，说明无法形成直线段，因此该区域被舍弃；如果生成的新的区域点数量大于等于2，则进行矩形拟合，同时计算矩形区域的密度，如果密度满足要求，则表示该矩形区域细化处理完成；否则，则尝试通过缩小区域半径的方法来调整该区域。

截止到步骤6，可以认为线段已经提取完毕，但是在论文中，又通过NFA公式对矩形进行修正，改进矩形的拟合精度

7.矩形修正(Rectangle Improvement)

该部分未能理解，公式等未搜到相关解释,仅将能够看到的解释附上，查看LSD算法的路上最不容易理解的点也是这里，仅主观上能够接受.

如果当前的R仍旧不能满足NFA，以下的方法将对其进行改进。考虑到在有些情况下，删除line-support region中的一个点会减少R的 length-1个点（想象为对角线）。对不满足NFA的R，采取以下迭代策略：
1.反复减小p=p/2
2.短边反复减少一行
3.长边反复减少一行
4.另一长边反复减少一行
5.反复减小p=p/2 直至满足NFA或R过小被拒或p为原来的1/32。

the Helmholtz principle：在完美噪声图像图像中不应该检测到目标。
contrario approach：一个不会检测到目标的噪声图像。对于本课题，contrario model是一个像素值为level-line angle的图像。其level-line angle随机分解为独立且服从平均分布于[0-2*PI]。
p-对齐点对于矩形判定很重要。这个概念是指设定一个阈值p，计算该矩形区域中的像素点，将矩形方向的夹角与该阈值比较，并将小于等于p的点进行统计。设其数量为k，通常p的值取为τ/π

NFA计算:
NFA stands for Number of False Alarms:
$$\mathrm{NFA}=NT\cdot B(n,k,p)$$

• NT - number of tests
• B(n,k,p) - p-对齐点所满足的公式。 tail of binomial distribution with parameters n,k and p:

而$ NT= (NM)^{5/2}$，其中N表示图像列数，M表示图像行。

$$B(n,k,p)=\sum _{j=k}^n\left(\begin{array}{c}n\\ j\end{array}\right)p^j(1-p{)}^{n-j}$$

$$\left(\begin{array}{c}n\\ k\end{array}\right)=\frac{\Gamma (n+1)}{\Gamma (k+1)\cdot \Gamma (n-k+1)}$$

最后得到的矩形区域，顶点坐标值需要额外加0.5，以弥补梯度计算时的误差;然后将顶点坐标放大回原图尺寸.

2个gamma公式(在计算NFA值时有用到):
$$\Gamma (x)=\frac{{\sum }_{n=0}^N{q}_n{x}^n}{{\Pi }_{n=0}^N(x+n)}(x+5.5{)}^{x+0.5}{e}^{-(x+5.5)}$$

$$\log \Gamma (x)=\log \left(\sum _{n=0}^N{q}_n{x}^n\right)+(x+0.5)\log (x+5.5)-(x+5.5)-\sum _{n=0}^N\log (x+n)$$

$$\begin{array}{r}q0=75122.6331530\\ q1=80916.6278952\\ q2=36308.2951477\\ q3=8687.24529705\\ q4=1168.92649479\\ q5=83.8676043424\\ q6=2.50662827511\end{array}$$

第2个:
$$\Gamma (x)=\sqrt{\frac{2\pi }{x}}{\left(\frac{x}{e}\sqrt{x\sinh (1/x)+\frac{1}{810x^6}}\right)}^x$$

因此:
$$\log \Gamma (x)=0.5\log (2\pi )+(x-0.5)\log (x)-x+0.5x\log \left(x\sinh (1/x)+\frac{1}{810x^6}\right)$$

优缺点

基于LSD的直线提取算法本质上是一种局部提取直线的算法（Hough算法为全局算法），因此其具有局部算法的缺点：
1）对于直线相交情况，因为设置了每个点是否USED，因此每个点只能属于一条直线，若有相交必有至少一条直线被割裂为两条。又因为其基于梯度变化，直线交点梯度值往往因较小而不被检测为边缘点，因此很有可能相交的两条直线在交点处被割裂为四条线段；
2）由于局部检测算法自增长的特点，对于长线段被遮挡、局部模糊等原因经常割裂为多条直线。因此，若实现具有普适性的直线段检测，需对其进行聚类和连接两方面的算法优化。

代码及注释

附上源码及注释,源码来源于OpenCV
头文件lsd.h:

//
// Created by liheng on 12/26/21.
//


#include <opencv2/imgproc.hpp>

//! @addtogroup imgproc_feature
//! @{

/** @brief Line segment detector class

following the algorithm described at @cite Rafael12 .

@note Implementation has been removed due original code license conflict

*/
class lLineSegmentDetector
{
public:

/** @brief Finds lines in the input image.

This is the output of the default parameters of the algorithm on the above shown image.

![image](pics/building_lsd.png)

@param _image A grayscale (CV_8UC1) input image. If only a roi needs to be selected, use:
`lsd_ptr-\>detect(image(roi), lines, ...); lines += Scalar(roi.x, roi.y, roi.x, roi.y);`
@param _lines A vector of Vec4i or Vec4f elements specifying the beginning and ending point of a line. Where
Vec4i/Vec4f is (x1, y1, x2, y2), point 1 is the start, point 2 - end. Returned lines are strictly
oriented depending on the gradient.
@param width Vector of widths of the regions, where the lines are found. E.g. Width of line.
@param prec Vector of precisions with which the lines are found.
@param nfa Vector containing number of false alarms in the line region, with precision of 10%. The
bigger the value, logarithmically better the detection.
- -1 corresponds to 10 mean false alarms
- 0 corresponds to 1 mean false alarm
- 1 corresponds to 0.1 mean false alarms
This vector will be calculated only when the objects type is #LSD_REFINE_ADV.
*/
virtual void detect(cv::InputArray _image, cv::OutputArray _lines,
                            cv::OutputArray width = cv::noArray(), cv::OutputArray prec = cv::noArray(),
                            cv::OutputArray nfa = cv::noArray()) = 0;

/** @brief Draws the line segments on a given image.
@param _image The image, where the lines will be drawn. Should be bigger or equal to the image,
where the lines were found.
@param lines A vector of the lines that needed to be drawn.
 */
virtual void drawSegments(cv::InputOutputArray _image, cv::InputArray lines) = 0;

/** @brief Draws two groups of lines in blue and red, counting the non overlapping (mismatching) pixels.

@param size The size of the image, where lines1 and lines2 were found.
@param lines1 The first group of lines that needs to be drawn. It is visualized in blue color.
@param lines2 The second group of lines. They visualized in red color.
@param _image Optional image, where the lines will be drawn. The image should be color(3-channel)
in order for lines1 and lines2 to be drawn in the above mentioned colors.
 */
virtual int compareSegments(const cv::Size& size, cv::InputArray lines1, cv::InputArray lines2, cv::InputOutputArray _image = cv::noArray()) = 0;

virtual ~lLineSegmentDetector() { }
};

/** @brief Creates a smart pointer to a LineSegmentDetector object and initializes it.

The LineSegmentDetector algorithm is defined using the standard values. Only advanced users may want
to edit those, as to tailor it for their own application.

@param _refine The way found lines will be refined, see #LineSegmentDetectorModes
            LSD_REFINE_NONE:没有改良的方式
            LSD_REFINE_STD:标准改良方式,将带弧度的线拆成多个可以逼近原线段的直线段
            LSD_REFINE_ADV:进一步改良方式,计算出错误警告数量,通过增加精度,减少尺寸进一步精确直线
@param _scale The scale of the image that will be used to find the lines. Range (0..1].
@param _sigma_scale Sigma for Gaussian filter. It is computed as sigma = _sigma_scale/_scale.
@param _quant Bound to the quantization error on the gradient norm.
@param _ang_th Gradient angle tolerance in degrees.
@param _log_eps Detection threshold: -log10(NFA) \> log_eps. Used only when advance refinement
is chosen.
@param _density_th Minimal density of aligned region points in the enclosing rectangle.
@param _n_bins Number of bins in pseudo-ordering of gradient modulus.

@note Implementation has been removed due original code license conflict
 */
std::shared_ptr<lLineSegmentDetector> createlLineSegmentDetector(
        int _refine = cv::LSD_REFINE_STD, double _scale = 0.8,
        double _sigma_scale = 0.6, double _quant = 2.0, double _ang_th = 22.5,
        double _log_eps = 0, double _density_th = 0.7, int _n_bins = 1024);

//! @} imgproc_feature

实现:

/*M///
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                           License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistributions of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistributions in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//   * The name of the copyright holders may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/

#include "lsd.h"
#include <vector>
#include <opencv2/opencv.hpp>

/
// Default LSD parameters
// SIGMA_SCALE 0.6    - Sigma for Gaussian filter is computed as sigma = sigma_scale/scale.
// QUANT       2.0    - Bound to the quantization error on the gradient norm.
// ANG_TH      22.5   - Gradient angle tolerance in degrees.
// LOG_EPS     0.0    - Detection threshold: -log10(NFA) > log_eps
// DENSITY_TH  0.7    - Minimal density of region points in rectangle.
// N_BINS      1024   - Number of bins in pseudo-ordering of gradient modulus.

#define M_3_2_PI    (3 * CV_PI) / 2   // 3/2 pi
#define M_2__PI     (2 * CV_PI)         // 2 pi

#ifndef M_LN10
#define M_LN10      2.30258509299404568402
#endif

#define NOTDEF      double(-1024.0) // Label for pixels with undefined gradient.

#define NOTUSED     0   // Label for pixels not used in yet.
#define USED        1   // Label for pixels already used in detection.

#define RELATIVE_ERROR_FACTOR 100.0

const double DEG_TO_RADS = CV_PI / 180;

#define log_gamma(x) ((x)>15.0?log_gamma_windschitl(x):log_gamma_lanczos(x))

struct edge
{
    cv::Point p;
    bool taken;
};

/

inline double distSq(const double x1, const double y1,
                     const double x2, const double y2)
{
    return (x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1);
}

inline double dist(const double x1, const double y1,
                   const double x2, const double y2)
{
    return sqrt(distSq(x1, y1, x2, y2));
}

// Signed angle difference
inline double angle_diff_signed(const double& a, const double& b)
{
    double diff = a - b;
    while(diff <= -CV_PI) diff += M_2__PI;
    while(diff >   CV_PI) diff -= M_2__PI;
    return diff;
}

// Absolute value angle difference
inline double angle_diff(const double& a, const double& b)
{
    return std::fabs(angle_diff_signed(a, b));
}

// Compare doubles by relative error.
inline bool double_equal(const double& a, const double& b)
{
    // trivial case
    if(a == b) return true;

    double abs_diff = fabs(a - b);
    double aa = fabs(a);
    double bb = fabs(b);
    double abs_max = (aa > bb)? aa : bb;

    if(abs_max < DBL_MIN) abs_max = DBL_MIN;

    return (abs_diff / abs_max) <= (RELATIVE_ERROR_FACTOR * DBL_EPSILON);
}

inline bool AsmallerB_XoverY(const edge& a, const edge& b)
{
    if (a.p.x == b.p.x) return a.p.y < b.p.y;
    else return a.p.x < b.p.x;
}

/**
 *   Computes the natural logarithm of the absolute value of
 *   the gamma function of x using Windschitl method.
 *   See http://www.rskey.org/gamma.htm
 */
inline double log_gamma_windschitl(const double& x)
{
    return 0.918938533204673 + (x-0.5)*log(x) - x
           + 0.5*x*log(x*sinh(1/x) + 1/(810.0*pow(x, 6.0)));
}

/**
 *   Computes the natural logarithm of the absolute value of
 *   the gamma function of x using the Lanczos approximation.
 *   See http://www.rskey.org/gamma.htm
 */
inline double log_gamma_lanczos(const double& x)
{
    static double q[7] = { 75122.6331530, 80916.6278952, 36308.2951477,
                           8687.24529705, 1168.92649479, 83.8676043424,
                           2.50662827511 };
    double a = (x + 0.5) * log(x + 5.5) - (x + 5.5);
    double b = 0;
    for(int n = 0; n < 7; ++n)
    {
        a -= log(x + double(n));
        b += q[n] * pow(x, double(n));
    }
    return a + log(b);
}
///

class LineSegmentDetectorImpl final : public lLineSegmentDetector
{
    public:

/**
 * Create a LineSegmentDetectorImpl object. Specifying scale, number of subdivisions for the image, should the lines be refined and other constants as follows:
 *
 * @param _refine       How should the lines found be refined?
 *                      LSD_REFINE_NONE - No refinement applied.
 *                      LSD_REFINE_STD  - Standard refinement is applied. E.g. breaking arches into smaller line approximations.
 *                      LSD_REFINE_ADV  - Advanced refinement. Number of false alarms is calculated,
 *                                    lines are refined through increase of precision, decrement in size, etc.
 * @param _scale        The scale of the image that will be used to find the lines. Range (0..1].
 * @param _sigma_scale  Sigma for Gaussian filter is computed as sigma = _sigma_scale/_scale.
 * @param _quant        Bound to the quantization error on the gradient norm.
 * @param _ang_th       Gradient angle tolerance in degrees.
 * @param _log_eps      Detection threshold: -log10(NFA) > _log_eps
 * @param _density_th   Minimal density of aligned region points in rectangle.
 * @param _n_bins       Number of bins in pseudo-ordering of gradient modulus.
 */
    LineSegmentDetectorImpl(int _refine = cv::LSD_REFINE_STD, double _scale = 0.8,
                            double _sigma_scale = 0.6, double _quant = 2.0, double _ang_th = 22.5,
                            double _log_eps = 0, double _density_th = 0.7, int _n_bins = 1024);

/**
 * Detect lines in the input image.
 *
 * @param _image    A grayscale(CV_8UC1) input image.
 *                  If only a roi needs to be selected, use
 *                  lsd_ptr->detect(image(roi), ..., lines);
 *                  lines += cv::Scalar(roi.x, roi.y, roi.x, roi.y);
 * @param _lines    Return: A vector of Vec4f elements specifying the beginning and ending point of a line.
 *                          Where Vec4f is (x1, y1, x2, y2), point 1 is the start, point 2 - end.
 *                          Returned lines are strictly oriented depending on the gradient.
 * @param width     Return: Vector of widths of the regions, where the lines are found. E.g. Width of line.
 * @param prec      Return: Vector of precisions with which the lines are found.
 * @param nfa       Return: Vector containing number of false alarms in the line region, with precision of 10%.
 *                          The bigger the value, logarithmically better the detection.
 *                              * -1 corresponds to 10 mean false alarms
 *                              * 0 corresponds to 1 mean false alarm
 *                              * 1 corresponds to 0.1 mean false alarms
 *                          This vector will be calculated _only_ when the objects type is REFINE_ADV
 */
    void detect(cv::InputArray _image, cv::OutputArray _lines,
                cv::OutputArray width = cv::noArray(), cv::OutputArray prec = cv::noArray(),
                cv::OutputArray nfa = cv::noArray()) override;

/**
 * Draw lines on the given canvas.
 *
 * @param image     The image, where lines will be drawn.
 *                  Should have the size of the image, where the lines were found
 * @param lines     The lines that need to be drawn
 */
    void drawSegments(cv::InputOutputArray _image, cv::InputArray lines) override;

/**
 * Draw both vectors on the image canvas. Uses blue for lines 1 and red for lines 2.
 *
 * @param size      The size of the image, where lines1 and lines2 were found.
 * @param lines1    The first lines that need to be drawn. Color - Blue.
 * @param lines2    The second lines that need to be drawn. Color - Red.
 * @param image     An optional image, where lines will be drawn.
 *                  Should have the size of the image, where the lines were found
 * @return          The number of mismatching pixels between lines1 and lines2.
 */
    int compareSegments(const cv::Size& size, cv::InputArray lines1, cv::InputArray lines2, cv::InputOutputArray _image = cv::noArray()) override;

    private:
    cv::Mat image;//输入图像CV_8UC1
    cv::Mat scaled_image;//原图:scale=1.or (高斯滤波+缩放)后的图像:scale!=1
    cv::Mat_<double> angles;// 梯度幅角  in rads
    cv::Mat_<double> modgrad;//梯度幅值
    cv::Mat_<uchar> used;//状态列表.该点是否已属于某线段

    int img_width;//缩放后的图像尺寸
    int img_height;
    double LOG_NT;

    bool w_needed;
    bool p_needed;
    bool n_needed;

    const double SCALE;
    const int doRefine;
    const double SIGMA_SCALE;
    const double QUANT;
    const double ANG_TH;
    const double LOG_EPS;
    const double DENSITY_TH;//对齐点 密度阈值
    const int N_BINS;

    struct RegionPoint {
        int x;
        int y;
        uchar* used;//此处是指针类型!!,指向状态列表cv::Mat_<uchar> used中(y,x)处的元素
        double angle;
        double modgrad;
    };

    struct normPoint
    {
        cv::Point2i p;//在图像中的位置
        int norm;//量化后的梯度幅值
    };

    std::vector<normPoint> ordered_points;//将图像中各点梯度伪排序后的点

    struct rect
    {
        double x1, y1, x2, y2;    // first and second point of the line segment
        double width;             // rectangle width
        double x, y;              // center of the rectangle
        double theta;             // angle
        double dx,dy;             // (dx,dy) is vector oriented as the line segment
        double prec;              // tolerance angle
        double p;                 // probability of a point with angle within 'prec'
    };

    LineSegmentDetectorImpl& operator= (const LineSegmentDetectorImpl&); // to quiet MSVC

/**
 * Detect lines in the whole input image.
 *
 * @param lines         Return: A vector of Vec4f elements specifying the beginning and ending point of a line.
 *                              Where Vec4f is (x1, y1, x2, y2), point 1 is the start, point 2 - end.
 *                              Returned lines are strictly oriented depending on the gradient.
 * @param widths        Return: Vector of widths of the regions, where the lines are found. E.g. Width of line.
 * @param precisions    Return: Vector of precisions with which the lines are found.
 * @param nfas          Return: Vector containing number of false alarms in the line region, with precision of 10%.
 *                              The bigger the value, logarithmically better the detection.
 *                                  * -1 corresponds to 10 mean false alarms
 *                                  * 0 corresponds to 1 mean false alarm
 *                                  * 1 corresponds to 0.1 mean false alarms
 */
    void flsd(std::vector<cv::Vec4f>& lines,
              std::vector<double>& widths, std::vector<double>& precisions,
              std::vector<double>& nfas);

/**
 * Finds the angles and the gradients of the image. Generates a list of pseudo ordered points.
 * 计算图像各点的梯度幅值与幅角,并按梯度幅值进行伪排序
 *
 * @param threshold      The minimum value of the angle that is considered defined, otherwise NOTDEF
 * @param n_bins         The number of bins with which gradients are ordered by, using bucket sort.
 * @param ordered_points Return: Vector of coordinate points that are pseudo ordered by magnitude.
 *                       Pixels would be ordered by norm value, up to a precision given by max_grad/n_bins.
 */
    void ll_angle(const double& threshold, const unsigned int& n_bins);

/**
 * Grow a region starting from point s with a defined precision,
 * returning the containing points size and the angle of the gradients.
 *
 * @param s         Starting point for the region.
 * @param reg       Return: Vector of points, that are part of the region
 * @param reg_angle Return: The mean angle of the region.
 * @param prec      The precision by which each region angle should be aligned to the mean.
 */
    void region_grow(const cv::Point2i& s, std::vector<RegionPoint>& reg,
                     double& reg_angle, const double& prec);

/**
 * Finds the bounding rotated rectangle of a region.
 *
 * @param reg       The region of points, from which the rectangle to be constructed from.
 * @param reg_angle The mean angle of the region.
 * @param prec      The precision by which points were found.
 * @param p         Probability of a point with angle within 'prec'.
 * @param rec       Return: The generated rectangle.
 */
    void region2rect(const std::vector<RegionPoint>& reg, const double reg_angle,
                     const double prec, const double p, rect& rec) const;

/**
 * Compute region's angle as the principal inertia axis of the region.
 * @return          Regions angle.
 */
    double get_theta(const std::vector<RegionPoint>& reg, const double& x,
                     const double& y, const double& reg_angle, const double& prec) const;

/**
 * An estimation of the angle tolerance is performed by the standard deviation of the angle at points
 * near the region's starting point. Then, a new region is grown starting from the same point, but using the
 * estimated angle tolerance. If this fails to produce a rectangle with the right density of region points,
 * 'reduce_region_radius' is called to try to satisfy this condition.
 */
    bool refine(std::vector<RegionPoint>& reg, double reg_angle,
                const double prec, double p, rect& rec, const double& density_th);

/**
 * Reduce the region size, by elimination the points far from the starting point, until that leads to
 * rectangle with the right density of region points or to discard the region if too small.
 */
    bool reduce_region_radius(std::vector<RegionPoint>& reg, double reg_angle,
                              const double prec, double p, rect& rec, double density, const double& density_th);

/**
 * Try some rectangles variations to improve NFA value. Only if the rectangle is not meaningful (i.e., log_nfa <= log_eps).
 * @return      The new NFA value.
 */
    double rect_improve(rect& rec) const;

/**
 * Calculates the number of correctly aligned points within the rectangle.
 * @return      The new NFA value.
 */
    double rect_nfa(const rect& rec) const;

/**
 * Computes the NFA values based on the total number of points, points that agree.
 * n, k, p are the binomial parameters.
 * @return      The new NFA value.
 */
    double nfa(const int& n, const int& k, const double& p) const;

    /**
     * Is the point at place 'address' aligned to angle theta, up to precision 'prec'?
     * 区域生长中,根据梯度方向,判断点(x,y)能否和种子点进行合并
     * (x,y),theta 待判断点的坐标及梯度方向
     * prec 梯度幅角合并阈值
     * @return      Whether the point is aligned.
     */
    bool isAligned(int x, int y, const double& theta, const double& prec) const;

    public:
    // Compare norm
    static inline bool compare_norm( const normPoint& n1, const normPoint& n2 )
    {
        return (n1.norm > n2.norm);
    }
};

/

std::shared_ptr<lLineSegmentDetector> createlLineSegmentDetector(
        int _refine, double _scale, double _sigma_scale, double _quant, double _ang_th,
        double _log_eps, double _density_th, int _n_bins)
{
    return std::make_shared<LineSegmentDetectorImpl>(
            _refine, _scale, _sigma_scale, _quant, _ang_th,
            _log_eps, _density_th, _n_bins);
}

/

LineSegmentDetectorImpl::LineSegmentDetectorImpl(int _refine, double _scale, double _sigma_scale, double _quant,
                                                 double _ang_th, double _log_eps, double _density_th, int _n_bins)
        : img_width(0), img_height(0), LOG_NT(0), w_needed(false), p_needed(false), n_needed(false),
          SCALE(_scale), doRefine(_refine), SIGMA_SCALE(_sigma_scale), QUANT(_quant),
          ANG_TH(_ang_th), LOG_EPS(_log_eps), DENSITY_TH(_density_th), N_BINS(_n_bins)
{
    CV_Assert(_scale > 0 && _sigma_scale > 0 && _quant >= 0 &&
              _ang_th > 0 && _ang_th < 180 && _density_th >= 0 && _density_th < 1 &&
              _n_bins > 0);
}

void LineSegmentDetectorImpl::detect(cv::InputArray _image, cv::OutputArray _lines,
                                     cv::OutputArray _width, cv::OutputArray _prec, cv::OutputArray _nfa)
{
    

    image = _image.getMat();
    CV_Assert(!image.empty() && image.type() == CV_8UC1);

    std::vector<cv::Vec4f> lines;
    std::vector<double> w, p, n;
    w_needed = _width.needed();
    p_needed = _prec.needed();
    if (doRefine < cv::LSD_REFINE_ADV)
        n_needed = false;
    else
        n_needed = _nfa.needed();

    flsd(lines, w, p, n);

    cv::Mat(lines).copyTo(_lines);
    if(w_needed) cv::Mat(w).copyTo(_width);
    if(p_needed) cv::Mat(p).copyTo(_prec);
    if(n_needed) cv::Mat(n).copyTo(_nfa);

    // Clear used structures
    ordered_points.clear();
}

void LineSegmentDetectorImpl::flsd(std::vector<cv::Vec4f>& lines,
                                   std::vector<double>& widths, std::vector<double>& precisions,
                                   std::vector<double>& nfas)
{
    // Angle tolerance
    const double prec = CV_PI * ANG_TH / 180;//角度阈值,待选点能够和已确定点进行合并的梯度幅角阈值
    const double p = ANG_TH / 180;//角度符合level-line角度阈值的概率
    const double rho = QUANT / sin(prec);    // gradient magnitude threshold.小梯度抑制中的梯度阈值

    if(SCALE != 1)
    {
        cv::Mat gaussian_img;
        const double sigma = (SCALE < 1)?(SIGMA_SCALE / SCALE):(SIGMA_SCALE);
        const double sprec = 3;
        const unsigned int h =  (unsigned int)(ceil(sigma * sqrt(2 * sprec * log(10.0))));
        cv::Size ksize(1 + 2 * h, 1 + 2 * h); // kernel size
        cv::GaussianBlur(image, gaussian_img, ksize, sigma);
        // Scale image to needed size
        cv::resize(gaussian_img, scaled_image, cv::Size(), SCALE, SCALE, cv::INTER_LINEAR_EXACT);//位精确双线性插值
        ll_angle(rho, N_BINS);//计算梯度幅值,幅角,并进行伪排序
    }
    else
    {
        scaled_image = image;
        ll_angle(rho, N_BINS);
    }

    //  11 * (X*Y)^(5/2)  whose logarithm value is
    //   log10(11) + 5/2 * (log10(X) + log10(Y)).
    LOG_NT = 5 * (log10(double(img_width)) + log10(double(img_height))) / 2 + log10(11.0);
    const size_t min_reg_size = size_t(-LOG_NT/log10(p)); // minimal number of points in region that can give a meaningful event

    // // Initialize region only when needed
    // cv::Mat region = cv::Mat::zeros(scaled_image.size(), CV_8UC1);
    used = cv::Mat_<uchar>::zeros(scaled_image.size()); // zeros = NOTUSED
    std::vector<RegionPoint> reg;

    // Search for line segments
    for(size_t i = 0, points_size = ordered_points.size(); i < points_size; ++i)
    {
        //cv::Mat modgradU;
        //cv::normalize(modgrad,modgradU,0,255,cv::NORM_MINMAX,CV_8UC1);
        //cv::cvtColor(modgradU,modgradU,cv::COLOR_GRAY2BGR);
        //cv::imshow("modgradU",modgradU);

        const cv::Point2i& point = ordered_points[i].p;
        if((used.at<uchar>(point) == NOTUSED) && (angles.at<double>(point) != NOTDEF))
        {
            //找出位置相邻接且方向相近的点构成的区域
            double reg_angle;
            region_grow(ordered_points[i].p, reg, reg_angle, prec);//对该点进行区域生长,得到和该点梯度方向一致的区域


            //区域点过少,舍去该区域
            // Ignore small regions
            if(reg.size() < min_reg_size)
            { continue; }

            //for(int k=0;k<reg.size();++k)
            //    modgradU.at<cv::Vec3b>(reg[k].y,reg[k].x) = cv::Vec3b(0,255,0);
            //modgradU.at<cv::Vec3b>(ordered_points[i].p.y,ordered_points[i].p.x) = cv::Vec3b(0,0,255);
            //cv::imshow("modgradURegion",modgradU);


            //检查矩形中点的密度,若小于阈值,则改进矩形表示.若仍旧不满足阈值条件,则清除该区域
            // Construct rectangular approximation for the region
            rect rec;
            region2rect(reg, reg_angle, prec, p, rec);//reg_angle:该区域内所有点的梯度幅角;prec:点属于同一区域的角度偏差容许值

            double log_nfa = -1;
            if(doRefine > cv::LSD_REFINE_NONE)
            {
                // At least REFINE_STANDARD lvl.
                //调整角度容差阈值和缩小区域半径的方法来对区域进行精炼.若通过两种方法均不能满足要求,
                //则继续下一轮遍历
                if(!refine(reg, reg_angle, prec, p, rec, DENSITY_TH))
                { continue; }

                //精炼满足要求后,对矩形进行改进.
                if(doRefine >= cv::LSD_REFINE_ADV)
                {
                    // Compute NFA
                    log_nfa = rect_improve(rec);
                    if(log_nfa <= LOG_EPS)
                    { continue; }
                }
            }
            // Found new line
            //找到一条新的线段


            //补偿梯度计算时的误差
            // Add the offset
            rec.x1 += 0.5; rec.y1 += 0.5;
            rec.x2 += 0.5; rec.y2 += 0.5;

            //尺寸恢复到原图
            // scale the result values if a sub-sampling was performed
            if(SCALE != 1)
            {
                rec.x1 /= SCALE; rec.y1 /= SCALE;
                rec.x2 /= SCALE; rec.y2 /= SCALE;
                rec.width /= SCALE;
            }

            //Store the relevant data
            lines.push_back(cv::Vec4f(float(rec.x1), float(rec.y1), float(rec.x2), float(rec.y2)));

            if(w_needed) widths.push_back(rec.width);
            if(p_needed) precisions.push_back(rec.p);
            if(n_needed && doRefine >= cv::LSD_REFINE_ADV) nfas.push_back(log_nfa);
        }
    }
}

void LineSegmentDetectorImpl::ll_angle(const double& threshold,
                                       const unsigned int& n_bins)
{
    //Initialize data
    angles = cv::Mat_<double>(scaled_image.size());
    modgrad = cv::Mat_<double>(scaled_image.size());

    img_width = scaled_image.cols;
    img_height = scaled_image.rows;

    // Undefined the down and right boundaries
    angles.row(img_height - 1).setTo(NOTDEF);
    angles.col(img_width - 1).setTo(NOTDEF);

    //求解各点梯度幅值和幅角.计算2x2的区域,因此实际计算的梯度是(x+0.5,y+0.5)处点的梯度,后续需要对该计算偏差进行补偿
    // Computing gradient for remaining pixels
    double max_grad = -1;//梯度最大值
    for(int y = 0; y < img_height - 1; ++y)
    {
        const uchar* scaled_image_row = scaled_image.ptr<uchar>(y);
        const uchar* next_scaled_image_row = scaled_image.ptr<uchar>(y+1);
        double* angles_row = angles.ptr<double>(y);
        double* modgrad_row = modgrad.ptr<double>(y);
        for(int x = 0; x < img_width-1; ++x)
        {
            int DA = next_scaled_image_row[x + 1] - scaled_image_row[x];
            int BC = scaled_image_row[x + 1] - next_scaled_image_row[x];
            int gx = DA + BC;    // gradient x component
            int gy = DA - BC;    // gradient y component
            double norm = std::sqrt((gx * gx + gy * gy) / 4.0); // gradient norm

            modgrad_row[x] = norm;    // store gradient

            //小梯度抑制
            //梯度过小,在幅角图中忽略该点
            if (norm <= threshold)  // norm too small, gradient no defined
            {
                angles_row[x] = NOTDEF;
            }
            else
            {
                angles_row[x] = cv::fastAtan2(float(gx), float(-gy)) * DEG_TO_RADS;  // gradient angle computation
                if (norm > max_grad) { max_grad = norm; }
            }

        }
    }

    //计算梯度直方图.
    //或者说,将梯度量化为[0,1024).量化后梯度最小为0,最大为1023.
    // Compute histogram of gradient values
    double bin_coef = (max_grad > 0) ? double(n_bins - 1) / max_grad : 0; // If all image is smooth, max_grad <= 0
    for(int y = 0; y < img_height - 1; ++y)//将图像中所有点的梯度进行量化,保存
    {
        const double* modgrad_row = modgrad.ptr<double>(y);
        for(int x = 0; x < img_width - 1; ++x)
        {
            normPoint _point;
            int i = int(modgrad_row[x] * bin_coef);
            _point.p = cv::Point(x, y);//该点在图像中的位置
            _point.norm = i;//量化后的梯度幅值
            ordered_points.push_back(_point);
        }
    }

    // Sort.量化后的梯度幅值从大到小进行排序
    std::sort(ordered_points.begin(), ordered_points.end(), compare_norm);
}

void LineSegmentDetectorImpl::region_grow(const cv::Point2i& s, std::vector<RegionPoint>& reg,
                                          double& reg_angle, const double& prec)
{
    reg.clear();

    // Point to this region
    RegionPoint seed;
    seed.x = s.x;
    seed.y = s.y;
    seed.used = &used.at<uchar>(s);
    reg_angle = angles.at<double>(s);//初值为种子点的梯度幅角
    seed.angle = reg_angle;
    seed.modgrad = modgrad.at<double>(s);
    reg.push_back(seed);//区域的第一个点来源于梯度最大值的点,作为种子点

    //sumdx,sumdy 区域各点梯度的x,y方向的分量
    float sumdx = float(std::cos(reg_angle));
    float sumdy = float(std::sin(reg_angle));
    *seed.used = USED;

    //对区域内的各个种子点邻域遍历,进行区域生长
    //Try neighboring regions
    for (size_t i = 0;i<reg.size();i++)
    {
        //8邻域点
        const RegionPoint& rpoint = reg[i];
        int xx_min = std::max(rpoint.x - 1, 0), xx_max = std::min(rpoint.x + 1, img_width - 1);
        int yy_min = std::max(rpoint.y - 1, 0), yy_max = std::min(rpoint.y + 1, img_height - 1);
        for(int yy = yy_min; yy <= yy_max; ++yy)
        {
            uchar* used_row = used.ptr<uchar>(yy);
            const double* angles_row = angles.ptr<double>(yy);
            const double* modgrad_row = modgrad.ptr<double>(yy);
            for(int xx = xx_min; xx <= xx_max; ++xx)
            {
                uchar& is_used = used_row[xx];
                if(is_used != USED &&
                   (isAligned(xx, yy, reg_angle, prec)))//该点未被标记且梯度幅角和和种子点的梯度幅角/方向接近
                {
                    const double& angle = angles_row[xx];
                    // Add point
                    is_used = USED;
                    RegionPoint region_point;
                    region_point.x = xx;
                    region_point.y = yy;
                    region_point.used = &is_used;
                    region_point.modgrad = modgrad_row[xx];
                    region_point.angle = angle;
                    reg.push_back(region_point);//将该点加入区域中,下次进行遍历

                    //根据区域内现有点的总梯度,计算该区域的梯度方向,并更新,以使下一个点进行判断
                    //根据矢量相加:总梯度=梯度1+梯度2+...
                    // Update region's angle
                    sumdx += cos(float(angle));
                    sumdy += sin(float(angle));
                    // reg_angle is used in the isAligned, so it needs to be updates?
                    reg_angle = cv::fastAtan2(sumdy, sumdx) * DEG_TO_RADS;
                }
            }
        }
    }
}

void LineSegmentDetectorImpl::region2rect(const std::vector<RegionPoint>& reg,
                                          const double reg_angle, const double prec, const double p, rect& rec) const
{
    //求解该矩形区域的重心(所有点的重心.每个点的质量为梯度幅值)
    double x = 0, y = 0, sum = 0;
    for(size_t i = 0; i < reg.size(); ++i)
    {
        const RegionPoint& pnt = reg[i];
        const double& weight = pnt.modgrad;
        x += double(pnt.x) * weight;//利用所有点的梯度规范化值平均计算重心.可以将梯度幅值理解为点的重量
        y += double(pnt.y) * weight;
        sum += weight;
    }

    // Weighted sum must differ from 0
    CV_Assert(sum > 0);

    x /= sum;
    y /= sum;

    double theta = get_theta(reg, x, y, reg_angle, prec);

    //求解最小外接矩形参数.'类似'于OpenCV的旋转矩形.
    //求解公式来历?
    // Find length and width
    double dx = cos(theta);
    double dy = sin(theta);
    double l_min = 0, l_max = 0, w_min = 0, w_max = 0;

    //l和w表示各点 沿着主方向 到矩形中心的距离.l_max-lmin表示旋转矩形的长度.
    for(size_t i = 0; i < reg.size(); ++i)
    {
        double regdx = double(reg[i].x) - x;
        double regdy = double(reg[i].y) - y;

        double l = regdx * dx + regdy * dy;
        double w = -regdx * dy + regdy * dx;

        if(l > l_max) l_max = l;
        else if(l < l_min) l_min = l;
        if(w > w_max) w_max = w;
        else if(w < w_min) w_min = w;
    }

    // Store values
    rec.x1 = x + l_min * dx;
    rec.y1 = y + l_min * dy;
    rec.x2 = x + l_max * dx;
    rec.y2 = y + l_max * dy;
    rec.width = w_max - w_min;
    rec.x = x;
    rec.y = y;
    rec.theta = theta;
    rec.dx = dx;
    rec.dy = dy;
    rec.prec = prec;
    rec.p = p;

    // Min width of 1 pixel
    if(rec.width < 1.0) rec.width = 1.0;
}

double LineSegmentDetectorImpl::get_theta(const std::vector<RegionPoint>& reg, const double& x,
                                          const double& y, const double& reg_angle, const double& prec) const
{
    double Ixx = 0.0;
    double Iyy = 0.0;
    double Ixy = 0.0;

    //惯量矩阵 or  协方差矩阵
    // Compute inertia matrix
    for(size_t i = 0; i < reg.size(); ++i)
    {
        const double& regx = reg[i].x;
        const double& regy = reg[i].y;
        const double& weight = reg[i].modgrad;
        double dx = regx - x;
        double dy = regy - y;
        Ixx += dy * dy * weight;//此处应该是Iyy?
        Iyy += dx * dx * weight;//此处应该是Ixx?.笔误?.对最终结果无影响
        Ixy -= dx * dy * weight;
    }

    //计算特征值.
    // Check if inertia matrix is null
    CV_Assert(!(double_equal(Ixx, 0) && double_equal(Iyy, 0) && double_equal(Ixy, 0)));

    // Compute smallest eigenvalue
    double lambda = 0.5 * (Ixx + Iyy - sqrt((Ixx - Iyy) * (Ixx - Iyy) + 4.0 * Ixy * Ixy));

    // Compute angle
    double theta = (fabs(Ixx)>fabs(Iyy))?
                   double(cv::fastAtan2(float(lambda - Ixx), float(Ixy))):
                   double(cv::fastAtan2(float(Ixy), float(lambda - Iyy))); // in degs
    theta *= DEG_TO_RADS;

    // Correct angle by 180 deg if necessary
    if(angle_diff(theta, reg_angle) > prec) { theta += CV_PI; }



//    //上面部分:也可以通过求解协方差矩阵,求解特征向量来代替,
//    //得到特征向量后,求特征向量与x轴的夹角.这样理解会比较简单.
//    //求得是最小特征值对应的特征向量
//    cv::Mat D(2,reg.size(),CV_32FC1);
//    for(size_t i = 0; i < reg.size(); ++i)
//    {
//        const double& weight = std::sqrt(reg[i].modgrad);
//        D.at<float>(0,i) = (reg[i].x-x)*weight;
//        D.at<float>(1,i) = (reg[i].y-y)*weight;
//    }
//
//    cv::Mat C = D*D.t();
//    cv::Mat eigenvalue,eigenVector;
//    cv::eigen(C,eigenvalue,eigenVector);
//
//    auto a = cv::fastAtan2(eigenVector.at<float>(0,0),eigenVector.at<float>(1,0));
//    auto my_theta = cv::fastAtan2(eigenVector.at<float>(0,1),eigenVector.at<float>(1,1));

    return theta;
}

bool LineSegmentDetectorImpl::refine(std::vector<RegionPoint>& reg, double reg_angle,
                                     const double prec, double p, rect& rec, const double& density_th)
{
    //计算区域内点的密度 点的数量/矩形的面积
    double density = double(reg.size()) / (dist(rec.x1, rec.y1, rec.x2, rec.y2) * rec.width);

    if (density >= density_th) { return true; }

    // Try to reduce angle tolerance
    double xc = double(reg[0].x);
    double yc = double(reg[0].y);//将第一个点作为种子点
    const double& ang_c = reg[0].angle;
    double sum = 0, s_sum = 0;
    int n = 0;

    for (size_t i = 0; i < reg.size(); ++i)
    {
        *(reg[i].used) = NOTUSED;//释放该点,后续被再次使用.该矩形区域内所有点状态均重置为NOTUSED.指针操作!!
        if (dist(xc, yc, reg[i].x, reg[i].y) < rec.width)//距离种子点的距离小于矩形的宽度
        {
            const double& angle = reg[i].angle;
            double ang_d = angle_diff_signed(angle, ang_c);//unit:rad
            sum += ang_d;
            s_sum += ang_d * ang_d;
            ++n;
        }
        //else 不满足if条件的点,其之后的点均不满足上面的判断条件,
        //可将该 for循环及if判断修改,减少计算
    }
    CV_Assert(n > 0);
    double mean_angle = sum / double(n);
    // 2 * standard deviation
    //求解 标准差*2
    double tau = 2.0 * sqrt((s_sum - 2.0 * mean_angle * sum) / double(n) + mean_angle * mean_angle);

    // Try new region
    //以该矩形区域的第一个点重新生长,角度容差利用上面的2*标准差.该角度容差比原设置值会稍大
    region_grow(cv::Point(reg[0].x, reg[0].y), reg, reg_angle, tau);

    if (reg.size() < 2)//此处判断条件,建议利用已经定义的变量 min_reg_size ,可能要合适点.
    { return false; }
    //上面选择2的原因是,region2rect 需要至少2个点才能求解一个矩形;2个点才能确定一个直线段

    region2rect(reg, reg_angle, prec, p, rec);
    density = double(reg.size()) / (dist(rec.x1, rec.y1, rec.x2, rec.y2) * rec.width);

    if (density < density_th)
    {
        //通过调整角度容差,仍不能满足要求,则 尝试通过 缩小区域半径 来满足要求.
        return reduce_region_radius(reg, reg_angle, prec, p, rec, density, density_th);
    }
    else
    {
        return true;
    }
}

bool LineSegmentDetectorImpl::reduce_region_radius(std::vector<RegionPoint>& reg, double reg_angle,
                                                   const double prec, double p, rect& rec, double density, const double& density_th)
{
    // Compute region's radius
    double xc = double(reg[0].x);
    double yc = double(reg[0].y);
    double radSq1 = distSq(xc, yc, rec.x1, rec.y1);
    double radSq2 = distSq(xc, yc, rec.x2, rec.y2);
    double radSq = radSq1 > radSq2 ? radSq1 : radSq2;

    while(density < density_th)
    {
        radSq *= 0.75*0.75; // Reduce region's radius to 75% of its value.由于radSq是平方形式,因此是 0.75*0.75
        // Remove points from the region and update 'used' map
        for (size_t i = 0; i < reg.size(); ++i)
        {
            if(distSq(xc, yc, double(reg[i].x), double(reg[i].y)) > radSq)
            {
                // Remove point from the region
                *(reg[i].used) = NOTUSED;
                std::swap(reg[i], reg[reg.size() - 1]);//移到最后,并删除
                reg.pop_back();
                --i; // To avoid skipping one point
            }
        }

        if(reg.size() < 2)//同上.此处判断条件,建议利用已经定义的变量 min_reg_size ,可能要合适点.
        { return false; }

        // Re-compute rectangle
        region2rect(reg ,reg_angle, prec, p, rec);

        // Re-compute region points density
        density = double(reg.size()) /
                  (dist(rec.x1, rec.y1, rec.x2, rec.y2) * rec.width);
    }

    return true;
}

double LineSegmentDetectorImpl::rect_improve(rect& rec) const
{
    double delta = 0.5;
    double delta_2 = delta / 2.0;

    double log_nfa = rect_nfa(rec);

    if(log_nfa > LOG_EPS) return log_nfa; // Good rectangle

    // Try to improve
    // Finer precision
    rect r = rect(rec); // Copy
    for(int n = 0; n < 5; ++n)
    {
        r.p /= 2;
        r.prec = r.p * CV_PI;
        double log_nfa_new = rect_nfa(r);
        if(log_nfa_new > log_nfa)
        {
            log_nfa = log_nfa_new;
            rec = rect(r);
        }
    }
    if(log_nfa > LOG_EPS) return log_nfa;

    // Try to reduce width
    r = rect(rec);
    for(unsigned int n = 0; n < 5; ++n)
    {
        if((r.width - delta) >= 0.5)
        {
            r.width -= delta;
            double log_nfa_new = rect_nfa(r);
            if(log_nfa_new > log_nfa)
            {
                rec = rect(r);
                log_nfa = log_nfa_new;
            }
        }
    }
    if(log_nfa > LOG_EPS) return log_nfa;

    // Try to reduce one side of rectangle
    r = rect(rec);
    for(unsigned int n = 0; n < 5; ++n)
    {
        if((r.width - delta) >= 0.5)
        {
            r.x1 += -r.dy * delta_2;
            r.y1 +=  r.dx * delta_2;
            r.x2 += -r.dy * delta_2;
            r.y2 +=  r.dx * delta_2;
            r.width -= delta;
            double log_nfa_new = rect_nfa(r);
            if(log_nfa_new > log_nfa)
            {
                rec = rect(r);
                log_nfa = log_nfa_new;
            }
        }
    }
    if(log_nfa > LOG_EPS) return log_nfa;

    // Try to reduce other side of rectangle
    r = rect(rec);
    for(unsigned int n = 0; n < 5; ++n)
    {
        if((r.width - delta) >= 0.5)
        {
            r.x1 -= -r.dy * delta_2;
            r.y1 -=  r.dx * delta_2;
            r.x2 -= -r.dy * delta_2;
            r.y2 -=  r.dx * delta_2;
            r.width -= delta;
            double log_nfa_new = rect_nfa(r);
            if(log_nfa_new > log_nfa)
            {
                rec = rect(r);
                log_nfa = log_nfa_new;
            }
        }
    }
    if(log_nfa > LOG_EPS) return log_nfa;

    // Try finer precision
    r = rect(rec);
    for(unsigned int n = 0; n < 5; ++n)
    {
        if((r.width - delta) >= 0.5)
        {
            r.p /= 2;
            r.prec = r.p * CV_PI;
            double log_nfa_new = rect_nfa(r);
            if(log_nfa_new > log_nfa)
            {
                rec = rect(r);
                log_nfa = log_nfa_new;
            }
        }
    }

    return log_nfa;
}

double LineSegmentDetectorImpl::rect_nfa(const rect& rec) const
{
    int total_pts = 0, alg_pts = 0;
    double half_width = rec.width / 2.0;
    double dyhw = rec.dy * half_width;
    double dxhw = rec.dx * half_width;

    edge ordered_x[4];
    edge* min_y = &ordered_x[0];
    edge* max_y = &ordered_x[0]; // Will be used for loop range

    /* build list of rectangle corners ordered
       in a circular way around the rectangle */
    ordered_x[0].p.x = int(rec.x1 - dyhw); ordered_x[0].p.y = int(rec.y1 + dxhw); ordered_x[0].taken = false;
    ordered_x[1].p.x = int(rec.x2 - dyhw); ordered_x[1].p.y = int(rec.y2 + dxhw); ordered_x[1].taken = false;
    ordered_x[2].p.x = int(rec.x2 + dyhw); ordered_x[2].p.y = int(rec.y2 - dxhw); ordered_x[2].taken = false;
    ordered_x[3].p.x = int(rec.x1 + dyhw); ordered_x[3].p.y = int(rec.y1 - dxhw); ordered_x[3].taken = false;

    std::sort(ordered_x, ordered_x + 4, AsmallerB_XoverY);//4个顶点,按照x从小到大排序

    // Find min y. And mark as taken. find max y.
    for(unsigned int i = 1; i < 4; ++i)
    {
        if(min_y->p.y > ordered_x[i].p.y) {min_y = &ordered_x[i]; }
        if(max_y->p.y < ordered_x[i].p.y) {max_y = &ordered_x[i]; }
    }
    min_y->taken = true;

    // Find leftmost untaken point;
    edge* leftmost = 0;
    for(unsigned int i = 0; i < 4; ++i)
    {
        if(!ordered_x[i].taken)
        {
            if(!leftmost) // if uninitialized
            {
                leftmost = &ordered_x[i];
            }
            else if (leftmost->p.x > ordered_x[i].p.x)
            {
                leftmost = &ordered_x[i];
            }
        }
    }
    CV_Assert(leftmost != NULL);
    leftmost->taken = true;

    // Find rightmost untaken point;
    edge* rightmost = 0;
    for(unsigned int i = 0; i < 4; ++i)
    {
        if(!ordered_x[i].taken)
        {
            if(!rightmost) // if uninitialized
            {
                rightmost = &ordered_x[i];
            }
            else if (rightmost->p.x < ordered_x[i].p.x)
            {
                rightmost = &ordered_x[i];
            }
        }
    }
    CV_Assert(rightmost != NULL);
    rightmost->taken = true;

    // Find last untaken point;
    edge* tailp = 0;
    for(unsigned int i = 0; i < 4; ++i)
    {
        if(!ordered_x[i].taken)
        {
            if(!tailp) // if uninitialized
            {
                tailp = &ordered_x[i];
            }
            //一共4个点,经过上面的步骤后,3个点 taken已经被标记为true,这里仅剩一个点未被标记
            //这里的else if应无必要!!
            else if (tailp->p.x > ordered_x[i].p.x)
            {
                tailp = &ordered_x[i];
            }
        }
    }
    CV_Assert(tailp != NULL);
    tailp->taken = true;

    double flstep = (min_y->p.y != leftmost->p.y) ?
                    (min_y->p.x - leftmost->p.x) / (min_y->p.y - leftmost->p.y) : 0; //first left step
    double slstep = (leftmost->p.y != tailp->p.x) ?
                    (leftmost->p.x - tailp->p.x) / (leftmost->p.y - tailp->p.x) : 0; //second left step

    double frstep = (min_y->p.y != rightmost->p.y) ?
                    (min_y->p.x - rightmost->p.x) / (min_y->p.y - rightmost->p.y) : 0; //first right step
    double srstep = (rightmost->p.y != tailp->p.x) ?
                    (rightmost->p.x - tailp->p.x) / (rightmost->p.y - tailp->p.x) : 0; //second right step

    double lstep = flstep, rstep = frstep;

    double left_x = min_y->p.x, right_x = min_y->p.x;

    // Loop around all points in the region and count those that are aligned.
    int min_iter = min_y->p.y;
    int max_iter = max_y->p.y;
    for(int y = min_iter; y <= max_iter; ++y)
    {
        if (y < 0 || y >= img_height) continue;

        for(int x = int(left_x); x <= int(right_x); ++x)
        {
            if (x < 0 || x >= img_width) continue;

            ++total_pts;/* total number of pixels counter */
            if(isAligned(x, y, rec.theta, rec.prec))
            {
                ++alg_pts;/* aligned points counter */
            }
        }

        if(y >= leftmost->p.y) { lstep = slstep; }
        if(y >= rightmost->p.y) { rstep = srstep; }

        left_x += lstep;
        right_x += rstep;
    }

    return nfa(total_pts, alg_pts, rec.p);
}

double LineSegmentDetectorImpl::nfa(const int& n, const int& k, const double& p) const
{
    // Trivial cases
    if(n == 0 || k == 0) { return -LOG_NT; }
    if(n == k) { return -LOG_NT - double(n) * log10(p); }

    /* probability term */
    double p_term = p / (1 - p);

    /* compute the first term of the series */
    /*
       binomial_tail(n,k,p) = sum_{i=k}^n bincoef(n,i) * p^i * (1-p)^{n-i}
       where bincoef(n,i) are the binomial coefficients.
       But
         bincoef(n,k) = gamma(n+1) / ( gamma(k+1) * gamma(n-k+1) ).
       We use this to compute the first term. Actually the log of it.
     */
    double log1term = (double(n) + 1) - log_gamma(double(k) + 1)
                      - log_gamma(double(n-k) + 1)
                      + double(k) * log(p) + double(n-k) * log(1.0 - p);
    double term = exp(log1term);

    //in some cases no more computations are needed
    if(double_equal(term, 0))/* the first term is almost zero */
    {

        if(k > n * p)  /* at begin or end of the tail?  */
            return -log1term / M_LN10 - LOG_NT;/* end: use just the first term  */
        else
            return -LOG_NT;/* begin: the tail is roughly 1  */
    }

    // Compute more terms if needed
    double bin_tail = term;
    double tolerance = 0.1; // an error of 10% in the result is accepted
    for(int i = k + 1; i <= n; ++i)
    {
        /*
           As
             term_i = bincoef(n,i) * p^i * (1-p)^(n-i)
           and
             bincoef(n,i)/bincoef(n,i-1) = n-1+1 / i,
           then,
             term_i / term_i-1 = (n-i+1)/i * p/(1-p)
           and
             term_i = term_i-1 * (n-i+1)/i * p/(1-p).
           1/i is stored in a table as they are computed,
           because divisions are expensive.
           p/(1-p) is computed only once and stored in 'p_term'.
         */
        double bin_term = double(n - i + 1) / double(i);
        double mult_term = bin_term * p_term;
        term *= mult_term;
        bin_tail += term;
        if(bin_term < 1)
        {
            /* When bin_term<1 then mult_term_j<mult_term_i for j>i.
               Then, the error on the binomial tail when truncated at
               the i term can be bounded by a geometric series of form
               term_i * sum mult_term_i^j.                            */
            double err = term * ((1 - pow(mult_term, double(n-i+1))) / (1 - mult_term) - 1);
            /* One wants an error at most of tolerance*final_result, or:
              tolerance * abs(-log10(bin_tail)-logNT).
              Now, the error that can be accepted on bin_tail is
              given by tolerance*final_result divided by the derivative
              of -log10(x) when x=bin_tail. that is:
              tolerance * abs(-log10(bin_tail)-logNT) / (1/bin_tail)
              Finally, we truncate the tail if the error is less than:
              tolerance * abs(-log10(bin_tail)-logNT) * bin_tail        */
            if(err < tolerance * fabs(-log10(bin_tail) - LOG_NT) * bin_tail) break;
        }

    }
    return -log10(bin_tail) - LOG_NT;
}

inline bool LineSegmentDetectorImpl::isAligned(int x, int y, const double& theta, const double& prec) const
{
    if(x < 0 || y < 0 || x >= angles.cols || y >= angles.rows) { return false; }
    const double& a = angles.at<double>(y, x);//候选点的梯度幅角
    if(a == NOTDEF) { return false; }

    // It is assumed that 'theta' and 'a' are in the range [-pi,pi]
    double n_theta = theta - a;//种子点梯度幅角-候选点梯度幅角
    if(n_theta < 0) { n_theta = -n_theta; }
    if(n_theta > M_3_2_PI)
    {
        n_theta -= M_2__PI;
        if(n_theta < 0) n_theta = -n_theta;
    }

    return n_theta <= prec;
}


void LineSegmentDetectorImpl::drawSegments(cv::InputOutputArray _image, cv::InputArray lines)
{

    CV_Assert(!_image.empty() && (_image.channels() == 1 || _image.channels() == 3));

    if (_image.channels() == 1)
    {
        cvtColor(_image, _image, cv::COLOR_GRAY2BGR);
    }

    cv::Mat _lines = lines.getMat();
    const int N = _lines.checkVector(4);

    CV_Assert(_lines.depth() == CV_32F || _lines.depth() == CV_32S);

    cv::RNG rng(cv::getTickCount());


    // Draw segments
    if (_lines.depth() == CV_32F)
    {
        for (int i = 0; i < N; ++i)
        {
            const cv::Vec4f& v = _lines.at<cv::Vec4f>(i);
            const cv::Point2f b(v[0], v[1]);
            const cv::Point2f e(v[2], v[3]);

            //auto color = cv::Scalar(rng.uniform(0,255),rng.uniform(0,255),rng.uniform(0,255));
            cv::Scalar color(0,0,255);
            line(_image, b, e, color, 1);
            cv::circle(_image,b,2,color,-1);
        }
    }
    else
    {
        for (int i = 0; i < N; ++i)
        {
            const cv::Vec4i& v = _lines.at<cv::Vec4i>(i);
            const cv::Point2i b(v[0], v[1]);
            const cv::Point2i e(v[2], v[3]);

            //auto color = cv::Scalar(rng.uniform(0,255),rng.uniform(0,255),rng.uniform(0,255));
            cv::Scalar color(0,0,255);
            line(_image, b, e, color, 1);
            cv::circle(_image,b,2,color,-1);
        }
    }
}


int LineSegmentDetectorImpl::compareSegments(const cv::Size& size, cv::InputArray lines1, cv::InputArray lines2, cv::InputOutputArray _image)
{

    cv::Size sz = size;
    if (_image.needed() && _image.size() != size) sz = _image.size();
    CV_Assert(!sz.empty());

    cv::Mat_<uchar> I1 = cv::Mat_<uchar>::zeros(sz);
    cv::Mat_<uchar> I2 = cv::Mat_<uchar>::zeros(sz);

    cv::Mat _lines1 = lines1.getMat();
    cv::Mat _lines2 = lines2.getMat();
    const int N1 = _lines1.checkVector(4);
    const int N2 = _lines2.checkVector(4);

    CV_Assert(_lines1.depth() == CV_32F || _lines1.depth() == CV_32S);
    CV_Assert(_lines2.depth() == CV_32F || _lines2.depth() == CV_32S);

    if (_lines1.depth() == CV_32S)
        _lines1.convertTo(_lines1, CV_32F);
    if (_lines2.depth() == CV_32S)
        _lines2.convertTo(_lines2, CV_32F);

    // Draw segments
    for(int i = 0; i < N1; ++i)
    {
        const cv::Point2f b(_lines1.at<cv::Vec4f>(i)[0], _lines1.at<cv::Vec4f>(i)[1]);
        const cv::Point2f e(_lines1.at<cv::Vec4f>(i)[2], _lines1.at<cv::Vec4f>(i)[3]);
        line(I1, b, e, cv::Scalar::all(255), 1);
    }
    for(int i = 0; i < N2; ++i)
    {
        const cv::Point2f b(_lines2.at<cv::Vec4f>(i)[0], _lines2.at<cv::Vec4f>(i)[1]);
        const cv::Point2f e(_lines2.at<cv::Vec4f>(i)[2], _lines2.at<cv::Vec4f>(i)[3]);
        line(I2, b, e, cv::Scalar::all(255), 1);
    }

    // Count the pixels that don't agree
    cv::Mat Ixor;
    bitwise_xor(I1, I2, Ixor);
    int N = countNonZero(Ixor);

    if (_image.needed())
    {
        CV_Assert(_image.channels() == 3);
        cv::Mat img = _image.getMatRef();
        CV_Assert(img.isContinuous() && I1.isContinuous() && I2.isContinuous());

        for (unsigned int i = 0; i < I1.total(); ++i)
        {
            uchar i1 = I1.ptr()[i];
            uchar i2 = I2.ptr()[i];
            if (i1 || i2)
            {
                unsigned int base_idx = i * 3;
                if (i1) img.ptr()[base_idx] = 255;
                else img.ptr()[base_idx] = 0;
                img.ptr()[base_idx + 1] = 0;
                if (i2) img.ptr()[base_idx + 2] = 255;
                else img.ptr()[base_idx + 2] = 0;
            }
        }
    }

    return N;
}

调用示例:

//
// Created by liheng on 12/26/21.
//

#include <opencv2/opencv.hpp>
#include "lsd.h"

int main()
{
    auto b = DBL_EPSILON;
    auto a=100*DBL_EPSILON;


    cv::Mat image = cv::imread("../Pictures/Door_Camera_000001.bmp",cv::IMREAD_GRAYSCALE);
    image = image(cv::Rect(640,360,640,360));
    cv::Mat imageWithLines = image.clone();



    cv::Canny(image,image,50,100,3,false);
    cv::imshow("canny",image);
    auto pLSD = createlLineSegmentDetector(cv::LSD_REFINE_ADV);
    std::vector<cv::Vec4f> lines_std;

    pLSD->detect(image,lines_std);

    pLSD->drawSegments(imageWithLines,lines_std);
    cv::imshow("imageWithLines",imageWithLines);

    cv::waitKey(0);


    return 0;
}

Canny和LSD运行结果如下：

参考资料

1.LSD论文
2.LSD:一种直线检测算法简介
3.直线段检测算法—LSD：a Line Segment Detector
4.直线段检测算法（LSD：a Line Segment Detector）
5.基于C+OpenCV4.0的LineSegmentDetector算法实现
6.图像LSD直线检测
7.LSD直线检测算法
8.基于LSD的直线提取算法

最先接触的图像算法就是LSD直线检测和FLD直线检测来实现车道线检测,一直希望能够看下LSD直线检测的源码，顺便整理下其实现流程，以便于对其进行修改和扩展，不过总是断断续续、零散的记录，终于趁着几个周末，将之前未完成的内容给补充下去。

---