概况由于最近看到一篇博客Steger算法实现结构光光条中心提取(python版本)-CSDN博客其中描述了关于Steger算法的一些问题刚好我也遇到相同的问题于是我参考了一篇名为《基于多尺度分析的激光光条中心点坐标提取方法》[1]光学学报的文章进行改进。该文章是北航在14年提出的该算法主要是通过优化方差的方式来优化Steger这个方法大量被学者借鉴如《Laser stripe center extraction method base on Hessian matrix improved by stripe width precise calculation》[2]该文章通过BP神经网络来预测方差从而优化Steger算法。基于文献[1]改进的代码为import cv2 import numpy as np import time from scipy.spatial import cKDTree # # 工具函数生成高斯导数核 # def get_gaussian_deriv_kernels(sigma, radius): size 2 * radius 1 x np.arange(-radius, radius 1, dtypenp.float32) y np.arange(-radius, radius 1, dtypenp.float32) xx, yy np.meshgrid(x, y) sigma2 sigma ** 2 sigma4 sigma2 ** 2 g np.exp(-(xx**2 yy**2) / (2 * sigma2)) / (2 * np.pi * sigma2) gx -xx / sigma2 * g gy -yy / sigma2 * g gxx (xx**2 - sigma2) / sigma4 * g gyy (yy**2 - sigma2) / sigma4 * g gxy xx * yy / sigma4 * g return gx, gy, gxx, gyy, gxy # # 单尺度整图导数计算加速核心 # def compute_scale_maps(gray_img, sigma): radius int(np.ceil(3 * sigma)) gx, gy, gxx, gyy, gxy get_gaussian_deriv_kernels(sigma, radius) ru cv2.filter2D(gray_img, cv2.CV_32F, gx, borderTypecv2.BORDER_REPLICATE) rv cv2.filter2D(gray_img, cv2.CV_32F, gy, borderTypecv2.BORDER_REPLICATE) ruu cv2.filter2D(gray_img, cv2.CV_32F, gxx, borderTypecv2.BORDER_REPLICATE) rvv cv2.filter2D(gray_img, cv2.CV_32F, gyy, borderTypecv2.BORDER_REPLICATE) ruv cv2.filter2D(gray_img, cv2.CV_32F, gxy, borderTypecv2.BORDER_REPLICATE) # 向量化计算特征值与C响应 trace ruu rvv delta (ruu - rvv)**2 4 * ruv**2 sqrt_delta np.sqrt(np.maximum(delta, 0)) lambda1 (trace sqrt_delta) / 2.0 lambda2 (trace - sqrt_delta) / 2.0 abs_l1, abs_l2 np.abs(lambda1), np.abs(lambda2) mask abs_l1 abs_l2 lambda_max np.where(mask, lambda1, lambda2) abs_lambda_max np.where(mask, abs_l1, abs_l2) C_val sigma * sigma * abs_lambda_max # 法线方向 nx np.zeros_like(ruu) ny np.ones_like(ruu) mask_ruv np.abs(ruv) 1e-6 nx[mask_ruv] 1.0 ny[mask_ruv] (lambda_max[mask_ruv] - ruu[mask_ruv]) / ruv[mask_ruv] norm_len np.sqrt(nx**2 ny**2) norm_len[norm_len 1e-6] 1.0 nx / norm_len ny / norm_len return { sigma: sigma, C: C_val, nx: nx, ny: ny, ru: ru, rv: rv, ruu: ruu, rvv: rvv, ruv: ruv } # # 自适应多尺度提取类 # class AdaptiveMultiScaleExtractor: def __init__(self, threshold_C15.0): self.threshold_C threshold_C self.k_search 1 # RON搜索半径 self.k_neighbor 2 # 拟合直线邻近点数 self.profile_range 20 # 法线方向剖面采样半长 def _get_skeleton(self, gray_img): 二值化最大连通域筛选骨架提取 _, binary cv2.threshold(gray_img, 0, 255, cv2.THRESH_BINARY cv2.THRESH_OTSU) kernel cv2.getStructuringElement(cv2.MORPH_RECT, (3, 3)) binary cv2.morphologyEx(binary, cv2.MORPH_OPEN, kernel) # 新增只保留最大连通域过滤杂散伪光条 num_labels, labels, stats, _ cv2.connectedComponentsWithStats(binary, connectivity8) if num_labels 1: max_idx np.argmax(stats[1:, cv2.CC_STAT_AREA]) 1 binary (labels max_idx).astype(np.uint8) * 255 # 骨架提取 try: from cv2 import ximgproc skeleton ximgproc.thinning(binary, thinningTypecv2.ximgproc.THINNING_ZHANGSUEN) except ImportError: dist cv2.distanceTransform(binary, cv2.DIST_L2, 5) skeleton (dist 2).astype(np.uint8) * 255 skeleton_pts np.column_stack(np.where(skeleton 0)) return skeleton_pts, binary def _estimate_width_per_point(self, gray_img, skeleton_pts, kdtree, idx): 逐点估计光条宽度2w 法线方向对应论文2.1节 h, w gray_img.shape py, px skeleton_pts[idx] # k近邻拟合直线求法线 _, indices kdtree.query([py, px], k2 * self.k_neighbor 1) near_pts skeleton_pts[indices[1:]] if len(near_pts) 2: return None, 0.0 near_xy near_pts[:, ::-1].astype(np.float32) [vx, vy, _, _] cv2.fitLine(near_xy, cv2.DIST_L2, 0, 0.01, 0.01) normal np.array([-vy[0], vx[0]]) # 切线→法线 normal normal / (np.linalg.norm(normal) 1e-8) # 沿法线采样灰度剖面0.2A法求宽度 profile [] for d in range(-self.profile_range, self.profile_range 1): sx int(round(px d * normal[0])) sy int(round(py d * normal[1])) if 0 sy h and 0 sx w: profile.append(float(gray_img[sy, sx])) else: profile.append(0.0) peak max(profile) if peak 20: return normal, 0.0 thr 0.2 * peak mid len(profile) // 2 left mid while left 0 and profile[left] thr: left - 1 right mid while right len(profile) - 1 and profile[right] thr: right 1 width_2w max(float(right - left), 1.0) return normal, width_2w def extract(self, gray_img): h, w gray_img.shape final_points [] # 1. 提取骨架初始点 skeleton_pts, _ self._get_skeleton(gray_img) if len(skeleton_pts) 0: return np.array([]) kdtree cKDTree(skeleton_pts) # 2. 逐点估算宽度统计全局宽度范围生成分档尺度 width_list [] normal_list [] for i in range(len(skeleton_pts)): norm, w2 self._estimate_width_per_point(gray_img, skeleton_pts, kdtree, i) normal_list.append(norm) if w2 1.0: width_list.append(w2) if not width_list: return np.array([]) # 生成5档覆盖全范围的σ从最细到最粗 min_w max(np.percentile(width_list, 10), 2.0) max_w np.percentile(width_list, 90) * 1.2 # 等比生成5个宽度档位 width_levels np.geomspace(min_w, max_w, 5) sigma_levels (width_levels / 2.0) / np.sqrt(3.0) # 3. 预计算所有档位的整图导数图一次卷积全程复用 scale_maps [] for sigma in sigma_levels: scale_maps.append(compute_scale_maps(gray_img, sigma)) # 4. 逐骨架点搜索最优尺度与亚像素中心 for idx in range(len(skeleton_pts)): py, px skeleton_pts[idx] w2 width_list[idx] if idx len(width_list) else min_w if w2 1.0: continue # 找到当前宽度最匹配的3个连续尺度档位 target_sigma (w2 / 2.0) / np.sqrt(3.0) scale_dists [abs(s - target_sigma) for s in sigma_levels] best_scale_idx np.argmin(scale_dists) # 取前后共3个尺度做精细比较 select_indices [] for offset in [-1, 0, 1]: s_idx best_scale_idx offset if 0 s_idx len(sigma_levels): select_indices.append(s_idx) # 不足3个则补全 while len(select_indices) 3: if select_indices[0] 0: select_indices.insert(0, select_indices[0] - 1) elif select_indices[-1] len(sigma_levels) - 1: select_indices.append(select_indices[-1] 1) else: break best_C -1.0 best_sdata None best_yx (py, px) # RON邻域搜索 for dy in range(-self.k_search, self.k_search 1): for dx in range(-self.k_search, self.k_search 1): cy, cx py dy, px dx if not (0 cy h and 0 cx w): continue for s_idx in select_indices: sdata scale_maps[s_idx] c_val sdata[C][cy, cx] if c_val best_C: best_C c_val best_sdata sdata best_yx (cy, cx) if best_C self.threshold_C or best_sdata is None: continue # 亚像素求解 cy, cx best_yx nx best_sdata[nx][cy, cx] ny best_sdata[ny][cy, cx] ru best_sdata[ru][cy, cx] rv best_sdata[rv][cy, cx] ruu best_sdata[ruu][cy, cx] rvv best_sdata[rvv][cy, cx] ruv best_sdata[ruv][cy, cx] denom nx*nx*ruu 2*nx*ny*ruv ny*ny*rvv if abs(denom) 1e-6: continue t -(nx*ru ny*rv) / denom if abs(t) 0.6: continue sub_x cx t * nx sub_y cy t * ny final_points.append([sub_x, sub_y]) return np.array(final_points) # # 主程序 # if __name__ __main__: img cv2.imread(test_3.bmp, cv2.IMREAD_GRAYSCALE) if img is None: print(图片读取失败请检查路径) exit() extractor AdaptiveMultiScaleExtractor(threshold_C25.0) start time.time() centers extractor.extract(img) cost time.time() - start print(f提取中心点数量{len(centers)}) print(f耗时{cost*1000:.2f} ms) # 可视化 result cv2.cvtColor(img, cv2.COLOR_GRAY2BGR) for pt in centers: cv2.circle(result, (int(round(pt[0])), int(round(pt[1]))), 1, (0, 0, 255), -1) cv2.imshow(Adaptive Multi-scale Result, result) cv2.imwrite(Result_fixed.png, result) cv2.waitKey(0) cv2.destroyAllWindows()提取效果博客Steger算法实现结构光光条中心提取(python版本)-CSDN博客的代码运行结果注两个图像我没有上传原图我只截一部分总结参考文献[1]描述该论文的方法能适用于光条宽度发生较大变化且较为反光的场景。但是由于本人实验条件有限未做相应实验进行验证。[1]李凤娇,李小菁,刘震.基于多尺度分析的激光光条中心点坐标提取方法[J].光学学报,2014,34(11):111-116.[2] Bo Q , Hou B , Miao X W Y .Laser stripe center extraction method base on Hessian matrix improved by stripe width precise calculation[J].Optics and Lasers in Engineering, 2024, 172(Jan.):1.1-1.9.DOI:10.1016/j.optlaseng.2023.107896.
