基于LSD的直线提取算法

LSD是一种局部提取直线的算法，速度比Hough要快。

但是有局部算法的缺点：

1.对于直线相交情况，因为设置了每个点是否USED，因此每个点只能属于一条直线，若有相交必有至少一条直线被割裂为两条。又因为其基于梯度，直线交点梯度值往往又较小（不被检测为边缘点），因此很有可能相交的两条直线在交点处被割裂为四条线段。

2.由于局部检测算法自增长的特点，对于长线段被遮挡、局部模糊等原因经常割裂为多条直线。这些缺点在Hough变换中不存在。

算法步骤概括：

1.以 s=0.8的尺度对输入图像进行高斯下采样。

2.计算每一个点的梯度值以及梯度方向（level-line orientation）。

3.根据梯度值对所有点进行伪排序（pseudo-ordered），建立状态列表，所有点设置为UNUSED。

4.将梯度值小于ρ的点状态表中相应位置设置为USED。

5.取出列表中梯度最大（伪排列的首位）的点作为种子点（seed），状态列表中设为USED。

   do：

    a.以seed为起点，搜索周围UNUSED并且方向在阈值[ -t, t]范围内的点，状态改为USED。

    b.生成包含所有满足点的矩形R。

    c.判断同性点（aligned pt）密度是否满足阈值D，若不满足，截断（cut）R变为多个矩形框，直至满足。

    d.计算NFA。

    e.改变R使NFA的值更小直至NFA <= ε ，R加入输出列表。

伪代码说明：

// 带有尺度参数和区域输出的 LSD 接口函数
double * lsd_scale_region( int * n_out,double * img, int X, int Y, double scale,
			int ** reg_img, int * reg_x, int * reg_y )
{
	double sigma_scale = 0.6  ; // Sigma for Gaussian filter is computed as sigma = sigma_scale/scale.
	double quant       = 2.0  ; // Bound to the quantization error on the gradient norm.
	double ang_th      = 22.5 ; // Gradient angle tolerance in degrees.
	double log_eps     = 0.0  ; // Detection threshold: -log10(NFA) > log_eps
	double density_th  = 0.7  ; // Minimal density of region points in rectangle.
	int    n_bins      = 1024 ; // Number of bins in pseudo-ordering of gradient modulus.

	return LineSegmentDetection( n_out, img, X, Y, scale, sigma_scale, quant,
		                    ang_th, log_eps, density_th, n_bins,
		                    reg_img, reg_x, reg_y );
}

