图像处理之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
       连通域属性列表，可用的属性在下一节
    """

属性及数据类型

{#********几何特征********#
    '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零阶空间矩(强度权重面积) 
}

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

其实例化了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)

area

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

使用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)])

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

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

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

coords、centroid

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

质心即连通域坐标取平均值即可，其中获取连通域坐标的函数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

这里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)

凸包图像计算，使用了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)

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)

偏心率 和连通域有相同二阶中心矩的椭圆的离心率。公式：焦距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)

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

@property
@_cached
def inertia_tensor(self):
    mu = self.moments_central
    return _moments.inertia_tensor(self.image, mu)

其调用了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	# 惯性张量特征值
# 然后根据惯性张量特征值计算离心率。

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

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)

其调用了moments_central

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

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

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

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

中心矩定义为：
$${\mu }_{pq}={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }(x-\bar{x}{)}^p(y-\bar{y}{)}^{q}f(x,y)dxdy$$
其中，$\bar{x}=\frac{M_{10}}{M_{00}}$，$\bar{y}=\frac{M_{01}}{M_{00}}$，是质心的分量

如果$f(x,y)$是数字图像，上面的公式变为：
$${\mu }_{pq}=\sum _x\sum _y(x-\bar{x}{)}^p(y-\bar{y}{)}^{q}f(x,y)$$
三阶中心矩为：
$$\begin{array}{rl}& {\mu }_{00}=M_{00},\\ & {\mu }_{01}=0,\\ & {\mu }_{10}=0,\\ & {\mu }_{11}=M_{11}-\bar{x}{M}_{01}=M_{11}-\bar{y}{M}_{10},\\ & {\mu }_{20}=M_{20}-\bar{x}{M}_{10},\\ & {\mu }_{02}=M_{02}-\bar{y}{M}_{01},\\ & {\mu }_{21}=M_{21}-2\bar{x}{M}_{11}-\bar{y}{M}_{20}+2{\bar{x}}^2{M}_{01},\\ & {\mu }_{12}=M_{12}-2\bar{y}{M}_{11}-\bar{x}{M}_{02}+2{\bar{y}}^2{M}_{10},\\ & {\mu }_{30}=M_{30}-3\bar{x}{M}_{20}+2{\bar{x}}^2{M}_{10},\\ & {\mu }_{03}=M_{03}-3\bar{y}{M}_{02}+2{\bar{y}}^2{M}_{01}.\end{array}$$
可以写成：
$${\mu }_{pq}=\sum _m^p\sum _n^q\left(\begin{array}{c}p\\ m\end{array}\right)\left(\begin{array}{c}q\\ n\end{array}\right)(-\bar{x}{)}^{(p-m)}(-\bar{y}{)}^{(q-n)}{M}_{mn}$$
中心矩是平移不变的。

利用二阶中心矩构造协方差矩阵，可以得到图像的方向信息。
$$\begin{array}{rl}& {\mu }_{20}^{\prime }={\mu }_{20}/{\mu }_{00}=M_{20}/M_{00}-{\bar{x}}^2\\ & {\mu }_{02}^{\prime }={\mu }_{02}/{\mu }_{00}=M_{02}/M_{00}-{\bar{y}}^2\\ & {\mu }_{11}^{\prime }={\mu }_{11}/{\mu }_{00}=M_{11}/M_{00}-\bar{x}\bar{y}\end{array}$$

图像$I(x,y)$的协方差矩阵变为：
$$\mathrm{cov}[I(x,y)]=\left[\begin{array}{ll}{\mu }_{20}^{\prime }& {\mu }_{11}^{\prime }\\ {\mu }_{11}^{\prime }& {\mu }_{02}^{\prime }\end{array}\right]$$
协方差矩阵的特征值为：
$${\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}$$
其和特征向量轴长度的平方成正比。因此，特征值大小的相对差异可以指示图像的偏心率，或它被拉长的程度。偏心率计算公式为：
$$\sqrt{1-\frac{{\lambda }_2}{{\lambda }_1}}$$

• 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)))

扩展 （矩）

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

原始矩

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

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

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

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​​}

中心矩

中心矩定义为：
$${\mu }_{pq}={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }(x-\bar{x}{)}^p(y-\bar{y}{)}^{q}f(x,y)dxdy$$
其中，$\bar{x}=\frac{M_{10}}{M_{00}}$，$\bar{y}=\frac{M_{01}}{M_{00}}$，是质心的分量

如果$f(x,y)$是数字图像，上面的公式变为：
$${\mu }_{pq}=\sum _x\sum _y(x-\bar{x}{)}^p(y-\bar{y}{)}^{q}f(x,y)$$
三阶中心矩为：
$$\begin{array}{rl}& {\mu }_{00}=M_{00},\\ & {\mu }_{01}=0,\\ & {\mu }_{10}=0,\\ & {\mu }_{11}=M_{11}-\bar{x}{M}_{01}=M_{11}-\bar{y}{M}_{10},\\ & {\mu }_{20}=M_{20}-\bar{x}{M}_{10},\\ & {\mu }_{02}=M_{02}-\bar{y}{M}_{01},\\ & {\mu }_{21}=M_{21}-2\bar{x}{M}_{11}-\bar{y}{M}_{20}+2{\bar{x}}^2{M}_{01},\\ & {\mu }_{12}=M_{12}-2\bar{y}{M}_{11}-\bar{x}{M}_{02}+2{\bar{y}}^2{M}_{10},\\ & {\mu }_{30}=M_{30}-3\bar{x}{M}_{20}+2{\bar{x}}^2{M}_{10},\\ & {\mu }_{03}=M_{03}-3\bar{y}{M}_{02}+2{\bar{y}}^2{M}_{01}.\end{array}$$
可以写成：
$${\mu }_{pq}=\sum _m^p\sum _n^q\left(\begin{array}{c}p\\ m\end{array}\right)\left(\begin{array}{c}q\\ n\end{array}\right)(-\bar{x}{)}^{(p-m)}(-\bar{y}{)}^{(q-n)}{M}_{mn}$$
中心矩是平移不变的。

利用二阶中心矩构造协方差矩阵，可以得到图像的方向信息。
$$\begin{array}{rl}& {\mu }_{20}^{\prime }={\mu }_{20}/{\mu }_{00}=M_{20}/M_{00}-{\bar{x}}^2\\ & {\mu }_{02}^{\prime }={\mu }_{02}/{\mu }_{00}=M_{02}/M_{00}-{\bar{y}}^2\\ & {\mu }_{11}^{\prime }={\mu }_{11}/{\mu }_{00}=M_{11}/M_{00}-\bar{x}\bar{y}\end{array}$$

图像$I(x,y)$的协方差矩阵变为：
$$\mathrm{cov}[I(x,y)]=\left[\begin{array}{ll}{\mu }_{20}^{\prime }& {\mu }_{11}^{\prime }\\ {\mu }_{11}^{\prime }& {\mu }_{02}^{\prime }\end{array}\right]$$
协方差矩阵的特征值为：
$${\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}$$
其和特征向量轴长度的平方成正比。因此，特征值大小的相对差异可以指示图像的偏心率，或它被拉长的程度。偏心率计算公式为：
$$\sqrt{1-\frac{{\lambda }_2}{{\lambda }_1}}$$

不变矩

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

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

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

平移不变

所有阶的中心矩${\mu }_{ij}$，通过构造，对于平移是不变的。

缩放不变

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

旋转不变

如Hu[1][2]的工作所示，可以构造关于平移、缩放和旋转的不变量:
$$\begin{array}{rl}& I_1={\eta }_{20}+{\eta }_{02}\\ & I_2={\left({\eta }_{20}-{\eta }_{02}\right)}^2+4{\eta }_{11}^2\\ & I_3={\left({\eta }_{30}-3{\eta }_{12}\right)}^2+{\left(3{\eta }_{21}-{\eta }_{03}\right)}^2\\ & I_4={\left({\eta }_{30}+{\eta }_{12}\right)}^2+{\left({\eta }_{21}+{\eta }_{03}\right)}^2\\ & I_5=\left({\eta }_{30}-3{\eta }_{12}\right)\left({\eta }_{30}+{\eta }_{12}\right)\left[{\left({\eta }_{30}+{\eta }_{12}\right)}^2-3{\left({\eta }_{21}+{\eta }_{03}\right)}^2\right]+\left(3{\eta }_{21}-{\eta }_{03}\right)\left({\eta }_{21}+{\eta }_{03}\right)\left[3{\left({\eta }_{30}+{\eta }_{12}\right)}^2-{\left({\eta }_{21}+{\eta }_{03}\right)}^2\right]\\ & I_6=\left({\eta }_{20}-{\eta }_{02}\right)\left[{\left({\eta }_{30}+{\eta }_{12}\right)}^2-{\left({\eta }_{21}+{\eta }_{03}\right)}^2\right]+4{\eta }_{11}\left({\eta }_{30}+{\eta }_{12}\right)\left({\eta }_{21}+{\eta }_{03}\right)\\ & I_7=\left(3{\eta }_{21}-{\eta }_{03}\right)\left({\eta }_{30}+{\eta }_{12}\right)\left[{\left({\eta }_{30}+{\eta }_{12}\right)}^2-3{\left({\eta }_{21}+{\eta }_{03}\right)}^2\right]-\left({\eta }_{30}-3{\eta }_{12}\right)\left({\eta }_{21}+{\eta }_{03}\right)\left[3{\left({\eta }_{30}+{\eta }_{12}\right)}^2-{\left({\eta }_{21}+{\eta }_{03}\right)}^2\right]\end{array}$$
这些就是众所周知的Hu矩不变量。

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

J. Flusser.[3]提出了获取旋转矩不变量完备独立集的一般理论。他证明了传统的胡矩不变量集既不是独立的，也不是完全的。I3不是很有用，因为它依赖于其他(如何?)在原Hu集合中缺失一个三阶独立不变量:
$$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)$$
和$I_7$一样，$I_8$也是反射反对称的。

后来，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