2019-07-08 15:38:02 Osean_li 阅读数 145
  • Matplotlib 数据分析可视化

    数据分析三剑客,NumPy、Pandas、Matplotlib,本课程是对Matplotlib的讲解,Matplotlib可以是分析的数据可视化,可以更直观的查看数据分析的结果,本课程独辟蹊径,不光教大家如何绘图,例如:饼图、柱状图、条形图、直方图等,而且深入剖析的Matplotlib的绘图原理,例如:如何操作图片,如何加载本地数据等。

    15306 人正在学习 去看看 郭宏志

引言

直方图广泛运用于很多计算机视觉运用当中。
参考整理的很好

什么是图像直方图?

通过标记帧与帧之间显著的边缘和颜色的统计变化,来检测视频中场景的变化。在每个兴趣点设置一个有相近特征的直方图所构成 “标签”,用以确定图像中的兴趣点。边缘,色彩,角度等直方图构成了可以被传递给目标识别分类器的一个通用特征类型

直方图就是数据进行统计的一种方法

  • 且具有图像平移、旋转、缩放不变性等众多优点。直方图在进行图像计算处理时代价较小,所以经常用于图像处理!
  • 直方图是一种对数据分布情况的图形表示,是一种二位统计图表,它的两个坐标分别是统计样本和该样本对应的某个属性的度量
  • 直方图是图像中像素强度分布的图形表达式
  • 它统计了每一个强度值所具有的像素个数

如何读懂直方图?

直方图的显示方式是左暗又亮,左边用于描述图像的暗度,右边用于描述图像的亮度
在这里插入图片描述

2019-09-09 15:48:53 weixin_41709536 阅读数 30
  • Matplotlib 数据分析可视化

    数据分析三剑客,NumPy、Pandas、Matplotlib,本课程是对Matplotlib的讲解,Matplotlib可以是分析的数据可视化,可以更直观的查看数据分析的结果,本课程独辟蹊径,不光教大家如何绘图,例如:饼图、柱状图、条形图、直方图等,而且深入剖析的Matplotlib的绘图原理,例如:如何操作图片,如何加载本地数据等。

    15306 人正在学习 去看看 郭宏志

知识点

图像直方图的解释
图像直方图是图像像素值的统计学特征、计算代价较小,具有图像平移、旋转、缩放不变性等众多优点,广泛地应用于图像处理的各个领域,特别是灰度图像的阈值分割、基于颜色的图像检索以及图像分类、反向投影跟踪。常见的分为

  • 灰度直方图
  • 颜色直方图

Bins是指直方图的大小范围, 对于像素值取值在0~255之间的,最少有256个bin,此外还可以有16、32、48、128等,256除以bin的大小应该是整数倍。

OpenCV中相关API
calcHist(&bgr_plane[0], 1, 0, Mat(), b_hist, 1, bins, ranges);
cv.calcHist([image], [i], None, [256], [0, 256])

python代码

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt

#灰度直方图
def custom_hist(gray):
    h, w = gray.shape
    hist = np.zeros([256], dtype=np.int32)
    for row in range(h):
        for col in range(w):
            pv = gray[row, col]
            hist[pv] += 1

    y_pos = np.arange(0, 256, 1, dtype=np.int32)
    plt.bar(y_pos, hist, align='center', color='r', alpha=0.5)
    plt.xticks(y_pos, y_pos)
    plt.ylabel('Frequency')
    plt.title('Histogram')

    # plt.plot(hist, color='r')
    # plt.xlim([0, 256])
    plt.show()

#颜色直方图
def image_hist(image):
    cv.imshow("input", image)
    color = ('blue', 'green', 'red')
    for i, color in enumerate(color):
        #三通道分离
        hist = cv.calcHist([image], [i], None, [256], [0, 256])
        plt.plot(hist, color=color)
        plt.xlim([0, 256])
    plt.show()


