激光条纹中心线提取算法总结和复现

滤波、分割等预处理过程省略。
输入图像为灰度图，激光条纹水平走向。

几何中心法

检测出光条边界 l、h 后，把两边界的中间线(l + h)/2作为激光条纹的中心线。

#include <iostream>
#include <opencv2/opencv.hpp>int main(int argc, char** argv)
{cv::Mat src_img = cv::imread("70.bmp", 0);cv::Mat dst_img = src_img.clone();cv::cvtColor(dst_img, dst_img, cv::COLOR_GRAY2RGB);uchar *p = src_img.data;std::vector<cv::Point> pts;for (size_t i = 0; i < src_img.cols; ++i){int y_min = INT_MIN, y_max = INT_MAX;for (size_t j = 0; j < src_img.rows; ++j){if (*(p + i + src_img.cols * j) != 0){y_min = j;break;}}for (size_t j = src_img.rows - 1; j > 0; --j){if (*(p + i + src_img.cols * j) != 0){y_max = j;break;}}if (y_min >= 0 && y_max < src_img.rows && y_min < y_max){pts.push_back(cv::Point(i, (y_min + y_max) / 2));}}for (size_t i = 0; i < pts.size(); ++i){cv::circle(dst_img, cv::Point(round(pts[i].x), round(pts[i].y)), 0.5, cv::Scalar(0, 0, 255), -1); }cv::imwrite("70.1.bmp", dst_img);system("pause");return EXIT_SUCCESS;
}

极值法

极值法是将激光条纹横截面上灰度值最大点作为激光条纹的中心。

#include <iostream>
#include <opencv2/opencv.hpp>int main(int argc, char** argv)
{cv::Mat src_img = cv::imread("70.bmp", 0);cv::Mat dst_img = src_img.clone();cv::cvtColor(dst_img, dst_img, cv::COLOR_GRAY2RGB);uchar *p = src_img.data;std::vector<cv::Point> pts;for (size_t i = 0; i < src_img.cols; ++i){int col_scalar_max = 0;int y = INT_MIN;for (size_t j = 0; j < src_img.rows; ++j){if (*(p + i + src_img.cols * j) > col_scalar_max){col_scalar_max = *(p + i + src_img.cols * j);y = j;}}pts.push_back(cv::Point(i, y));}for (size_t i = 0; i < pts.size(); ++i){cv::circle(dst_img, cv::Point(round(pts[i].x), round(pts[i].y)), 0.5, cv::Scalar(0, 0, 255), -1);}cv::imwrite("70.2.bmp", dst_img);system("pause");return EXIT_SUCCESS;
}

细化法

骨架细化法是重复地剥掉二值图像的边界像素，在剥离的过程中必须保持目标的连通性，直到得到图像的骨架。
具体原理介绍可见Zhang-Suen 图像骨架提取算法的原理和OpenCV实现