<p>// LSD 完整的接口函数
double * LineSegmentDetection( int * n_out,
         double * img, int X, int Y,
         double scale, double sigma_scale, double quant,
         double ang_th, double log_eps, double density_th,
         int n_bins,
         int ** reg_img, int * reg_x, int * reg_y )
{
 int ls_count = 0;     /* line segments are numbered 1,2,3,... */</p><p>
  // 角度阈值
 prec = M_PI * ang_th / 180.0;
 p = ang_th / 180.0;  // contrario model中角度符合level-line角度阈值的概率
 rho = quant / sin(prec); // 梯度阈值</p><p> // 对图像进行尺度变换
 image = new_image_double_ptr( (unsigned int) X, (unsigned int) Y, img );
 if( scale != 1.0 )
 {
  scaled_image = gaussian_sampler( image, scale, sigma_scale );
  // 计算每个点的方向
  angles = ll_angle( scaled_image, rho, &list_p, &mem_p, &modgrad, (unsigned int) n_bins );
  free_image_double(scaled_image);
 }
 else
  // 计算每个点的方向
  angles = ll_angle( image, rho, &list_p, &mem_p, &modgrad, (unsigned int) n_bins );
 xsize = angles->xsize;
 ysize = angles->ysize;</p><p> /* Number of Tests - NT</p><p> improve rectangle step根据p的衰减看出运行11次，共(XY)^(5/2)种线段，得到Ntest</p><p> The theoretical number of tests is Np.(XY)^(5/2) where X and Y are number of columns and rows of the image.
 Np corresponds to the number of angle precisions considered. As the procedure 'rect_improve' tests 5 times
 to halve the angle precision, and 5 more times after improving other factors, 11 different precision values
 are potentially tested. Thus, the number of tests is 
   11 * (X*Y)^(5/2)
 whose logarithm value is
   log10(11) + 5/2 * (log10(X) + log10(Y)).                                       */
 logNT = 5.0 * ( log10( (double) xsize ) + log10( (double) ysize ) ) / 2.0 + log10(11.0);
 min_reg_size = (int) (-logNT/log10(p)); /* minimal number of points in region that can give a meaningful event */</p><p> // 搜索线段
 for(; list_p != NULL; list_p = list_p->next )
 {
  if( used->data[ list_p->x + list_p->y * used->xsize ] == NOTUSED &&
   angles->data[ list_p->x + list_p->y * angles->xsize ] != NOTDEF )
   /* there is no risk of double comparison problems here
   because we are only interested in the exact NOTDEF value */
  {
   // 找出位置相邻接且方向相近的点构成的区域
   region_grow( list_p->x, list_p->y, angles, reg, &reg_size, &reg_angle, used, prec );</p><p>   // 清除较小的区域
   if( reg_size < min_reg_size ) continue;</p><p>   // 构造矩形，近似该区域
   region2rect(reg,reg_size,modgrad,reg_angle,prec,p,&rec);</p><p>   // 检查矩形中点的密度，若小于阈值，则改进矩形表示。若仍旧不满足阈值条件，则清除该区域
   if( !refine( reg, &reg_size, modgrad, reg_angle,
    prec, p, &rec, used, angles, density_th ) ) continue;</p><p>   // 计算 NFA 值
   log_nfa = rect_improve(&rec,angles,logNT,log_eps);
   if( log_nfa <= log_eps ) continue;</p><p>   // 找到一条新的线段
   ++ ls_count;  // 线段数目 増一</p><p>   // 补偿梯度计算时的误差
   rec.x1 += 0.5; rec.y1 += 0.5;
   rec.x2 += 0.5; rec.y2 += 0.5;</p><p>   // 尺度变换
   if( scale != 1.0 )
   {
    rec.x1 /= scale; rec.y1 /= scale;
    rec.x2 /= scale; rec.y2 /= scale;
    rec.width /= scale;
   }</p><p>   // 将找到的直线加入到输出变量中
   add_7tuple( out, rec.x1, rec.y1, rec.x2, rec.y2,
    rec.width, rec.p, log_nfa );</p><p>   // 在 'region' 图像中标记出该直线
   if( region != NULL )
    for(i=0; i<reg_size; i++)
     region->data[ reg[i].x + reg[i].y * region->xsize ] = ls_count;
  }
 }</p><p>  // 返回结果
  return return_value;
}</p>

1.尺度缩放:

按照原文，这一步可选择跳过（s =1.0）。算法的默认选择是s=0.8。目的是为了消除锯齿效应。然后用高斯下采样的方式对输入图像进行操作，其中σ =∑ /s。∑=0.6。

2.梯度和方向计算：

梯度计算中论文用到的是2*2的模板。主要是为了用contrario model的时候保证相邻点方向分布的独立性。

gx(x,y) = ( i(x + 1, y) + i(x + 1, y + 1) − i(x, y) − i(x, y + 1)) / 2

gy(x, y) =( i(x, y + 1) + i(x + 1, y + 1) − i(x, y) − i(x + 1, y)) / 2

规范化后的梯度：

g(x,y) = (gx^2 + gy^2)^0.5

