感知哈希算法

在数字图像处理领域中，认为图像在本质上是一种信号，表达图像的方式有两种，一种是时域方式，它由点阵构成，里面的每一个点存储着该位置对应的颜色值称为像素，整体就形成了人眼中的图像， Windows操作系统的位图格式(Bitmap)就是一种典型的时域信号形式；另一种是频域方式，图像由一组不同频率的正弦波和余弦波叠加而成。时域信号与频域信号可以相互转换，常见的有傅立叶变换FFT、离散余弦变换DCT及他们的逆变换iFFT、iDCT等，目前广为使用的JEPG格式使用的就是DCT/iDCT变换。

频域信号的特点是低频信号表达了图像的轮廓信息，高频信号则表达了细节信息。pHash算法正是利用这一特点，通过对比图像间低频信号特征来计算图像的相似度。以下为不同频率的信号对图像的影响情况：

DCT

$$ \begin{align} F(0, 0) &= \frac{1}{n}\sum^{n-1}\sum^{n-1}f(x,y) \\ F(0, v) &= \frac{\sqrt{2}}{n}\sum^{n-1}\sum^{n-1}f(x,y) \cdot \cos\frac{(2y+1)v\pi}{2n} \\ F(u, 0) &= \frac{\sqrt{2}{n}}\sum^{n-1}\sum^{n-1}f(x,y) \cdot \cos\frac{(2x+1)u\pi}{2n} \\ F(u, v) &= \frac{\sqrt{2}{n}}\sum^{n-1}\sum^{n-1}f(x,y) \cdot \cos\frac{(2x+1)u\pi}{2n}\cos\frac{(2y+1)v\pi}{2n} \end{align} $$

其中，$u,v$是频域系数矩阵元素的座标，$F(u,v)$是元素的值； $x, y$ 是时域像素矩阵元素的座标，$f(x,y)$是元素的值；$n$是两个矩阵的阶。

算法描述

  • 步骤1图像处理：首先将图片转换成灰阶并缩小至固定的$M_{32×32}$
  • 步骤2 信号变换：使用DCT将图片$M_{32×32}$转换为频域系数矩阵，记作$C_{32×32}$
  • 步骤3 低频截取：从$C_{32×32}$左上角第二行第二列开始截取获得低频矩阵，记作$N_{8×8}$
  • 步骤4 二值化：计算矩阵$N_{8×8}$的均值，记为$m$，遍历矩阵$N_{8×8}$的每个元素e，若$e \gt m$则记为1，否则记0。得到表征图像特征的二进制串哈希值。
  • 步骤5 相似度计算：计算两个图像哈希值的海明距离来表征图像间的相似度，距离越大则相似度低，反之则相似度越高

实践中也有在步骤1图像处理的时候使用模糊滤镜进行去噪，以及在步骤4二值化中使用中值代替均值的做法，需要进行性能测试验证其有效性。

算法实现

https://github.com/klesh/qt-phash

性能测试

测试目标

测量在感知哈希算法中，对图像处理时使用不同的去噪算法及二值时使用不同的阈值对最终表现的影响。根据测试结果选择最优组合。

测试设计

测试数据：使用一组随机图片作为样本，并对应生成五组对照图片，分别为缩略组、旋转组、模糊组、裁切组、高对比组。即，同一张图片分别对应有五张相似图片，不同的图片间则互为不相似图片。

测试组合：图像处理和二值化是算法的不同步骤，因此需要分别取他们的可能值形成所有可能组合，并对这些组合进行对比。图像处理时的去噪算法有三种可能取值：不去噪、盒模糊去噪及高斯模糊；二值化时阈值有两种取法：均值、中值。共有六种可能组合：不去噪-均值、盒模糊-均值、高斯模糊-均值、不去噪-中值、盒模糊-中值、高斯模糊-中值。

测试方法：图像哈希算法应反映出图片间真实的差异程度，因此，它应在同时满足两个条件：一、 对相似的两张图片相似输出较低的海明距离，二、 对不相似的两张图片输出较高的海明距离。因此，首先要计算出这些组合对于相似图片及不相似图片输出的距离数值，对比这些距离的均值将反映出组合对性能的影响，距离的方差反映出对算法稳定性的影响。

测试实现

图像卷积操作实现：

 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

