图像特征Vol.1：计算机视觉特征度量|第二弹：【统计区域度量】

一、前言

在前篇图像特征Vol.1：计算机视觉特征度量【纹理区域特征】中，我们说到计算机视觉度量的三类方法，在那篇博客中，我们介绍了纹理区域度量的各个方法。在本篇博客中，我们将继续介绍纹理区域特征的第二类方法：统计区域度量。Let’s Go！

二、统计区域度量

统计区域度量，是利用统计方法，描述图像区域中像素值分布和变化情况的一种方法。 它主要应用在图像处理和计算机视觉领域。对于图像处理方面，在前面的博客【cv】图像预处理技术——从特征检测讲述图像预处理理论、实践、应用|01中，利用全局直方图和局部直方图，对图像进行光照校正、增强对比度，就是统计区域度量在图像处理方面应用的一个体现。这篇博客中，我们将介绍统计区域度量在计算机视觉特征度量中的应用体现。

常见的图像像素统计度量包括以下八种

1. 图像矩特征
2. 点度量特征
3. 全局直方图
4. 局部直方图
5. 散点图和3D直方图
6. 多尺度直方图
7. 径向直方图
8. 轮廓或边缘直方图

这些统计度量可以通过对图像像素值进行数学运算和分析来得到。它们可以为图像处理和计算机视觉算法提供重要的特征和信息，从而提高算法的性能和鲁棒性。

接下来，我会基于每种方法给出基本原理，和步骤，以及一些简单的应用举例~Go

2.1：图像矩特征

👧🏻 什么是矩和图像矩？
 矩是概率与统计中的一个概念，是随机变量的一种数字特征。在数字图像处理、计算机视觉与相关领域中，图像矩(Image moments)是指图像的某些特定像素灰度的加权平均值（矩），或者是图像具有类似功能或意义的属性。

 图像矩通常用来描述 分割后的图像对象。可以通过图像的矩来获得图像的部分性质，包括面积(或总体亮度)，以及有关几何中心和方向的信息 。

💁🏻‍♀️基本思想和原理？
 实际上，对于矩，我们实际上可以看作是二维图像在基空间上的投影。 为什么这么说呢？接着看完它的基本思想和原理你应该就懂了。

根据加权求和的方法不同，有不同种类的矩，且每一种矩有各自不同的表示图像特征的特征值。

2.1.1：原始矩/几何矩

对于灰度图像f,其(p+q)阶矩(有时称为"原始矩"),定义如下:

原始矩阵中包含以下可以描述图像特征的信息：

1. 二值图像或者或度图像的像素总和,可以表示为:
   
2. 图像的几何中心也叫质心 (Cx,Cy) ,计算公式如下
   

2.1.2：中心距

当图像平移时，原始矩会发生变化（质心肯定会因为平移而改变嘛）。但是中心矩具有平移不变性（但在图像旋转时会发生变化；归一化中心距不仅具有平移不变性，而且具有比例不变性（尺度不变性））。相比于原始矩，它添加了一个基于质心的平移，对于灰度图像f,其(p+q)阶中心矩的定义如下：

ps.述公式中（cx, cy)是用原始矩求出的图像的几何质心

中心矩可以用来计算图像的方向和形状等特征，不管是一维的序列分布还是二维的图像，不同的阶次均具备下述视觉上的规律性性质：

1. 零阶矩：表示一维均值或二维质心；
2. 一阶中心矩：包含二维面积、质心和物体/目标大小等相关信息；
3. 二阶中心矩：与方差和2D椭圆度量相关；
4. 三阶中心矩：提供了二维形状（或偏度）的对称信息；
5. 四阶中心矩：用来度量二维分布，如高、矮、细、短、胖等形态；
6. 更高阶的矩：可由多个矩的比值组成，比如协方差。

2.1.3：归一化的中心矩

当图像发生平移时，几何距 mpq ​也会发生变化；中心矩 mupq具有平移不变性，但在图像旋转时会发生变化；归一化中心距不仅具有平移不变性，而且具有比例不变性（尺度不变性）。归一化的中心矩定义如下：

2.1.4：不变矩——Hu矩

Hu 利用二阶和三阶归一化中心距构造了 7 个不变矩 M1~M7， 在连续图像下具有平移、灰度、尺度、旋转不变性， 是高度浓缩的图像特征。不变矩能够描述图像的整体性质，从而在边缘提取、图像匹配及目标识别中得到了广泛的应用。

在实际的图像识别中，只有 M1 和 M2 的不变性比较好，而 M3~M7 的误差较大，识别率比较低。Hu不变矩对于物体的形状描述比较稳定，比较适合识别较大尺寸的物体，如水果的形状、车牌的字符等。

2.1.5：OpenCv实现矩特征及其应用

OpenCV 提供了函数 cv2.moments() 计算图像矩 Mu，函数 cv2.HuMoments() 计算目标轮廓的 Hu 不变矩。

函数说明：

cv.moments(array[, binaryImage]) → Mu
cv.HuMoments(Mu[, hu]) → Hu

参数说明：

1. array：是一幅单通道 8-bits 图像，或一个二维浮点数组（如轮廓列表 contours 中的一个轮廓）
2. binaryImage：指示输入图像是否二值图像

