直方图均衡化算法、直方图匹配算法 C++ 代码

这两天一直在研究匀光匀色算法才了解到了直方图匹配算法，想要了解这个算法又要先了解直方图均衡化算法，通过网上查找了很多资料，没有现成C++代码，经过仔细思考和实验后大概复现了该算法。特此记录，以备查阅
参考链接如下：
1、匀光匀色–直方图匹配算法实现与应用
2、基于OpenCV的直方图匹配
3、直方图均衡化的数学原理

先介绍一下基本概念

测试数据

原始数据

底图影像

直方图匹配匀光匀色效果

从上图可以看出原始影像与模板颜色不一致，经过直方图匹配匀光匀色后颜色一致性较好。

直方图均衡算法

先看看直方图均衡算法原理知识：

先在根据上面的算法步骤，计算一幅影像的直方图均衡影像

1、获取一张影像Mat的直方图，即计算hs

void getImgHistogram(Mat image, map<int, double> & map)
{//获取一张影像Mat的直方图，即计算hsint w = image.cols;int h = image.rows;double totalPixel = w * h;double rate = 1 / (image.channels() * totalPixel);//所有像素分之一for (int i = 0; i < 256; i++) {//map.put(i, 0.0);// 通过循环，往集合里面填充0~255个位置,初始值都为0map[i] = 0.0;}//分别统计图像上0~255上分布总数for (int row = 0; row < h; row++) {uchar* uc_pixel = image.data + row * image.step;for (int col = 0; col < w; col++) {//uc_pixel[0] = 255 - uc_pixel[0];//uc_pixel[1] = 255 - uc_pixel[1];//uc_pixel[2] = 255 - uc_pixel[2];for (int BandID = 0; BandID < image.channels(); BandID++){map[uc_pixel[BandID]] += rate;}uc_pixel += image.channels();//map[uc_pixel[0]] += rate;//map[uc_pixel[1]] += rate;//map[uc_pixel[2]] += rate;//uc_pixel += 3;}}return;
}

2、计算累计直方图，即hp

 double SrcSum = 0;map<int, double> SrcImg_AccMap;//计算累计直方图，即hpfor (int i = 0; i < 256; i++){SrcSum += SrcImg_map[i];SrcImg_AccMap[i] = SrcSum;}

3、对直方图中所有像素乘以255即可得到最终直方图均衡的影像，不单独列出代码，最后看所有代码里面可以找到。

所有代码如下：

