简 介： 利用对所获得的圆环边缘亮度平均变化曲线，利用Sigmoid函数进行建模逼近，可以对原来经由HoughCircle所获得的半径进行补偿。对于利用扫描仪所获得的图片进行实际处理，可以看到它可以获得与利用轮廓面积所测量的结果精度接近的效果。但精度略差一些。

关键词： 亮度曲线建模，Sigmoid函数，抑菌圈面积测量

 

§01 问题提出

---

1.1 背景介绍

  在 抑菌圈金属模板 圆环直径提取过程中，使用HoughCircles可以比较方便求取图像中符合在一定直径范围内的圆形区域，但限于 Hough Transforms 本身的缺陷，对于圆形直径的分辨率比较低。通过 轮廓面积 来求取圆形半径，可以提升圆形直径计算精度，但是会发现这样所获取的圆形的直径受到图像在二值化过程中的阈值影响比较大。

  下图显示了对于 模板在扫描以上 不同位置下（实际上是模板沿着扫描仪的对角线直线移动）所拍摄的同样模板通过轮廓反过来计算四个圆形直径对应的变化。

▲ 图1.1 利用轮廓面积获取圆形直径的数据

  在博文 轮廓面积中通过手动调整阈值，可以看到所计算出来的圆形直径会产生大约1 ~ 2像素的变化。这是算法对于图像亮度所不适应性带来的影响。反过来，这种收到的影响也会受到扫描仪在不同横向和总线位置上的补光灯的亮度差也会带来对应的影响。

  那么是否可以找到一种更加有效的方式，来弥补这种由于阈值（或者亮度）对于测量结果所带来的差异吗?

1.2 边缘建模

  一种解决方式就是对于模板图像中圆环的边缘灰度图像进行建模，找到其中对于亮度不敏感，而是以来在于实际物体物理反射特性的参数，用于估算圆形直径。

1.2.1 分析变化变化特性

  下图反映了图像中圆环边缘处的图像细节：

• 图像存在一个亮度变化过渡带；
• 像素离散化带来的马赛克效应明显；

▲ 图1.2.1 圆形边缘像素细节

  变化亮度变化的过渡带表明，随着对该图像进行二值化，阈值选取不同，就会对于变化的像素在分割上出现差异，从而带来所的带到的圆形直径发生偏转。

  下图给出了在圆形边缘处的灰度变化情况。

▲ 图1.2.2 圆形灰度图边缘变化

  由于像素离散所带来的灰度变化的台阶，可以采用 数据插值的方式 对数据进行连续平滑处理。考虑到我们是对于圆周（半径在90， 110）所对应的大约 628个像素进行处理，所以一开始先不进行插值进行统计。

1.2.2 圆环灰度统计

  具体处理方法如下：

  首先根据最初通过 HoughCircles获得圆环的中心 $\left(x_0,y_0\right)$ 以及半径 $r_0$ 。选择合适的扫描区间： $\left(r_0-\Delta r,r_0+\Delta r\right)$ 和扫描步数 $N_0$ ，确定一系列的圆：

• 圆心： $\left(x_0,y_0\right)$
• 第n个圆半径： $r_n\in \left[r_0-\Delta r,r_0+{\Delta }_r\right]$ 。

r n = 2 ⋅ Δ r N ⋅ n + r 0 − Δ r ,     n ∈ { 1 , 2 , ⋯   , N − 1 } r_n = {{2 \cdot \Delta _r } \over N} \cdot n + r_0 - \Delta _r ,\,\,\,n \in \left\{ {1,2, \cdots ,N - 1} \right\} rn​=N2⋅Δr​​⋅n+r0​−Δr​,n∈{1,2,⋯,N−1}

▲ 图1.2.3 进行圆环灰度统计

  对于第 $n$ 圆进行圆周采样，采样点 $\left(a_{ni},b_{ni}\right)$ 为：
$$a_{ni}=\mathrm{int}\left(x_0+r_n\cdot \cos \left({\theta }_i\right)\right),b_{ni}=\mathrm{int}\left(y_0+r_n\cdot \sin \left({\theta }_i\right)\right)$$

  其中角度 ${\theta }_n=\frac{n\cdot 2\pi }{N}$ ，N是每个圆采样点数。

  计算采样点对应的弧度图像平均亮度作为采样圆周的平均亮度：$$L_n=\frac{1}{N}\sum _{i=1}^{N}img\left[a_{ni},b_{ni}\right]$$

  最终绘制这所有采样圆亮度曲线，然后进行建模估计圆的半径。

 

§02 数据处理

---

2.1 获取边缘亮度变化曲线

2.1.1 处理程序

