计算质心：原理、方法与编程实例详解

如何计算质心

原始文档：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​=∑i​wi​∑i​wi​xi​​

C

y

=

∑

i

w

i

y

i

∑

i

w

i

C_y = \frac{\sum_i w_i y_i}{\sum_i w_i}

Cy​=∑i​wi​∑i​wi​yi​​

这里要注意，公式中表示各个点的方式与前面直接基于坐标值的方式有所不同，而是通过一个额外的点索引

i

i

i来对不同的点进行编码排序，从而构建了公式。

对于不同的点有着不同的对应权重，我们可以理解为质量或者其他的度量形式，所以对应相同的点序号

i

i

i，其各个轴向上的坐标权重也是一样的

w

i

w_i

wi​，且

W

=

∑

i

w

i

W = \sum_i w_i

W=∑i​wi​可以表示为图像对应的整体质量或者其他的度量。如果各点权重均为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​=∑i​Wi​∑i​Ci,x​Wi​​=∑i​∑j​wi,j​∑i​∑j​wi,j​∑j​xi,j​wi,j​​∑j​wi,j​​=∑i​∑j​wi,j​∑i​∑j​xi,j​wi,j​​=∑l​wl​∑l​xl​wl​​

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​=∑i​Ai​∑i​Ci,y​Ai​​

除了直接基于定义的形式进行计算，还可以利用图像的

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​=m00​m10​​=∫−∞∞​∫−∞∞​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​=m00​m01​​=∫−∞∞​∫−∞∞​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=1N​xip​yjq​f(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​=m00​m10​​=∑i=1N​∑j=1N​f(xi​,yj​)∑i=1N​∑j=1N​xip​f(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​=m00​m01​​=∑i=1N​∑j=1N​f(xi​,yj​)∑i=1N​∑j=1N​yjq​f(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)

直接计算

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/Centroidhttps://baike.baidu.com/item/%E8%B4%A8%E5%BF%83