精华内容
参与话题
问答
  • CLAHE算法

    万次阅读 2019-04-08 18:28:52
    限制对比度自适应直方图均衡(Contrast Limited Adaptive Histogram Equalization,CLAHE) 自适应直方图均衡(AHE)算法,对于图像中存在明显比其他区域亮或者暗的地方时,普通的直方图均衡算法就不能将该处的细节信息...

    限制对比度自适应直方图均衡(Contrast Limited Adaptive Histogram Equalization,CLAHE)

    自适应直方图均衡(AHE)算法,对于图像中存在明显比其他区域亮或者暗的地方时,普通的直方图均衡算法就不能将该处的细节信息描述出来。AHE算法通过在当前处理像素周边的一个矩形区域内进行直方图均衡,来达到扩大局部对比度,显示平滑区域细节的作用。
       AHE算法的2个属性:1、AHE算法处理的局部领域,矩形领域小,局部对比度强,矩形领域大,局部对比度弱。2、如果矩形区域内的图像块信息比较平坦,灰度接近,其灰度直方图呈尖状,在直方图均衡的过程中就可能会出现过度放大噪声的情况。

       CLAHE,对比度受限的自适应直方图均衡算法就能够有效的限制噪声放大的情形。下图表示的就是局部矩形领域内的灰度直方图,由于对比度放大的程度与像素点的概率分布直方图的曲线斜度成比例,所以为了限制对比度,将大于一定阈值的部分平均分配到直方图的其他地方,如下图所示,这样的话,通过限制CDF(累积分布函数)的斜率来一定程度限制对比度。
     

    插值过程,得到了CDF函数,也就获得了对应的亮度变换函数,在计算变换函数的时候可以通过插值过程来降低计算量。其中红色块(图像角点处)的变换函数是完全按照定义获得的,绿色块(图像边缘)的变换函数是通过旁边两个图像块的变换函数线性插值得到的,蓝色部分图像块的变换函数则是通过双线性插值得到。
     

    Matlab和OpenCV中都已经集成了CLAHE函数,在Matlab中,就是函数J = adapthisteq(I);

    全局映射算子
    每一个像素点将会根据它的全图特征和亮度信息进行映射,不管其空间位置几何。全局算子一个比较典型的例子就是色调曲线。
    导致图像色调映射过后看起来很平坦,失去了其局部的细节信息。

    局部映射算子
    像素点所在的空间位置会被考虑,具有相同亮度值的两个像素点会被映射成不同的值,因为它们的空间位置周边的亮度信息可能不一样。
    局部色调映射会很好的保护高亮和阴影部分的局部对比度和细节信息。
    https://www.cnblogs.com/whw19818/p/5766005.html

     算法主要分以下三部:
     1.图像分块,以块为单位,先计算直方图,然后修剪直方图,最后均衡;
     

     2.块间线性插值,这里需要遍历、操作各个图像块,处理起来复杂一些;
     

    3.与原图做图层滤色混合操作:f(a, b) = 1 - (1 - a)*(1 - b)。
    https://blog.csdn.net/grafx/article/details/53311915
     

    展开全文
  • 1. 概述 本篇文章是基于这篇博文写的,然后经过粗略查看之后,将其运用到课题中,这里将此...bool CLAHE_Algorithm::CLAHE_Process(cv::Mat& src_img) { if (!src_img.data) { std::cout 无图像数据" ; retur

    1. 概述

    本篇文章是基于这篇博文写的,然后经过粗略查看之后,将其运用到课题中,这里将此记录下来作为记录,文章中若有错误的地方敬请谅解。

    2. 实现代码

    bool CLAHE_Algorithm::CLAHE_Process(cv::Mat& src_img)
    {
    	if (!src_img.data)
    	{
    		std::cout << "无图像数据" << endl;
    		return false;
    	}
    
    	int rows(src_img.rows);
    	int cols(src_img.cols);
    
    	int img_type = src_img.type();	//图像的类型
    	if (CV_8UC1 != img_type && CV_8UC3 != img_type)
    	{
    		//cout << "输入8位的图像" << endl;
    		return false;
    	}
    	int img_channels = src_img.channels();
    	unsigned char* data_array = nullptr;
    	unsigned char* data = nullptr;
    
    	if (3 == img_channels)
    	{
    		std::vector<cv::Mat> splited_img;
    		cv::Mat src_img_temp = src_img.clone();
    		cv::resize(src_img_temp, src_img_temp, cv::Size(512,512));
    		rows=src_img_temp.rows;
    		cols=src_img_temp.cols;
    		cv::split(src_img_temp, splited_img);
    		for (int i = 0; i < 3; i++)
    		{
    			delete[] data_array;
    			data_array = nullptr;
    
    			data_array = new unsigned char[rows*cols];
    
    			cv::Mat& temp = splited_img[i];
    			for (int i = 0; i < rows; i++)
    			{
    				data = temp.ptr<unsigned char>(i);
    				for (int j = 0; j < cols; j++)
    				{
    					data_array[i*rows + j] = *data++;
    				}
    			}
    			this->CLAHE(data_array, cols, rows, 0, 255, 8, 8, 128, 10.0);
    
    			//cv::Mat result(temp.rows, temp.cols, temp.type(), cv::Scalar::all(0));
    			for (int i = 0; i < rows; i++)
    			{
    				data = temp.ptr<unsigned char>(i);
    				for (int j = 0; j < cols; j++)
    				{
    					*data++ = data_array[i*rows + j];
    				}
    			}
    		}
    		cv::Mat& blue_img = splited_img[0];	//blue channel
    		cv::Mat& green_img = splited_img[1];	//green channel
    		cv::Mat& red_img = splited_img[2];	//red channel
    		cv::Mat result;
    		cv::merge(splited_img, result);
    		result.copyTo(src_img);
    	}
    	else if (1== img_channels)
    	{
    		data_array = new unsigned char[rows*cols];
    
    		for (int i = 0; i < rows; i++)
    		{
    			data = src_img.ptr<unsigned char>(i);
    			for (int j=0; j<cols; j++)
    			{
    				data_array[i*rows+j] = *data++;
    			}
    		}
    		this->CLAHE(data_array, cols, rows, 0, 255, 8, 8, 128, 10.0);
    
    		for (int i = 0; i < rows; i++)
    		{
    			data = src_img.ptr<unsigned char>(i);
    			for (int j = 0; j < cols; j++)
    			{
    				*data++ = data_array[i*rows + j];
    			}
    		}
    	}
    
    	if (data_array)
    	{
    		delete[] data_array;
    		data_array = nullptr;
    	}
    
    	return true;
    }
    
    
    //************************************************************************
    // 函数名称:    	CLAHE
    // 访问权限:    	private 
    // 创建日期:		2016/11/27
    // 创 建 人:		
    // 函数说明:		CLAHE算法的主函数
    // 函数参数: 	kz_pixel_t * pImage		输入的图像数据
    // 函数参数: 	unsigned int uiXRes		图像X轴上的分辨率
    // 函数参数: 	unsigned int uiYRes		图像Y轴上的分辨率
    // 函数参数: 	kz_pixel_t Min			输入图像也是输出图像的最小像素值
    // 函数参数: 	kz_pixel_t Max			输入图像也是输出图像的最大像素值
    // 函数参数: 	unsigned int uiNrX		输入图像在X轴方向上进行分块的数目(min 2, max uiMAX_REG_X)
    // 函数参数: 	unsigned int uiNrY		输入图像在Y轴方向上进行分块的数目(min 2, max uiMAX_REG_X)
    // 函数参数: 	unsigned int uiNrBins	输出图像的有效灰度级数,输出图像同样拥有之前设定好的最大和最小灰度值,通常选择128既能加快处理速度并且效果不错
    // 函数参数: 	float fCliplimit		图像规范化区间限制,数值越大图像对比度越高(当为1时直接返回输入图像)
    // 返 回 值:   	int $
    //************************************************************************
    int CLAHE_Algorithm::CLAHE (kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes, kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrX,
    	unsigned int uiNrY, unsigned int uiNrBins, float fCliplimit)
    {
        unsigned int uiX, uiY;							/* counters */
        unsigned int uiXSize, uiYSize, uiSubX, uiSubY;	/* size of context. reg. and subimages */
        unsigned int uiXL, uiXR, uiYU, uiYB;			/* auxiliary variables interpolation routine */
        unsigned long ulClipLimit, ulNrPixels;			/* clip limit and region pixel count */
        kz_pixel_t* pImPointer;							/* pointer to image */
        kz_pixel_t aLUT[uiNR_OF_GREY];					/* lookup table used for scaling of input image */
        unsigned long* pulHist, *pulMapArray;			/* pointer to histogram and mappings*/
        unsigned long* pulLU, *pulLB, *pulRU, *pulRB;	/* auxiliary pointers interpolation */
    
        if (uiNrX > uiMAX_REG_X) return -1;				/* # of regions x-direction too large */
        if (uiNrY > uiMAX_REG_Y) return -2;				/* # of regions y-direction too large */
        if (uiXRes % uiNrX) return -3;					/* x-resolution no multiple of uiNrX */
        if (uiYRes & uiNrY) return -4;					/* y-resolution no multiple of uiNrY */
        if (Max >= uiNR_OF_GREY) return -5;				/* maximum too large */
        if (Min >= Max) return -6;						/* minimum equal or larger than maximum */
        if (uiNrX < 2 || uiNrY < 2) return -7;			/* at least 4 contextual regions required */
        if (fCliplimit == 1.0) return 0;				/* is OK, immediately returns original image. */
        if (uiNrBins == 0) uiNrBins = 128;				/* default value when not specified */
    
        pulMapArray=(unsigned long *)malloc(sizeof(unsigned long)*uiNrX*uiNrY*uiNrBins);	//全部块的直方图查询表
        if (pulMapArray == 0) return -8;				/* Not enough memory! (try reducing uiNrBins) */
    
        uiXSize = uiXRes/uiNrX; uiYSize = uiYRes/uiNrY;  /* Actual size of contextual regions */
        ulNrPixels = (unsigned long)uiXSize * (unsigned long)uiYSize;
    
        if(fCliplimit > 0.0) {							/* Calculate actual cliplimit     */
           ulClipLimit = (unsigned long) (fCliplimit * (uiXSize * uiYSize) / uiNrBins);
           ulClipLimit = (ulClipLimit < 1UL) ? 1UL : ulClipLimit;
        }	//计算直方图限制值		
        else ulClipLimit = 1UL<<14;						/* Large value, do not clip (AHE) */
        MakeLut(aLUT, Min, Max, uiNrBins);				/* Make lookup table for mapping of greyvalues */
        /* Calculate greylevel mappings for each contextual region */
    	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);							//计算最后的直方图的映射lookup表
    		}
    		pImPointer += (uiYSize - 1) * uiXRes;          /* skip lines, set pointer */
    	}
    
        /* Interpolate greylevel mappings to get CLAHE image */
    	for (pImPointer = pImage, uiY = 0; uiY <= uiNrY; uiY++)
    	{
    		if (uiY == 0) {                      /* special case: top row */
    			uiSubY = uiYSize >> 1;  uiYU = 0; uiYB = 0;
    		}
    		else
    		{
    			if (uiY == uiNrY) {                  /* special case: bottom row */
    				uiSubY = uiYSize >> 1;    uiYU = uiNrY - 1;     uiYB = uiYU;
    			}
    			else {                      /* default values */
    				uiSubY = uiYSize; uiYU = uiY - 1; uiYB = uiYU + 1;
    			}
    		}
    		for (uiX = 0; uiX <= uiNrX; uiX++)
    		{
    			if (uiX == 0)
    			{                  /* special case: left column */
    				uiSubX = uiXSize >> 1; uiXL = 0; uiXR = 0;
    			}
    			else
    			{
    				if (uiX == uiNrX)
    				{              /* special case: right column */
    					uiSubX = uiXSize >> 1;  uiXL = uiNrX - 1; uiXR = uiXL;
    				}
    				else
    				{              /* default values */
    					uiSubX = uiXSize; uiXL = uiX - 1; uiXR = uiXL + 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 */
    	return 0;									/* return status OK */
    }
    
    //************************************************************************
    // 函数名称:    	ClipHistogram
    // 访问权限:    	private static 
    // 创建日期:		2016/11/27
    // 创 建 人:		
    // 函数说明:		直方图修剪This function performs clipping of the histogram and redistribution of bins.
    //				The histogram is clipped and the number of excess pixels is counted.Afterwards
    //				the excess pixels are equally redistributed across the whole histogram(providing
    //				the bin count is smaller than the cliplimit
    // 函数参数: 	unsigned long * pulHistogram	统计直方图数据
    // 函数参数: 	unsigned int uiNrGreylevels		输出灰度等级
    // 函数参数: 	unsigned long ulClipLimit		限制值
    // 返 回 值:   	void $
    //************************************************************************
    void CLAHE_Algorithm::ClipHistogram (unsigned long* pulHistogram, unsigned int uiNrGreylevels, unsigned long ulClipLimit)
    {
        unsigned long* pulBinPointer, *pulEndPointer, *pulHisto;
        unsigned long ulNrExcess, 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;
    			//}
    			//应修正为
    			if (pulHistogram[i] > ulUpper)
    			{       /* high bin count */
    				ulNrExcess -= (ulClipLimit - pulHistogram[i]); pulHistogram[i] = ulClipLimit;
    			}
    			else
    			{                    /* low bin count */
    				ulNrExcess -= ulBinIncr;
    				pulHistogram[i] += ulBinIncr;
    			}
    		}
    	}
    
    	while (ulNrExcess)
    	{   /* Redistribute remaining excess  */
    		pulEndPointer = &pulHistogram[uiNrGreylevels]; pulHisto = pulHistogram;
    
    		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 */
    		}
    	}
    }
    
    //************************************************************************
    // 函数名称:    	MakeHistogram
    // 访问权限:    	private static 
    // 创建日期:		2016/11/27
    // 创 建 人:		
    // 函数说明:		
    // 函数参数: 	kz_pixel_t * pImage			输入的图像数据
    // 函数参数: 	unsigned int uiXRes			输入图像X轴方向的分辨率
    // 函数参数: 	unsigned int uiSizeX		输入图像经过分割之后拆分得到的单个块X轴的分辨率
    // 函数参数: 	unsigned int uiSizeY		输入图像经过分割之后拆分得到的单个块Y轴的分辨率
    // 函数参数: 	unsigned long * pulHistogram 输出统计直方图
    // 函数参数: 	unsigned int uiNrGreylevels	输出图像的灰度级数
    // 函数参数: 	kz_pixel_t * pLookupTable	查询表
    // 返 回 值:   	void $
    //************************************************************************
    void CLAHE_Algorithm::MakeHistogram (kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiSizeX, unsigned int uiSizeY,
            unsigned long* pulHistogram, unsigned int uiNrGreylevels, kz_pixel_t* pLookupTable)
    /* This function classifies the greylevels present in the array image into
     * a greylevel histogram. The pLookupTable specifies the relationship
     * between the greyvalue of the pixel (typically between 0 and 4095) and
     * the corresponding bin in the histogram (usually containing only 128 bins).
     */
    {
        kz_pixel_t* pImagePointer;
        unsigned int i;
    
        for (i = 0; i < uiNrGreylevels; i++) 
    		pulHistogram[i] = 0L; /* clear histogram */
    
    	for (i = 0; i < uiSizeY; i++)
    	{
    		pImagePointer = &pImage[uiSizeX];
    		while (pImage < pImagePointer)
    			pulHistogram[pLookupTable[*pImage++]]++;	//统计图像的灰度直方图
    		pImagePointer += uiXRes;						//将图像的指针转换到下一行像素数据
    		//pImage = &pImagePointer[-uiSizeX];
    		pImage = &(pImagePointer-uiSizeX)[0];
    	}
    }
    
    //************************************************************************
    // 函数名称:    	MapHistogram
    // 访问权限:    	private static 
    // 创建日期:		2016/11/27
    // 创 建 人:		
    // 函数说明:		计算映射的直方图This function calculates the equalized lookup table (mapping) by
    //				cumulating the input histogram.Note: lookup table is rescaled in range[Min..Max].
    // 函数参数: 	unsigned long * pulHistogram	直方图
    // 函数参数: 	kz_pixel_t Min					图像的最小灰度值
    // 函数参数: 	kz_pixel_t Max					图像的最大灰度值
    // 函数参数: 	unsigned int uiNrGreylevels		输出图像的灰度级数
    // 函数参数: 	unsigned long ulNrOfPixels		单个块的总像素点个数
    // 返 回 值:   	void $
    //************************************************************************
    void CLAHE_Algorithm::MapHistogram (unsigned long* pulHistogram, kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrGreylevels,
    	unsigned long ulNrOfPixels)
    {
        unsigned int i;  unsigned long ulSum = 0;
        const float fScale = ((float)(Max - Min)) / ulNrOfPixels;
        const unsigned long ulMin = (unsigned long) Min;
    
        for (i = 0; i < uiNrGreylevels; i++) 
    	{
    		ulSum += pulHistogram[i]; 
    		pulHistogram[i]=(unsigned long)(ulMin+ulSum*fScale);
    		if (pulHistogram[i] > Max) 
    			pulHistogram[i] = Max;
        }
    }
    
    //************************************************************************
    // 函数名称:    	MakeLut
    // 访问权限:    	private static 
    // 创建日期:		2016/11/27
    // 创 建 人:		
    // 函数说明:		To speed up histogram clipping, the input image [Min,Max] is scaled down to*[0, uiNrBins - 1].
    //				This function calculates the LUT.
    // 函数参数: 	kz_pixel_t * pLUT	查询表
    // 函数参数: 	kz_pixel_t Min		输入图像的最小值
    // 函数参数: 	kz_pixel_t Max		输出图像的最大值
    // 函数参数: 	unsigned int uiNrBins	输出图像指定的灰度级数
    // 返 回 值:   	void $
    //************************************************************************
    void CLAHE_Algorithm::MakeLut (kz_pixel_t * pLUT, kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrBins)
    {
        int i;
        const kz_pixel_t BinSize = (kz_pixel_t) (1 + (Max - Min) / uiNrBins);
    
        for (i = Min; i <= Max; i++)
    		pLUT[i] = (i - Min) / BinSize;
    }
    
    void CLAHE_Algorithm::Interpolate (kz_pixel_t * pImage, int uiXRes, unsigned long * pulMapLU, unsigned long * pulMapRU, unsigned long * pulMapLB,
    	unsigned long * pulMapRB, unsigned int uiXSize, unsigned int uiYSize, kz_pixel_t * pLUT)
    /* pImage      - pointer to input/output image
     * uiXRes      - resolution of image in x-direction
     * pulMap*     - mappings of greylevels from histograms
     * uiXSize     - uiXSize of image submatrix
     * uiYSize     - uiYSize of image submatrix
     * pLUT           - lookup table containing mapping greyvalues to bins
     * This function calculates the new greylevel assignments of pixels within a submatrix
     * of the image with size uiXSize and uiYSize. This is done by a bilinear interpolation
     * between four different mappings in order to eliminate boundary artifacts.
     * It uses a division; since division is often an expensive operation, I added code to
     * perform a logical shift instead when feasible.
     */
    {
        const unsigned int uiIncr = uiXRes-uiXSize; /* Pointer increment after processing row */
        kz_pixel_t GreyValue; unsigned int uiNum = uiXSize*uiYSize; /* Normalization factor */
    
        unsigned int 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++ = (kz_pixel_t)((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++ = (kz_pixel_t)((uiYInvCoef* (uiXInvCoef * pulMapLU[GreyValue]
    					+ uiXCoef * pulMapRU[GreyValue])
    					+ uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
    					+ uiXCoef * pulMapRB[GreyValue])) >> uiShift);
    			}
    		}
    	}
    }

    3. 实现效果



    展开全文
  • CLAHE算法学习

    千次阅读 2019-12-30 19:49:12
    图像识别工程开发中需要增强图像对比度,便于后续处理,接触到了CLAHE(Contrast Limited Adaptive Histogram Equalization),记录一下其中的学习过程。 1.直方图均衡 1.1灰度直方图 灰度图中像素值的分布为0-...

    0.前言

    图像识别工程开发中需要增强图像对比度,便于后续处理,接触到了CLAHE(Contrast Limited Adaptive Histogram Equalization),记录一下其中的学习过程。

    1.直方图均衡

    1.1灰度直方图

    灰度图中像素值的分布为0-255,以灰度值为横坐标,纵坐标为该灰度值对应的像素点数目/比例,则得到了灰度图像的直方图,体现的是图像中灰度的整体分布情况。

    OpenCV中提供了相应的计算函数calcHist(),可以得到相应分布范围的像素点数;将其绘制出来观察,更为直观。下图为某输入图像和其灰度直方图分布。

    一般性规律:对于灰度图像,以个人观察意愿将其分为前景和背景区域,前景是我们所感兴趣的区域,背景则反之。若前景和背景的灰度差异大,则易于人们观察,否则不易观察。一般来说若图像的灰度直方图集中在某个区域,其整体对比度较差,不易于我们区分前景和背景;反之若灰度直方图分布较为均匀,则整体对比度较好,易于我们区分前景和背景。

    1.2直方图均衡化(Histogram Equalization, HE)

    1)效果

    直方图均衡化就是一种使图像的灰度直方图分布更佳均匀的变换方法。OpenCV中的equalizeHist()实现了该功能

    左侧为原始灰度图,中间上方为其灰度分布;右侧为HE处理后得到的图像,中间下方为其灰度分布;对比灰度分布可知右方图像的灰度分布更为均匀,图像效果上来看结果图像对比度更强。

    2)数学原理

    先说概率分布与概率密度函数,看一下下图中的题目回想一下以前的数学知识。

    接下来对概率密度函数的变换做一些数学推导,输入公式比较麻烦,截图为图像粘贴进来。

    这里可以看到转换函数T()为单调递增函数,不过结果图的直方图分布并不是理想的一条水平线;原因在于我们是在连续分布上推导得到转换函数,作用于离散分布,且不允许产生非整数形式的新结果。不过整体看来该转换函数使新的分布具有展开直方图分布的效果,增强了整体的对比度。

    2.CLAHE

    上述的HE算法得到的结果存在一些问题,主要包括:1)部分区域由于对比度增强过大,成为噪点;2)一些区域调整后变得更暗/更亮,丢失细节信息。

    2.1CLHE

    针对问题1),有人提出了CLHE,即加入对比度限制的HE算法。转换函数T受像素灰度概率分布的影响:其变化率与概率密度成比例;概率密度过大,会导致原始像素经过变换后灰度值过高,对比度增强过大。

    因此需要限制对比度,即限制灰度的分布;同时还需要确保概率密度的积分仍然为1。

    如下图所示,一般操作为设置直方图分布的阈值,将超过该阈值的分布“均匀”分散至概率密度分布上,由此来限制转换函数(累计直方图)的增幅。

    下图是使用CLHE和HE的对比,左侧图像为HE结果,中间上方图像为其灰度直方图;右侧图像为CLHE结果,中间下方图像为其灰度直方图;观察灰度直方图可看到上述的限制效果;同时CLHE结果中噪点的分布相对于HE结果要小一些。

    2.2 AHE

    针对问题2),有人提出自适应直方图均衡算法(AHE),基本思想为:对原图中每个像素,计算其周围一个邻域内的直方图,并使用HE映射得到新的像素值。为了能够处理图像边缘的像素,一般先对原始图像做镜像扩边处理。

    不过这样做计算量偏大,为了减少计算量,一般先把原始图像分为MXN个子区域,分别对每个子区域进行HE变换。但是这样又产生了新的问题,子区域相接处像素值分布不连续,产生很强的分割线,如下图所示。

    为了解决该问题,有人提出加入双线性插值的AHE,也就形成了完整的CLAHE。算法流程如下:

    1)将图像分为多个矩形块大小,对于每个矩形块子图,分别计算其灰度直方图和对应的变换函数(累积直方图)

    2)将原始图像中的像素按照分布分为三种情况处理;对比下图,红色区域中的像素按照其所在子图的变换函数进行灰度映射,绿色区域中的像素按照所在的两个相邻子图变换函数变换后进行线性插值得到;紫色区域中的像素按照其所在的四个相邻子图变换函数变换后双线性插值得到

    Clahe-tileinterpol.svg

    3. OpenCV中的CLAHE

    OpenCV中提供了CLAHE的实现和接口,调用也十分简单,其提供的两个设置接口分别用于设置自适应处理中的阈值,和子图的大小。下面是调用OpenCV的HE和CLAHE处理的示例代码和结果图像。

    void opencvEqualizeHist()
    {
    	Mat img, img_gray, result;
    	Mat clahe_result;
    	img = imread("test.jpg");
    	cvtColor(img, img_gray, CV_BGR2GRAY);
    	equalizeHist(img_gray, result);
    
    	cv::Ptr<CLAHE> clahe = cv::createCLAHE();
    	clahe->setClipLimit(4);
    	clahe->setTilesGridSize(cv::Size(10, 10));
    	clahe->apply(img_gray, clahe_result);
    
    	imshow("src", img_gray);
    	imshow("result", result);
    	imshow("clahe_result", clahe_result);
    
    	waitKey();
    }
    

    对比可发现CLAHE的结果相对于HE的结果,未出现结果像素相比于原图因调整得更亮/更暗而丢失了细节的情况,且局部对比度更强。

    4.参考

    https://www.cnblogs.com/aoru45/p/10633313.html

    https://blog.csdn.net/superjunenaruto/article/details/52431941

    https://www.jianshu.com/p/5e8392ec1940

    https://www.cnblogs.com/Imageshop/archive/2013/04/07/3006334.html

    https://blog.csdn.net/eternity1118_/article/details/51492105

    https://www.cnblogs.com/jsxyhelu/p/6435601.html

    展开全文
  • 包含一篇对CLAHE的原理阐述论文,以及具体实现的MATLAB代码
  • 基于直方图的图像增强算法(HE、CLAHE)之(二)

    万次阅读 多人点赞 2016-02-15 10:29:38
    作为图像增强算法系列的第二篇文章,下面我们将要介绍功能强大、用途广泛、影响深远的对比度有限的自适应直方图均衡(CLAHE,Contrast Limited Adaptive Histogram Equalization)算法。尽管最初它仅仅是被当作一种...

    作为图像增强算法系列的第二篇文章,下面我们将要介绍功能强大、用途广泛、影响深远的对比度有限的自适应直方图均衡(CLAHE,Contrast Limited Adaptive Histogram Equalization)算法。尽管最初它仅仅是被当作一种图像增强算法被提出,但是现今在图像去雾、低照度图像增强,水下图像效果调节、以及数码照片改善等方面都有应用。这个算法的算法原理看似简单,但是实现起来却并不那么容易。我们将结合相应的Matlab代码来对其进行解释。希望你在阅读本文之前对朴素的直方图均衡算法有所了解,相关内容可以参见本系列的前一篇文章:基于直方图的图像增强算法(HE、CLAHE)之(一)(http://blog.csdn.net/baimafujinji/article/details/50614332)。

     

    先来看一下待处理的图像效果:

     

    下面是利用CLAHE算法处理之后得到的两个效果(后面我们还会具体介绍我们所使用的策略):

     

       

                                 效果图A                                                                         效果图B

    对于一幅图像而言,它不同区域的对比度可能差别很大。可能有些地方很明亮,而有些地方又很暗淡。如果采用单一的直方图来对其进行调整显然并不是最好的选择。于是人们基于分块处理的思想提出了自适应的直方图均衡算法AHE。维基百科上说的也比较明白:AHE improves on this by transforming each pixel with a transformation function derived from a neighbourhood region. 但是这种方法有时候又会将一些噪声放大,这是我们所不希望看到的。于是荷兰乌得勒支大学的Zuiderveld教授又引入了CLAHE,利用一个对比度阈值来去除噪声的影响。特别地,为了提升计算速度以及去除分块处理所导致的块边缘过渡不平衡效应,他又建议采用双线性插值的方法。关于算法的介绍和描述,下面这两个资源已经讲得比较清楚。

    [1] https://en.wikipedia.org/wiki/Adaptive_histogram_equalization#Contrast_Limited_AHE

    [2] K. Zuiderveld: Contrast Limited Adaptive Histogram Equalization. In: P. Heckbert: Graphics Gems IV, Academic Press 1994 (http://www.docin.com/p-119164091.html)

    事实上,尽管这个算法原理,然而它实现起来却仍然有很多障碍。但在此之前,笔者还需说明的是,Matlab中已经集成了实现CLAHE的函数adapthisteq(),如果你仅仅需要一个结果,其实直接使用这个函数就是最好的选择。我给出一些示例代码用以生成前面给出之效果。函数adapthisteq()只能用来处理灰度图,如果要处理彩色图像,则需要结合自己编写的代码来完成。上一篇文章介绍了对彩色图像进行直方图均衡的两种主要策略:一种是对R、G、B三个通道分别进行处理;另外一种是转换到另外一个色彩空间中再进行处理,例如HSV(转换后只需对V通道进行处理即可)。

    首先,我们给出对R、G、B三个通道分别使用adapthisteq()函数进行处理的示例代码:

    img = imread('space.jpg');
    rimg = img(:,:,1);
    gimg = img(:,:,2);
    bimg = img(:,:,3);
    resultr = adapthisteq(rimg);
    resultg = adapthisteq(gimg);
    resultb = adapthisteq(bimg);
    result = cat(3, resultr, resultg, resultb);
    imshow(result);

    上述程序之结果效果图A所示。

     

    下面程序将原图像的色彩空间转换到LAB空间之后再对L通道进行处理。

    clear;
    img = imread('space.jpg');
    cform2lab = makecform('srgb2lab');
    LAB = applycform(img, cform2lab);
    L = LAB(:,:,1);
    LAB(:,:,1) = adapthisteq(L);
    cform2srgb = makecform('lab2srgb');
    J = applycform(LAB, cform2srgb);
    imshow(J);
    


    上述程序所得之结果如图B所示。

     

    如果你希望把这个算法进一步提升和推广,利用用于图像去雾、低照度图像改善和水下图像处理,那么仅仅知其然是显然不够的,你还必须知其所以然。希望我下面一步一步实现的代码能够帮你解开这方面的困惑。鉴于前面所列之文献已经给出了比较详细的算法描述,下面将不再重复这部分内容,转而采用Matlab代码来对其中的一些细节进行演示。

     

    首先来从灰度图的CLAHE处理开始我们的讨论。为此清理一下Matlab的环境。然后,读入一张图片(并将其转化灰度图),获取图片的长、宽、像素灰度的最大值、最小值等信息。

    clc;
    clear all;
    Img = rgb2gray(imread('space.jpg'));
    [h,w] = size(Img);
    minV = double(min(min(Img)));
    maxV = double(max(max(Img)));
    imshow(Img);
    
    


    图像的初始状态显示如下。此外该图的 Height = 395,Width = 590,灰度最大值为255,最小值为8。

     

     

    我们希望把原图像水平方向分成8份,把垂直方向分成4份,即原图将被划分成4 × 8 = 32个SubImage。然后可以算得每个块(tile)的height = 99,width = 74。注意,由于原图的长、宽不太可能刚好可被整除,所以我在这里的处理方式是建立一个稍微大一点的图像,它的宽和长都被补上了deltax和deltay,以保证长、宽都能被整除。

    NrX = 8;
    NrY = 4;
    HSize = ceil(h/NrY);
    WSize = ceil(w/NrX);
    
    deltay = NrY*HSize - h;
    deltax = NrX*WSize - w;
    
    tmpImg = zeros(h+deltay,w+deltax);
    tmpImg(1:h,1:w) = Img;
    

     

    对长和宽进行填补之后,对新图像的一些必要信息进行更新。

    new_w = w + deltax;
    new_h = h + deltay;
    
    NrPixels = WSize * WSize;


    然后指定图像中直方图横坐标上取值的计数(也就指定了统计直方图上横轴数值的间隔或计数的精度),对于色彩比较丰富的图像,我们一般都要求这个值应该大于128。

    % NrBins - Number of greybins for histogram ("dynamic range")
    NrBins = 256;

    然后用原图的灰度取值范围重新映射了一张Look-Up Table(当然你也可以直接使用0~255这个范围,这取决你后续建立直方图的具体方法),并以此为基础为每个图像块(tile)建立直方图。

    LUT = zeros(maxV+1,1);
    
    for i=minV:maxV
        LUT(i+1) = fix(i - minV);%i+1
    end
    
    Bin = zeros(new_h, new_w);
    for m = 1 : new_h
        for n = 1 : new_w
            Bin(m,n) = 1 + LUT(tmpImg(m,n) + 1);
        end
    end
    
    Hist = zeros(NrY, NrX, 256);
    for i=1:NrY
        for j=1:NrX
            tmp = uint8(Bin(1+(i-1)*HSize:i*HSize, 1+(j-1)*WSize:j*WSize));
            %tmp = tmpImg(1+(i-1)*HSize:i*HSize,1+(j-1)*WSize:j*WSize);
            [Hist(i, j, :), x] = imhist(tmp, 256);
        end
    end
    
    Hist = circshift(Hist,[0, 0, -1]);


    注意:按通常的理解,上面这一步我们应该建立的直方图(集合)应该是一个4×8=32个长度为256的向量(你当然也可以这么做)。但由于涉及到后续的一些处理方式,我这里是生成了一个长度为256的4×8矩阵。Index = 1的矩阵其实相当于是整张图像各个tile上灰度值=0的像素个数计数。例如,我们所得的Hist(:, :, 18)如下。这就表明图像中最左上角的那个tile里面灰度值=17的像素有零个。同理,它右边的一个tile则有46个灰度值=17的像素。

    Hist(:,:,18) =
    
         0    46   218    50    14    55    15     7
         0     0    21    18   114    15    74    73
         0     1     0     0     2    67   124    82
         0     0     0     0     0     1     9     2

     

    然后来对直方图进行裁剪。Matlab中内置的函数adapthisteq()中cliplimit参数的取值范围是0~1。这里我们所写的方法则要求该值>1。当然这完全取决你算法实现的策略,它们本质上并没有差异。然后我们将得到新的(裁剪后的)映射直方图。

    ClipLimit = 2.5;
    ClipLimit = max(1,ClipLimit * HSize * WSize/NrBins);
    
    Hist = clipHistogram(Hist,NrBins,ClipLimit,NrY,NrX);
    Map=mapHistogram(Hist, minV, maxV, NrBins, NrPixels, NrY, NrX);

    因为这里没有具体给出clipHistogram函数的实现,所以此处我希望插入一部分内容来解释一下我的实现策略(也就是说,在实际程序中并不需要包含这部分)。我们以图像最左上角的一个tile为例,它的原直方图分布可以用下面代码来绘出:

    tmp_hist = reshape(Hist(1,1,:), 1, 256);
    plot(tmp_hist)

    输出结果下图中的左图所示。

     

    如果我们给ClipLimit赋初值为2.5,则经过语句ClipLimit = max(1,ClipLimit * HSize * WSize/NrBins);计算之后,ClipLimit将变成71.54。然后我们再用上述代码绘制新的直方图,其结果将如上图中的右图所示。显然,图中大于71.54的部分被裁剪掉了,然后又平均分配给整张直方图,所以你会发现整张图都被提升了。这就是我们这里进行直方图裁剪所使用的策略。但是再次强调,matlab中的内置函数adapthisteq()仅仅是将这个参数进行了归一化,这与我们所使用的方法并没有本质上的区别。

     

    继续回到程序实现上的讨论。最后,也是最关键的步骤,我们需要对结果进程插值处理。这也是Zuiderveld设计的算法中最复杂的部分。

    yI = 1;
    for i = 1:NrY+1
        if i == 1
            subY = floor(HSize/2);
            yU = 1;
            yB = 1;
        elseif i == NrY+1
            subY = floor(HSize/2);
            yU = NrY;
            yB = NrY;
        else
            subY = HSize;
            yU = i - 1;
            yB = i;
        end
        xI = 1;
        for j = 1:NrX+1
            if j == 1
                subX = floor(WSize/2);
                xL = 1;
                xR = 1;
            elseif j == NrX+1
                subX = floor(WSize/2);
                xL = NrX;
                xR = NrX;
            else
                subX = WSize;
                xL = j - 1;
                xR = j;
            end
            UL = Map(yU,xL,:);
            UR = Map(yU,xR,:);
            BL = Map(yB,xL,:);
            BR = Map(yB,xR,:);
            subImage = Bin(yI:yI+subY-1,xI:xI+subX-1);
    
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            sImage = zeros(size(subImage));
            num = subY * subX;
            for i = 0:subY - 1
                inverseI = subY - i;
                for j = 0:subX - 1
                    inverseJ = subX - j;
                    val = subImage(i+1,j+1);
                    sImage(i+1, j+1) = (inverseI*(inverseJ*UL(val) + j*UR(val)) ...
                                    + i*(inverseJ*BL(val) + j*BR(val)))/num;
                end
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            output(yI:yI+subY-1, xI:xI+subX-1) = sImage;
            xI = xI + subX;
        end
        yI = yI + subY;
    end


    这个地方,作者原文中已经讲得比较清楚了,我感觉我也没有必要狗尾续貂,班门弄斧了。下面截作者原文中的一段描述,足以说明问题。

    最后来看看我们处理的效果如何(当然,这里还需要把之前我们填补的部分裁掉)。

    output = output(1:h, 1:w);
    figure, imshow(output, []);


    来看看结果吧~可以对比一下之前的灰度图,不难发现,图像质量已有大幅改善。

    更多有趣有用的图像处理算法还可以参考我的《数字图像中的数学修炼》

     

    展开全文
  • CLAHE

    千次阅读 2017-06-29 09:30:21
    实现了基于OpenCV的CLAHE实现和研究。从最基本的开始做,分别实现HE算法,AHE算法,CLHE算法和CLAHE算法。同时调用OpenCV生成代码和自己编写代码进行比对。以别人的灰度图像例子,自己研究彩色处理。int _tmain(int ...
  • CLAHE

    2020-07-08 23:51:20
    直方图均衡 目标:学习直方图均衡化的概念,并将其用于改善图像的对比度。 理论: 考虑一个图像,其像素值仅限于特定的值范围。例如,较亮的图像会将所有像素限制在较高的值。但是,好的图像将具有来自所有区域的...
  • %% tic %% 清空工作区与变量 clc; clear; for image_number=1:1 imageName=strcat(num2str(image_number),'.jpg'); img = imread(imageName); %% 在LAB空间进行去雾 % RGB转LAB transform = makecform('... ...
  • CLAHE

    千次阅读 2019-06-12 10:18:37
    CLAHE同普通的自适应直方图均衡不同的地方主要是其对比度限幅。这个特性也可以应用到全局直方图均衡化中,即构成所谓的限制对比度直方图均衡(CLHE),但这在实际中很少使用。在CLAHE中,对于每个小区域都必须使用...
  • CLAHE

    千次阅读 2013-11-28 16:08:25
    #include #include "opencv2/core/core.hpp"#include "opencv2/core/utility.hpp"#include "opencv2/imgproc/imgproc.hpp"#include "opencv2/highgui/highgui.hpp"#include "opencv2/ocl/ocl.hpp"using namespace cv
  • CLAHE

    千次阅读 2015-03-30 16:54:57
     CLAHE同普通的自适应直方图均衡不同的地方主要是其对比度限幅。这个特性也可以应用到全局直方图均衡化中,即构成所谓的限制对比度直方图均衡(CLHE),但这在实际中很少使用。在CLAHE中,对于每个小区域都必须使用...
  • CLAHE 代码

    热门讨论 2015-11-02 13:45:59
    Matlab编写的CLAHE,非常值得借鉴。
  • CLAHE

    2020-09-17 22:36:00
    https://blog.csdn.net/jingbo18/article/details/81707181
  • CLAHE代码

    千次阅读 2016-01-14 18:59:19
    为了提高运行速度,把CLAHE的代码进行了修改,方便运行。 CLAHE.h class CLAHE { public: CLAHE(void); ~CLAHE(void); int m_nGridX; /* CLAHE 宽网格的个数 */ int m_nGridY; /* CLAHE 高网格的个数 ...
  • CLAHE

    2018-07-18 16:42:00
    CLAHE 读音:clay Contrast Limiting Adaptive Histogram Equalization 用来处理灰度图像非常好用。 一种理解,相当于样本均衡,对于不要求绝对回归的任务中可以使用。 转载于:...
  • CLAHE算法对于医学图像,特别是医学红外图像的增强效果非常明显。 CLAHE https://en.wikipedia.org/wiki/Adaptive_histogram_equalization 中文方面非常好的资料 限制对比度自适应直方图均衡化算法原理、实现及效果 ...
  • 论文:Contrast limited adaptive histogram equalization Zuiderveld, Karel. "Contrast limited adaptive histogram equalization." Graphics gems IV. Academic Press Professional, Inc., 1994. ...
  • OpenCV源码的本地路径: %OPENCV%\opencv\sources\modules\imgproc\src\clahe.cpp clahe.cpp // ---------------------------------------------------------------------- // CLAHE namespace { class C
  • python实现clahe对比度增强

    千次阅读 2019-12-24 21:51:41
    python实现clahe对比度增强一、获取直方图二、累积分布函数三、双线性插值四、主函数五、使用的库六、与官方函数的对比 一、获取直方图 def hist_cal(img): img = np.array(img) heigh,width=img.shape histogram...
  • cv::CLAHE 使用CLAHE算法例子

    万次阅读 2014-04-01 20:08:38
    具体算法可以看源码 #include using namespace cv; int _tmain(int argc, _TCHAR* argv[]) ... Mat src = imread("d:/1.... Ptr clahe = createCLAHE(); Mat dst; clahe ->apply(src,dst); imshow("src",src);

空空如也

1 2 3 4 5 ... 20
收藏数 1,096
精华内容 438
关键字:

clahe