图像哈希2：基于环形分区和NMF的鲁棒感知图像哈希(CCF A),复现论文

图像哈希系列：
图像哈希1：基于四元数离散傅里叶变换和对数极坐标的鲁棒图像哈希算法
图像哈希2：基于环形分区和NMF的鲁棒感知图像哈希(CCF A)
图像哈希3：基于四元数离散余弦变换的鲁棒感知图像哈希研究（性能优于近年一些较好的研究）
图像哈希4：基于四元数 SVD 的奇异值的鲁棒图像哈希
未完待续…
欢迎交流# 基于环形分区和NMF的鲁棒感知图像哈希

论文：Robust Perceptual Image Hashing Based on Ring Partition and NMF, 论文链接
doi:10.1109/TKDE.2013.45
摘要：本文利用环形分割和非负矩阵因式分解（NMF）设计了一种高效的图像散列方法，它既具有旋转鲁棒性，又具有良好的判别能力。其主要贡献在于首次在图像散列中使用了旋转不变二次图像的新构造，这有助于使图像散列具有抗旋转性。此外，NMF 系数通过内容保护操作近似线性变化，从而用相关系数衡量哈希相似性。用 346 幅图像进行了实验，以说明其效率。实验结果表明，所提出的哈希算法对图像旋转、JPEG 压缩、水印嵌入、高斯低通滤波、伽玛校正、亮度调整、对比度调整和图像缩放等内容保护操作具有鲁棒性。此外，还将哈希算法与最先进的算法进行了ROC比较，结果表明，在鲁棒性和区分度方面，所提出的哈希算法在分类性能上远远优于所有这些算法。
关键词： 鲁棒图像哈希，多媒体安全，非负矩阵分解，环形分区

介绍

互联网和多媒体设备，如数码相机、扫描仪和智能手机，使我们能够更方便、更快捷地捕捉、存储、分享和传递成千上万的图像、歌曲和视频。在我们越来越享受多媒体产品带来的生活的同时，学术界和产业界仍然面临着许多具有挑战性的问题。例如，在个人电脑中，一幅图像可能有多个副本，这些副本可以是不同的数字表现形式，但视觉内容与原始图像相同。如何从大规模多媒体数据中有效搜索图像的所有相似版本（包括原始版本及其副本），显然是一个现实而又具有挑战性的问题。特别是，虽然强大的多媒体工具使得编辑和篡改等数字操作比以往任何时候都更加容易，但确保多媒体内容的安全性一直是多媒体界关注的重要问题。这些现实需求催生了一种新兴的多媒体技术，即图像哈希。在本文中研究了一种新的稳健图像哈希技术。

图像散列不仅能让我们在大型数据库中快速找到图像副本，还能确保数字图像的内容安全。图像哈希将输入图像映射成一个短字符串，称为图像哈希，已被广泛应用于图像检索、图像认证、数字水印、图像副本检测、篡改检测、图像索引、多媒体取证和还原参考图像质量评估等领域。

有一些经典的加密哈希函数，如 SHA-1 和 MD5，可以将输入信息（文档、图像、图形、文本、视频等）转换为固定大小的字符串。但是，它们对比特级变化很敏感，不适合图像哈希。这是因为在实际应用中，数字图像通常会经过正常的数字处理，如 JPEG 压缩和图像增强，而不会改变图像的视觉内容。这是图像散列函数两个基本特性中的第一个特性，即感知鲁棒性，也就是散列函数应该对几何变换、格式转换和 JPEG 压缩等保留内容的操作具有鲁棒性。换句话说，图像的哈希值和处理过的版本应该是相同或非常相似的。只有当视觉内容被恶意操作（如删除对象和插入对象）改变时，哈希值才会发生重大变化。另一个基本特性是鉴别能力，即不同内容的图像应具有不同的哈希值。这意味着不同图像之间的哈希值距离应该足够大。此外，图像哈希函数在应用于某些应用时还应满足其他属性。例如，应用于数字取证时，图像哈希值必须依赖于一个或多个密钥。

