【第二部分图像处理】第3章 Opencv图像处理进阶【3 直方图与匹配 C】

3.4直方图对比

3.4.1直方图对比概述

要比较两个直方图( and ), 首先必须要选择一个衡量直方图相似度的对比标准。OpenCV 函数 compareHist 执行了具体的直方图对比的任务。该函数提供了4种对比标准来计算相似度：
 相关：Correlation ( CV_COMP_CORREL )
这里写图片描述
其中

是直方图中bin的数目。
 卡方：Chi-Square ( CV_COMP_CHISQR )

 直方图相交：Intersection ( CV_COMP_INTERSECT )

 Bhattacharyya 距离( CV_COMP_BHATTACHARYYA )

3.4.2直方图对比相关API及源码

 compareHist()函数原型

C++: double compareHist(InputArray H1, 
                        InputArray H2, 
                        int method)

C++: double compareHist(const SparseMat& H1, 
                        const SparseMat& H2, 
                        int method)

【参数】
第一个参数，H1 – First compared histogram.
第二个参数，H2 – Second compared histogram of the same size as H1 .
第三个参数，method -Comparison method that could be one of the following:
 CV_COMP_CORREL Correlation
 CV_COMP_CHISQR Chi-Square
 CV_COMP_INTERSECT Intersection
 CV_COMP_BHATTACHARYYA Bhattacharyya distance

 compareHist()函数源代码

/*【compareHist ( )源代码】***********************************************************
 * @Version：OpenCV 3.0.0（Opnencv2和Opnencv3差别不大，Linux和PC的对应版本源码完全一样，均在对应的安装目录下）  
 * @源码路径：…\opencv\sources\modules\imgproc\src\histogram.cpp
 * @起始行数：2272行   
********************************************************************************/
double cv::compareHist( InputArray _H1, InputArray _H2, int method )
{
    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, len = (int)it.size;

    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() );

#if CV_SSE2
    bool haveSIMD = checkHardwareSupport(CV_CPU_SSE2);
