• 图像处理之Bolb分析（一）

    文章目录

    前言

    Bolb分析，即连通域分析。开箱即用的开源库有OpenCV中的SimpleBlobDetector类、findContours，以及skimage.measue的regionprops等。这里介绍skimage.measue的regionprops。

    skimage中的regionprops

    regionprops

    以skimage 0.19.1为例：

    源码位置：https://github.com/scikit-image/scikit-image/blob/main/skimage/measure/_regionprops.py

    regionprops函数

    def regionprops(label_image, intensity_image=None, cache=True,
                coordinates=None, *, extra_properties=None):
    r"""测量标记图像连通域的属性.
    Parameters
    ----------
    label_image : (M, N[, P]) ndarray
        输入标签图像. 0标签忽略
    intensity_image : (M, N[, P][, C]) ndarray, optional
        和标签图像大小相同的强度图像, 可选参数。
        可选的多通道数据的额外维度。目前,这个额外的通道尺寸，如果存在，必须是最后一个轴。默认是None
    cache : bool, optional
        确定是否缓存计算的属性。缓存属性的计算速度要快得多，但内存消耗会增加。
    coordinates : DEPRECATED
        This argument is deprecated and will be removed in a future version
        of scikit-image.
    extra_properties : Iterable of callables
        添加skimage中没有的额外属性计算函数。属性名和函数名函数名一致，dtype是通过在一个小样本上调用函数来推断的。
        如果额外属性的名称与现有属性的名称冲突这个额外的属性将不可见，并使用UserWarning警告。
        一个属性计算函数必须接受一个区域掩码作为它的第一个参数。如果属性需要一个强度图像，它必须接受强度图像作为第二个参数。

    Returns
    -------
    properties : list of RegionProperties
       连通域属性列表，可用的属性在下一节
    """

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25

    属性及数据类型

    {#********几何特征********#
    'area': int,	# 连通域面积
    'area_bbox': int,	# 外接矩形面积
    'area_convex': int,	# 最小凸包面积
    'area_filled': int,	# 孔填充后面积
    'axis_major_length': float,	# 长轴 和连通域有相同标准第二中心矩的椭圆长轴
    'axis_minor_length': float,	# 短轴
    'bbox': int,	# 边界框，即外接矩形坐标 (min_row, min_col, max_row, max_col)
    'centroid': float,	# 质心坐标	(row, col)
    'centroid_local': float,	# 质心-局部坐标
    'centroid_weighted': float,	# 使用强度图像加权的质心坐标
    'centroid_weighted_local': float,# 使用强度图像加权的质心-局部坐标
    'eccentricity': float,	# 偏心率 和连通域有相同二阶中心矩的椭圆的离心率。公式：焦距c/长轴a [0, 1)，0时为圆。
    'equivalent_diameter_area': float,	# 等效直径 和连通域面积相同的圆的直径
    'euler_number': int,	# 欧拉数	非零像素集合的欧拉特征。公式:连通域数目 - 洞数目。3D里面：连通域数目*洞数目-隧道数目
    'extent': float,	# 范围 连通域面积/外接矩形面积 area / (rows * cols)
    'feret_diameter_max': float,	# 最大Feret直径，由计算通过find_contours获取的连通域凸包轮廓点之间的最长距离
    'orientation': float,	# 与连通域有相同二阶矩的椭圆的第0轴(rows)与主轴之间的夹角，从' -pi/2 '到' pi/2 '，逆时针方向。
    'perimeter': float,	# 连通域周长，其轮廓近似为通过边界轮廓像素中心的线，使用4连通分析。
    'perimeter_crofton': float,	# 用四方向的克罗夫顿公式逼近物体的周长。
    'solidity': float,	# 连通域像素与凸包像素之比
 #********连通域********#
    'slice': object,	# 从源图像中提取物体的切片。
    'coords': object,	# 连通域的坐标 (N, 2) ndarray
    'image': object,	# 与边界框大小相同的二值区域图像切片
    'image_convex': object,	# 与边界框大小相同的二值凸包图像。
    'image_filled': object,	# 与边界框大小相同的填充孔的二值图像
 #********灰度特征********#
    'intensity_max': float,	# 连通域强度最大值（需提供二值图像对应的强度图像）
    'intensity_mean': float,	# 连通域强度平均值
    'intensity_min': float,	# 连通域强度最小值 
    'label': int,	# 标记图像中连通域的label
 #********矩特征********#
    'image_intensity': object,	# 边界框内的强度图像
    'inertia_tensor': float,	# 绕质量旋转的连通域区域惯性张量。
    'inertia_tensor_eigvals': float,	# 惯性张量特征值，按降序排列
    'moments': float,	# 三阶空间矩 m_ij = sum{ array(row, col) * row^i * col^j }，其中sum在连通域行列上求和 (3, 3) ndarray
    'moments_central': float,	# 三阶的中心矩(平移不变)， mu_ij = sum{ array(row, col) * (row - row_c)^i * (col - col_c)^j }，' row_c '和' col_c '是连通域的坐标重心。
    'moments_hu': float,	# hu矩(平移、缩放和旋转不变量)。
    'moments_normalized': float,	# 三阶归一化矩(平移和缩放不变) nu_ij = mu_ij / m_00^[(i+j)/2 + 1], 其中m_00是零阶空间矩
    'moments_weighted': float,	# 强度图像的三阶空间矩，wm_ij = sum{ array(row, col) * row^i * col^j }，其中sum在连通域行列上求和
    'moments_weighted_central': float,	# 强度图像的三阶中心空间矩， wmu_ij = sum{ array(row, col) * (row - row_c)^i * (col - col_c)^j }，其中sum在连通域行列上求和， ' row_c '和' col_c '是连通域加权的坐标重心。
    'moments_weighted_hu': float,	# 强度图像的三阶hu矩(平移、缩放和旋转不变量) 
    'moments_weighted_normalized': float,# 强度图像的三阶归一化矩(平移和缩放不变) wnu_ij = wmu_ij / wm_00^[(i+j)/2 + 1]，：“wm_00”是0零阶空间矩(强度权重面积) 
}

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35
    • 36
    • 37
    • 38
    • 39
    • 40
    • 41
    • 42
    • 43
    • 44
    • 45

    References
    [1] Wilhelm Burger, Mark Burge. Principles of Digital Image Processing:
    Core Algorithms. Springer-Verlag, London, 2009.
    [2] B. Jähne. Digital Image Processing. Springer-Verlag,
    Berlin-Heidelberg, 6. edition, 2005.
    [3] T. H. Reiss. Recognizing Planar Objects Using Invariant Image
    Features, from Lecture notes in computer science, p. 676. Springer,
    Berlin, 1993.
    [4] https://en.wikipedia.org/wiki/Image_moment
    [5] W. Pabst, E. Gregorová. Characterization of particles and particle
    systems, pp. 27-28. ICT Prague, 2007.
    https://old.vscht.cz/sil/keramika/Characterization_of_particles/CPPS%20_English%20version_.pdf

    源码

    这里介绍一些重要的源码，数据类型检测，参数检查之类的这里不再赘述。

    regionprops函数

    def regionprops(label_image, intensity_image=None, cache=True,
                coordinates=None, *, extra_properties=None):
    if label_image.ndim not in (2, 3):	# 函数只支持2维和3维图片
        raise TypeError('Only 2-D and 3-D images supported.')

    if not np.issubdtype(label_image.dtype, np.integer):	# 定义支持的数据类型 即integer的子类型，比如np.int16 np.uint8
        if np.issubdtype(label_image.dtype, bool):
            raise TypeError(
                    'Non-integer image types are ambiguous: '
                    'use skimage.measure.label to label the connected'
                    'components of label_image,'
                    'or label_image.astype(np.uint8) to interpret'
                    'the True values as a single label.')
        else:
            raise TypeError(
                    'Non-integer label_image types are ambiguous')

    if coordinates is not None:	# 坐标系 rc坐标系即row-column坐标系
        if coordinates == 'rc':
            msg = ('The coordinates keyword argument to skimage.measure.'
                   'regionprops is deprecated. All features are now computed '
                   'in rc (row-column) coordinates. Please remove '
                   '`coordinates="rc"` from all calls to regionprops before '
                   'updating scikit-image.')
            warn(msg, stacklevel=2, category=FutureWarning)
        else:
            msg = ('Values other than "rc" for the "coordinates" argument '
                   'to skimage.measure.regionprops are no longer supported. '
                   'You should update your code to use "rc" coordinates and '
                   'stop using the "coordinates" argument, or use skimage '
                   'version 0.15.x or earlier.')
            raise ValueError(msg)

    regions = []

    objects = ndi.find_objects(label_image)	# 返回每一个连通域的切片 忽略0标签
    for i, sl in enumerate(objects):
        if sl is None:
            continue

        label = i + 1

        props = RegionProperties(sl, label, label_image, intensity_image,
                                 cache, extra_properties=extra_properties)	# 构造RegionProperties类
        regions.append(props)	# regions里面存的是一个个类

    return regions


    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35
    • 36
    • 37
    • 38
    • 39
    • 40
    • 41
    • 42
    • 43
    • 44
    • 45
    • 46
    • 47
    • 48

    其实例化了RegionProperties类

    class RegionProperties:
    """Please refer to `skimage.measure.regionprops` for more information
    on the available region properties.
    """
	# 初始化
    def __init__(self, slice, label, label_image, intensity_image,
                 cache_active, *, extra_properties=None):

        if intensity_image is not None:	# 判断是否有强度图像，用于计算某些需要强度图像的属性，比如intensity_max等
            ndim = label_image.ndim
            if not (
                    intensity_image.shape[:ndim] == label_image.shape
                    and intensity_image.ndim in [ndim, ndim + 1]
                    ):
                raise ValueError('Label and intensity image shapes must match,'
                                 ' except for channel (last) axis.')
            multichannel = label_image.shape < intensity_image.shape
        else:
            multichannel = False
		# 初始化参数
        self.label = label

        self._slice = slice
        self.slice = slice
        self._label_image = label_image
        self._intensity_image = intensity_image

        self._cache_active = cache_active
        self._cache = {}
        self._ndim = label_image.ndim
        self._multichannel = multichannel
        self._spatial_axes = tuple(range(self._ndim))
		
        self._extra_properties = {}
        if extra_properties is not None:	# 如果有自定义的属性函数，在这里获取到其函数名
            for func in extra_properties:
                name = func.__name__
                if hasattr(self, name):	# 如果自定义属性函数名字与内置属性函数名重叠，则忽略并警告
                    msg = (
                        f"Extra property '{name}' is shadowed by existing "
                        f"property and will be inaccessible. Consider "
                        f"renaming it."
                    )
                    warn(msg)
            self._extra_properties = {
                func.__name__: func for func in extra_properties
            }
    # 重载__getattr__函数，在里面做一些判断  
    def __getattr__(self, attr):
        if attr in self._extra_properties:
            func = self._extra_properties[attr]
            n_args = _infer_number_of_required_args(func)
            # 确定函数是否需要强度图像
            if n_args == 2:
                if self._intensity_image is not None:
                    if self._multichannel:
                        multichannel_list = [func(self.image,
                                                  self.image_intensity[..., i])
                                             for i in range(
                            self.image_intensity.shape[-1])]
                        return np.stack(multichannel_list, axis=-1)
                    else:
                        return func(self.image, self.image_intensity)
                else:
                    raise AttributeError(
                        f'intensity image required to calculate {attr}'
                    )
            elif n_args == 1:
                return func(self.image)
            else:
                raise AttributeError(
                    f'Custom regionprop function\'s number of arguments must '
                    f'be 1 or 2, but {attr} takes {n_args} arguments.'
                )
        elif attr in PROPS and attr.lower() == attr:	# 为了向后兼容，有些过时的名字
            # 检索已弃用的属性(不包括旧的CamelCase属性) 
            return getattr(self, PROPS[attr])	
        else:
            raise AttributeError(
                f"'{type(self)}' object has no attribute '{attr}'"
            )
	# 重载__setattr__函数，在里面做一些判断
    def __setattr__(self, name, value):
        if name in PROPS:
            super().__setattr__(PROPS[name], value)
        else:
            super().__setattr__(name, value)

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35
    • 36
    • 37
    • 38
    • 39
    • 40
    • 41
    • 42
    • 43
    • 44
    • 45
    • 46
    • 47
    • 48
    • 49
    • 50
    • 51
    • 52
    • 53
    • 54
    • 55
    • 56
    • 57
    • 58
    • 59
    • 60
    • 61
    • 62
    • 63
    • 64
    • 65
    • 66
    • 67
    • 68
    • 69
    • 70
    • 71
    • 72
    • 73
    • 74
    • 75
    • 76
    • 77
    • 78
    • 79
    • 80
    • 81
    • 82
    • 83
    • 84
    • 85
    • 86
    • 87

    area

    @property