• ⭐函数 cv2.moments() 以字典（Dict）形式返回图像的矩。
  返回值 Mu 是字典格式，包括 24个键值对。[‘m00’,‘m10’,‘m01’,‘m20’,‘m11’,‘m02’,‘m30’,‘m21’,‘m12’,‘m03’,‘mu20’,‘mu11’,‘mu02’,‘mu30’,‘mu21’,‘mu12’,‘mu03’,‘nu20’,‘nu11’,‘nu02’,‘nu30’,‘nu21’,nu12’,‘nu03’]
  1. 前 10个键值 [‘m00’,‘m10’,‘m01’,‘m20’,‘m11’,‘m02’,‘m30’,‘m21’,‘m12’,‘m03’] 代表几何矩（p+q阶矩），也称原点矩
  2. 中间 7 个键值 [‘mu20’,‘mu11’,‘mu02’,‘mu30’,‘mu21’,‘mu12’,‘mu03’] 代表中心矩
  3. 最后 7 个键值 [‘nu20’,‘nu11’,‘nu02’,‘nu30’,‘nu21’,nu12’,‘nu03’] 代表归一化的中心矩
• ⭐ 函数 cv2.HuMoments() 以列表（List）形式返回目标轮廓的 Hu 不变矩。
  返回值 Hu 是一个列表（List），包括 7 个不变矩 M1~M7，浮点数。

🧑🏻‍🚀实例1、基于轮廓计算矩及其面积：

import cv2
import numpy as np

o = cv2.imread('2.jpg')
cv2.imshow("original",o)
gray = cv2.cvtColor(o,cv2.COLOR_BGR2GRAY)
ret, binary = cv2.threshold(gray,127,255,cv2.THRESH_BINARY) # 转成二值图
contours, hierarchy = cv2.findContours(binary,
 cv2.RETR_LIST,
 cv2.CHAIN_APPROX_SIMPLE)
n=len(contours)
contoursImg=[]
for i in range(n):         # 分别展示轮廓
     temp=np.zeros(o.shape,np.uint8)
     contoursImg.append(temp)
     contoursImg[i]=cv2.drawContours(contoursImg[i],contours,i,255,3)
     cv2.imshow("contours[" + str(i)+"]",contoursImg[i])
print("观察各个轮廓的矩（moments）:")
for i in range(n):
     print("轮廓"+str(i)+"的矩:\n",cv2.moments(contours[i]))
print("观察各个轮廓的面积:")
for i in range(n):
     print("轮廓"+str(i)+"的面积:%d" %cv2.moments(contours[i])['m00'])

cv2.waitKey()
cv2.destroyAllWindows()

输出结果：

观察各个轮廓的矩（moments）:
轮廓0的矩:
 {'m00': 1456.0, 'm10': 58968.0, 'm01': 242400.8333333333, 'm20': 2504655.333333333, 'm11': 9817233.75, 'm02': 40623554.166666664, 'm30': 110871099.0, 'm21': 416981784.45, 'm12': 1645253943.75, 'm03': 6852302097.35, 'mu20': 116451.33333333302, 'mu11': 0.0, 'mu02': 267672.2980578467, 'mu30': 4.470348358154297e-08, 'mu21': -3476.5477868951857, 'mu12': -8.568167686462402e-08, 'mu03': 339.670916557312, 'nu20': 0.05493153252425217, 'nu11': 0.0, 'nu02': 0.1262643297051641, 'nu30': 5.526344190010775e-16, 'nu21': -4.297785792992532e-05, 'nu12': -1.0592159697520654e-15, 'nu03': 4.199087511397351e-06}
轮廓1的矩:
 {'m00': 3443.5, 'm10': 420127.3333333333, 'm01': 306481.8333333333, 'm20': 53186976.416666664, 'm11': 37392173.125, 'm02': 27783857.25, 'm30': 6959813898.700001, 'm21': 4733709260.066667, 'm12': 3389725532.9, 'm03': 2562927246.9500003, 'mu20': 1928960.9632681608, 'mu11': -420.26935005933046, 'mu02': 506054.385658171, 'mu30': -540.4464931488037, 'mu21': 11304.976172447205, 'mu12': 5699.499477289617, 'mu03': -140.53871870040894, 'nu20': 0.16267591725263073, 'nu11': -3.5442760800216457e-05, 'nu02': 0.042677308112644856, 'nu30': -7.766985329226633e-07, 'nu21': 1.6246859807910376e-05, 'nu12': 8.190991964093441e-06, 'nu03': -2.0197414178314225e-07}
轮廓2的矩:
 {'m00': 368.0, 'm10': 80776.0, 'm01': 7728.0, 'm20': 17746554.666666664, 'm11': 1696296.0, 'm02': 170138.66666666666, 'm30': 3902490500.0, 'm21': 372677648.0, 'm12': 37345437.333333336, 'm03': 3902640.0, 'mu20': 16222.666666664183, 'mu11': 0.0, 'mu02': 7850.666666666657, 'mu30': 1.430511474609375e-06, 'mu21': 5.21540641784668e-08, 'mu12': 4.6566128730773926e-09, 'mu03': 4.656612873077393e-10, 'nu20': 0.11979166666664831, 'nu11': 0.0, 'nu02': 0.05797101449275355, 'nu30': 5.50645125550301e-13, 'nu21': 2.0075603535688056e-14, 'nu12': 1.792464601400719e-15, 'nu03': 1.7924646014007192e-16}
观察各个轮廓的面积:
轮廓0的面积:1456
轮廓1的面积:3443
轮廓2的面积:368

🧑🏻‍🚀实例2、Hu矩进行形状匹配

下面这个代码是将原图像进行缩放和旋转等操作，再分别计算各自的Hu矩，根据输出结果可以看到，hu矩是具有平移不变性、旋转不变性和尺度不变性的。

# # 14.14 特征描述之不变矩
import cv2
import numpy as np
import matplotlib.pyplot as plt

gray = cv2.imread("Hu1.png", flags=0)  # 灰度图像
height, width = gray.shape  # (568, 568)

