【C++ 图像处理 OpenCV】20230913-HOG特征提取

豆得儿不是猫

已于 2023-12-19 20:17:43 修改

阅读量370

点赞数

分类专栏： OpenCV C++ 文章标签：计算机视觉 opencv 人工智能

于 2023-09-14 09:46:24 首次发布

本文链接：https://blog.csdn.net/qq_37245700/article/details/132869609

版权

C++ 同时被 2 个专栏收录

9 篇文章 0 订阅

订阅专栏

OpenCV

6 篇文章 0 订阅

订阅专栏

HOG特征提取

参考博文：
1.Python提取图像HOG特征-CSDN博客
2.图像之HOG特征描述算子-行人检测

一、步骤

1）灰度化

将图像看做一个x,y,z（灰度）的三维图像，python读取图像时直接灰度处理

img = cv2.imread('C:/zhaokai/software/DataSet/images/1.png', cv2.IMREAD_GRAYSCALE)

2）Gamma校正

采用Gamma校正法对输入图像进行颜色空间的标准化（归一化）；目的是调节图像的对比度，降低图像局部的阴影和光照变化所造成的影响，同时可以抑制噪音的干扰；
$outImg = inImg^{gamma}$
简单一点的就是，直接实现上述公式。这种方式适用于图片不多，手动调节gamma参数。如下

def gamma(inImg):
    outImg = np.power( inImg / 255.0, 2.2 )
    return outImg

gamma一般取值【2.2，2.4】。

当复杂一些的情况下，比如需要批量传入图像时，

def gamma_batch(img):
    mean = np.mean(img)
    gamma_val = math.log10(0.5) / math.log10(mean / 255)  # 公式计算gamma
    gamma_table = [np.power(x / 255.0, gamma_val) * 255.0 for x in range(256)]  # 建立映射表
    gamma_table = np.round(np.array(gamma_table)).astype(np.uint8)  # 颜色值为整数
    gamma_val =cv2.LUT(img, gamma_table)
    return gamma_val  # 图片颜色查表。另外可以根据光强（颜色）均匀化原则设计自适应算法。

图片多的情况下，手动设置gamma值显得过于麻烦，这时候就需要采用自动Gamma矫正，将RGB图片转成灰度图（上述方法中没有进行转换灰度图，需要自己转换），计算灰度图的数据均值，通过下面的计算公式计算gamma值。

gamma_val = math.log10(0.5) / math.log10(mean / 255)

3）计算图像每个像素的梯度（包括大小和方向）

主要是为了捕获轮廓信息，同时进一步弱化光照的干扰。计算图像横坐标和纵坐标方向的梯度，并据此计算每个像素位置的梯度方向值；求导操作不仅能够捕获轮廓，人影和一些纹理信息，还能进一步弱化光照的影响。

公式如下图

在这里插入图片描述

在python中如下

gradient_values_x = cv2.Sobel(img, cv2.CV_64F, 1, 0, ksize=5)  # x方向梯度
gradient_values_y = cv2.Sobel(img, cv2.CV_64F, 0, 1, ksize=5)  # y方向梯度
gradient_magnitude = np.sqrt(np.power(gradient_values_x, 2) + np.power(gradient_values_y, 2))#梯度幅值
gradient_angle = np.arctan2(gradient_values_x, gradient_values_y)#梯度方向

做完上述步骤后，还需要进行一些修改，是为了之后的获取直方图做准备。现在gradient_angle的取值范围是
$gradient\_angle \in [\frac{-\pi}{2},\frac{\pi}{2}]$
需要将其调整到
$gradient\_angle \in [0,{\pi}]$
所以增加下面的代码：

# 角度转换至（0-180）
    gradient_angle[gradient_angle > 0] *= 180 / 3.14
    gradient_angle[gradient_angle < 0] = (gradient_angle[gradient_angle < 0] + 3.14) * 180 / 3.14

经过处理后的图像区别如下。

在这里插入图片描述

4）将图像划分成小cells

例如8*8像素/cell；