QImage convolve(const QImage &source, const float (&kernel)[KS][KS]) {
    assert(source.format() == QImage::Format_Grayscale8);
    const int w = source.width(), h = source.height();
    QImage target(w, h, source.format());
    const int kcx = KS/2, kcy = kcx; // center position of the kernel
    #pragma omp parallel for
    for (int y = 0; y < h; y++) { // fill out target pixel one by one.
      for (int x = 0; x < w; x++) {
        int ksx = x - kcx, ksy = y - kcy, kex = x + kcx + 1, key = y + kcy + 1;
        assert ( kex - ksx == KS );
        assert ( key - ksy == KS );
        ksx = ksx < 0 ? 0 : ksx; // kernel crop. (called Neumann method in CImg)
        ksy = ksy < 0 ? 0 : ksy; // ksx = kernel starting x, kex = kernel ending x
        kex = kex > w ? w : kex;
        key = key > h ? h : key;
        // calculate pixel value by image kernel
        float totalValue = 0.0, totalWeight = 0.0;
        for (int ky = ksy; ky < key; ky++) {
          for (int kx = ksx; kx < kex; kx++) {
            int i = kx - ksx, j = ky - ksy;
            float weight = kernel[i][j];
            totalWeight += weight;
            float value = source.scanLine(ky)[kx];
            totalValue += weight * value;
          }
        }
        target.scanLine(y)[x] = totalValue / totalWeight;
      }
   }
    return target;
}

求第k个最小值算法实现：

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

template<typename T>
T kthSmallest(T *array, int index, int low, int high) {
    int l = low, r = high;
    while (true) {
      int m = (l + r + 1) >> 1;
      if (array[l] > array[m]) std::swap(array[l], array[m]);
      if (array[r] > array[m]) std::swap(array[r], array[m]);
      if (array[l] > array[r]) std::swap(array[l], array[r]);
      int  p = r--;
      while (true) {
        while (array[l] < array[p]) l++;
        while (array[r] > array[p]) r--;
        if (l > r) break;
        std::swap(array[l], array[r]);
      }
      std::swap(array[l], array[p]);
      if (l == index) break;
      else if (index < l) r = high = l - 1, l = low; 
      else if (index > l) l = low = l + 1, r = high;
    }
    return array[index];
};

求取矩阵中值算法实现：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

template<int N, int M, typename T>
T matrixMedian(QGenericMatrix<N, M, T> matrix) {
    const int S = N * M, m = S >> 1, n = m - 1, e = S - 1;
    T array[S];
    matrix.copyDataTo(array);
    const T mv = kthSmallest(array, m, 0, e);
    if (S % 2) return mv;
    const T nv = kthSmallest(array, n, 0, n);
    return (mv + nv) / 2;
}

测试数据获取与生成：unsplash.com 是一个提供免费图片的网站，同时提供了随机图片获取的接口。利用 UNIX 系统下的命令和工具可以很方便地获得所需的测试数据。之后再通过开源 imagemagick 工具软件，可以很方便地对图片进行变换生成我们需要的相似图片：

获取随机图片脚本：

1
2
3
4
5

#!/usr/bin/env fish

for i in (seq 99)
  curl -L 'https://source.unsplash.com/random' -o original/(printf %02d $i).jpg
end

生成相似图片脚本：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

#!/usr/bin/env fish

rm -rf samples testsets
mkdir -p samples \
  testsets/thumbnail \
  testsets/rotate \
  testsets/blur \
  testsets/crop \
  testsets/contrast
for img in (ls original)
  convert original/$img -resize 640x480 samples/$img
  convert samples/$img -blur 2x2 testsets/blur/$img
  convert samples/$img -rotate 1 testsets/rotate/$img
  convert samples/$img -resize 120x80 testsets/thumbnail/$img
  convert samples/$img -gravity Center -crop 98%x+1%+1% testsets/crop/$img
  convert samples/$img -level 10% testsets/contrast/$img
end

测试结果

本次测试使用了100张随机下载的图片，分别抽取20张/50张和全量数量进行测量。结果如下：

20 张图片测试结果

组合 相似图片 不相似图片
Mean Variance Mean Variance
none-mean 2.83 11.18 20.92 106.18
none-median 2.64 12.39 20.93 106.15
box-mean 2.48 9.67 21.26 114.47
box-median 2.82 7.85 20.97 105.99
gaussian-mean 2.34 8.06 20.99 108.08
gaussian-median 2.8 10.24 21.05 108.16

50 张图片测试结果

组合 相似图片 不相似图片
Mean Variance Mean Variance
none-mean 2.8 10.72 20.9 85.56
none-median 2.85 11.76 21.12 86.13
box-mean 3.05 13.3 21.06 84.17
box-median 3.26 9.96 21.09 86.39
gaussian-mean 2.44 8.06 21.07 84.48
gaussian-median 2.85 10.1 21.07 85.44

100 张图片测试结果

组合 相似图片 不相似图片
Mean Variance Mean Variance
none-mean 2.77 11.33 20.6 110.74
none-median 2.85 11.05 21.37 114.27
box-mean 2.83 12.28 21.73 124.97
box-median 3.28 10.03 21.38 117.19
gaussian-mean 2.43 8.45 21.62 120.99
gaussian-median 2.71 9.28 21.31 116.77

从测试结果可以看出，在不同数量的数据集上高斯模糊和均值的组合性能表现相对较好，同时也较为稳定。