# 图像的平移，缩放，旋转和镜像
grayList = []
grayList.append(gray)  # [0]，原始图像
mat = np.float32([[1, 0, 50], [0, 1, 50]])
grayList.append(cv2.warpAffine(gray, mat, (height, width)))  # [1]，图像平移
top, bottom, left, right = height // 4, height // 4, width // 4, width // 4
grayResize = cv2.resize(gray, (width // 2, height // 2))  # 图像缩放 (284, 284)
replicate = cv2.copyMakeBorder(grayResize, top, bottom, left, right, cv2.BORDER_CONSTANT, value=0)
grayList.append(replicate)  # [2]，图像缩放并填充至原来尺寸 (568, 568)
grayList.append(cv2.flip(gray, 1))  # [3]，图像镜像，水平翻转
mar = cv2.getRotationMatrix2D((width // 2, height // 2), angle=45, scale=1)  # 图像中心作为旋转中心
rotate = cv2.warpAffine(gray, mar, (height, width))  # 旋转变换，默认为黑色填充
grayList.append(rotate)  # [4]，图像旋转 45度
grayList.append(cv2.rotate(gray, cv2.ROTATE_90_COUNTERCLOCKWISE))  # [5]，图像逆时针旋转90度
print(gray.shape, grayResize.shape, replicate.shape, rotate.shape, len(grayList))

plt.figure(figsize=(9, 6))
for i in range(len(grayList)):
    moments = cv2.moments(grayList[i])  # 返回几何矩 mpq, 中心矩 mupq 和归一化矩 nupq
    huM = cv2.HuMoments(moments)  # 计算 Hu 不变矩
    plt.subplot(2, 3, i + 1), plt.axis('off'), plt.imshow(grayList[i], 'gray')
    # print("Moments of gray:\n", moments)
    print("HuMoments of gray:\n", np.log10(np.abs(huM.T)).round(4))

plt.show()

输出结果：

HuMoments of gray:
 [[ -2.7619  -6.5682 -10.2914  -9.8954 -20.4914 -13.7315 -20.0114]]
HuMoments of gray:
 [[ -2.7198  -7.1721  -9.9113  -9.5982 -19.8569 -13.185  -19.3754]]
HuMoments of gray:
 [[ -2.7612  -6.5649 -10.2883  -9.8929 -20.4913 -13.7241 -20.0055]]
HuMoments of gray:
 [[ -2.7619  -6.5682 -10.2914  -9.8954 -20.4914 -13.7315 -20.0114]]
HuMoments of gray:
 [[ -2.7092  -6.5701  -9.6241 -10.073  -20.0139 -13.7624 -20.1516]]
HuMoments of gray:
 [[ -2.7619  -6.5682 -10.2914  -9.8954 -20.4914 -13.7315 -20.0114]]

OpenCV提供了一个名为matchShapes的实用程序函数，它可以获取两幅图像(或轮廓)，并使用Hu矩查找它们之间的距离。所以，你不需要明确地计算Hu力矩。只需要简单地对图像进行二值化并使用matchShapes。

方法介绍：

retval = cv2.matchShapes( contour1, contour2, method, parameter )

• contour1 第一个轮廓或者灰度图像

• contour2 第二个灰度或轮廓图像 ，不管是什么，这个函数都会提前计算好对象的Hu矩。

• method 比较两个对象的Hu 矩的方法(距离函数），ℎ𝑖𝐴和ℎ𝑖𝐵分别是对象 A 和对象 B 的 Hu 矩。

  
  以下代码实现的功能是分别将Hu1和Hu1自身比、Hu1和Hu2比、Hu1和Hu3比
  

import cv2

o1 = cv2.imread('Hu4.png')
o2 = cv2.imread('Hu5.png')
o3 = cv2.imread('Hu6.png')
gray1 = cv2.cvtColor(o1,cv2.COLOR_BGR2GRAY)
gray2 = cv2.cvtColor(o2,cv2.COLOR_BGR2GRAY)
gray3 = cv2.cvtColor(o3,cv2.COLOR_BGR2GRAY)
ret, binary1 = cv2.threshold(gray1,127,255,cv2.THRESH_BINARY)
ret, binary2 = cv2.threshold(gray2,127,255,cv2.THRESH_BINARY)
ret, binary3 = cv2.threshold(gray3,127,255,cv2.THRESH_BINARY)
contours1, hierarchy = cv2.findContours(binary1,
 cv2.RETR_LIST,
 cv2.CHAIN_APPROX_SIMPLE)
contours2, hierarchy = cv2.findContours(binary2,
 cv2.RETR_LIST,
 cv2.CHAIN_APPROX_SIMPLE)
contours3, hierarchy = cv2.findContours(binary3,
 cv2.RETR_LIST,
 cv2.CHAIN_APPROX_SIMPLE)
cnt1 = contours1[0]
cnt2 = contours2[0]
cnt3 = contours3[0]
ret0 = cv2.matchShapes(cnt1,cnt1,1,0.0)
ret1 = cv2.matchShapes(cnt1,cnt2,1,0.0)
ret2 = cv2.matchShapes(cnt1,cnt3,1,0.0)
print("相同图像的 matchShape=",ret0)
print("相似图像的 matchShape=",ret1)
print("不相似图像的 matchShape=",ret2)

输出结果：

相同图像的 matchShape= 0.0
相似图像的 matchShape= 0.2135632086586597
不相似图像的 matchShape= 1.6939812402126888

🧑🏻‍🚀实例3、前景物体的中心点以及长轴的角度

# 利用矩求出前景物体的中心点以及长轴的角度,我们可以计算出大雁群的质心和长轴角度
import cv2
import numpy as np

# 读入图像
img_file = "1.jpg"
img = cv2.imread(img_file,cv2.IMREAD_COLOR)
# 将图像转换为灰度图像
gray_img = cv2.cvtColor(img,cv2.COLOR_RGB2GRAY)

# 二值化，分割前景物体
retval, bin_img = cv2.threshold(gray_img, 50, 1, cv2.THRESH_BINARY)
# show_img = bin_img * 255
show_img = (1 - bin_img) * 255

moments = cv2.moments(show_img)

cx = moments['m10']/moments['m00'] # m10/m00
cy = moments['m01']/moments['m00'] # m01/m00
print('Centriod: ({0:.2f}, {1:.2f})'.format(cx, cy))

mu00 = moments['m00']
mu11 =moments['mu11']
mu20 = moments['mu20']
mu02 = moments['mu02']
theta = 1/2*np.arctan2(2*mu11/mu00, (mu20 - mu02)/mu00)
angle_degrees = theta * 180 / np.pi
print('Angle: {0:.2f}'.format(angle_degrees))

# visual
rho = 800
dx_major = rho * np.cos(theta)
dy_major = rho * np.sin(theta)
dx_minor = 0.3 * rho * np.cos(theta - np.pi / 2)
dy_minor = 0.3 * rho * np.sin(theta - np.pi / 2)

# short
short_axis=[(int(cx-dx_minor),int(cy-dy_minor)),(int(cx),int(cy)),(int(cx+dx_minor),int(cy+dy_minor))]
for i in range(len(short_axis)-1):
    cv2.line(img,short_axis[i],short_axis[i+1],color=(255,0,0),thickness=2)
for pt in short_axis:
    cv2.circle(img,pt,radius=5,color=(255,0,0),thickness=3,lineType=8)
# long
long_axis = [(int(cx - dx_major), int(cy - dy_major)), (int(cx), int(cy)), (int(cx + dx_major), int(cy + dy_major))]
for i in range(len(long_axis) - 1):
    cv2.line(img, long_axis[i], long_axis[i + 1], color=(0, 0, 255), thickness=2)
for pt in long_axis:
    cv2.circle(img,pt,radius=5,color=(0,0,255),thickness=3,lineType=8)
# center
cv2.circle(img, (int(cx),int(cy)), radius=5, color=(0, 255, 0), thickness=3, lineType=8)
# show
cv2.imshow("img",img)
cv2.waitKey(0)

运行结果：

Centriod: (773.86, 497.30)
Angle: 36.76

2.2：点度量特征

点度量特征主要就是指描述图像中单个像素点的特征的方法。对于像素点，特征主要包括以下几种

• 灰度值
• 梯度
• 方向
• 角点
• 纹理特征

其中，前三个较为常见，前面也一直有用到。对于角点，属于兴趣点特征，涉及到特征提取方面的知识，这个在后续的会深入讲解。纹理特征呢，就不是单个像素点，而是多个像素点组合的特征度量，这个在前一篇纹理特征度量中也有很好的分析。

点特征度量的核心下游任务就是图像分割。

这里就简单对灰度、梯度、方向三个点特征进行阐述：

• 灰度
  $$I(x,y)$$

• 梯度

$$▽I(x,y)=(G_x,G_y)=\left(\begin{array}{c}\frac{\partial I}{\partial x}\\ \frac{\partial I}{\partial y}\end{array}\right)$$

• 方向
  $$\theta (x,y)=\arctan \frac{G_x}{G_y}$$

🧑🏻‍🚀实例1、基于HSV颜色阈值的图像分割

import cv2
import numpy as np
import matplotlib.pyplot as plt
# 1. 加载图像
img_file = "./lan.jpg"
image = cv2.imread(img_file,cv2.IMREAD_COLOR)

# 2. 颜色空间转换（从RGB到HSV）
hsv_image = cv2.cvtColor(image, cv2.COLOR_BGR2HSV)

# 3. 颜色阈值分割（这里选择蓝色的范围）
lower_blue = np.array([90, 50, 50])  # HSV中的下限
upper_blue = np.array([130, 255, 255])  # HSV中的上限

# 根据颜色阈值创建一个二进制掩码
mask = cv2.inRange(hsv_image, lower_blue, upper_blue)

# 4. 创建二值掩码
result = cv2.bitwise_and(image, image, mask=mask)

# 显示分割结果
plt.subplot(121)
plt.imshow(image)
plt.title('origin image')
plt.axis('off')
plt.subplot(122)
plt.imshow(result)
plt.title('Segmented Image')
plt.axis('off')

plt.show()

2.3：全局直方图和局部区域直方图

图像的全局直方图和局部区域直方图都是描述图像中像素值分布的方法，可以用来表示图像中每个灰度级别的像素数量。

全局直方图是一个一维向量，长度为图像的灰度级别数，向量元素值通常是该灰度值出现的像素个数或者像素频率。局部区域直方图就是加窗之后的窗内小图像的全局直方图而已。

下面是SAR图像的全局直方图示例，其中右边列为对应Image的全局直方图，直方图向量的分布变化可以作为诸如SAR图像变化检测的依据
a

2.4：散点图和3D直方图

散点图和3D直方图都是描述图像特征之间关系的方法。它们可以用来表示图像中不同特征之间的相关性和分布情况，从而提供图像的多维特征信息。

（1）散点图
散点图是二维坐标体系，其中每个点表示图像中两个特征之间的一个数据点。
$$S(i,j)=(f_i(x,y),f_j(x,y))|I(x,y)\in \Omega$$

其中S(i,j)表示特征fi和特征fj之间的散点图。比如我们要表征RGB三通道之间的相互关系，这里图像特征可以定义为不同通道下的像素值，I(x,y)为图像在(x,y)位置的像素点，Ω为图像的像素点集合。

(2) 3D直方图

图像3D直方图同样也是用来描述图像中不同特征之间的相关性和分布情况。

下面就以一个表示图像中颜色特征分布的示例来更好得了解一下

下面创建的3D直方图用于显示图像中红色通道和绿色通道的像素分布情况。这意味着x轴表示红色通道的像素值，y轴表示绿色通道的像素值，而z轴表示在这两个通道上像素值组合的频率。从图中可以了解哪些颜色组合在图像中较为常见，哪些较为罕见。这对于分析图像中的颜色特征和结构非常有帮助

# 这个示例演示了如何创建一个表示图像颜色分布的3D直方图
import cv2
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
# 读取图像
img = cv2.imread('./4.jpg')

# 将图像从BGR颜色空间转换为RGB
img_rgb = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)

# 获取图像的三个颜色通道
r, g, b = cv2.split(img_rgb)

# 创建3D直方图
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

hist, xedges, yedges = np.histogram2d(r.ravel(), g.ravel(), bins=(256, 256))

x, y = np.meshgrid(xedges[:-1], yedges[:-1])

# 绘制3D直方图
ax.plot_surface(x, y, hist.T, cmap='viridis')

# 设置坐标轴标签
ax.set_xlabel('Red Channel')
ax.set_ylabel('Green Channel')
ax.set_zlabel('Frequency')

# 显示图形
plt.show()

结果：

• 直方图的箱子（Bins）： 在x轴和y轴上，像素值被划分为不同的分箱（bins），每个区间代表一个像素值范围。在这个例子中，x轴和y轴各有256个分箱，分别表示0到255的像素值范围。

• Z轴的高度： Z轴上的高度（颜色越亮）表示在对应的红色和绿色通道像素值组合中出现的频率。较高的柱子表示该像素值组合在图像中出现的次数更多。

• 颜色映射： 3D直方图中的颜色映射（cmap）通常用于表示频率的大小。在这个示例中，使用的是’viridis’颜色映射，颜色从暗紫色到明黄色表示频率从低到高。

• 坐标轴标签： x轴标签是“Red Channel”（红色通道），y轴标签是“Green Channel”（绿色通道），z轴标签是“Frequency”（频率）。

2.5：多尺度直方图

将图像进行尺度变换后，得到不同尺度下的直方图，这有利于尺度不变的特征分析。

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取图像
image = cv2.imread('./5.jpg', cv2.IMREAD_GRAYSCALE)

# 定义不同尺度的比例因子
scales = [0.5, 1.0, 1.5]

# 初始化一个空的多分辨率直方图列表
multi_resolution_histograms = []

for scale in scales:
    # 缩放图像
    scaled_image = cv2.resize(image, None, fx=scale, fy=scale)
    
    # 计算直方图
    hist = cv2.calcHist([scaled_image], [0], None, [256], [0, 256])
    
    # 归一化直方图
    hist = cv2.normalize(hist, hist)
    
    # 添加到多分辨率直方图列表
    multi_resolution_histograms.append(hist)

# 显示多分辨率直方图
for i, hist in enumerate(multi_resolution_histograms):
    plt.plot(hist, label=f'Scale {scales[i]}')

plt.title('Multi-Resolution Histograms')
plt.xlabel('Pixel Value')
plt.ylabel('Frequency')
plt.legend()
plt.show()

2.6：径向直方图

径向直方图是一种描述图像中像素到中心点的距离和方向的关系方法，它可以用来表示图像中径向纹理和特征信息。径向直方图是由一系列的环形区域直方图组成的，每个环形区域直方图对应着图像中不同距离和方向下的像素分布情况。

径向直方图的计算步骤：

• 首先计算出选定图像或者区域的质心，作为起点
• 定义从质心发出投影射线的数量以及步长
• 遍历投影射线（也就是不同方向上）计算距离质心距离为r的像素值

在直方图中，横坐标是像素点距离质心的距离，竖坐标是累计像素值

# 径向直方图绘制
# 我们将使用形状的质心作为起点，绘制多个直线到形状的边界，并记录每条线上的像素值。然后，我们将这些像素值分箱，得到径向直方图。
import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取包含形状的二值图像
shape_image = cv2.imread('./5.jpg', cv2.IMREAD_GRAYSCALE)

# 寻找形状的质心
M = cv2.moments(shape_image)
cx = int(M['m10'] / M['m00'])
cy = int(M['m01'] / M['m00'])

# 定义投影的角度范围和步长
num_projections = 360  # 投影线的数量
angle_step = 360 / num_projections

# 初始化径向直方图
hist_size = 100  # 直方图的箱数
hist_range = (0, 255)  # 像素值范围
radial_histogram = np.zeros(hist_size, dtype=np.float32)

# 计算径向直方图
for angle in range(0, 360, int(angle_step)):
    angle_rad = np.deg2rad(angle)
    for r in range(hist_size): # r是像素点距离质心的距离
        x = int(cx + r * np.cos(angle_rad))
        y = int(cy + r * np.sin(angle_rad))

        if x >= 0 and x < shape_image.shape[1] and y >= 0 and y < shape_image.shape[0]:
            pixel_value = shape_image[y, x]
            radial_histogram[r] += pixel_value

# 可视化径向直方图
plt.figure(figsize=(8, 6))
plt.plot(radial_histogram)
plt.title('Radial Histogram')
plt.xlabel('Radius from Centroid')
plt.ylabel('Sum of Pixel Values')
plt.grid(True)
plt.show()

输出结果：

径向直方图一般用于具备多个环形特征或者极坐标图像下目标特征描述、检测和定位。

2.6：轮廓或边缘直方图

在计算机视觉特征度量中，关于轮廓或边缘直方图描述符主要有三种,分别是：边缘直方图（Edge Histogram， EHD）、边缘方向直方图（Edge Direction Histogram,EDH)、梯度直方图（Histogram of Oriented Gradients, HOG)。 它们在特征提取和图像识别任务中有不同的应用和特点。

它们的核心都是：对特征进行定位并且得到灰度图、计算得到各个像素点的梯度大小和边缘方向角（边缘检测）、对灰度图进行分块并且计算每一块子图的灰度边缘方向直方图特征向量。

下面展示了绘制图像边缘方向直方图并且利用计算的每个像素点的梯度信息进行了模板匹配

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 加载原始图像
image = cv2.imread('tuxiang.jpg', 0)
origin_image = image.copy()

# 计算边缘梯度
sobel_x = cv2.Sobel(image, cv2.CV_64F, 1, 0, ksize=3)
sobel_y = cv2.Sobel(image, cv2.CV_64F, 0, 1, ksize=3)
gradient_magnitude = np.sqrt(sobel_x**2 + sobel_y**2)
gradient_direction = np.arctan2(sobel_y, sobel_x)


# 将梯度幅值转换为8位无符号整数类型
gradient_magnitude = cv2.convertScaleAbs(gradient_magnitude)

# 定义边缘梯度特征的分箱和直方图参数
num_bins = 9  # 分箱数目
bin_range = (0, 2*np.pi)  # 方向范围

# 计算边缘梯度特征直方图
hist, bin_edges = np.histogram(gradient_direction, bins=num_bins, range=bin_range)

# 可视化边缘梯度特征直方图
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
plt.figure()
plt.bar(bin_centers, hist, width=2*np.pi/num_bins)
plt.xlabel('Gradient Direction')
plt.ylabel('Frequency')
plt.show()

# 加载边缘模板
template = cv2.imread('moshi2.jpg', 0)

# 计算边缘梯度
sobel_x_template = cv2.Sobel(template, cv2.CV_64F, 1, 0, ksize=3)
sobel_y_template = cv2.Sobel(template, cv2.CV_64F, 0, 1, ksize=3)
gradient_magnitude_template = np.sqrt(sobel_x_template**2 + sobel_y_template**2)
gradient_direction_template = np.arctan2(sobel_y_template, sobel_x_template)


# 将梯度幅值转换为8位无符号整数类型
gradient_magnitude_template = cv2.convertScaleAbs(gradient_magnitude_template)


# 进行模板匹配
result = cv2.matchTemplate(gradient_magnitude, gradient_magnitude_template, cv2.TM_CCOEFF_NORMED)

# 设置阈值，获取匹配结果的位置
threshold = 0.8
locations = np.where(result >= threshold)

# 在原始图像上绘制匹配结果的边框
for loc in zip(*locations[::-1]):
    cv2.rectangle(image, loc, (loc[0] + template.shape[1], loc[1] + template.shape[0]), (0, 255, 0), 2)

# 显示结果图像
plt.subplot(131)
plt.imshow(origin_image,cmap='gray')
plt.title('origin image')
plt.axis('off')
plt.subplot(132)
plt.imshow(template,cmap='gray')
plt.title('template Image')
plt.axis('off')
plt.subplot(133)
plt.imshow(image,cmap='gray')
plt.title('Result')
plt.axis('off')

plt.show()
# cv2.imshow('Result', image)
# cv2.waitKey(0)
# cv2.destroyAllWindows()

输出结果：

上面关于轮廓或边缘直方图的案例种种其实没有真正用到边缘直方图作为特征描述符。在我看到的一篇发明专利CN104680164A确实提到了用边缘方向直方图特征实现车标特征提取的方法，它还结合了前面纹理特征的Hu不变矩特征，它将这两个特征结合，形成特征向量送入神经网络进行训练。这其实启发了一点就是不一定就是只用一种方法作为特征描述符，这样有很大局限性，比如如果仅仅只用灰度直方图，那显然弱化了图像中的纹理特征，而边缘直方图恰恰就比只用灰度直方图强，并且可以结合多个特征度量方法一同来构成特征向量。

其实在写到轮廓或边缘直方图时，我觉得仅仅用它来作为特征向量在特征匹配方面有一个很大的局限性，虽然它确实能体现纹理特征，但是正是因为它和方向密切相关，导致它不具有旋转不变性！并且在学习这一特征度量方法时，我联想到了SIFT算法，它恰恰就是应用边缘方向直方图来计算得到特征的主方向，从而让特征具有旋转不变性。

但是也只是说在特征匹配上，边缘方向直方图确实不是很好，但是正是因为它的这种方向性，我们可以看到，直方图的主导方向捕捉了人的形状，特别是在躯干和腿部周围。

所以它的方向信息，在做运动估计行为识别可能大有用处！
🧑🏻‍🚀实例1、方向梯度直方图（Histogram of Oriented Gradient, HOG）

方向梯度直方图，就是基于图像的梯度和方向构建的特征描述符，梯度的幅值是边缘和角点检测的基础、梯度的方向包含着丰富的图像特征。HOG的基本思想：将图像的局部特征用梯度幅值和方向分布来描述。

基本方法：

1. 首先，将图像划分为多个单元格(cell)8x8；（这种表示方法，不仅使得特征更紧凑，而且对单个像素值的变换不敏感，能够抵抗噪声干扰）

2. 计算单元格的方向梯度直方图

3. 将2x2的cell组成一个子块block,通过滑动窗口，设置子块滑动步长blockStride来得到。
   

4. 比如上面这张图，整个图片为8x16个cell.子块滑动步长为1，则2x2个cell组成一个子块block的花，有(16-1)x(8-1)一共105个子块，每个子块有4x9=36长度的特征向量。最后将这105个block合并，得到了整个图像的特征描述符，长度为105x36为3780

HOG描述符的向量维数取决于检测图像大小和单元格大小。

HOG描述符不具有尺度和旋转不变性，但是具有良好的 几何不变性和光学不变性，适合于人体检测。

💁🏻‍♀️OpenCV中的HOGDescriptor类：
OpenCV提供了cv::HOGDescriptor类。在python中，提供了HOG类的接口函数cv.HOGDescriptor

• 函数原型

  cv.HOGDescriptor(_winSize, _blockSize, _blockStride, _cellSize, _nbins) → retval
  hog.compute(img[, _winStride=Size(), _padding=Size()]) → descriptors
  

• 参数说明

  • winSize:检测窗口的大小，形如(w,h)的元组，默认值：(64,128)
  • blockSize:子块的大小，形如(w,h)的元组，默认值(16,16)
  • blockStride:子块的滑动步长，形如(w,h)的元组，默认值为(8,8)
  • cellSize:单元格大小，将子块划分为多个单元格。形如(w,h)的元组，默认值(8,8)
  • nbins:直方图的条数，默认值是9 【方向梯度取值范围为0~180°,（这被称为“无符号”梯度，认为两个完全相反的方向是相同的）等分为nbins个区间，单元格像素的梯度方向分配到nbins】, 默认值为9
  • img:输入图像，单通道，类型CV_8U
  • winStride:窗口大小，可选项，必须是blockStride的整数倍
  • descriptors：Hog描述符，Numpy数组，数据类型CV_32F

• 函数说明

  ⑴ 计算每个单元格cell的HOG：方向梯度的取值范围0～180度，等分为nbins个区间，单元格像素的梯度方向分配到nbins个扇形区间，累加每个区间内的像素数，得到nbins位的HOG向量。
  ⑵ 构造子块block的HOG：多个单元格cell组合为子块，子块的HOG描述符就是多个单元格HOG向量的串联，长度为nbins*blockSize/cellSize。
  ⑶ 整个检测窗口的HOG：子块block以步长blockStride在检测窗口内滑动，遍历检测窗口，检测窗口的HOG就是每个子块block的HOG的串联。
  因此，检测窗口的HOG的向量维数是：

  lenHOG = nbins * (blockSize[0]/cellSize[0]) * (blockSize[1]/cellSize[1])x((winSize[0]-blockSize[0])/blockStride[0] + 1)x((winSize[1]-blockSize[1])/blockStride[1] + 1)
  

• 注意问题

  • 函数cv.HOGDescriptor实例化HOGDescriptor类，定义一个HOGDescriptor类对象。成员函数hog.compute计算给定图像的HOG描述符。

  # 构造 HOG 检测器
  winSize = (40, 40)
  blockSize = (20, 20)
  blockStride = (10, 10)
  cellSize = (10, 10)
  nbins = 8
  hog = cv.HOGDescriptor(winSize, blockSize, blockStride, cellSize, nbins)
  # hog = cv.HOGDescriptor(_winSize=(40,40), _blockSize=(20,20), 
  #      _blockStride=(10,10), _cellSize=(10,10), _nbins=8)
  
  

  • 检测窗口大小winSize是子块大小blockSize的整数倍，子块大小blockSize是单元格大小cellSize的整数倍，子块大小blockSize是滑动步长blockStride的整数倍
  • 函数中方向梯度的取值范围是0～180度，而不是0～360度。
  • cv::HOGDescriptor类的功能丰富，参数和成员函数很多，例如可以实现尺度不变性检测。更多使用方法可以参见OpenCV官方文档

• 例程：本例程示例HOG描述符的使用。为了便于解释HOG原理和绘图，例程中将检测窗口、子块和单元格设为相同的尺寸，实际应用时可以参考函数默认值来设置。

# 【1609】特征描述之 HOG 描述符
import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt

def drawHOG(image, descriptors, cx, cy, rad):
    angles = np.arange(0, 180, 22.5).astype(np.float32)  # start, stop, step
    normGrad = descriptors/np.max(descriptors).astype(np.float32)
    gx, gy = cv.polarToCart(normGrad*rad, angles, angleInDegrees=True)
    for i in range(angles.shape[0]):
        px, py = int(cx+gx[i]), int(cy+gy[i])
        cv.arrowedLine(image, (cx,cy), (px, py), 0, tipLength=0.1)  # 黑色
    return image

if __name__ == '__main__':
    # (1) 读取样本图像，构造样本图像集合
    img = cv.imread("4.jpg", flags=0)  # 灰度图像
    height, width, wCell, d = 200, 200, 20, 10
    img = cv.resize(img, (width, height))  # 调整为统一尺寸

    # (2) 构造 HOG 检测器
    winSize = (20, 20)
    blockSize = (20, 20)
    blockStride = (20, 20)
    cellSize = (20, 20)
    nbins = 8
    hog = cv.HOGDescriptor(winSize, blockSize, blockStride, cellSize, nbins)
    lenHOG = nbins * (blockSize[0]/cellSize[0]) * (blockSize[1]/cellSize[1]) \
            * ((winSize[0]-blockSize[0])/blockStride[0] + 1) \
            * ((winSize[1]-blockSize[1])/blockStride[1] + 1)
    print("length of descriptors:", lenHOG)

    # (3) 计算检测区域的 HOG 描述符
    xt, yt = 80, 80  # 检测区域位置
    cell = img[xt:xt+wCell, yt:yt+wCell]
    cellDes = hog.compute(cell)  # HOG 描述符，(8,)
    normGrad = cellDes/np.max(cellDes).astype(np.float32)
    print("shape of descriptors:{}".format(cellDes.shape))
    print(cellDes)

    # (4) 绘制方向梯度示意图
    imgGrad = cv.resize(cell, (wCell*10, wCell*10), interpolation=cv.INTER_AREA)
    Gx = cv.Sobel(img, cv.CV_32F, 1, 0, ksize=5)  # X 轴梯度 Gx
    Gy = cv.Sobel(img, cv.CV_32F, 0, 1, ksize=5)  # Y 轴梯度 Gy
    magG, angG = cv.cartToPolar(Gx, Gy, angleInDegrees=True)  # 极坐标求幅值与方向 (0~360)
    print(magG.min(), magG.max(), angG.min(), angG.max())
    angCell = angG[xt:xt+wCell, yt:yt+wCell]
    box = np.zeros((4, 2), np.int32)  # 计算旋转矩形的顶点, (4, 2)
    for i in range(wCell):
        for j in range(wCell):
            cx, cy = i*10+d, j*10+d
            rect = ((cx,cy), (8,1), angCell[i,j])  # 旋转矩形类
            box = np.int32(cv.boxPoints(rect))  # 计算旋转矩形的顶点, (4, 2)
            cv.drawContours(imgGrad, [box], 0, (0,0,0), -1)

    # (5) 绘制检测区域的方向梯度直方图
    cellHOG = np.ones((201,201), np.uint8)  # 白色
    cellHOG = drawHOG(cellHOG, cellDes, xt+d, yt+d, 40)

    # (6) 绘制图像的方向梯度直方图
    imgHOG = np.ones(img.shape, np.uint8)*255  # 白色
    for i in range(10):
        for j in range(10):
            xc, yc = 20*i, 20*j
            cell = img[xc:xc+wCell, yc:yc+wCell]
            descriptors = hog.compute(cell)  # HOG 描述符，(8,)
            imgHOG = drawHOG(imgHOG, descriptors, xc+d, yc+d, 8)
    imgWeight = cv.addWeighted(img, 0.5, imgHOG, 0.5, 0)

    plt.figure(figsize=(9, 6.2))
    plt.subplot(231), plt.title("1. Original")
    cv.rectangle(img, (xt,yt), (xt+wCell,yt+wCell), (0,0,0), 2)  # 绘制 block
    plt.axis('off'), plt.imshow(img, cmap='gray')
    plt.subplot(232), plt.title("2. Oriented gradient")
    angNorm = np.uint8(cv.normalize(angG, None, 0, 255, cv.NORM_MINMAX))
    plt.axis('off'), plt.imshow(angNorm, cmap='gray')
    plt.subplot(233), plt.title("3. Image with HOG")
    cv.rectangle(imgWeight, (xt,yt), (xt+wCell,yt+wCell), (0,0,0), 2)  # 绘制 block
    plt.axis('off'), plt.imshow(imgWeight, cmap='gray')
    plt.subplot(234), plt.title("4. Grad angle of cell")
    plt.axis('off'), plt.imshow(imgGrad, cmap='gray')
    plt.subplot(235), plt.title("5. HOG of cell")
    strAng = ("0", "22", "45", "67", "90", "112", "135", "157")
    plt.bar(strAng, cellDes*wCell*wCell)
    plt.subplot(236), plt.title("6. HOG diagram of cell")
    plt.axis('off'), plt.imshow(cellHOG, cmap='gray')
plt.tight_layout()
plt.show()

⑴ 子图1是原始图像，图中黑色方框是一个单元格cell。子图2是原始图像的梯度方向图，像素值的大小反映梯度方向的角度。
⑵ 子图4是子图1中方框位置单元格cell的梯度方向图，图中的线段表示像素点的梯度方向。注意例程中梯度方向的范围是0～180度。
⑶ 子图5是对子图4单元格中的所有像素点，按8个方向区间绘制的方向梯度直方图。子图6是子图5的单元格方向梯度直方图的空间矢量表示。
⑷ 子图3是整个图像的可视化方向梯度直方图.将图像划分为10*10个单元格，计算每个单元格的HOG，表示为如子图6的空间矢量形式。
⑸ 虽然例程给出了HOG处理过程和结果的各种图像，这是为了便于理解HOG的思路和计算步骤。在实际应用中，检测图像的HOG是维数为lenHOG的特征向量，而不是二维图像。

---

💐本节主要讲解了计算机视觉特征度量的第二类方法：统计区域度量；关于第三类基空间变换的方法，将在下篇进行讲解~
因为博主个人专业能力和知识还在提升中，若文章内容有错误或有失偏颇的地方，还请大家多多指出🍊

本文的部分内容参考和整合自：
[1]一文弄懂图像的矩和相关应用
[2]使用Hu矩进行形状匹配
[3]opencv入门：矩特征，Hu矩
[4]【youcans 的 OpenCV 例程200篇】196.图像的矩和不变矩（cv.moments）
[5]【OpenCV 例程 300篇】233. 区域特征之矩不变量
[6]第十三章 图像特征Vol.1：全局特征与区域特征
[7]【youcans 的 OpenCV 例程200篇】194.寻找图像轮廓（cv.findContours）
[8]边缘直方图描述符（EHD）Edge Histogram Descriptor
[9]openCV【实践系列】2——OpenCV方向梯度直方图
[10]基于灰度直方图和边界方向直方图的图像检索
[11]一文讲解梯度直方图Hog

---

💐在此表示感谢！