图片相似度——opencv折腾记

0x01 背景

  业务需要,得对所有图片去重,单纯的md5去重肯定是不满足的,就用到了感知哈希。先前也只是单纯的听说,代码也是同时给的,我就单纯的部署了下。然而,考虑到效率、环境部署等方面的问题,需要将原本计算感知哈希的python脚本换成c的。
  先前也算没有怎么接触过Opencv吧,说有,可能还真有,大学毕业那会儿做人脸识别好像还用到过,但也是许久没碰过了。可能学习方法的问题吧,要解决这事,并没有认真去学习phash是啥、opencv的有关函数等等~怎么说呢,有好有坏吧,折腾过程中的确可以学到更多零碎知识,对比看书学习是有优势。但另一方面,他的劣势(也不一定算是劣势)在于,折腾的时间不可控,学习到的只是不整体。
  我也是完整折腾完了之后,才弄明白什么是pHash,python计算转c++后,有哪些坑等等,这篇文,也就是来总结一下学习到的东西吧。

0x02 图片相似度算法

  计算图片相似度其实有蛮多算法,比如ahash、dhash、phash,由于同事先前给的代码就是python版本的phash。所以,这里也直接按按他的代码来吧。
  网上找到的介绍各算法原理内容如下:

1.ahash算法

  在向下采样的过程中,只保留图像的低频信息。
  工作过程:
  ①缩小尺寸,简化色彩:8*8灰度
  ②计算灰度均值:64个灰度值的平均
  ③生成哈希码:将64个灰度值与上一步生成的平均值进行比较。比均值大,则置1;比均值小,则置0。从而得到一个包含64个元素的“指纹”。
  ④计算汉明距离:将两个“指纹”进行比较,计算它们的汉明距离。(==0):非常相似;(<5): 相似; (>10): 不同

2.phash算法

  均值哈希虽然简单,但是受均值影响大。如果对图像进行伽马校正或者进行直方图均值化都会影响均值,从而影响哈希值的计算。所以就有人提出更健壮的方法,通过离散余弦(DCT)进行低频提取。
  离散余弦变换(DCT)是种图像压缩算法,它将图像从像素域变换到频率域。DCT变换属于变换压缩方法(TransformCompression),变换压缩的一个特点是将从前密度均匀的信息分布变换为密度不同的信息分布。在图像中,低频部分的信息量要大于高频部分的信息量,尽管低频部分的数据量比高频部分的数据量要小的多。例如删除掉占50%存储空间的高频部分,信息量的损失可能还不到5%。
  一般图像都存在很多冗余和相关性的,所以转换到频率域之后,只有很少的一部分频率分量的系数才不为0,大部分系数都为0(或者说接近于0)。下图的右图是对lena图进行离散余弦变换(DCT)得到的系数矩阵图。从左上角依次到右下角,频率越来越高,由图可以看到,左上角的值比较大,到右下角的值就很小很小了。换句话说,图像的能量几乎都集中在左上角这个地方的低频系数上面了。
  工作过程:
  ①缩小尺寸,简化色彩:32*32灰度
  ②计算DCT:对图像进行离散余弦变换,得到32*32的DCT系数矩阵
  ③截取DCT:因为只有左上角的部分呈现图像的最低频部分,所以我们截取左上角的8*8矩阵
  ④计算均值:64个DCT系数的均值
  ⑤生成哈希码:将64个DCT系数与上一步生成的平均值进行比较。之后过程同均值哈希

Alt text
  再看一个实例,下图很好的证实了DCT的压缩性,即,其可逆算回原来的图片,试想,取DCT低频部分,复原后,原图大部分信息还将保存有。

Alt text

3.dhash算法

  这是一种通过渐变实现的哈希算法,相比phash,速度上有较大的优势。
  工作过程:
   ①缩小尺寸、简化色彩:9*8灰度
   ②计算差异:每行9个像素做差得到8个差异值,总共64个差异值
   ③生成哈希码:差异值与0比较。大于0,置1;小于0,置0。之后过程同均值哈希

  python版本的phash主要计算代码如下,与上文引用中的phash介绍基本一致:

