上阈值—面积曲线的拟合与最佳阈值选取

书名：数字图像处理实战
作者：杨坦 张良均
出版社：人民邮电出版社有限公司
出版时间：2023-11-01
ISBN：9787115623850

---

第9章　钢轨表面缺陷检测

9.4 基于区域生长算法的钢轨表面缺陷检测

9.4.2 上阈值的自适应选择

3.上阈值—面积曲线的拟合与最佳阈值选取

• 计算得到上阈值—面积曲线后，本章采用有a、b、c和d共4个参数的Sigmoid函数对上阈值—面积曲线进行拟合，如式(9-1)所示。

  
  

• 这4个参数联合在一起可以控制曲线两个平滑段的高度、跳跃段的水平位置和陡峭程度。该Sigmoid函数的一阶导数的定义如式(9-2)所示。

  
  

• 对样本rail_1中缺陷的上阈值—面积曲线进行拟合，得到的Sigmoid函数及其导数的示意如图9-15所示。

  
  图9-15 Sigmoid函数及其导数的示意

• 观察图9-15中的函数曲线，当横坐标x在(0，20)区间时Sigmoid曲线处于第一个平滑段，这一段的y值由参数a确定。这个区间对应于钢轨表面坑洞陡峭的内壁中的像素灰度范围，所以此时生长出的连通区域的面积随上阈值的增加只产生微小的变化。
  在x从20增加到25的过程中，曲线有一个显著的跳跃，阶跃的位置由参数c确定，阶跃的陡峭程度由参数d确定。参数c的值大致对应于坑洞边缘像素的灰度值，上阈值超过c后连通区域将迅速扩展到坑洞周围的平坦的背景区域，所以生长出的连通区域的面积随着上阈值的微小增加而急剧增大。
  当x大于25后，曲线再次变得平滑，因为连通区域已经扩展到背景区域了，所以此时增大上阈值很难再引起y值的增大。
  根据上面的分析，一个常规的思路是将上阈值设置为Sigmoid曲线上阶跃结束位置附近的点的横坐标，对应于缺陷边缘的像素的灰度值。但这个思路也存在问题，由于不同种子点的上阈值—面积曲线拟合得到的Sigmoid曲线千差万别，很难从曲线本身得到可以定量处理的上阈值选择策略。

• 图9-15中也展示了Sigmoid函数的一阶导数，它在Sigmoid曲线的两个平滑段的值接近于0，在Sigmoid曲线的阶跃段表现为一个向上的波峰。通过在Sigmoid函数一阶导数曲线上选择波峰右侧下降段上特定y值对应的x，可以定量选取区域生长算法的上阈值。
  经过实践验证，本章中取导数值为0.1时的上阈值作为区域生长算法的最终上阈值。

• 在得到拟合Sigmoid函数需要的数据后，使用op.curve_fit函数来估计Sigmoid函数的参数a、b、c和d，完整的实现过程如代码9-10所示。

• 代码9-10 估计Sigmoid函数的参数

#定义带参数的sigmoid函数
def sigmoid(x,a,b,c,d):
  y=a+1/(b+np.exp(c-d*x))
  return y

x_data = np.arange(len(data))
y_data = data
curve =op.curve_fit(sigmoid,x_data,y_data, bounds=([-5,00，-1]，[5，0.1,20,1]))
#使用op.curve fit函数估计sigmoid函数的参数
a,b,c,d= curve[0]
x=np.arange(0，len(data)，0.1)
y=a+1/(b+np.exp(c-d*x))    # 最终得到的sigmoid函数表达

• 得到拟合曲线的表达式后，绘图并使用SciPy库的root函数求解Sigmoid函数的导数值为0.1时的横坐标，完整的实现过程如代码9-11所示，结果如图9-16所示。
  图中使用竖线表示求得的上阈值与Sigmoid函数的关系。
• 代码9-11 根据Sigmoid函数计算上阈值

from scipy.optimize import root,fsolve
# 定义sigmoid函数的导数
def f1(x):
  return(d*np.exp(c-d*x))/(np.exp(c-d*x)+ b)**2

#定义导数与0.1的差值，用于生成方程
def f2(x):
  return f1(x)-0.1

ans_fsolve =fsolve(f2,[25])[0]
print(ans_fsolve)
# 绘制导数的图像
plt.plot(x,y,color='red',label='curve_fit\na={0:.2f} b={1:.2f}\nc={2:.2f}d={3:.2f}'.format(a,b,c,d))
plt.plot(x，f1(x),color='blue',label='F\'(x)')
plt.legend(bbox_to_anchor=(0,1),loc=2,borderaxespad=0)
plt.title('sigmoid:Rail {θ},F\'(x)'.format(picture_label))
plt.savefig('../tmp/'+'sigmoid Rail {o} daoshu.png'.format(picture label),dpi=400)
plt.show()

图9-16 上阈值位置示意图（基于样本rail_1）

• 本小节展示了如何一步步地确定单个种子的生长算法中的上阈值，将上阈值—面积曲线的计算、Sigmoid函数拟合以及上阈值计算封装成函数Calculation_of_threshold，如代码9-12所示。调用该函数可以得到针对一幅图像中所有种子对应的上阈值列表。
• 代码9-12 定义计算上阈值的函数Calculation_of_threshold

def Calculation of threshold(zhongzi, img):
  data = area_of_gray(zhongzi[e],zhongzi[1],img.copy(),yuzhi=70)
  x data =np.arange(len(data))
  y_data = data
  curve = op.curve fit(sigmoid,x data, y data, bounds=([-5，0，0，-1]，[5，0.1,20，1]))
  a,b,c,d= curve[0]
  # ans=np.log((2*np.exp(c)+3**(1/2)*np.exp(c))/b)/d
  f1=lambda x:(d*np.exp(c-d*x))/(np.exp(c-d*x)+b)**2-0.1
  ans fsolve =fsolve(f1,[25])[0]
  return ans fsolve

threshold_of_seed =[]
for i in range(len(zhongzi)):
  threshold=Calculation_of_threshold(zhongzi[i],img=imggray4.copy())
  threshold_of_seed .append(threshold)