#include <iostream>
#include <opencv2/opencv.hpp>int main(int argc, char** argv)
{cv::Mat src_img = cv::imread("70.bmp", 0);cv::Mat dst_img = src_img.clone();cv::cvtColor(dst_img, dst_img, cv::COLOR_GRAY2RGB);//zhang细化算法（中轴变换法)cv::Mat copy_img = src_img.clone();while (1){bool stop = false;//step1for (int i = 1; i < src_img.cols - 1; i++)for (int j = 0; j < src_img.rows; j++){if (src_img.at<uchar>(j, i)>0){int p1 = int(src_img.at<uchar>(j, i))>0 ? 1 : 0;int p2 = int(src_img.at<uchar>(j - 1, i))>0 ? 1 : 0;int p3 = int(src_img.at<uchar>(j - 1, i + 1))>0 ? 1 : 0;int p4 = int(src_img.at<uchar>(j, i + 1))>0 ? 1 : 0;int p5 = int(src_img.at<uchar>(j + 1, i + 1))>0 ? 1 : 0;int p6 = int(src_img.at<uchar>(j + 1, i))>0 ? 1 : 0;int p7 = int(src_img.at<uchar>(j + 1, i - 1))>0 ? 1 : 0;int p8 = int(src_img.at<uchar>(j, i - 1))>0 ? 1 : 0;int p9 = int(src_img.at<uchar>(j - 1, i - 1))>0 ? 1 : 0;int np1 = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9;int sp2 = (p2 == 0 && p3 == 1) ? 1 : 0;int sp3 = (p3 == 0 && p4 == 1) ? 1 : 0;int sp4 = (p4 == 0 && p5 == 1) ? 1 : 0;int sp5 = (p5 == 0 && p6 == 1) ? 1 : 0;int sp6 = (p6 == 0 && p7 == 1) ? 1 : 0;int sp7 = (p7 == 0 && p8 == 1) ? 1 : 0;int sp8 = (p8 == 0 && p9 == 1) ? 1 : 0;int sp9 = (p9 == 0 && p2 == 1) ? 1 : 0;int sp1 = sp2 + sp3 + sp4 + sp5 + sp6 + sp7 + sp8 + sp9;if (np1 >= 2 && np1 <= 6 && sp1 == 1 && ((p2*p4*p6) == 0) && ((p4*p6*p8) == 0)){stop = true;copy_img.at<uchar>(j, i) = 0;}}}//step2for (int i = 1; i < copy_img.cols - 1; i++){for (int j = 0; j < src_img.rows; j++){if (src_img.at<uchar>(j, i)>0){int p2 = int(src_img.at<uchar>(j - 1, i))>0 ? 1 : 0;int p3 = int(src_img.at<uchar>(j - 1, i + 1)) > 0 ? 1 : 0;int p4 = int(src_img.at<uchar>(j, i + 1)) > 0 ? 1 : 0;int p5 = int(src_img.at<uchar>(j + 1, i + 1)) > 0 ? 1 : 0;int p6 = int(src_img.at<uchar>(j + 1, i)) > 0 ? 1 : 0;int p7 = int(src_img.at<uchar>(j + 1, i - 1)) > 0 ? 1 : 0;int p8 = int(src_img.at<uchar>(j, i - 1)) > 0 ? 1 : 0;int p9 = int(src_img.at<uchar>(j - 1, i - 1)) > 0 ? 1 : 0;int np1 = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9;int sp2 = (p2 == 0 && p3 == 1) ? 1 : 0;int sp3 = (p3 == 0 && p4 == 1) ? 1 : 0;int sp4 = (p4 == 0 && p5 == 1) ? 1 : 0;int sp5 = (p5 == 0 && p6 == 1) ? 1 : 0;int sp6 = (p6 == 0 && p7 == 1) ? 1 : 0;int sp7 = (p7 == 0 && p8 == 1) ? 1 : 0;int sp8 = (p8 == 0 && p9 == 1) ? 1 : 0;int sp9 = (p9 == 0 && p2 == 1) ? 1 : 0;int sp1 = sp2 + sp3 + sp4 + sp5 + sp6 + sp7 + sp8 + sp9;if (np1 >= 2 && np1 <= 6 && sp1 == 1 && (p2*p4*p8) == 0 && (p2*p6*p8) == 0){stop = true;copy_img.at<uchar>(j, i) = 0;}}}}copy_img.copyTo(src_img);if (!stop){break;}}   uchar *p = copy_img.data;std::vector<cv::Point> pts;for (int i = 0; i < copy_img.cols; i++)for (int j = 0; j < copy_img.rows; j++)if (*(p + i + src_img.cols * j) != 0)pts.push_back(cv::Point(i, j));for (size_t i = 0; i < pts.size(); ++i){cv::circle(dst_img, cv::Point(round(pts[i].x), round(pts[i].y)), 0.5, cv::Scalar(0, 0, 255), -1);}cv::imwrite("70.3.bmp", dst_img);system("pause");return EXIT_SUCCESS;
}

灰度重心法

