一、自适应直方图均衡化(Adaptive histgram equalization/AHE)
1.简述
自适应直方图均衡化(AHE)用来提升图像的对比度的一种计算机图像处理技术。和普通的直方图均衡算法不同，AHE算法通过计算图像的局部直方图，然后重新分布亮度来来改变图像对比度。因此，该算法更适合于改进图像的局部对比度以及获得更多的图像细节。
不过，AHE有过度放大图像中相同区域的噪音的问题，另外一种自适应的直方图均衡算法即限制对比度直方图均衡（CLAHE）算法能有限的限制这种不利的放大。
2. 算法的解释
普通的直方图均衡算法对于整幅图像的像素使用相同的直方图变换，对于那些像素值分布比较均衡的图像来说，算法的效果很好。然后，如果图像中包括明显比图像其它区域暗或者亮的部分，在这些部分的对比度将得不到有效的增强。
AHE算法通过对局部区域执行响应的直方图均衡来改变上述问题。该算法首先被开发出来适用于改进航天器驾驶舱的显示效果。其最简单的形式，就是每个像素通过其周边一个矩形范围内的像素的直方图进行均衡化。均衡的方式则完全同普通的均衡化算法：变换函数同像素周边的累积直方图函数（CDF）成比例。
图像边缘的像素需要特殊处理，因为边缘像素的领域不完全在图像内部。这个通过镜像图像边缘的行像素或列像素来解决。直接复制边缘的像素进行扩充是不合适的。因为这会导致带有剑锋的领域直方图。
3. AHE的属性
- 领域的大小是该方法的一个参数。领域小，对比度得到增强，领域大，则对比度降低。
- 当某个区域包含的像素值非常相似，其直方图就会尖状化，此时直方图的变换函数会将一个很窄范围内的像素映射到整个像素范围。这将使得某些平坦区域中的少量噪音经AHE处理后过度放大。
二、限制对比度自适应直方图均衡（Contrast Limited Adaptive histgram equalization/CLAHE)
1.简述
CLAHE同普通的自适应直方图均衡不同的地方主要是其对比度限幅。这个特性也可以应用到全局直方图均衡化中，即构成所谓的限制对比度直方图均衡（CLHE），但这在实际中很少使用。在CLAHE中，对于每个小区域都必须使用对比度限幅。CLAHE主要是用来克服AHE的过度放大噪音的问题。
这主要是通过限制AHE算法的对比提高程度来达到的。在指定的像素值周边的对比度放大主要是由变换函数的斜度决定的。这个斜度和领域的累积直方图的斜度成比例。CLAHE通过在计算CDF前用预先定义的阈值来裁剪直方图以达到限制放大幅度的目的。这限制了CDF的斜度因此，也限制了变换函数的斜度。直方图被裁剪的值，也就是所谓的裁剪限幅，取决于直方图的分布因此也取决于领域大小的取值。
通常，直接忽略掉那些超出直方图裁剪限幅的部分是不好的，而应该将这些裁剪掉的部分均匀的分布到直方图的其他部分。如下图所示。
这个重分布的过程可能会导致那些倍裁剪掉的部分由重新超过了裁剪值（如上图的绿色部分所示）。如果这不是所希望的，可以不带使用重复不的过程指导这个超出的部分已经变得微不足道了。
2. 通过插值加快计算速度
如上所述的直接的自适应直方图，不管是否带有对比度限制，都需要对图像中的每个像素计算器领域直方图以及对应的变换函数，这使得算法及其耗时。
而插值使得上述算法效率上有极大的提升，并且质量上没有下降。首先，将图像均匀分成等份矩形大小，如下图的右侧部分所示（8行8列64个块是常用的选择）。然后计算个块的直方图、CDF以及对应的变换函数。这个变换函数对于块的中心像素（下图左侧部分的黑色小方块）是完全符合原始定义的。而其他的像素通过哪些于其临近的四个块的变换函数插值获取。位于图中蓝色阴影部分的像素采用双线性查插值，而位于便于边缘的（绿色阴影）部分采用线性插值，角点处（红色阴影处）直接使用块所在的变换函数。
这样的过程极大的降低了变换函数需要计算的次数，只是增加了一些双线性插值的计算量。
clahe.h
//*****************************************************************************
// Contrast Limited Adaptive Histogram Equalization (CLAHE) for OpenCV
//-----------------------------------------------------------------------------
// Original CLAHE implementation by Karel Zuiderveld, karel@cv.ruu.nl
// in "Graphics Gems IV", Academic Press, 1994.
//-----------------------------------------------------------------------------
// Converted to OpenCV format by Toby Breckon, toby.breckon@cranfield.ac.uk
// Copyright (c) 2009 School of Engineering, Cranfield University
// License : LGPL - http://www.gnu.org/licenses/lgpl.html
//-----------------------------------------------------------------------------
// Improved by Shervin Emami on 17th Nov 2010, shervin.emami@gmail.com
// http://www.shervinemami.co.cc/
//*****************************************************************************
#include <cv.h> // open cv general include file
// *****************************************************************************
// CLAHE input/output range flag definitions
#define CV_CLAHE_RANGE_FULL 0
#define CV_CLAHE_RANGE_INPUT 1
// *****************************************************************************
// cvAdaptEqualize(src, dst, xdivs, ydivs, bins)
//
// src - pointer to source image (must be single channel 8-bit)
// dst - pointer to destination image (must be single channel 8-bit)
// xdivs - number of cell divisions to use in vertical (x) direction (MIN=2 MAX = 16)
// ydivs - number of cell divisions to use in vertical (y) direction (MIN=2 MAX = 16)
// bins - number of histogram bins to use per division
// range - either of CV_CLAHE_RANGE_INPUT or CV_CLAHE_RANGE_FULL to limit the output
// pixel range to the min/max range of the input image or the full range of the
// pixel depth (8-bit in this case)
void cvAdaptEqualize(IplImage *src, IplImage *dst,
unsigned int xdivs, unsigned int ydivs, unsigned int bins, int range);
// cvCLAdaptEqualize(src, dst, xdivs, ydivs, bins, limit)
//
// src - pointer to source image (must be single channel 8-bit)
// dst - pointer to destination image (must be single channel 8-bit)
// xdivs - number of cell divisions to use in vertical (x) direction (MIN=2 MAX = 16)
// ydivs - number of cell divisions to use in vertical (y) direction (MIN=2 MAX = 16)
// bins - number of histogram bins to use per division
// limit - contrast limit for localised changes in contrast
// range - either of CV_CLAHE_RANGE_INPUT or CV_CLAHE_RANGE_FULL to limit the output
// pixel range to the min/max range of the input image or the full range of the
// pixel depth (8-bit in this case)
void cvCLAdaptEqualize(IplImage *src, IplImage *dst,
unsigned int xdivs, unsigned int ydivs, unsigned int bins, float limit,
int range);
// *****************************************************************************
/*
redefine : CV_ERROR_LOCAL macro unconditionally raises error with passed code and message.
After raising error, control will be transferred to the exit label.
*/
#undef CV_ERROR
#define CV_ERROR( Code, Msg ) \
{ \
cvError( (Code), "CLAHE code", Msg, __FILE__, __LINE__ ); \
exit(1); \
}
// *****************************************************************************
clahe.cpp
//*****************************************************************************
// Contrast Limited Adaptive Histogram Equalization (CLAHE) for OpenCV
//-----------------------------------------------------------------------------
// Original CLAHE implementation by Karel Zuiderveld, karel@cv.ruu.nl
// in "Graphics Gems IV", Academic Press, 1994.
//-----------------------------------------------------------------------------
// Converted to OpenCV format by Toby Breckon, toby.breckon@cranfield.ac.uk
// Copyright (c) 2009 School of Engineering, Cranfield University
// License : LGPL - http://www.gnu.org/licenses/lgpl.html
//-----------------------------------------------------------------------------
// Improved by Shervin Emami on 17th Nov 2010, shervin.emami@gmail.com
// http://www.shervinemami.co.cc/
//*****************************************************************************
#include <cv.h> // open cv general include file
#include "clahe.h" // opencv CLAHE include file
#include <stdio.h>
#include <stdlib.h>
// *****************************************************************************
// type defs. for Graphic Gemms Code - see later
#define BYTE_IMAGE
#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 Graphic Gemms CLAHE function. *********************/
static 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);
// *****************************************************************************
// General Notes:
//
// The number of "effective" greylevels in the output image is set by bins; 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 (?? is this correct ) results in
// standard (non-contrast limited) AHE.
// cvAdaptEqualize(src, dst, xdivs, ydivs, bins)
//
// perform adaptive histogram equalization (AHE)
//
// src - pointer to source image (must be single channel 8-bit)
// dst - pointer to destination image (must be single channel 8-bit)
// xdivs - number of cell divisions to use in vertical (x) direction (MIN=2 MAX = 16)
// ydivs - number of cell divisions to use in vertical (y) direction (MIN=2 MAX = 16)
// bins - number of histogram bins to use per division
// range - either of CV_CLAHE_RANGE_INPUT or CV_CLAHE_RANGE_FULL to limit the output
// pixel range to the min/max range of the input image or the full range of the
// pixel depth (8-bit in this case)
void cvAdaptEqualize(IplImage *src, IplImage *dst,
unsigned int xdivs, unsigned int ydivs, unsigned int bins, int range)
{
// call CLAHE with negative limit (as flag value) to perform AHE (no limits)
cvCLAdaptEqualize(src, dst, xdivs, ydivs, bins, -1.0, range);
}
// cvCLAdaptEqualize(src, dst, xdivs, ydivs, bins, limit)
//
// perform contrast limited adaptive histogram equalization (CLAHE)
//
// src - pointer to source image (must be single channel 8-bit)
// dst - pointer to destination image (must be single channel 8-bit)
// xdivs - number of cell divisions to use in vertical (x) direction (MIN=2 MAX = 16)
// ydivs - number of cell divisions to use in vertical (y) direction (MIN=2 MAX = 16)
// bins - number of histogram bins to use per division (MIN=2 MAX = 256)
// limit - contrast limit for localised changes in contrast
// (limit >= 0 gives standard AHE without contrast limiting)
// range - either of CV_CLAHE_RANGE_INPUT or CV_CLAHE_RANGE_FULL to limit the output
// pixel range to the min/max range of the input image or the full range of the
// pixel depth (8-bit in this case)
void cvCLAdaptEqualize(IplImage *src, IplImage *dst,
unsigned int xdivs, unsigned int ydivs, unsigned int bins,
float limit, int range)
{
double min_val, max_val;
unsigned char min, max;
// CV_FUNCNAME( "cvCLAdaptEquualize" );
// __BEGIN__;
// check the inputs to the function
int type;
if ((src == NULL) || (dst == NULL))
CV_ERROR( CV_StsUnsupportedFormat, "NULL value passed as image to function" );
CV_CALL( type = cvGetElemType( src ));
if( type != CV_8UC1 )
CV_ERROR( CV_StsUnsupportedFormat, "Only 8uC1 images are supported" );
CV_CALL( type = cvGetElemType( src ));
if( type != CV_8UC1 )
CV_ERROR( CV_StsUnsupportedFormat, "Only 8uC1 images are supported" );
//if( !CV_ARE_SIZES_EQ( src, dst )) // Modified by Shervin Emami, 17Nov2010.
if (src->width != dst->width || src->height != dst->height)
CV_ERROR( CV_StsUnmatchedSizes, "src and dst images must be equal sizes" );
if (((xdivs < 2) || (xdivs > 16)) || ((ydivs < 2) || (ydivs > 16)))
CV_ERROR( CV_StsBadFlag, "xdivs and ydivs must in range (MIN=2 MAX = 16)" );
if ((bins < 2) || (bins > 256))
CV_ERROR( CV_StsBadFlag, "bins must in range (MIN=2 MAX = 256)" );
// copy src to dst for in-place CLAHE.
cvCopy(src, dst);
// If the dimensions of the image are not a multiple of the xdivs and ydivs, then enlarge the image to be a correct size and then shrink the image.
// Also make sure the image is aligned to 8 pixels width, so that OpenCV won't add extra padding to the image.
// Added by Shervin Emami, 17Nov2010.
int enlarged = 0;
int origW = dst->width;
int origH = dst->height;
IplImage *tmpDst = 0;
if ((dst->width & (8-1)) || (dst->height & (8-1)) || (dst->width % xdivs) || (dst->height % ydivs)) {
int newW = ((dst->width + 8-1) & -8); // Align to 8 pixels, so that widthStep hopefully equals width.
newW = ((newW + xdivs-1) & -xdivs); // Also align for CLAHE.
int newH = ((dst->height + ydivs-1) & -ydivs);
//std::cout << "w=" << dst->width << ", h=" << dst->height << ". new w = " << newW << ", h = " << newH << std::endl;
IplImage *enlargedDst = cvCreateImage(cvSize(newW, newH), dst->depth, dst->nChannels);
cvResize(dst, enlargedDst, CV_INTER_CUBIC);
//cvReleaseImage(&dst);
tmpDst = dst;
dst = enlargedDst; // Use the enlarged image instead of the original dst image.
enlarged = 1; // signal that we need to shrink later!
}
// Check if OpenCV adds padding bytes on each row. Added by Shervin Emami, 17Nov2010.
if (dst->width != dst->widthStep)
CV_ERROR( CV_StsBadFlag, "dst->widthStep should be the same as dst->width. The IplImage has padding, so use a larger image." );
// check number of xdivs and ydivs is a multiple of image dimensions
if (dst->width % xdivs)
CV_ERROR( CV_StsBadFlag, "xdiv must be an integer multiple of image width " );
if (dst->height % ydivs)
CV_ERROR( CV_StsBadFlag, "ydiv must be an integer multiple of image height " );
// get the min and max values of the image
if (range == CV_CLAHE_RANGE_INPUT) {
cvMinMaxLoc(dst, &min_val, &max_val);
min = (unsigned char) min_val;
max = (unsigned char) max_val;
} else {
min = 0;
max = 255;
}
// call CLHAHE for in-place CLAHE
//int rcode =
CLAHE((kz_pixel_t*) (dst->imageData), (unsigned int) dst->width, (unsigned int)
dst->height, (kz_pixel_t) min, (kz_pixel_t) max, (unsigned int) xdivs, (unsigned int) ydivs,
(unsigned int) bins, (float) limit);
//printf("RCODE %i\n", rcode);
// If the dst image was enlarged to fit the alignment, then shrink it back now.
// Added by Shervin Emami, 17Nov2010.
if (enlarged) {
//std::cout << "w=" << dst->width << ", h=" << dst->height << ". orig w=" << origW << ", h=" << origH << std::endl;
cvResize(dst, tmpDst, CV_INTER_CUBIC); // Shrink the enlarged image back into the original dst image.
cvReleaseImage(&dst); // Free the enlarged image.
}
}
// *****************************************************************************
/*
* 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)
*/
/*
EULA: The Graphics Gems code is copyright-protected. In other words, you cannot
claim the text of the code as your own and resell it. Using the code is permitted
in any program, product, or library, non-commercial or commercial. Giving credit
is not required, though is a nice gesture. The code comes as-is, and if there are
any flaws or problems with any Gems code, nobody involved with Gems - authors,
editors, publishers, or webmasters - are to be held responsible. Basically,
don't be a jerk, and remember that anything free comes with no guarantee.
- http://tog.acm.org/resources/GraphicsGems/ (August 2009)
*/
/*********************** Local prototypes ************************/
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 static unsigned int uiMAX_REG_X = 16; /* max. # contextual regions in x-direction */
const static unsigned int uiMAX_REG_Y = 16; /* max. # contextual regions in y-direction */
/************************** main function CLAHE ******************/
static 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 #TPB FIX */
#ifndef BYTE_IMAGE /* #TPB FIX */
if (Max >= uiNR_OF_GREY) return -5; /* maximum too large */
#endif
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);
}
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 */
}
void ClipHistogram (unsigned long* pulHistogram, unsigned int
uiNrGreylevels, unsigned long ulClipLimit)
/* 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* pulBinPointer, *pulEndPointer, *pulHisto;
unsigned long ulNrExcess, ulUpper, ulBinIncr, ulStepSize, i;
unsigned long ulOldNrExcess; // #IAC Modification
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;
}
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 */
// }
//}
/* ####
IAC Modification:
In the original version of the loop below (commented out above) it was possible for an infinite loop to get
created. If there was more pixels to be redistributed than available space then the
while loop would never end. This problem has been fixed by stopping the loop when all
pixels have been redistributed OR when no pixels where redistributed in the previous iteration.
This change allows very low clipping levels to be used.
*/
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));
/* Finish loop when we have no more pixels or we can't redistribute any more pixels */
}
void 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];
}
}
void MapHistogram (unsigned long* pulHistogram, kz_pixel_t Min, kz_pixel_t Max,
unsigned int uiNrGreylevels, unsigned long ulNrOfPixels)
/* This function calculates the equalized lookup table (mapping) by
* cumulating the input histogram. Note: lookup table is rescaled in range [Min..Max].
*/
{
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;
}
}
void MakeLut (kz_pixel_t * pLUT, kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrBins)
/* To speed up histogram clipping, the input image [Min,Max] is scaled down to
* [0,uiNrBins-1]. This function calculates the LUT.
*/
{
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 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);
}
}
}
}
// *****************************************************************************
adapthistequal.cpp
// Example : adaptive histogram equalise grayscale image
// usage: prog {<image_name> | <video_name>}
// Author : Toby Breckon, toby.breckon@cranfield.ac.uk
// Copyright (c) 2009 School of Engineering, Cranfield University
// License : LGPL - http://www.gnu.org/licenses/lgpl.html
#include "cv.h" // opencv general include file
#include "highgui.h" // opencv GUI include file
#include "clahe.h" // opencv CLAHE include file
#include "useful_macros.h" // opencv macros include file
#include <stdio.h>
/******************************************************************************/
// setup the camera index properly based on OS platform
// 0 in linux gives first camera for v4l
//-1 in windows gives first device or user dialog selection
#ifdef linux
#define CAMERA_INDEX 0
#else
#define CAMERA_INDEX -1
#endif
/******************************************************************************/
// function that takes a gray scale image and draws a histogram
// image for it in a pre-allocated image
void create_histogram_image(IplImage* grayImg, IplImage* histogramImage){
CvHistogram *hist = NULL; // pointer to histogram object
float max_value = 0; // max value in histogram
int hist_size = 256; // size of histogram (number of bins)
int bin_w = 0; // initial width to draw bars
float range_0[]={0,256};
float* ranges[] = { range_0 };
hist = cvCreateHist(1, &hist_size, CV_HIST_ARRAY, ranges, 1);
cvCalcHist( &grayImg, hist, 0, NULL );
cvGetMinMaxHistValue( hist, 0, &max_value, 0, 0 );
cvScale( hist->bins, hist->bins, ((double)histogramImage->height)/max_value, 0 );
cvSet( histogramImage, cvScalarAll(255), 0 );
bin_w = cvRound((double)histogramImage->width/hist_size);
for(int i = 0; i < hist_size; i++ )
{
cvRectangle( histogramImage, cvPoint(i*bin_w, histogramImage->height),
cvPoint((i+1)*bin_w, histogramImage->height
- cvRound(cvGetReal1D(hist->bins,i))),
cvScalarAll(0), -1, 8, 0 );
}
cvReleaseHist (&hist);
}
/******************************************************************************/
int main( int argc, char** argv )
{
IplImage* img = NULL; // image object
CvCapture* capture = NULL; // capture object
IplImage* eqHistogramImage = NULL; // histogram images
IplImage* grayHistogramImage = NULL;
const char* windowName = "CL Adaptive Histogram Equalization"; // window name
const char* windowName1 = "Grayscale"; // window name
const char* windowNameH1 = "Adaptive Equalised Histogram"; // window name
const char* windowNameH2 = "Original Histogram"; // window name
bool keepProcessing = true; // loop control flag
char key; // user input
int EVENT_LOOP_DELAY = 40; // delay for GUI window
// 40 ms equates to 1000ms/25fps = 40ms per frame
int xdivs = 2;
int ydivs = 2;
int bins = 256;
int limit_counter = 14;
// if command line arguments are provided try to read image/video_name
// otherwise default to capture from attached H/W camera
if(
( argc == 2 && (img = cvLoadImage( argv[1], 1)) != 0 ) ||
( argc == 2 && (capture = cvCreateFileCapture( argv[1] )) != 0 ) ||
( argc != 2 && (capture = cvCreateCameraCapture( 0 )) != 0 )
)
{
// create window object (use flag=0 to allow resize, 1 to auto fix size)
cvNamedWindow(windowName, 1); // flag set to 1 by Shervin Emami, 17Nov2010.
cvNamedWindow(windowName1, 1); // flag set to 1 by Shervin Emami, 17Nov2010.
cvNamedWindow(windowNameH1, 1); // flag set to 1 by Shervin Emami, 17Nov2010.
cvNamedWindow(windowNameH2, 1); // flag set to 1 by Shervin Emami, 17Nov2010.
cvNamedWindow("Simple Histogram Equalization", 1); // Added by Shervin Emami, 17Nov2010.
cvCreateTrackbar("X cells", windowName, &xdivs, 16, NULL);
cvCreateTrackbar("Y cells", windowName, &ydivs, 16, NULL);
cvCreateTrackbar("bins", windowName, &bins, 256, NULL);
cvCreateTrackbar("limit (x 0.1)", windowName, &limit_counter, 30, NULL);
// define required images for intermediate processing
// (if using a capture object we need to get a frame first to get the size)
if (capture) {
// cvQueryFrame s just a combination of cvGrabFrame
// and cvRetrieveFrame in one call.
img = cvQueryFrame(capture);
if(!img){
if (argc == 2){
printf("End of video file reached\n");
} else {
printf("ERROR: cannot get next fram from camera\n");
}
exit(0);
}
}
IplImage* grayImg =
cvCreateImage(cvSize(img->width,img->height), img->depth, 1);
grayImg->origin = img->origin;
IplImage* eqImg =
cvCreateImage(cvSize(img->width,img->height), img->depth, 1);
eqImg->origin = img->origin;
eqHistogramImage = cvCreateImage(cvSize(255,200), 8, 1);
grayHistogramImage = cvCreateImage(cvSize(255,200), 8, 1);
// start main loop
while (keepProcessing) {
// if capture object in use (i.e. video/camera)
// get image from capture object
if (capture) {
// cvQueryFrame s just a combination of cvGrabFrame
// and cvRetrieveFrame in one call.
img = cvQueryFrame(capture);
if(!img){
if (argc == 2){
printf("End of video file reached\n");
} else {
printf("ERROR: cannot get next fram from camera\n");
}
exit(0);
}
} else {
// if not a capture object set event delay to zero so it waits
// indefinitely (as single image file, no need to loop)
EVENT_LOOP_DELAY = 0;
}
// *** Histogram Equalisation Processing
// if input is not already grayscale, convert to grayscale
if (img->nChannels > 1){
cvCvtColor(img, grayImg, CV_BGR2GRAY);
} else {
grayImg = img;
}
// histogram equalize it (checking first for valid cells numbers)
// if ((img->width % xdivs) || (img->height % ydivs))
if (0) // Check has been removed by Shervin Emami, 17Nov2010.
{
printf("X cells and Y cells must be multiples of image height and image width\n\n");
} else {
CV_TIMER_START(X)
cvCLAdaptEqualize(grayImg, eqImg, (unsigned int) xdivs, (unsigned int) ydivs,
(unsigned int) bins, (float) limit_counter * 0.1, CV_CLAHE_RANGE_FULL);
CV_TIMER_STOP(X, "CLAHE")
}
// ***
// *** draw histograms
create_histogram_image(grayImg, grayHistogramImage);
create_histogram_image(eqImg, eqHistogramImage);
// display image in window
cvShowImage( windowName, eqImg );
cvShowImage( windowName1, grayImg );
cvShowImage( windowNameH1, eqHistogramImage );
cvShowImage( windowNameH2, grayHistogramImage );
// Simple Histogram Equalization. Added by Shervin Emami, 17Nov2010.
IplImage *histEqImage = cvCreateImage(cvGetSize(grayImg), grayImg->depth, grayImg->nChannels);
cvEqualizeHist(grayImg, histEqImage);
cvShowImage( "Simple Histogram Equalization", histEqImage );
cvReleaseImage(&histEqImage);
// start event processing loop (very important,in fact essential for GUI)
// 4 ms roughly equates to 100ms/25fps = 4ms per frame
key = cvWaitKey(EVENT_LOOP_DELAY);
if (key == 'x' || key == 27) { // 'Esc' key added by Shervin Emami, 17Nov2010.
// if user presses "x" or 'Esc' then exit
printf("Keyboard exit requested : exiting now - bye!\n");
keepProcessing = false;
}
}
// destroy window objects
// (triggered by event loop *only* window is closed)
cvDestroyAllWindows();
// destroy image object (if it does not originate from a capture object)
if (!capture){
cvReleaseImage( &img );
}
cvReleaseImage( &grayImg );
cvReleaseImage( &eqImg );
// release histogram images
cvReleaseImage( &eqHistogramImage );
cvReleaseImage( &grayHistogramImage );
// all OK : main returns 0
return 0;
}
// not OK : main returns -1
return -1;
}
/******************************************************************************/
3478
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