# 获取梯度值cell图像，梯度方向cell图像
def div(img, cell_x, cell_y, cell_w):
    # 创建的是一个四维矩阵，cell[i][j]表示第i行第j个cell的像素矩阵值，每个cell为8*8的矩阵
    cell = np.zeros(shape=(cell_x, cell_y, cell_w, cell_w))
    # axis = 0 横向切分img 平均切分cell_x次。
    img_x = np.split(img, cell_x, axis=0)
    for i in range(cell_x):
        # axis = 1 纵向切分img_x[i] 平均切分cell_y次。
        img_y = np.split(img_x[i], cell_y, axis=1)
        # print("len(img_y) == ",len(img_y))
        for j in range(cell_y):
            cell[i][j] = img_y[j]
    return cell

5）统计每个cell的梯度直方图

（不同梯度的个数），即可形成每个cell的descriptor；

def get_bins(grad_cell, ang_cell):
    bins = np.zeros(shape=(grad_cell.shape[0], grad_cell.shape[1], 9))
    for i in range(grad_cell.shape[0]):
        for j in range(grad_cell.shape[1]):
            binn = np.zeros(9) # 9个柱的直方图，统计9个方向的梯度幅值总和 
            grad_list = np.int8(grad_cell[i, j].flatten())  # 每个cell中的64个梯度值展平，并转为整数
            ang_list = ang_cell[i, j].flatten()  # 每个cell中的64个梯度方向展平)
            ang_list = np.int8(ang_list / 20.0)  # 0-9(梯度方向（0-180°）除20，归化到0-9)
            ang_list[ang_list >= 9] = 0
            for m in range(len(ang_list)):
                binn[ang_list[m]] += int(grad_list[m])  # 不同角度对应的梯度值相加，为直方图的幅值

            # 每个cell的梯度方向直方图可视化
            N = 9
            x = np.arange( N )
            str1 = ( '0-20', '20-40', '40-60', '60-80', '80-100', '100-120', '120-140', '140-160', '160-180' )
            plt.bar( x, height = binn, width = 0.8, label = 'cell histogram', tick_label = str1 )
            for a, b in zip(x, binn):
                plt.text( a, b+0.05, '{}'.format(b), ha = 'center', va = 'bottom', fontsize = 10 )
            plt.show()

            bins[i][j] = binn
    return bins

其中一个cell的结果如下图所示
在这里插入图片描述

6）将每几个cell组成一个block

（例如2*2个cell/block），一个block内所有cell的特征descriptor串联起来便得到该block的HOG特征descriptor。

feature = []
    # 2*2个cell组成一个block，并计算特征值，这里定义的特征值计算方式为feature = x-E(x)
    for i in range(cell_x - 1):
        for j in range(cell_y - 1):
            tmp = []
            tmp.append(bins[i, j])
            tmp.append(bins[i + 1, j])
            tmp.append(bins[i, j + 1])
            tmp.append(bins[i + 1, j + 1])
            tmp -= np.mean(tmp)

7）获取features特征

将图像image内的所有block的HOG特征descriptor串联起来就可以得到该image（你要检测的目标）的HOG特征descriptor了。这个就是最终的可供分类使用的特征向量了。总共的特征个数为：
$Total_{features} = cellnums_{block}*bins.size*(cell\_x-block.x+1)*(cell\_y-block.y+1)$
其中

cellnums_block表示一个block中的cell个数,为4
binn.size为直方图个列数，为9。
cell_x为80，cell_y为160，
block.x == block.y == 2

所以本本文测试的最后特征个数为4*9*79*159=452196个。

feature.append(tmp.flatten())

二、完整代码

