关于直方图规范化的C++编程实现_2015_7_24

直方图规范化原理简介

07.24.0代码版本

1.直方图的规范化

直方图的规范化主要完成的任务是提高图像的对比度，实现过程主要如下，将灰度图像的的直方图计算出来，然后通过将不均匀的直方图进行拉伸（即对原图像中的灰度图进行调整）使得整个直方图灰度分布较为平均，这样得到的图像的对比度较高。是预处理阶段较为常见的而有效的图像处理手段。

最简单的直方图是线性变换通过对密集区域的拉伸实现灰度的平均分布。但有时难以达到理想的效果，所以一种改进的的方法可以对原图像的灰度直方图进行分段，分为拉伸端和压缩段已达到更好的效果。另一种改进的方法是先设计一种灰度拉伸的结果然建立原图像的灰度直方图和拉升后的灰度直方图的映射，可以不同于原图像分布的可以订制的直方图分布。

本文参考一篇论文的形式将进行原直方图进行分段拉伸，主要的难点在于选取原直方图的分段阈值点。这里采用最小误差法的方式。

最小误差法将原图像的灰度直方图分成对象和背景两个部分，并且将两个部分认为是两个高斯分布的拟合。用一个评价函数来确定阈值的值，当评价函数最小的时候说明设定的阈值越好，拟合效果越好，本文通过两次最小误差法的计算，第一次阈值用于将图像分成对象和非两个部分，第二次用阈值非对象分成背景和过度部分。然后对对象部分进行拉伸，对过度部分保持不变，对背景进行压缩。用将原图像分成三段并进行不同处理的方式完成提高对比度的任务。

PS:本人度这套算法有两点不同意见。

首先，我认为第二次阈值计算的方式有些不太合理，因为虽然可以将对象和背景视为两个高斯分布的拟合，但是将过度部分和背景也视为两个高斯部分的拟合有些不太合理，而且实际的效果也说明只能针对特定的图像完成两侧阈值的设置，有些图像在进行第二次阈值判断时会出现将阈值设置在整个灰度范围的两端，也就是说评价函数认为这个带判断的区域为可以视为一个高斯函数的拟合。

其次，对于第二次阈值判断时，是建立在第一次阈值的基础上的，但是如果是对阈值之前的灰度进行第二次阈值判断还是对阈值之后的灰度直方图进行第二次阈值判断，文章中采用的是根据图像的直方图直观感觉，这样的方式难以实现自动化的阈值判断实用性有待提高。

2.07.24.0版本

07.24.0目前完成编程实现，但是依然还是C语言的风格，接下来的任务把程序改写成类的形式。目前代码如下：

// thread1.cpp : 定义控制台应用程序的入口点。
//

#include <stdio.h>
#include <tchar.h>
#include <opencv2\opencv.hpp>
#include <highgui.h>
#include <cv.h>
#include <math.h>

using namespace std;
using namespace cv;

//定义一个直方图的类
class Histogram1D{
private:
int histSize[1];    //项的数目
float hranges[2];   //像素的最小值和最大值
const float * ranges[1];
int channels[1];    //仅用到一个通道
public:
Histogram1D(){
//准备1D直方图的参数
histSize[0] = 256;
hranges[0] = 0.0;
hranges[1] = 255.0;
ranges[0] = hranges;
channels[0] = 0;    //默认情况下我们考察0通道。
}

//计算一个图像的直方图，并返回一个列表
MatND getHistogram(const Mat &image){
MatND hist;
calcHist(&image,
1,          //计算单张图片的直方图
channels,   //通道数量
Mat(),      //不使用图像作为掩码
hist,       //返回直方图
1,          //这是1D的直方图
histSize,   //项的数量
ranges      //像素值的范围
);
return hist;
}

//计算1D直方图，并返回一幅图像
Mat getHistogramImage(const Mat &image, int minValue = 0){
//首先计算直方图
MatND hist = getHistogram(image);
//获取最大值和最小值
double maxVal = 0;
double minVal = 0;
minMaxLoc(hist, &minVal, &maxVal, 0, 0);
//显示直方图的图像
Mat histImg(histSize[0], histSize[0], CV_8U, Scalar(255));
//设置最高点为nbing的90%
int hpt = static_cast<int>(0.9*histSize[0]);
//每个条目都绘制一条垂直线
for (int h = 0; h < histSize[0]; h++){
float binVal = hist.at<float>(h);
int intensity = static_cast<int>(binVal*hpt / maxVal);
//两点之间绘制一条直线
line(histImg, Point(h, histSize[0]),
Point(h, histSize[0] - intensity),
Scalar::all(0));
}
return histImg;
}
};

//matlab联调接口函数

