图像直方图比较

对于直方图的比较,我们可以使用 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种:

  1. HISTCMP_CORREL:相关性
  2. HISTCMP_CHISQR:卡方
  3. HISTCMP_INTERSECT:交集
  4. HISTCMP_BHATTACHARYYA:巴氏距离

而根据 HistCompMethods[2],OpenCV的 compareHist 函数计算直方图匹配度所提供的指标有以下7种(实际为6种,巴氏和海林格算一种):

  1. HISTCMP_CORREL:相关性(Correlation)

    d ( H 1 , H 2 ) = ∑ I ( H 1 ( I ) − H 1 ˉ ) ( H 2 ( I ) − H 2 ˉ ) ∑ I ( H 1 ( I ) − H 1 ˉ ) 2 ∑ I ( H 2 ( I ) − H 2 ˉ ) 2 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}} d(H1,H2)=I(H1(I)H1ˉ)2I(H2(I)H2ˉ)2 I(H1(I)H1ˉ)(H2(I)H2ˉ)

    其中,

    H 1 ˉ = 1 N ∑ I H 1 ( I ) \bar{H_1} = \frac{1}{N}\sum_I H_1(I) H1ˉ=N1IH1(I)

    N N N是直方图桶的数量。

  2. 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

  3. 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))

  4. 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)=1H1ˉH2ˉN2 1IH1(I)H2(I)

  5. HISTCMP_HELLINGER:海林格距离(Hellinger distance),同巴氏距离。

  6. 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

    常被用于材质比较。

  7. 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

参考资料

  1. Histogram Comparison
  2. HistCompMethods
  • 12
    点赞
  • 11
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Humbunklung

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值