图形处理中有一种对比度变换，像显示器就有对比度调节，PhotoShop也有图片的对比度修改，对比度的提高可以使图像细节清晰，相反，对比度的减小可以隐藏图像的细节，在一定程度上使图像柔和。
对比度变换其中一种比较简单的方法是直方图均衡化。
所谓直方图就是在某一灰度级的象素个数占整幅图像的象素比 h=nj/N,其中nj是灰度级在j的象素数，N是总象素数，扫描整幅图像得出的h的离散序列就是图像的直方图，h求和必然＝1，所以直方图可以看成是象素对于灰度的概率分布函数。
直方图是高低不齐的，因为象素灰度是随机变化的，直方图均衡化就是用一定的算法使直方图大致平和。
算法如下：
对于一个直方图
设 Pr(r)是原始图像直方图，Ps(s)是均衡化的直方图，
由于其是一个概率分布函数
所以有 Ps(s)ds=Pr(r)dr (编辑关系，ds,dr是积分变量)
因为要进行均衡化，令 Ps(s)=1,
得 ds=Pr(r)dr/1
两边积分得 s=F Pr(r)dr (因为编辑关系，左边F表示积分符号….-_-++)
数字图像是离散的，因此离散化上式得
sk=E{j=0,k}(nj/N) 左式k,j是离散量下标，因为编辑关系，E{0,k}表示下标0到k的连加符号，N是象素总数
由此得出每一象素的sk为均衡化后的正规化灰度（即灰度正规化到[0,1])，统计sk即可得出均衡化后的直方图。
在均衡化过程中可以对每一象素映射到新的实际灰度值sk*255，就实现了图像的变换
（严格理论中应该是灰度正规化到[0,1]区间，然后均衡化后的sk还要量化到原始的正规灰度以实现灰度合并，下面的BCB程序并没有量化，而且255是固定灰度级，因为256色BMP的彩色表就是256个表项)
现在开始实践
用BCB对一BMP灰度图像进行直方图均衡化处理，代码如下
//—————————-BCB6代码
#include <vcl.h>
#pragma hdrstop
#include<stdio.h>
#include "Unit1.h"
#include"File1.h"
#pragma pack(1)
//BMP文件头
struct BITMAPFILEHEADER
{
short type;
int bfSize;
short re1,re2;
int Offbits;
};
//BMP信息头
struct BITMAPINFO
{
long size;
long width,height;
short planes,bitCount;
long comp,sizeImg;
long xpels,ypels;
long used,important;
};
//BMP彩色表项
struct COLOR
{
char blue,green,red,re;
};
//——将BMP彩色表的数据校正到BCB TColor的数据。
void SwitchColor(long &c)
{
long blue=c& 0x000000ff;
long green=c& 0x0000ff00;
long red=c& 0x00ff0000;
c=(blue<<16) | green | (red>>16);
}
void xxx()
{
FILE f=fopen("f:\bbs_prev2.bmp","rb");
if(f==NULL) /判断文件是否打开成功/
{
ShowMessage("File open error");
return;
}
fseek(f,0,0);//移动到开头
//———-读BMP文件头
BITMAPFILEHEADER_ bmph=new BITMAPFILEHEADER();
if(fread((char*)bmph,sizeof(BITMAPFILEHEADER),1,f)==NULL)
{
ShowMessage("File read error");
return;
}
//———–读BMP信息头
BITMAPINFO *bmpi=new BITMAPINFO();
if(fread((char)bmpi,sizeof(BITMAPINFO_),1,f)==NULL)
{
ShowMessage("File read error2");
return;
}
//————–读彩色表
long c=new long[bmph->Offbits-sizeof(BITMAPFILEHEADER)-sizeof(BITMAPINFO)];
fread((char)c,bmph->Offbits-sizeof(BITMAPFILEHEADER)-sizeof(BITMAPINFO),1,f);
//———-显示一些信息
Form1->Edit1->Text=IntToStr(bmph->bfSize);
Form1->Edit2->Text=IntToStr(bmpi->width);
Form1->Edit3->Text=IntToStr(bmpi->height);
Form1->Edit4->Text=IntToStr(bmpi->comp);
Form1->Edit5->Text=IntToStr(bmpi->used);
int i,j,k,wc;
long N=bmph->bfSize- bmph->Offbits;//象素总数
unsigned char image=new char[N]; //位图矩阵
unsigned char newimage=new char[N];//变换后的位图矩阵
fread(image,N,1,f);//读入位图矩阵
//———直方图数列初始化
//———直方图数列用来存储正规化后的灰度
double h=new double[255];//255个灰度级，保存原始图像正规化灰度直方图数据
for(i=0;i<255;i++)
h[i]=0.0;
double nh=new double[255];//255个灰度级，保存变换后的图像正规化灰度直方图
for(i=0;i<255;i++)
nh[i]=0.0;
long count=new long[255]; //每一灰度级的象素数量统计
for(i=0;i<255;i++)
count[i]=0;
for(i=0;i<N;i++)
{
count[image[i]]++;
}
//－－－－－正规化灰度概率统计
for(i=0;i<255;i++)
{
h[i]=count[i]/(double)N;
}
//——正规化新灰度图
double hc;
for(i=0;i<N;i++)
{
hc=0;
for(j=0;j<image[i];j++)
hc+=h[j];
nh[image[i]]+=hc; //保存新正规化灰度图
newimage[i]=hc255; //保存新图像灰度索引
}
//———-显示直方图
for(i=0;i<255;i++)
{
//原始直方图
Form1->Canvas->MoveTo(10+i,200);
Form1->Canvas->LineTo(10+i,200+h[i]N);
//新直方图
Form1->Canvas->MoveTo(300+i,200);
Form1->Canvas->LineTo(300+i,200+nh[i]255);
}
//——显示图形
TColor tc;
if(bmpi->width%4==0)//———–因为BMP图像4字节对齐
wc=bmpi->width/44;
else
wc=(bmpi->width/4+1)4;
long pos=0;
for( i=0;i<bmpi->height;i++)
{
for(j=0;j<wc;j++)
{
//—–原始图形
a= c[image[pos]];
SwitchColor(a);
Form1->Canvas->Pixels[10+j][600-i]=a;
//——新图形
a= c[newimage[pos]];
SwitchColor(a);
Form1->Canvas->Pixels[300+j][600-i]=a;
pos++;
}
}
fclose(f);
}
这个程序使用256色BMP文件，但程序代码是针对灰度图像的，用于彩色图像时得出一些古怪色彩配合而已。
在对灰度图像均衡化时
如果原始图像对比度本来就很高，如果再均衡化则灰度调和，对比度降低。
在泛白缓和的图像中，由于均衡化过程中会合并一些象素灰度，则会增大对比度，这里255灰度级太多，合并不明显。
5.9 灰度直方图均衡化http://www-scf.usc.edu/~flv/ipbook/chap05.htm
在介绍灰度直方图均衡化(histogram equalization)之前，先讲讲直方图修正。所谓直方图修正，就是通过一个灰度映射函数Gnew=F(Gold)，将原灰度直方图改造成你所希望的直方图。所以，直方图修正的关键就是灰度映射函数。我们刚才介绍的阈值化、削波、灰度窗口变换等等，都是灰度映射函数。
直方图均衡化是一种最常用的直方图修正。它是把给定图象的直方图分布改造成均匀直方图分布。由信息学的理论来解释，具有最大熵(信息量)的图象为均衡化图象。直观地讲，直方图均衡化导致图象的对比度增加。
由于直方图均衡化涉及到很多概率和数学的知识，具体的细节这里就不介绍了，只给出算法。通过下面的例子，就很容易明白了。
有一幅图象，共有16级灰度，其直方图分布为Pi, i=0,1,…,15，求经直方图均衡化后，量化级别为10级的灰度图象的直方图分布Qi，其中Pi和Qi为分布的概率，即灰度i出现的次数与总的点数之比。
Pi： 0.03,0,0.06,0.10,0.20,0.11,0,0,0,0.03,0,0.06,0.10,0.20,0.11,0
步骤1：用一个数组s记录Pi，即
s[0]=0.03,s[1]=0,s[2]=0.06,…,s[14]=0.11,s[15]=0
步骤2：i从1开始，令s[i]=s[i]+s[i-1]+ … + s[0]，得到的结果是
s： 0.03,0.03,0.09,0.19,0.39,0.50,0.50,0.50,0.50,0.53,0.53,0.59,0.69,0.89,1.0,1.0
步骤3：用一个数组L记录新的调色板索引值，即令L[i]=s[i]×(10-1)，得到的结果是L：0,0,1,2,4,5,5,5,5,5,5,5,6,8,9,9
这样就找到了原来的调色板索引值和新的调色板索引值之间的对应关系，即
0→0,1→0,2→1,3→2,4→4,5→5,6→5,7→5,8→5,9→5,10→5,11→5,12→6,
13→8,14→9,15→9。
步骤4：将老的索引值对应的概率合并，作为对应的新的索引值的概率。例如，原来的索引值0,1都对应了新的索引值0，则灰度索引值为0的概率为P0+P1=0.03；新的索引值3和7找不到老的索引值与之对应，所以令Q3和Q7为0。最后得到的结果是Qi：0.03,0.06,0.10,0,0.20,0.20,0.10,0,0.20,0.11
图5.17为Pi的分布，图5.18为Qi的分布，对照一下，不难发现图5.18的分布比图5.17要均匀一些。
图5.17 Pi的分布
图5.18 Qi的分布
要注意的是，均衡化处理后的图象只能是近似均匀分布。均衡化图象的动态范围扩大了，但其本质是扩大了量化间隔，而量化级别反而减少了，因此，原来灰度不同的象素经处理后可能变的相同，形成了一片的相同灰度的区域，各区域之间有明显的边界，从而出现了伪轮廓。
图5.19为图5.13经直方图均衡化处理后，量化为128级灰度的结果；图5.20为它的直方图分布。为什么天亮了起来呢？分析一下就明白了：因为原图中低灰度的点太多了，所以s数组前面的元素很大。经过L[i]=s[i]×(128-1)的处理后，原图中低灰度的点的灰度值提高了不少，所以那片暗区变亮了。同时可以看出，天空中出现了伪轮廓。
图5.20 图5.19的灰度直方图
图5.21为图5.1直方图均衡化后的结果(128级灰度)，暗的区域(手)变亮了，看起来更清楚一些。
图5.21 图5.1直方图均衡化后的结果
下面给出直方图均衡化的源程序：
int EquaScale; //为新的灰度级别
BOOL HistogramEqua(HWND hWnd)
{
DLGPROC dlgInputBox = NULL;
DWORD BufSize,OffBits;
LPBITMAPINFOHEADER lpImgData;
LPSTR lpPtr;
HLOCAL hTempImgData;
LPBITMAPINFOHEADER lpTempImgData;
LPSTR lpTempPtr;
HDC hDc;
HFILE hf;
LONG x,y;
LOGPALETTE pPal;
HPALETTE hPrevPalette;
HLOCAL hPal;
WORD i;
int Gray;
DWORD GrayHits[256];
int GrayIndex[256];
float s[256];
if( NumColors!=256){ //必须是256级灰度图
MessageBox(hWnd,"Must be a 256 grayscale bitmap!","Error Message",
MB_OK|MB_ICONEXCLAMATION);
return FALSE;
}
//出现对话框，输入新的灰度级别
dlgInputBox = (DLGPROC) MakeProcInstance ( (FARPROC)InputBox, ghInst );
DialogBox (ghInst, "INPUTBOX", hWnd, dlgInputBox);
FreeProcInstance ( (FARPROC) dlgInputBox );
if( EquaScale >=255){ //量化级别不能大于255
MessageBox(hWnd,"The new scale can not be larger than 255",
"Error Message",MB_OK|MB_ICONEXCLAMATION);
return FALSE;
}
//OffBits为到实际位图数据的偏移值
OffBits=bf.bfOffBits-sizeof(BITMAPFILEHEADER);
//BufSize为缓冲区的大小
BufSize=OffBits+bi.biHeightLineBytes; if((hTempImgData=LocalAlloc(LHND,BufSize))==NULL)
{
MessageBox(hWnd,"Error alloc memory!","Error Message",MB_OK|
MB_ICONEXCLAMATION);
return FALSE;
}
//lpImgData指向原图，//lpTempImgData指向新开的缓冲区
lpImgData=(LPBITMAPINFOHEADER)GlobalLock(hImgData);
lpTempImgData=(LPBITMAPINFOHEADER)LocalLock(hTempImgData);
//拷贝头信息
memcpy(lpTempImgData,lpImgData,OffBits);
//ColorHits为记录颜色使用频率的数组，ColorIndex为记录颜色索引值的
//数组
//先清零
memset(GrayHits,0,256sizeof(DWORD));
memset(GrayIndex,0,256sizeof(WORD));
for(y=0;y<bi.biHeight;y++){
lpPtr=(unsigned char )lpImgData+(BufSize-LineBytes-yLineBytes);
for(x=0;x<bi.biWidth;x++){
Gray=(unsigned char )(lpPtr++);
GrayHits[Gray]++; //统计该颜色用到的次数
}
}
for(i=0;i<256;i++)
//次数除以总点数得到频率
s[i]=(float)GrayHits[i]/((float)bi.biWidth(float)bi.biHeight);
for(i=1;i<256;i++)
s[i]+=s[i-1]; //每一项都是前面所有项的累加
for(i=0;i<256;i++)
//根据新的量化级别，计算新的灰度索引值
GrayIndex[i]=(int)(s[i](EquaScale-1));
//为新的调色板分配内存
hPal=LocalAlloc(LHND,sizeof(LOGPALETTE)+
256sizeof(PALETTEENTRY));
pPal =(LOGPALETTE )LocalLock(hPal);
//先将调色板内存全部清零
memset(pPal,0,sizeof(LOGPALETTE) + 256 sizeof(PALETTEENTRY));
pPal->palNumEntries =(WORD) 256;
pPal->palVersion = 0x300;
lpTempPtr=(char )lpTempImgData+sizeof(BITMAPINFOHEADER);
for (i = 0; i < EquaScale; i++) {
Gray=(int)(i255.0/(EquaScale-1)); //根据新的量化级别，计算灰度值
pPal->palPalEntry[i].peRed=(BYTE)Gray;
pPal->palPalEntry[i].peGreen=(BYTE)Gray;
pPal->palPalEntry[i].peBlue=(BYTE)Gray;
pPal->palPalEntry[i].peFlags=(BYTE)0;
(lpTempPtr++)=(unsigned char)Gray;
(lpTempPtr++)=(unsigned char)Gray;
(lpTempPtr++)=(unsigned char)Gray;
(lpTempPtr++)=0;
}
if(hPalette!=NULL)
DeleteObject(hPalette);
//产生新的逻辑调色板
hPalette=CreatePalette(pPal);
LocalUnlock(hPal);
LocalFree(hPal);
hDc=GetDC(hWnd);
if(hPalette){
hPrevPalette=SelectPalette(hDc,hPalette,FALSE);
RealizePalette(hDc);
}
for(y=0;y<bi.biHeight;y++){
lpPtr=(char )lpImgData+(BufSize-LineBytes-yLineBytes);
lpTempPtr=(char )lpTempImgData+(BufSize-LineBytes-yLineBytes);
for(x=0;x<bi.biWidth;x++){
Gray=(unsigned char )(lpPtr++); //原灰度索引值
Gray=GrayIndex[Gray]; //对应的新的灰度索引值
(lpTempPtr++)=(unsigned char)Gray;
}
}
if(hBitmap!=NULL)
DeleteObject(hBitmap);
//产生新的位图
hBitmap=CreateDIBitmap(hDc,(LPBITMAPINFOHEADER)lpTempImgData,
(LONG)CBM_INIT,
(LPSTR)lpTempImgData+
sizeof(BITMAPINFOHEADER)+
256*sizeof(RGBQUAD),
(LPBITMAPINFO)lpTempImgData,
DIB_RGB_COLORS);
if(hPalette && hPrevPalette){
SelectPalette(hDc,hPrevPalette,FALSE);
RealizePalette(hDc);
}
hf=_lcreat("c:\equa.bmp",0);
_lwrite(hf,(LPSTR)&bf,sizeof(BITMAPFILEHEADER));
_lwrite(hf,(LPSTR)lpTempImgData,BufSize);
_lclose(hf);
//释放内存和资源
ReleaseDC(hWnd,hDc);
LocalUnlock(hTempImgData);
LocalFree(hTempImgData);
GlobalUnlock(hImgData);
return TRUE;
}