| 如何计算质心 | 您所在的位置：网站首页 › 数学二质心和形心公式面积 › 如何计算质心 |
|
如何计算质心
原始文档：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) 直接计算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 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 |
| CopyRight 2018-2019 实验室设备网 版权所有 |