对于直方图的比较,我们可以使用 OpenCV 提供的函数 compareHist()
进行比较,从而得到一个数值,表示两个直方图的匹配程度(相似性)。
原理
对于两个直方图( H 1 H_{1} H1 和 H 2 H_{2} H2),我们首先选用一个指标( d ( H 1 , H 2 ) d(H_{1}, H_{2}) d(H1,H2)),表示两个直方图的匹配度。
根据 Histogram Comparison 参考资料[1],OpenCV的 compareHist
函数计算直方图匹配度所提供的指标有以下4种:
HISTCMP_CORREL
:相关性HISTCMP_CHISQR
:卡方HISTCMP_INTERSECT
:交集HISTCMP_BHATTACHARYYA
:巴氏距离
而根据 HistCompMethods
[2],OpenCV的 compareHist
函数计算直方图匹配度所提供的指标有以下7种(实际为6中,巴氏和海林格算一种):
-
HISTCMP_CORREL
:相关性(Correlation)$
d(H_1,H_2) = \frac{\sum_I (H_1(I) - \bar{H_1}) (H_2(I) - \bar{H_2})}{\sqrt{\sum_I(H_1(I) - \bar{H_1})^2 \sum_I(H_2(I) - \bar{H_2})^2}}
$其中,
$
\bar{H_1} = \frac{1}{N}\sum_I H_1(I)
$且 N N N是直方图桶的数量。
-
HISTCMP_CHISQR
:卡方(Chi-Square)d ( H 1 , H 2 ) = ∑ I ( H 1 ( I ) − H 2 ( I ) ) 2 H 1 ( I ) d(H_1,H_2) = \sum _I \frac{\left(H_1(I)-H_2(I)\right)^2}{H_1(I)} d(H1,H2)=∑IH1(I)(H1(I)−H2(I))2
-
HISTCMP_INTERSECT
:交集(Intersection)d ( H 1 , H 2 ) = ∑ I min ( H 1 ( I ) , H 2 ( I ) ) d(H_1,H_2) = \sum _I \min (H_1(I), H_2(I)) d(H1,H2)=∑Imin(H1(I),H2(I))
-
HISTCMP_BHATTACHARYYA
:巴氏距离(Bhattacharyya distance),事实上 OpenCV 计算的是与巴氏距离系数相关的海林格距离(Hellinger distance)。d ( H 1 , H 2 ) = 1 − 1 H 1 ˉ H 2 ˉ N 2 ∑ I H 1 ( I ) ⋅ H 2 ( I ) d(H_1,H_2) = \sqrt{1 - \frac{1}{\sqrt{\bar{H_1} \bar{H_2} N^2}} \sum_I \sqrt{H_1(I) \cdot H_2(I)}} d(H1,H2)=1−H1ˉH2ˉN2 1∑IH1(I)⋅H2(I)
-
HISTCMP_HELLINGER
:海林格距离(Hellinger distance),同巴氏距离。 -
HISTCMP_CHISQR_ALT
:修正卡方(Alternative Chi-Square)d ( H 1 , H 2 ) = 2 × ∑ I ( H 1 ( I ) − H 2 ( I ) ) 2 H 1 ( I ) + H 2 ( I ) d(H_1,H_2) = 2 × \sum _I \frac{\left(H_1(I)-H_2(I)\right)^2}{H_1(I)+H_2(I)} d(H1,H2)=2×∑IH1(I)+H2(I)(H1(I)−H2(I))2
常被用于材质比较。
-
HISTCMP_KL_DIV
:库尔巴克-莱布勒散度(Kullback-Leibler divergence)d ( H 1 , H 2 ) = ∑ I H 1 ( I ) log ( H 1 ( I ) H 2 ( I ) ) d(H_1,H_2) = \sum _I H_1(I) \log \left(\frac{H_1(I)}{H_2(I)}\right) d(H1,H2)=∑IH1(I)log(H2(I)H1(I))
OpenCV 中的声明与定义(以 4.10 版为例)
-
OpenCV 中
compareHist
函数声明:/** @brief Compares two histograms. The function cv::compareHist compares two dense or two sparse histograms using the specified method. The function returns \f$d(H_1, H_2)\f$ . While the function works well with 1-, 2-, 3-dimensional dense histograms, it may not be suitable for high-dimensional sparse histograms. In such histograms, because of aliasing and sampling problems, the coordinates of non-zero histogram bins can slightly shift. To compare such histograms or more general sparse configurations of weighted points, consider using the #EMD function. @param H1 First compared histogram. @param H2 Second compared histogram of the same size as H1 . @param method Comparison method, see #HistCompMethods */ CV_EXPORTS_W double compareHist( InputArray H1, InputArray H2, int method ); /** @overload */ CV_EXPORTS double compareHist( const SparseMat& H1, const SparseMat& H2, int method );
-
OpenCV 中
HistCompMethods
枚举:enum HistCompMethods { /** Correlation \f[d(H_1,H_2) = \frac{\sum_I (H_1(I) - \bar{H_1}) (H_2(I) - \bar{H_2})}{\sqrt{\sum_I(H_1(I) - \bar{H_1})^2 \sum_I(H_2(I) - \bar{H_2})^2}}\f] where \f[\bar{H_k} = \frac{1}{N} \sum _J H_k(J)\f] and \f$N\f$ is a total number of histogram bins. */ HISTCMP_CORREL = 0, /** Chi-Square \f[d(H_1,H_2) = \sum _I \frac{\left(H_1(I)-H_2(I)\right)^2}{H_1(I)}\f] */ HISTCMP_CHISQR = 1, /** Intersection \f[d(H_1,H_2) = \sum _I \min (H_1(I), H_2(I))\f] */ HISTCMP_INTERSECT = 2, /** Bhattacharyya distance (In fact, OpenCV computes Hellinger distance, which is related to Bhattacharyya coefficient.) \f[d(H_1,H_2) = \sqrt{1 - \frac{1}{\sqrt{\bar{H_1} \bar{H_2} N^2}} \sum_I \sqrt{H_1(I) \cdot H_2(I)}}\f] */ HISTCMP_BHATTACHARYYA = 3, HISTCMP_HELLINGER = HISTCMP_BHATTACHARYYA, //!< Synonym for HISTCMP_BHATTACHARYYA /** Alternative Chi-Square \f[d(H_1,H_2) = 2 * \sum _I \frac{\left(H_1(I)-H_2(I)\right)^2}{H_1(I)+H_2(I)}\f] This alternative formula is regularly used for texture comparison. See e.g. @cite Puzicha1997 */ HISTCMP_CHISQR_ALT = 4, /** Kullback-Leibler divergence \f[d(H_1,H_2) = \sum _I H_1(I) \log \left(\frac{H_1(I)}{H_2(I)}\right)\f] */ HISTCMP_KL_DIV = 5 };
-
OpenCV 中
compareHist
函数定义:
// C O M P A R E H I S T O G R A M S
double cv::compareHist( InputArray _H1, InputArray _H2, int method )
{
CV_INSTRUMENT_REGION();
Mat H1 = _H1.getMat(), H2 = _H2.getMat();
const Mat* arrays[] = {&H1, &H2, 0};
Mat planes[2];
NAryMatIterator it(arrays, planes);
double result = 0;
int j;
CV_Assert( H1.type() == H2.type() && H1.depth() == CV_32F );
double s1 = 0, s2 = 0, s11 = 0, s12 = 0, s22 = 0;
CV_Assert( it.planes[0].isContinuous() && it.planes[1].isContinuous() );
for( size_t i = 0; i < it.nplanes; i++, ++it )
{
const float* h1 = it.planes[0].ptr<float>();
const float* h2 = it.planes[1].ptr<float>();
const int len = it.planes[0].rows*it.planes[0].cols*H1.channels();
j = 0;
if( (method == CV_COMP_CHISQR) || (method == CV_COMP_CHISQR_ALT))
{
for( ; j < len; j++ )
{
double a = h1[j] - h2[j];
double b = (method == CV_COMP_CHISQR) ? h1[j] : h1[j] + h2[j];
if( fabs(b) > DBL_EPSILON )
result += a*a/b;
}
}
else if( method == CV_COMP_CORREL )
{
#if (CV_SIMD_64F || CV_SIMD_SCALABLE_64F)
v_float64 v_s1 = vx_setzero_f64();
v_float64 v_s2 = vx_setzero_f64();
v_float64 v_s11 = vx_setzero_f64();
v_float64 v_s12 = vx_setzero_f64();
v_float64 v_s22 = vx_setzero_f64();
for ( ; j <= len - VTraits<v_float32>::vlanes(); j += VTraits<v_float32>::vlanes())
{
v_float32 v_a = vx_load(h1 + j);
v_float32 v_b = vx_load(h2 + j);
// 0-1
v_float64 v_ad = v_cvt_f64(v_a);
v_float64 v_bd = v_cvt_f64(v_b);
v_s12 = v_muladd(v_ad, v_bd, v_s12);
v_s11 = v_muladd(v_ad, v_ad, v_s11);
v_s22 = v_muladd(v_bd, v_bd, v_s22);
v_s1 = v_add(v_s1, v_ad);
v_s2 = v_add(v_s2, v_bd);
// 2-3
v_ad = v_cvt_f64_high(v_a);
v_bd = v_cvt_f64_high(v_b);
v_s12 = v_muladd(v_ad, v_bd, v_s12);
v_s11 = v_muladd(v_ad, v_ad, v_s11);
v_s22 = v_muladd(v_bd, v_bd, v_s22);
v_s1 = v_add(v_s1, v_ad);
v_s2 = v_add(v_s2, v_bd);
}
s12 += v_reduce_sum(v_s12);
s11 += v_reduce_sum(v_s11);
s22 += v_reduce_sum(v_s22);
s1 += v_reduce_sum(v_s1);
s2 += v_reduce_sum(v_s2);
#elif CV_SIMD && 0 //Disable vectorization for CV_COMP_CORREL if f64 is unsupported due to low precision
v_float32 v_s1 = vx_setzero_f32();
v_float32 v_s2 = vx_setzero_f32();
v_float32 v_s11 = vx_setzero_f32();
v_float32 v_s12 = vx_setzero_f32();
v_float32 v_s22 = vx_setzero_f32();
for (; j <= len - VTraits<v_float32>::vlanes(); j += VTraits<v_float32>::vlanes())
{
v_float32 v_a = vx_load(h1 + j);
v_float32 v_b = vx_load(h2 + j);
v_s12 = v_muladd(v_a, v_b, v_s12);
v_s11 = v_muladd(v_a, v_a, v_s11);
v_s22 = v_muladd(v_b, v_b, v_s22);
v_s1 += v_a;
v_s2 += v_b;
}
s12 += v_reduce_sum(v_s12);
s11 += v_reduce_sum(v_s11);
s22 += v_reduce_sum(v_s22);
s1 += v_reduce_sum(v_s1);
s2 += v_reduce_sum(v_s2);
#endif
for( ; j < len; j++ )
{
double a = h1[j];
double b = h2[j];
s12 += a*b;
s1 += a;
s11 += a*a;
s2 += b;
s22 += b*b;
}
}
else if( method == CV_COMP_INTERSECT )
{
#if (CV_SIMD_64F || CV_SIMD_SCALABLE_64F)
v_float64 v_result = vx_setzero_f64();
for ( ; j <= len - VTraits<v_float32>::vlanes(); j += VTraits<v_float32>::vlanes())
{
v_float32 v_src = v_min(vx_load(h1 + j), vx_load(h2 + j));
v_result = v_add(v_result, v_add(v_cvt_f64(v_src), v_cvt_f64_high(v_src)));
}
result += v_reduce_sum(v_result);
#elif CV_SIMD
v_float32 v_result = vx_setzero_f32();
for (; j <= len - VTraits<v_float32>::vlanes(); j += VTraits<v_float32>::vlanes())
{
v_float32 v_src = v_min(vx_load(h1 + j), vx_load(h2 + j));
v_result = v_add(v_result, v_src);
}
result += v_reduce_sum(v_result);
#endif
for( ; j < len; j++ )
result += std::min(h1[j], h2[j]);
}
else if( method == CV_COMP_BHATTACHARYYA )
{
#if (CV_SIMD_64F || CV_SIMD_SCALABLE_64F)
v_float64 v_s1 = vx_setzero_f64();
v_float64 v_s2 = vx_setzero_f64();
v_float64 v_result = vx_setzero_f64();
for ( ; j <= len - VTraits<v_float32>::vlanes(); j += VTraits<v_float32>::vlanes())
{
v_float32 v_a = vx_load(h1 + j);
v_float32 v_b = vx_load(h2 + j);
v_float64 v_ad = v_cvt_f64(v_a);
v_float64 v_bd = v_cvt_f64(v_b);
v_s1 = v_add(v_s1, v_ad);
v_s2 = v_add(v_s2, v_bd);
v_result = v_add(v_result, v_sqrt(v_mul(v_ad, v_bd)));
v_ad = v_cvt_f64_high(v_a);
v_bd = v_cvt_f64_high(v_b);
v_s1 = v_add(v_s1, v_ad);
v_s2 = v_add(v_s2, v_bd);
v_result = v_add(v_result, v_sqrt(v_mul(v_ad, v_bd)));
}
s1 += v_reduce_sum(v_s1);
s2 += v_reduce_sum(v_s2);
result += v_reduce_sum(v_result);
#elif CV_SIMD && 0 //Disable vectorization for CV_COMP_BHATTACHARYYA if f64 is unsupported due to low precision
v_float32 v_s1 = vx_setzero_f32();
v_float32 v_s2 = vx_setzero_f32();
v_float32 v_result = vx_setzero_f32();
for (; j <= len - VTraits<v_float32>::vlanes(); j += VTraits<v_float32>::vlanes())
{
v_float32 v_a = vx_load(h1 + j);
v_float32 v_b = vx_load(h2 + j);
v_s1 += v_a;
v_s2 += v_b;
v_result += v_sqrt(v_a * v_b);
}
s1 += v_reduce_sum(v_s1);
s2 += v_reduce_sum(v_s2);
result += v_reduce_sum(v_result);
#endif
for( ; j < len; j++ )
{
double a = h1[j];
double b = h2[j];
result += std::sqrt(a*b);
s1 += a;
s2 += b;
}
}
else if( method == CV_COMP_KL_DIV )
{
for( ; j < len; j++ )
{
double p = h1[j];
double q = h2[j];
if( fabs(p) <= DBL_EPSILON ) {
continue;
}
if( fabs(q) <= DBL_EPSILON ) {
q = 1e-10;
}
result += p * std::log( p / q );
}
}
else
CV_Error( cv::Error::StsBadArg, "Unknown comparison method" );
}
if( method == CV_COMP_CHISQR_ALT )
result *= 2;
else if( method == CV_COMP_CORREL )
{
size_t total = H1.total();
double scale = 1./total;
double num = s12 - s1*s2*scale;
double denom2 = (s11 - s1*s1*scale)*(s22 - s2*s2*scale);
result = std::abs(denom2) > DBL_EPSILON ? num/std::sqrt(denom2) : 1.;
}
else if( method == CV_COMP_BHATTACHARYYA )
{
s1 *= s2;
s1 = fabs(s1) > FLT_EPSILON ? 1./std::sqrt(s1) : 1.;
result = std::sqrt(std::max(1. - result*s1, 0.));
}
return result;
}
应用示例
比较两个图像的灰度图匹配度。
#include "opencv2/imgcodecs.hpp"
#include "opencv2/highgui.hpp"
#include "opencv2/imgproc.hpp"
#include <iostream>
using namespace cv;
using std::cout;
using std::endl;
using std::string;
// 命令行参数
const char* keys =
"{ help h| | Print help message. }"
"{ @input1 |wukong1.jpg | Path to input image 1. }"
"{ @input2 |wukong2.jpg | Path to input image 2. }";
const string CMP_METHODS[]{"Correlation", "Chi-Squared", "Intersection", "Bhattacharyya", "Hellinger", "Alternative Chi-Squared", "Kullback-Leibler Divergence"};
int main(int argc, char** argv)
{
// 命令行解析
CommandLineParser parser(argc, argv, keys);
// 加载图像
Mat img1 = imread(parser.get<String>("@input1"));
Mat img2 = imread(parser.get<String>("@input2"));
// 检查图像是否加载成功
if (img1.empty() || img2.empty())
{
cout << "Could not open or find the image!\n" << endl;
cout << "Usage: " << argv[0] << " <Input_Image1> <Input_Image2>" << endl;
parser.printMessage();
return EXIT_FAILURE;
}
// 转换为灰度图
Mat gray1, gray2;
cvtColor(img1, gray1, COLOR_BGR2GRAY);
cvtColor(img2, gray2, COLOR_BGR2GRAY);
// 计算直方图
int histSize = 256; // 256个灰度级,桶数
// 灰度级范围
float range[] = { 0, 256 };
const float* histRange = { range };
// 其它直方图计算参数
bool uniform = true, accumulate = false;
Mat hist1, hist2;
calcHist(&gray1, 1, 0, Mat(), hist1, 1, &histSize, &histRange, uniform, accumulate);
calcHist(&gray2, 1, 0, Mat(), hist2, 1, &histSize, &histRange, uniform, accumulate);
// 归一化直方图
normalize(hist1, hist1, 0, 1, NORM_MINMAX, -1, Mat());
normalize(hist2, hist2, 0, 1, NORM_MINMAX, -1, Mat());
// 比较直方图
double correlation = compareHist(hist1, hist2, HISTCMP_CORREL);
cout << "Similarity between the two images with " << CMP_METHODS[0] << ": " << correlation << endl;
double chiSquared = compareHist(hist1, hist2, HISTCMP_CHISQR);
cout << "Similarity between the two images with " << CMP_METHODS[1] << ": " << chiSquared << endl;
double intersection = compareHist(hist1, hist2, HISTCMP_INTERSECT);
cout << "Similarity between the two images with " << CMP_METHODS[2] << ": " << intersection << endl;
double bhattacharyya = compareHist(hist1, hist2, HISTCMP_BHATTACHARYYA);
cout << "Similarity between the two images with " << CMP_METHODS[3] << ": " << bhattacharyya << endl;
double hellinger = compareHist(hist1, hist2, HISTCMP_HELLINGER);
cout << "Similarity between the two images with " << CMP_METHODS[4] << ": " << hellinger << endl;
double alternativeChiSquared = compareHist(hist1, hist2, HISTCMP_CHISQR_ALT);
cout << "Similarity between the two images with " << CMP_METHODS[5] << ": " << alternativeChiSquared << endl;
double kullbackLeibler = compareHist(hist1, hist2, HISTCMP_KL_DIV);
cout << "Similarity between the two images with " << CMP_METHODS[6] << ": " << kullbackLeibler << endl;
// 显示图像
imshow("Image 1", img1);
imshow("Image 2", img2);
waitKey(0);
return EXIT_SUCCESS;
}
import cv2
import numpy as np
# 读取图像
image1_path = '../data/Histogram_Comparison_Source_0.jpg'
image2_path = '../data/Histogram_Comparison_Source_1.jpg'
image1 = cv2.imread(image1_path)
image2 = cv2.imread(image2_path)
# 转换为灰度图像
gray_image1 = cv2.cvtColor(image1, cv2.COLOR_BGR2GRAY)
gray_image2 = cv2.cvtColor(image2, cv2.COLOR_BGR2GRAY)
# 计算灰度图像的直方图
hist1 = cv2.calcHist([gray_image1], [0], None, [256], [0, 256])
hist2 = cv2.calcHist([gray_image2], [0], None, [256], [0, 256])
# 归一化直方图
cv2.normalize(hist1, hist1, alpha=0, beta=1, norm_type=cv2.NORM_MINMAX)
cv2.normalize(hist2, hist2, alpha=0, beta=1, norm_type=cv2.NORM_MINMAX)
# 使用compareHist函数比较直方图
comparison_method = cv2.HISTCMP_CORREL # 相关性比较方法
corr_similarity = cv2.compareHist(hist1, hist2, comparison_method)
chi_squared_similarity = cv2.compareHist(hist1, hist2, cv2.HISTCMP_CHISQR)
intersection_similarity = cv2.compareHist(hist1, hist2, cv2.HISTCMP_INTERSECT)
bhattacharyya_similarity = cv2.compareHist(hist1, hist2, cv2.HISTCMP_BHATTACHARYYA)
hellinger_similarity = cv2.compareHist(hist1, hist2, cv2.HISTCMP_HELLINGER)
alternative_chi_squared_similarity = cv2.compareHist(hist1, hist2, cv2.HISTCMP_CHISQR_ALT)
kbl_divergence_similarity = cv2.compareHist(hist1, hist2, cv2.HISTCMP_KL_DIV)
print(f"Similarity between the two images with Correlation: {corr_similarity}")
print(f"Similarity between the two images with Chi-Squared: {chi_squared_similarity}")
print(f"Similarity between the two images with Intersection: {intersection_similarity}")
print(f"Similarity between the two images with Bhattacharyya: {bhattacharyya_similarity}")
print(f"Similarity between the two images with Hellinger: {hellinger_similarity}")
print(f"Similarity between the two images with Alternative Chi-Squared: {alternative_chi_squared_similarity}")
print(f"Similarity between the two images with Kullback-Leibler Divergence: {kbl_divergence_similarity}")
# 显示结果
cv2.imshow('Image 1', image1)
cv2.imshow('Image 2', image2)
cv2.waitKey(0)
cv2.destroyAllWindows()
Similarity between the two images with Correlation: 0.5185838942269291
Similarity between the two images with Chi-Squared: 53.894400613305585
Similarity between the two images with Intersection: 25.52862914837897
Similarity between the two images with Bhattacharyya: 0.358873990739656
Similarity between the two images with Hellinger: 0.358873990739656
Similarity between the two images with Alternative Chi-Squared: 34.90277478480459
Similarity between the two images with Kullback-Leibler Divergence: 71.12863163014705