1 引言
最近使用传统方法进行石头检测和圆的检测的文章,受到了大家一致的好评.
嗯嗯,应该是一致的好评!
今天我们来研究一个新的方向,如何使用传统方法来进行大米的分割计数和大米长轴方向的标记?
如上图所示,我们拍摄了一张包含大米的图像,问题是如何使用图像处理的方法来求出图中共有多少颗大米?
嗯嗯,学过图像处理的同学,好多都碰到过类似问题吧,
碰巧我也接触了一些图像处理的技巧,那么我们就来尝试解决这个问题吧.
2 解决方案
2.1 读取图像
这里我们直接以灰度方式读取大米的图像,代码如下:
# Step 1: Read image
image_file = './images/rice.png'
image = cv2.imread(image_file, cv2.IMREAD_GRAYSCALE)
结果如下:
2.2 生成背景图
注意观察上图,可以看出图像中的亮度不均匀,即上图中背景的上半部分的亮度相比下半部分的亮度偏亮,因此我们使用大kernel的模板来进行开运算生成背景图,代码如下:
selem = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (31, 31))
background = cv2.morphologyEx(image, cv2.MORPH_OPEN, selem)
结果如下:
2.3 获取前景图
我们有了原图,有了生成的前景图,将二者相减,即可得到前景图.代码实现如下:
foreground = cv2.subtract(image, background)
运行结果如下:
可以看出,经过处理后的图像背景亮度基本趋于一致,不会再有亮度不均匀的现象.
2.4 二值化
为了进一步区分前景大米和背景,这里我们采用OTSU的方法来进行二值化操作,方法如下:
thresh_otsu, binary_image = cv2.threshold(foreground.astype(np.uint8), 0, 255,
cv2.THRESH_BINARY+cv2.THRESH_OTSU)
结果如下:
可以看出,二值化后前景物体大米都被分割出来了.
2.5 连通域标记
接下来我们使用opencv中的connectedComponents 函数来进行连通域的标记,代码如下:
num_labels, label_image = cv2.connectedComponents(binary_image, connectivity=4)
解释如下:
num_labels：所有连通域的数目
labels_image：图像上每一像素的标记，用数字1、2、3…表示（不同的数字表示不同的连通域）
将label_image 进行可视化,结果如下:
2.6 查看单个米粒
经过上述操作,我们已经将每颗米粒使用不同的颜色标记出来,接下来我们来挑选单个米粒来进行标记结果校验.这里我们挑选第43颗米粒来进行展示,代码实现如下:
# First, we select the grain labeled 43, and finds its contour
single_grain_full = (label_image == 43)
cnt_list, _ = cv2.findContours(single_grain_full.astype(np.uint8),
cv2.RETR_EXTERNAL,
cv2.CHAIN_APPROX_NONE)
# As there is only one contour, we select it, and use it to compute a
# bounding box around the grain
grain_contour = cnt_list[0]
y_start, x_start, dy, dx = cv2.boundingRect(grain_contour)
grain_min_bbox = label_image[x_start:x_start+dx, y_start:y_start+dy].astype(np.uint8)
# grain_min_bbox is the minimal bounding box around the grain, so we pad it a
# bit before we display it.
grain_bbox = cv2.copyMakeBorder(grain_min_bbox, 10, 10, 10, 10,
cv2.BORDER_CONSTANT, value=0)
主要流程为:
- 首先将label_image中标记为43号的区域取出来,得到图像single_grain_full
- 接着使用cv2.findContours 找到该米粒的外轮廓
- 然后利用函数cv2.boundingRect 求出该外轮廓的外接矩形
- 最后将外接矩形扩充10个像素,得到最后的结果图
得到第43颗米粒的效果图如下:
2.7 图像的矩
观察我们标记后的图像,可以看出所有米粒均被我们处理成一个个独立的封闭区域,接下来我们来求对应每个目标米粒的矩,进而确定每个目标相应的性质,关于矩的定义和性质参考前文.代码实现如下:
moment_dict = moments.get_moment_list(label_image)
我们查看单个目标相关性质的输出:
print(moment_dict[1])
输出如下:
{'m00': 61.0, 'm01': 204.0, 'm10': 1392.0, 'm11': 4502.0,
'm02': 1020.0, 'm12': 22208.0, 'm21': 100212.0, 'm20': 32118.0, 'm22': 486906.0, 'mu00': 61.0, 'mu01': 0, 'mu10': 0, 'mu11': -153.21311475409857,
'mu20': 353.01639344262367, 'mu02': 337.7704918032787,
'mu21': -206.4703036818537, 'mu12': -43.29588820209756,
'nu00': 1, 'nu01': 0, 'nu10': 0, 'nu11': -0.0411752525541786, '
nu02': 0.09077411765742506, 'nu20': 0.09487137689938825,
'nu12': -0.0014897797471970042, 'nu21': -0.007104491664128479,
'area': 61.0, 'cx': 22.81967213114754, 'cy': 3.3442622950819674,
'cov_mat': array([[ 5.78715399, -2.51169041],
[-2.51169041, 5.53722118]]), 'theta': -0.7605417074824052}
相关解释:
- m00…m22 为单个目标对应的0阶,1阶,2阶,3阶,4阶原始矩
- mu00…mu12 为单个目标对应的0阶,1阶,2阶,3阶中心矩
- nu00…nu12 为单个目标对应的0阶,1阶,2阶,3阶正则中心矩
- area 为目标的面积
- cx,cy 为目标的质心
- cov_mat 为协方差矩阵
- theta 为目标长轴方向
2.8 结果可视化
经过上述运算,我们得到了每粒大米的面积,质心和长轴方向,这里我们将其可视化表示,代码如下:
# Superimplose labels, centriods and orientation over this image
fig = plt.figure(fig_num)
plt.imshow(image, cmap='gray', interpolation='none')
plt.xticks([]), plt.yticks([])
plt.tight_layout()
obj_properties = {}
for ind, obj_moments in moment_dict.items():
if obj_moments['m00'] > 0:
# Area
area = obj_moments['m00']
theta = obj_moments['theta']
cx = int(obj_moments['m10'] / obj_moments['m00'])
cy = int(obj_moments['m01'] / obj_moments['m00'])
rho = 10
dx = rho*np.cos(theta)
dy = rho*np.sin(theta)
if abs(theta*180/np.pi) <= 20:
plt.plot([cx-dx, cx, cx+dx], [cy-dy, cy, cy+dy], 'go-')
else:
plt.plot([cx-dx, cx, cx+dx], [cy-dy, cy, cy+dy], 'bo-')
plt.text(cx, cy, str(ind), color='red', fontsize=14)
# Store properties
props = {}
props['area'] = area
props['cx'] = cx
props['cy'] = cy
props['theta'] = theta*180/np.pi
obj_properties[ind] = props
print(obj_properties[1])
axes = plt.gca()
axes.set_xlim([0, label_image.shape[0]])
axes.set_ylim([label_image.shape[1], 0])
上述代码中的运行结果如下:
上图中:
- 每个大米目标的长轴标记出来,长轴的中间点表示大米的质心位置;
- 长轴角度介于[-20,20]度的使用绿色标记,其他角度的使用蓝色标记出来
- 每个大米的id使用红色字体标记出来.
2.9 统计分析
上述2.8节的代码中,将每个大米的面积,质心位置和朝向都记录下来,并存储在obj_properties中,接下来,我们对角度和面积进行展示并计算相应的直方分布图,代码如下:
areas = []
orientations = []
print('Object Area Orientation Centroid')
for ind, props in obj_properties.items():
areas.append(props['area'])
orientations.append(props['theta'])
print('{0:>6} {1:>9,.0f} {2:>11,.2f} ({3:>3}, {4:>3})'.format(ind,
props['area'],props['theta'],
props['cx'],props['cy']))
fig_num = plot_line_and_histogram(areas, title='Areas',marker='*',linestyle='',
fig_num=fig_num, num_bins=30,
normed=False, label='Area (pixels)')
fig_num = plot_line_and_histogram(orientations, title='Orientations',marker='*', linestyle='',
fig_num=fig_num, num_bins=30,
normed=False, label='Angle (degrees)')
# Index + 1 because of zero indexing in list, while the object labels are
# enumerated from 1 (since label 0 is reserved to backround).
print('Grain {} has the largest area of {} pixels'.format(np.argmax(areas)+1,
np.max(areas)))
print('Mean area: {0:.2f}'.format(np.mean(areas)))
运行结果如下:
Object Area Orientation Centroid
1 61 -43.58 ( 22, 3)
2 130 -40.53 ( 58, 5)
3 186 -80.62 ( 79, 11)
4 5 5.15 ( 88, 0)
5 17 -9.17 (149, 0)
6 119 -64.62 (163, 7)
7 187 -54.87 (188, 9)
8 46 -18.47 (204, 2)
9 91 -86.96 (217, 6)
10 118 -24.79 (238, 4)
我们一共标记出101个大米,这里仅展示前10个大米相应的属性
面积直方图如下:
上图中,上半部分,长轴为标记出的101个大米,每一个大米对应的面积取值.下半部分表示面积分布的直方图.
角度直方图如下:
上图中,上半部分,长轴为标记出的101个大米,每一个大米对应的角度取值.下半部分表示长轴角度分布的直方图.
3 总结
本文介绍了使用传统图像处理方法进行大米图像中大米的计数和方向标记的基本方法,并给出了完整的代码实现.
您学废了吗?
关注公众号《AI算法之道》，获取更多AI算法资讯。
注: 关注公众号,后台回复 大米 , 可获取完整代码
315
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