因此，近十年来有许多图像散列算法被成功设计出来。从应用的角度来看，图像哈希仍存在一些局限性。例如，图像旋转是一种有用的操作。在照片的后期处理中，人们经常利用图像旋转来修正地面不平整的照片。然而，大多数算法对旋转都很敏感。这意味着它们会错误地将这些经过校正的照片识别为不同的图像。有些方法对旋转具有鲁棒性，但它们的鉴别能力不够好。在开发高性能哈希算法时，如何同时满足旋转鲁棒性和良好判别能力的要求是一项具有挑战性的任务。在本文中，提出了一种高效的鲁棒图像哈希算法，以同时满足旋转鲁棒性和良好的判别能力。该算法基于环形分割和非负矩阵因式分解（NMF）。其关键技术是通过环形分割来构建新颖的二次图像。二次图像不受图像旋转的影响，从而使算法具有抗旋转能力。此外，NMF 的使用还为所提出的算法提供了良好的判别能力。这是因为 NMF 是一种高效的基于部件表示的学习技术，而图像哈希包含的局部图像内容信息越多，图像哈希的判别能力就越强。用 346 幅图像（包括 USC-SIPI 图像数据库中的 146 幅图像）和 200 幅不同颜色的图像（其中 100 幅图像来自地面实况数据库，67 幅图像从互联网下载，33 幅图像由数码相机拍摄）进行了实验，以说明该算法的效率。实验证明，所提出的算法在旋转鲁棒性和判别能力之间达到了理想的平衡。

一般来说，图像哈希的构建基于三个基本步骤，即预处理、图像特征提取和哈希的构建。其中，最关键的当然是图像特征提取步骤。

提出的图像哈希算法

NMF

NMF 是一种高效的降维技术，与主成分分析（PCA）和矢量量化（VQ）相比，它在学习基于部件的表示时表现出更好的性能。事实上，NMF 已成功应用于人脸识别、图像表示、图像分析、信号分离、数据聚类等领域。一个非负矩阵$V=(V_{ij}{)}_{M\times N}$通常被视为大小为$M\times 1$的$N$个向量的组合。V 的 NMF 结果是两个非负矩阵因子，即$B=(B_{ij}{)}_{(}K\times N)$和$C=(C_{ij}{)}_{K\times N}$，其中 K 是 NMF 的秩，且$K<min(M,N)$。矩阵B和C分别称为基矩阵和系数矩阵（或编码矩阵）。它们可以用来近似表示 V，即
$$V\approx BC$$
许多文献提出了各种算法来实现NMF,具体方法请看论文，这里就不具体展开了。

环形分区

一般来说，旋转操作以图像中心为坐标原点。下图(a)显示了图像的中心部分，下图(b)显示了旋转图像的相应部分。显然，(a)和(b)对应环中的视觉内容在图像旋转后保持不变。因此，为了使图像哈希对旋转具有鲁棒性，我们可以将图像划分为不同的环，并使用它们来形成对旋转不变的二次图像。下图为二次图像构造示意图，其中(a)为正方形图像，分为7个环，(b)为这些环形成的二次图像。二次图像构造的详细方案如下：

设正方形图像的大小为$m\times m$ ；n是环的总数，$R_k$是第k个环中的一组像素值（k=1,2,…,n）。在这里，我们只使用正方形图像的内切圆中的像素，并将内切圆划分为面积相等的环，这是因为每个环都应该是二次图像的一列。环划分可以通过计算圆半径和每个像素与图像中心之间的距离来完成。除最内层环外，每个环的像素可以由两个相邻半径决定。假设$r_k$是第 k 个半径（k=1,2,…,n），从小值标记到大值。因此，$r_1$和 $r_n$分别是最内层和最外层圆的半径。显然，$r_n=\lfloor m/2\rfloor$对于$m\times m$图像，其中$\lfloor \cdot \rfloor$表示向下舍入。为了确定其他半径，首先计算内切圆A的面积和每个环${\mu }_A$的平均面积如下:
$$A=\pi r_n^2,\mu =\lfloor A/n\rfloor$$
$r_1$可以通过以下方式计算。
$$r_1=\sqrt{\frac{{\mu }_A}{\pi }}$$
因此其他半径$r_k(k=2,3,...,n-1)$可以通过以下等式获取：
$$r_k=\sqrt{\frac{{\mu }_A+\pi r_{k-1}^2}{\pi }}$$
令 p(x,y)为图像第 y 行和第 x 列中的像素值$1\le x,y\le m$。假设$(x_c,y_c)$是图像中心的坐标。因此，如果 m 是偶数，$x_c=m/2+0.5$和$y_c=m/2+0.5$ 。否则，$x_c=(m+1)/2$和$y_c=(m+1)/2$。因此p(x,y)与图像中心$(x_c,y_c)$之间的距离可以用欧氏距离来衡量，如下所示:
$$d_{x,y}=\sqrt{(x-x_c{)}^2+(y-y_c{)}^2}$$
在得到圆半径和像素距离后，我们可以将这些像素值分为n个集合，如下所示:
R 1 = { p ( x , y ) ∣ d x , y ⩽ r 1 } , R k = { p ( x , y ) ∣ r k − 1 < d x , y ⩽ r k } ( k = 2 , 3 , … , n ) \begin{array}{c} R_{1}=\left\{p(x, y) \mid d_{x, y} \leqslant r_{1}\right\}, \\ R_{k}=\left\{p(x, y) \mid r_{k-1}<d_{x, y} \leqslant r_{k}\right\}(k=2,3, \ldots, n) \end{array} R1​={p(x,y)∣dx,y​⩽r1​},Rk​={p(x,y)∣rk−1​<dx,y​⩽rk​}(k=2,3,…,n)​
接下来，我们重新排列$R_k(k=1,2,...,n)$的元素，以升序排列除一个排列序列 ，这一操作保证了$u_k$与旋转无关。由于像素坐标是离散的，每组像素的个数并不总是等于${\mu }_A$。由于每个环的像素都要构成二次图像的一列，因此要通过线性插值将${\mu }_k$映射到大小为${\mu }_A\times 1$的新向量$v_k$中。因此，二次图像 V 是通过对这些新向量进行如下排列得到的：
$$V=[v_1,v_2,...,v_n]$$
由于$v_k$与旋转无关，V对该操作也是不变的。除了旋转不变的优点外，二次图像还有一个优点，即它的列比原始图像少。具体可以参考原论文。