@_cached
def area(self):
    return np.sum(self.image)

    • 1
    • 2
    • 3
    • 4

    使用np.sum()统计连通域面积

    关于装饰器@property和@_cached：

    @property装饰器就是负责把一个方法变成属性调用的。是python内置装饰器

    参考：https://blog.csdn.net/qq_42350970/article/details/84943534

    @_cached装饰器是_reginoprops.py里面定义的，用于判断计算的属性是否缓存。

    后面所有的属性函数都有@property装饰器，用于将该方法变成属性调用，一些可能会复用的属性，有@_cached装饰器，用于控制缓存该属性值，加快计算速度。

    bbox

    @property
def bbox(self):
    """
    Returns
    -------
    A tuple of the bounding box's start coordinates for each dimension,
    followed by the end coordinates for each dimension
    """
    return tuple([self.slice[i].start for i in range(self._ndim)] +
                 [self.slice[i].stop for i in range(self._ndim)])

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10

    这里的连通域slice切片存储格式为：

    ((lf_top_r, rg_bot_r, None), (lf_top_c, lf_top_c, None))

    • 1

    所以边界框直接取出来就行。

    coords、centroid

    @property
def centroid(self):
    return tuple(self.coords.mean(axis=0))

    • 1
    • 2
    • 3

    质心即连通域坐标取平均值即可，其中获取连通域坐标的函数coords()定义如下：

    @property