src = cv.imread("C:/Users/qqxd/Desktop/opencvcode/images/flower.png")
cv.namedWindow("input", cv.WINDOW_AUTOSIZE)
gray = cv.cvtColor(src, cv.COLOR_BGR2GRAY)
cv.imshow("input", gray)
custom_hist(gray)
image_hist(src)
cv.waitKey(0)
cv.destroyAllWindows()

c++代码

#include <opencv2/opencv.hpp>
#include <iostream>

using namespace cv;
using namespace std;

const int bins = 256;
Mat src;
const char* winTitle = "input image";

void drawHistogram(Mat &image);
int main(int argc, char** argv) {
	src = imread("C:/Users/qqxd/Desktop/opencvcode/images/flower.png");
	if (src.empty()) {
		printf("could not load image...\n");
		return 0;
	}
	namedWindow(winTitle, WINDOW_AUTOSIZE);
	imshow(winTitle, src);
	drawHistogram(src);
	waitKey(0);
	return 0;
}

void drawHistogram(Mat &image) {
	// 定义参数变量
	const int channels[1] = { 0 };
	const int bins[1] = { 256 };
	float hranges[2] = { 0,255 };
	const float* ranges[1] = { hranges };
	int dims = image.channels();
	if (dims == 3) {
		vector<Mat> bgr_plane;
		split(src, bgr_plane);
		Mat b_hist;
		Mat g_hist;
		Mat r_hist;
		// 计算Blue, Green, Red通道的直方图
		calcHist(&bgr_plane[0], 1, 0, Mat(), b_hist, 1, bins, ranges);
		calcHist(&bgr_plane[1], 1, 0, Mat(), g_hist, 1, bins, ranges);
		calcHist(&bgr_plane[2], 1, 0, Mat(), r_hist, 1, bins, ranges);
		// 显示直方图
		int hist_w = 512;
		int hist_h = 400;
		int bin_w = cvRound((double)hist_w / bins[0]);
		Mat histImage = Mat::zeros(hist_h, hist_w, CV_8UC3);
		// 归一化直方图数据
		normalize(b_hist, b_hist, 0, histImage.rows, NORM_MINMAX, -1, Mat());
		normalize(g_hist, g_hist, 0, histImage.rows, NORM_MINMAX, -1, Mat());
		normalize(r_hist, r_hist, 0, histImage.rows, NORM_MINMAX, -1, Mat());
		// 绘制直方图曲线
		for (int i = 1; i < bins[0]; i++) {
			line(histImage, Point(bin_w*(i - 1), hist_h - cvRound(b_hist.at<float>(i - 1))),
				Point(bin_w*(i), hist_h - cvRound(b_hist.at<float>(i))), Scalar(255, 0, 0), 2, 8, 0);
			line(histImage, Point(bin_w*(i - 1), hist_h - cvRound(g_hist.at<float>(i - 1))),
				Point(bin_w*(i), hist_h - cvRound(g_hist.at<float>(i))), Scalar(0, 255, 0), 2, 8, 0);
			line(histImage, Point(bin_w*(i - 1), hist_h - cvRound(r_hist.at<float>(i - 1))),
				Point(bin_w*(i), hist_h - cvRound(r_hist.at<float>(i))), Scalar(0, 0, 255), 2, 8, 0);
		}
		// 显示直方图
		namedWindow("Histogram Demo", WINDOW_AUTOSIZE);
		imshow("Histogram Demo", histImage);
	}
	else {
		Mat hist;
		// 计算Blue, Green, Red通道的直方图
		calcHist(&image, 1, 0, Mat(), hist, 1, bins, ranges);
		// 显示直方图
		int hist_w = 512;
		int hist_h = 400;
		int bin_w = cvRound((double)hist_w / bins[0]);
		Mat histImage = Mat::zeros(hist_h, hist_w, CV_8UC3);
		// 归一化直方图数据
		normalize(hist, hist, 0, histImage.rows, NORM_MINMAX, -1, Mat());
		// 绘制直方图曲线
		for (int i = 1; i < bins[0]; i++) {
			line(histImage, Point(bin_w*(i - 1), hist_h - cvRound(hist.at<float>(i - 1))),
				Point(bin_w*(i), hist_h - cvRound(hist.at<float>(i))), Scalar(0, 255, 0), 2, 8, 0);
		}
		// 显示直方图
		namedWindow("Histogram Demo", WINDOW_AUTOSIZE);
		imshow("Histogram Demo", histImage);
	}
}