阐述所提出的方法

下面为提出的图像哈希具体步骤：
1.Preprocessing.输入图像首先通过双线性插值映射到归一化大小$m\times m$。这确保了我们的哈希序列对缩放操作具有弹性，并且不同大小的图像哈希具有相同的长度。对于 RGB 彩色图像将其转换为 YCbCr 颜色空间。其中 R、G 和 B 表示像素的红色、绿色和蓝色分量，Y 、Cb 和 Cr 是亮度、蓝差色度和色差色度。在颜色空间转换之后，我们采用亮度分量Y来表示。
2.Ring partition(环形分区).将Y分成n个环，并利用利用它们通过上一节描述的方案产生二次图像V。这一步的目的是构造一个旋转不变矩阵进行降维。
3.NMF.对 V 应用 NMF，然后得到系数矩阵 C。将矩阵条目合并，就能得到一个紧凑的图像哈希值。因此，哈希长度为 L=nK，其中 n 是环数，K 是 NMF的秩。
为了衡量两个图像哈希之间的相似性，我们将相关系数作为度量。设$h^{(1)}=[h_1^{(1)},h_2^{(1)},...,h_L^{(1)}]$ 和$h^{(2)}=[h_1^{(2)},h_2^{(2)},...,h_L^{(2)}]$是两个图像散列。因此，相关系数定义为:
$$S=\frac{{\sum }_{i=1}^L[h_i^{(1)}-{\mu }_1][h_i^{(2)}-{\mu }_2]}{\sqrt{{\sum }_{i=1}^L[h_i^{(1)}-{\mu }_1{]}^2}\sqrt{{\sum }_{i=1}^L[h_i^{(2)}-{\mu }_2{]}^2}+\epsilon }$$
输入图像越相似，S的值越大。如果 S 大于预定义的阈值 T，则图像被视为视觉上相同的图像。否则，它们是不同的图像，或者一个是另一个被篡改的版本。选择相关系数作为度量的原因是基于观察到 NMF 系数通过内容保存操作近似线性变化的观察结果。

实验参数

这里复现使用的参数基本和论文使用的一致，具体内容可以看论文和后面的复现代码。

实验结果

这里使用的数据集不是论文的数据集。随着互联网的发展，一张图片或许会受到多种不同的攻击，我们利用柯达数据库来评估我们的方案的鲁棒性，VOC2007数据库用于检查我们的方案的唯一性。其中柯达数据集有24张彩色图像。对于每个图像，有100个相似的版本，由16种多重攻击类型组成，包括旋转图像变形(IDR)(角度:30°)+椒盐噪声(SPN)(密度:0.01)+BA、IDR+SPN+对比度调整(CA)、IDR+SPN+Gamma校正(GC)、IDR+SPN+高斯低通滤波(GLF)、IDR+SPN+JC等。在实验中，相似图像总数为2400对。总之，使用了2400+24 =2424 张图像。VOC2007数据库中有 5,011张彩色图像，生成该数据库中每对哈希码之间的相关系数。因此，距离的总数可以由组合数C(5011，2)=5011(5011−1)/2=11255255来建立。