梯度方向的计算：
ang(x,y) = arctan（gx(x, y) / （−gy(x, y））

这种方法计算的梯度其实是（x+0.5，y+0.5）处的，在最后的步骤会对这里的计算偏差进行补偿。

3.梯度伪排序： 

梯度值越大，越是显著的边缘点，因此更适合作为种子点。但是对梯度值进行完全排序是一个时耗性很高的工作。因此简单的将梯度值划分为1024个等级（bins），这1024个等级涵盖了梯度由0~255的变化范围，这种排序是一个线性的时耗。种子点从梯度值最高的bin开始搜索，依次往下，直至所有点标记为USED。

小梯度值抑制：

小梯度值点往往出现在平滑区域，或者仅仅是噪声。不在关注的范围内，但是他们的存在往往会严重影响直线角度的计算。考虑到有角度承受范围[-t,t]。因此小梯度“脏点”的最大影响角度需要小于t。假设“脏点”的阈值为q。那么|bad ang|=q / sin（t）。这里，t被坐着设置为22.5度。因此在这里q=2。

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

由种子点搜索角度满足并且状态为UNUSED的点（八邻域）形成的区域称为line-support region。邻域内在角度承受范围t内满足整个区域的方向角region-ang的点加入到这个区域。每次增加一个点，

reigon-ang会更新：

region-ang=arctan（∑（sin（pt(i)）/ ∑（cos（pt(i)）），i为region内所有的点。 

直线（矩形）勾勒：

上一步后line-support region一系列相邻的离散点，需要将他们包含在一个矩形框内R，这个可以看做宽度为R的宽，长度为R的长的候选直线。R内，上一步所选择的点即为类内点（aligned pt），显然，矩形框内还很有可能包含很多上一步没有选择近line-support region的outline pt。R选择能包含line-support region的最小矩形，所有点的梯度规范化值平均计算重心。R长轴的方向设置为R的方向。

5.NFA的计算： 

The validation step is based on the a contrario approachandthe Helmholtz principle proposed by Desolneux.

the Helmholtz principle：在完美噪声图像图像中不应该检测到目标。

       contrario approach：一个不会检测到目标的噪声图像。对于本课题，contrario model是一个像素值为level-line angle的图像。其level-line angle随机分解为独立且服从平均分布于[0-2*PI]。

这里用NFA（Number of False Alarms）来评判observe img中某个候选Rec少于contrario model中相同位置Rec里aligned pt点的数量的概率，NFA越大，表明当前Rec与contrario model中相同位置越相似，相反的，当前rec越有可能是“真正的目标”。

其中:

 NFA = N * Ph0[ k(r,I) >= k(r,i) ]

上式中，N为当前大小（m*n）图像中直线（矩形框）的数量。在m*n的图像中直线的起点和终点分别有m*n种选择，所以一共有(m*n）*(m*n)种起点和终点搭配。线段的宽度为[0,m*n ^0.5]。因此在m*n大小的图像中有（m*n）^2.5 种不同直线。

N=（m*n）^2.5

k(r,I) 为contrario model ，I 中 r 矩形里aligned pt的数量。

k(r,i) 为observe img，i 中 r 矩形里aligned pt的数量。

根据contrario model的建立规则，每个像素都服从值域为[0,2*PI]的二项分布。

设k(r,i)=k。

其中：

角度正负容忍误差为t，总容忍误差为2t。那么在contrario model中某个点为aligned pt的概率为 p=2*t / 2*Pi = t / PI。

那么，在 I 中的 r 矩形里，总像素个数为 n。I 中的 r 矩形里aligned pt个数k(r,I)大于等于k的话，可选择的值为k(r,i)、k(r,i)+1、k(r,i)......n。

概率为:   ∑ p^i * (1-p)^(n-i) ，i=k，k+1，k+2……n。

因此：

NFA = （m*n）^2.5  *  ∑ p^i * (1-p)^(n-i) ，i=k，k+1，k+2……n

如果NFA的值很大，认为当前事件出现在contrario model的概率很大，将其认为是背景中的一部分。相反的，认为目标是相对突出（rare）的，是一个合适的“直线”。

即：NFA <= ε     

ε-meaningful rectangles :

E(ho)[ ∑NFA(r) <= ε ] <= ε ,  r∈ R。R为图像中所有直线的集合。E为期望运算。ho为contrario model。

根据论文，ε =1

 类内点（aligned pt）密度计算

fig 1

 如fig 1所示，因为t的存在且有正负号容忍范围，容易出现图中line-support region增长不合理，此时应将R截断为两个更小的R更为合适。若如图分割的区域包含的直线形成的角为钝角，本次操作只留下包含seed的一段，舍弃另一段。舍弃的点重置为UNUSED。

6.类内点密度： 

d = k / n。k为类内点个数，n为R的length*width。若d > D。accepted。否则，需要将R截断。在这里设置D=0.7。文章认为这个数字既能保证同一个R中的类内点属性相近，也能保证R不会被过分的分割为小的矩形。

若没满足d > D 分割方法分为“缩小角度容忍阈值”与“缩小区域半径”的方法。

缩小角度容忍阈值：简单的将t值缩小，再次从当前R的seed开始搜寻。如果这一步仍为满足d > D，将逐步舍弃离seed较远的点，直至满足不等式未知。

类内点的计算将会解决曲线因为t的存在而误识别为直线的状况。

7.Rectangle Improvement： 

如果当前的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。

 Results：

  

Reference:

1. <LSD :A Fast Line Segment Detector with a False Detection Control>. Rafael Grompone von Gioi.etc

2.<LSD: a Line Segment Detector>.Rafael Grompone von Gioi.etc

代码下载：
http://www.ipol.im/pub/art/2012/gjmr-lsd/