如何计算质心

原始文档：https://www.yuque.com/lart/idh721/gpbigm

概念

质心，即质量中心的简称。质点系的质心是质点系质量分布的平均位置。指物质系统上被认为质量集中于此的一个假想点，与重心不同的是质心不一定要在有重力场的系统中，值得注意的是除非重力场是均匀的，否则同一物质系统的质心与重心通常不在同一假想点上。

计算

质心坐标等于所有点关于每个坐标的以质量为权重的加权平均值。一般主要在二维空间讨论，尤其是图像数据，但是这里直接按照更一般的形式进行定义。首先对于任意$n$维空间中的连续形式的子集$P$的质心可以定义为：

$$C=\frac{\int pg(p)dp}{\int g(p)dp}$$

其中：

• $p\in {\mathbb{R}}^n$表示该子集中的一点；
• $g$表示该子集的特征函数（indicator function or a characteristic function）。一个比较实际的场景是，它可以用来表示各个位置对应的质量。

也可以看到，这里的分母是对这个集合的一个度量，因而如果度量为0，那么就不可以被计算质心。

而对于第$k$个坐标$C_k$的计算，可以通过如下形式：

$$C_k=\frac{\int z{S}_k(z)dz}{\int g(x)dz}$$

这里$S_k(z)$表示的是对应的$P$与由$p_k=z$定义的超平面（hyperplane）的交集的度量，在这个超平面上，会涉及到其他所有的坐标轴。

对于一个平面图，即二维情形，上面的式子可以用于求取$C_x$和$C_y$：

• $C_x=\frac{\int x{S}_x(x)dx}{A}$
• $C_y=\frac{\int y{S}_y(y)dy}{A}$

即对单一轴向上的坐标积分，每个坐标对应乘以一个与之关联的量$S$（该量会涉及到另一个轴），$A=\int S_x(x)dx=\int S_y(y)dy$表示图形的面积。一般情况下，我们可以简单的理解为这是过点$(x,0)$或是$(0,y)$的垂线与图像区域的相交构成的线段的长度。

对于更实际的离散且有限点集的情形下，前面二维的形式可以转化为如下形式：

• $C_x=\frac{{\sum }_i{w}_i{x}_i}{{\sum }_i{w}_i}$
• $C_y=\frac{{\sum }_i{w}_i{y}_i}{{\sum }_i{w}_i}$

这里要注意，公式中表示各个点的方式与前面直接基于坐标值的方式有所不同，而是通过一个额外的点索引$i$来对不同的点进行编码排序，从而构建了公式。

对于不同的点有着不同的对应权重，我们可以理解为质量或者其他的度量形式，所以对应相同的点序号$i$，其各个轴向上的坐标权重也是一样的$w_i$，且$W={\sum }_i{w}_i$可以表示为图像对应的整体质量或者其他的度量。如果各点权重均为1，则这里的$W$则实际上便是点的数量，对于数字图像而言，就是图形面积了。

更一般的，这里的点$i$实际上还可以替换为有着面积（或者说离散点的数量）$A_i$的区域$P_i$。计算过程中将其质心作为这里的点。

• $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}$
• $C_y=\frac{{\sum }_i{C}_{i,y}{A}_i}{{\sum }_i{A}_i}$

除了直接基于定义的形式进行计算，还可以利用图像的$p+q$阶矩（空间矩/几何矩/原点矩）$m_{pq}$和中心矩${\mu }_{pq}$来定义。

对于一幅二维连续图像，$f(x,y)\ge 0$，两个矩的定义为：

• $m_{pq}={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{x}^p{y}^{q}f(x,y)dxdy$
• ${\mu }_{pq}={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }(x-x_c{)}^p(y-y_c{)}^{q}f(x,y)dxdy$

这里的$(x_c,y_c)$即为质心坐标：

• $x_c=\frac{m_{10}}{m_{00}}=\frac{{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }xf(x,y)dxdy}{{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)dxdy}$
• $y_c=\frac{m_{01}}{m_{00}}=\frac{{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }yf(x,y)dxdy}{{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)dxdy}$

对于离散情形，可以定义为：

• $m_{pq}={\sum }_{i=1}^N{\sum }_{j=1}^N{x}_i^p{y}_j^{q}f(x_i,y_j)$
• ${\mu }_{pq}={\sum }_{i=1}^N{\sum }_{j=1}^N(x_i-x_c{)}^p(y_j-y_c{)}^{q}f(x_i,y_j)$

对应的执行可以计算为：

• $x_c=\frac{m_{10}}{m_{00}}=\frac{{\sum }_{i=1}^N{\sum }_{j=1}^N{x}_i^{p}f(x_i,y_j)}{{\sum }_{i=1}^N{\sum }_{j=1}^{N}f(x_i,y_j)}$
• $y_c=\frac{m_{01}}{m_{00}}=\frac{{\sum }_{i=1}^N{\sum }_{j=1}^N{y}_j^{q}f(x_i,y_j)}{{\sum }_{i=1}^N{\sum }_{j=1}^{N}f(x_i,y_j)}$

编程实现

网上有很多的实现方式，有基于定义的，也有基于矩的形式的，这里找到了几个进行一下简单的分析。

基于定义

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)

直接计算

https://github.com/lartpang/PySODMetrics/blob/4aa253a59aff71507f92daf2dffe539c5c97ce46/py_sod_metrics/sod_metrics.py#L277-L282

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

参考

• https://en.wikipedia.org/wiki/Centroid
• https://baike.baidu.com/item/%E8%B4%A8%E5%BF%83