#!/usr/local/bin/python
# -*- coding: gbk -*-
#============================================================
# PROC1.PY                     -- by Dr. ZhuoQing 2022-01-27
#
# Note:
#============================================================
from headm import *                 # =
import cv2
npzfile = '/home/aistudio/work/Scanner/ScanDiagBlock.npz'
imagedir = '/home/aistudio/work/Scanner/ScanDiagBlock'
#filedim = sorted([s for s in os.listdir(imagedir) if s.find("jpg") > 0])
#printt(filedim:)
alldata = load(npzfile, allow_pickle=True)
alldim = alldata['alldim']
allfile = alldata['allfile']
#printt(alldim:, allfile:)
#------------------------------------------------------------
CIRCLE_NUM          = 100
DELTA_RATIO         = 7
SAMPLE_NUM          = 300
theta = linspace(0, 2*pi, SAMPLE_NUM)
def circleEdges(img, circles):
    clight = []
    for c in circles:
        rdim = linspace(c[-1]-DELTA_RATIO, c[-1]+DELTA_RATIO, CIRCLE_NUM)
        rmean = []
        for r in rdim:
            thetadim = []
            for a in theta:
                x = int(c[0] + r*cos(a))
                y = int(c[1] + r*sin(a))
                thetadim.append(img[y,x])
            rmean.append(mean(thetadim))
        clight.append(rmean)
    return clight
#------------------------------------------------------------
for id,f in enumerate(allfile):
    cc = alldim[id]
    img = cv2.imread(f)
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    ce = circleEdges(gray, cc)
    plt.clf()
    plt.figure(figsize=(10,10))
    for i in range(4):
        plt.subplot(2,2,i+1)
        plt.plot(ce[i])
        plt.xlabel("n")
        plt.ylabel("light")
        plt.grid(True)
        plt.tight_layout()
    plt.savefig('/home/aistudio/stdout.jpg')
    plt.show()
    break
#------------------------------------------------------------
#        END OF FILE : PROC1.PY
#============================================================

2.1.2 处理结果

处理参数：

M：100
DeltaR：7
N：300

▲ 图2.1.1 获得四个圆环变换亮度变化曲线

  这是对于所采集到的100张模板图片进行处理之后所获得边缘亮度变化曲线图。

▲ 图2.1.2 100个图片的数据

2.2 边缘建模

2.2.1 选择Sigmoid函数

  对于边缘亮度变化使用一个Sigmoid函数进行建模。

  下面是一个普通的Sigmoid函数。 如果使用该函数去逼近圆形边缘亮度曲线，还需要有四个参数。

$$f\left(x\right)=\frac{1}{1+e^{-\lambda x}}$$

▲ 图2.2.1 Lambda=2对应的Sigmoid函数

  这四个参数分表代表了对于曲线在横行和纵向方向上的拉伸和平移。

待定系数：

a：纵向上的尺度，初始化为：140
b：纵向上的平移，初始化为：20
c：横向方面的平移：50
d：横向是好的尺度：初始化为：0.2

$$g\left(x\right)=\frac{a}{1+e^{-d\left(x-c\right)}}+b$$

  下面是使用上述初始化值绘制出的Sigmoid曲线。针对每一个圆环的亮度曲线，使用scipy的曲线拟合，获取最佳参数。

▲ 图2.2.2 应用初始化值绘制的图像

  其中参数 $c$ 代表了曲线的实际上中心平移。在根据前面对于圆环扫描的参数，可以将原来HoughCircles估计的圆环半径进行修正。

$$r_m=r_0+\left(c-0.5N\right)\frac{2\Delta r}{N}$$

2.2.2 对实际边缘亮度数据进行建模

▲ 图2.2.3 利用Sigmoid对于测量曲线进行建模后的曲线

2.3 边缘处理后的曲线

2.3.1 处理后的结果

  下面是对于曲线进行亮度估计之后所得到得到的圆环直径变化曲线。

▲ 图2.3.1 经过建模之后所获得的模板圆环直径变化

maxvalue: [112.49083917 112.39641417  94.72715966  94.83673293]
minvalue: [111.7984239  111.65632129  93.62277949  93.76408117]
deltavalue: [0.69241527 0.74009288 1.10438017 1.07265176]

  作为对比，下面是直接利用HoughCircle所得到的各个圆的半径曲线：

▲ 图2.3.2 利用HoughCircle所得到的各个圆的半径

maxvalue: [113.3 114.4  96.   96.1]
minvalue: [110.  110.8  92.5  92.8]
deltavalue: [3.300003  3.5999985 3.5       3.2999954]

  通过上述对比，可以看到通过对于边界建模补偿之后，大大提高了圆环半径估计的精度。

2.3.2 利用轮廓面积处理结果

  下面是在 利用圆圈轮廓面积求取圆环半径：cv2.findContours, contourArea 利用轮廓面积反过来计算圆环直径所得到的结果。

▲ 图2.3.3 二值化阈值为100时，采用轮廓面积计算圆环半径的结果