$$柯达数据集(多种攻击)$$
利用(ROC)对不同参数下的鲁棒性和判别能力进行可视化分类性能。因此，真阳性率 (TPR) $P_{TPR}$和假阳性率 (FPR) $P_{FPR}$都可以计算如下：
$$\begin{gathered} P_{TPR}=\frac{n_1}{N_1}, \\ {P}_{TPR}=\frac{n_2}{N_2} \end{gathered}$$
其中 n1 是被视为相似图像的视觉上相同图像对的数量，N1 是视觉相同图像的总对，n2 是被视为相似图像的不同图像对的数量，N2 是不同图像的总对。显然，TPR 和 FPR 分别表示感知鲁棒性和判别能力。实验数据集的ROC曲线如下：

主要代码

import numpy as np    #1.20.3
import cv2
import matplotlib.pyplot as plt
import os
import pandas as pd
import sklearn.decomposition
from scipy.interpolate import interp1d
from sklearn.decomposition import NMF
import imageio

def readImg(im_fn):
    im = cv2.imread(im_fn)
    if im is None:
        # print('{} cv2.imread failed'.format(im_fn))
        tmp = imageio.mimread(im_fn)
        if tmp is not None:
            imt = np.array(tmp)
            imt = imt[0]
            im = imt[:, :, 0:3]
        im=cv2.cvtColor(im, cv2.COLOR_RGB2BGR)
    return im
def read_images_from_folder(folder_path):
    """
    读取文件夹中所有的图片
    """
    image_list = []
    image_paths = []  # 用于保存图片的文件路径
    # 遍历文件夹中的所有文件
    for filename in os.listdir(folder_path):
        file_path = os.path.join(folder_path, filename)
        # 检查文件是否为图片
        if file_path.lower().endswith(('.png', '.jpg', '.jpeg', '.gif', '.bmp','.tif')):
            try:
                # 使用OpenCV读取图片
                image = readImg(file_path)
                image_list.append(image)
                image_paths.append(file_path)
            except Exception as e:
                # 处理异常，例如无法读取的文件
                print(f"Error reading image {file_path}: {e}")

    return image_list, image_paths


def divide_into_rings(image, n):
    """
    image为彩色图像,n为划分环的个数,row为图像大小
    """
    ############计算每个环形分区的半径############
    row, col= image.shape
    Rn = row // 2
    A = np.pi * Rn ** 2
    u = A // n
    R1 = np.sqrt(u / np.pi)
    R = [R1, ]
    for k in range(1, n):
        Rk = np.sqrt((u + np.pi * R[k - 1] ** 2) / np.pi)
        R.append(Rk)
    # print(f"图片大小为：{row, col}")
    # print(f"内切圆的面积A={A},每个环的面积u={u},每个环对应的半径{R}")

    # 计算图像中心
    if row % 2 == 0:
        x_center, y_center = (row / 2 + 0.5 - 1, col / 2 + 0.5 - 1)
    else:
        x_center, y_center = ((row + 1) / 2 - 1, (col + 1) / 2) - 1
    # print(f"原点为{x_center, y_center}")

    ######划分环#######
    # 计算每个像素到中心的距离
    y, x = np.ogrid[:col, :row]
    distances = np.sqrt((x - x_center) ** 2 + (y - y_center) ** 2)
    # 初始化每个环的像素集合
    rings = [set() for _ in range(n)]
    # 划分像素到每个环
    # 对于R[0]的情况
    mask = (0 < distances) & (distances <= R[0])
    rings[0] = set(zip(np.where(mask)[0], np.where(mask)[1]))

    for k in range(1, n):
        mask = (R[k - 1] < distances) & (distances <= R[k])
        rings[k] = set(zip(np.where(mask)[0], np.where(mask)[1]))
    # u=int(u//n)*n
    # 初始化矩阵结果
    V = np.zeros((int(u), n), dtype=np.uint8)
    # 对每个环进行线性插值
    for i, ring in enumerate(rings):
        pixel_values = [image[x, y] for x, y in ring]
        sort_pixel_values = sorted(pixel_values)
        size = len(sort_pixel_values)
        target_size = int(u)
        # 线性插值
        # interpolated_values = np.interp(np.linspace(0, 1, target_size),
        #                                 np.linspace(0, 1, size), sort_pixel_values)
        #
        #
        # 插值函数
        interpolator = interp1d(np.linspace(0, 1, size), sort_pixel_values, kind='linear', fill_value='extrapolate')

        new_positions = np.linspace(0, 1, target_size)
        # 使用线性插值
        interpolated_values = interpolator(new_positions)
        V[:, i] = interpolated_values

        # 小于target_size在序列后面用0来填补，大于则截取
        # if size<target_size:
        #     padding_size=target_size-size
        #     V[:, i, channel]=np.concatenate([np.array(sort_pixel_values),np.zeros(padding_size)])
        # else:
        #     V[:, i, channel]=sort_pixel_values[:target_size]
    return V