//计算各灰度值出现的概率
void ComProArr(float probabilityGray[], //输出的各灰度值概率的数组
const MatND  histout_arr,           //直方图各灰度值出现概率
const int number,                   //计算灰度值范围
const int rows,                     //行数
const int cols,                     //列数
const int channels                  //通道数
){
float TotalPointNumber = rows*cols*channels;
float testsum = 0;
for (int i = 0; i < 256; i++)
{
probabilityGray[i] = histout_arr.at<float>(i) / TotalPointNumber;
//                      cout <<"proability:  "<<i<<"    "<< probabilityGray[i] << endl;

}
}

//根据输入阈值来计算两个前景和背景的概率
void ComProSum(float  probabilityphase[],//输出概率存放数组地址
const int thread,                   //阈值相对位置
const float  probabilityGray[], //各灰度级的概率
const int number                    //灰度级总数
){
if (thread > number)
{
cout << "阈值设置错误";
}
float total = 0.0;      //用于计算实际灰度分布区域用灰度分割的概率的
float sumPro = 0.0;     //用于计算实际灰度分布区域的总概率
for (int i = 0; i < thread; i++)
total += probabilityGray[i];
probabilityphase[0] = total;
sumPro = total;
total = 0;
for (int i = thread; i < number; i++)
total += probabilityGray[i];
sumPro += total;
probabilityphase[0] = probabilityphase[0] / sumPro;
probabilityphase[1] = total / sumPro;
//      cout << "proabilityphase:0   " << thread<<"    "<< probabilityphase[0] << endl;
//      cout << "proabilityphase:1   " << thread << "    " << probabilityphase[1] << endl;
//      cout <<"sumproabilityphase:   "<< probabilityphase[0] + probabilityphase[1]<<endl;
}

//根据阈值输入来计算两个均值
void ComMean(float meanPhase[],     //输出计算均值
const float  probabilityGray[],//直方图各灰度出现概率
const int number,               //计算灰度的范围
const int thread,               //阈值相对位置
const float probailityphase[]   //输入阶段概率地址
){
float sumpro = 0.0;
for (int i = 0; i < number; i++)
{
sumpro += probabilityGray[i];
}
//      cout << "sumpro:" << sumpro;
//      cout <<"probabilityPhase0:"<< probailityphase[0]<<endl;
//      cout << "probabilityPhase1:" << probailityphase[1]<<endl;
float temp = 0;
for (int i = 0; i < thread; i++){
temp += i*probabilityGray[i] / sumpro;
}
temp = temp / (probailityphase[0] + 0.000000001);
meanPhase[0] = temp;
//      cout<<"meanPhase[0]  "<<thread<<"   "<< temp << endl;

temp = 0.0;
for (int i = thread; i < number; i++){
temp += i *probabilityGray[i] / sumpro;
//      cout << "probability:" << i << "    " << probabilityGray[i] << endl;
}
temp = temp / probailityphase[1];
meanPhase[1] = temp;
//      cout <<thread << "meanPhase[1]" << temp << endl;
}

void ComVariance(float variancePhase[],     //输出计算方差
const float  probabilityGray[],//直方图各灰度出现概率
const int number,               //计算灰度的范围
const int absothread,           //绝对阈值
const int thread,               //相对阈值
const float probailityphase[],  //阶段概率存放地址
const float meanPhase[] //阶段概率存放地址
){
float sumpro = 0.0;
for (int i = 0; i < number; i++)
{
sumpro += probabilityGray[i];
}
float temp = 0.0;
//  for (int i = 0; i < thread; i++){
//      temp += (((i + absothread) - meanPhase[0]) *((i + absothread) - meanPhase[0])) / sumpro* probabilityGray[i];
//  }
//  temp = temp / (probailityphase[0] + 0.000000001);
//  variancePhase[0] = temp;
for (int i = 0; i < thread; i++)
{
temp += probabilityGray[i] / sumpro  * (i + absothread)* (i + absothread);
}
temp = temp / probailityphase[0] - meanPhase[0] * meanPhase[0];
variancePhase[0] = temp;
//      cout <<thread <<" "<<"variancePhase0""  "<< variancePhase[0]<<endl;

temp = 0.0;
//  for (int i = thread; i < number; i++){
//      temp += (((i + absothread) - meanPhase[1]) *((i + absothread) - meanPhase[1])) / sumpro* probabilityGray[i];
//  }
//  temp = temp / probailityphase[1];
//  variancePhase[1] = temp;

for (int i = thread; i < number; i++)
{
temp += probabilityGray[i] / sumpro*(i + absothread)*(i + absothread);
}
temp = temp / probailityphase[1] - meanPhase[1] * meanPhase[1];
variancePhase[1] = temp;
//      cout <<thread<<"  "<<"variancePhase1   "<< variancePhase[1]<<endl;
}