void showHistogram() {
	// 三通道分离
	vector<Mat> bgr_plane;
	split(src, bgr_plane);
	// 定义参数变量
	const int channels[1] = { 0 };
	const int bins[1] = { 256 };
	float hranges[2] = { 0,255 };
	const float* ranges[1] = { hranges };
	Mat b_hist;
	Mat g_hist;
	Mat r_hist;
	// 计算Blue, Green, Red通道的直方图
	calcHist(&bgr_plane[0], 1, 0, Mat(), b_hist, 1, bins, ranges);
	calcHist(&bgr_plane[1], 1, 0, Mat(), g_hist, 1, bins, ranges);
	calcHist(&bgr_plane[2], 1, 0, Mat(), r_hist, 1, bins, ranges);
	// 显示直方图
	int hist_w = 512;
	int hist_h = 400;
	int bin_w = cvRound((double)hist_w / bins[0]);
	Mat histImage = Mat::zeros(hist_h, hist_w, CV_8UC3);
	// 归一化直方图数据
	normalize(b_hist, b_hist, 0, histImage.rows, NORM_MINMAX, -1, Mat());
	normalize(g_hist, g_hist, 0, histImage.rows, NORM_MINMAX, -1, Mat());
	normalize(r_hist, r_hist, 0, histImage.rows, NORM_MINMAX, -1, Mat());
	// 绘制直方图曲线
	for (int i = 1; i < bins[0]; i++) {
		line(histImage, Point(bin_w*(i - 1), hist_h - cvRound(b_hist.at<float>(i - 1))),
			Point(bin_w*(i), hist_h - cvRound(b_hist.at<float>(i))), Scalar(255, 0, 0), 2, 8, 0);
		line(histImage, Point(bin_w*(i - 1), hist_h - cvRound(g_hist.at<float>(i - 1))),
			Point(bin_w*(i), hist_h - cvRound(g_hist.at<float>(i))), Scalar(0, 255, 0), 2, 8, 0);
		line(histImage, Point(bin_w*(i - 1), hist_h - cvRound(r_hist.at<float>(i - 1))),
			Point(bin_w*(i), hist_h - cvRound(r_hist.at<float>(i))), Scalar(0, 0, 255), 2, 8, 0);
	}
	// 显示直方图
	namedWindow("Histogram Demo", WINDOW_AUTOSIZE);
	imshow("Histogram Demo", histImage);
}

python运行结果如下:
在这里插入图片描述
在这里插入图片描述
c++运行结果如下:
在这里插入图片描述
在这里插入图片描述

注:之所以python与c++得到的三通道颜色直方图大小不一致,因为c++版本做了归一化处理。

2019-01-16 21:12:31 qq_42887760 阅读数 174
  • Matplotlib 数据分析可视化

    数据分析三剑客,NumPy、Pandas、Matplotlib,本课程是对Matplotlib的讲解,Matplotlib可以是分析的数据可视化,可以更直观的查看数据分析的结果,本课程独辟蹊径,不光教大家如何绘图,例如:饼图、柱状图、条形图、直方图等,而且深入剖析的Matplotlib的绘图原理,例如:如何操作图片,如何加载本地数据等。

    15306 人正在学习 去看看 郭宏志

直方图概念