def coords(self):
    indices = np.nonzero(self.image)
    return np.vstack([indices[i] + self.slice[i].start
                          for i in range(self._ndim)]).T

    • 1
    • 2
    • 3
    • 4
    • 5

    这里np.nonzero()获取的坐标是局部坐标，所以需要加上连通域左上角的坐标得到的才是全局坐标。

    image_convex、area_convex

    @property
@_cached
def image_convex(self):
    from ..morphology.convex_hull import convex_hull_image
    return convex_hull_image(self.image)

    • 1
    • 2
    • 3
    • 4
    • 5

    凸包图像计算，使用了skimage.morphology.convex_hull_image函数，具体算法实现可参考：

    https://github.com/scikit-image/scikit-image/blob/main/skimage/morphology/convex_hull.py

    @property
@_cached
def area_convex(self):
    return np.sum(self.image_convex)

    • 1
    • 2
    • 3
    • 4

    area_convex凸包图像面积使用np.sum()计算

    eccentricity

    @property
@only2d
def eccentricity(self):
    l1, l2 = self.inertia_tensor_eigvals
    if l1 == 0:
        return 0
    return sqrt(1 - l2 / l1)

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7

    偏心率 和连通域有相同二阶中心矩的椭圆的离心率。公式：焦距c/长轴a [0, 1)，0时为圆。

    其使用到了inertia_tensor_eigvals()函数：

    @property
@_cached
def inertia_tensor_eigvals(self):
    return _moments.inertia_tensor_eigvals(self.image,
                                           T=self.inertia_tensor)

    • 1
    • 2
    • 3
    • 4
    • 5

    其调用了_moments.inertia_tensor_eigvals()函数：

    vb6的156种图像处理。
    rar 5星 超过95%的资源 92KB
    下载
    @property
@_cached
def inertia_tensor(self):
    mu = self.moments_central
    return _moments.inertia_tensor(self.image, mu)

    • 1
    • 2
    • 3
    • 4
    • 5

    其调用了moments_central和_moments.inertia_tensor()

    而moments_central又调用了_moments.moments_central(self.image.astype(np.uint8),
    self.centroid_local, order=3)

    self.centroid_local调用了moments()方法，moments()调用了 _moments.moments

    所以离心率的计算先后顺序为：

     _moments.moments # 矩