max(rdim1)-min(rdim1): 0.45511065871777134	
max(rdim2)-min(rdim2): 0.5779160789730753
max(rdim3)-min(rdim3): 0.9834303156866042
max(rdim4)-min(rdim4): 0.8077211123899559

  对比使用变换灰度建模方法所得到的半径结果与通过边缘轮廓面积所得到的结果来看，它们所得到的结果很相似。但利用轮廓面积所得到的结果要更好一些。

 

※ 测试总结 ※

---

  利用对所获得的圆环边缘亮度平均变化曲线，利用Sigmoid函数进行建模逼近，可以对原来经由HoughCircle所获得的半径进行补偿。对于利用扫描仪所获得的图片进行实际处理，可以看到它可以获得与利用轮廓面积所测量的结果精度接近的效果。但精度略差一些。

  那么剩下一个问题：是否可以将边缘建模与轮廓面积合二为一，这样就可以提高圆环半径测量精度了呢？

---

■ 相关文献链接:

• 扫描仪标准模板滑动采集图像及其处理
• Hough Transforms
• 利用圆圈轮廓面积求取圆环半径：cv2.findContours, contourArea
• 获取棋盘格与标准模板在扫描仪上不同位置图片
• 测试scipy.interpolate.interp2d对于低分辨率图像插值处理结果

● 相关图表链接:

• 图1.1 利用轮廓面积获取圆形直径的数据
• 图1.2.1 圆形边缘像素细节
• 图1.2.2 圆形灰度图边缘变化
• 图1.2.3 进行圆环灰度统计
• 图2.1.1 获得四个圆环变换亮度变化曲线
• 图2.1.2 100个图片的数据
• 图2.2.1 Lambda=2对应的Sigmoid函数
• 图2.2.2 应用初始化值绘制的图像
• 图2.2.3 利用Sigmoid对于测量曲线进行建模后的曲线
• 图2.3.1 经过建模之后所获得的模板圆环直径变化
• 图2.3.2 利用HoughCircle所得到的各个圆的半径
• 图2.3.3 二值化阈值为100时，采用轮廓面积计算圆环半径的结果

测量程序

#!/usr/local/bin/python
# -*- coding: gbk -*-
#============================================================
# PROC2.PY                     -- by Dr. ZhuoQing 2022-01-27
#
# Note:
#============================================================
from headm import *                 # =
from scipy.optimize        import curve_fit
#------------------------------------------------------------
cealldiag = '/home/aistudio/work/Scanner/cealldiag.npz'
npzdiag = '/home/aistudio/work/Scanner/ScanDiagBlock.npz'
alldata = load(npzdiag, allow_pickle=True)
alldim = alldata['alldim']
#printt(alldim:)
cealldata = load(cealldiag, allow_pickle=True)
ceall = cealldata['ceall']
#printt(len(ceall))
#------------------------------------------------------------
def sigmoidx(x, a, b, c, d):
    return a/(1+exp(-(x-c)*d)) + b
#------------------------------------------------------------
#printt(ceall[0])
delta = 7
N = 100
ceallmod = []
for id,ce in enumerate(ceall):
#    plt.clf()
#    plt.figure(figsize=(10,10))
    rdim = alldim[id]
#    printt(rdim)
    x = linspace(0, 100, 100, endpoint=False)
    cmod = []
    for i,c in enumerate(ce):
        y = c
        param = (140,20,50,0.2)
        param, conv = curve_fit(sigmoidx, x, y, p0=param)
        r0 = rdim[i][-1]
        rmod = r0+(param[2]-N/2)*(2*delta/N)
        cmod.append(rmod)
#        printf(param)
#        yfit = sigmoidx(x, *param)
#        plt.subplot(2,2,i+1)
#        plt.plot(c, linewidth=4, label='Origin')
#        plt.plot(yfit, linewidth=2, label='Fit')
#        plt.legend(loc='upper left')
#        plt.xlabel("n")
#        plt.ylabel("Value")
#        plt.grid(True)
#        plt.tight_layout()
#    plt.savefig('/home/aistudio/stdout.jpg')
#    plt.show()
    ceallmod.append(cmod)
#    break
printt(ceallmod:)
#------------------------------------------------------------
#ceallmod = [cc.T[2] for cc in alldim]
 #------------------------------------------------------------
carray = array(ceallmod).T
plt.clf()
plt.figure(figsize=(10,8))
plt.plot(carray[0], label='Circle1')
plt.plot(carray[1], label='Circle2')
plt.plot(carray[2], label='Circle3')
plt.plot(carray[3], label='Circle4')
plt.xlabel("n")
plt.ylabel("Ratio")
plt.grid(True)
plt.tight_layout()
plt.savefig('/home/aistudio/stdout.jpg')
plt.show()
#------------------------------------------------------------
maxvalue = array([max(c) for c in carray])
minvalue = array([min(c) for c in carray])
deltavalue = maxvalue-minvalue
printt(maxvalue:, minvalue:, deltavalue:)
#------------------------------------------------------------
#        END OF FILE : PROC2.PY
#============================================================