#include <opencv.hpp>
#include <map>
#include <vector>using namespace std;
using namespace cv;
void getImgHistogram(Mat image, map<int, double> & map);
Mat GetEquHistogramMat(Mat SrcImage, map<int, double> SrcImg_AccMap);//获取均衡化的影像Matint main()
{const char* SrcPath = "D:/datas/ColorMap/back/Test.png";const char* outPutPath = "D:/datas/ColorMap/back/Test_ZFTJH.png";Mat SrcImage = imread(SrcPath);map<int, double> SrcImg_map;//定义直方图map,直方图可以使用map进行定义，也可以使用一个一维数组定义getImgHistogram(SrcImage, SrcImg_map);//获取影像直方图//计算累计直方图double SrcSum = 0;map<int, double> SrcImg_AccMap;//累计直方图for (int i = 0; i < 256; i++){SrcSum += SrcImg_map[i];SrcImg_AccMap[i] = SrcSum;}Mat SrcImage_Equ = GetEquHistogramMat(SrcImage, SrcImg_AccMap);//获取直方图均衡化的Matimwrite(outPutPath, SrcImage_Equ);printf("Success save img:%s\n", outPutPath);imshow("SrcImage", SrcImage);imshow("SrcImage_Equ", SrcImage_Equ);waitKey(0);return 0;
}void getImgHistogram(Mat image, map<int, double> & map)
{//获取影像直方图int w = image.cols;int h = image.rows;double totalPixel = w * h;double rate = 1 / (image.channels() * totalPixel);//所有像素分之一for (int i = 0; i < 256; i++) {// 通过循环，往集合里面填充0~255个位置,初始值都为0map[i] = 0.0;}//分别统计图像上0~255上分布总数for (int row = 0; row < h; row++) {uchar* uc_pixel = image.data + row * image.step;for (int col = 0; col < w; col++) {//每个波段像素for (int BandID = 0; BandID < image.channels(); BandID++){map[uc_pixel[BandID]] += rate;}uc_pixel += image.channels();//map[uc_pixel[0]] += rate;//map[uc_pixel[1]] += rate;//map[uc_pixel[2]] += rate;//uc_pixel += 3;}}return;
}Mat GetEquHistogramMat(Mat SrcImage, map<int, double> SrcImg_AccMap)//获取直方图均衡化的Mat
{//累计直方图，然后再每个直方图像素乘以255即得到均衡的直方图Mat resultImg = SrcImage.clone();//均衡化的src图for (int row = 0; row < resultImg.rows; row++) {uchar* uc_pixel = resultImg.data + row * resultImg.step;for (int col = 0; col < resultImg.cols; col++) {for (int BandID = 0; BandID < SrcImage.channels(); BandID++){uc_pixel[BandID]=255* SrcImg_AccMap[uc_pixel[BandID]];}uc_pixel += SrcImage.channels();//uc_pixel[0] = 255 * SrcImg_AccMap[uc_pixel[0]];//uc_pixel[1] = 255 * SrcImg_AccMap[uc_pixel[1]];//uc_pixel[2] = 255 * SrcImg_AccMap[uc_pixel[2]];//uc_pixel += 3;}}return resultImg;}

直方图匹配算法

直方图匹配算法是在直方图均衡算法基础上进行的衍生
直方图匹配：将一张图片的直方图匹配到目标图上，使两张图的视觉感觉接近
关于理论公式有很多，如果喜欢看理论知识，这里还有一篇博客：直方图匹配的数学原理。写了很多，但是公式太复杂了，看不是很懂。直到看了下面这幅图我才懂了直方图匹配的流程步骤

上图说的是待匹配影像A和参考底图B，分别对其进行直方图均衡化（说简单点就是累计直方图乘以255），输出的直方图均衡的影像称为：A_Equ、B_Equ影像，看A_Equ和B_Equ那个像素值距离最近，则将B上的该直方图值映射给A,形成一个映射表。例如通过上面最下面这副图可以得到映射表（映射表可以用map<int ,uchar>数据格式表达，或者uchar*的一维数组来表达，以数组下标表示直方图值），为了值观表达，上图的映射表{0：0，1：1，2：1，3：2，4：4}。
我认为上图映射表中其实有点问题（该表我是从参考链接三中搬过来的）：A图的像素3，转换成均衡化影像像素值为18，18这个值在B_Equ中距离15是3（fabs(18-15)），距离20是2（fabs(18-20)），依照取出B_Equ最近的像素作为映射像素，所以该值应该映射到B_Equ中的像素20，对应B中的的是3，所以映射表中3映射成3，而不是2。
最后一步就是根据引射表{0：0，1：1，2：1，3：2，4：4}进行修正A影像得到的就是右下角的直方图匹配后的影像，将该影像存出来则是最终结果。

注意这是直方图匹配大概流程，如果想要颜色与底图颜色更接近，则需要对待匹配影像的每个波段都与底图影像每个波段做匹配，效果则最好。

上代码：


#include <opencv.hpp>
#include <map>
#include <vector>using namespace std;
using namespace cv;
void getSigBandImgHistogram(Mat SrcImage, map<int, double> &SrcImg_SigMap, int BandId);int main(int argc, char* argv[])
{//const char* SrcPath = "D:/datas/test/ColorCorrection/Src.png";//const char* BaseMap = "D:/datas/test/ColorCorrection/baseMap.png";//const char* outPutPath= "D:/datas/test/ColorCorrection/outPut.png";const char* SrcPath = argv[1];const char* BaseMap = argv[2];const char* outPutPath = argv[3];int IsDislay = atoi(argv[4]);//是否显示匹配后影像效果，0为不显示，其他的都显示Mat SrcImage = imread(SrcPath);Mat BaseMapImg = imread(BaseMap);cout << SrcImage.channels() << endl;vector<uchar*>resultMultMap;//存储每个波段最终的映射关系for (int BandId = 0; BandId < SrcImage.channels(); BandId++){//单独对影像每个波段进行直方图均衡，然后找到最终的映射关系printf("Start Get BandID:%d Histogram!\n", BandId);uchar* SigResultBandMap = new uchar[256];//首先定义一下单个波段的映射结果//1、获取单波段直方图map<int, double> SrcImg_SigMap;//map<int, double> BaseMapImg_SigMap;getSigBandImgHistogram(SrcImage, SrcImg_SigMap, BandId);getSigBandImgHistogram(BaseMapImg, BaseMapImg_SigMap, BandId);//计算累计直方图double BaseSum_SigMap = 0;double SrcSum_SigMap = 0;map<int, double> SrcImg_AccSigMap;//累计直方图map<int, double> BaseMapImg_AccSigMap;//累计直方图for (int i = 0; i < 256; i++){SrcSum_SigMap += SrcImg_SigMap[i];SrcImg_AccSigMap[i] = SrcSum_SigMap;BaseSum_SigMap += BaseMapImg_SigMap[i];BaseMapImg_AccSigMap[i] = BaseSum_SigMap;}printf("Start Build resultMap :%d!\n", BandId);for (int i = 0; i < 256; i++){double MinValue = 255;//最小值初始化for (int j = 0; j < 256; j++){double SubVelue = fabs(SrcImg_AccSigMap[i] * 255 - BaseMapImg_AccSigMap[j] * 255);//取绝对值if (MinValue > SubVelue){//找到最靠近直方图均衡化的值，作为映射值MinValue = SubVelue;SigResultBandMap[i] = j;}}}resultMultMap.push_back(SigResultBandMap);}//最后进行直方图像素替换Mat resultImg = SrcImage.clone();//均衡化的src图for (int row = 0; row < resultImg.rows; row++) {uchar* uc_pixel = resultImg.data + row * resultImg.step;for (int col = 0; col < resultImg.cols; col++){for (int BID = 0; BID < SrcImage.channels(); BID++){uc_pixel[BID] = resultMultMap[BID][uc_pixel[BID]];}//uc_pixel[0] = resultMap[uc_pixel[0]];//uc_pixel[2] = resultMap[uc_pixel[2]];//uc_pixel[1] = resultMap[uc_pixel[1]];//uc_pixel += 3;uc_pixel += SrcImage.channels();}}imwrite(outPutPath, resultImg);printf("Success save img:%s\n", outPutPath);if (IsDislay != 0){//是否显示影像imshow("SrcImage", SrcImage);//imshow("SrcImage_Equ", SrcImage_Equ);imshow("BaseMapImg", BaseMapImg);//imshow("BaseMapImg_Equ", BaseMapImg_Equ);imshow("resultImg", resultImg);waitKey(0);}return 0;
}void getSigBandImgHistogram(Mat SrcImage, map<int, double> &SrcImg_SigMap, int BandId)
{//获取指定波段的直方图信息int w = SrcImage.cols;int h = SrcImage.rows;double totalPixel = w * h;double rate = 1 / totalPixel;for (int i = 0; i < 256; i++) {//map.put(i, 0.0);// 通过循环，往集合里面填充0~255个位置,初始值都为0SrcImg_SigMap[i] = 0.0;}//分别统计图像上0~255上分布总数for (int row = 0; row < h; row++) {uchar* uc_pixel = SrcImage.data + row * SrcImage.step;for (int col = 0; col < w; col++) {SrcImg_SigMap[uc_pixel[BandId]] += rate;uc_pixel += 3;}}return;}

结果如下图