void ComEvaluate(double & evaluationPhase,  //输出计算评价函数的值
const float probabilityphase[],             //输入各阶段的概率
const float variancePhase[]             //输出各阶段的方差

){
//  cout << endl;
//  cout << endl;
//  cout << endl;
//  cout << endl;
//  cout << probabilityphase[0] << endl;
//  cout << probabilityphase[1] << endl;
//  cout << variancePhase[0] << endl;
//  cout << variancePhase[1] << endl;
evaluationPhase = 1 + 2 * ((double)probabilityphase[0] * log(sqrt((double)variancePhase[0])) + (double)probabilityphase[1] * log(sqrt((double)variancePhase[1])))
- 2 * ((double)probabilityphase[0] * log((double)probabilityphase[0]) + (double)probabilityphase[1] * log((double)probabilityphase[1]));
//      evaluationPhase = probabilityphase[0] * log(variancePhase[0] / probabilityphase[0] / probabilityphase[0]) +
//          probabilityphase[1] * log(variancePhase[1] / probabilityphase[1] / probabilityphase[1]);
//          cout <<"  evaluation  "<<"    "<< evaluationPhase << endl;
};

void findMinScript(int & minSubscript,      //输出最小坐标
const double evaluatelist[],                            //输入数组首地址
const int number                                //输入数组数量
){
double temp = 500000;
int subscript = 0;
//  cout << "find start:" << endl;
//  cout << temp;

for (int i = 0; i < number; i++)
{
if (temp > abs(evaluatelist[i]))
{
temp = evaluatelist[i];
minSubscript = i;
}
//      cout << "J list:" << i << "    " << evaluatelist[i]<<endl;
//          cout << evaluatelist[i]<<endl;
//          cout << temp << endl;
}
//  cout << "find end:" << endl;
}

void Fcomspan(int & firstvalue,     //灰度值开始点
int &lastvalue,                 //灰度值结束点
const MatND histout_arr)        //histout_arr
{
for (int i = 0; i < 256; i++)
if (histout_arr.at<float>(i)>9)
{
firstvalue = i;
break;
}
for (int i = 255; i > 0; i--)
if (histout_arr.at<float>(i) > 9)
{
lastvalue = i;
break;
}
//cout << endl << firstvalue;
//cout << endl << lastvalue;
}

void scanImage(Mat &out_img,                    //输出拉伸后图像
const int thread1,              //输入阈值1
const int thread2,              //输入阈值2
const float k0,                 //输入斜率0
const float k2,                 //输入斜率2
const int firstvalue,           //输入实际范围最小值
const int lastvalue,            //输出实际范围最大值
const int c,                    //调用c
const int d                     //调用d
){
/*  Mat_<Vec3b>::iterator it = out_img.begin<Vec3b>();
Mat_<Vec3b>::iterator itend = out_img.end<Vec3b>();
for (; it != itend; ++it)
{
(*it)[0] = 256-(*it)[0];

if ((*it)[0] < thread1)
(*it)[0] = k0*((*it)[0] - firstvalue);
else if ((*it)[0] < thread2)
(*it)[0] = (*it)[0] - thread1 + c;
else
(*it)[0] = k2*((*it)[0] - thread2) + d;
*/
int temp = 0;
for (int j = 0; j < out_img.rows; j++)
{
uchar *data = out_img.ptr<uchar>(j);
for (int i = 0; i < out_img.cols; i++)
{

if (data[i] < thread1)
temp = k0*(data[i] - firstvalue);
else if (data[i] < thread2)
temp = data[i] - thread1 + c;
else
temp = k2*(data[i] - thread2) + d;
if (temp>255)
temp = 254;
data[i] = temp;
}
}
}