self.centroid_local	# 局部质心
_moments.moments_central	# 中心矩
_moments.inertia_tensor	# 惯性张量
_moments.inertia_tensor_eigvals	# 惯性张量特征值
# 然后根据惯性张量特征值计算离心率。

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6

    关于中心矩、惯性张量、惯性张量特征值算法实现可参考：

    https://github.com/scikit-image/scikit-image/blob/main/skimage/measure/_moments.py

    • moments 矩
    def moments(image, order=3):
    """计算指定阶数的图像原始矩.
    以下属性可以通过图像原始矩计算得到:
     * Area as: ``M[0, 0]``.
     * Centroid as: {``M[1, 0] / M[0, 0]``, ``M[0, 1] / M[0, 0]``}.
    注意：图像原始矩不是平移、缩放、旋转不变的。
    Parameters
    ----------
    image : nD double or uint8 array
        Rasterized shape as image.
    order : int, optional
        Maximum order of moments. Default is 3.
    Returns
    -------
    m : (``order + 1``, ``order + 1``) array
        Raw image moments.
    References
    ----------
    .. [1] Wilhelm Burger, Mark Burge. Principles of Digital Image Processing:
           Core Algorithms. Springer-Verlag, London, 2009.
    .. [2] B. Jähne. Digital Image Processing. Springer-Verlag,
           Berlin-Heidelberg, 6. edition, 2005.
    .. [3] T. H. Reiss. Recognizing Planar Objects Using Invariant Image
           Features, from Lecture notes in computer science, p. 676. Springer,
           Berlin, 1993.
    .. [4] https://en.wikipedia.org/wiki/Image_moment
    Examples
    --------
    >>> image = np.zeros((20, 20), dtype=np.double)
    >>> image[13:17, 13:17] = 1
    >>> M = moments(image)
    >>> centroid = (M[1, 0] / M[0, 0], M[0, 1] / M[0, 0])
    >>> centroid
    (14.5, 14.5)
    """
    return moments_central(image, (0,) * image.ndim, order=order)

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35
    • 36

    其调用了moments_central

    对于二维连续函数f(x,y)，(p + q)阶矩(有时称为原始矩)定义为：
    M p q = ∫ − ∞ ∞ ∫ − ∞ ∞ x p y q f ( x , y ) d x d y M_{p q}=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} x^{p} y^{q} f(x, y) d x d y Mpq​=∫−∞∞​∫−∞∞​xpyqf(x,y)dxdy
    对于p,q = 0,1,2，…将其适应于像素强度为I(x,y)的标量(灰度)图像，原始图像矩Mij计算为：
    M i j = ∑ x ∑ y x i y j I ( x , y ) M_{i j}=\sum_{x} \sum_{y} x^{i} y^{j} I(x, y) Mij​=x∑​y∑​xiyjI(x,y)
    在某些情况下，这可以通过将图像考虑为概率密度函数来计算，即，通过除以上面的
    ∑ x ∑ y I ( x , y ) \sum_{x} \sum_{y} I(x, y) x∑​y∑​I(x,y)
    一个唯一性定理(Hu[1962])指出，如果f(x,y)是分段连续的，且仅在xy平面的有限部分上具有非零值，则所有阶矩都存在，且矩序列(Mpq)是由f(x,y)唯一确定的。相反，(Mpq)唯一确定f(x,y)。在实际应用中，图像是用几个低阶矩的函数来概括的。

    通过原始矩可计算的简单图像属性有：

    1. 面积 (用于二值图像)或灰度总和(用于灰度图像)： M 00 M_{00} M00​

    2. 中心： { x ˉ , y ˉ } = { M 10 M 00 , M 01 M 00 } \{\bar{x}, \bar{y}\}=\left\{\frac{M_{10}}{M_{00}}, \frac{M_{01}}{M_{00}}\right\} {xˉ,yˉ​}={M00​M10​​,M00​M01​​}

    • moments_central
    def moments_central(image, center=None, order=3, **kwargs):
    """Calculate all central image moments up to a certain order.
    The center coordinates (cr, cc) can be calculated from the raw moments as:
    {``M[1, 0] / M[0, 0]``, ``M[0, 1] / M[0, 0]``}.
    Note that central moments are translation invariant but not scale and
    rotation invariant.
    Parameters
    ----------
    image : nD double or uint8 array
        Rasterized shape as image.
    center : tuple of float, optional
        Coordinates of the image centroid. This will be computed if it
        is not provided.
    order : int, optional
        The maximum order of moments computed.
    Returns
    -------
    mu : (``order + 1``, ``order + 1``) array
        Central image moments.
    References
    ----------
    .. [1] Wilhelm Burger, Mark Burge. Principles of Digital Image Processing:
           Core Algorithms. Springer-Verlag, London, 2009.
    .. [2] B. Jähne. Digital Image Processing. Springer-Verlag,
           Berlin-Heidelberg, 6. edition, 2005.
    .. [3] T. H. Reiss. Recognizing Planar Objects Using Invariant Image
           Features, from Lecture notes in computer science, p. 676. Springer,
           Berlin, 1993.
    .. [4] https://en.wikipedia.org/wiki/Image_moment
    Examples
    --------
    >>> image = np.zeros((20, 20), dtype=np.double)
    >>> image[13:17, 13:17] = 1
    >>> M = moments(image)
    >>> centroid = (M[1, 0] / M[0, 0], M[0, 1] / M[0, 0])
    >>> moments_central(image, centroid)
    array([[16.,  0., 20.,  0.],
           [ 0.,  0.,  0.,  0.],
           [20.,  0., 25.,  0.],
           [ 0.,  0.,  0.,  0.]])
    """
    if center is None:	# 如果没有指定中心，使用1阶原始矩计算中心
        center = centroid(image)
    float_dtype = _supported_float_type(image.dtype)
    calc = image.astype(float_dtype, copy=False)
    # 开始计算中心矩
    for dim, dim_length in enumerate(image.shape):
        delta = np.arange(dim_length, dtype=float_dtype) - center[dim]
        powers_of_delta = (
            delta[:, np.newaxis] ** np.arange(order + 1, dtype=float_dtype)
        )
        calc = np.rollaxis(calc, dim, image.ndim)
        calc = np.dot(calc, powers_of_delta)
        calc = np.rollaxis(calc, -1, dim)
    return calc

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35
    • 36
    • 37
    • 38
    • 39
    • 40
    • 41
    • 42
    • 43
    • 44
    • 45
    • 46
    • 47
    • 48
    • 49
    • 50
    • 51
    • 52
    • 53
    • 54
    • 55

    中心矩定义为：
    μ p q = ∫ − ∞ ∞ ∫ − ∞ ∞ ( x − x ˉ ) p ( y − y ˉ ) q f ( x , y ) d x d y \mu_{p q}=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty}(x-\bar{x})^{p}(y-\bar{y})^{q} f(x, y) d x d y μpq​=∫−∞∞​∫−∞∞​(x−xˉ)p(y−yˉ​)qf(x,y)dxdy
    其中， x ˉ = M 10 M 00 \bar{x}=\frac{M_{10}}{M_{00}} xˉ=M00​M10​​， y ˉ = M 01 M 00 \bar{y}=\frac{M_{01}}{M_{00}} yˉ​=M00​M01​​，是质心的分量

    如果 f ( x , y ) f(x, y) f(x,y)是数字图像，上面的公式变为：
    μ p q = ∑ x ∑ y ( x − x ˉ ) p ( y − y ˉ ) q f ( x , y ) \mu_{p q}=\sum_{x} \sum_{y}(x-\bar{x})^{p}(y-\bar{y})^{q} f(x, y) μpq​=x∑​y∑​(x−xˉ)p(y−yˉ​)qf(x,y)
    三阶中心矩为：
    μ 00 = M 00 , μ 01 = 0 , μ 10 = 0 , μ 11 = M 11 − x ˉ M 01 = M 11 − y ˉ M 10 , μ 20 = M 20 − x ˉ M 10 , μ 02 = M 02 − y ˉ M 01 , μ 21 = M 21 − 2 x ˉ M 11 − y ˉ M 20 + 2 x ˉ 2 M 01 , μ 12 = M 12 − 2 y ˉ M 11 − x ˉ M 02 + 2 y ˉ 2 M 10 , μ 30 = M 30 − 3 x ˉ M 20 + 2 x ˉ 2 M 10 , μ 03 = M 03 − 3 y ˉ M 02 + 2 y ˉ 2 M 01 .

    μ00=M00,μ01=0,μ10=0,μ11=M11−x¯M01=M11−y¯M10,μ20=M20−x¯M10,μ02=M02−y¯M01,μ21=M21−2x¯M11−y¯M20+2x¯2M01,μ12=M12−2y¯M11−x¯M02+2y¯2M10,μ30=M30−3x¯M20+2x¯2M10,μ03=M03−3y¯M02+2y¯2M01." role="presentation" style="position: relative;">μ00=M00,μ01=0,μ10=0,μ11=M11−x¯M01=M11−y¯M10,μ20=M20−x¯M10,μ02=M02−y¯M01,μ21=M21−2x¯M11−y¯M20+2x¯2M01,μ12=M12−2y¯M11−x¯M02+2y¯2M10,μ30=M30−3x¯M20+2x¯2M10,μ03=M03−3y¯M02+2y¯2M01.
    ​μ00​=M00​,μ01​=0,μ10​=0,μ11​=M11​−xˉM01​=M11​−yˉ​M10​,μ20​=M20​−xˉM10​,μ02​=M02​−yˉ​M01​,μ21​=M21​−2xˉM11​−yˉ​M20​+2xˉ2M01​,μ12​=M12​−2yˉ​M11​−xˉM02​+2yˉ​2M10​,μ30​=M30​−3xˉM20​+2xˉ2M10​,μ03​=M03​−3yˉ​M02​+2yˉ​2M01​.​
    可以写成：
    μ p q = ∑ m p ∑ n q ( p m ) ( q n ) ( − x ˉ ) ( p − m ) ( − y ˉ ) ( q − n ) M m n \mu_{p q}=\sum_{m}^{p} \sum_{n}^{q}\left(
    pm" role="presentation" style="position: relative;">pm
    \right)\left(
    qn" role="presentation" style="position: relative;">qn
    \right)(-\bar{x})^{(p-m)}(-\bar{y})^{(q-n)} M_{m n}     μpq​=m∑p​n∑q​(pm​)(qn​)(−xˉ)(p−m)(−yˉ​)(q−n)Mmn​
    中心矩是平移不变的。

    利用二阶中心矩构造协方差矩阵，可以得到图像的方向信息。
    μ 20 ′ = μ 20 / μ 00 = M 20 / M 00 − x ˉ 2 μ 02 ′ = μ 02 / μ 00 = M 02 / M 00 − y ˉ 2 μ 11 ′ = μ 11 / μ 00 = M 11 / M 00 − x ˉ y ˉ

    μ20′=μ20/μ00=M20/M00−x¯2μ02′=μ02/μ00=M02/M00−y¯2μ11′=μ11/μ00=M11/M00−x¯y¯" role="presentation" style="position: relative;">μ20′=μ20/μ00=M20/M00−x¯2μ02′=μ02/μ00=M02/M00−y¯2μ11′=μ11/μ00=M11/M00−x¯y¯
    ​μ20′​=μ20​/μ00​=M20​/M00​−xˉ2μ02′​=μ02​/μ00​=M02​/M00​−yˉ​2μ11′​=μ11​/μ00​=M11​/M00​−xˉyˉ​​

    图像 I ( x , y ) I(x, y) I(x,y)的协方差矩阵变为：
    cov ⁡ [ I ( x , y ) ] = [ μ 20 ′ μ 11 ′ μ 11 ′ μ 02 ′ ] \operatorname{cov}[I(x, y)]=\left[

    μ20′μ11′μ11′μ02′" role="presentation" style="position: relative;">μ20′μ11′μ11′μ02′
    \right] cov[I(x,y)]=[μ20′​μ11′​​μ11′​μ02′​​]
    协方差矩阵的特征值为：
    λ i = μ 20 ′ + μ 02 ′ 2 ± 4 μ 11 ′ 2 + ( μ 20 ′ − μ 02 ′ ) 2 2 \lambda_{i}=\frac{\mu_{20}^{\prime}+\mu_{02}^{\prime}}{2} \pm \frac{\sqrt{4 \mu_{11}^{\prime 2}+\left(\mu_{20}^{\prime}-\mu_{02}^{\prime}\right)^{2}}}{2} λi​=2μ20′​+μ02′​​±24μ11′2​+(μ20′​−μ02′​)2 ​​
    其和特征向量轴长度的平方成正比。因此，特征值大小的相对差异可以指示图像的偏心率，或它被拉长的程度。偏心率计算公式为：
    1 − λ 2 λ 1 \sqrt{1-\frac{\lambda_{2}}{\lambda_{1}}} 1−λ1​λ2​​ ​

    • centroid_local
    @property