# coding:utf-8
# *********************************************************************************************************
'''
说明：利用python/numpy/opencv实现图像HOG特征的提取
算法思路：
算法思路:
        1)以灰度图的方式加载图片，resize到（128,64）;
        2）灰度图像gamma校正;
		3)利用一阶微分算子Sobel函数，分别计算出灰度图像X方向和Y方向上的一阶微分/梯度图像，根据得到的两幅
        梯度图像(X方向上的梯度图像和Y方向上的梯度图像)，计算出这两幅梯度图像所对应的梯度幅值图像gradient_magnitude、
        梯度方向图像gradient_angle
		4)构造(cell_x = 128/8 =16, cell_y= 64/8 =8)大小的cell图像----梯度幅值的grad_cell图像，梯度方向的ang_cell图像，
        每个cell包含有8*8 = 64个值；
		5)将每个cell根据角度值（0-180）分为9个bin，并计算每个cell中的梯度方向直方图,每个cell有9个值；
		6)每（2*2）个cell为一个block，总共15*7个block,计算每个block的梯度方向直方图，并进行归一化处理，每个block中有9*4=36个值；
		7)计算整幅图像的梯度方向直方图HOG:将计算出来的所有的Block的HOG梯度方向直方图的特征向量首尾相接组成一个维度很大的向量
        长度为：15*7*36 = 3780，
        这个特征向量就是整幅图像的梯度方向直方图特征，这个特征可用于SVM分类。
'''
import math

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


# 灰度图像gamma校正
def gamma(img):
    gamma_val = np.power( img / 255.0, 2.2 )
    return gamma_val

def gamma_batch(img):
    mean = np.mean(img)
    gamma_val = math.log10(0.5) / math.log10(mean / 255)  # 公式计算gamma
    gamma_table = [np.power(x / 255.0, gamma_val) * 255.0 for x in range(256)]  # 建立映射表
    gamma_table = np.round(np.array(gamma_table)).astype(np.uint8)  # 颜色值为整数
    gamma_val =cv2.LUT(img, gamma_table)
    return gamma_val  # 图片颜色查表。另外可以根据光强（颜色）均匀化原则设计自适应算法。




# 获取梯度值cell图像，梯度方向cell图像
def div(img, cell_x, cell_y, cell_w):
    # 创建的是一个四维矩阵，cell[i][j]表示第i行第j个cell的像素矩阵值，每个cell为8*8的矩阵
    cell = np.zeros(shape=(cell_x, cell_y, cell_w, cell_w))
    # axis = 0 横向切分img 平均切分cell_x次。
    img_x = np.split(img, cell_x, axis=0)
    for i in range(cell_x):
        # axis = 1 纵向切分img_x[i] 平均切分cell_y次。
        img_y = np.split(img_x[i], cell_y, axis=1)
        # print("len(img_y) == ",len(img_y))
        for j in range(cell_y):
            cell[i][j] = img_y[j]
    return cell


# 获取梯度方向直方图图像，每个像素点有9个值
def get_bins(grad_cell, ang_cell):
    bins = np.zeros(shape=(grad_cell.shape[0], grad_cell.shape[1], 9))
    for i in range(grad_cell.shape[0]):
        for j in range(grad_cell.shape[1]):
            binn = np.zeros(9)
            grad_list = np.int8(grad_cell[i, j].flatten())  # 每个cell中的64个梯度值展平，并转为整数
            ang_list = ang_cell[i, j].flatten()  # 每个cell中的64个梯度方向展平)
            ang_list = np.int8(ang_list / 20.0)  # 0-9
            ang_list[ang_list >= 9] = 0
            for m in range(len(ang_list)):
                binn[ang_list[m]] += int(grad_list[m])  # 不同角度对应的梯度值相加，为直方图的幅值

            # 每个cell的梯度方向直方图可视化
            N = 9
            x = np.arange( N )
            str1 = ( '0-20', '20-40', '40-60', '60-80', '80-100', '100-120', '120-140', '140-160', '160-180' )
            plt.bar( x, height = binn, width = 0.8, label = 'cell histogram', tick_label = str1 )
            for a, b in zip(x, binn):
                plt.text( a, b+0.05, '{}'.format(b), ha = 'center', va = 'bottom', fontsize = 10 )
            plt.show()

            bins[i][j] = binn
    return bins