int _tmain(int argc, _TCHAR* argv[])
{
Mat src_img;                //原始图像
src_img = imread("src.jpg");

if (src_img.empty())
{
cout << "打开图像失败！" << endl;
return -1;
};

//调用窗口显示
namedWindow("source image", CV_WINDOW_AUTOSIZE);

imshow("source image", src_img);

//转换成灰度图
Histogram1D hist_arr;
MatND histout_arr = hist_arr.getHistogram(src_img);
//可视化显示
namedWindow("Histogram");
imshow("Histogram", hist_arr.getHistogramImage(src_img));

//计算各点出现的概率
float probabilityGray[256];
ComProArr(probabilityGray, histout_arr, 256, src_img.rows, src_img.cols, 1);

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
int firstvalue = 0;
int lastvalue = 0;

Fcomspan(firstvalue, lastvalue, histout_arr);

//  for (lastvalue = 1; lastvalue < 70; lastvalue++)

//      cout << " lastvalue = " << lastvalue<<"---";

//  lastvalue = 120;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//设定阈值

cout << "thread compute start:" << endl;
int thread = 0;
const int size = 256;
double evaluatelist[256];
int minSubscript[2];
minSubscript[1] = lastvalue;
minSubscript[0] = firstvalue;
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

for (; thread < lastvalue; thread++)
{
//根据设定的值计算P0和P1的总概率
float probabilityphase[2];
//      cout << "sum probability compute start" << endl;
ComProSum(probabilityphase, thread - minSubscript[0], &probabilityGray[minSubscript[0]], minSubscript[1] - minSubscript[0]);
/*
float total = 0;
for (int i = 68; i < 78; i++)
{
total += probabilityGray[i];
cout << probabilityGray[i] << endl;
}
cout << total;
*/
//根据设定的阈值计算U0和U1的均值
float meanPhase[2] = { 0.0, 0.0 };
ComMean(meanPhase, &probabilityGray[minSubscript[0]], minSubscript[1] - minSubscript[0], thread - minSubscript[0], probabilityphase);
meanPhase[0] = minSubscript[0] + meanPhase[0];
meanPhase[1] = minSubscript[0] + meanPhase[1];

//根据设定的阈值计算sigma0和sigma的方差
float variancePhase[2] = { 0.0, 0.0 };
ComVariance(variancePhase, &probabilityGray[minSubscript[0]], minSubscript[1] - minSubscript[0], minSubscript[0], thread - minSubscript[0], probabilityphase, meanPhase);

//根据设定的阈值计算评价函数
ComEvaluate(evaluatelist[thread], probabilityphase, variancePhase);
}
int thread1 = 0;
findMinScript(thread1, &evaluatelist[minSubscript[0]], minSubscript[1] - minSubscript[0]);
thread1 = thread1 + minSubscript[0];
cout << thread1 << endl;

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
minSubscript[0] = thread1;
thread = thread1;
for (; thread < lastvalue; thread++)
{
//根据设定的值计算P0和P1的总概率
float probabilityphase[2];
//      cout << "sum probability compute start" << endl;
ComProSum(probabilityphase, thread - minSubscript[0], &probabilityGray[minSubscript[0]], minSubscript[1] - minSubscript[0]);

//      float total = 0;
//      for (int i = 68; i < 78; i++)
//      {
//          total += probabilityGray[i];
//          cout << probabilityGray[i] << endl;
//      }
//      cout << total;

//根据设定的阈值计算U0和U1的均值
float meanPhase[2] = { 0.0, 0.0 };
ComMean(meanPhase, &probabilityGray[minSubscript[0]], minSubscript[1] - minSubscript[0], thread - minSubscript[0], probabilityphase);
meanPhase[0] = minSubscript[0] + meanPhase[0];
meanPhase[1] = minSubscript[0] + meanPhase[1];

//根据设定的阈值计算sigma0和sigma的方差
float variancePhase[2] = { 0.0, 0.0 };
ComVariance(variancePhase, &probabilityGray[minSubscript[0]], minSubscript[1] - minSubscript[0], minSubscript[0], thread - minSubscript[0], probabilityphase, meanPhase);

//根据设定的阈值计算评价函数
ComEvaluate(evaluatelist[thread], probabilityphase, variancePhase);
}
int thread2 = 0;
findMinScript(thread2, &evaluatelist[minSubscript[0]], minSubscript[1] - minSubscript[0]);
thread2 = thread2 + minSubscript[0];
cout << thread2 << endl;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//下面计算cd的范围
float k0 = 0;
float k1 = 1;
float k2 = 0;
float c = 0;
float d = 0;

k0 = (float)(256 - lastvalue + thread1) / (float)(thread1 - firstvalue);
k2 = (256 - thread2 + thread1 - k0*(thread1 - firstvalue)) / (thread1 - firstvalue);
c = k0*(thread1 - 0);
d = thread2 - thread1 + k0*(thread1 - 0);

cout << endl << "parameter display:" << endl;
cout << "minvalue:" << firstvalue << endl;
cout << "maxvalue:" << lastvalue << endl;
cout << "thread0:" << thread1 << endl;
cout << "thread1:" << thread2 << endl;
cout << "k0:" << k0 << endl;
cout << "k1:" << k1 << endl;
cout << "k2:" << k2 << endl;
cout << "c:" << c << endl;
cout << "d:" << d << endl;

Mat out_img;
out_img.create(src_img.size(), src_img.type());
cvtColor(src_img, out_img, CV_BGR2GRAY);

scanImage(out_img, thread1, thread2, k0, k2, firstvalue, lastvalue, c, d);
namedWindow("gray image", CV_WINDOW_AUTOSIZE);
imshow("gray image", out_img);

Histogram1D hist_arrout1d;
MatND hist_arrout = hist_arrout1d.getHistogram(out_img);
//可视化显示
namedWindow("HistogramOUT");
imshow("HistogramOUT", hist_arrout1d.getHistogramImage(out_img));

waitKey();
return 0;
}