def preprocessing(image,M=512,m=32):
    """
    :param image:输入BGR彩色图片
    :param M: 调整的大小（M×M）
    :param m:环形分区环的个数
    :return:
    """
    #1.使用双线性插值将图片大小调整为标准大小M×M
    img_resized = cv2.resize(image, (M,M))
    #2.对图片使用3×3的高斯低通滤波，减轻JPEG压缩和噪声污染等微小操作队图像哈希的影响。
    kernel=np.array([[1,2,1],
                     [2,4,2],
                     [1,2,1]])/16.0
    #对每一个颜色通道应用滤波
    img_Gaussian=cv2.filter2D(img_resized,-1,kernel)
    #应用7x7的平均滤波器
    # k_size = 14
    # img_blurred = cv2.blur(img_Gaussian, (k_size, k_size))

    #3.将BGR转为RGB
    # RGB_image = cv2.cvtColor(img_Gaussian, cv2.COLOR_BGR2RGB)

    # #4.将BGR转LAB
    # Lab_image=cv2.cvtColor(img_Gaussian,cv2.COLOR_BGR2LAB)
    #
    # #4.将BGR转HSV
    # HSV_image=cv2.cvtColor(img_Gaussian,cv2.COLOR_BGR2HSV)
    #
    # #4.将BGR转HLS
    # HLS_image=cv2.cvtColor(img_Gaussian,cv2.COLOR_BGR2HSV)

    YCbCr_image=cv2.cvtColor(img_Gaussian,cv2.COLOR_BGR2YCrCb)[:,:,0]
    ring_image=divide_into_rings(YCbCr_image,32)
    return ring_image
#
# def NMF(image,k):
#     # can be used for example for dimensionality reduction, source separation or topic extraction
#     # 个人认为最重要的参数是n_components、alpha、l1_ratio、solver
#     # nmf = sklearn.decomposition.NMF(n_components=k, # k value,默认会保留全部特征
#     #           init='nndsvdar',  # W H 的初始化方法，包括'random' | 'nndsvd'(默认) |  'nndsvda' | 'nndsvdar' | 'custom'.
#     #           solver='mu',  # 'cd' | 'mu'
#     #           beta_loss='kullback-leibler', # {'frobenius', 'kullback-leibler', 'itakura-saito'}，一般默认就好
#     #           max_iter=500   # 最大迭代次数
#     #           # random_state=None,
#     #           # alpha=0.,  # 正则化参数
#     #           # l1_ratio=0.,  # 正则化参数
#     #           # verbose=0,  # 冗长模式
#     #           # shuffle=False  # 针对"cd solver"
#     #           )
#     # W = nmf.fit_transform(image)
#     # H = nmf.components_
#     print(0)
#     n_components=k
#     model = sklearn.decomposition.NMF(n_components=n_components,
#                                       init='random',
#                                       solver = "mu",
#                                       beta_loss = "kullback-leibler",
#                                       random_state=0,
#                                       max_iter = 1000)
#     print(0)
#     W = model.fit_transform(image)
#     print(0)
#     H = model.components_
#     H = H.reshape(1, -1)[0]
#     return H

def NMF(image, k):
    try:
        n_components = k
        model = sklearn.decomposition.NMF(n_components=n_components,
                                          init='random',
                                          solver="mu",
                                          beta_loss="kullback-leibler",
                                          random_state=0,
                                          max_iter=1000
                                          )
        W = model.fit_transform(image)
        H = model.components_
        H = H.reshape(1, -1)[0]
        return H
    except Exception as e:
        print(f"An error occurred in the NMF function: {e}")
        traceback.print_exc()  # 添加这一行以打印详细的堆栈跟踪信息
        return None