程序:
/*
pImage:输入图像/输出图像
uiXRes:图像的宽,彩色图像为宽*3
uiYRes:图像高
Min:图像的像素的最小值
Max:图像像素的最大值
uiNrX:最小值大于2,最大值小于16,分割块的宽度
uiNrY:最小值大于2,最大值小于16,分割块的高度
uiNrBins:色阶数
fCliplimit:<=1.0不做剪裁,>1.0剪裁
*/
void CLAHE(uchar* pImage,uint uiXRes,uint uiYRes,uchar Min,
uchar Max,uint uiNrX,uint uiNrY,uint uiNrBins,float fCliplimit)
{
uint uiX, uiY;
uint uiXSize, uiYSize; //分割块的宽高的数量
uint uiSubX, uiSubY;
uint uiXL, uiXR, uiYU, uiYB;//线性插值的变量
ulong ulClipLimit, ulNrPixels; //剪裁区域和区域的像素数量
uchar* pImPointer;
uchar aLUT[256];
ulong* pulHist;//直方图
ulong *pulMapArray; //映射函数
ulong* pulLU, *pulLB, *pulRU, *pulRB;//线性插值的数据,分别表示左上、左下、右上、右下
pulMapArray=(ulong *)malloc(sizeof(ulong)*uiNrX*uiNrY*uiNrBins);
uiXSize = uiXRes/uiNrX;
uiYSize = uiYRes/uiNrY;
ulNrPixels = (ulong)uiXSize * (ulong)uiYSize;
if(fCliplimit > 0.0)
{
//计算实际剪裁
ulClipLimit = (ulong)(fCliplimit*(uiXSize*uiYSize)/uiNrBins);
ulClipLimit = (ulClipLimit<1UL)?1UL:ulClipLimit;
}
else
{
ulClipLimit = 1UL<<14;
}
MakeLut(aLUT,Min,Max,uiNrBins);//生成查找表
//每个块计算映射
for (uiY = 0, pImPointer = pImage; uiY < uiNrY; uiY++)
{
for (uiX = 0; uiX < uiNrX; uiX++, pImPointer += uiXSize)
{
pulHist = &pulMapArray[uiNrBins * (uiY * uiNrX + uiX)];
MakeHistogram(pImPointer,uiXRes,uiXSize,uiYSize,pulHist,uiNrBins,aLUT);
ClipHistogram(pulHist, uiNrBins, ulClipLimit);
MapHistogram(pulHist, Min, Max, uiNrBins, ulNrPixels);
}
pImPointer += (uiYSize - 1) * uiXRes; //skip lines, set pointer
}
//插值
for (pImPointer = pImage, uiY = 0; uiY < uiNrY; uiY++)
{
uiSubY = uiYSize;
if (uiY == uiNrY-1)
{
uiYU = uiNrY-1;
uiYB = uiNrY-1;
}
else
{
uiYU = uiY ;
uiYB = uiY + 1;
}
for (uiX = 0; uiX < uiNrX; uiX++)
{
uiSubX = uiXSize;
if (uiX == uiNrX-1)
{
uiXL = uiNrX-1;
uiXR = uiNrX-1;
}
else
{
uiXL = uiX;
uiXR = uiX + 1;
}
pulLU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXL)];
pulRU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXR)];
pulLB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXL)];
pulRB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXR)];
Interpolate(pImPointer,uiXRes,pulLU,pulRU,pulLB,pulRB,uiSubX,uiSubY,aLUT);
pImPointer += uiSubX; //set pointer on next matrix
}
pImPointer += (uiSubY - 1) * uiXRes;
}
free(pulMapArray); //free space for histograms
}
void MakeLut(uchar* pLUT,uchar Min,uchar Max,uint uiNrBins)
{
uchar BinSize = (uchar)(1+(Max-Min)/uiNrBins);
for (int i=Min;i<=Max;i++)
pLUT[i] = (i - Min) / BinSize;
}
void MakeHistogram (uchar* pImage,uint uiXRes,uint uiSizeX,uint uiSizeY,
ulong* pulHistogram,uint uiNrGreylevels,uchar* pLookupTable)
{
uchar* pImagePointer;
uint i;
for (i = 0; i < uiNrGreylevels; i++)
pulHistogram[i]=0L;
for (i = 0; i < uiSizeY; i++)
{
pImagePointer = &pImage[uiSizeX];
while (pImage < pImagePointer)
pulHistogram[pLookupTable[*pImage++]]++;
pImagePointer += uiXRes;
pImage = pImagePointer-uiSizeX;
}
}
void MapHistogram (ulong* pulHistogram,uchar Min,uchar Max,
uint uiNrGreylevels,ulong ulNrOfPixels)
{
uint i;
ulong ulSum = 0;
float fScale = ((float)(Max - Min)) / ulNrOfPixels;
ulong ulMin = (ulong) Min;
for (i = 0; i < uiNrGreylevels; i++)
{
ulSum += pulHistogram[i];
pulHistogram[i]=(ulong)(ulMin+ulSum*fScale);
if (pulHistogram[i] > Max)
pulHistogram[i] = Max;
}
// int I, Sum = 0, Amount = 0;
// for (I = 0; I < 256; I++) Amount += pulHistogram[I];
// for (I = 0; I < 256; I++)
// {
// Sum += pulHistogram[I];
// pulHistogram[I] = Sum * 255/ Amount; // 计算分布
// }
}
void ClipHistogram (ulong* pulHistogram,uint uiNrGreylevels,ulong ulClipLimit)
{
ulong* pulBinPointer, *pulEndPointer, *pulHisto;
ulong ulNrExcess, ulOldNrExcess, ulUpper, ulBinIncr, ulStepSize, i;
long lBinExcess;
ulNrExcess = 0; pulBinPointer = pulHistogram;
for (i = 0; i < uiNrGreylevels; i++)
{ /* calculate total number of excess pixels */
lBinExcess = (long) pulBinPointer[i] - (long) ulClipLimit;
if (lBinExcess > 0) ulNrExcess += lBinExcess; /* excess in current bin */
};
/* Second part: clip histogram and redistribute excess pixels in each bin */
ulBinIncr = ulNrExcess / uiNrGreylevels; /* average binincrement */
ulUpper = ulClipLimit - ulBinIncr; /* Bins larger than ulUpper set to cliplimit */
for (i = 0; i < uiNrGreylevels; i++)
{
if (pulHistogram[i] > ulClipLimit)
pulHistogram[i] = ulClipLimit; /* clip bin */
else
{
if (pulHistogram[i] > ulUpper)
{ /* high bin count */
//ulNrExcess -= (pulHistogram[i] - ulUpper); pulHistogram[i]=ulClipLimit;
ulNrExcess -= (ulClipLimit -pulHistogram[i]);
pulHistogram[i]=ulClipLimit;
}
else
{ /* low bin count */
ulNrExcess -= ulBinIncr;
pulHistogram[i] += ulBinIncr;
}
}
}
do { /* Redistribute remaining excess */
pulEndPointer = &pulHistogram[uiNrGreylevels]; pulHisto = pulHistogram;
ulOldNrExcess = ulNrExcess; /* Store number of excess pixels for test later. */
while (ulNrExcess && pulHisto < pulEndPointer)
{
ulStepSize = uiNrGreylevels / ulNrExcess;
if (ulStepSize < 1)
ulStepSize = 1; /* stepsize at least 1 */
for (pulBinPointer=pulHisto; pulBinPointer < pulEndPointer && ulNrExcess;
pulBinPointer += ulStepSize)
{
if (*pulBinPointer < ulClipLimit)
{
(*pulBinPointer)++; ulNrExcess--; /* reduce excess */
}
}
pulHisto++; /* restart redistributing on other bin location */
}
} while ((ulNrExcess) && (ulNrExcess < ulOldNrExcess));
}
void Interpolate(uchar* pImage,int uiXRes,ulong* pulMapLU,ulong* pulMapRU,
ulong* pulMapLB,ulong* pulMapRB,uint uiXSize,uint uiYSize,uchar* pLUT)
{
uint uiIncr = uiXRes-uiXSize; /* Pointer increment after processing row */
uchar GreyValue;
uint uiNum = uiXSize*uiYSize; /* Normalization factor */
uint uiXCoef, uiYCoef, uiXInvCoef, uiYInvCoef, uiShift = 0;
if (uiNum & (uiNum - 1)) /* If uiNum is not a power of two, use division */
for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;
uiYCoef++, uiYInvCoef--,pImage+=uiIncr)
{
for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;
uiXCoef++, uiXInvCoef--)
{
GreyValue = pLUT[*pImage]; /* get histogram bin value */
*pImage++ = (uchar ) ((uiYInvCoef * (uiXInvCoef*pulMapLU[GreyValue]
+ uiXCoef * pulMapRU[GreyValue])+ uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
+ uiXCoef * pulMapRB[GreyValue]))/ uiNum);
}
}
else
{ /* avoid the division and use a right shift instead */
while (uiNum >>= 1) uiShift++; /* Calculate 2log of uiNum */
for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;
uiYCoef++, uiYInvCoef--,pImage+=uiIncr)
{
for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;
uiXCoef++, uiXInvCoef--)
{
GreyValue = pLUT[*pImage]; /* get histogram bin value */
*pImage++ = (uchar)((uiYInvCoef*(uiXInvCoef * pulMapLU[GreyValue]
+ uiXCoef * pulMapRU[GreyValue])+uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
+ uiXCoef * pulMapRB[GreyValue])) >> uiShift);
}
}
}
//=====================================测试=========================================
#include "clahe.h"
#include "opencv/cv.h"
#include "opencv/cxcore.h"
#include "opencv/highgui.h"
#include "stdio.h"
#include <omp.h>
IplImage* Process1(IplImage *pImage)
{
IplImage *HIS = cvCreateImage( cvGetSize(pImage), 8, 3 );
IplImage *H = cvCreateImage( cvGetSize(pImage), 8, 1 );
IplImage *S = cvCreateImage( cvGetSize(pImage), 8, 1 );
IplImage *I = cvCreateImage( cvGetSize(pImage), 8, 1 );
cvCvtColor( pImage, HIS, CV_RGB2HLS);//改为CV_RGB2HLS可以转hls空间
cvCvtPixToPlane(HIS,H,I,S,0);
IplImage *resimage = cvCreateImage( cvGetSize(HIS), 8, 3 );
long Y_Height,X_Width,Step;
Y_Height=I->height;
X_Width=I->width;
Step=I->widthStep;
unsigned char *data;
data=(unsigned char*)I->imageData;
unsigned int NrX,NrY,uiNrBins;
NrX=2;
NrY=2;
uiNrBins=256;
double Min,Max;
cvMinMaxLoc(I, &Min, &Max);
IplImage *Image = cvCreateImage(cvGetSize(HIS), HIS->depth, HIS->nChannels);
IplImage *pImageChannel[4] = { 0, 0, 0, 0 };
IplImage *pImage_1;
if( I )
{
pImage_1 = cvCreateImage(cvGetSize(HIS), HIS->depth, 1);
cvCopy(I, pImage_1,NULL);
/*S = ContrastExtend(S);*/
CLAHE( (unsigned char *)(pImage_1->imageData),X_Width,Y_Height,0,255,NrX,NrY,uiNrBins,1.2);
CvScalar pixel,pixel_1;
for (int m = 0; m < I->height; m++)
{
for (int n = 0; n < I->width; n++)
{
pixel_1 = cvGet2D(pImage_1, m, n);
pixel.val[0] = pixel_1.val[0];
cvSet2D(I, m, n, pixel);
}
}
// 信道组合
cvMerge(H, I, S, 0, Image);
cvCvtColor(Image, resimage, CV_HLS2RGB);//改为CV_BGR2HLS可以转hls空间
}
// 释放资源
cvReleaseImage( &Image);
Image = 0;
return resimage;
}
void test1()
{
IplImage* test=cvLoadImage("4.jpg");
cvNamedWindow("test1",CV_WINDOW_AUTOSIZE);
cvShowImage("test1",test);
IplImage* res=Process1(test);
cvNamedWindow("res1",CV_WINDOW_AUTOSIZE);
cvShowImage("res1",res);
}
void test2()
{
IplImage* test=cvLoadImage("1.bmp");
int width=test->width;
int height=test->height;
int len=width*height;
unsigned char* dataR=(unsigned char*)calloc(width*height,sizeof(unsigned char));
unsigned char* dataG=(unsigned char*)calloc(width*height,sizeof(unsigned char));
unsigned char* dataB=(unsigned char*)calloc(width*height,sizeof(unsigned char));
for (int i=0;i<height;i++)
{
unsigned char* dataR_ptr=dataR+i*width;
unsigned char* dataG_ptr=dataG+i*width;
unsigned char* dataB_ptr=dataB+i*width;
char* test_ptr=test->imageData+i*test->widthStep;
for (int j=0;j<width;j++)
{
dataR_ptr[j]=test_ptr[j*3];
dataG_ptr[j]=test_ptr[j*3+1];
dataB_ptr[j]=test_ptr[j*3+2];
}
}
clock_t start, finish;
double duration;
start = clock();
#pragma omp parallel sections
{
#pragma omp section
CLAHE(dataR,width,height,0,255,2,8,256,2.0);
#pragma omp section
CLAHE(dataG,width,height,0,255,2,8,256,2.0);
#pragma omp section
CLAHE(dataB,width,height,0,255,2,8,256,2.0);
}
finish = clock();
duration = (double)(finish - start) / CLOCKS_PER_SEC;
printf("%lf\n",duration);
CvSize size = cvSize(width,height);
IplImage* res=cvCreateImage(size,IPL_DEPTH_8U,3);
for (int i=0;i<height;i++)
{
unsigned char* dataR_ptr=dataR+i*width;
unsigned char* dataG_ptr=dataG+i*width;
unsigned char* dataB_ptr=dataB+i*width;
char* res_ptr=res->imageData+i*res->widthStep;
for (int j=0;j<width;j++)
{
res_ptr[j*3]=dataR_ptr[j];
res_ptr[j*3+1]=dataG_ptr[j];
res_ptr[j*3+2]=dataB_ptr[j];
}
}
cvNamedWindow("res2",CV_WINDOW_AUTOSIZE);
cvShowImage("res2",res);
}
void test3()
{
IplImage* test=cvLoadImage("1.jpg");
int width=test->width;
int height=test->height;
cvNamedWindow("test3",CV_WINDOW_AUTOSIZE);
cvShowImage("test3",test);
unsigned char* data=(unsigned char*)calloc(width*height*3,sizeof(unsigned char));
for (int i=0;i<height;i++)
{
unsigned char* inputData_ptr=data+i*width*3;
char* test_ptr=test->imageData+i*test->widthStep;
for (int j=0;j<width;j++)
{
inputData_ptr[j*3]=test_ptr[j*3];
inputData_ptr[j*3+1]=test_ptr[j*3+1];
inputData_ptr[j*3+2]=test_ptr[j*3+2];
}
}
CLAHE(data,width*3,height,0,255,2,2,256,2.0);
CvSize size = cvSize(width,height);
IplImage* res=cvCreateImage(size,IPL_DEPTH_8U,3);
for (int i=0;i<height;i++)
{
unsigned char* outputData_ptr=data+i*width*3;
char* res_ptr=res->imageData+i*res->widthStep;
for (int j=0;j<width;j++)
{
res_ptr[j*3]=outputData_ptr[j*3];
res_ptr[j*3+1]=outputData_ptr[j*3+1];
res_ptr[j*3+2]=outputData_ptr[j*3+2];
}
}
cvNamedWindow("res3",CV_WINDOW_AUTOSIZE);
cvShowImage("res3",res);
}
int main()
{
test2();
cvWaitKey(0);
return 0;
}