本文记录《机器视觉》 第三章第二节 —— 简单几何性质，一些学习笔记和个人理解，其中核心内容为二值图的转动惯量求解。
我们已经有了一组二值图，我们可以根据二值图来确定其表示物体的简单几何性质。
特征函数
二值图的特征函数 b(x, y)比较简单，当[x, y]处有物体时值为1，否则为0
面积
- 可以用特征函数求得二值图的面积A:
可以认为是二值图的 0 阶矩的物理意义。
质心
空间位置按照密度加权平均即是质心的位置 (bar{x}, bar{y}) ：
可以认为是二值图的 1 阶矩（静力矩）物理意义。
朝向
- 如果我们想要知道二值图物体表示的朝向，则需要用到转动惯量的概念。如果找到了使得物体转动惯量最小的轴，那么这个轴向就是物体的朝向。
- 在当前图像为二维的情况下，转动惯量是物体针对某条直线，将物体上的每个点到直线距离的平方按照密度计算积分，即得到了图像关于该轴向的转动惯量值。
- 我们的任务是为给定的二值图物体找到使得其转动惯量最小的直线。
- 转动惯量计算方法：
- 其中 r 表示二值图上的点到直线的距离，虽然还没有这条直线
直线建模
- 为我们的目标直线建模，取2个参数 原点到直线的距离 rho 直线和 x 轴之间（沿逆时针方向）的夹角 theta
- 这种建模方式有一些方便之处：
- 当坐标系平移或旋转时，这两个参数的变化是连续的
- 当直线平行（或近似平行）于某个坐标轴时，用这两个参数来表示直线也不会产生问题（相比于：使用斜率和截距来表示直线的情况）
- 使用这两个参数，可以将直线方程写为如下形式：
- 在直线上，距离原点最近的点 (-rho sin theta, rho cos theta) ，通过这个点，沿着夹角 theta 运动任意距离s 的点仍在直线上，因此可以将直线上任意一点(x_0,y_0)表示为：
最短距离
- 回到我们的二值图，在给定直线方程的情况下，二值图上一点(x,y)，直线上距离其最近的点(x_0,y_0)，二者距离显然可以表示为：
- 将直线上点公式eqref{6}代入，得到：
- 对于每个点，我们都需要解eqref{8}这样的优化方程
- 即给定了x,y,rho,theta，求解使得r最小的s,我们在eqref{8}中对s求导，可得：
- 将eqref{9}带入eqref{6}，二值图上(x,y)到直线上距离最近的点(x_0,y_0)可得到关系为：
- 将eqref{10}带入 eqref{7}可得：
- 此处可以看到，将某点(x,y)带入eqref{5}，得到值的绝对值即为该点到直线的垂线（最短）距离。
二阶矩轴向通过质心
- 我们已经得到了二值图上一点到任意直线的距离计算方法，将eqref{11}带入eqref{4}，得到：
- 对rho求导，并令倒数为0，得到：
- 其中，A是区域面积，而是区域质心。因此，我们得到了结论：
确定轴向倾角
我们已经确定该轴经过一个确定的点 (bar{x}, bar{y}) 了，仅需要再确定直线倾角即可。
- 将二值图平移到原点与质心重合的位置，那么我们要求得的就是一条穿过原点的直接倾角
- 也就是直接去除 rho 参数的影响
- 转动惯量计算方式如下：
其中，我们定义平移后的二值图I’上点的坐标为 ({x’}, {y’}) 。
- 可以表示为：
- 即：
- 上式对theta求导，并令求导结果等于零
- 假设a≠c，我们可以得到：
- 因此除非出现 b=0 并且 a=c 的情况, 否则, 我们最终可以得到：
- 至此我们已经求出了使得该二值图转动惯量最小和最大的两个轴
- E 的的最小值和最大值的比值，给出了一些关于物体有“多么圆”的信息。对于直线，这个比值是0对于圆，这个比值是1。
拉格朗日
从式eqref{15}开始，事实上我们要解的就是一个带约束的优化方程组，可以使用拉格朗日乘数法求解：
- 将E设为f(x,y)，约束条件设为g(x,y)=0，构建拉格朗日方程：
- 构造拉格朗日方程组：
- 重新令 x = sintheta, y=contheta可得：
- 即推出eqref{17}相同结论。
特征向量
可以将eqref{19}看作是一个二次型优化问题，原带约束的方程可以写成：
- 那么拉格朗日方程可以写成：
- L对bf{s}的偏导数为0：
- 而式eqref{25}就是在寻找矩阵bf{A}的特征向量和特征值。
- 也就是说，对于给定的二值图，求解其对应的a,b,c，构造出矩阵bf{A}，求解bf{A}的特征向量即是寻找最大、最小转动惯量的方向。
- 二者大小的比值也类似于特征值的比值，也就是矩阵的条件数。
参考示例
- 示例图像：
- 示例程序
import math
import cv2
import numpy as np
import matplotlib.pyplot as plt
from numpy.lib.function_base import iterable
def vvd_round(num):
if iterable(num):
return np.round(np.array(num)).astype('int32').tolist()
return int(round(num))
def show_image(image):
plt.imshow(image.astype('uint8'))
plt.show()
pass
def gravity_center(mask):
Ys, Xs = mask.nonzero()
A = (mask > 0).sum()
C_X = (Xs).sum() / A
C_Y = (Ys).sum() / A
return C_X, C_Y
def load_gray_image(image_path):
image = cv2.imread(image_path)
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
image = (image == 0).astype('uint8') * 255
return image
def moment_of_inertia(mask, center):
temp_image = mask.copy().astype('uint8') * 128
C_X, C_Y = center
Ys, Xs = mask.nonzero()
Ys = (Ys - C_Y) / 100
Xs = (Xs - C_X) / 100
a = (Xs * Xs).sum()
b = 2 * (Xs * Ys).sum()
c = (Ys * Ys).sum()
if b == 0:
theta = 0
elif a == c:
theta = - np.pi * 0.5 * 0.5
else:
theta = - math.atan(b / (a - c)) / 2
point_1 = vvd_round([C_X math.cos(theta) * 200, C_Y - math.sin(theta) * 200])
point_2 = vvd_round([C_X - math.cos(theta) * 200, C_Y math.sin(theta) * 200])
temp_image = cv2.line(temp_image.astype('uint8'), point_1, point_2, 255, 2)
point_1 = vvd_round([C_X math.cos(theta 0.5 * np.pi) * 200, C_Y - math.sin(theta 0.5 * np.pi) * 200])
point_2 = vvd_round([C_X - math.cos(theta 0.5 * np.pi) * 200, C_Y math.sin(theta 0.5 * np.pi) * 200])
temp_image = cv2.line(temp_image.astype('uint8'), point_1, point_2, 200, 2)
theta_1 = theta
theta_2 = theta 0.5 * np.pi
E_1 =(math.sin(theta_1)) ** 2 * a - b * math.sin(theta_1) * math.cos(theta_1) c * (math.cos(theta_1)) ** 2
E_2 =(math.sin(theta_2)) ** 2 * a - b * math.sin(theta_2) * math.cos(theta_2) c * (math.cos(theta_2)) ** 2
return temp_image, min(E_2, E_1) / max(E_2, E_1, 1)
if __name__ == '__main__':
image_path = 'test.png'
image = load_gray_image(image_path)
center = gravity_center(image)
temp_image, rate = moment_of_inertia(image, center)
show_image(temp_image)
pass
- 示例效果
参考资料
- 伯特霍尔德・霍恩著BERTHOLDKLAUSPAULHORN. 机器视觉[M]. 中国青年出版社, 2014.