在这里插入图片描述
上述直方图概念是基于图像像素值,其实对图像梯度、每个像素的角度、等一切图像的属性值,我们都可以建立直方图。这个才是直方图的概念真正意义,不过是基于图像像素灰度直方图是最常见的。
直方图最常见的几个属性:

  • dims:要收集数据的参数数量。 在我们的示例中,dims = 1,因为我们只计算每个像素的强度值(在灰度图像中)。
  • bin:它是每个暗淡的细分数量。 在我们的示例中,bin = 16
  • range:要测量的值的限制。 在这种情况下:范围= [0,255]

API学习

  • split(// 把多通道图像分为多个单通道图像
    const Mat &src, //输入图像
    Mat* mvbegin)// 输出的通道图像数组

  • calcHist(
    const Mat* images,//输入图像指针
    int images,// 图像数目
    const int* channels,// 通道数,要计算的通道数的下标,可以传一个数组 {0, 1} 表示计算第0通道与第1通道的直方图,此数组长度要与histsize ranges 数组长度一致
    InputArray mask,//输入mask,可选。如有,则表示只计算mask元素值为255的位置的直方图
    OutputArray hist,//输出的直方图数据
    int dims,// 维数
    const int* histsize,// 直方图级数 ,对应 bins
    const float* ranges,// 值域范围
    bool uniform,// true by default 是否归一化到 0-1 之间
    bool accumulate// false by defaut
    )

程序代码

#include<opencv2/opencv.hpp>
#include<iostream>
using namespace std;
using namespace cv;

