python图像处理教程：初步📷插值变换📷形态学处理📷滤波

光斑是工程中经常出现的图像数据，其特点是目标明确，分布清晰。对光斑图像的分析，主要包括质心定位、目标截取以及半径拟合等。

质心

所谓质心，就是质量的中心，连续物体的质心可表示为$\frac{\int xf(x)\text{d}x}{\int f(x)\text{d}x}$，则图像质心可表示为

$$\bar{x}=\frac{\sum i{p}_{ij}}{\sum p_{ij}}\bar{y}=\frac{\sum j{p}_{ij}}{\sum p_{ij}}$$

由于图像是二维数组，有两个可以操作的方向，故而除了全局的质心外，其每一行或每一列均有一个质心，行质心和列质心的交叉点，就是图像的质心，图像如下，左侧为原始光斑，右侧为其行质心和列质心。

绘图代码如下，考虑到求质心时，有一个求和项需要做除数，所以在求质心之前，将图像加上$10^{-9}$，以确保不会出现除0错误，同时不会对质心造成干扰。

from numpy as np
import matplotlib.pyplot as plt

# 创建一副光斑图像
size = [256,256]
x, y = np.indices(size)
d2 = (x-80)**2+(y-80)**2
img = 255*np.exp(-d2/32)

# 求行质心和列质心
arr = np.arange(256).reshape(-1,1)

cImg = img + 1e-9
cenRow = (cImg @ arr).reshape(-1)/np.sum(cImg, 0)
cenCol = (cImg.T @ arr).reshape(-1)/np.sum(cImg, 1)

# 绘图
fig = plt.figure()
ax1 = fig.add_subplot(121)
ax1.imshow(img, cmap='gray')

ax2 = fig.add_subplot(122, sharey=ax1)
ax2.tick_params(axis="y", labelleft=False)
arr = arr.reshape(-1)
plt.imshow(img)
plt.plot(arr, cenRow)
plt.plot(cenCol,arr)
plt.show()

目标截取

由于光斑只占据图片中很小一部分，大量的冗余信息等同于噪声，会影响数据处理，故需截取感兴趣的区域。【plt.ginput】函数提供一种交互操作方法，可返回鼠标点击的位置，其输入参数为选取点数，输出为点击的点的坐标。根据返回的点的位置，即可对图像进行截取。

plt.imshow(img)
p1, p2 = plt.ginput(2)

roi = img[int(p1[1]):int(p2[1]), int(p1[0]):int(p2[0])]

在本例中，$p_1,p_2$两点分别在左上角和右下角，从而无论$x$还是$y$，$p_1$的值均小于$p_2$。另一方面，ginput的返回值是$(x,y)$坐标，而矩阵截取过程则是先行后列，行对应的是$y$，这就是roi在截取时的选区逻辑。

最终截取结果如下。

fig = plt.figure()
ax = fig.add_subplot(121)
plt.imshow(roi)

ax = fig.add_subplot(122, projection='3d')
yGrid, xGrid = np.indices(roi.shape)
ax.plot_surface(yGrid, xGrid, roi)
plt.show()

光斑半径

一般来说，Gauss光束以其强度最大值的$\frac{1}{e}$处为边界，也就是说，最简单的求解光斑半径的方法，就是统计所有边界内部的点的个数，然后再根据圆面积公式求取其半径。

area = np.sum(roi>255/np.e)
r = np.sqrt(area/np.pi)
# 5.56

另一种方法，则需对光斑的某一维度进行高斯拟合。具体方法是对roi按列求最大值，然后对得到的数组进行高斯拟合，最终得到半径值$r=5.65$，显然更接近于我们最开始设定的$\sqrt{32}$，其拟合结果如下。

from scipy.optimize import curve_fit
def gauss(x, a, b, c):
    return a*np.exp(-(x-b)**2/c**2)

y = np.max(roi,0)
x = np.arange(len(y))

abc,  para = curve_fit(gauss,x,y)
print(abc)							#即上式中的a,b,c
# [255.          19.           5.65685425]
Y = gauss(x,abc[0],abc[1],abc[2])	#拟合值
plt.scatter(x,y, label="origin")    #绘制原始数据的散点图
plt.plot(x,Y, label="curve_fit")	#绘制拟合数据的曲线图
plt.legend()
plt.show()