这两天一直在研究匀光匀色算法才了解到了直方图匹配算法，想要了解这个算法又要先了解直方图均衡化算法，通过网上查找了很多资料，没有现成C++代码，经过仔细思考和实验后大概复现了该算法。特此记录，以备查阅
参考链接如下：
1、匀光匀色–直方图匹配算法实现与应用
2、基于OpenCV的直方图匹配
3、直方图均衡化的数学原理
先介绍一下基本概念
原始数据
底图影像
直方图匹配匀光匀色效果
从上图可以看出原始影像与模板颜色不一致，经过直方图匹配匀光匀色后颜色一致性较好。
先看看直方图均衡算法原理知识：
先在根据上面的算法步骤，计算一幅影像的直方图均衡影像
1、获取一张影像Mat的直方图，即计算hs
void getImgHistogram(Mat image, map<int, double> & map)
{//获取一张影像Mat的直方图，即计算hs
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++) {
//map.put(i, 0.0);// 通过循环，往集合里面填充0~255个位置,初始值都为0
map[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;//计算累计直方图，即hp
for (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);//获取均衡化的影像Mat
int 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);//获取直方图均衡化的Mat
imwrite(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个位置,初始值都为0
map[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个位置,初始值都为0
SrcImg_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;
}
结果如下图
