主要流程
分为三个步骤
1、颜色通道分离,引导线(红色)颜色提取(详情参考前面的博文);
2、霍夫直线检测;
3、霍夫直线聚类;
霍夫直线
检测函数
- 第一个检测函数
void HoughLines(InputArray image, OutputArray lines,
double rho, double theta, int threshold,
double srn=0, double stn=0,
double min_theta = 0, double max_theta = CV_PI )
参数解释:
第一个参数:InputArray类型的image,输入图像,即源图像,需为8位的单通道二进制图像,可以将任意的源图载入进来后由函数修改成此格式后,再填在这里。
第二个参数:InputArray类型的lines,经过调用HoughLines函数后储存了霍夫线变换检测到线条的输出矢量。每一条线由具有两个元素的矢量表示,其中,是离坐标原点((0,0)(也就是图像的左上角)的距离。 是弧度线条旋转角度(0垂直线,π/2水平线)。
第三个参数:double类型的rho,以像素为单位的距离精度。另一种形容方式是直线搜索时的进步尺寸的单位半径。PS:Latex中/rho就表示 。
第四个参数:double类型的theta,以弧度为单位的角度精度。另一种形容方式是直线搜索时的进步尺寸的单位角度。
第五个参数:int类型的threshold,累加平面的阈值参数,即识别某部分为图中的一条直线时它在累加平面中必须达到的值。大于阈值threshold的线段才可以被检测通过并返回到结果中。
第六个参数:double类型的srn,有默认值0。对于多尺度的霍夫变换,这是第三个参数进步尺寸rho的除数距离。粗略的累加器进步尺寸直接是第三个参数rho,而精确的累加器进步尺寸为rho/srn。
第七个参数:double类型的stn,有默认值0,对于多尺度霍夫变换,srn表示第四个参数进步尺寸的单位角度theta的除数距离。且如果srn和stn同时为0,就表示使用经典的霍夫变换。否则,这两个参数应该都为正数。
第八个参数:double类型的 min_theta,对于标准和多尺度Hough变换,检查线条的最小角度。必须介于0和max_theta之间。
第九个参数:double类型的 max_theta, 对于标准和多尺度Hough变换,检查线条的最大角度。必须介于min_theta和CV_PI之间.
- 第二个检测函数
直线描述
霍夫直线参数含义
霍夫直线参数位rho,theta。
其中,
rho代表图像原点至直线的距离;
theta代表直线的角度,垂直为0度,顺时针旋转逐渐增大,水平位90度。
theta的取值范围:大于等于0度, 小于180度。
数学方程表示直线
已知: rho,theta;
直线系数:a = cos(theta), b = sin(theta);
垂点坐标:x0 = a * rho, y0 = b * rho;
直线方程:x = x0 - bt, y = y0 + at;
直线反侧变换
标准直线参数的theta范围位[0,PI),在直线之间进行运算时会用到theta在此范围外的值,例如将1度的直线变换为181度,或者将179度的直线变换为-1度,这种变换成为反侧变换,反侧变换前后的直线等效,为同一直线。
- 反侧变换一:
new_theta = old_theta - PI;
new_rho = - old_rho;
- 反侧变换二:
new_theta = old_theta + PI;
new_rho = - old_rho;
直线规范化操作
将直线参数的角度值控制在[0,PI)范围内。
while (center_theta > PI)
{
center_theta -= PI;
center_rho = - center_rho;
}
while (center_theta < 0)
{
center_theta += PI;
center_rho = - center_rho;
}
霍夫直线间求差
在霍夫直线聚类时,需要将直线参数看作二维点,在此计算二维点之间的距离。
设两直线参数为(rho1, theta1)、(rho2, theta2),求取绝对误差:
error_theta = abs(theta1-theta2);
if error_theta <= PI/2 时:
error_rho = abs(rho1 - rho2);
if error_theta > PI/2 时:
error_theta = PI - error_theta;
error_rho = abs(rho1 + rho2);
计算归一化直线误差(距离):
d0 = drho/Rho_max;
d1 = dtheta/(CV_PI*0.5);
error_cj = sqrtf(d0*d0 + d1*d1);
霍夫直线求和与均值
在霍夫直线聚类求取中心点时,需要计算多直线参数累加和、均值。
- 计算累加和:
当样本直线距离该聚类中心直线的角度差大于90度时,需要做反侧变换,伪代码描述过程如下。
sum_theta = 0;
sum_rho = 0;
this_cluster_num = 0
for every lines in this cluster center:
this_cluster_num ++;
# 直线角度距离大于90度,则做反侧变换
if (this_theta - center_theta) > PI/2:
this_theta = this_theta - PI;
this_rho = - this_rho;
else if (this_theta - center_theta) < - PI/2:
this_theta = this_theta + PI;
this_rho = - this_rho;
sum_theta += this_theta;
sum_rho += this_rho;
- 求取均值(中心点)
center_theta = sum_theta/this_cluster_num;
center_rho = sum_rho/this_cluster_num;
- 结果规范化操作
为了保持中心直线参数在多伦的聚类迭代运算时,始终保持在一定的范围内,需要对其进行常规化操作,将其范围控制在[0,PI)之间。
霍夫直线滤波
在直线跟踪时,为了保证目标直线的平滑行,需要加入滤波处理。
普通一阶低通滤波的方法: y = y_1 + k*(x - y_1);
霍夫直线滤波时需要考虑直线反侧变换,过程用伪代码表示:
# 偏差检测与反侧变换
error_theta = x_theta - y_theta;
if (error_theta > PI/2):
# 直线反侧变换
x_u_theta = x_theta - PI;
x_u_rho = - x_rho;
else if (error_theta < -PI/2):
# 直线反侧变换
x_u_theta = x_theta + PI;
x_u_rho = - x_rho;
else:
x_u_theta = x_theta;
x_u_rho = -x_rho;
# 滤波处理
y_line = y_line + k*(x_u_line - y_line);
霍夫直线聚类
直线聚类的算法流程
聚类中心数量cluster_center_num = 1;
while 聚类未完成:
for 聚类N次数:
聚类初始化工作;
随机分布聚类中心点:在数据样本集中随机选取cluster_center_num 个点,作为数据中心点;
while 本次聚类未达到收敛状态:
根据当前聚类中心值,计算每一个点归属于哪一个聚类中心;
根据归属关系,计算新的聚类中心值;
当新的聚类中心与上次聚类中心一致时,本次聚类达到收敛状态;
计算本次聚类完成后的样本点均方差,并记录;
选取N次聚类结果中,均方差最小的那一次数据,作为聚类结果。
若聚类均方差小于设定阈值,则聚类完成;否则,cluster_center_num ++,开启新的N轮聚类。
随机分布聚类中心点的方法:
(1)在取值空间中随机取值坐标点,产生cluster_center_num 个数据点;
(2)在数据样本集中随机选取cluster_center_num 个点,作为数据中心点;
如果样本点在空间中分布非常不均匀的情况下,方法(1)产生的数据中心点有些偏离数据较远,则很难达到较好的收敛状态,即某些聚类中心可能没有隶属样本点。
需要考虑的问题:
- 聚类的总类别个数不确定
因为事先不知道图形中存在多少真实的直线,所以要自适应寻找聚类的类别数,判断依据是聚类后的均方差是否足够小; - 霍夫直线参数误差计算时的特殊处理
(1) 角度处理
由霍夫直线参数物理含义可知,1度的直线与179度的直线,在斜率上仅仅差别2度。所以,
当两直线的误差error = | theta1 - theta2 |大于90度时,
其真实误差是| PI - error|。
(2) 距离处理
当两个直线的角度偏差计算经过折算处理时,其斜率计算也需要进行折算处理。
效果展示
图1
原图
霍夫直线检测与聚类
原图添加效果
图2
原图
霍夫直线检测与聚类
原图添加效果
图3
原图
霍夫直线检测与聚类
原图添加效果
程序示例
程序主体
#include <opencv2/opencv.hpp>
#include <opencv2/xfeatures2d.hpp>
#include <iostream>
using namespace cv;
using namespace cv::xfeatures2d; // 不要忘了导入扩展模块
using namespace std;
Mat src_img, gray_img;
const string output_name = "二值化";
int minHessian = 100; // 定义SURF中的hessian阈值特征点检测算子
int max_value = 250;
Mat bool_img;
// =================================================
// 颜色通道分离 + 霍夫直线检测 + 直线绘制
// =================================================
int Rho_max = 500;
void cluster_lines(vector<Vec2f> lines, vector<Vec2f> *lines_center);
string Name = "6";
int main()
{
// Mat src = imread("/home/parobot/Pictures/rgb.jpeg");
// resize(src, src_img, Size(round(1.0*src.cols), round(1.0*src.rows)));
Mat src = imread("/home/parobot/Pictures/line/" + Name + ".jpg");
resize(src, src_img, Size(round(0.15*src.cols), round(0.15*src.rows)));
// resize(src, src_img, Size(round(1.0*src.cols), round(1.0*src.rows)));
if (src_img.empty())
{
printf("could not load the image...\n");
return -1;
}
namedWindow("原图", WINDOW_AUTOSIZE);
imshow("原图", src_img);
cvtColor(src_img, gray_img, COLOR_BGR2GRAY);
vector<Mat> channels;
split(src_img, channels);
// imshow("red", channels.at(2));
// imshow("green", channels.at(1));
// imshow("blue", channels.at(0));
cout << channels.at(0).size << endl;
Rho_max = sqrtf(src_img.rows*src_img.rows + src_img.cols*src_img.cols);
for (int i = 0; i < channels.at(0).rows; i ++)
{
for (int j = 0; j < channels.at(0).cols; j ++)
{
// 提取红色
// 红色通道值,应该大于其它通道值
int red = channels.at(2).at<uchar>(i,j) - 30;
if ((red < channels.at(0).at<uchar>(i,j))
||(red < channels.at(1).at<uchar>(i,j)))
{
channels.at(0).at<uchar>(i,j) = 0;
channels.at(1).at<uchar>(i,j) = 0;
channels.at(2).at<uchar>(i,j) = 0;
}
}
}
Mat red_line;
merge(channels, red_line);
// 霍夫直线检测
vector<Vec2f> lines;
HoughLines(channels.at(2), lines, 1, CV_PI/180, 200, 0, 0);
cout << lines.size() << endl;
for (int i = 0; i < lines.size(); i ++)
{
float rho = lines[i][0], theta = lines[i][1];
Point p1, p2;
float a = cos(theta), b = sin(theta);
float x0 = a*rho, y0 = b*rho;
p1.x = cvRound(x0 + 1000*(-b));
p1.y = cvRound(y0 + 1000*a);
p2.x = cvRound(x0 - 1000*(-b));
p2.y = cvRound(y0 - 1000*a);
line(red_line, p1, p2, Scalar(255, 0, 0), 1, LINE_AA);
}
// 概率霍夫直线检测
// vector<Vec4i> lines;
// HoughLinesP(channels.at(2), lines, 1, CV_PI/180, 80, 50, 10);
// cout << lines.size() << endl;
// for (int i = 0; i < lines.size(); i ++)
// {
// Vec4i l = lines[i];
// line(red_line, Point(l[0],l[1]), Point(l[2], l[3]), Scalar(255, 0, 0), 1, LINE_AA);
// }
imwrite("/home/parobot/Pictures/line/orgin" + Name + ".jpg", src_img);
// 聚类
vector<Vec2f> lines_center;
cluster_lines(lines, &lines_center);
cout << "lines_center.size : " << lines_center.size() << endl;
for (int i = 0; i < lines_center.size(); i ++)
{
float rho = lines_center[i][0], theta = lines_center[i][1];
Point p1, p2;
float a = cos(theta), b = sin(theta);
float x0 = a*rho, y0 = b*rho;
p1.x = cvRound(x0 + 1000*(-b));
p1.y = cvRound(y0 + 1000*a);
p2.x = cvRound(x0 - 1000*(-b));
p2.y = cvRound(y0 - 1000*a);
line(red_line, p1, p2, Scalar(0, 255, 0), 2, LINE_AA);
line(src_img, p1, p2, Scalar(0, 255, 0), 2, LINE_AA);
}
imshow("cluster", red_line);
imshow("result", src_img);
imwrite("/home/parobot/Pictures/line/cluster" + Name + ".jpg", red_line);
imwrite("/home/parobot/Pictures/line/result" + Name + ".jpg", src_img);
waitKey(0);
return 0;
}
聚类算法
void cluster_lines(vector<Vec2f> lines, vector<Vec2f> *lines_center)
{
int cluster_num = 1;
bool cluster_done = false;
vector<Vec2f> center_this;
vector<Vec2f> center_best;
vector<vector<Vec2f>> center_episode;
center_episode.resize(3);
while (!cluster_done)
{
cout << " ================================================ " << endl;
cout << "新的一次聚类中心点数量, cluster_num = " << cluster_num << endl;
center_this.resize(cluster_num);
center_best.resize(cluster_num);
float error_min_rec = 0.0;
float error_points[lines.size()];
float error_center[cluster_num];
float error_center_sum_min;
// 多次聚类取最优
int episode_best_index = 0;
for (int episode = 0; episode < 3; episode ++)
{
cout << " ------------------------------------------ " << endl;
cout << "cluster begin, episode = " << episode << endl;
// 初始化:给中心点赋随机值
// for (int i = 0; i < cluster_num; i ++)
// {
// unsigned seed = time(0);
// srand(seed);
// int rho = rand()%(Rho_max + 1);
// int angle = rand()%180;
// if (angle > 90)
// {
// center_this[i][0] = -(float)rho;
// }
// else
// {
// center_this[i][0] = (float)rho;
// }
// center_this[i][1] = ((float)angle)*CV_PI/180.0;
// }
unsigned seed = time(0);
srand(seed);
for (int i = 0; i < cluster_num; i ++)
{
int index = rand()%lines.size();
center_this[i] = lines[index];
cout << "center_this[" << i << "]" << center_this[i] << endl;
}
// 循环多次,直到聚类中心收敛到稳态
bool center_stable = false;
int Point_Belong[lines.size()];
int Center_Points_Num[cluster_num];
//while (!center_stable)
for (int time = 0; time < 3; time ++)
{
cout << "update the cluster center ... ... ..." << endl;
// 新一轮聚类的初始化
for (int cj = 0; cj < cluster_num; cj ++)
{
Center_Points_Num[cj] = 0.0;
}
// 分类:遍历所有点,判断点归属于哪个中心
Vec2f Center_Points_Vaule_Sum[cluster_num];
for (int pi = 0; pi < lines.size(); pi ++)
{
float moni_d0_min , moni_d1_min;
float error_min = 1.0e10;
// 计算中心距离:遍历所有中心点
for (int cj = 0; cj < cluster_num; cj ++)
{
// 折算最近方向的距离
float dtheta = fabsf(lines[pi][1] - center_this[cj][1]);
float drho;
if (dtheta > CV_PI*0.5)
{
dtheta = fabsf(CV_PI - dtheta);
drho = fabsf(lines[pi][0] + center_this[cj][0]);
}
else
{
drho = fabsf(lines[pi][0] - center_this[cj][0]);
}
// 计算归一化距离
float d0 = drho/Rho_max;
float d1 = dtheta/(CV_PI*0.5);
float error_cj = sqrtf(d0*d0 + d1*d1);
if (error_min > error_cj)
{
error_min = error_cj;
Point_Belong[pi] = cj;
moni_d0_min = d0;
moni_d1_min = d1;
}
}
error_points[pi] = error_min;
cout << "d0 = " << moni_d0_min << ", d1 = " << moni_d1_min << endl;
// 计算累加和、个数
Center_Points_Num[Point_Belong[pi]] ++;
// 计算累加和:需要折算
Vec2f re_line = lines[pi];
float theta_err = lines[pi][1] - center_this[Point_Belong[pi]][1];
if (theta_err > CV_PI*0.5)
{
re_line[1] -= CV_PI;
re_line[0] = - lines[pi][0];
}
else if (theta_err < -CV_PI*0.5)
{
re_line[1] += CV_PI;
re_line[0] = - lines[pi][0];
}
Center_Points_Vaule_Sum[Point_Belong[pi]] += re_line;
// cout << "lines[" << pi << "] " << lines[pi] << endl;
// cout << "Center_Points_Vaule_Sum[" << Point_Belong[pi]<< "] " << Center_Points_Vaule_Sum[Point_Belong[pi]] << endl;
}
// 更新中心点坐标
bool center_change = false;
for (int cj = 0; cj < cluster_num; cj ++)
{
cout << " center_this[" << cj << "] " << center_this[cj] << endl;
cout << "Center_Points_Vaule_Sum[" << cj << "] " << Center_Points_Vaule_Sum[cj] << endl;
cout << " Center_Points_Num[" << cj << "] " << Center_Points_Num[cj] << endl;
Vec2f center_temp = Center_Points_Vaule_Sum[cj]/Center_Points_Num[cj];
if (center_this[cj] != center_temp)
{
center_change = true;
center_this[cj] = center_temp;
}
cout << " center_temp " << center_temp << endl;
}
center_stable = !center_change;
}
cout << "this elipse cluster is done " << center_stable << endl;
// 记录本次
center_episode[episode] = center_this;
// 更新方差
for (int cj = 0; cj < cluster_num; cj ++)
{
error_center[cj] = 0;
}
for (int pi = 0; pi < lines.size(); pi ++)
{
error_center[Point_Belong[pi]] += error_points[pi]*error_points[pi]/Center_Points_Num[Point_Belong[pi]];
}
// 选择最小方差
float error_center_sum = 0;
for (int cj = 0; cj < cluster_num; cj ++)
{
error_center_sum += error_center[cj];
cout << "error_center[" << cj << "] " << error_center[cj] << endl;
}
if (episode == 0)
{
error_center_sum_min = error_center_sum;
episode_best_index = 0;
}
else
{
if (error_center_sum_min < error_center_sum)
{
error_center_sum_min = error_center_sum;
episode_best_index = episode;
}
}
}
center_best = center_episode[episode_best_index];
if((error_center_sum_min/cluster_num) < 0.05)
{
cluster_done = true;
}
else
{
cluster_num ++;
}
cout << "error_center_sum_min = " << error_center_sum_min << ", avg = " << error_center_sum_min/cluster_num << endl;
}
lines_center->resize(center_best.size());
(*lines_center) = center_best;
for (int i = 0; i < cluster_num; i ++)
{
cout << "center_best[i] " << center_best[i] << endl;
}
}