如何计算质心
本文介绍了质心的概念,以及基于Numpy、Scipy、OpenCV等工具的多种实现方式。
如何计算质心
原始文档:https://www.yuque.com/lart/idh721/gpbigm
概念
质心,即质量中心的简称。质点系的质心是质点系质量分布的平均位置。指物质系统上被认为质量集中于此的一个假想点,与重心不同的是质心不一定要在有重力场的系统中,值得注意的是除非重力场是均匀的,否则同一物质系统的质心与重心通常不在同一假想点上。
计算
质心坐标等于所有点关于每个坐标的以质量为权重的加权平均值。一般主要在二维空间讨论,尤其是图像数据,但是这里直接按照更一般的形式进行定义。首先对于任意 n n n维空间中的连续形式的子集 P P P的质心可以定义为:
C = ∫ p g ( p ) d p ∫ g ( p ) d p C = \frac{\int p g(p)dp}{\int g(p)dp} C=∫g(p)dp∫pg(p)dp
其中:
- p ∈ R n p \in \mathbb{R}^n p∈Rn表示该子集中的一点;
- g g g表示该子集的特征函数(indicator function or a characteristic function)。一个比较实际的场景是,它可以用来表示各个位置对应的质量。
也可以看到,这里的分母是对这个集合的一个度量,因而如果度量为0,那么就不可以被计算质心。
而对于第 k k k个坐标 C k C_k Ck的计算,可以通过如下形式:
C k = ∫ z S k ( z ) d z ∫ g ( x ) d z C_k = \frac{\int z S_k(z) dz}{\int g(x) dz} Ck=∫g(x)dz∫zSk(z)dz
这里 S k ( z ) S_k(z) Sk(z)表示的是对应的 P P P与由 p k = z p_k=z pk=z定义的超平面(hyperplane)的交集的度量,在这个超平面上,会涉及到其他所有的坐标轴。
对于一个平面图,即二维情形,上面的式子可以用于求取 C x C_x Cx和 C y C_y Cy:
- C x = ∫ x S x ( x ) d x A C_x = \frac{\int x S_x(x) dx}{A} Cx=A∫xSx(x)dx
- C y = ∫ y S y ( y ) d y A C_y = \frac{\int y S_y(y) dy}{A} Cy=A∫ySy(y)dy
即对单一轴向上的坐标积分,每个坐标对应乘以一个与之关联的量 S S S(该量会涉及到另一个轴), A = ∫ S x ( x ) d x = ∫ S y ( y ) d y A = \int S_x(x) dx = \int S_y(y) dy A=∫Sx(x)dx=∫Sy(y)dy表示图形的面积。一般情况下,我们可以简单的理解为这是过点 ( x , 0 ) (x, 0) (x,0)或是 ( 0 , y ) (0, y) (0,y)的垂线与图像区域的相交构成的线段的长度。
对于更实际的离散且有限点集的情形下,前面二维的形式可以转化为如下形式:
- C x = ∑ i w i x i ∑ i w i C_x = \frac{\sum_i w_i x_i}{\sum_i w_i} Cx=∑iwi∑iwixi
- C y = ∑ i w i y i ∑ i w i C_y = \frac{\sum_i w_i y_i}{\sum_i w_i} Cy=∑iwi∑iwiyi
这里要注意,公式中表示各个点的方式与前面直接基于坐标值的方式有所不同,而是通过一个额外的点索引 i i i来对不同的点进行编码排序,从而构建了公式。
对于不同的点有着不同的对应权重,我们可以理解为质量或者其他的度量形式,所以对应相同的点序号 i i i,其各个轴向上的坐标权重也是一样的 w i w_i wi,且 W = ∑ i w i W = \sum_i w_i W=∑iwi可以表示为图像对应的整体质量或者其他的度量。如果各点权重均为1,则这里的 W W W则实际上便是点的数量,对于数字图像而言,就是图形面积了。
更一般的,这里的点 i i i实际上还可以替换为有着面积(或者说离散点的数量) A i A_i Ai的区域 P i P_i Pi。计算过程中将其质心作为这里的点。
- C x = ∑ i C i , x W i ∑ i W i = ∑ i ∑ j x i , j w i , j ∑ j w i , j ∑ j w i , j ∑ i ∑ j w i , j = ∑ i ∑ j x i , j w i , j ∑ i ∑ j w i , j = ∑ l x l w l ∑ l w l C_x = \frac{\sum_i C_{i,x} W_i}{\sum_i W_i} = \frac{\sum_i \frac{\sum_j x_{i,j} w_{i,j}}{\sum_j w_{i,j}} \sum_j w_{i,j} }{\sum_i \sum_j w_{i,j}} = \frac{\sum_i \sum_j x_{i,j} w_{i,j}}{\sum_i \sum_j w_{i,j}} = \frac{\sum_l x_l w_l}{\sum_l w_l} Cx=∑iWi∑iCi,xWi=∑i∑jwi,j∑i∑jwi,j∑jxi,jwi,j∑jwi,j=∑i∑jwi,j∑i∑jxi,jwi,j=∑lwl∑lxlwl
- C y = ∑ i C i , y A i ∑ i A i C_y = \frac{\sum_i C_{i,y}A_i}{\sum_i A_i} Cy=∑iAi∑iCi,yAi
除了直接基于定义的形式进行计算,还可以利用图像的 p + q p+q p+q阶矩(空间矩/几何矩/原点矩) m p q m_{pq} mpq和中心矩 μ p q \mu_{pq} μpq来定义。
对于一幅二维连续图像, f ( x , y ) ≥ 0 f(x, y) \ge 0 f(x,y)≥0,两个矩的定义为:
- m p q = ∫ − ∞ ∞ ∫ − ∞ ∞ x p y q f ( x , y ) d x d y m_{pq} = \int^{\infty}_{-\infty}\int^{\infty}_{-\infty} x^p y^q f(x,y) dxdy mpq=∫−∞∞∫−∞∞xpyqf(x,y)dxdy
- μ p q = ∫ − ∞ ∞ ∫ − ∞ ∞ ( x − x c ) p ( y − y c ) q f ( x , y ) d x d y \mu_{pq} = \int^{\infty}_{-\infty}\int^{\infty}_{-\infty} (x-x_c)^p (y-y_c)^q f(x,y) dxdy μpq=∫−∞∞∫−∞∞(x−xc)p(y−yc)qf(x,y)dxdy
这里的 ( x c , y c ) (x_c, y_c) (xc,yc)即为质心坐标:
- x c = m 10 m 00 = ∫ − ∞ ∞ ∫ − ∞ ∞ x f ( x , y ) d x d y ∫ − ∞ ∞ ∫ − ∞ ∞ f ( x , y ) d x d y x_c = \frac{m_{10}}{m_{00}} = \frac{\int^{\infty}_{-\infty}\int^{\infty}_{-\infty} x f(x,y) dxdy}{\int^{\infty}_{-\infty}\int^{\infty}_{-\infty} f(x,y) dxdy} xc=m00m10=∫−∞∞∫−∞∞f(x,y)dxdy∫−∞∞∫−∞∞xf(x,y)dxdy
- y c = m 01 m 00 = ∫ − ∞ ∞ ∫ − ∞ ∞ y f ( x , y ) d x d y ∫ − ∞ ∞ ∫ − ∞ ∞ f ( x , y ) d x d y y_c = \frac{m_{01}}{m_{00}} = \frac{\int^{\infty}_{-\infty}\int^{\infty}_{-\infty} y f(x,y) dxdy}{\int^{\infty}_{-\infty}\int^{\infty}_{-\infty} f(x,y) dxdy} yc=m00m01=∫−∞∞∫−∞∞f(x,y)dxdy∫−∞∞∫−∞∞yf(x,y)dxdy
对于离散情形,可以定义为:
- m p q = ∑ i = 1 N ∑ j = 1 N x i p y j q f ( x i , y j ) m_{pq} = \sum^{N}_{i=1}\sum^{N}_{j=1} x_i^p y_j^q f(x_i,y_j) mpq=∑i=1N∑j=1Nxipyjqf(xi,yj)
- μ p q = ∑ i = 1 N ∑ j = 1 N ( x i − x c ) p ( y j − y c ) q f ( x i , y j ) \mu_{pq} = \sum^{N}_{i=1}\sum^{N}_{j=1} (x_i - x_c)^p (y_j-y_c)^q f(x_i,y_j) μpq=∑i=1N∑j=1N(xi−xc)p(yj−yc)qf(xi,yj)
对应的执行可以计算为:
- x c = m 10 m 00 = ∑ i = 1 N ∑ j = 1 N x i p f ( x i , y j ) ∑ i = 1 N ∑ j = 1 N f ( x i , y j ) x_c = \frac{m_{10}}{m_{00}} = \frac{\sum^{N}_{i=1}\sum^{N}_{j=1} x_i^p f(x_i,y_j)}{\sum^{N}_{i=1}\sum^{N}_{j=1} f(x_i,y_j)} xc=m00m10=∑i=1N∑j=1Nf(xi,yj)∑i=1N∑j=1Nxipf(xi,yj)
- y c = m 01 m 00 = ∑ i = 1 N ∑ j = 1 N y j q f ( x i , y j ) ∑ i = 1 N ∑ j = 1 N f ( x i , y j ) y_c = \frac{m_{01}}{m_{00}} = \frac{\sum^{N}_{i=1}\sum^{N}_{j=1} y_j^q f(x_i,y_j)}{\sum^{N}_{i=1}\sum^{N}_{j=1} f(x_i,y_j)} yc=m00m01=∑i=1N∑j=1Nf(xi,yj)∑i=1N∑j=1Nyjqf(xi,yj)
编程实现
网上有很多的实现方式,有基于定义的,也有基于矩的形式的,这里找到了几个进行一下简单的分析。
基于定义
scipy.ndimage.center_of_mass
文档可见https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.center_of_mass.html
从实现中我们可以直观看出来,这是直接基于定义实现的:
normalizer = sum(input, labels, index)
grids = numpy.ogrid[[slice(0, i) for i in input.shape]]
results = [sum(input * grids[dir].astype(float), labels, index) / normalizer
for dir in range(input.ndim)]
if numpy.isscalar(results[0]):
return tuple(results)
return [tuple(v) for v in numpy.array(results).T]
这里提供了同时对多个不同区域的质心计算的支持。
质心计算时,首先使用numpy.ogrid
构造了坐标系网格,之后针对不同的坐标轴遍历,分别对图形区域内的坐标使用原始数据加权求和并归一化。
numpy.argwhere
对于特殊情况,即我们针对二值图计算质心时,可以考虑使用这一方法。
考虑到此时质心的计算实际上仅仅是图形内部坐标的平均,所以可以直接利用argwhere
获得图像区域的像素坐标,对其平均即可。
# https://stackoverflow.com/a/38933601
np.argwhere(x).mean(0)
直接计算
area_object = np.sum(matrix)
row_ids = np.arange(h)
col_ids = np.arange(w)
x = np.round(np.sum(np.sum(matrix, axis=0) * col_ids) / area_object)
y = np.round(np.sum(np.sum(matrix, axis=1) * row_ids) / area_object)
但是这里的代码存在溢出的风险。
numpy不同于python自身的数据表示形式,本身是存在类型限制的,尤其是整型数组容易出现溢出问题:
这里如果对于特别大的图像进行计算,会出现与前两种方式明显不同的结果:
In [53]: def print_centroid(x, h, w):
...: print(np.sum(np.sum(x, axis=1) * np.arange(h)) / np.count_nonzero(x), np.sum(np.sum(x
...: , axis=0) * np.arange(w)) / np.count_nonzero(x))
...: print(center_of_mass(x))
...: print(np.argwhere(x).mean(0))
In [59]: print_centroid(np.random.random((1024, 1024)) > 0.7, 1024, 1024)
511.3934109266679 511.69975831584304
(511.3934109266679, 511.69975831584304)
[511.39341093 511.69975832]
In [60]: print_centroid(np.random.random((2*1024, 2*1024)) > 0.7, 2*1024, 2*1024)
1023.3879048154441 1022.7496662402068
(1023.3879048154441, 1022.7496662402068)
[1023.38790482 1022.74966624]
In [61]: print_centroid(np.random.random((3*1024, 3*1024)) > 0.7, 3*1024, 3*1024)
18.466653739911568 18.801060064556072 <--此时输出已经开始异常
(1535.4062819085118, 1535.7406882331563)
[1535.40628191 1535.74068823]
In [62]: print_centroid(np.random.random((4*1024, 4*1024)) > 0.7, 4*1024, 4*1024)
341.39504553161055 341.632512348338
(2047.3403782063924, 2047.5778450231198)
[2047.34037821 2047.57784502]
In [63]: print_centroid(np.random.random((8*1024, 8*1024)) > 0.7, 8*1024, 8*1024)
42.42419332351165 42.500630084653515
(4095.6871022111004, 4095.7635389722423)
[4095.68710221 4095.76353897]
参考前面scipy的实现,我们可以这样修改:
In [73]: def print_centroid(x, h, w):
...: print(np.sum(np.sum(x, axis=1) * np.arange(h).astype(float)) / np.count_nonzero(x), n
...: p.sum(np.sum(x, axis=0) * np.arange(w).astype(float)) / np.count_nonzero(x))
...: print(center_of_mass(x))
...: print(np.argwhere(x).mean(0))
In [74]: print_centroid(np.random.random((8*1024, 8*1024)) > 0.7, 8*1024, 8*1024)
4095.778456882133 4095.5295789226266
(4095.778456882133, 4095.5295789226266)
[4095.77845688 4095.52957892]
In [75]: print_centroid(np.random.random((3*1024, 3*1024)) > 0.7, 3*1024, 3*1024)
1535.384796905072 1535.7708201396363
(1535.384796905072, 1535.7708201396363)
[1535.38479691 1535.77082014]
此时便不再容易溢出了。
基于矩的方式
cv2.moments
关于不同矩的介绍可见中的介绍。
从这篇文章https://www.geeksforgeeks.org/python-opencv-find-center-of-contour/中我们可以注意到,opencv提供了计算图像矩的功能。
该函数有两种使用方式:
核心代码如下:
# 直接处理图像
m = cv2.moments(image)
cx = m['m10']/m['m00']
cy = m['m01']/m['m00']
# 提取轮廓
contours, hierarchies = cv.findContours(thresh, cv.RETR_LIST, cv.CHAIN_APPROX_SIMPLE)
# 根据轮廓计算不同轮廓对应的质心
for i in contours:
m = cv.moments(i)
cx = m['m10'] / m['m00']
cy = m['m01'] / m['m00']
使用该函数计算最好使用真实图像,结果和前三种是一致的。
In [41]: from skimage import data, img_as_float
...: img = img_as_float(data.camera()) > 0.5
In [44]: def print_centroid(x, h, w):
...: print(np.sum(np.sum(x, axis=1) * np.arange(h).astype(float)) / np.count_nonzero(x), n
...: p.sum(np.sum(x, axis=0) * np.arange(w).astype(float)) / np.count_nonzero(x))
...: print(center_of_mass(x))
...: print(np.argwhere(x).mean(0))
...: mu = cv2.moments(x.astype(np.uint8))
...: print(mu['m01'] / mu['m00'], mu['m10'] / mu['m00'])
...:
In [45]: np.count_nonzero(img)
Out[45]: 168559
In [46]: print_centroid(img, 512, 512)
231.7689117756987 309.7281604660683
(231.7689117756987, 309.7281604660683)
[231.76891178 309.72816047]
231.7689117756987 309.7281604660683
参考
更多推荐
所有评论(0)