#endif

    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>();
        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_SSE2
            if (haveSIMD)
            {
                __m128d v_s1 = _mm_setzero_pd(), v_s2 = v_s1;
                __m128d v_s11 = v_s1, v_s22 = v_s1, v_s12 = v_s1;

                for ( ; j <= len - 4; j += 4)
                {
                    __m128 v_a = _mm_loadu_ps(h1 + j);
                    __m128 v_b = _mm_loadu_ps(h2 + j);

                    // 0-1
                    __m128d v_ad = _mm_cvtps_pd(v_a);
                    __m128d v_bd = _mm_cvtps_pd(v_b);
                    v_s12 = _mm_add_pd(v_s12, _mm_mul_pd(v_ad, v_bd));
                    v_s11 = _mm_add_pd(v_s11, _mm_mul_pd(v_ad, v_ad));
                    v_s22 = _mm_add_pd(v_s22, _mm_mul_pd(v_bd, v_bd));
                    v_s1 = _mm_add_pd(v_s1, v_ad);
                    v_s2 = _mm_add_pd(v_s2, v_bd);

                    // 2-3
                    v_ad = _mm_cvtps_pd(_mm_castsi128_ps(_mm_srli_si128(_mm_castps_si128(v_a), 8)));
                    v_bd = _mm_cvtps_pd(_mm_castsi128_ps(_mm_srli_si128(_mm_castps_si128(v_b), 8)));
                    v_s12 = _mm_add_pd(v_s12, _mm_mul_pd(v_ad, v_bd));
                    v_s11 = _mm_add_pd(v_s11, _mm_mul_pd(v_ad, v_ad));
                    v_s22 = _mm_add_pd(v_s22, _mm_mul_pd(v_bd, v_bd));
                    v_s1 = _mm_add_pd(v_s1, v_ad);
                    v_s2 = _mm_add_pd(v_s2, v_bd);
                }

                double CV_DECL_ALIGNED(16) ar[10];
                _mm_store_pd(ar, v_s12);
                _mm_store_pd(ar + 2, v_s11);
                _mm_store_pd(ar + 4, v_s22);
                _mm_store_pd(ar + 6, v_s1);
                _mm_store_pd(ar + 8, v_s2);

                s12 += ar[0] + ar[1];
                s11 += ar[2] + ar[3];
                s22 += ar[4] + ar[5];
                s1 += ar[6] + ar[7];
                s2 += ar[8] + ar[9];
            }
            #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_NEON
            float32x4_t v_result = vdupq_n_f32(0.0f);
            for( ; j <= len - 4; j += 4 )
                v_result = vaddq_f32(v_result, vminq_f32(vld1q_f32(h1 + j), vld1q_f32(h2 + j)));
            float CV_DECL_ALIGNED(16) ar[4];
            vst1q_f32(ar, v_result);
            result += ar[0] + ar[1] + ar[2] + ar[3];
            #elif CV_SSE2
            if (haveSIMD)
            {
                __m128d v_result = _mm_setzero_pd();
                for ( ; j <= len - 4; j += 4)
                {
                    __m128 v_src = _mm_min_ps(_mm_loadu_ps(h1 + j),
                                              _mm_loadu_ps(h2 + j));
                    v_result = _mm_add_pd(v_result, _mm_cvtps_pd(v_src));
                    v_src = _mm_castsi128_ps(_mm_srli_si128(_mm_castps_si128(v_src), 8));
                    v_result = _mm_add_pd(v_result, _mm_cvtps_pd(v_src));
                }

                double CV_DECL_ALIGNED(16) ar[2];
                _mm_store_pd(ar, v_result);
                result += ar[0] + ar[1];
            }
            #endif
            for( ; j < len; j++ )
                result += std::min(h1[j], h2[j]);
        }
        else if( method == CV_COMP_BHATTACHARYYA )
        {
            #if CV_SSE2
            if (haveSIMD)
            {
                __m128d v_s1 = _mm_setzero_pd(), v_s2 = v_s1, v_result = v_s1;
                for ( ; j <= len - 4; j += 4)
                {
                    __m128 v_a = _mm_loadu_ps(h1 + j);
                    __m128 v_b = _mm_loadu_ps(h2 + j);

                    __m128d v_ad = _mm_cvtps_pd(v_a);
                    __m128d v_bd = _mm_cvtps_pd(v_b);
                    v_s1 = _mm_add_pd(v_s1, v_ad);
                    v_s2 = _mm_add_pd(v_s2, v_bd);
                    v_result = _mm_add_pd(v_result, _mm_sqrt_pd(_mm_mul_pd(v_ad, v_bd)));

                    v_ad = _mm_cvtps_pd(_mm_castsi128_ps(_mm_srli_si128(_mm_castps_si128(v_a), 8)));
                    v_bd = _mm_cvtps_pd(_mm_castsi128_ps(_mm_srli_si128(_mm_castps_si128(v_b), 8)));
                    v_s1 = _mm_add_pd(v_s1, v_ad);
                    v_s2 = _mm_add_pd(v_s2, v_bd);
                    v_result = _mm_add_pd(v_result, _mm_sqrt_pd(_mm_mul_pd(v_ad, v_bd)));
                }

                double CV_DECL_ALIGNED(16) ar[6];
                _mm_store_pd(ar, v_s1);
                _mm_store_pd(ar + 2, v_s2);
                _mm_store_pd(ar + 4, v_result);
                s1 += ar[0] + ar[1];
                s2 += ar[2] + ar[3];
                result += ar[4] + ar[5];
            }
            #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_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;
}