def pHash(imgfile):
"""get image pHash value"""
#加载并调整图片为32x32灰度图片
img=cv2.imread(imgfile, cv2.CV_LOAD_IMAGE_GRAYSCALE)
img=cv2.resize(img,(32,32),interpolation=cv2.INTER_CUBIC)
#创建二维列表
h, w = img.shape[:2]
vis0 = np.zeros((h,w), np.float32)
vis0[:h,:w] = img #填充数据
#二维Dct变换
vis1 = cv2.dct(cv2.dct(vis0))
#cv.SaveImage('a.jpg',cv.fromarray(vis0)) #保存图片
vis1.resize(8,8)
#把二维list变成一维list
img_list=flatten(vis1.tolist())
#计算均值
avg = sum(img_list)*1./len(img_list)
avg_list = ['0' if i<avg else '1' for i in img_list]
#得到哈希值
return [int(''.join(avg_list[x:x+4]),2) for x in range(0,64)]

0x02 转C++

  在没有了解phash原理的情况下,我是直接去网上找了一些phash的c代码,然后打算对着改,编译通过后,对着一步步去调试,看两边矩阵的值的情况。
  从网上找到了一段ahash的c++代如下,代码来自http://www.itnose.net/detail/6249373.html
  按照他的代码,以及后来了解到的phash的原理,整理main函数后得到可用的c++代码大致如下:

int main(int argc, char* argv[]){
Mat img = imread(argv[1], 1);
if(!img.data){
cout << "the image is not exist" << endl;
return 0;
}
int size = 32; // 图片缩放后大小
cvtColor(img, img, COLOR_BGR2GRAY); // 灰度化
resize(img, img, Size(size,size), 0,0, CV_INTER_CUBIC); // 缩放到32*32
double **iMatrix = new double*[size];
for(int i = 0; i < size; i++)
iMatrix[i] = new double[size];
DCT(img, size, iMatrix); // 离散余弦变换
DCT1(iMatrix, size, iMatrix); // 第二次余弦变换
char phash[64]={0};
double averagePix = calcAverage(iMatrix, 8); // 计算均值
fingerPrint(iMatrix, 8, phash, averagePix); // 计算hash值
// 为了输出保持一致
cout << "[";
for (int i=0; i<64;i++)
{
char tmp[4]={0};
memcpy(tmp, phash+i, 4);
tmp[4] = '\0';
int re = bin2int(tmp);
if (i != 63) cout << re << ", ";
else cout << re << "]" << endl;
}
return 0;
}

  也就是说,用上面这段c++代码,即可得到与python版本一致的结果,当然,还改了其他的函数。
  

0x04 遇到的坑

  其实这里才是本文的主要内容。
  1、编译
  在centos上搭建好opencv的编译环境后,就开始来折腾了,编译就直接报错:

phash.cpp:(.text+0xe23): undefined reference to `cv::imread(std::string const&, int)'
phash.cpp:(.text+0xe96): undefined reference to `cv::_OutputArray::_OutputArray(cv::Mat&)'
phash.cpp:(.text+0xeaf): undefined reference to `cv::_InputArray::_InputArray(cv::Mat const&)'
phash.cpp:(.text+0xecf): undefined reference to `cv::cvtColor(cv::_InputArray const&, cv::_OutputArray const&, int, int)'
phash.cpp:(.text+0xf05): undefined reference to `cv::_OutputArray::_OutputArray(cv::Mat&)'
phash.cpp:(.text+0xf1e): undefined reference to `cv::_InputArray::_InputArray(cv::Mat const&)'
phash.cpp:(.text+0xf48): undefined reference to `cv::resize(cv::_InputArray const&, cv::_OutputArray const&, cv::Size_<int>, double, double, int)'

  google到编译参数如下:

g++ phash.cpp -lopencv_objdetect -lm -o phash

  于是,能成功编译出可用elf了。
  2、函数参数默认值不同
  各位看官觉得,img初始值相同的情况下,下面两行代码得到的img结果是否会一致?

img=cv2.resize(img,(32,32),interpolation=cv2.INTER_CUBIC)  #python
cv::resize(img, img, cv::Size(32,32), CV_INTER_CUBIC);     //c++

  当然,我也觉得是一致的,可现实打了我的脸:

python脚本处理后,img数组最后几行:
Alt text

c++处理完后,img数组最后几行:
Alt text

  看完我简直惊呆了,怎么都想不通,然后我在测测interpolation参数,你猜当interpolation=CV2.INTER_LINEAR时,python和c++得到的结果是否一直?
  对的,是一致的,然后我就百思不得其解了,开始怀疑是否是精度问题,但想想又不太对劲,python最后还是要调用c库的。然后就去找了找resize的原型。

void resize(InputArray src, OutputArray dst, Size dsize, double fx=0, double fy=0, int interpolation=INTER_LINEAR )

  猛然发现,咦,还有俩参数:fx和fy。是不是他们的影响,于是,改了下上面的代码,固定好fx和fy均为0,如下:

img=cv2.resize(img,(32,32),fx=0,fy=0,interpolation=cv2.INTER_CUBIC)  #python
cv::resize(img, img, cv::Size(32,32),0,0, CV_INTER_CUBIC);     //c++

  是的,这样得到img完全一致。也就是说:

同一幅图片,用c++里的resize和python的resize默认fx和fy的默认参数值,居然是不相同的!!!

  3、dct结果
  是的,上面resize得到的结果终于一致后,同样放到dct下面,得到的结果又不同了:
  python:

Alt text
  c++:

Alt text
  最后一行数组的值总是相差一点,至于为啥,我现在都还没弄明白,这次似乎不是参数的问题,感觉还是精度问题。好在python脚本里做了两次dct,两次dct后,取平均值,再计算phash得到的结果还好是一致的,而计算一次得到的结果却有一定几率不同。
  
  4、根据phash算法描述中,需要取左上角8*8做均值计算。
  然而在python脚本中,np.resize(8,8)后,得到的却是32*2的矩阵,并不是左上角8*8,至于为啥,我也不太清楚,我只有跟着python该c++(与历史问题保持一致)。
  大概,问题就是这么多了吧,因为先前没怎么接触过opencv,在进行mat转换的时候,真的,有点奔溃,不知道其具体的内存布局,这是个很大问题。
  一会儿需要mat,一会儿需要二维数组。
  对了,还有个问题,opencv自带的dct,需要浮点型的Mat,就需要用convertTo来做一次转换,由于这样得到的精度更低,所以,还是用了ahash代码中自己写的DCT算法。
  最后,为了可移植性,想编译出static版本的elf,然后,我放弃了,用ldd找出引用,然后再把lib拷贝过去,在需要运行的机器上写个sh来临时export LD_LIBRARY_PATH,当然,结果是好的,的确可用。

0x05 参考

http://blog.csdn.net/zouxy09/article/details/17471401

代码,其中是有蛮多我的调试记录(注释部分):

#include <iostream>
#include <bitset>
#include <string>
#include <iomanip>
#include <cmath>
#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"

using namespace std;
using namespace cv;

#define PI 3.1415926
#define hashLength 64

/*
功能:获取DCT系数
n:矩阵大小
quotient: 系数
quotientT: 系数转置
*/

void coefficient(const int &n, double **quotient, double **quotientT){
double sqr = 1.0/sqrt(n+0.0);
for(int i = 0; i < n; i++){
quotient[0][i] = sqr;
quotientT[i][0] = sqr;
}

for(int i = 1; i < n; i++){
for(int j = 0; j < n; j++){
quotient[i][j] = sqrt(2.0/n)*cos(i*(j+0.5)*PI/n); // 由公式得到
quotientT[j][i] = quotient[i][j];
}
}

}
/*
功能:两矩阵相乘
A和B:源输入矩阵
result:输出矩阵
*/

void matrixMultiply(double **A, double **B, int n, double **result){
double t = 0;
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
t = 0;
for(int k = 0; k < n; k++)
t += A[i][k]*B[k][j];
result[i][j] = t;
}
}
}


void DCT(Mat_<uchar> image, const int &n, double **iMatrix){
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
iMatrix[i][j] = (double)image(i,j);
}
}

// 为系数分配空间
double **quotient = new double*[n];
double **quotientT = new double*[n];
double **tmp = new double*[n];
for(int i = 0; i < n; i++){
quotient[i] = new double[n];
quotientT[i] = new double[n];
tmp[i] = new double[n];
}
// 计算系数矩阵
coefficient(n, quotient, quotientT);
matrixMultiply(quotient, iMatrix, n, tmp); // 由公式成绩结果
matrixMultiply(tmp, quotientT, n, iMatrix);

for(int i = 0; i < n; i++){
delete []tmp[i];
delete []quotient[i];
delete []quotientT[i];
}
delete []tmp;
delete []quotient;
delete []quotientT;
}

// 计算8*8图像的平均灰度
float calcAverage(double **iMatrix, const int &size){
float sum = 0;
for(int i = 0 ; i < 2; i++){
for(int j = 0; j < 32; j++){
sum += iMatrix[i][j];
}
}
return sum/(size*size);
}

/* 计算hash值
image:8*8的灰度图像
size: 图像大小 8*8
phash:存放64位hash值
averagePix: 灰度值的平均值
*/

void fingerPrint(double **iMatrix, const int &size, char* phash, const float &averagePix){
for(int i = 0; i < 2; i++){
int pos = i * 32;
for(int j = 0; j < 32; j++){
phash[pos+j] = iMatrix[i][j] >= averagePix ? '1':'0';
}
}
//cout << endl;
}

/*计算汉明距离*/
int hammingDistance(const bitset<hashLength> &query, const bitset<hashLength> &target){
int distance = 0;
for(int i = 0; i < hashLength; i++){
distance += (query[i] == target[i] ? 0 : 1);
}
return distance;
}

string bitTohex(const bitset<hashLength> &target){
string str;
for(int i = 0; i < hashLength; i=i+4){
int sum = 0;
string s;
sum += target[i] + (target[i+1]<<1) + (target[i+2]<<2) + (target[i+3]<<3);
stringstream ss;
ss << hex <<sum; // 以十六进制保存
ss >> s;
str += s;
}
return str;
}

void DCT1(double **image, const int &n, double **iMatrix){
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
iMatrix[i][j] = image[i][j];
}
}
double **quotient = new double*[n];
double **quotientT = new double*[n];
double **tmp = new double*[n];
for(int i = 0; i < n; i++){
quotient[i] = new double[n];
quotientT[i] = new double[n];
tmp[i] = new double[n];
}
coefficient(n, quotient, quotientT);
matrixMultiply(quotient, iMatrix, n, tmp); // 由公式成绩结果
matrixMultiply(tmp, quotientT, n, iMatrix);

for(int i = 0; i < n; i++){
delete []tmp[i];
delete []quotient[i];
delete []quotientT[i];
}
delete []tmp;
delete []quotient;
delete []quotientT;
}

int bin2int(char* bin){
int i,sum;
for(sum=i=0;bin[i];(sum*=2)+=bin[i++]-'0');
return sum;
}


int main(int argc, char* argv[]){
Mat img = imread(argv[1], 1);
if(!img.data){
cout << "the image is not exist" << endl;
return 0;
}
int size = 32; // 图片缩放后大小


cvtColor(img, img, COLOR_BGR2GRAY); // 灰度化

resize(img, img, Size(size,size), 0,0, CV_INTER_CUBIC); // 缩放到32*32

double **iMatrix = new double*[size];
for(int i = 0; i < size; i++)
iMatrix[i] = new double[size];

DCT(img, size, iMatrix); // 离散余弦变换
DCT1(iMatrix, size, iMatrix);
/*
Mat tmp;
img.convertTo(tmp, CV_32F);
dct(tmp, tmp);
dct(tmp, tmp);
for (int i=0; i<8; i++)
for (int j=0; j<8; j++)
cout << iMatrix[i][j] << " ";
cout<< "\n" << endl;
*/

char phash[64]={0};
double averagePix = calcAverage(iMatrix, 8);
//cout << averagePix << endl;
//bitset<hashLength> phash;
fingerPrint(iMatrix, 8, phash, averagePix);
//cout << phash << endl;
//string str = bitTohex(phash);
cout << "[";
for (int i=0; i<64;i++)
{
char tmp[4]={0};
memcpy(tmp, phash+i, 4);
tmp[4] = '\0';
//cout << tmp << " ";
int re = bin2int(tmp);
if (i != 63) cout << re << ", ";
else cout << re << "]" << endl;
}
return 0;
}

——Tracy_梓朋
2017年3月29日15:09:01