def centroid_local(self):
    M = self.moments
    M0 = M[(0,) * self._ndim]

    def _get_element(axis):
        return (0,) * axis + (1,) + (0,) * (self._ndim - 1 - axis)

    return np.asarray(
        tuple(M[_get_element(axis)] / M0 for axis in range(self._ndim)))

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10

    扩展 （矩）

    参考: https://en.wikipedia.org/wiki/Image_moment

    原始矩
    图像处理与分析方法介绍
    ppt 0星 超过10%的资源 24.33MB
    下载

    对于二维连续函数f(x,y)，(p + q)阶矩(有时称为原始矩)定义为：
    M p q = ∫ − ∞ ∞ ∫ − ∞ ∞ x p y q f ( x , y ) d x d y M_{p q}=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} x^{p} y^{q} f(x, y) d x d y Mpq​=∫−∞∞​∫−∞∞​xpyqf(x,y)dxdy
    对于p,q = 0,1,2，…将其适应于像素强度为I(x,y)的标量(灰度)图像，原始图像矩Mij计算为：
    M i j = ∑ x ∑ y x i y j I ( x , y ) M_{i j}=\sum_{x} \sum_{y} x^{i} y^{j} I(x, y) Mij​=x∑​y∑​xiyjI(x,y)
    在某些情况下，这可以通过将图像考虑为概率密度函数来计算，即，通过除以上面的
    ∑ x ∑ y I ( x , y ) \sum_{x} \sum_{y} I(x, y) x∑​y∑​I(x,y)
    一个唯一性定理(Hu[1962])指出，如果f(x,y)是分段连续的，且仅在xy平面的有限部分上具有非零值，则所有阶矩都存在，且矩序列(Mpq)是由f(x,y)唯一确定的。相反，(Mpq)唯一确定f(x,y)。在实际应用中，图像是用几个低阶矩的函数来概括的。

    通过原始矩可计算的简单图像属性有：

    1. 面积 (用于二值图像)或灰度总和(用于灰度图像)： M 00 M_{00} M00​

    2. 中心： { x ˉ , y ˉ } = { M 10 M 00 , M 01 M 00 } \{\bar{x}, \bar{y}\}=\left\{\frac{M_{10}}{M_{00}}, \frac{M_{01}}{M_{00}}\right\} {xˉ,yˉ​}={M00​M10​​,M00​M01​​}

    中心矩

    中心矩定义为：
    μ p q = ∫ − ∞ ∞ ∫ − ∞ ∞ ( x − x ˉ ) p ( y − y ˉ ) q f ( x , y ) d x d y \mu_{p q}=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty}(x-\bar{x})^{p}(y-\bar{y})^{q} f(x, y) d x d y μpq​=∫−∞∞​∫−∞∞​(x−xˉ)p(y−yˉ​)qf(x,y)dxdy
    其中， x ˉ = M 10 M 00 \bar{x}=\frac{M_{10}}{M_{00}} xˉ=M00​M10​​， y ˉ = M 01 M 00 \bar{y}=\frac{M_{01}}{M_{00}} yˉ​=M00​M01​​，是质心的分量

    如果 f ( x , y ) f(x, y) f(x,y)是数字图像，上面的公式变为：
    μ p q = ∑ x ∑ y ( x − x ˉ ) p ( y − y ˉ ) q f ( x , y ) \mu_{p q}=\sum_{x} \sum_{y}(x-\bar{x})^{p}(y-\bar{y})^{q} f(x, y) μpq​=x∑​y∑​(x−xˉ)p(y−yˉ​)qf(x,y)
    三阶中心矩为：
    μ 00 = M 00 , μ 01 = 0 , μ 10 = 0 , μ 11 = M 11 − x ˉ M 01 = M 11 − y ˉ M 10 , μ 20 = M 20 − x ˉ M 10 , μ 02 = M 02 − y ˉ M 01 , μ 21 = M 21 − 2 x ˉ M 11 − y ˉ M 20 + 2 x ˉ 2 M 01 , μ 12 = M 12 − 2 y ˉ M 11 − x ˉ M 02 + 2 y ˉ 2 M 10 , μ 30 = M 30 − 3 x ˉ M 20 + 2 x ˉ 2 M 10 , μ 03 = M 03 − 3 y ˉ M 02 + 2 y ˉ 2 M 01 .

    μ00=M00,μ01=0,μ10=0,μ11=M11−x¯M01=M11−y¯M10,μ20=M20−x¯M10,μ02=M02−y¯M01,μ21=M21−2x¯M11−y¯M20+2x¯2M01,μ12=M12−2y¯M11−x¯M02+2y¯2M10,μ30=M30−3x¯M20+2x¯2M10,μ03=M03−3y¯M02+2y¯2M01." role="presentation" style="position: relative;">μ00=M00,μ01=0,μ10=0,μ11=M11−x¯M01=M11−y¯M10,μ20=M20−x¯M10,μ02=M02−y¯M01,μ21=M21−2x¯M11−y¯M20+2x¯2M01,μ12=M12−2y¯M11−x¯M02+2y¯2M10,μ30=M30−3x¯M20+2x¯2M10,μ03=M03−3y¯M02+2y¯2M01.
    ​μ00​=M00​,μ01​=0,μ10​=0,μ11​=M11​−xˉM01​=M11​−yˉ​M10​,μ20​=M20​−xˉM10​,μ02​=M02​−yˉ​M01​,μ21​=M21​−2xˉM11​−yˉ​M20​+2xˉ2M01​,μ12​=M12​−2yˉ​M11​−xˉM02​+2yˉ​2M10​,μ30​=M30​−3xˉM20​+2xˉ2M10​,μ03​=M03​−3yˉ​M02​+2yˉ​2M01​.​
    可以写成：
    μ p q = ∑ m p ∑ n q ( p m ) ( q n ) ( − x ˉ ) ( p − m ) ( − y ˉ ) ( q − n ) M m n \mu_{p q}=\sum_{m}^{p} \sum_{n}^{q}\left(
    pm" role="presentation" style="position: relative;">pm
    \right)\left(
    qn" role="presentation" style="position: relative;">qn
    \right)(-\bar{x})^{(p-m)}(-\bar{y})^{(q-n)} M_{m n}     μpq​=m∑p​n∑q​(pm​)(qn​)(−xˉ)(p−m)(−yˉ​)(q−n)Mmn​
    中心矩是平移不变的。

    利用二阶中心矩构造协方差矩阵，可以得到图像的方向信息。
    μ 20 ′ = μ 20 / μ 00 = M 20 / M 00 − x ˉ 2 μ 02 ′ = μ 02 / μ 00 = M 02 / M 00 − y ˉ 2 μ 11 ′ = μ 11 / μ 00 = M 11 / M 00 − x ˉ y ˉ

    μ20′=μ20/μ00=M20/M00−x¯2μ02′=μ02/μ00=M02/M00−y¯2μ11′=μ11/μ00=M11/M00−x¯y¯" role="presentation" style="position: relative;">μ20′=μ20/μ00=M20/M00−x¯2μ02′=μ02/μ00=M02/M00−y¯2μ11′=μ11/μ00=M11/M00−x¯y¯
    ​μ20′​=μ20​/μ00​=M20​/M00​−xˉ2μ02′​=μ02​/μ00​=M02​/M00​−yˉ​2μ11′​=μ11​/μ00​=M11​/M00​−xˉyˉ​​

    图像 I ( x , y ) I(x, y) I(x,y)的协方差矩阵变为：
    cov ⁡ [ I ( x , y ) ] = [ μ 20 ′ μ 11 ′ μ 11 ′ μ 02 ′ ] \operatorname{cov}[I(x, y)]=\left[

    μ20′μ11′μ11′μ02′" role="presentation" style="position: relative;">μ20′μ11′μ11′μ02′
    \right] cov[I(x,y)]=[μ20′​μ11′​​μ11′​μ02′​​]
    协方差矩阵的特征值为：
    λ i = μ 20 ′ + μ 02 ′ 2 ± 4 μ 11 ′ 2 + ( μ 20 ′ − μ 02 ′ ) 2 2 \lambda_{i}=\frac{\mu_{20}^{\prime}+\mu_{02}^{\prime}}{2} \pm \frac{\sqrt{4 \mu_{11}^{\prime 2}+\left(\mu_{20}^{\prime}-\mu_{02}^{\prime}\right)^{2}}}{2} λi​=2μ20′​+μ02′​​±24μ11′2​+(μ20′​−μ02′​)2 ​​
    其和特征向量轴长度的平方成正比。因此，特征值大小的相对差异可以指示图像的偏心率，或它被拉长的程度。偏心率计算公式为：
    1 − λ 2 λ 1 \sqrt{1-\frac{\lambda_{2}}{\lambda_{1}}} 1−λ1​λ2​​ ​

    不变矩

    矩因其在图像分析中的应用而闻名，因为它们可以用于获得关于特定变换类的不变量。

    在这种情况下，术语 不变矩经常出现。然而，当不变量由矩构成时，唯一的本身是不变量的是中心矩。

    注意，下面详述的不变量仅在连续域中是完全不变的。在离散域中，缩放和旋转都没有很好的定义:以这种方式变换的离散图像通常是一种近似，并且变换是不可逆的。因此，当描述离散图像中的形状时，这些不变量只是近似不变的。

    平移不变

    所有阶的中心矩 μ i j μ_{ij} μij​，通过构造，对于平移是不变的。

    缩放不变

    平移和尺度上的不变量 η i j η_{ij} ηij​可以由中心矩构成，通过除以适当比例的零点中心矩:
    η i j = μ i j μ 00 ( 1 + i + j 2 ) \eta_{i j}=\frac{\mu_{i j}}{\mu_{00}^{\left(1+\frac{i+j}{2}\right)}} ηij​=μ00(1+2i+j​)​μij​​
    其中， i + j ≥ 2 i+j \geq 2 i+j≥2。注意，只有在使用中心矩的情况下试平移不变的。

    旋转不变

    如Hu[1][2]的工作所示，可以构造关于平移、缩放和旋转的不变量:
    I 1 = η 20 + η 02 I 2 = ( η 20 − η 02 ) 2 + 4 η 11 2 I 3 = ( η 30 − 3 η 12 ) 2 + ( 3 η 21 − η 03 ) 2 I 4 = ( η 30 + η 12 ) 2 + ( η 21 + η 03 ) 2 I 5 = ( η 30 − 3 η 12 ) ( η 30 + η 12 ) [ ( η 30 + η 12 ) 2 − 3 ( η 21 + η 03 ) 2 ] + ( 3 η 21 − η 03 ) ( η 21 + η 03 ) [ 3 ( η 30 + η 12 ) 2 − ( η 21 + η 03 ) 2 ] I 6 = ( η 20 − η 02 ) [ ( η 30 + η 12 ) 2 − ( η 21 + η 03 ) 2 ] + 4 η 11 ( η 30 + η 12 ) ( η 21 + η 03 ) I 7 = ( 3 η 21 − η 03 ) ( η 30 + η 12 ) [ ( η 30 + η 12 ) 2 − 3 ( η 21 + η 03 ) 2 ] − ( η 30 − 3 η 12 ) ( η 21 + η 03 ) [ 3 ( η 30 + η 12 ) 2 − ( η 21 + η 03 ) 2 ]

    I1=η20+η02I2=(η20−η02)2+4η112I3=(η30−3η12)2+(3η21−η03)2I4=(η30+η12)2+(η21+η03)2I5=(η30−3η12)(η30+η12)[(η30+η12)2−3(η21+η03)2]+(3η21−η03)(η21+η03)[3(η30+η12)2−(η21+η03)2]I6=(η20−η02)[(η30+η12)2−(η21+η03)2]+4η11(η30+η12)(η21+η03)I7=(3η21−η03)(η30+η12)[(η30+η12)2−3(η21+η03)2]−(η30−3η12)(η21+η03)[3(η30+η12)2−(η21+η03)2]" role="presentation" style="position: relative;">I1=η20+η02I2=(η20−η02)2+4η112I3=(η30−3η12)2+(3η21−η03)2I4=(η30+η12)2+(η21+η03)2I5=(η30−3η12)(η30+η12)[(η30+η12)2−3(η21+η03)2]+(3η21−η03)(η21+η03)[3(η30+η12)2−(η21+η03)2]I6=(η20−η02)[(η30+η12)2−(η21+η03)2]+4η11(η30+η12)(η21+η03)I7=(3η21−η03)(η30+η12)[(η30+η12)2−3(η21+η03)2]−(η30−3η12)(η21+η03)[3(η30+η12)2−(η21+η03)2]
    ​I1​=η20​+η02​I2​=(η20​−η02​)2+4η112​I3​=(η30​−3η12​)2+(3η21​−η03​)2I4​=(η30​+η12​)2+(η21​+η03​)2I5​=(η30​−3η12​)(η30​+η12​)[(η30​+η12​)2−3(η21​+η03​)2]+(3η21​−η03​)(η21​+η03​)[3(η30​+η12​)2−(η21​+η03​)2]I6​=(η20​−η02​)[(η30​+η12​)2−(η21​+η03​)2]+4η11​(η30​+η12​)(η21​+η03​)I7​=(3η21​−η03​)(η30​+η12​)[(η30​+η12​)2−3(η21​+η03​)2]−(η30​−3η12​)(η21​+η03​)[3(η30​+η12​)2−(η21​+η03​)2]​
    这些就是众所周知的Hu矩不变量。

    第一个 I 1 I_1 I1​，类似于围绕图像质心的转动惯量，其中像素的强度类似于物理密度。前六个 I 1 , . . . I 6 I_1,...I_6 I1​,...I6​反射对称，比如，如果将图像更改为镜像，则这些参数不变。最后一个 I 7 I_7 I7​，是反射反对称（反射情况下符号改变），这使它能够区分其他相同图像的镜像。

    J. Flusser.[3]提出了获取旋转矩不变量完备独立集的一般理论。他证明了传统的胡矩不变量集既不是独立的，也不是完全的。I3不是很有用，因为它依赖于其他(如何?)在原Hu集合中缺失一个三阶独立不变量:
    I 8 = η 11 [ ( η 30 + η 12 ) 2 − ( η 03 + η 21 ) 2 ] − ( η 20 − η 02 ) ( η 30 + η 12 ) ( η 03 + η 21 ) I_{8}=\eta_{11}\left[\left(\eta_{30}+\eta_{12}\right)^{2}-\left(\eta_{03}+\eta_{21}\right)^{2}\right]-\left(\eta_{20}-\eta_{02}\right)\left(\eta_{30}+\eta_{12}\right)\left(\eta_{03}+\eta_{21}\right) I8​=η11​[(η30​+η12​)2−(η03​+η21​)2]−(η20​−η02​)(η30​+η12​)(η03​+η21​)
    和 I 7 I_7 I7​一样， I 8 I_8 I8​也是反射反对称的。

    后来，J. Flusser和T. Suk[4]将n个旋转对称形状的理论专门化。

    应用

    Zhang等人应用Hu矩不变量来解决病理脑检测(PBD)问题Doerr和Florence利用与二阶中心矩相关的目标方向信息，从微x射线层析成像图像数据中有效提取平移和旋转不变的目标截面[6]。

    equivalent_diameter_area

    参考：W. Pabst, E. Gregorová. Characterization of particles and particle
    systems, pp. 27-28. ICT Prague, 2007.
    https://old.vscht.cz/sil/keramika/Characterization_of_particles/CPPS%20_English%20version_.pdf

  • 相关阅读:
    德克萨斯大学奥斯汀分校自然语言处理硕士课程汉化版(第九周) - 可解释性
    flex布局原理，常见的父项属性,flex布局子项常见属性，利用flex布局和float制作的有轮播图的携程旅行首页
    创意电子学-小知识：如何使用面包板
    webservice接口测试
    北京化工大学数据结构2022/11/17作业 题解
    JAVA生成图片缩略图、JAVA截取图片局部内容
    DolphinScheduler 集群部署
    2020年数维杯数学建模A题舆情监测情感倾向分析建模求解全过程文档及程序
    【Compose】Desktop Application 初尝试
    测试一下 Baichuan2-7B-Chat 的性能
  • 原文地址：https://blog.csdn.net/weixin_44456692/article/details/126134136