/*【compareHist ( )源代码】***********************************************************
 * @Version：OpenCV 3.0.0（Opnencv2和Opnencv3差别不大，Linux和PC的对应版本源码完全一样，均在对应的安装目录下）  
 * @源码路径：…\opencv\sources\modules\imgproc\src\histogram.cpp
 * @起始行数： 2478行   
********************************************************************************/
double cv::compareHist( const SparseMat& H1, const SparseMat& H2, int method )
{
    double result = 0;
    int i, dims = H1.dims();

    CV_Assert( dims > 0 && dims == H2.dims() && H1.type() == H2.type() && H1.type() == CV_32F );
    for( i = 0; i < dims; i++ )
        CV_Assert( H1.size(i) == H2.size(i) );

    const SparseMat *PH1 = &H1, *PH2 = &H2;
    if( PH1->nzcount() > PH2->nzcount() && method != CV_COMP_CHISQR && method != CV_COMP_CHISQR_ALT && method != CV_COMP_KL_DIV )
        std::swap(PH1, PH2);

    SparseMatConstIterator it = PH1->begin();
    int N1 = (int)PH1->nzcount(), N2 = (int)PH2->nzcount();

    if( (method == CV_COMP_CHISQR) || (method == CV_COMP_CHISQR_ALT) )
    {
        for( i = 0; i < N1; i++, ++it )
        {
            float v1 = it.value<float>();
            const SparseMat::Node* node = it.node();
            float v2 = PH2->value<float>(node->idx, (size_t*)&node->hashval);
            double a = v1 - v2;
            double b = (method == CV_COMP_CHISQR) ? v1 : v1 + v2;
            if( fabs(b) > DBL_EPSILON )
                result += a*a/b;
        }
    }
    else if( method == CV_COMP_CORREL )
    {
        double s1 = 0, s2 = 0, s11 = 0, s12 = 0, s22 = 0;

        for( i = 0; i < N1; i++, ++it )
        {
            double v1 = it.value<float>();
            const SparseMat::Node* node = it.node();
            s12 += v1*PH2->value<float>(node->idx, (size_t*)&node->hashval);
            s1 += v1;
            s11 += v1*v1;
        }

        it = PH2->begin();
        for( i = 0; i < N2; i++, ++it )
        {
            double v2 = it.value<float>();
            s2 += v2;
            s22 += v2*v2;
        }

        size_t total = 1;
        for( i = 0; i < H1.dims(); i++ )
            total *= H1.size(i);
        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_INTERSECT )
    {
        for( i = 0; i < N1; i++, ++it )
        {
            float v1 = it.value<float>();
            const SparseMat::Node* node = it.node();
            float v2 = PH2->value<float>(node->idx, (size_t*)&node->hashval);
            if( v2 )
                result += std::min(v1, v2);
        }
    }
    else if( method == CV_COMP_BHATTACHARYYA )
    {
        double s1 = 0, s2 = 0;

        for( i = 0; i < N1; i++, ++it )
        {
            double v1 = it.value<float>();
            const SparseMat::Node* node = it.node();
            double v2 = PH2->value<float>(node->idx, (size_t*)&node->hashval);
            result += std::sqrt(v1*v2);
            s1 += v1;
        }

        it = PH2->begin();
        for( i = 0; i < N2; i++, ++it )
            s2 += it.value<float>();

        s1 *= s2;
        s1 = fabs(s1) > FLT_EPSILON ? 1./std::sqrt(s1) : 1.;
        result = std::sqrt(std::max(1. - result*s1, 0.));
    }
    else if( method == CV_COMP_KL_DIV )
    {
        for( i = 0; i < N1; i++, ++it )
        {
            double v1 = it.value<float>();
            const SparseMat::Node* node = it.node();
            double v2 = PH2->value<float>(node->idx, (size_t*)&node->hashval);
            if( !v2 )
                v2 = 1e-10;
            result += v1 * std::log( v1 / v2 );
        }
    }
    else
        CV_Error( CV_StsBadArg, "Unknown comparison method" );

    if( method == CV_COMP_CHISQR_ALT )
        result *= 2;

    return result;
}