灰度重心法就是对图像中的每一列(行)提取灰度重心作为激光条纹的中心位置。若某行的非零区间为［p，q］，则该行的灰度重心位置为:

式中，I i 是第 i 个像素点的灰度值。

#include <iostream>
#include <opencv2/opencv.hpp>int main(int argc, char** argv)
{cv::Mat src_img = cv::imread("70.bmp", 0);cv::Mat dst_img = src_img.clone();cv::cvtColor(dst_img, dst_img, cv::COLOR_GRAY2RGB);uchar *p = src_img.data;std::vector<cv::Point2f> pts;for (size_t i = 0; i < src_img.cols; ++i){int sum = 0; float y = 0; for (size_t j = 0; j < src_img.rows; ++j){int s = *(p + i + src_img.cols * j);if (s){sum += s;y += j*s;}}if (sum){y /= sum;pts.push_back(cv::Point2f(i, y));}}for (size_t i = 0; i < pts.size(); ++i){cv::circle(dst_img, cv::Point(round(pts[i].x), round(pts[i].y)), 0.5, cv::Scalar(0, 0, 255), -1);}cv::imwrite("70.4.bmp", dst_img);system("pause");return EXIT_SUCCESS;
}

法向质心法

#include <iostream>
#include <opencv2/opencv.hpp>//亚像素灰度值计算
float ijpixel(float x, float y, cv::Mat& img)
{int x_0 = int(x);int x_1 = int(x + 1);int y_0 = int(y);int y_1 = int(y + 1);int px_0y_0 = int(img.at<uchar>(y_0, x_0));int px_0y_1 = int(img.at<uchar>(y_1, x_0));int px_1y_0 = int(img.at<uchar>(y_0, x_1));int px_1y_1 = int(img.at<uchar>(y_1, x_1));float x_y0 = px_0y_0 + (x - float(x_0))*(px_1y_0 - px_0y_0);float x_y1 = px_0y_1 + (x - float(x_0))*(px_1y_1 - px_0y_1);float x_y = x_y0 + (y - float(y_0))*(x_y1 - x_y0);return x_y;
}int main(int argc, char** argv)
{cv::Mat src_img = cv::imread("70.bmp", 0);cv::Mat dst_img = src_img.clone();cv::cvtColor(dst_img, dst_img, cv::COLOR_GRAY2RGB);uchar *p = src_img.data;std::vector<cv::Point2f> pts;//灰度重心法//for (size_t i = 0; i < src_img.cols; ++i){int sum = 0;float y = 0;for (size_t j = 0; j < src_img.rows; ++j){int s = *(p + i + src_img.cols * j);if (s){sum += s;y += j*s;}}if (sum){y /= sum;pts.push_back(cv::Point2f(i, y));}}//法向迭代质心法//int iter = 0;int max_iter = 10;float distance = 0;float thre_distance = 0.05;std::vector<cv::Point2f> pts_new = pts;do{int pts_tmp_size = 5;assert((pts_tmp_size - 1) % 2 == 0);for (size_t i = (pts_tmp_size - 1) % 2; i < pts.size() - (pts_tmp_size - 1) % 2; ++i){std::vector<cv::Point2f> pts_tmp;for (size_t j = 0; j < pts_tmp_size; ++j){pts_tmp.push_back(pts[i + j - (pts_tmp_size - 1) % 2]);}cv::Vec4f line_para;cv::fitLine(pts_tmp, line_para, cv::DIST_L2, 0, 1e-2, 1e-2);    //最小二乘拟合直线方程系数float k = line_para[1] / line_para[0];float sin_theta =  1 / sqrt(k * k + 1);float cos_theta =  - k / sqrt(k * k + 1);float sum = ijpixel(pts[i].x, pts[i].y, src_img);float sumx = ijpixel(pts[i].x, pts[i].y, src_img)* pts[i].x;float sumy = ijpixel(pts[i].x, pts[i].y, src_img)* pts[i].y;int range = 10;for (size_t j = 1; j < range; ++j){float x_cor_left = pts[i].x - j*cos_theta;float y_cor_left = pts[i].y - j*sin_theta;float x_cor_right = pts[i].x + j*cos_theta;float y_cor_right = pts[i].y + j*sin_theta;if (x_cor_left >= 0 && x_cor_left<src_img.cols && y_cor_left >= 0 && y_cor_left<src_img.rows &&x_cor_right >= 0 && x_cor_right<src_img.cols && y_cor_right >= 0 && y_cor_right< src_img.rows){if (ijpixel(x_cor_left, y_cor_left, src_img) || ijpixel(x_cor_right, y_cor_right, src_img)){sum += ijpixel(x_cor_left, y_cor_left, src_img) + ijpixel(x_cor_right, y_cor_right, src_img);sumx += ijpixel(x_cor_left, y_cor_left, src_img)*x_cor_left + ijpixel(x_cor_right, y_cor_right, src_img)*x_cor_right;sumy += ijpixel(x_cor_left, y_cor_left, src_img)*y_cor_left + ijpixel(x_cor_right, y_cor_right, src_img)*y_cor_right;}}}pts_new[i].x = (float)sumx / sum;pts_new[i].y = (float)sumy / sum;}distance = 0;for (int i = 0; i < pts.size(); i++){distance += pow(pts_new[i].x - pts[i].x, 2) + pow(pts_new[i].y - pts[i].y, 2);}distance = sqrt(distance/ pts.size());iter++;std::cout << iter << "\t" << distance << std::endl;pts = pts_new;} while (iter < max_iter && distance > thre_distance);for (size_t i = 0; i < pts.size(); ++i){cv::circle(dst_img, cv::Point(round(pts[i].x), round(pts[i].y)), 0.5, cv::Scalar(0, 0, 255), -1);}cv::imwrite("70.5.bmp", dst_img);system("pause");return EXIT_SUCCESS;
}

