对⽐度受限⾃适应直⽅图均衡化⽅法
图像增强技术可具体为时域和频域。 时域增强往往被应⽤于提⾼图像的对⽐度且改进其灰度级， 其原理是在灰度映射转化的基础上， 对像素进⾏直接性地处理； 频域增强的作⽤是以强化图像的低、 ⾼频来达到改善图像的平滑性和边缘为主，其原理是通过傅⾥叶转换的⽅式提升兴趣区的频率分量。
图像灰度分布不均的问题对于图像处理技术⽽⾔是⼀个影响较⼤的问题， 为了能解决这个问题， 基于时域增强理论， 采⽤映射函数当中的⾮线性函数， 通过转化后使其像素分布均匀。
直⽅图均衡化⽅法通过不同的灰度分布⽅案使转变为灰度值均匀化分布的图像，以眼底图为例，利⽤直⽅图提取灰度映射曲线，再对其进⾏转换，从⽽提升亮度。但该⽅法仍有瑕疵，即⽆法增强其对⽐度，且在此过程当中，⼀并将噪声信号也同时放⼤了。以 AHE (Adaptive Histogtam Equalization)为例，该⽅法将图像利⽤⽹格线切割为数量众多的⼩格⼦区域， 并对每⼀个格⼦进⾏均衡化处理。但该⽅法会导致图像失真， 且其微噪也会因此⽽放⼤。 因此， 这⾥提出了 对⽐度受限制⾃适应直⽅图 (Contrast limited Adaptive Histogtam Equalization， CLAHE) ，即对每⼀个划分单元进⾏对⽐度的限制处理，这也是该⽅法与传统的AHE的不同之处。
CLAHE 通过限制局部直⽅图的⾼度来限制噪声放⼤和局部对⽐度增强。该⽅法将图像划分为多个⼦区域； 然后对每个⼦区域的直⽅图进⾏分类。 再对每个⼦区域分别进⾏直⽅图均衡化， 最后通过对每个
像素进⾏插值运算来获得变换后的灰度值，以此实现对⽐度受限⾃适应直⽅图均衡化图像增强。
CLAHE ⽅法原理：
CLAHE 与 AHE 不同的地⽅是增加对⽐度限幅，这就可以克服 AHE 的过度放⼤噪声的问题；
综上所述，改变最⼤的映射函数斜率Smax 及相应的最⼤直⽅图⾼度Hmax，可获得不同增强效果的图像。也就是说限制CDF的斜率就相当于限制Hist的幅度。
因此我们需要对⼦块中统计得到的直⽅图进⾏裁剪，使其幅值低于某个上限，当然裁剪掉的部分⼜不能扔掉，我们还需要将这部分裁剪值均匀地分布在整个灰度区间上，以保证直⽅图总⾯积不变，如下图：
可以看到，这时直⽅图⼜会整体上升了⼀个⾼度，貌似会超过我们设置的上限。其实在具体实现的时候有很多解决⽅法，你可以多重复⼏次裁剪过程，使得上升的部分变得微不⾜道，或是⽤另⼀种常⽤的⽅法：
设裁剪值为ClipLimit，求直⽅图中⾼于该值的部分的和totalExcess，此时假设将totalExcess均分给所有灰度级， 求出这样导致的直⽅图整体上升的⾼度L=totalExcess/N，以upper= ClipLimit-L为界限对直⽅图进⾏如下处理：
(1)若幅值⾼于ClipLimit，直接置为ClipLimit；
(2)若幅值处于Upper和ClipLimit之间，将其填补⾄ClipLimit；
(3)若幅值低于Upper，直接填补L个像素点；
经过上述操作，⽤来填补的像素点个数通常会略⼩于totalExcess，也就是还有⼀些剩余的像素点没分出去，这个剩余来⾃于(1)(2)两处。这时我们可以再把这些点均匀地分给那些⽬前幅值仍然⼩于ClipLimit的灰度值。
% total number of pixels overflowing clip limit in each bin
totalExcess = sum(max(imgHist - clipLimit,0));
% clip the histogram and redistribute the excess pixels in each bin
avgBinIncr = floor(totalExcess/numBins);
upperLimit = clipLimit - avgBinIncr; % bins larger than this will be
% set to clipLimit
% this loop should speed up the operation by putting multiple pixels
% into the "obvious" places first
for k=1:numBins
if imgHist(k) > clipLimit
imgHist(k) = clipLimit;
else
if imgHist(k) > upperLimit % high bin count
totalExcess = totalExcess - (clipLimit - imgHist(k));
imgHist(k) = clipLimit;
else
totalExcess = totalExcess - avgBinIncr;
imgHist(k) = imgHist(k) + avgBinIncr;
end
end
end
% this loops redistributes the remaining pixels, one pixel at a time
k = 1;
while (totalExcess ~= 0)
%keep increasing the step as fewer and fewer pixels remain for
%the redistribution (spread them evenly)
stepSize = max(floor(numBins/totalExcess),1);
for m=k:stepSize:numBins
if imgHist(m) < clipLimit
imgHist(m) = imgHist(m)+1;
totalExcess = totalExcess - 1; %reduce excess
if totalExcess == 0
break;
end
end
end
k = k+1; %prevent from always placing the pixels in bin #1
if k > numBins % start over if numBins was reached
k = 1;
end
end
CLAHE和AHE中另⼀个重要的问题：插值。
将图像进⾏分块处理，若每块中的像素点仅通过该块中的映射函数进⾏变换，则会导致最终图像呈块状效应：
为了解决这个问题，我们需要利⽤插值运算，也就是每个像素点出的值由它周围4个⼦块的映射函数值进⾏双线性插值得到。
上图中，为了求蓝⾊像素点处的值，需要利⽤它周围四个⼦块的映射函数分别做变换得到四个映射值，再对这四个值做双线性插值即可。
当然对于边界处的像素点则不是通过四个⼦块进⾏插值，如上图红⾊像素点直接以⼀个⼦块的映射函数做变换，绿⾊像素则以两个⼦块做映射函数做线性插值。这⾥讲的边界处像素是指落在图像左上⾓，左下⾓、右上⾓，右下⾓的四个⼦块中⼼像素点围成的四边形之外的像素。如下图，将图像分为8x8⼦块，边界像素即落在灰⾊区域的像素点。
CLAHE算法的源代码参考:
/*
* ANSI C code from the article
* "Contrast Limited Adaptive Histogram Equalization"
* by Karel Zuiderveld, karel@cv.ruu.nl
* in "Graphics Gems IV", Academic Press, 1994
*
*
* These functions implement Contrast Limited Adaptive Histogram Equalization.
* The main routine (CLAHE) expects an input image that is stored contiguously in
* memory; the CLAHE output image overwrites the original input image and has the
* same minimum and maximum values (which must be provided by the user).
* This implementation assumes that the X- and Y image resolutions are an integer
* multiple of the X- and Y sizes of the contextual regions. A check on various other
* error conditions is performed.
*
* #define the symbol BYTE_IMAGE to make this implementation suitable for
* 8-bit images. The maximum number of contextual regions can be redefined
* by changing uiMAX_REG_X and/or uiMAX_REG_Y; the use of more than 256
* contextual regions is not recommended.
*
* The code is ANSI-C and is also C++ compliant.
*
* Author: Karel Zuiderveld, Computer Vision Research Group,
* Utrecht, The Netherlands (karel@cv.ruu.nl)
*/
#ifdef BYTE_IMAGE
typedef unsigned char kz_pixel_t; /* for 8 bit-per-pixel images */
#define uiNR_OF_GREY (256)
#else
typedef unsigned short kz_pixel_t; /* for 12 bit-per-pixel images (default) */
# define uiNR_OF_GREY (4096)
#endif
/******** Prototype of CLAHE function. Put this in a separate include file. *****/
int 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);
/*********************** Local prototypes ************************/
static void ClipHistogram (unsigned long*, unsigned int, unsigned long);
static void ClipHistogram (unsigned long*, unsigned int, unsigned long);
static void MakeHistogram (kz_pixel_t*, unsigned int, unsigned int, unsigned int,
unsigned long*, unsigned int, kz_pixel_t*);
static void MapHistogram (unsigned long*, kz_pixel_t, kz_pixel_t,
unsigned int, unsigned long);
局部直方图均衡化static void MakeLut (kz_pixel_t*, kz_pixel_t, kz_pixel_t, unsigned int);
static void Interpolate (kz_pixel_t*, int, unsigned long*, unsigned long*,
unsigned long*, unsigned long*, unsigned int, unsigned int, kz_pixel_t*);
/************** Start of actual code **************/
#include <stdlib.h> /* To get prototypes of malloc() and free() */
const unsigned int uiMAX_REG_X = 16; /* max. # contextual regions in x-direction */
const unsigned int uiMAX_REG_Y = 16; /* max. # contextual regions in y-direction */
/************************** main function CLAHE ******************/
int 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)
/* pImage - Pointer to the input/output image
* uiXRes - Image resolution in the X direction
* uiYRes - Image resolution in the Y direction
* Min - Minimum greyvalue of input image (also becomes minimum of output image)
* Max - Maximum greyvalue of input image (also becomes maximum of output image)
* uiNrX - Number of contextial regions in the X direction (min 2, max uiMAX_REG_X)
* uiNrY - Number of contextial regions in the Y direction (min 2, max uiMAX_REG_Y)
* uiNrBins - Number of greybins for histogram ("dynamic range")
* float fCliplimit - Normalized cliplimit (higher values give more contrast)
* The number of "effective" greylevels in the output image is set by uiNrBins; selecting
* a small value (eg. 128) speeds up processing and still produce an output image of
* good quality. The output image will have the same minimum and maximum value as the input * image. A clip limit smaller than 1 results in standard (non-contrast limited) AHE.
*/
{
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) */
发表评论