# 计算图像HOG特征向量
def hog(img, cell_x, cell_y, cell_w):
    height, width = img.shape
    gradient_values_x = cv2.Sobel(img, cv2.CV_64F, 1, 0, ksize=5)  # x方向梯度
    gradient_values_y = cv2.Sobel(img, cv2.CV_64F, 0, 1, ksize=5)  # y方向梯度
    gradient_magnitude = np.sqrt(np.power(gradient_values_x, 2) + np.power(gradient_values_y, 2))
    gradient_angle = np.arctan2(gradient_values_x, gradient_values_y)
    # plt.figure()
    # plt.subplot(3,2,1)
    # plt.title("gradient_values_x")
    # plt.imshow(gradient_values_x)
    # plt.subplot(3,2,2)
    # plt.title("gradient_values_y")
    # plt.imshow(gradient_values_y)
    # plt.subplot(3,2,3)
    # plt.title("gradient_magnitude")
    # plt.imshow(gradient_magnitude)
    # plt.subplot(3,2,4)
    # plt.title("gradient_angle")
    # plt.imshow(gradient_angle)
    print(gradient_magnitude.shape, gradient_angle.shape)
    # 角度转换至（0-180）
    gradient_angle[gradient_angle > 0] *= 180 / 3.14
    gradient_angle[gradient_angle < 0] = (gradient_angle[gradient_angle < 0] + 3.14) * 180 / 3.14
    # plt.subplot( 3, 2, 5 )
    # plt.title("gradient_angle transform")
    # plt.imshow( gradient_angle )
    # plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.5, hspace=0.5)
    # plt.show()

    # 获取梯度幅值的cell
    grad_cell = div(gradient_magnitude, cell_x, cell_y, cell_w)
    # 获取梯度角度的cell
    ang_cell = div(gradient_angle, cell_x, cell_y, cell_w)
    bins = get_bins(grad_cell, ang_cell)
    feature = []
    # 2*2个cell组成一个block，并计算特征值，这里定义的特征值计算方式为feature = x-E(x)
    for i in range(cell_x - 1):
        for j in range(cell_y - 1):
            tmp = []
            tmp.append(bins[i, j])
            tmp.append(bins[i + 1, j])
            tmp.append(bins[i, j + 1])
            tmp.append(bins[i + 1, j + 1])
            tmp -= np.mean(tmp)
            feature.append(tmp.flatten())
    return np.array(feature).flatten()


if __name__ == '__main__':
    # 修改为自己的图片路径
    img = cv2.imread('your own image path', cv2.IMREAD_GRAYSCALE)
    if (img is None):
        print('Not read image.')
    print(img.shape)
    resizeimg = cv2.resize(img, (1280, 640), interpolation=cv2.INTER_CUBIC)
    cell_w = 8
    cell_x = int(resizeimg.shape[0] / cell_w)  # cell行数
    cell_y = int(resizeimg.shape[1] / cell_w)  # cell列数
    print('The size of cellmap is {}*{} '.format(cell_x, cell_y))
    # gammaimg = gamma(resizeimg) * 255
    gammaimg = gamma_batch(resizeimg) * 255
    print("cell_x == ",cell_x,"cell_y == ", cell_y)
    feature = hog(gammaimg, cell_x, cell_y, cell_w)
    print(feature.shape)

三、和pytorch提供的比较

3.1 比较

from skimage.feature import hog, local_binary_pattern
import cv2
ppc = (8, 8)
cpb = (2, 2)


def get_hog_feature(img):
    res = hog(img, orientations=9, pixels_per_cell=ppc, cells_per_block=cpb, block_norm='L1', transform_sqrt=False, feature_vector=True)
    print("get_hog_feature == ", res)
    return res
img = cv2.imread(r'your own image path', cv2.IMREAD_GRAYSCALE)
if (img is None):
    print('Not read image.')
cv2.imshow("img",img)
print(img.shape)
resizeimg = cv2.resize(img, (64, 128), interpolation=cv2.INTER_CUBIC)
cv2.imshow("resizeimg",resizeimg)
print(resizeimg.shape)

res = get_hog_feature(resizeimg)
print("res.shape() == ", res.size)
cv2.waitKey(0)

结果展示为：
在这里插入图片描述

自己手动计算HOG的结果为
在这里插入图片描述