Steger算法

#include <iostream>
#include <opencv2/opencv.hpp>int main(int argc, char** argv)
{cv::Mat src_img = cv::imread("14.bmp", 0);cv::Mat dst_img = src_img.clone();cv::cvtColor(dst_img, dst_img, cv::COLOR_GRAY2RGB);std::vector<cv::Point2f> pts;//src_img.convertTo(src_img, CV_32FC1);cv::GaussianBlur(src_img, src_img, cv::Size(0, 0), 6, 6);//一阶偏导数cv::Mat m1, m2;m1 = (cv::Mat_<float>(1, 3) << 1, 0, -1);  //x偏导m2 = (cv::Mat_<float>(3, 1) << 1, 0, -1);  //y偏导cv::Mat dx, dy;cv::filter2D(src_img, dx, CV_32FC1, m1); //卷积cv::filter2D(src_img, dy, CV_32FC1, m2);//二阶偏导数cv::Mat m3, m4, m5;m3 = (cv::Mat_<float>(1, 3) << 1, -2, 1);   //二阶x偏导m4 = (cv::Mat_<float>(3, 1) << 1, -2, 1);   //二阶y偏导m5 = (cv::Mat_<float>(2, 2) << 1, -1, -1, 1);   //二阶xy偏导cv::Mat dxx, dyy, dxy;cv::filter2D(src_img, dxx, CV_32FC1, m3);cv::filter2D(src_img, dyy, CV_32FC1, m4);cv::filter2D(src_img, dxy, CV_32FC1, m5);//hessian矩阵for (size_t i = 0; i < src_img.cols; ++i){for (size_t j = 0; j < src_img.rows; ++j){if (src_img.at<uchar>(j, i) > 50){cv::Mat hessian(2, 2, CV_32FC1);hessian.at<float>(0, 0) = dxx.at<float>(j, i);hessian.at<float>(0, 1) = dxy.at<float>(j, i);hessian.at<float>(1, 0) = dxy.at<float>(j, i);hessian.at<float>(1, 1) = dyy.at<float>(j, i);cv::Mat eValue, eVectors;cv::eigen(hessian, eValue, eVectors); //求特征值与特征向量double nx, ny;if (fabs(eValue.at<float>(0, 0)) >= fabs(eValue.at<float>(1, 0)))  //求特征值最大时对应的特征向量{nx = eVectors.at<float>(0, 0);ny = eVectors.at<float>(0, 1);}else{nx = eVectors.at<float>(1, 0);ny = eVectors.at<float>(1, 1);}double t = -(nx*dx.at<float>(j, i) + ny*dy.at<float>(j, i)) /(nx*nx*dxx.at<float>(j, i) + 2 * nx*ny*dxy.at<float>(j, i) + ny*ny*dyy.at<float>(j, i));if (fabs(t*nx) <= 0.5 && fabs(t*ny) <= 0.5){pts.push_back(cv::Point2f(i + t*nx, j + t*ny));}}}}for (size_t i = 0; i < pts.size(); ++i){cv::circle(dst_img, cv::Point(round(pts[i].x), round(pts[i].y)), 0.5, cv::Scalar(0, 0, 255), -1);}cv::imwrite("70.6.bmp", dst_img);system("pause");return EXIT_SUCCESS;
}