int main(int argc, char** argv)
{
    Mat src;
	// 1. 加载源图像
    src=imread("E:/Experiment/OpenCV/Pictures/dog2.jpg");
	imshow("input Image",src);

    // 2. 在R、G、B平面中分离源图像,把多通道图像分为多个单通道图像。使用OpenCV函数cv::split。
    vector<Mat> bgr_planes;
    split(src, bgr_planes);// 把多通道图像分为多个单通道图像

    printf("channels=%d\n", bgr_planes.size());//3通道,所以size也是3
    imshow("channels_b", bgr_planes[0]);
    imshow("channels_g", bgr_planes[1]);
    imshow("channels_r", bgr_planes[2]);

    // 3. 现在我们准备开始为每个平面配置直方图。 由于我们正在使用B,G和R平面,我们知道我们的值将在区间[0,255]范围内
    int histBins = 256;//建立箱数(5,10 ......)
    float range[] = { 0, 256 };//设置值的范围(在0到255之间)
    const float * histRanges = range;//注意:函数形参 float ** 与 const float ** 是两种不同数据类型。
	bool uniform = true, accumulate = false;//我们希望我们的箱子具有相同的尺寸(均匀)并在开头清除直方图
    Mat b_hist, g_hist, r_hist;//calcHist计算出来的Mat中元素的最大值可能上几千,所以最好归一化后再绘制直方图
    //使用OpenCV函数cv::calcHist计算直方图:
	calcHist(&bgr_planes[0], 1, 0, Mat(), b_hist, 1, &histBins, &histRanges, uniform, accumulate);//计算直方图
    calcHist(&bgr_planes[1], 1, 0, Mat(), g_hist, 1, &histBins, &histRanges, uniform, accumulate);
    calcHist(&bgr_planes[2], 1, 0, Mat(), r_hist, 1, &histBins, &histRanges, uniform, accumulate);

    // 4. 归一化
    int hist_cols = 400;
    int hist_rows = 512;
    int bin_w = hist_rows / histBins;
    /*
        normalize(  // normalize函数作用为 归一化数据
            InputArray src, // 输入数组
            InputOutputArray dst, // 输出数组,数组的大小和原数组一致
            double alpha = 1, // 1,用来规范值,2.规范范围,并且是下限
            double beta = 0, // 只用来规范范围并且是上限
            int norm_type = NORM_L2, // 归一化选择的数学公式类型
            int dtype = -1, // 当为负,输出在大小深度通道数都等于输入,当为正,输出只在深度与输如不同,不同的地方由dtype决定
            InputArray mask = noArray() // 掩码。选择感兴趣区域,选定后只能对该区域进行操作
        );

        值归一化举例说明:   参考博客 https://blog.csdn.net/cosmispower/article/details/64457406
            src={10,23,71},参数 beta 设为0,值的范围为 0-参数alpha 设置的值,例子中alpha设置为1
            NORM_L1运算后得到   dst={0.096,0.221,0.683}       NORM_* 的归一化公式参考上述博客
            NORM_INF运算后得到  dst={0.141,0.324,1}
            NORM_L2运算后得到   dst={0.133,0.307,0.947}
            NORM_MINMAX运算得到 dst={0,0.377,1}             P = Ak / (max(Ai)-min(Ai))  Ak等于最大最小Ai时,不按此公式计算,P直接等于1, 0

        范围归一化时,beta必有值不等于0,范围为 alpha-beta ,alpha为下限(可为0也可非0),beta为上限
    */
	//请注意,在绘制之前,我们首先对直方图进行cv :: normalize,使其值落在输入参数指示的范围内:
    normalize(b_hist, b_hist, 0, hist_cols, NORM_MINMAX, -1, Mat());//b_hist中元素的值转换到 0-hist_cols 之间
    normalize(g_hist, g_hist, 0, hist_cols, NORM_MINMAX, -1, Mat());
    normalize(r_hist, r_hist, 0, hist_cols, NORM_MINMAX, -1, Mat());//传参 0, hist_cols 或 hist_cols, 0 结果一致

    // 5. 绘制直方图
    Mat histImage(hist_rows, hist_cols, CV_8UC3, Scalar(0, 0, 0));
    for (int i = 1; i < histBins; i++)
    {
        // cvRound 四舍五入,返回整型值
        line(histImage, Point((i - 1)*bin_w, hist_cols - cvRound(b_hist.at<float>(i - 1))),
            Point(i*bin_w, hist_cols - cvRound(b_hist.at<float>(i))), Scalar(255, 0, 0), 2, LINE_AA);
        line(histImage, Point((i - 1)*bin_w, hist_cols - cvRound(g_hist.at<float>(i - 1))),
            Point(i*bin_w, hist_cols - cvRound(g_hist.at<float>(i))), Scalar(0, 255, 0), 2, LINE_AA);
        line(histImage, Point((i - 1)*bin_w, hist_cols - cvRound(r_hist.at<float>(i - 1))),
            Point(i*bin_w, hist_cols - cvRound(r_hist.at<float>(i))), Scalar(0, 0, 255), 2, LINE_AA);
    }
	// 6. 最后,我们显示直方图并等待用户退出:
    imshow("histImage", histImage);

    waitKey(0);
	return 0;
}

运行结果

在这里插入图片描述

参考代码

  1. https://blog.csdn.net/huanghuangjin/article/details/81173895
  2. https://blog.csdn.net/LYKymy/article/details/83189170
  3. https://blog.csdn.net/u011574296/article/details/73381576
2012-08-14 17:07:47 masikkk 阅读数 7327
  • Matplotlib 数据分析可视化

    数据分析三剑客,NumPy、Pandas、Matplotlib,本课程是对Matplotlib的讲解,Matplotlib可以是分析的数据可视化,可以更直观的查看数据分析的结果,本课程独辟蹊径,不光教大家如何绘图,例如:饼图、柱状图、条形图、直方图等,而且深入剖析的Matplotlib的绘图原理,例如:如何操作图片,如何加载本地数据等。

    15306 人正在学习 去看看 郭宏志

    直方图就是对数据进行统计,将统计值组织到一系列事先定义好的bin(直方图中的柱子)中。bin中的数值是从数据中计算出的特征的统计量,这些数据可以是诸如梯度、方向、色彩或任何其他特征。无论如何,直方图获得的是数据分布的统计图。

