• Steger算法实现结构光光条中心提取(python版本)


    Steger算法原理

    对结构光进行光条中心提取时,Steger算法是以Hessian矩阵为基础的。它的基础步骤如下所示:

    1. 从Hessian矩阵中求出线激光条纹的法线方向
    2. 在光条纹法线方向上将其灰度分布按照泰勒多项式展开,求取的极大值即为光条在该法线方向上的亚像素坐标。对于二维离散图像I(u,v)来说,Hessian矩阵可以表示为:

    H(u,v) = \begin{bmatrix} I_{uu} &I_{uv} \\ I_{uv} & I_{vv} \end{bmatrix}

    这里面的u,v表示的就是像素的行坐标和列坐标,I_{uv}代表像素(u,v)灰度也可以称之为灰度分布函数。而I_{uu}I_{uv}I_{vv}都可以通过I_{uv}与二维高斯函数G(u,v)卷积运算得到。

    G(u,v) = \frac{1}{\sigma\sqrt{2 \pi}}\cdot e^{\left ( -\frac{u^{2}+v^{2}}{2\sigma^{2}} \right )}

    I_{uu}=\frac{\partial ^{2}}{\partial u\partial u}G(u,v)\bigotimes I_{uv}

    I_{uv}=\frac{\partial ^{2}}{\partial u\partial v}G(u,v)\bigotimes I_{uv}

    I_{vv}=\frac{\partial ^{2}}{\partial v\partial v}G(u,v)\bigotimes I_{uv} 

    在这里,二维高斯函数,其主要作用是为了让光条灰度分布特性更加明显,在G(u,v)表达式中,\sigma是标准差,一般取\sigma \geq \frac{W}{\sqrt{3}},W代表光条宽度。在像素(u,v)处的Hessian矩阵有两个特征向量,其中一个绝对值较大,为该像素处的法线方向向量,而另外一个则是切向方向向量。因此可以通过求取Hessian矩阵的特征向量来计算出法线方向。某一像素点H(u_{0},v_{0}),在二阶泰勒展开式Hessian矩阵为:

    H(u0,v0)=[IuuIuvIuvIvv]

    由该点的Hessian矩阵所求出的特征值和特征向量分别与该点的法线方向和该方向的二阶方向导数相对应。其中法线方向的单位向量为:e=[eu,ev],并且在光条纹法线方向上的像素点(u0+teu,v0+tev),并且在光条纹法线方向上的像素点I=(u0+teu,v0+tev)可以由像素(u0,v0)的灰度I(u0,v0)和二阶泰勒展开多项式表示为:

    I(u0+teu,v0+tev)=I(u0,v0)+te[Iu,Iv]T+teH(u,v)eT

    t=euIu+evIve2uIuu+2euevIuv+e2vIvv

    再将t(即泰勒展开式)代入其中即可求出光条纹中心的亚像素坐标。

    Steger算法Python实现

    这一部分,我在网上找到了许多C++版本的Steger算法实现,本人参考了现有的C++版本Steger算法实现,在此基础上进行了改进,用Python重新去实现该算法。

    计算图像的一阶导数和二阶导数

    这里我采用了三种方法,分别是自定义卷积、Scharr 滤波器、Sobel 滤波器来计算图像的导数和二阶导数。

    1. import cv2
    2. import numpy as np
    3. def _derivation_with_Filter(Gaussimg):
    4. dx = cv2.filter2D(Gaussimg,-1, kernel=np.array([[1], [0], [-1]]))
    5. dy = cv2.filter2D(Gaussimg,-1, kernel=np.array([[1, 0, -1]]))
    6. dxx = cv2.filter2D(Gaussimg,-1, kernel=np.array([[1], [-2], [1]]))
    7. dyy = cv2.filter2D(Gaussimg,-1, kernel=np.array([[1, -2, 1]]))
    8. dxy = cv2.filter2D(Gaussimg,-1, kernel=np.array([[1, -1], [-1, 1]]))
    9. return dx, dy, dxx, dyy, dxy
    10. def _derivation_with_Scharr(Gaussimg):
    11. dx = cv2.Scharr(Gaussimg, cv2.CV_32F, 1, 0)
    12. dy = cv2.Scharr(Gaussimg, cv2.CV_32F, 0, 1)
    13. dxx = cv2.Scharr(dx, cv2.CV_32F, 1, 0)
    14. dxy = cv2.Scharr(dx, cv2.CV_32F, 0, 1)
    15. dyy = cv2.Scharr(dy, cv2.CV_32F, 0, 1)
    16. return dx, dy, dxx, dyy, dxy
    17. def _derivation_with_Sobel(Gaussimg):
    18. dx = cv2.Sobel(Gaussimg, cv2.CV_32F, 1, 0, ksize=3)
    19. dy = cv2.Sobel(Gaussimg, cv2.CV_32F, 0, 1, ksize=3)
    20. dxx = cv2.Sobel(dx, cv2.CV_32F, 1, 0, ksize=3)
    21. dxy = cv2.Sobel(dx, cv2.CV_32F, 0, 1, ksize=3)
    22. dyy = cv2.Sobel(dy, cv2.CV_32F, 0, 1, ksize=3)
    23. return dx, dy, dxx, dyy, dxy
    24. if __name__=="__main__":
    25. from pyzjr.dlearn import Runcodes
    26. sigmaX, sigmaY = 1.1, 1.1
    27. img = cv2.imread(r"D:\PythonProject\net\line\linedata\Image_1.jpg")
    28. gray_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    29. Gaussimg = cv2.GaussianBlur(gray_img, ksize=(0, 0), sigmaX=sigmaX, sigmaY=sigmaY)
    30. with Runcodes("Filter"):
    31. dx, dy, dxx, dyy, dxy = _derivation_with_Filter(Gaussimg)
    32. print(type(dx[0][0]))
    33. with Runcodes("Scharr"):
    34. dx1, dy1, dxx1, dyy1, dxy1 = _derivation_with_Scharr(Gaussimg)
    35. print(type(dx1[0][0]))
    36. with Runcodes("Sobel"):
    37. dx2, dy2, dxx2, dyy2, dxy2 = _derivation_with_Sobel(Gaussimg)
    38. print(type(dx2[0][0]))

    在控制台中输出为:


    Filter: 0.00602 sec

    Scharr: 0.00770 sec

    Sobel: 0.01025 sec

    从运算时间上来看,自定义的Filter速度最快,scharr滤波器居中,sobel滤波器最慢;从元素类型上来看,第一种为uint8,另外两种为float32类型。

    通过选择计算方法可以去创建Hessian矩阵。

    1. def Magnitudefield(dxx, dyy):
    2. """计算幅度和相位"""
    3. dxx = dxx.astype(float)
    4. dyy = dyy.astype(float)
    5. mag = cv2.magnitude(dxx, dyy)
    6. phase = mag*180./np.pi
    7. return mag, phase
    8. def derivation(gray_img, sigmaX, sigmaY, method="Scharr", nonmax=False):
    9. """
    10. 计算图像的一阶导数 dx 和 dy,以及二阶导数 dxx、dyy 和 dxy
    11. :param gray_img: 灰度图
    12. :param sigmaX: 在水平方向的高斯核标准差,用于激光线提取建议取1-2
    13. :param sigmaY: 在垂直方向上的高斯核标准差,用于激光线提取建议取1-2
    14. :param method:"Scharr" or "Filter" or "Sobel"
    15. 选择什么方式获取dx, dy, dxx, dyy, dxy,提供了卷积与Scharr和Sobel滤波器三种方式计算,
    16. Scharr滤波器通常会产生更平滑和准确的结果,所以这里默认使用"Scharr"方法,虽然
    17. "Sobel"运行比另外两种要慢,但在使用的时候,建议自己试试
    18. :return: dx, dy, dxx, dyy, dxy
    19. """
    20. Gaussimg = cv2.GaussianBlur(gray_img, ksize=(0, 0), sigmaX=sigmaX, sigmaY=sigmaY)
    21. dx, dy, dxx, dyy, dxy = _derivation_with_Scharr(Gaussimg)
    22. if method == "Filter":
    23. dx, dy, dxx, dyy, dxy = _derivation_with_Filter(Gaussimg)
    24. elif method == "Sobel":
    25. dx, dy, dxx, dyy, dxy =_derivation_with_Sobel(Gaussimg)
    26. if nonmax:
    27. normal, phase = Magnitudefield(dxx, dyy)
    28. dxy = nonMaxSuppression(normal, phase)
    29. return dx, dy, dxx, dyy, dxy
    30. def nonMaxSuppression(det, phase):
    31. """非最大值抑制"""
    32. gmax = np.zeros(det.shape)
    33. # thin-out evry edge for angle = [0, 45, 90, 135]
    34. for i in range(gmax.shape[0]):
    35. for j in range(gmax.shape[1]):
    36. if phase[i][j] < 0:
    37. phase[i][j] += 360
    38. if ((j+1) < gmax.shape[1]) and ((j-1) >= 0) and ((i+1) < gmax.shape[0]) and ((i-1) >= 0):
    39. # 0 degrees
    40. if (phase[i][j] >= 337.5 or phase[i][j] < 22.5) or (phase[i][j] >= 157.5 and phase[i][j] < 202.5):
    41. if det[i][j] >= det[i][j + 1] and det[i][j] >= det[i][j - 1]:
    42. gmax[i][j] = det[i][j]
    43. # 45 degrees
    44. if (phase[i][j] >= 22.5 and phase[i][j] < 67.5) or (phase[i][j] >= 202.5 and phase[i][j] < 247.5):
    45. if det[i][j] >= det[i - 1][j + 1] and det[i][j] >= det[i + 1][j - 1]:
    46. gmax[i][j] = det[i][j]
    47. # 90 degrees
    48. if (phase[i][j] >= 67.5 and phase[i][j] < 112.5) or (phase[i][j] >= 247.5 and phase[i][j] < 292.5):
    49. if det[i][j] >= det[i - 1][j] and det[i][j] >= det[i + 1][j]:
    50. gmax[i][j] = det[i][j]
    51. # 135 degrees
    52. if (phase[i][j] >= 112.5 and phase[i][j] < 157.5) or (phase[i][j] >= 292.5 and phase[i][j] < 337.5):
    53. if det[i][j] >= det[i - 1][j - 1] and det[i][j] >= det[i + 1][j + 1]:
    54. gmax[i][j] = det[i][j]
    55. return gmax

    Magnitudefield(dxx, dyy)函数计算梯度的幅度和相位,使用给定的 x 和 y 导数 (dxx 和 dyy)。幅度代表梯度的强度,相位代表梯度的方向;nonMaxSuppression(det, phase)函数根据梯度的相位 (phase) 对梯度幅度 (det) 执行非最大值抑制。它有助于在梯度方向上仅保留局部最大值,从而细化检测到的边缘;然后应用于derivation函数中去调控。

    二阶泰勒展开式Hessian矩阵

    1. def HessianMatrix(self, dx, dy, dxx, dyy, dxy, threshold=0.5):
    2. """
    3. HessianMatrix = [dxx dxy]
    4. [dxy dyy]
    5. compute hessian:
    6. [dxx dxy] [00 01]
    7. ====>
    8. [dxy dyy] [10 11]
    9. """
    10. point=[]
    11. direction=[]
    12. value=[]
    13. for x in range(0, self.col):
    14. for y in range(0, self.row):
    15. if dxy[y,x] > 0:
    16. hessian = np.zeros((2,2))
    17. hessian[0,0] = dxx[y,x]
    18. hessian[0,1] = dxy[y,x]
    19. hessian[1,0] = dxy[y,x]
    20. hessian[1,1] = dyy[y,x]
    21. # 计算矩阵的特征值和特征向量
    22. _, eigenvalues, eigenvectors = cv2.eigen(hessian)
    23. if np.abs(eigenvalues[0,0]) >= np.abs(eigenvalues[1,0]):
    24. nx = eigenvectors[0,0]
    25. ny = eigenvectors[0,1]
    26. else:
    27. nx = eigenvectors[1,0]
    28. ny = eigenvectors[1,1]
    29. # Taylor展开式分子分母部分,需要避免为0的情况
    30. Taylor_numer = (dx[y, x] * nx + dy[y, x] * ny)
    31. Taylor_denom = dxx[y,x]*nx*nx + dyy[y,x]*ny*ny + 2*dxy[y,x]*nx*ny
    32. if Taylor_denom != 0:
    33. T = -(Taylor_numer/Taylor_denom)
    34. # Hessian矩阵最大特征值对应的特征向量对应于光条的法线方向
    35. if np.abs(T*nx) <= threshold and np.abs(T*ny) <= threshold:
    36. point.append((x,y))
    37. direction.append((nx,ny))
    38. value.append(np.abs(dxy[y,x]+dxy[y,x]))
    39. return point, direction, value

    这一部分我定义在steger类下了

    绘制中心线在背景与原图中

    1. def centerline(self,img, sigmaX, sigmaY, method="Scharr", usenonmax=True, color=(0, 0, 255)):
    2. dx, dy, dxx, dyy, dxy = derivation(self.gray_image, sigmaX, sigmaY, method=method, nonmax=usenonmax)
    3. point, direction, value = self.HessianMatrix(dx, dy, dxx, dyy, dxy)
    4. for point in point:
    5. self.newimage[point[1], point[0]] = 255
    6. img[point[1], point[0], :] = color
    7. return img, self.newimage

    代码以及效果展示

    1. import cv2
    2. import numpy as np
    3. def _derivation_with_Filter(Gaussimg):
    4. dx = cv2.filter2D(Gaussimg,-1, kernel=np.array([[1], [0], [-1]]))
    5. dy = cv2.filter2D(Gaussimg,-1, kernel=np.array([[1, 0, -1]]))
    6. dxx = cv2.filter2D(Gaussimg,-1, kernel=np.array([[1], [-2], [1]]))
    7. dyy = cv2.filter2D(Gaussimg,-1, kernel=np.array([[1, -2, 1]]))
    8. dxy = cv2.filter2D(Gaussimg,-1, kernel=np.array([[1, -1], [-1, 1]]))
    9. return dx, dy, dxx, dyy, dxy
    10. def _derivation_with_Scharr(Gaussimg):
    11. dx = cv2.Scharr(Gaussimg, cv2.CV_32F, 1, 0)
    12. dy = cv2.Scharr(Gaussimg, cv2.CV_32F, 0, 1)
    13. dxx = cv2.Scharr(dx, cv2.CV_32F, 1, 0)
    14. dxy = cv2.Scharr(dx, cv2.CV_32F, 0, 1)
    15. dyy = cv2.Scharr(dy, cv2.CV_32F, 0, 1)
    16. return dx, dy, dxx, dyy, dxy
    17. def _derivation_with_Sobel(Gaussimg):
    18. dx = cv2.Sobel(Gaussimg, cv2.CV_32F, 1, 0, ksize=3)
    19. dy = cv2.Sobel(Gaussimg, cv2.CV_32F, 0, 1, ksize=3)
    20. dxx = cv2.Sobel(dx, cv2.CV_32F, 1, 0, ksize=3)
    21. dxy = cv2.Sobel(dx, cv2.CV_32F, 0, 1, ksize=3)
    22. dyy = cv2.Sobel(dy, cv2.CV_32F, 0, 1, ksize=3)
    23. return dx, dy, dxx, dyy, dxy
    24. def Magnitudefield(dxx, dyy):
    25. """计算幅度和相位"""
    26. dxx = dxx.astype(float)
    27. dyy = dyy.astype(float)
    28. mag = cv2.magnitude(dxx, dyy)
    29. phase = mag*180./np.pi
    30. return mag, phase
    31. def derivation(gray_img, sigmaX, sigmaY, method="Scharr", nonmax=False):
    32. """
    33. 计算图像的一阶导数 dx 和 dy,以及二阶导数 dxx、dyy 和 dxy
    34. :param gray_img: 灰度图
    35. :param sigmaX: 在水平方向的高斯核标准差,用于激光线提取建议取1-2
    36. :param sigmaY: 在垂直方向上的高斯核标准差,用于激光线提取建议取1-2
    37. :param method:"Scharr" or "Filter" or "Sobel"
    38. 选择什么方式获取dx, dy, dxx, dyy, dxy,提供了卷积与Scharr和Sobel滤波器三种方式计算,
    39. Scharr滤波器通常会产生更平滑和准确的结果,所以这里默认使用"Scharr"方法,虽然
    40. "Sobel"运行比另外两种要慢,但在使用的时候,建议自己试试
    41. :return: dx, dy, dxx, dyy, dxy
    42. """
    43. Gaussimg = cv2.GaussianBlur(gray_img, ksize=(0, 0), sigmaX=sigmaX, sigmaY=sigmaY)
    44. dx, dy, dxx, dyy, dxy = _derivation_with_Scharr(Gaussimg)
    45. if method == "Filter":
    46. dx, dy, dxx, dyy, dxy = _derivation_with_Filter(Gaussimg)
    47. elif method == "Sobel":
    48. dx, dy, dxx, dyy, dxy =_derivation_with_Sobel(Gaussimg)
    49. if nonmax:
    50. normal, phase = Magnitudefield(dxx, dyy)
    51. dxy = nonMaxSuppression(normal, phase)
    52. return dx, dy, dxx, dyy, dxy
    53. def nonMaxSuppression(det, phase):
    54. """非最大值抑制"""
    55. gmax = np.zeros(det.shape)
    56. # thin-out evry edge for angle = [0, 45, 90, 135]
    57. for i in range(gmax.shape[0]):
    58. for j in range(gmax.shape[1]):
    59. if phase[i][j] < 0:
    60. phase[i][j] += 360
    61. if ((j+1) < gmax.shape[1]) and ((j-1) >= 0) and ((i+1) < gmax.shape[0]) and ((i-1) >= 0):
    62. # 0 degrees
    63. if (phase[i][j] >= 337.5 or phase[i][j] < 22.5) or (phase[i][j] >= 157.5 and phase[i][j] < 202.5):
    64. if det[i][j] >= det[i][j + 1] and det[i][j] >= det[i][j - 1]:
    65. gmax[i][j] = det[i][j]
    66. # 45 degrees
    67. if (phase[i][j] >= 22.5 and phase[i][j] < 67.5) or (phase[i][j] >= 202.5 and phase[i][j] < 247.5):
    68. if det[i][j] >= det[i - 1][j + 1] and det[i][j] >= det[i + 1][j - 1]:
    69. gmax[i][j] = det[i][j]
    70. # 90 degrees
    71. if (phase[i][j] >= 67.5 and phase[i][j] < 112.5) or (phase[i][j] >= 247.5 and phase[i][j] < 292.5):
    72. if det[i][j] >= det[i - 1][j] and det[i][j] >= det[i + 1][j]:
    73. gmax[i][j] = det[i][j]
    74. # 135 degrees
    75. if (phase[i][j] >= 112.5 and phase[i][j] < 157.5) or (phase[i][j] >= 292.5 and phase[i][j] < 337.5):
    76. if det[i][j] >= det[i - 1][j - 1] and det[i][j] >= det[i + 1][j + 1]:
    77. gmax[i][j] = det[i][j]
    78. return gmax
    79. class Steger():
    80. def __init__(self, gray_image):
    81. self.gray_image = gray_image
    82. self.row, self.col = self.gray_image.shape[:2]
    83. self.newimage = np.zeros((self.row, self.col), np.uint8)
    84. def centerline(self,img, sigmaX, sigmaY, method="Scharr", usenonmax=True, color=(0, 0, 255)):
    85. dx, dy, dxx, dyy, dxy = derivation(self.gray_image, sigmaX, sigmaY, method=method, nonmax=usenonmax)
    86. point, direction, value = self.HessianMatrix(dx, dy, dxx, dyy, dxy)
    87. for point in point:
    88. self.newimage[point[1], point[0]] = 255
    89. img[point[1], point[0], :] = color
    90. return img, self.newimage
    91. def HessianMatrix(self, dx, dy, dxx, dyy, dxy, threshold=0.5):
    92. """
    93. HessianMatrix = [dxx dxy]
    94. [dxy dyy]
    95. compute hessian:
    96. [dxx dxy] [00 01]
    97. ====>
    98. [dxy dyy] [10 11]
    99. """
    100. point=[]
    101. direction=[]
    102. value=[]
    103. for x in range(0, self.col):
    104. for y in range(0, self.row):
    105. if dxy[y,x] > 0:
    106. hessian = np.zeros((2,2))
    107. hessian[0,0] = dxx[y,x]
    108. hessian[0,1] = dxy[y,x]
    109. hessian[1,0] = dxy[y,x]
    110. hessian[1,1] = dyy[y,x]
    111. # 计算矩阵的特征值和特征向量
    112. _, eigenvalues, eigenvectors = cv2.eigen(hessian)
    113. if np.abs(eigenvalues[0,0]) >= np.abs(eigenvalues[1,0]):
    114. nx = eigenvectors[0,0]
    115. ny = eigenvectors[0,1]
    116. else:
    117. nx = eigenvectors[1,0]
    118. ny = eigenvectors[1,1]
    119. # Taylor展开式分子分母部分,需要避免为0的情况
    120. Taylor_numer = (dx[y, x] * nx + dy[y, x] * ny)
    121. Taylor_denom = dxx[y,x]*nx*nx + dyy[y,x]*ny*ny + 2*dxy[y,x]*nx*ny
    122. if Taylor_denom != 0:
    123. T = -(Taylor_numer/Taylor_denom)
    124. # Hessian矩阵最大特征值对应的特征向量对应于光条的法线方向
    125. if np.abs(T*nx) <= threshold and np.abs(T*ny) <= threshold:
    126. point.append((x,y))
    127. direction.append((nx,ny))
    128. value.append(np.abs(dxy[y,x]+dxy[y,x]))
    129. return point, direction, value
    130. def threshold(mask, std=127.5):
    131. """阈值法"""
    132. mask[mask > std] = 255
    133. mask[mask < std] = 0
    134. return mask.astype("uint8")
    135. def Bessel(xi: list):
    136. """贝塞尔公式"""
    137. xi_array = np.array(xi)
    138. x_average = np.mean(xi_array)
    139. squared_diff = (xi_array - x_average) ** 2
    140. variance = squared_diff / (len(xi)-1)
    141. bessel = np.sqrt(variance)
    142. return bessel
    143. def test_Steger():
    144. img = cv2.imread(r"D:\PythonProject\net\line\data\Image_20230215160728679.jpg")
    145. gray_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    146. gray_img = threshold(gray_img)
    147. steger = Steger(gray_img)
    148. img,newimage=steger.centerline(img,sigmaX=1.1,sigmaY=1.1,method="Sobel")
    149. # point, direction, value = HessianMatrix(gray_img, 1.1, 1.1, usenonmax=True)
    150. cv2.imwrite("result/main3.png", newimage)
    151. cv2.imwrite("result/main3_img.png", img)
    152. def test_derivation_methods_time():
    153. from pyzjr.dlearn import Runcodes
    154. img = cv2.imread(r"D:\PythonProject\net\line\data\Image_20230215160728411.jpg")
    155. img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
    156. gray_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    157. with Runcodes(description="Filter"):
    158. print(derivation(gray_img,1.1,1.1,"Filter"))
    159. # Filter: 0.01344 sec
    160. with Runcodes(description="Scharr"):
    161. print(derivation(gray_img, 1.1, 1.1))
    162. # Scharr: 0.00959 sec
    163. with Runcodes(description="Sobel"):
    164. print(derivation(gray_img, 1.1, 1.1,"Sobel"))
    165. # Sobel: 0.01820 sec
    166. if __name__=="__main__":
    167. # test_derivation_methods_time()
    168. test_Steger()

    我目前也没有弄明白为什么会出现三条线,似乎是中心线和轮廓,后续如果有改进会再次的更新

  • 相关阅读:
    simulink基础学习笔记
    Cadence23学习笔记(三)
    模型训练中的常见超参数解析
    【探索C语言中VS调试技巧】:提高效率和准确性
    python求多叉树任意两点之间的距离
    led台灯哪个牌子最好?2022最新的台灯牌子排名
    计算机毕业设计之java+ssm家校通网站
    角速度变化时四元数和旋转矩阵微分方程的证明
    与创新者同行!Apache Doris 首届线下峰会即将开启,最新议程公开!|即刻预约
    ubuntu14.04改静态ip
  • 原文地址:https://blog.csdn.net/m0_62919535/article/details/133245716