后续文章：中心线提取–GPU加速

激光条纹中心线提取算法总结和复现相关推荐

【必备知识】：线激光条纹中心线提取算法导读
线激光条纹特性线激光器是由点激光器和前置透镜组成的.点激光器可以为He-Ne激光器或半导体激光器.相比较He-Ne激光器,半导体激光器因其输出光源具有发散性,更适合用于制作线激光器.需要说明的是,半 ...
传统激光条纹中心提取算法研究现状
传统激光条纹中心提取算法研究现状前言一.边缘法二.中心法三.阈值法四. 细化法五.极值法六.灰度重心法七.方向模板八.曲线拟合法九.Steger 前言光条中心提取是将宽度大于1的 ...
激光条纹中心提取——方法总结
激光条纹中心提取--方法总结算法优势缺点边缘法处理速度快:适用于精度要求低的大型物体测量存在很大误差:要求图像质量较好且结构光特性较高中心法适用于条纹质量好且形状规则的物体测量:精度高 ...
激光条纹中心提取——灰度中心法python
激光条纹中心提取--灰度中心法python 灰度中心法 python代码灰度中心法灰度重心法是根据每行光条纹横截面内的灰度分布特征逐行进行处理,通过在行坐标的方向上,逐行计算提取光条纹区域的灰度重 ...
激光条纹中心提取——Zhang-Suen法python
Zhang-Suen法原理-- Zhang-Suen法代码--python代码原理-- Zhang-Suen法细化法(又称形态学骨架法)是通过对光条纹不断地进行腐蚀操作,剥离光条纹边界,得到单 ...
代码：法向质心法提取线激光条纹中心线（CPP+OpenCV）
激光条纹中心提取——ZhangSuen法python
ZhangSuen法: 论文连接:A fast parallel algorithm for thinning digital patterns 代码连接:https://github.com/bsd ...
Steger算法（Line_Gauss）-光条中心线提取（基于Hessian矩阵）
Steger算法(Line_Gauss)-光条中心线提取(基于Hessian矩阵) 算法背景介绍 Hessian 矩阵与泰勒多项式关于求t 导数与中心点.亚像素点高斯函数作用文献算法背景介绍 ...
中线提取算法_综述|线结构光中心提取算法研究发展
摘要: 线结构光扫描是三维重建领域的关键技术.光条纹中心提取算法是决定线结构光三维重建精度以及光条纹轮廓定位准确性的重要因素.本文详细阐述了光条纹中心提取算法的理论基础及发展历程,将现有算法分为三类 ...
基于灰度质心法和骨架的激光中心线提取
之前博主一直在做线结构光成像,硬件比较垃圾,相机加镜头和线激光器总共成本在1000以内,精度在0.1mm左右,感觉这种成本做出来还是不错的,其实主要大部分时间花在了分析上来达到最好的效果. 一般对于激 ...