灰度图像的直方图的性质:

(1) 直方图是一幅图像中各像素灰度出现频次的统计结果,它只反映图像中不同灰度值出现的次数,而没反映某一灰度所在的位置。也就是说,它只包含了该图像的某一灰度像素出现的概率,而丢失了其所在的位置信息。

(2) 任一幅图像,都有惟一确定一幅与它对应的直方图,但不同的图像可能有相同的直方图。即图像与直方图之间是多对一的映射关系。

(3) 由于直方图是对具有相同灰度值的像素统计得到的,因此,一幅图像各子区的直方图之和就等于该图像全图的直方图。

通过直方图均衡化进行图像增强

       直方图均衡化是灰度变换的一个重要应用,广泛应用在图像增强处理中,它是以累计分布函数变换为基础的直方图修正法,可以产生一幅灰度级分布具有均匀概率密度的图像,扩展了像素的取值动态范围。许多图像的灰度值是非均匀分布的,其中灰度值集中在一个小区间内的图像是很常见的,直方图均衡化是一种通过重新均匀地分布各灰度值来增强图像对比度的方法,经过直方图均衡化的图像对二值化阈值选取十分有利。一般来说,直方图修正能提高图像的主观质量,因此在处理艺术图像时非常有用。直方图均衡化处理的中心思想是把原始图像的灰度直方图从比较集中的某个灰度区间变成在全部灰度范围内的均匀分布。

OpenCV中提供了现成的直方图均衡化函数:void cvEqualizeHist( const CvArr* src, CvArr* dst );

实现如下:

  IplImage * pImage = cvLoadImage("长虹大厦20120810.jpg", 1);//原始图像
	IplImage * pMergeImage = cvCreateImage(cvGetSize(pImage),pImage->depth,pImage->nChannels);//处理后的图像
	IplImage * pImageChannel[4] = {0,0,0,0};//分别保存4个通道的灰度图像
	//创建各个灰度图像
	for(int i=0; i<pImage->nChannels; i++)
		pImageChannel[i] = cvCreateImage(cvGetSize(pImage),pImage->depth,1);
	//分割通道
	cvSplit(pImage,pImageChannel[0],pImageChannel[1],pImageChannel[2],pImageChannel[3]);
	//对每个信道分别做直方图均衡化
	for(i=0; i<pImage->nChannels; i++)
		cvEqualizeHist(pImageChannel[i],pImageChannel[i]);
	//合并通道
	cvMerge(pImageChannel[0],pImageChannel[1],pImageChannel[2],pImageChannel[3],pMergeImage);


	cvNamedWindow("原始图像",1);
	cvNamedWindow("处理后的图像",1);
	cvShowImage("原始图像",pImage);
	cvShowImage("处理后的图像",pMergeImage);

处理前:

处理前图像的直方图:

处理后:

处理后图像的直方图:

 

可以看到处理后图像的直方图分布扩散了,说明色彩分布更宽广了,不集中在一起了。

 

源码下载:

http://download.csdn.net/detail/masikkk/4499482

 

 

2018-03-03 11:32:34 u013162035 阅读数 333
  • Matplotlib 数据分析可视化

    数据分析三剑客,NumPy、Pandas、Matplotlib,本课程是对Matplotlib的讲解,Matplotlib可以是分析的数据可视化,可以更直观的查看数据分析的结果,本课程独辟蹊径,不光教大家如何绘图,例如:饼图、柱状图、条形图、直方图等,而且深入剖析的Matplotlib的绘图原理,例如:如何操作图片,如何加载本地数据等。

    15306 人正在学习 去看看 郭宏志

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;
}

3.4.3直方图对比实例

代码参看附件【demo1】。

这里写图片描述

图1

这里写图片描述

图2

这里写图片描述

图3

参考:
中文
英文

本章参考代码

点击进入

没有更多推荐了,返回首页