在介绍灰度直方图均衡化(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.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.19 图5.13经直方图均衡化处理后的结果
图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.biHeight*LineBytes; 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,256*sizeof(DWORD));
memset(GrayIndex,0,256*sizeof(WORD));
for(y=0;y<bi.biHeight;y++){
lpPtr=(unsigned char *)lpImgData+(BufSize-LineBytes-y*LineBytes);
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)+
256*sizeof(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)(i*255.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-y*LineBytes);
lpTempPtr=(char *)lpTempImgData+(BufSize-LineBytes-y*LineBytes);
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;
}