可以看到最后的feature的大小是一致的。feature值的不同是因为自己的没有做归一化处理。

使用L2-norm，对feature归一化处理后，

在这里插入图片描述

不同的归一化方式，其特征值会有所不一样。

3.2 替换代码

替换第一节中第6）步骤的代码，增加了对block特征值的归一化处理

def hog(img, cell_x, cell_y, cell_w):
    height, width = img.shape
    gradient_values_x = cv2.Sobel(img, cv2.CV_64F, 1, 0, ksize=5)  # x方向梯度
    gradient_values_y = cv2.Sobel(img, cv2.CV_64F, 0, 1, ksize=5)  # y方向梯度
    gradient_magnitude = np.sqrt(np.power(gradient_values_x, 2) + np.power(gradient_values_y, 2))
    gradient_angle = np.arctan2(gradient_values_x, gradient_values_y)
    plt.figure()
    plt.subplot(3,2,1)
    plt.title("gradient_values_x")
    plt.imshow(gradient_values_x)
    plt.subplot(3,2,2)
    plt.title("gradient_values_y")
    plt.imshow(gradient_values_y)
    plt.subplot(3,2,3)
    plt.title("gradient_magnitude")
    plt.imshow(gradient_magnitude)
    plt.subplot(3,2,4)
    plt.title("gradient_angle")
    plt.imshow(gradient_angle)
    print(gradient_magnitude.shape, gradient_angle.shape)
    # 角度转换至（0-180）
    gradient_angle[gradient_angle > 0] *= 180 / 3.14
    gradient_angle[gradient_angle < 0] = (gradient_angle[gradient_angle < 0] + 3.14) * 180 / 3.14
    plt.subplot( 3, 2, 5 )
    plt.title("gradient_angle transform")
    plt.imshow( gradient_angle )
    plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.5, hspace=0.5)
    plt.show()

    # 获取梯度幅值的cell
    grad_cell = div(gradient_magnitude, cell_x, cell_y, cell_w)
    # 获取梯度角度的cell
    ang_cell = div(gradient_angle, cell_x, cell_y, cell_w)
    bins = get_bins(grad_cell, ang_cell)
    feature = []
    # 2*2个cell组成一个block，并计算特征值，这里定义的特征值计算方式为feature = x-E(x)
    for i in range(cell_x - 1):
        for j in range(cell_y - 1):
            tmp = []
            tmp.append(bins[i, j])
            tmp.append(bins[i + 1, j])
            tmp.append(bins[i, j + 1])
            tmp.append(bins[i + 1, j + 1])
            tmp -= np.mean(tmp)
            feature.append(tmp.flatten())
    # # 3*3个cell组成一个block，并计算特征值，这里定义的特征值计算方式为feature = x-E(x)
    # for i in range(cell_x - 2):
    #     for j in range(cell_y - 2):
    #         tmp = []
    #         tmp.append(bins[i, j])
    #         tmp.append(bins[i + 1, j])
    #         tmp.append(bins[i + 2, j])
    #         tmp.append(bins[i, j+1])
    #         tmp.append(bins[i + 1, j+1])
    #         tmp.append(bins[i + 2, j+1])
    #         tmp.append(bins[i, j])
    #         tmp.append(bins[i + 1, j+2])
    #         tmp.append(bins[i + 2, j+2])
    #         tmp -= np.mean(tmp)
    #         feature.append(tmp.flatten())
    
    # 增加
    res = np.array(feature).flatten()
    tensor = torch.from_numpy(res)
    re = torch.nn.functional.normalize(tensor, p=2, dim=-1)
    print("tensor == ",re)
    return re

#         tmp.append(bins[i + 2, j+1])
#         tmp.append(bins[i, j])
#         tmp.append(bins[i + 1, j+2])
#         tmp.append(bins[i + 2, j+2])
#         tmp -= np.mean(tmp)
#         feature.append(tmp.flatten())

# 增加
res = np.array(feature).flatten()
tensor = torch.from_numpy(res)
re = torch.nn.functional.normalize(tensor, p=2, dim=-1)
print("tensor == ",re)
return re